Aerodynamic analysis and multi-objective optimization of wings in ground effect

Ocean Engineering 68 (2013) 1–13

Contents lists available at SciVerse ScienceDirect

Ocean Engineering journal homepage: www.elsevier.com/locate/oceaneng

Aerodynamic analysis and multi-objective optimization of wings in ground effect Sang-Hwan Lee a,1, Juhee Lee b,n a b

Department of Mechanical Engineering, Hanyang University, Seoul 133-791, South Korea Department of Mechatronics Engineering, Hoseo University, Chungnam 336-795, South Korea

art ic l e i nf o

a b s t r a c t

Article history: Received 19 January 2012 Accepted 6 April 2013 Available online 21 May 2013

Optimization of the sectional shapes of wings in ground effect (WIG) has been performed in this study by using computational fluid dynamics (CFD) and multi-objective optimization technology. The primary factors of the aerodynamic characteristics of the wings in ground effect are the lift force, the static height stability, and the lift-to-drag ratio. The strong trade-off among aerodynamic characteristics makes it difficult to simultaneously satisfy the design requirements of high aerodynamic performance and high stability. In this study, three characteristics – the lift coefficient, the aerodynamic center of height, and the lift-to-drag ratio – are chosen as the objective functions. The aerodynamic center of height is selected, rather than static height stability, because it is more suitable for single airfoil optimization and can be directly extended to a WIG vehicle. Also, 18 coordinates in the sectional airfoil are adopted as design variables. In multi-objective optimization, the optimal solutions are not unique, but are a set of nondominated and potential optima, called Pareto optima or Pareto sets. The Pareto optima (113 individuals) include various airfoil shapes such as a flat lower surface and a convex lower surface next to the trailing edge, which show high lift and high stability, respectively. & 2013 Elsevier Ltd. All rights reserved.

Keywords: Multi-objective optimization Ground effect Airfoil section Shape optimization Computational fluid dynamics Aerodynamic characteristics

1. Introduction Wing-in-ground (WIG) vehicles operate close to the water or ground surface by utilizing the air cushion of relatively highly pressurized air created between the airfoil and the surface. The air cushion augments lift and reduces drag considerably compared to an out-of-ground effect vehicle; these phenomena further enhance the lift-to-drag ratio, and have been observed during takeoff and landing by a number of researchers since the 1900s (e. g., Rozhdestvensky, 2006; Kikuchi et al., 1997; Joh and Kim, 2004; Kim and Joh, 2004). Studies have shown that the ground has a significant influence on the pressure distributions along the wing surface. As a WIG vehicle moves forward, the speed of the oncoming air gradually decreases under the lower wing surface, and dynamic pressure changes to static pressure. This increased pressure is called an air cushion or a ram effect, and it necessitates a longer runway for landing. However, this aerodynamic enhancement makes a WIG vehicle extremely efficient relative to other transportation options. No fast ship can match the speed of a WIG vehicle, and no economical aircraft can improve on the low operational expenses of a WIG n

Corresponding author. Tel.: +82 41 540 5956; fax: +82 41 540 5808. E-mail addresses: [email protected] (S.-H. Lee), [email protected], [email protected] (J. Lee). 1 Tel./fax: +82 2 2220 0445. 0029-8018/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.oceaneng.2013.04.018

vehicle. Therefore, WIG vehicles satisfy the niche between ships and aircrafts and may become an alternative means of transportation in the near future. Along with their advantages, WIG vehicles are challenged by technical difficulties, such as hump drag, which impedes the high speed needed to take off above the surface, and the instability problem, which is not generally observed in a typical airplane. In particular, a cambered airfoil, while able to achieve additional lift via the ground effect, does not easily satisfy static height stability (H.S.) when it is in ground effect. Attempts to improve the stability and controllability of a WIG effect vehicle have included employing a large horizontal tail, a tandem configuration (Rozhdestvensky, 2006), a sophisticated control system (Kornev and Matveev, 2003), and an S-shaped main wing (Shin et al., 2000), but the inherent challenges have generally hindered the development of WIG-vehicle industries. Until recently, experimental and numerical research on WIG vehicles has concentrated on the aerodynamic characteristics. For instance, Kikuchi et al. (1997) measured the forces on an airfoil from the interactions between airfoil and ground. Joh and Kim (2004) numerically studied the pressure change, drag, and stability around airfoils of a WIG vehicle for two- and three-dimensional configurations. Kornev and Matveev (2003) analyzed the static height stability using the vortex lattice method (VLM). Their study examined three important factors for static height stability: the tail unit, the profiles of the wing sections, and the main wing profile and they found that static height stability of a WIG vehicle could not be

2

S.-H. Lee, J. Lee / Ocean Engineering 68 (2013) 1–13

Nomenclature Bij c Cd Cl Cm Cp di, dj of f
line fe on
line fe F(X) h H:S: k

blending function of Bézier curve chord length drag coefficient lift coefficient moment coefficient at aerodynamic center pressure coefficient niche distance offline performance online performance objective function height at trailing edge static height stability turbulent kinetic energy

achieved by moving the center of gravity. They found that the favorable range of height stability for a stable flight with the ground effect was between −0.15 and −0.05. When the height stability was less than −0.15, the static stability became excessive and led to dynamic instability. When the static height stability was between 0 and −0.05, there was insufficient stability, which led to long-period (phugoid) instability. As a result, Kornev and Matveev concluded that, for a stable flight, the center of gravity should be located between the aerodynamic centers of height and pitch and close to the center of height. The optimal design of WIG airfoils has been studied by only a few researchers. Most have employed a single-objective optimization with a local optimization technique. In one study, Kim and Chun (1998) performed computational optimization for an airfoil shape. They chose the pressure distributions (inverse design) and lift coefficient as the objective functions and obtained their optimal solutions using a sequential quadratic programming (SQP) method, which is a gradient-based local optimization technique. At present, most of the WIG vehicles used a very large and high T-tail and sophisticated dynamic control to improve stability. However, these designs have high costs and lose the WIG vehicle's advantages. Designing effective and stable wings for WIG vehicles remains a great challenge. More recently, Park and Lee (2010) performed numerical optimization by considering the lift coefficient, static height stability, and lift-to-drag ratio as objectives. They obtained 64 non-dominated shapes for the airfoil section and compared the aerodynamic characteristics of the base model (NACA 0012) and Pareto individuals. Based on their analysis, they suggested that Pareto individuals with a flat lower wing surface have advantages in high lift and low drag when the wing is in ground effect. Static height stability consists of two aerodynamic centers: the aerodynamic center of height and center of pitch angle. The main and trailing wings play important roles affecting the aerodynamic centers of height and pitch, respectively. A numerical optimization was performed by Lee and Lee (2011) of a three-dimensional wing in ground effect. They classified four different groups of Pareto individuals according to the wing configuration: the drag and static height stability group, the static height stability group, the lift group, and the intermediate group. Interestingly, the intermediate group did not show any advantages over the other three groups, but it was still in the Pareto frontier and one of the global solution groups. They found that a slight forward sweep and a small chord ratio could improve the static height stability. However, both Park and Lee (2010) and Lee and Lee (2011) did not clearly investigate the relationship between wing shape and the static height stability. Considering only the

p ! Pi Re Sij ui , uj xi X ac Xh Xα α ε ρ μt vi

pressure control point of Bézier curve Reynolds number mean strain rates velocity component coordinate system aerodynamic center aerodynamic center of height aerodynamic center of pitch angle angle of attack or pitch angle turbulent dissipation rate density turbulent viscosity design variables

main wing for optimization, it is difficult to find a clear relationship between sectional shapes and static height stability. This study has two purposes: finding multiple optimal solutions under the cruise condition, and analysis of the aerodynamic characteristics and ground effect of the optimal set. In order to obtain stable and high-performance airfoils under the influence of the ground effect, we numerically performed shape optimization using computational fluid dynamics (CFD) and the multi-objective optimization technique. The lift coefficient, aerodynamic center of height, and lift-to-drag ratio, all of which significantly influence the performance of a WIG vehicle, are considered as the objective functions. The aerodynamic center of height is used for this optimization instead of the static height stability used in previous studies (Park and Lee, 2010; Lee and Lee, 2011). The airfoil shape is parameterized using Bézier curves, and their control points are used as the design variables. Due to tradeoffs between the conflicting objectives, the optimal solutions are a set of individuals (i.e., designs) that are not dominated by the other individuals within the design space. The 113 non-dominated optimal solutions, known as the Pareto optima (or sets), can be obtained by using a multi-objective genetic algorithm (MOGA). The Pareto optima include various shapes of the airfoil, and the aerodynamic characteristics of these Pareto optima have been analyzed in detail.

2. Physical model 2.1. Flows around airfoil A schematic configuration and the coordinate system of an airfoil in the WIG vehicle considered in this study are shown in Fig. 1. The flow around the airfoil is assumed to be two-dimensional, turbulent, and steady state with incompressible fluid. The turbulent flow of the air can be described by the Reynolds-averaged Navier–Stokes (RANS)

Fig. 1. Physical model and coordinate system of airfoil.

S.-H. Lee, J. Lee / Ocean Engineering 68 (2013) 1–13

3

equations, which can be expressed in tensor notation for mass and momentum as follows: ∂ ðρuj Þ ¼ 0 ∂xj

ð1Þ

∂ ∂p ðρuj ui −τij Þ ¼ − þ si ∂xj ∂xi

ð2Þ

2 ∂u τij ¼ μt Sij − μt k δij −ρu′i u′j 3 ∂xk

ð3Þ

where xj , j¼ 1,2 are the Cartesian coordinate vectors; uj represents the mean velocity components; −ρu′j u′i is the Reynolds stress tensor; and μt and Sij are the turbulent viscosity and the modulus of the mean strain rate tensor, respectively. These variables are defined as
 2 ρk ∂ui ∂uj μt ¼ f μ ⋅C μ ; Sij ¼ þ ð4Þ ε ∂xj ∂xi In this study, the RNG k−ε model proposed by Yakhot et al. (1992) is applied to turbulent flow around the airfoil. It is known that the RNG k−ε model with an additional term in its ε-equation can significantly improve the accuracy of turbulent flow (Yakhot et al., 1992). In this study the flow is assumed to be fully turbulent, and the boundary layer transition on the airfoil surface is ignored. 2.2. Flow analysis Two-dimensional incompressible Navier–Stokes equations are solved in the present computational analysis. To avoid the blockage effects and disturbances of outer boundaries, the computational domains are extended as follows: the inlet region is separated from the leading edge by 10c and the exit in the downstream direction is located at a distance of 20c from the trailing edge. In addition, the top outer boundary in the y-direction is extended by 10c from the upper crest location. To calculate the aerodynamic characteristics such as the lift, the lift-to-drag ratio, and the aerodynamic center of height, which are the objective functions for shape optimization, we used a commercial CFD package, STAR-CD (2001). Most of the computational cost for optimization is due to the evaluation of the objective functions. Therefore, a cost-effective and high-fidelity RANS solver for turbulent flow, such as STAR-CD S/W, needs to be used. For the treatment of the convective and diffusion terms in the governing equations, we employ the second-order upwind and central differencing schemes, respectively. The pressure–velocity coupling is resolved using the SIMPLE algorithm (Patankar, 1980). The accuracy and efficiency of the computational results are also dependent on how the grid system is constructed. For accurate and effective computation, this study used H-type grids with nonlinear distribution (refined next to the airfoil) as shown in Fig. 2. The computation was carried out for a Reynolds number of 2:8  106 based on the airfoil chord (c), and a corresponding inlet velocity of 34.28 m/s was applied. No-slip boundary conditions on the airfoil surface, a moving wall boundary at the ground, and a pressure boundary condition at the exit (x-direction) were employed. The slip-wall boundary condition was prescribed at the top-outer boundary. The solutions were considered to have converged when all the residuals for the continuity and momentum equations are less than or equal to 10−6. An AMD AthlonTM processor (2.2 GHz) with Linux was used for the computations. One evaluation required about 3 min to obtain converged solutions for the selected design variables. For the calculation of the aerodynamic center of height, which was one of the objectives, derivatives with respect to height were required; thus, the CFD calculation had to run three times to evaluate one design. The overall computation time to complete one optimization

Fig. 2. Non-uniformly distributed grid for two-dimensional wing in ground effect. (a) Computational grid and extended computational domain. (b) Close-up of twodimensional grid around wing.

process was about 180 h (about 8 days). To minimize the computation time, previous computational results were stored, so that when a new evaluation was required, the optimizer searched the storage first and then began a new evaluation.

3. Airfoil shape optimization Many optimization technologies have been developed to obtain optimal solutions for various industrial applications. Among them, gradient-based methods, such as SQP and sequential linear programming (SLP), are well-known techniques that seek to find the optimum using local gradient information so that they obtain a local optimum. In particular, to treat multi-objective optimization problems, gradient-based methods convert a number of the objectives into a single objective using a weighted sum of all the objective functions with a normalization technique. However, this approach may distort the original design space and lead to a solution of local optima. In addition, when the orders of magnitude among the objectives are significantly different, gradientbased optimizers tend to focus on the largest objective function only. Normalizing the objectives can resolve the multi-objective problems, but the distortion of the design space may still remain. To avoid these disadvantages of using the gradient-based methods, a multi-objective genetic algorithm (MOGA) is introduced to find the global optima; this is an attractive approach for airfoil shape optimization with multiple objectives. The optimal solutions of the multi-objective problem cannot be defined uniquely except for a linearly dependent system. All designs can be classified into dominated and non-dominated individuals. The latter is called a Pareto set, which is located on the front line of the design space. 3.1. Concept of Pareto frontier (or set) A multi-objective optimization does not look for a unique solution, but rather seeks a set of potential solutions (i.e., Pareto set).

4

S.-H. Lee, J. Lee / Ocean Engineering 68 (2013) 1–13

Thus, all such solutions belong to a hyper-surface: a locus of optimal solutions that is the trade-off surface between conflicting objective functions, called the Pareto frontier. From the Pareto frontier point of view, none of the optima are dominated. In other words, none of the objectives can be improved without worsening at least one of the other objectives. To illustrate the Pareto frontier, we consider a multiobjective problem with m objective functions to be minimized. Individual A dominates with respect to individual B if, for at least one of the objectives, the values of the objective functions for A are strictly better than those of B and if, for all other objectives, A is not worse than B. According to this definition, an individual will be considered Pareto optimal if it is non-dominated. Mathematically, B is Pareto optimal if no individual A exists such that A 4 B⇔½∀i∈f1; …mgF i ðAÞ ≤F i ðBÞ∩½∃i∈f1; …mgF i ðAÞ o F i ðBÞ p

ð5Þ

3.2. Multi-objective genetic algorithm (MOGA) The genetic algorithm (GA) is a unique global optimum algorithm based on “survival of the fittest,” or mechanism of natural selection and reproduction (Goldberg, 1989; Booker, 1994; Poloni et al., 1999). In this algorithm, each individual has a number of strings representing the characteristics of a particular design, similar to a gene in natural life. Each generation consists of a group of individuals. To evolve a new solution, new individuals for the next generation are obtained by using genetic operations such as selection, crossover, mutation, and niche. The selection process governs the direction of the evolution toward survival of the fittest and determines parents from the genetic pool. Crossover produces offspring from the superior individuals to evolve the generation. However, since crossover and selection are deterministic processes that only rearrange the genes, they do not enable introducing completely new characteristics. This deficiency of searching the design space causes the solution to converge to a local optimum. Mutations, however, can prevent this deficiency, and balanced searching between crossovers and mutations can effectively lead the GA to the global optima. In contrast to a random search, these genetic operations continually improve the average fitness of a generation. Note that the properties of adjacent individuals in the design space are similar to each other. To extend exploration, it is necessary to control the number of individuals inside the niche radius. Before sufficient exploration, at an early stage of the evolution, the existence of local optima is very attractive. The offspring tends to gather around the local optimum points, and this leads to early mature convergence. However, appropriate selection of the niche radius is able to control and prevent a GA from early mature convergence. In this study, the binary distance between two individuals – instead of the n-dimensional norm – is employed as the niche distance, and it can be written as follows:     L L L h mi r ij jdi −dj j m⋅Δx ¼ ∑ ¼ ∑ ¼ ∑ ð6Þ n⋅Δx R R k¼1 k¼1 k¼1 n k k k where jdi −dj jk is the distance between i and j individuals in the k variable, m is the binary distance, and n is the niche binary distance. It is difficult to present all niches simultaneously in the real design space with a single niche radius, so a niche radius for each design variable is required. For the binary distance, however, one niche radius is sufficient. In GA processes, a set of the solutions (population) evolves generation by generation. After the random first generation, the next generation can be obtained by the genetic operations of selection, crossover, and mutation. Any new individual with an identical twin in the new generation is ignored in order to maintain the diversity of the system. In this study, due to the multi-objectives, tournament selection is used. Pairs of individuals

are randomly selected from the population. The one with highest fitness is copied into the mating pool. We then employ tournament competitions n times, using the value of n to manage the selection pressure. The steps constituting the selection of individuals for the mating pool are as follows: (1) Is there a dominant individual between the selected two individuals? (2) Is there a Pareto individual in a pair of individuals in the same generation? (3) If both individuals are Pareto individuals, select the less crowded individual. (4) The superior individual of the two is selected. These processes are close to the mathematical definition of a Pareto optimum. Based on these operations, the authors of this study developed SMOGA (simple multi-objective genetic algorithm). Details can be found in the references (Lee and Lee, 2011; Lee et al., 2006).

4. Objective function and validation To evaluate the objectives of the individual designs such as the aerodynamic forces and aerodynamic center of height, a CFD analysis is automatically performed under the optimizer. A series of pre- and post-processes such as grid generation, CFD model setting, solving, and post-process are performed by a few batch scripts and small executable codes developed by the authors. A SMOGA is used as an optimizer so that normalization and weighting for the multi-objective functions are not required. Optimization is performed under normal cruise conditions such as h/c¼0.25 (non-dimensional height) and α ¼21 (angle of attack) at Re ¼ 2:8  106 . 4.1. Mathematical formula for optimization Optimization involves finding the best or potential value of the design variables that minimize and/or maximize the objective functions while satisfying the constraints and bounds. Optimization problems comprise design variable(s), objective function(s), constraint(s), and side constraint(s). 4.1.1. Mathematical formula for multi-objective optimization The ground effect can simultaneously improve the lift and drag of a typical airfoil when the airfoil is in the vicinity of the ground. However, these favorable phenomena can give rise to a stability problem for a WIG vehicle. Therefore, in this study, the lift coefficient, F1(X), the aerodynamic center of height, F2(X), and the lift-to-drag ratio, F3(X), are considered to be the objective functions, simultaneously. The static height stability condition consists of two aerodynamic centers, height (X h ) and pitch (X α ), given by H:S: ¼ X h −X a ¼

C m;h C m;a − C l;h C l;a

ð7Þ

The primary function of a horizontal tail is to provide longitudinal control and trim when the vehicle cruises steadily. The size of the horizontal tail must be sufficient to provide the necessary stability and control of the WIG vehicle; thus, a large T-tail is generally used for sufficient stability. The moment at the center of gravity (CG) is proportional to the value of the horizontal tail volume ratio that consists of the horizontal distance between the CG of the vehicle, the aerodynamic center of the horizontal tail, and the platform area of the horizontal tail. The horizontal tail is placed out of ground effect and has only a minor effect on X h .

S.-H. Lee, J. Lee / Ocean Engineering 68 (2013) 1–13

However, the main wing, which is in ground effect, is important to X h . As a result, the main wing and the horizontal tail are important factors for X h and X α , respectively. In this study, considering only a main wing, we used the X h as an objective function instead of H:S:, including both X h and X α . The value of C m is measured at the aerodynamic center (AC). Note that, if C m;h is zero, then X h is at the same location as the AC, which is set as a quarter chord in this study. The strict stability condition (Kornev and Matveev, 2003) of X h oX cg o X α or X h ∼X cg o X α can easily be satisfied when a wing with X h  0 or X h  X ac is installed on a WIG vehicle with a small horizontal tail. The mathematical formulation of optimization is defined as follows: Find the control points; X ¼ fx1 ; x2 ; …; x18 gT

ð8Þ

to maximize the lif t; F 1 ðXÞ ¼ C l

ð9aÞ

C   m;h  to minimize the aerodynamic center of height; F 2 ðXÞ ¼ jX h j ¼   C l;h ð9bÞ to maximize the lif t
to
drag ratio; F 3 ðXÞ ¼

Cl Cd

subject to the bounds; X Li ≤X i ≤X Ui ; for i ¼ 1; 18

ð9cÞ ð10Þ

4.1.2. Design variables and constraints A wide variety of methods to represent the airfoil shape have been used in the shape optimization of airfoils. One of the most popular methods is the Bézier curve, using control points to define airfoil shapes. The Bézier curve generally offers the following advantages: (1) there is direct control of the coordinates of control points, which provides an easy way to control airfoil profiles with a small number of points; and (2) the curve always passes through the first and last control points and lies within the convex hull of the control points which leads to smooth changes of airfoil shapes. In the present work, the airfoil is parameterized using four Bézier polynomials: two fourth-order Bézier curves for the leading edge and two third-order Bézier curves for the trailing edge. A Bézier curve can be represented as follows: k ! ! P ðtÞ ¼ ∑ P i Bi;k ðtÞ;

t∈½0; 1

ð11Þ

i¼0

k! t i ð1−tÞk−i ð12Þ i!ðk−iÞ! ! where P i (t) are the coordinates of control points, the parameter t is bounded from 0 to 1, and Bi;k is a blending function. Fig. 3 shows a schematic diagram of a WIG airfoil using a Bézier-curve. As shown in Fig. 3, the airfoil shape can be completed using six control points on each side of the airfoil. Except for continuous or fixed points on the airfoil, such as the maximum thickness (v5), the smooth leading edge (v1, v10) and the trailing edge, the number of design variables is 18. The three control points near the maximum thickness have the same height, which ensures a smooth airfoil profile. The only explicit constraints are the upper and lower boundaries of the design variables. Generally, the choices of range and resolution in design variables are the most important factors in the optimization processes. Therefore, the pre-optimizations are performed four times in order to properly decide the upper and lower limits of the design variables. Because there is insufficient information for the design space after first pre-optimization, the upper and lower bounds are extended arbitrarily until the intersection of control points is avoided. To save computational cost in this study, eight-bit resolution is used for every design variable, and the design space reaches to Bi;k ðtÞ ¼

5

Fig. 3. Airfoil geometry parameterization.

25618. In order to calculate cost effectively, only 20 evolutions and 25 individuals are used for preliminary optimization. Once the Pareto set is obtained from the preliminary optimization, the boundaries of the design variables are modified. Note that when optima are concentrated in the vicinity of one point, the range of the design variables is reduced. In contrast, the range of the design variables increases when the optima are placed around the bounds. Using this methodology, relatively high resolution can be obtained with a relatively low number of bits. After the final pre-optimization, 10-bit resolution is used and the design space becomes 102418 for the optimization. This high-resolution binary representation of the properly defined upper and lower bounds can almost resolve the real values of the design variables. All the bounds of design variables are listed in Table 1. In Table 1, v2–v4 and v11–v13 represent the leading edge profiles, v5 and v14 are the maximum thicknesses, and v6 and v15 are the locations of maximum thickness. 4.1.3. Optimization parameters of MOGA Table 2 provides the optimization parameters used in this study. Forty individuals for one population are used and the selection pressure value is chosen to be 2 to enhance the convergence rate. In each tournament, two candidates are randomly selected from the current generation, and through two competitions, the winner has a chance to become a parent for the reproduction. The number of cutting lines for the crossover operation that exchanges genes is important. In this study, two cutting lines are used to maximize the life of a schema, which is a useful pattern in the gene. Mutation is the occasional random alteration of the value of a string position with a small probability. By itself, mutation is a random walk through the string space. When the value of a mutation is as high as a few percent, the GA cannot converge to proper solutions with evolutions, and the search may become completely random. To prevent this and maintain the balance between exploitation and exploration, a 0.5% mutation rate is chosen. When a new offspring individual is found to be a genetic twin of any other individuals in the next generation, the individual is ignored, and a different individual is generated. 4.2. CFD model validation To validate the CFD models and evaluation processes, lift coefficients (Cl) of the NACA0015 airfoil are compared with the experimental results for two different conditions: in and out of ground effect. Fig. 4 shows Cl according to the non-dimensional height (h/c), considering the ground effect at α ¼ 01 and Re¼2:8  106 . The values of Cl from the CFD model are also compared with the experimental results measured by Ranzenbach and Barlow (1994). Their experiment was conducted with a fixed ground (i.e., a no-slip wall). The hollow circle (−∘−) in Fig. 4 represents the results under the same conditions as their experiment, and the hollow triangle represents the results using the moving-wall boundary condition.

6

S.-H. Lee, J. Lee / Ocean Engineering 68 (2013) 1–13

Table 1 Upper and lower limits of design variables. Design variables

1st

2nd

3rd

4th

Final

0.036 0.045 0.08 0.22 0.06 0.32 0.55 0.025 0.85 −0.055 −0.065 0.05 0.15 −0.07 0.28 0.45 −0.02 0.65

0.035 0.045 0.085 0.2 0.075 0.31 0.6 0.02 0.9 −0.06 −0.07 0.05 0.16 −0.05 0.25 0.4 −0.035 0.6

0.04 0.05 0.085 0.2 0.09 0.31 0.65 0.02 0.9 −0.06 −0.03 0.07 0.16 −0.035 0.25 0.35 −0.03 0.6

0.04 0.05 0.085 0.2 0.09 0.31 0.65 0.02 0.9 −0.06 −0.03 0.07 0.16 −0.035 0.25 0.35 −0.03 0.6

0.048 0.065 0.09 0.26 0.08 0.35 0.65 0.035 0.95 −0.035 −0.05 0.14 0.2 −0.045 0.32 0.55 −0.01 0.75

0.05 0.065 0.09 0.24 0.09 0.33 0.7 0.04 0.96 −0.04 −0.05 0.14 0.19 −0.04 0.3 0.5 −0.005 0.8

0.05 0.065 0.09 0.24 0.1 0.34 0.7 0.04 0.98 −0.038 −0.015 0.1 0.2 −0.015 0.3 0.45 −0.002 0.75

0.05 0.065 0.09 0.24 0.01 0.034 0.07 0.04 0.98 −0.38 −0.015 0.1 0.2 −0.015 0.3 0.45 −0.02 0.75

Lower boundaries vari01 vari02 vari03 vari04 vari05 vari06 vari07 vari08 vari09 vari10 vari11 vari12 vari13 vari14 vari15 vari16 vari17 vari18

0.025 0.04 0.05 0.15 0.05 0.27 0.5 0.01 0.7 −0.045 −0.06 0.05 0.15 −0.07 0.27 0.5 −0.03 0.7 Upper boundaries

vari01 vari02 vari03 vari04 vari05 vari06 vari07 vari08 vari09 vari10 vari11 vari12 vari13 vari14 vari15 vari16 vari17 vari18

0.045 0.06 0.13 0.25 0.07 0.33 0.6 0.03 0.9 −0.025 −0.04 0.13 0.25 −0.05 0.33 0.6 −0.01 0.9

The lift coefficients obtained in this study are in good agreement with the experimental data. In addition, as the NACA0015 airfoil approaches the ground, the oncoming air accelerates at a leading edge and then decelerates after a maximum thickness (x/c¼0.3) because of a convergence–divergence passage. Thus, the increased velocity on the lower surface produces negative lift, a phenomenon called the Venturi effect. Therefore, a symmetric airfoil without a camber, such as NACA0015, is not appropriate for a WIG vehicle. The sudden increase in lift in the experiments was mainly a consequence of the flow separation in front of a narrow passage between the airfoil and ground, which does not occur in a real flight. From the comparison, it is clear that the present CFD model is capable of predicting the characteristics of the flow past the airfoil. On the other hand, comparisons of C l with respect to the angles of attack for three Reynolds numbers (Re¼3:31  105 , 1:27  106 , and 3:26  106 are calculated but not shown in this paper) show good agreement with the experimental results. The detailed results can be found in the authors’ previous study (Park and Lee, 2010).

5. Results and discussion 5.1. Optimal solutions To obtain optimal performance of the WIG vehicle, it is evident that both Cl and Cl/Cd must be maximized while X h is simultaneously minimized. However, for a multi-objective optimization,

the optimal solution is not unique; rather, it is a set of nondominated solutions that encompass the trade-off among objectives. Figs. 5 and 6 show the online and offline performances according to the generation, to demonstrate the convergence histories for Cl and X h . DeJong (1975) devised these two performance measures (i.e., online and offline) to quantitatively evaluate the performance of GAs. They are defined as follows: on
line

fe

ðiÞ ¼

of f
line

fe

ðiÞ ¼

1 i ∑ f ðjÞ i j¼1 e

ð13Þ

1 i n ∑ f ðjÞ i j¼1 e

ð14Þ

where f e ðjÞ is the value of the objective functions at generation j for n environment e, and f e ðjÞ represents the best values of the objective functions until a given generation i, for j ¼ 1; 2; …; i. The online on
line performance, f e ðiÞ, is an average fitness of all trials up to of f
line the current generation and the offline performance,f e ðiÞ, is a running average of the best individuals up to a particular generation. In addition, the online performance is used to measure the developing performance, whereas the offline performance was originally devised to gauge the convergence of the optimization process for a single-objective optimization problem. As shown in Figs. 5 and 6, the objective functions gradually converge as the generation is advanced. It can be seen that a moderate convergence is achieved for all the objective functions after 24 generations.

S.-H. Lee, J. Lee / Ocean Engineering 68 (2013) 1–13

7

Table 2 Parameters for genetic algorithm. Parameters

Values

Population Generation Crossover rate Mutation rate Tournament level Niche

40 30 0.8 0.5% 2 1

Fig. 5. Online and offline performance for convergence history of C l .

Fig. 4. Lift coefficient according to height (h/c) for slip and no-slip conditions (filled: experiments, Ranzenbach and Barlow, 1994; hollow: present study). Fig. 6. Online and offline performance for convergence history of X h .

In this study, the minimized or maximized objective functions are Cl, Cl/Cd, and X h ; and, these optimal solutions should be displayed in three-dimensional space (i.e., Cl, Cl/Cd, and X h ). However, it is hard to explain the relationship and features of each function in one threedimensional figure, so all individuals have been projected onto a twodimensional design space, as shown in Figs. 7–10. In the figures, the hollow squares represent the Pareto individuals and the filled squares represent the dominated solutions. The simple multi-objective genetic algorithm (SMOGA), which is based on the GA and developed by the authors, is different from a random search but cannot obtain exactly the same Pareto set each time, because of its random works in the mutation and selection operations. To confirm the computational feasibility of SMOGA, the optimization is performed three times in a row. The three Pareto sets are plotted in Fig. 7. Approximately 100 potential solutions (Pareto individuals) are obtained from each optimization. When all Pareto individuals are plotted in the same figure, each set of Pareto individuals cannot be different from the others, so only one-fourth of the individuals are plotted in Fig. 7. The three lines in Fig. 7 are linear regressions of each Pareto set. Every line has a slight deviation, but they show a similar tendency and each Pareto set is placed in a small band—the Pareto frontier. The three bands overlap each other, and SMOGA can find the Pareto set properly despite the random works in operations. In Figs. 8–10, some of the Pareto optima can be observed in the area behind the frontier line from which the dominated individuals exit. This inconsistency comes from the projection of the three-dimensional design space onto the two-dimensional space. The comparison for the Pareto set was performed between entire sets of individuals at the end of the optimization instead of during the last generation. Consequently, 113 Pareto individuals are selected from 1200 individuals (30 evolutions 40 individuals). In Fig. 8, Pareto individuals are numbered according to their lift coefficient from p001 (the lowest lift coefficient) to p113 (the highest lift coefficient). In order to observe the relationship between the design variables and the objectives, seven individuals

are randomly listed among the Pareto frontier, as shown in Fig. 8. A multi-objective optimization does not seek a unique solution but rather a set of potential individuals (Pareto optima). The trade-off between objectives (C l and X h ) can be observed in Fig. 8. When the lift is improved, X h is degraded necessarily, and vice versa. The actual values of X h are negative, and thus the locations of the aerodynamic center of height are located downstream of a quarter chord (the aerodynamic center). If the stability margin (or static height stability, X h −X α ) becomes lesser than −0.15, the static stability becomes excessive, which leads to dynamic instability. If the stability margin is negative but larger than −0.05, then there is insufficient stability, which leads to a long-period (phugoid) instability (Kornev and Matveev, 2003). Therefore, to cruise continuously under the ground effect, the center of gravity should be located between the aerodynamic centers of height and pitch and close to the center of height:X h 4 X CG 4 X α , X CG  X h . Generally, since the center of gravity of a vehicle is placed slightly behind the AC, the p001 individual whose X h is placed near the AC can easily satisfy the stability condition: −0.05 4H:S: 4−0.15. However, C l and X h conflict with one another, as shown in Fig. 8; p001 with high stability (small X h ) shows the lowest value of C l . This value indicates that the p001 individual has the lowest endurance among the Pareto set. On the other hand, the p113 individual, which has the largest C l among the Pareto set, is difficult to design in a WIG vehicle because its aerodynamic center of pitch angle is downstream of X h . This situation implies that the WIG vehicle should have a large horizontal tail to move X α further downstream, and a complex control system to satisfy stability when it cruises in ground effect. Fig. 9 shows the Pareto set and dominated individuals with C l =C d and C l . Since C l is proportional to C l =C d , there exists a unique solution (single optimum) instead of a Pareto set, such as p113, which is located at the extreme top right and has the largest C l , as well as C l =C d in the given environment. However, Pareto individuals are shown in the area behind the p113

8

S.-H. Lee, J. Lee / Ocean Engineering 68 (2013) 1–13

Fig. 7. Comparison of Pareto optima: line—first-order regression line for Pareto optima, symbol—Pareto optima.

Fig. 8. Pareto set and dominated individuals with respect to X h and C l .

Fig. 9. Pareto set and dominated individuals with respect to C l =C d and C l .

individual and in the other projected objective plots in Figs. 8 and 10 (C l vs X h ), and p113 is one of the Pareto optima because of the trade-offs between the other objectives: C l vs X h and C l =C d vs X h . Note that two objectives, C l and C l =C d , are linearly dependent, and thus, an individual high C l generally shows high C l =C d . However, a few individuals have high C l but low C l =C d because of high C d . The rank of Pareto individuals can be changed slightly according to the value of C d . Therefore, C l and C l =C d should be considered as objectives for the optimization. In Fig. 10, C l =C d is maximizing, X h is minimizing, and Pareto optima are placed along the Pareto frontier line at the bottom right side. Since SMOGA is based on GAs, the individuals are randomly and uniformly distributed across the objective space. According to the evolution, the individuals move gradually toward a Pareto

Fig. 10. Pareto set and dominated individuals with respect X h to C l =C d .

frontier, and thus, many individuals – including Pareto optima – are placed next to the Pareto frontier line, as shown in Fig. 10. By projecting the Pareto set from the three-dimensional objective space to the two-dimensional space, a Pareto set can also be observed behind the frontier line, as shown in Fig. 8. SMOGA, the optimizer used in this study, can find these potential solutions without weighting functions, which highlights the capabilities of the optimizer. Consequently, we understand that the two objective functions, C l =C d and X h , are inversely related to each other so that improving one necessarily deteriorates the other. The trade-off between C l =C d and X h can be observed in Fig. 10, showing that there are no individuals that improve the values of X h and C l =C d , simultaneously. To apply the Pareto set obtained in this study to designing a WIG effect vehicle, a designer decides the acceptable range of X h and then selects candidates from the Pareto individuals that satisfy the requirements of the lift or lift-to-drag ratio. Through a series of choices, a designer can achieve the optimal and global design under a given environments. Table 3 shows the correlation and probability value results from the correlation analysis comparing the design variables and objective functions. In the table, p-values less than 0.05 imply that the p-coefficient (correlation) is reliable. Note that there are a direct correlation between the design variables and objectives, except for v07, v08, v12, and v13, which are not listed in Table 3. The six design variables (in bold) that show the three largest values of the p-coefficient at each surface are listed. The design variables of the leading edge and the maximum thickness, which guide the oncoming air into the gap between the airfoil and the ground, are very important for the aerodynamics and optimization. The variables v05, v06, v14, and v15 represent the value of the maximum thickness, while v04 and v11 indicate the leading edges of the upper and lower surfaces, respectively, as shown in Fig. 3. In particular, v05 and v14, which represent the thickness of the airfoil, have the largest values for C l and C l =C d on each side except v05 for C l =C d . It is well known that thickness has a strong effect on both C l and C l =C d . However, increased thickness, especially v05, improves C l and tends to increase C d . Thus, the values are relatively small in C l =C d as shown in Table 3. However, increasing v05 has a negative effect on X h . It also highlights a strong conflict between lift and stability. The lower surface directly interferes with the ground surface; thus, it significantly influences the lift and the stability. Because the normal cruise condition at a low angle of attack (i.e., a twodimensional flow) and a constant cruise velocity are considered in this study, there is a slight flow separation. Accordingly, the drag is mainly a result of friction with the wetted surface; lift comes from stagnated air flow on the lower wing surface; and stability with respect to the AC comes from pressure distributions on the wing surface. Therefore, as shown in Table 3, v14 and v17 have the

S.-H. Lee, J. Lee / Ocean Engineering 68 (2013) 1–13

Table 3 Correlations and p-value between variables and objectives. Variable

Cl

Pareto

p coef

p-value

p coef

p-value

p coef

p-value

v01 v02 v04 v05 v06 v09 v10 v11 v14 v17

0.347 0.291 −0.321 0.406 0.121 0.189 0.201 0.360 0.644 0.414

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

−0.078 0.095 0.077 −0.128 −0.134 −0.137 −0.157 −0.127 0.285 0.512

0.007 0.001 0.008 0.000 0.000 0.000 0.000 0.000 0.000 0.000

0.248 0.244 −0.275 0.236

0.000 0.000 0.000 0.000

0.150 0.291 0.673 0.440

0.000 0.000 0.000 0.000

Xh

Cl/Cd

largest p-coefficient values among the variables. These values imply that C l and X h are sensitive to variations in v14 and v17, respectively. To maximize C l and minimize C d , the oncoming air should be stagnated on the lower wing surface, causing the dynamic pressure to change to static pressure. A lower wing surface that is parallel to the ground surface (e.g. following Carter, 1961; Fink and Lastinger, 1961; Park and Lee, 2010) or a thinner shape next to the trailing edge (such as NACA4412) are favorable designs for C l and C l =C d . The stagnation air under the lower wing surface leads to high lift and the low drag at the same time; yet, these airfoils tend to increase in pressure next to the trailing edge because of the stagnation (Park and Lee, 2010). This shape results in a high deviation of C m with respect to height, which leads to an unfavorable condition for the stability (X h ). On the other hand, v17 has the largest value among the p-coefficients for X h and a clear correlation with X h . Since X h is minimized, an increase in v17 leads to an increase in X h and results in an unfavorable stability condition. Recall that v17 represents the thickness of the lower wing surface next to the trailing edge, and its sign is negative (−). A larger absolute value of v17 reduces the gap between the lower wing surface and the ground and vice versa. The thick airfoil next to the trailing edge generates the Venturi effect locally and reduces the lift slightly. Nevertheless, the pressure next to the trailing edge decreases according to the Venturi effect, and C m with respect to the AC reduces so that the deviation of C m against height is reduced, which improves the stability. The highest stable airfoil, such as p001, with a small v17 is similar to an S-shaped airfoil (Chun and Chang, 2002; Rozhdestvensky, 2006), which is generally known as a high stability airfoil.

9

ground and the airfoil. At the same time, the Venturi effect will increase at the location of the convergence–divergence passage. Consequently, the derivative of C m with respect to the aerodynamic center will be constant when the height is changed. Thus, this convergence–divergence passage might sacrifice lift slightly, but it will improve X h . Since the center of the pitching moment is located at a quarter chord (i.e., at the aerodynamic center), X h is located at the same quarter chord when the value of X h is zero. However, for optimization purposes, the absolute value of X h is employed in this study. As shown in Table 4, all the values of X h are negative, and thus, the location of X h should be placed the downstream of the AC. The p001 individual, which has the smallest C l among the Pareto optima shows the smallest value of X h , which is favorable for stability. In contrast, Pareto individuals with larger values of C l have X h located further downstream. The X h of p001 is placed at the AC under the given environment, whereas the X h of p113 is 0.366. The actual location of X h of p113 is placed at x=c ¼ 0.616 as measured from the leading edge; this is slightly downstream from the center of a chord. In this case, the strict stability condition X h o X cg o X α cannot be easily satisfied without a large horizontal trailing wing. Consequently, p001 is more favorable for designing a WIG vehicle that is aerodynamically and structurally stable with a minor trade-off in C l . As shown in Table 4, C m;h has a more dominant effect on X h than C l;h does. The reason for the negative value of X h is that C l decreases as the height increases and thus, C l;h in the numerator becomes negative. This is one of the key features of a groundeffect vehicle. On the other hand, the reason for a positive value of C m;h is that the C p on the trailing edge with a large moment arm increases with respect to the decreasing height, which increases the negative pitching moment at the aerodynamic center. However, p001 has a convex lower wing surface next to the trailing edge, which is similar to an S-shaped airfoil. The convex surface can minimize the effect of a negative value of C m with height and can improve stability. Since the optimization in this study is performed at a low angle of attack (α ¼ 21), no flow separation is observed. Accordingly, the drag mainly comes from the friction of the surface area. Note that the C d values of the Pareto individuals listed in Table 4 are all nearly identical. Consequently, it is clear that C l =C d is determined by C l , but still slightly affected by C d according to the shape of the airfoil. Fig. 12(a) and (b) shows the lift coefficient vs angle of attack and the lift coefficient vs height (h=c), respectively. The results of three Pareto individuals, p001, p052, and p113, are included.

5.2. Aerodynamic characteristics of Pareto individuals The aerodynamic analysis of the Pareto optima is performed under the same conditions as the optimization (i.e., the cruise conditions including h/c¼0.25, and/or α ¼ 21 at Re ¼ 2:8  106 ). In Fig. 11, the sectional profile and its pressure distributions along the surface are plotted. According to individual characteristics, different shapes on the upper and lower wing surface can be obtained by optimization. The straight lower surface can utilize the ground effect and minimize the Venturi effect, which enforces a negative lift force as shown in various studies (Mahon and Zhang, 2005; Zerihan and Zhang, 2000). Park and Lee (2010) and Carter (1961) also presented the effect of the straight lower surface in detail. As illustrated in Fig. 11, p113 (which has the highest lift among the Pareto individuals) shows a straight lower wing surface that improves lift, whereas p001 has a convex lower wing surface at x=c ¼ 0:8, which implies that a convergence–divergence passage exists next to the trailing edge. The pressure distribution on the lower surface will increase with decreasing distance between the

Fig. 11. Airfoil profiles and pressure distributions of three Pareto individuals: p001, p052 and p113 (h/c ¼ 0.25, α ¼ 21, and Re ¼ 2:8  106 ).

10

S.-H. Lee, J. Lee / Ocean Engineering 68 (2013) 1–13

The value of C l in ground effect increases nonlinearly with respect to the angles of attack, as shown in Fig. 12(a). The general features of the ground effect increasing lift as the gap between the ground and the lower wing surface decreases is shown clearly in Fig. 12(b), except for α o 21. However, in the case of α ¼ 01, the value of C l at h=x ¼0.15 is smaller than that of C l at h=c ¼ 0.4. This difference is mainly due to the optimized shape of the lower wing surface, which is nearly parallel to the ground at α ¼21. At α ¼01, the value of C l decreases because of the Venturi effect via the convergence– divergence passage. In this condition (α ¼01), it is possible that the vehicle could suddenly lose lift or become out-of-trim with a small disturbance. This lift discrepancy is large at lower heights, such as h=c ¼ 0.15, as shown in Fig. 12(a). Comparing lift coefficients between h=c ¼0.4 and 0.15 for α ¼01, the lift reduction of p113 is about 5.6%, while the lift reduction of p001 is about 10%. These values imply that an S-shaped wing surface can improve stability when in ground effect, but there is a high probability that this shape could lose lift when its ground condition is not proper. To compare the aerodynamic characteristics of p001 (high stability), p052, and p113 (high lift), Fig. 13 depicts a drag polar. Aerodynamically, a favorable airfoil should show both high lift and low drag. The ratio of lift to drag, C l =C d , is a good measure of aerodynamic efficiency (Anderson, 1998). Lift is required to sustain the vehicle constantly in the air, but lift-to-drag is an important factor for cruising efficiently and reducing cruise cost. As shown in Fig. 13, lift and drag are both improved when the airfoil is in ground effect, as demonstrated by Wieselsberger (1922) using image vortex methods; however, drag reduction occurs in this study for a slightly different reason. Wieselsberger's research used three-dimensional analysis (image vortex methods), and a key factor for his drag reduction was a decrease in the induced drag in ground effect. In contrast, this study uses a two-dimensional analysis, so the main drag reduction comes from a decrease in the dynamic pressure on the lower surface. Thus, the increase in dynamic pressure eliminates the drag reduction caused by the ground effect, as shown in the circle in Fig. 13. Moreover, the

Venturi effect increases both the dynamic pressure and the drag at the lower wing surface. It is expected that p052 has a higher lift coefficient and a higher lift-to-drag ratio compared to p001, but it has high lift as well as the high drag. Therefore, p052 shows similar aerodynamic characteristics to p001, and p052 does not improve C l =C d as much as expected. Recall that X h is composed of C m;h (the numerator) and C l;h (the denominator). Each term of X h (i.e., C m and C l ) is closely investigated in Fig. 14, which compares the C m of p001 (stable), p052, and p113 (lift) with respect to heights and angles of attack. For this figure, the non-dimensional height between the wing and the ground is measured at the trailing edge, the center of the pitching moment is a quarter chord (i.e., the AC), and the nose-up pitching moment is positive. The moment coefficients of three individuals are negative, and the pitching moments are nose-down. In Fig. 14, the p113 individual shows the pitching moment's largest value, followed by p052 and p001. Variations of the moment with decreasing height are large at high angles of attack. Therefore, the largest variation of the pitching moment can be observed at α ¼ 101. These variations require the location of X h to be further behind the AC so that the CG can be placed in the strictly narrow range to satisfy the H:S: condition: X h ≤X cg o X α . Consequently, this makes a WIG vehicle difficult to operate in ground effect. On the other hand, the overall slopes of the pitching moment coefficients with respect to height are positive except at α ¼ 01, and the negative pitching moment becomes stronger as height decreases. Even though the forces (C p ) around the leading edge are large, the effect on the pitching moment is not significant because of its short moment arm. In contrast, the force next to the trailing edge is relatively small, but the effect on the pitching moment is large because of its large moment arm. Generally, when the airfoil is in ground effect, oncoming air is stagnated on the lower wing surface. Thus, the pressure distribution on the lower wing surface tends to be constant (Ahmed et al., 2007; Hsiun and Chen, 1996) or slightly increased at the trailing edge (Park and Lee, 2010). Therefore, the negative pitching moment is increased as the height

Table 4 Comparison of objective functions amon Pareto individuals.

p001 p014 p029 p052 p093 P104 p113

abs (X h )

Xh

C m;h

C l;h

Cm

Cl

Cd

C l =C d

0 0.005809 0.025581 0.037183 0.069385 0.087285 0.125957

0 −0.005809 −0.025581 −0.037183 −0.069385 −0.087285 −0.125957

0 0.001121 0.004552 0.007673 0.015453 0.018847 0.031094

−0.21168 −0.193024 −0.177956 −0.206369 −0.222707 −0.215924 −0.246865

−0.082265 −0.091658 −0.095846 −0.099626 −0.101641 −0.108605 −0.111067

0.647526 0.673585 0.682074 0.693821 0.713493 0.722933 0.740048

0.010559 0.010671 0.010652 0.010764 0.010538 0.010562 0.010257

61.322483 63.120823 64.029587 64.457047 67.708511 68.448288 72.152702

Fig. 12. Comparison of lift coefficients from three Pareto individuals: (a) C l vs. α and (b) C l vs. h=c.

S.-H. Lee, J. Lee / Ocean Engineering 68 (2013) 1–13

decreases as shown in Fig. 14. For α ¼01, the reverse tendency – an increase in the moment with a decrease in the height – is observed. The pitching moment increases because the airfoil shape is optimized at α ¼ 21. At a lower angle of attack, such as α ¼01, the oncoming air is insufficiently stagnated on the lower wing surface. In particular, both p001 and p052 have an S-shaped lower surface next to the trailing edge—thus, they have Venturi effects at the same location and show larger variations of pitching moment with height changes than does p113. As a result, the pitching moment for p001 and p052 is constant at α ¼21, whereas the pitching moment for p113 is constant at α ¼ 01, as shown in Fig. 14. Fig. 15 shows the pitching moment variations with respect to height for three Pareto individuals. As expected, p001 and p052 exhibit similar tendencies and are somewhat smaller than p113, which is a favorable individual for C l . On the other hand, for p113, which has a straight lower wing surface, C m;h increases for all angles of attack as height decreases, whereas for p001 and p052, C m;h increase at α≥41. As shown in Fig. 14, the straight lower wing surface continues to increase the pitching moments via the stagnated air under the lower wing surface. It is interesting to note that the variations of C m;h for the three Pareto individuals is similar at a high angles of attack, such as α ¼101. This similarity stems from the decreased Venturi effect (S-shape) at the trailing edge. Additionally, the pressure distributions under the lower wing surface are constantly increasing without a local Venturi effect because of the stagnation of oncoming air at a high angle of attack. Again recall that, to minimize X h , a larger value of C l;h is favorable. Because C l;h has a close relationship with C m;h , increasing one of them generally causes the other to increase, and vice versa. The sign of C l;h is usually negative except at α ¼01. In most cases, the lift should increase with decreasing height due to the ground effect and the static height stability condition (H:S:). For instance, when a WIG vehicle encounters a disturbance and becomes out-of-trim, it decreases in height, which

Fig. 13. Drag polar with respect to h=c ¼0.15 and 0.40.

11

in turn increases the lift and brings the vehicle back to its original height. However, the value of C l;h at α ¼01 is positive for the three individuals, as shown in Fig. 16(a) and (b), which implies that decrease in height leads to lower the lift for a WIG vehicle. Consequently, the vehicle could be exposed to hazardous conditions, such as contacting the water or ground surface. As a result, the condition of α ¼01 should be avoided while cruising in ground effect. As shown in Fig. 16, the value of C l;h at α ¼01 is positive; thus, operation at this angle of attack is not acceptable and is excluded from this analysis. Furthermore, since this analysis considered only an airfoil without the other compartments of the vehicle, the values of X h can be either negative or positive. A more important factor in the successful design of a WIG vehicle is the location of X h ; i.e., the location of X h should be as close to the AC as possible. In this configuration, a sufficient stability margin (H:S:, X h −X α ) can be obtained, and CG is conveniently located between X h and X α , which is a strict condition for static height stability (Kornev and Matveev, 2003). As shown in Fig. 17, the locations of X h according to angle of attack move backward, or downstream from the AC. The aerodynamic centers of height (X h ) for the Pareto individuals are located at 10% (upstream) to −25% (downstream) of the chord length from the AC. The distance of X h from the leading edge is therefore about 15–50% of the chord length. The CG is located downstream of X h and the stability margin is reduced for higher value of X h . To compensate for the reduced stability margin and to satisfy the static height stability, it is necessary to have a large horizontal tail at high angles of attack, such as α ¼101. Such a high angle of attack should be avoided when the vehicle is in ground effect. Except at α ¼101, X h is bounded within 0.2, as shown in Fig. 17(a) and (b). The optimization is performed at α ¼21 and h=c ¼0.25, and the smallest absolute value of X h is observed for the three cases. The location of X h for p001 and p052 moves upstream as the height decreases, whereas the location of X h for p113 is stationary. For α ¼ 41 and 61, X h moves slightly forward for p001 and p052 and is constant for p113 as the height changes. One can see from Fig. 17 that the X h for p001 and p052 with an S-shaped lower surface is located upstream (i.e., more stable) than the X h for p113, which has a straight lower wing surface and the highest lift coefficient. Fig. 18 shows the C p variations of the p001 and p113 individuals with respect to angles of attack (α ¼01 and 61) and heights (h=c ¼0.15 and 0.4). The C p variations of height for both cases increase when the angle of attack is large (α ¼61). Consequently, the C p distributions on the upper wing surface are almost constant for decreasing heights. This small variation on the upper wing surface can also be found in a few studies (Ahmed and Sharma, 2005; Hsiun and Chen, 1996; Park and Lee, 2010); however, the increase in C p on the lower wing surface is considerable because the stagnation on the lower surface has a significant effect on the

Fig. 14. Comparison of C m with respect to heights: (a) p001 (filled) and p052 (hollow and marked), and (b) p001 (filled) and p113 (hollow).

12

S.-H. Lee, J. Lee / Ocean Engineering 68 (2013) 1–13

Fig. 15. Comparison of C m;h with respect to heights: (a) p001 (filled) and p052 (hollow and marked) and (b) p001 (filled) and p113 (hollow).

Fig. 16. Comparison of C l;h with respect to heights: (a) p001 (filled) and p052 (hollow and marked), and (b) p001 (filled) and p113 (hollow).

Fig. 17. Comparison of X h with respect to heights: (a) p001 (filled) and p052 (hollow and marked), and (b) p001 (filled) and p113 (hollow).

pressure distributions. Therefore, the increase in pressure distributions is not significant, and it does not decrease with height when the angle of attack is small. In this study, optimization is performed at α ¼21. The convergence–divergence passage effect is significant at α ¼01 for both p001 and p113. However, p001 – with a convex shape next to the trailing edge of the lower wing surface – has lower C p distributions than does p113. Although the shape change of the trailing edge is small, the convex shape has a significant effect on C m because the moment arm is large, resulting in a nearly constant C m;h at α ¼21. When the angle of attack is low, the lower wing surface of p001 does not appear to contribute to C m .

6. Conclusions A multi-objective optimization of a two-dimensional airfoil in ground effect was carried out using CFD and a multi-objective genetic algorithm. The design variables were parameterized by using Bézier curves, and a pre-optimization process with lowresolution design variables was performed four times in order to determine the lower and upper limits of the design variables effectively and guarantee sufficient resolution of their real values, The multi-objective optimization problem was solved using SMOGA without weighting factors and without normalization processes. Comparing Pareto sets from three optimizations with the reduced

S.-H. Lee, J. Lee / Ocean Engineering 68 (2013) 1–13

13

Fig. 18. Pressure distributions of two Pareto individuals with respect to heights: (a) p001 and (b) p113.

final bounds of design variables, the three Pareto frontiers are found to be very similar. This demonstrates that SMOGA can find potential solutions without considering weighting factors, despite employing indeterminacy operations such as selection and mutation. Based on the analysis of the Pareto optima including potential solutions, aerodynamic characteristics and stability are explored. Two objectives, C l and C l =C d , are projected onto a twodimensional plot. These objectives are globally dependent on each other, but locally independent. In addition, a unique solution (p113) is found with a straight lower wing surface. This surface can prevent the Venturi effect, which is generally observed from symmetric airfoils in ground effect, and can improve the ram effect. This shape may also reduce the drag and increase lift simultaneously, but it can degrade the stability provided by the constant pressure distributions on the lower wing surface. Furthermore, since lift (C l and C l =C d ) and stability (X h ) are traded off in designs, a Pareto individual such as p113 with high lift shows low stability. On the other hand, an individual, such as p001, with a convex surface next to the trailing edge similar to an S-shaped airfoil shows high stability because a local Venturi effect is present at the convex surface. However, such an individual loses lift slightly. Clearly, the shape of the lower wing surface is a key factor in designing a WIG vehicle. A straight lower surface improves the aerodynamic forces, and an S-shape next to the lower wing surface moves the aerodynamic center of the height upstream and improves the stability. Pareto optima obtained in this study can be applied to real WIG vehicle designs. Acknowledgments The authors would like to thank to J. Park for hours of discussion and for helping with the statistical analysis. References Ahmed, M.R., Sharma, S.D., 2005. An investigation on the aerodynamics of a symmetrical airfoil in ground effect. Exp. Thermal Fluid Sci. 29, 633–647. Ahmed, M.R., Takasaki, T., Kohama, Y., 2007. Aerodynamics of a NACA4412 Airfoil in Ground Effect. AIAA J. 45 (1), 37–47. Anderson, J.D., 1998. Aircraft Performance and Design. McGraw-Hill, New York. Booker, L.B., 1994. Improving Search in Genetic Algorithms. In: Davis, L. (Ed.), Genetic Algorithms and Simulated Annealing. Morgan Kaufmann Publishers, Los Altos, CA.

Carter, A.W., 1961. Effect of Ground Proximity on the Aerodynamic Characteristics of Aspect-Ratio-1 Airfoils with and without End Plate. NASA TN D-970. Chun, H.H., Chang, C.H., 2002. Longitudinal stability and dynamic motions of a small passenger WIG craft. Ocean Eng. 29, 1145 1126. DeJong, K.A., 1975. An Analysis of the Behavior of a Class of Genetic Adaptive Systems. Doctoral Thesis, Department of Computer and Communication Sciences University of Michigan, Ann Arbor. Fink, M.P., Lastinger, J.L., 1961. Aerodynamic Characteristics of Low-Aspect-Ratio Wings in Close Proximity to the Ground. NASA TN D-926. Goldberg, D., 1989. Genetic Algorithms in Search, Optimization and Machine Learning. Addision-Wesle. Hsiun, C.M., Chen, C.K., 1996. Aerodynamic characteristics of a two-dimensional airfoil with ground effect. J. Aircr. 33 (2), 386–392. Joh, C.Y., Kim, Y.J., 2004. Computational aerodynamic analysis of airfoils for WIG (airfoil-in-ground-effect)-craft. J. Korean Soc. Aeronaut. Space Sci. 32 (8), 37–46. Kikuchi, K., Motoe, F., Yanagizawa, M., 1997. Numerical simulation of the ground effect using the boundary element method. Int. J. Numer. Methods Fluids 25, 1043–1056. Kim, H.J., Chun, H.H., 1998. Design of 2-dimensional WIG section by a nonlinear optimization method. J. Soc. Naval Archit. Korea 35 (3), 50–59. Kim, Y.J., Joh, C.Y., 2004. Aerodynamic design optimization of airfoils for WIG vehicle using response surface method. J. Korean Soc. Aeronaut. Space Sci. 33 (5), 18–27. Kornev, N., Matveev, K., 2003. Complex numerical modeling of dynamics and crashes of wing-in-ground vehicles. AIAA2003–2600. Lee, J., Kim, B., Park, K., 2006. Aerodynamics characteristics and shape optimization of airfoils in WIG vehicle. Trans. Korean Soc. Mech. Eng. B 30 (11), 1084–1092. Lee, S.H., Lee, J.H., 2011. Optimization of three-dimensional wings in ground effect using multiobjective genetic algorithm. J. Aircr. 48 (5), 1633–1645. Mahon, S., Zhang, X., 2005. Computational analysis of pressure and wake characteristics of an aerofoil in ground effect. Trans. ASME 127, 290–298. Park, K.W., Lee, J.H., 2010. Optimal design of two-dimensional wings in ground effect using multi-objective genetic algorithm. Ocean Eng. 37, 902–912. Patankar, S.V., 1980. Numerical Heat Transfer and Fluid Flow. Hemisphere Publication Corporation. Poloni, A.C., Giurgevich, A., Onesti, L., Pediroda, V., 1999. Hybridization of a MultiObjective Genetic Algorithm, a Neural Network and a Classical Optimizer for a Complex Design Problem in Fluid Dynamics. Dipartimento di Energetica Universita di Trieste, Italy. Ranzenbach, R., Barlow, J.B., 1994. Two-dimensional airfoil in ground effect, an experimental and computational study. G.L.M Wind Tunnel 1 (287), 241–250. Rozhdestvensky, K.V., 2006. Wing-in-ground effect vehicles. Prog. Aerospace Sci. 42, 211–283. STAR-CD v3.15, 2001. Methodology. Computational Dynamics, Co, London, UK. Shin, M.S., Yang, J.Y., Yang, S.H., Wang, G.Q., 2000. Numerical simulation of viscous flow around a three-dimensional wing in ground effect with endplates. J. Ships Ocean Eng. 29, 29–39. Wieselsberger, C., 1922. Wing Resistance Near the Ground. NACATM77. Yakhot, V., Orszag, S.A., Thangam, S., Gatski, T.B., Spezoale, C.G., 1992. Development of turbulent models for shear flows by a double expansion technique. Phys. Fluids A 4 (7), 1510–1520. Zerihan, J., Zhang, X., 2000. Aerodynamics of a single element wing in ground effect. J. Aircr. 37 (6), 1058–1064.