Introduction

At present, three primary methods are applied for calculating planing vessel motions, i.e. model experiment (Wencai Dong et al., 2004, Yang and Gao, 2008), 2.5D method, and CFD method. Except the model experiment, most of the theoretical/numerical studies are based on the slender body assumption (Tulin, 1956, Savander, 1997, Zhao et al., 1997). The 2.5D method is a kind of slender body theory which provides a means of simplifying a 3D shape into a series of 2D sections. The problem of a 2D section with knuckles was analyzed by Arai et al. (1994). The hydrodynamic properties of prismatic planing hulls were investigated through using the 2D+T method by Battistin and Iafrati. (2003). In these years, the CFD method has been widely used in ship design, and many pre-designs completed through it. Although the experimental approach is still useful it has its own restrictions and the model test is also very expensive. For planing hull, a large systematic series of experimental tests have been done. Depending on the experimental tests results, some attempts have been made and several relations were developed for the calculation of hydrodynamic forces acting on the planing hull. Through these experimental results a series of semi-empirical methods were set up and widely employed for the performance estimations of planing crafts. Nowadays, with the developments of computer and computational fluid dynamics, numerical simulations are becoming a common way for ship design. Many researchers have made a large number of comparisons between the CFD method and the experimental test to prove the accuracy of the numerical method. In order to prove the accuracy of CFD codes in the prediction of planing surfaces hydrodynamic characteristics, Brizzolara and Serra(2007) have made the comparison between the CFD results and the experimental results. The comparative results showed that the CFD method is able to reproduce the physical phenomena of the free sur face flow in proximity of a planing hull in a sufficient and accurate way. Except from a few particular cases, the error between lift, drag, and trim moment predicted by the CFD code remain well within the measuring error that affects the experimental results.

If we want to obtain the hydrodynamic performance of planing hull using numerical methods, the most important work is that we must confirm the running attitude for different speeds beforehand. As mentioned earlier, for most situations we confirm the running attitude through experiments. In this paper, however, Goal Driven Optimization (GDO) was applied for confirming the attitude and the hydrodynamic characteristics of planing hull. It was calculated under the obtained equilibrium attitude. It is obvious that the lift and lift center are the main conditions for keeping the planing hull in the equilibrium state. Therefore, we define lift and lift center as the equilibrium conditions. At the same time, response surface method will also be applied, and then we will use Multiple Objective Genetic Algorithm (MOGA) to find the equilibrium positions of every speeds which meet the equilibrium conditions well. The results obtained from the MOGA method will also be recalculated using CFX to reconfirm the MOGA results. At last, the hydrodynamic characteristics of planing hull were calculated by using ANSYS-CFX, then the calculated results and the experiment results were compared to improve the accuracy of this method. In order to illustrate the procedure clearly, the series 62 4667-1 planing hull model and the experiment results (Savitsky, 1992) are applied here. The geometry of the hull is built in Rhino 3D and the mesh of the fluid field is generated in ANSYS-MESH.

Materials And Method

Governing equations

The continuity equation for the mass conservation and the Navier–Stokes equations for the momentum transport govern the 3D motion of the incompressible and viscous fluid. These equations are expressed in terms by the set of partial differential equations as follows:

where ρis the density, P is the pressure, u_i is the velocity of the xi direction, τ_ij is the stress tensor, and is the component of the gravitational acceleration in the direction of the Cartesian coordinate x_i. In the case of constant density and gravity, the term ρg can be written as grad (ρg·r), where r is the position vector, r=x_i, i_i(usually, gravity is assumed to act in the negative z-direction, i.e. g_z=gZ_z, g_z being negative; in this case g_z·r=gz_z). Then -ρg·z_z is the hydrostatic pressure, which makes it convenient for numerical solutions to more efficiently define = P-ρg·z_z as the head and use it in place of the pressure. The term ρg_i then disappears from the above equation. If the actual pressure is needed, we only have to add ρgz to

The homogenous multiphase Eulerian fluid approach utilizes the Volume of Fluid (VOF) method to describe the free surface flow problem mathematically. The VOF method developed by Hirt and Nichols (1981) is a fixed mesh technique designed for two or more fluids, where in each cell of a mesh it is necessary to use only one value for each dependent variable defining the fluid state. The variables can be defined by a function α_q (where q represents the fraction of a volume; q = 1,2… ). The value of this variable is one at any point occupied by the fluid and is zero in other cases, it is written as:

the average value of in a cell represents the position of the interface of the fluid. In particular, a unity value of corresponds to a cell full of fluid, while a zero value indicates that the cell is empty. So, the cells with α values between zero and one must contain a free surface. The tracking of the interface is obtained by solving the continuity equation of the volume fraction. For the q liquid, troughs of the density are as follows:

For a Newtonian fluid the viscous stress tensor can be expressed as:

where is the strain-rate tensor and is the dynamic viscosity. For an incompressible fluid, Equation 5 reduces to:

the Navier-Stokes equations are a highly non-linear system. The strong non-linearity of the equations produces high frequency oscillations when the Reynolds number is increased, and the flow becomes unstable and turbulent. It is computationally very expensive to solve the equations directly, which makes it so that presently, only in very simple geometry configurations it is possible to solve the Navier-Stokes equations using direct methods (DNS). The most common approach at the moment in the hydraulic engineering practice is to solve the Reynolds Averaged Navier-Stokes (Versteeg HK and Malalasekera. W 1995; Anderson John David and Wendt J 1995). In simple cases the continuity and Navier–Stokes equations can be solved analytically. More complex flows can be tackled numerically with CFD techniques such as the finite volume method without additional approximations. Many, if not most, flows of engineering significance are turbulent, so the turbulent flow regime is not just of theoretical interest. Fluid engineers need access to viable tools capable of representing the effects of turbulence. In order to be able to compute turbulent flows with the RANS equations it is necessary to develop turbulence models to predict the Reynolds stresses and the scalar transport terms and close the system of mean flow equations. One of the most prominent turbulence models, the (k-epsilon) model, has been implemented in most general purpose CFD codes and is considered the industry standard model. It has proven to be stable and numerically robust and has a well established regime of predictive capability. For general purpose simulations, the k-e model offers a good compromise in terms of accuracy and robustness. Therefore, the k-e model is also applied in this research. Within CFX, the k-e turbulence model uses the scalable wall-function approach to improve robustness and accuracy when the near-wall mesh is very fine. The scalable wall functions allow solution on arbitrarily fine near wall grids, which is a significant improvement over standard wall functions

Goal Driven Optimization

Goal Driven Optimization (GDO) is a set of constrained multi-objective optimization (MOO) techniques, in which the "best" possible designs are obtained from a sample set, given the goals you set for parameters. There are three available optimization methods: Screening, MOGA, and NLPQL. MOGA and NLPQL can only be used when all input parameters are continuous. The GDO process allows you to determine the effect on input parameters with certain objectives applied for the output parameters. For example, in a structural engineering design problem, you may want to determine which set of designs (in terms of geometric problem dimensions and material types) best satisfy minimum mass, maximum natural frequency, maximum buckling and shear strengths, and minimum cost, with maximum value constraints on the von Mises stress and maximum displacement.

GDO can be used for design optimization in three ways: the Screening approach, the MOGA approach, or the NLPQL approach. The Screening approach is a non-iterative direct sampling method by a quasi-random number generator based on the Hammersley algorithm. The MOGA approach is an iterative Multi-Objective Genetic Algorithm, which can optimize problems with continuous input parameters. NLPQL is a gradient-based single objective optimizer that is based on quasi-Newton methods.

The concept of Genetic Algorithm (GA) was developed by Holland and his colleagues in the 1960s and 1970s. GA are inspired by the evolutionist theory explaining the origin of species. In nature, weak and unfit species within their environment are faced with extinction by natural selection. The strong ones have greater opportunity to pass their genes to future generations via reproduction. In the long run, species carrying the correct combination in their genes become dominant in their population. Sometimes, during the slow process of evolution, random changes may occur in genes. If these changes provide additional advantages in the challenge for survival, new species evolve from the old ones. Unsuccessful changes are eliminated by natural selection.

In GA terminology, a solution vector x∈X is called an individual or a chromosome. Chromosomes are made of discrete units called genes. Each gene controls one or more features of the chromosome. In the original implementation of GA by Holland, genes are assumed to be binary digits. In later implementations, more varied gene types have been introduced. Normally, a chromosome corresponds to a unique solution x in the solution space. This requires a mapping mechanism between the solution space and the chromosomes. This mapping is called an encoding. In fact, GA work on the encoding of a problem, not on the problem itself.

GA operate with a collection of chromosomes, called a population. The population is normally randomly initialized. As the search evolves, the population includes fitter and fitter solutions, and eventually it converges, meaning that it is dominated by a single solution. Holland also presented a proof of convergence (the schema theorem) to the global optimum where chromosomes are binary vectors.

GA use two operators to generate new solutions from existing ones: crossover and mutation. The crossover operator is the most important operator of GA. In crossover, generally two chromosomes, called parents, are combined together to form new chromosomes, called offspring. The parents are selected among existing chromosomes in the population with preference towards fitness so that offspring is expected to inherit good genes which make the parents fitter. By iteratively applying the crossover operator, genes of good chromosomes are expected to appear more frequently in the population, eventually leading to convergence to an overall good solution.

The mutation operator introduces random changes into characteristics of chromosomes. Mutation is generally applied at the gene level. In typical GA implementations, the mutation rate (probability of changing the properties of a gene) is very small and depends on the length of the chromosome. Therefore, the new chromosome produced by mutation will not be very different from the original one. Mutation plays a critical role in GA. As discussed earlier, crossover leads the population to converge by making the chromosomes in the population alike. Mutation reintroduces genetic diversity back into the population and assists the search escape from local optima. (Konak et al., 2006)

The use of genetic algorithms for multi-objective optimization has been growing considerably in the past few years. Gas were originally used for maximization or minimization of an unconstrained function. However, there has been increasing interest in optimizing two or more criteria simultaneously, especially if it is difficult to represent one criteria in terms of another. These problems are often referred to as multi-objective (or vector-valued) optimizations problems. One such field of study utilizing this concept is the aerospace industry, where an effort has been made to incorporate cost directly into the design process. This methodology can lead to high performance designs that can be built with available materials and manufacturing techniques. Furthermore, such studies can be used to formulate trade o studies between cost and weight which may aid in the selection of a design that minimizes cost and/or weight, two of the most important considerations in aerospace applications.

The goal of single objective optimization problems is straightforward: find the maximum or minimum value of a function for a given set of parameters. The optimization concept is less clear for multi-objective problems, since the best value for one objective usually does not imply that the other objective(s) is simultaneously optimized. Thus, the concept of Pareto-optimality is often used in multi-objective problems to help determine the best way to simultaneously satisfy all objectives to the greatest extent possible. (Soremekun Grant AE, 1997)

Analysis procedure

In this part, we will make the geometry model of the series 62 4667-1 using Rhino 3D and generate the mesh of the fluid field in the ANSYS Mesh, and then the calculations for different speeds were done using the ANSYS- CFX software. In order to prove the accuracy of CFD method for calculation of planing hull performance, the comparison between CFX results and experiment results is expected in the last part.

Geometry

Firstly, the 3D model is built using Rhino 3D. With the help of the main principal dimensions from Blount, (2011), the hull model is made in Rhino 3D and the gravity center of the hull is placed at the origin point of Rhino. The final model of the hull is shown in Fig.1 below.

Before starting the analysis, the accuracy of this model must first be confirmed. The results of the comparison between series 62 4667-1, mentioned by Blount, (2011), and Rhino model mentioned in Fig. 1 are shown in Table 1.

Where Ap is the projected area of the hull bottom bounded by the chine and transom, Bpx is the maximum projected chine beam and Bpt is the chine beam at transom. As shown in the Table 1, the main principal dimensions of the Rhino model have been modeled exactly which proves that it could be used for the next analysis.

The Rhino model was imported into the ANSYS Design Modeler (DM) and it was modified. Before we generate the fluid domain the hull model must be transformed into a solid body in the DM, then the domain was cut with the hull model to generate the fluid domain around the hull. More details are shown in Fig. 2.

Here, we divide the domain into two bodies so that we could generate a relatively better mesh region to catch the free surface motion using "Inflation" at interface. More details are explained in the following parts. The domain size is defined in the DM shown in Fig. 3.

Meshing

In order to avoid the influence of tetrahedron mesh character on the free surface which could make the free surface very rough, we inflate 20 layers downward on the interface for achieving an enough space to show the free surface transformation as shown in Fig. 4. Now, the mesh type around the interface has been transformed into triangular prism.

Furthermore, the first layer’s thickness of mesh around the hull should especially be considered. According to the existing experience and theory, the flow could be divided into core flow and near-wall flow in the normal direction on the wall. In the core areas, flow is a fully developed turbulence and we can put it aside for the moment. In the near-wall flow, the fluid motion is obviously influenced by the wall, and because of that we use a set of experience formulas to calculate the near-wall flow, which is called the "Wall Function". Therefore, the first layer node of mesh must be placed into the core region (Anderson John David and Wendt J, 1995).

The first layer’s thickness is calculated as follows (ITTC 2011):

where y⁺ is defined as 80 and the Δy is obtained as 1mm on the hull surface. The mesh of computational domain is shown below in Fig. 5.

Boundary conditions

The boundary conditions are defined after the fluid field mesh is imported into the CFX-Pre.

Boundary condition definitions are shown in Table 2 above. For hull surface the wall condition with no slip was used. Here ‘no slip’ means the effect of friction on the boundary is active in each direction which is also the effect to wall function. “Free slip” means that there is on friction on the boundary. For Inflow, the inlet velocities are constant at inlet boundary plane. For “opening”, the boundary condition allows the fluid to cross the boundary surface in either direction. For example, all of the fluid might flow into the domain at the opening, or all of the fluid might flow out of the domain, or a mixture of the two might occur. An opening boundary condition might be used where it is known that the fluid flows in both directions across the boundary. At “Top” the “Mass and Momentum” is specified as “Opening Pres. and Dir” which means that when the flow direction is into the domain, the pressure value is taken to be total pressure based on the normal component of velocity and when it is leaving the domain, it is taken to be relative static pressure. This is the most robust and stable setting because it puts a constraint on the momentum transported through the boundary condition. In this case, the VOF method is applied for describing the free surface which could be explained in CFX as follows:

CFX Expressions Language(CEL):

Calculation

In this part the Goal Driven Optimization (GDO) is applied to calculate the corresponding resistance and trim for the speeds from 4.87 to 13.74 knot. It is obvious that the weight and the gravity center of the hull are the two most important conditions for keeping the planing hull under the steady state. If the lift and lift center, which come from the hydrodynamic forces, could meet the weight and gravity center of the hull, we could affirm that the hull has been in an equilibrium state. And only in this circumstance are the results we obtained from CFD are receivable. Therefore, the most important work now is to find the equilibrium position of each speed. This work was done using the Goal Driven Optimization (GDO) module in ANSYS Workbench.

Besides, the "Response Surface" method is applied in this procedure to catch the relationship between "Input parameters" and "Output parameters". As we know, trim and draft are the main parameters affecting the planing hull running attitude. Therefore, these two parameters were defined as the "Input parameters". At the same time, as equilibrium conditions the lift and lift center were defined as "Output parameters". Now, let us take 5.86 knot for example to explain this procedure.

In order to generate the response surface we must calculate some feature points previously called "Design of Experiments" (DOE) in Workbench. In this step, the Central Composite Design-Face Centered-Standard is used, which could produce 9 DOE. After we produce 9 DOE, we will calculate these 9 DOE one by one using CFX.

Firstly, we should predict the approximate range of trim and draft. Because for most situations we have no experiment results in advance, the ranges of trim and draft have to be estimated according to experience. Therefore, we usually need to calculate from a large range and then narrow the range down step by step. But now, in order to show the process easily, the range is assumed around the experiment results. Thus, the range of trim and draft are defined as below:

0.9 * experiment trim≤Trim≤1.1 * experment trim

0.9 * experiment draft≤Draft≤1.1 * experiment draft

Where Trim is the degree between keel and horizontal line, (+) means stern trim in this paper. According to the results obtained by Savitsky (Savitsky, D, 1992) the trim and draft of 5.86 knot is 2.91 Deg and -0.06899 m. (The draft data has already been transformed in Rhino 3D depending on Blount, 2011) Finally, the lower bound and higher bound are defined

2.619≤Trim≤3.201

–0.075889≤–0.062091

At last, the DOE results obtained by CFX are shown in Fig. 6.

The response surface is generated through these DOE results by using Standard Response Surface-Full 2-nd Order Polynomials shown in Fig. 7.

In the last step, we use Multi-Objective Genetic Algorithms (MOGA) mentioned in “Goal Driven Optimization” to find the trim and draft of the equilibrium position, which is shown below in Fig. 8:

The equilibrium conditions are:

Lift = planing hull weight = 685.27 N

Lift center = planing hull gravity center = 0

Here, we calculate the three candidates one by one through CFX using the corresponding obtained trim and draft, and then select the result which could meet the equilibrium conditions most well as the final result and record the corresponding total resistance. The wave pattern for 5.86 knot is given below in Fig. 9.

The other velocities are calculated in the same way as mentioned earlier and the final results are provided in the last part.

Results of analysis

For planing hull, we are particularly concerned about the high speed situation. Thus, we just choose the speeds from 4.87 to 13.74 knot for the CFD analysis here. The experiment results, MOGA results and CFX calculated results based on MOGA are given in Table 3.

Where is Displacement-Froude number. Through comparison between experiment results and MOGA results we could find that the errors of trim, draft and total resistance for each velocity are all lower than 4%. For the CFX errors, they also did not exceed 5%. (Fig. 10 ~ Fig. 15) Fig. 11, 12, 13, 14

Results

In this paper, we calculated the hydrodynamic performance of the series 62 4667-1 planing hull through using the CFD and optimization methods in the ANSYS Workbench. The Goal Driven Optimization (GDO) caught the trim and draft of the planing hull for each velocity, and finally the equilibrium positions were confirmed. Then, we calculated the resistance under the equilibrium position from 4.87 to 13.74 knot. The results show that the MOGA errors of Rt in different speeds are lower than 4%. On the other hand, the CFD errors based on corresponding trim and draft in MOGA results show a slight increase than MOGA Rt in high speed. But we could also find that the verification results using CFD still show a good accuracy which proves that it is feasible to estimate the running attitude and corresponding resistance of planing hull through MOGA. Here, although we did not do any designs about planing hull using this method, the validity of it has been already proved through the results obtained earlier. Differ from other displacement ship, planing hull always changes attitude along with the speed, and we have to spend much more time and energy estimating the equilibrium position for calculating corresponding resistance. However, it is still a complex job for us. For traditional method, although the model test is considered as a relatively accurate method, the high-cost of it makes us feel embarrassed. And the estimation based on empirical data is always restricted by ship types. However, the method shown here combines CFD and MOGA to predict the performance of planing hull. Although CFD has been widely used in ship design, it is also week for searching the equilibrium position only by itself. Because it is impossible to calculate all likelihood. But the addition of MOGA raises efficiency sharply which gives us a mathematical guidance for acquiring the attitude of a planing hull. We could find some points met the equilibrium condition around response surface. It is obvious that it has an economic cost and also get rid of constrain of empirical method for completely different ship types. For the procedure of a new de sign of planing hull based on mother ship, this method could give us the credible predictions of running attitude and corresponding resistance. On the other hand, we could also apply it instead of an expensive model experiment when we want to carry out a brand new design. That could be considered as the most important practical value of this method.