Twenty-Fourth Symposium on Naval Hydrodynamics (2003)

Below is the uncorrected machine-read text of this chapter, intended to provide our own search engines and external engines with highly rich, chapter-representative searchable text of each book. Because it is UNCORRECTED material, please consider the following text as a useful but insufficient proxy for the authoritative book pages.

On the Role played by Turbulence Closures in Hull Shape Optimization at Model and Full Scale R.Duvigneau, M.Visonneau and G.B. Deng (Ecole Centrale de Nantes, France) Abstract The practical use of automated CFD-based design tools in ship-building industry requires powerful flow solvers, able to take into account realistic geometries as well as complex physical phenomena, such as turbulence. A shape optimization tool is developed in this framework. A derivative-free optimizer, yielding both flexibility and robustness, is preferred to the classical gradient-based method, more fastidious to implement and still limited to moderately complex problems. The flow solver in- cluded in the design procedure solves the incompressible Reynolds-averaged Navier-Stokes equations on unstruc- tured grids using a finite-volume formulation, involving several near-wall low-Reynolds number turbulence mod- els. The design tool is then applied to optimize the stern of a modern hull shape at model and full scale, different purposes being considered. Our interest is particularly fo- cused on the influence of turbulence modeling in the de- sign process. The effects of a two-equation model based on the eddy-viscosity assumption and a second-order clo- sure relying on the Reynolds stress transport equations are compared. ~ Introduction \ Computational Fluid Dynamics (CFD) may be considered today as practical analysis tools for the ship-building in- dustry and is integrated in fully automated optimization procedures during the design phase. Since the calcula- tions are computationally expensive, gradient-based opti- mization algorithms are usually employed in the design procedure [13. However, this approach is still limited to moderately complex problems, since the evaluation of the derivatives of the cost function requires to differentiate the flow solver, representing a cumbersome work when so- phisticated flow solvers are used. Moreover, these meth- ods are particularly sensitive to the noisy errors arising from the simulation process during the evaluation of the cost function, generating irregularities and spurious local optima [21. All these remarks make this strategy question- able. Genetic algorithms are an alternative approach L3l, but the high number of evaluations required still discards these methods at the present time. Furthermore, a so pow- erful and expensive approach is not always needed for some mild optimization problems. The present paper is devoted to the study of shape optimization procedures, in the particular framework of highly complex CFD problems, with a special emphasis on the influence of turbulence closures in the design pro- cess. In order to take into account realistic geometries as well as complex flows, the code ISIS developed in our lab- oratory is implemented in the procedure. Incompressible Reynolds-Averaged Navier-Stokes Equations (RANSE) are solved on unstructured grids, using a finite volume formulation. Turbulence modeling is achieved by near- wall low-Reynolds number models, ranging from the one- equation model of Spalart-Allmaras [4], k—w closures [5] to a full Reynolds stress transport Rij—w model [61. Multi-fluids calculations using surface tracking or surface capturing methods are also possible. The flow solver is in- cluded in a derivative-free optimization procedure, based on linear or quadratic interpolations. This approach pro- vide the capability to solve optimization problems involv- ing complex flow solvers, since the solver is considered as a black box, yielding both flexibility and robustness. In order to test the capabilities of the design tool, the optimization of the stern of a modern ship, the KVLCC2 tanker, is presented. Two goals are successively consid- ered, at model and full scale. First, a reduction of the drag is performed. Then, the characteristics of the flow at the propeller location are improved, homogenizing the longi- tudinal velocity field. However, it is now well admitted that turbulence closures play a central role in the simula- tion of flows in this area. Therefore, the influence of the turbulence closures on the design process is investigated, comparing the effects of a two-equation model based on the eddy-viscosity assumption and a second-order closure relying on the Reynolds stress transport equations. The first part of this paper is devoted to the description of the numerical methods involved in the present work. The whole design procedure, which is the framework of

this study, is first presented. Then, the tools used inside are detailed, including the flow solver, the optimizer, the parameterization tool and the mesh update strategy em- ployed. In a second part, the applications on hull shape optimization are presented, the results analyzed and dis- cussed. 2 Design Procedure 2.l Problem statement Shape optimization consists in determining the design pa- rameters D to minimize a cost function f (X, U), depend- ing on metric elements X(D) and flow variables UtD), both entities being related to the shape. The state equa- tions R(X, U) = 0 are assumed to be verified at each step of the optimization procedure. Some geometrical constraints, bounding the variation domain of the design parameters, are usually imposed. Physical constraints de- pending on the problem studied may also be taken into account. Finally, the problem may be expressed as: find D to minimize constraint to 2.2 Optimization loop f(X(D), U(D)) R(X(D), U(D)) = 0 Ci(X(D), U(D)) > 0 Ce(X(D), U(D)) = 0 A shape optimization procedure corresponds to an itera- tive search in the design parameters space, moving from a shape to an other, until an optimum is reached. For each optimization step k and each shape Dk, a mesh adapted to the shape should be built, providing the metric elements X(Dk). Then, the state equations are solved by the flow solver, providing the flow variables U(D k) corresponding to the current shape. Then, the fitness f (X (D k ), U (Dk ~ ~ of the shape Dk can be evaluated. Finally, this informa- tion is used by the optimizer to predict a next shape D k+t, expected to be better. The whole optimization loop is: determine the initial shape D ken 1 2. if stopping criteria reached then STOP 3. generate a grid adapted to D k and evaluate the metric elements X(Dk) 4. solve the state equations by the flow solver, provid- ing the flow variables U(Dk) 5. evaluate the fitness f (X(Dk), U(Dk)) of the current shape 6. predict an improved shape Dk+t by the optimizer 7. update k ~ k+1 goto 2 As described above, three main tools are included in the shape optimization procedure: the flow solver, the mesh update tool and the optimizer. The following sections de- tail the strategies and techniques employed for each of them. 3 Flow Solver The ISIS flow solver, developed by the DMN (Divi- sion Modelisation Numerique i.e. CFD Department of the Fluid Mechanics Laboratory), uses the incompress- ible unsteady Reynolds-averaged Navier Stokes equa- tions (RANSE). The solver is based on the finite vol- ume method to build the spatial discretization of the transport equations. The face-based method is gen- eralized to two-dimensional, rotationally-symmetric, or three-dimensional unstructured meshes for which non- overlapping control volumes are bounded by an arbitrary number of constitutive faces. (1 ) The velocity field is obtained from the momentum con- servation equations and the pressure field is extracted from the mass conservation constraint, or continuity equa- tion, transformed into a pressure-equation. In the case of turbulent flows, additional transport equa- tions for modeled variables are solved in a form similar to the momentum equations and they can be discretised and solved using the same principles. Incompressible and non-miscible flow phases are mod- elized through the use of conservation equations for each volume fraction of phase. 3.1 Conservation equations l~he flow solver can deal with multi-phase flows and mov- ing grids. Considering incompressible flow of viscous fluid under isothermal conditions, mass, momentum and volume fraction conservation equations can be written as (using the generalized form of Gauss' theorem): bt Iv pdV + Is P(t—td) · ~dS = 0 ~t |V pUidV + Is PUi(t - td) · ~dS = | (7ijIj—pIi) ~dS + | P9idV s V 0t |V CidV + Is Ci(t—td) · ~dS = 0 (2a) (2b)

where V is the domain of interest, or control volume, bounded by the closed surface S with a unit normal vector directed outward. The effective flow physical properties (viscosity and density) are obtained from each phase physical properties (pi and Pi) with the following constitutive relations: p=Ecipi; p=ECipi; 1= When the grid is moving, the so-called space conserva- tion law must also be satisfied: by; fv dV Is Ice Adds = 0 (4) 3.2 Numerical framework All the flow variables are stored at geometric centers of the arbitrary shaped cells. Surface and volume integrals are evaluated according to second-order accurate approx- imations by using the values of integrand that prevail at the center of the face /, or cell C, and neighbor cells Cnb (Fig. 1~. ~ ~ C nb no Figure 1: Arbitrary control volume The various fluxes appearing in the discretised equa- tions are built using centered and/or upwind schemes. For where example, the convective fluxes are obtained by two kinds of upwind schemes. A first scheme available in the flow solver, (HD) for Hybrid differencing, is a combination of upwind (UD) and centered (CD) schemes. Contrary to a 0 = C kS practical approach t7, 8] where CD/UD blending is fixed t' with a global blending factor for all faces of the mesh, the HD scheme results from a local blending factor based on the signed Peclet number at the face. An other upwind scheme which is implemented in ISIS is the Gamma Dif- ferencing Scheme (GDS) [93. Through a normalized vari- able diagram (NVD) analysis [101, this scheme enforces local monotonicity and convection boundedness criterium (CBC) [1 13. A pressure equation is obtained in the spirit of the Rhie and Chow [12] procedure. Special attention is given to the all way unsteady terms, 0/0t, and pseudo-unsteady terms, 0/~r, are interpolated so that the overall solution does not depend on the time step At and on the local fictitious time step A\. When running on parallel machines, the computational domain is splitted into multiple-connected domains hav- ing approximately the same number of unknowns. This is done with the help of the MeTiS fly, 14] partitioning algorithm. Communication of faces data between domains is per- formed according to the Message Passing Interface stan- dard El53. 4 Turbulence models 4.l Reynolds Stress Mode! The Reynolds stress transport equations can be written as: Dl; = Pig + Oij—Aid— ~~ (Ciik —fax UiUj) The production terms are given by: Pij (UiUk 0~ + fink ~ ) An isotropic model for the dissipation rate is used: cij = 35ij~ A linear model with wall reflection terms is retained in the present study for the pressure-strain terms: (pij = C(~)ij + ¢)(2)ij + (~(W)i `~(~)ij = -Cabin + C3k (bikSjk + bjkSik — 3bmrlSm~5ij) + C4 k (bikWjk + bjkWik)

and the wall reflection terms proposed by Gibson & Launder (1978) [16] are given by: · To predict a correct behavior of the law of wall and a good estimation of the wall friction velocity for sim- ple wall flow. This requirement is guaranteed by the ~ ~ 3 good prediction of the shear stress. 0(W)ij = Cut k ~ Ukum~kNm~ii — 2UkUi~knj \ 3 ~` k3/2 2UkUjNkni) C sy + Cw2 (~(2)km~knm~ij- 2°(2)kinknj 3 ~ k3/2 2°(2)kjnkniJ C ey In the above formula, Sij and Wij are the strain rate and the rotation rate tensors respectively: 1 {0Ui bUj ~ sij = - <_ + _J W ~ (bUi bUj) bid is the Reynolds stress anisotropy tensor defined as bid = 2k —3§ij ni is the wall normal vector, and Yw is the wall nor- mal distance. The Daly & Harlow model is used for the diffusion correlation terms: C C k buin ~ Axe A variety of linear pressure-strain correlation mod- els have been proposed. The most popular models are the Launder, Reece and Rodi (1975) t17] LRR model, the isotropization-of-production (IP) model, and more re- cently, the Speziale, Sarkar and Gatski (1991) t18] (SSG) model. It may be noticed that, when using the IP model, the 0(2)ij can be simplified as: ( ) 2 ( 3 ) The above model is only valid for high Reynolds num- ber Bows. A near wall modification needs to be included in order to be able to apply the model in the wall region. When calibrating a near wall model, the most important requirements that should be taken into account are the fol- where the function fC remains to be defined. lowing: The first two terms representing a linear eddy viscos- ity model are treated implicitly, while the last term is · To satisfy the numerical constraints so that all equa- treated explicitly. The fit,, function is calibrated such that lions can be integrated to the wall. · To predict correctly the normal stress anisotropy. · To predict a correct wall limiting behavior of the Reynolds stresses. In the present near wall model, the Car coefficient is just replaced by: Car ~ 2 + (C~—2)tar7,h ( ) The remaining coefficients are given below: = A — = A — 3 — 0.4 mirl(O 3, A) 2.5 0.22 1—gA2 + ~,A3 aija aijajka The Cw2 formula proposed by Hanjalic & Jakirlic (1998) [19] is used here since it is found useful to avoid any unphysical separation. The Reynolds stresses determined from transport equa- tions are not directly used in the momentum transport equations. Actually, the effective Reynolds stresses in the momentum equations are computed as: + fc (2L,tS:; +UiUi—35iJk) where k2 Z/t = 0 09/i1— fly = 1-e- 20 Ry = man (Y~4,0.006YU'~ ~) ~ V L/ J the implicit term gives a good approximation of the shear

stress for simple wall boundary layer flow. The choice of the correction function fC iS flexible. With fC = 0, the present model is reduced to an eddy-viscosity model even Reynolds stress transport equations are solved to de- termine k and £. Although this choice is not interesting for turbulence modelization, a first run with fC = 0 can give a good initialization of the Reynolds stress. With fC = 1, the above implementation can be considered as a defect correction approach. However, even with this choice, the implicit term does not cancel the correspond- ing explicit term. In one-dimension for example, the implicit discretization of the implicit term involves grid points (i-l,i,i+1), while the explicit discretization of the source term involves grid points (i-2,i,i+2~. The differ- ence between the implicit and the explicit discretization provides a numerical dissipation which is beneficial for the numerical stabilization. For the computation of complex flows such as the flow around the HSVA tanker [20], it is found that the com- putation with fC = 1 fails to converge. For such a com- plex flow, the validation of the wall reflection model cal- ibrated from simple wall boundary layer flow becomes questionable. However, this term has no direct contri- bution to the prediction of the turbulent kinetic energy. Therefore, we may expect that reasonable prediction of the kinetic energy can still be obtained when solving the Reynolds stress transport equation. Consequently, an eddy-viscosity model based on the predicted kinetic en- ergy can provide a reasonable approximation. To obtain a stable solution, we have chosen the correction fc as fc = 1 -e- 15 In this case, the present model is reduced to eddy- viscosity model in the near wall region. 4.2 Turbulent Frequency Equation for the Reynolds Stress Mode! The turbulent dissipation £ iS determined from a turbulent frequency w equation rather than from a turbulent dissi- pation equation. When y ~ 0, the w equation is decou- pled from any other equations. Consequently, it is less stiff compared with the ~ equation. In addition, a near wall model for w is easier to implement than for a. The transport equation for w used in the present study is given below. Model Cal c~2 COtw LRR 1.55 1.90 0.63 IF 1.53 1.92 0.60 SSG 1.53 1.92 ~ 0.60 Table 1: Model coefficients for the w transport equation Dw Dt Where = or ~ = = Cw = fw = C - ~ — aW ph _ ~W2 1~ ~ ~ k dw bw ~ CW—UkU! + ~ ark ~ ~X' Oak 2 2 ( ~ + C k ) 0w Irk (CC! - 1) fw + COW (1 - fw) Cp (Ce2 - 1) 0.18 tank ~ ( 0.002~yW ) = C`,kw 0.09 The coefficients Cat, Cot and C,,LW given in the follow- ing table are calibrated for different pressure-strain model by using the experimental data for a boundary layer under zero pressure gradient. Away from the wall, the w equation returns to the clas- sical ~ equation with the constants Cal and Cog given in table 2. Unlike in a near wall model for ~ where a damping function is often used to change the order of the dissipa- tion term, the damping function fin in the above model is used to change the value of the or coefficient as well as to cancel the cross diffusion term in the near wall region. Therefore, such a model is more universal. 5 Parameterizations The choice of the parameterization technique employed is considerable since it exerts influence on the search and the shapes obtained. Furthermore, some techni- cal requirements are leading the choice, as discussed by Samareh [211. The parameterization should be able to represent complex three-dimensional shapes, involving a number of variables as small as possible, in order to fa- cilitate the search. Then, the modified shapes have to

remain smooth, whatever the perturbations imposed by the optimizer. Moreover, the use of existing grids is ex- pected, which forbid the use of Computer Aided Design (CAD) softwares interacting with the procedure. A usual approach consists in representing the shape by B-spline curves or surfaces, the design parameters being the co- ordinates of the control points. However, the application of this method is difficult for three-dimensional problems in the particular framework of unstructured grids, since a curvilinear representation of the shape is not available. The nodes coordinates on the shape and connectivities are the only knowledge which could be used. \S ' \ <~ _ \ ~ An', '' '. 1. ~ tr~~~ I ~I.~~ :~ Figure 2: Parameterization of the stern of the ship Therefore, a Free Form Deformation (FFD) tech- nique [22] is used to control the shape perturbations dur- G Mesh update ing the design process. It consists in embedding in a box the object to be deformed, and then modifying the space in the box and the object inside by deforming the box, rather than modifying the object itself. In that manner, the shape of the object can be modified without even identifying its nature. Actually, a local coordinate system is imposed on a parallelepipedic volume including the object to be de- formed. Any point P in this volume has coordinates in this system such that: P = bet + 7'e2 + (e3, (5) where en, e2 and e3 are vectors of length unity, parallel to the leading directions of the parallelepipedic volume. A deformation function is then applied to the space inside the volume, by the use of a trivariate B-spline product. The new position of the previous P point is provided by the algebraic relation: PFFD = ~ Ni CONS Clank (~)PBS . . . 3,}.~ where Hi, Nj and Nk are the B-spline basis functions and . . . PBS are the Cartesian coordinates of the control points. By moving some of these control points, the space in- side the parallelepipedic volume is perturbed, yielding a smooth deformation of the shape inside. This approach provides an easy-to-use manipulation tool, since no curvi- linear representation of the object is required. The example of the parameterization of the stern of a ship is given, illustrated by the figure 2. A volume is de- fined around the part of the ship to be modified, mapped by a trivariate B-spline product of 5 x 6 x 3 control points. The shape of the ship may then be modified by moving some control points, for instance six control points in fig- ure 2. They are located strictly inside the volume in order to ensure a smooth transition between the deformed part of the shape and the part which remains unaltered. For each new shape proposed by the optimizer, a new grid should be built around it, before the flow solver determine the solution fields. As mentioned previously, the flow solver ISIS may take into account complex grids, with- out any assumptions concerning the nature of the control volumes. This point of view is expected for the design tool as well, to be able to include an automated mesh re- finement procedure for instance. Moreover, the ability to use existing grids is expected. Therefore, the possi- bility to use automated grids generation softwares is re- jected. As a consequence, the mesh update is performed by deformation of the initial grid. Thus, the topology of the grid is maintained and only the nodes coordinates are modified. Different grid deformation techniques for un- structured meshes exist, usually based on structural analo- gies t213. However, their application may be tedious in practice, when three-dimensional unstructured meshes adapted to viscous calculations are employed. Because (6) of the use of near-wall turbulence models, the distance between the first nodes and the wall is really low, creat- ing highly stretched volumes which generate problems to the deformation algorithms. As a consequence, the PPL' method is chosen to deform at once the shape and the grid around. Indeed, the relation 6 establish a spatial deforma- tion, which may be used to move the nodes of the grid. Since the volumic deformation is controlled by B-spline functions, the modification of the mesh is always smooth. However, the quality of the resulting grid is not assessed. Nevertheless, satisfactory results are obtained in practice, assuming perturbations are usually moderate. As example the figure 3 and 4 show the deformation of a two-dimensional grid around a NACA 0012 airfoil, composed of a structured grid near the wall and a trian- gular mesh in the outer domain. As expected, a smooth deformation 1S performed, with no crossing line although a severe perturbation is imposed. Meanwhile, a stretching of the grid is observed in some areas near the wall and the

orthogonality of the initial mesh in the structured part is not maintained. Figure 3: Initial mesh Figure 4: Final mesh 7 Optimization algorithm 7.l Strategy The optimizer plays a crucial role in the design procedure, since it predicts at each optimization step an improved shape, from the information collected previously. The use of gradient-based optimizers is usually motivated by their efficiency from a numerical point of view, since they can reach a minimum of the cost function in a low number of evaluations. However, some difficulties are arising when they are faced with complex realistic problems. The eval- uation of the derivatives of the cost function, calculated through a sophisticated simulation process, is rather prob- lematic. They are usually evaluated using an adjoins for- mulation, relying on the differentiation of the flow solver [23, 241. This task is tedious when high-order discretiza- tion schemes on unstructured grids are used, or complex turbulence models are employed. This approach often needs an a-priori simplification of the problem, neglecting turbulence and free surface variables for instance, or us- ing first-order discretization schemes, which provides an approximated gradient. It is shown in E251 that these sim- plifications often lead to erroneous gradient values. Thus, this approach is still limited to moderately complex prob- lems. Moreover, these algorithms are very sensitive to the noisy errors arising from the evaluation of the cost function. Indeed, the fitness of a shape is not determined exactly, some errors corrupting the exact value f en. Actu- ally, an approximated fitness fap is only evaluated: fap = fee + Elf + If + If (7) where Elf represents the error of discretization, If the error of convergence and If the error of modeling. The magnitude of the first term depends on the grid size and the discretization scheme employed. The second term comes from the fact that only partially converged solu- tions are computed. These two sources of errors are stud- ied by Madsen [2], who underlines the high-frequency errors introduced by the use of high-order discretization schemes and low converged solutions. The last term cor- responds to the error introduced by the turbulence models and will be examined later. Gradient-based methods are particularly sensitive to the high-frequency errors coming from the incomplete convergence, which generate irreg- ularities and spurious local minima. Thus, strongly con- verged and expensive flow analysis are required by this approach. To overcome these limitations, a derivative-free algo- rithm is employed, which is easier to implement in a com- plex numerical framework. It may be associated with so- phisticated flow solvers, since the solver is considered as a black box, which has not to be modified in order to be included in the design procedure. Furthermore, this ap- proach is less sensitive to the noise, because no informa- tion about the derivatives is needed to predict the opti- mization path, and only moderately converged solutions of the state equations may be used, reducing the calcula- tion costs. 7.2 Algorithm The optimizer chosen for this work rely on an iteratively updated linear or quadratic interpolating model of the cost function, developed by Marazzi and Nocedal t261. As- suming that the value of the cost function is known at a set of sample points including the best point ok, a local interpolating model is built at each optimization step k. It may be linear, including n + 1 points for a problem of . . . Amen son n: Ik (Xk + S) = f (ok ) + gk S S ~ urn (8) or quadratic, for which 2 (in + 1) (in + 2)—1 sample points are needed: ~kfXk+S) = f(xk)+9kS+S Irks s ~ Urn (9)

The unknown coefficients Ok and Hk are determined by solving a linear system. Assuming that the model chosen mk is valid in a sphere of radius Pk centered in Xk, called the trust-region, a subproblem is solved at each iteration to minimize mk in this region. This cheap operation may be easily performed using a simple gradient-based algorithm, for instance. The new point obtained amok may then be included in the set of sample points to replace X far the furthest point from Xk. However, an other constraint is added to the subproblem in order to ensure that the re- placement of afar by amok will not deteriorate the model at the next iteration. Indeed, it could be easily shown that, if amok belongs to a certain linear or quadratic subspace Ok, the resulting set of sample points will not be adequate to determine the model magi. To prevent from this case, the following subproblem is solved: minimize mktx) x constraint to ~ ~ X—Xk ~ ~ < Pk X ~ (9k- where Ok iS a prescribed region embedding Ok. In that manner, the accuracy of the interpolating model is im- proved during the update. In practice, the previous sub- problem is first solved without taking into account the last constraint. If the point obtained violates the constraint, a correction is made to move away from the region ~ k. The details of this operation are found in [261. After a solution Amos of the subproblem 10 is determined, the value of the cost function at this point is calculated. Amok is then in- cluded in the set of sample points if the current iteration is a success f (Amos) < f (Xk) orif Amok is closer to Xk than Of ar. The second case promote the inclusion of points in the vicinity of ok in the set of sample points, improving the quality of the model when unsuccessful trials are per- formed. The radius of the trust region is then updated. In case of success, it is increased since the model seems to be accurate and larger step may be performed, else it is reduced. Finally, the procedure goes on using the updated model and trust region. convection terms is achieved by the hybrid differencing scheme described previously. The influence of the turbu- lence closure is tested by comparing the results obtained by the two-equation SST k—~ model of Menter t5] and the second-order closure presented above. For all the cal- culations performed, the computational time required by the second-order closure is about twice more important than by the two-equation model. Although no experimen- tal data is available at full scale to validate the calcula- tions, the design test cases are performed at both model and full scale, for Reynolds numbers Rem = 4.6 106 and Ref = 109, since different behaviours are observed. A systematic study of the respective influence of mod- eling and discretization errors performed on a similar hull, indicated that the modeling errors were far more impor- tant than the discretization errors for this moderately re- fined grid. These conclusions were indirectly confirmed (10) by the similar results obtained by the CFD teams partic- ipating in the Gothenburg 2000 workshop. This is why it has been decided to keep such a moderately refined grid, considering the very large amount of computations required by hull shape optimizations. The successful prediction of the hook-shaped mean streamwise velocity contours at the propeller location x/l = 0.485, observed in the experimental measure- ments, is considered as a key criterion for assessing the performance of turbulence models designed for ship flows. A detailed discussion on the different behaviours of the turbulence closures compared to experiments for the KVLCC2 ship is given in [61. The isowakes at the pro- peller location are shown on figures 5 and 6 at model and full scales. As seen, hook-shaped contours are captured for both models, but the vortex predicted by the second- order closure is more intense and closer to experimental results than the one given by the two-equation model. By comparing full scale to model scale results, it is observed that the overall intensity of the longitudinal vorticity is drastically reduced at full scale, since the convection ef- fect is more important. The differences between the mod- els are lower but still exist. ~ Application to hull shape opti- ' · ~ mlzatlon 8.l Study of the flow around the KVECC2 tanker Computations are performed on a domain covering the whole ship, using a structured mesh with 121 x 81 x 41 nodes in the streamwise, radial and girthwise directions respectively. The free-surface plane is considered as a horizontal symmetry plane. The discretization of the As mentionned previously, the turbulence closures have a major influence on the prediction of the flow at the stern of the ship, providing different values of the drag and dif- ferent characteristics of the velocity at the propeller loca- tion. Meanwhile, these two criteria are usually considered as goals for ship design optimization, trying to reduce the drag and to homogenize the velocity field on the propeller. Therefore, an influence of turbulence modeling may be expected during the design procedure. This point is ana- lyzed in the next sections.

\\\\\~Rij-m model k-m model /// ~ ~ \~,~j-m model ~ ~~/ Figure 5: Isowakes at model scale 8.2 Drag Reduction The first exercise of hull shape optimization consists in trying to reduce the drag by modifying only the stern of the KVLCC2 ship. To achieve this goal, 6 control points are positionned on a vertical face of the FED control box as shown in figure 2. Actually, two complementary ob- jectives are pursued in this part of the present study. The respective influence of the turbulence closure on the final optimized shape has first to be assessed independently for model and full scales in order to evaluate the influence of the "modeling noise" on the optimization procedure. Then, the optimized solutions at model and full scales will be compared in order to check if the overall design proce- dure depends strongly on the Reynolds number scaling. X.2.1 Model scale For the model scale, the drag is reduced by 2% in about 60 RANSE evaluations with a moderate loss of hull volume. The figure 7 shows the evolution of the drag for the k—~ SST and Rij—w turbulence closures. Despite different values, one can observe that the rate of reduction of the drag looks very similar and does not depend strongly on the "modeling noise". However, one must notice that the reduction of drag is smaller than the difference between the drags predicted by different turbulence closures. The shapes of the final optimized hulls for each tur- bulence closure are shown in figures 11 to 13 for differ- ent x-stations. As expected, the reduction of drag results globally in an evolution of the hull from a U-shape to- wards a V-shape for both turbulence closures. However, the shape resulting from the RSM computations appears slightly more complicated than the one provided by the k—w SST computations. The difference between shapes are correlated to different flow fields as shown in the fig- ures 14 related to the model scale. One can observe that, even if the global evolution of the stern shape is similar during the optimization, the stern flow modeled by the k—`v SST closure loses a major part of its longitudi- nal vorticity whereas the flow modeled by the Rij—w Figure 6: Isowakes at full scale closure keeps a noticeable amount of longitudinal vor- ticity although concentrated in an narrower region near the vertical plane of symmetry. Actually, during the opti- mization process, the shape of the hull converges towards a V-shape, and it is well known that such a shape does not produce intense longitudinal vortices, contrary to U- shape. However, the strong dependence of the flow fields with respect to turbulence closures which was observed on the initial shape, is not significantly reduced as long as the sections of the KVLCC2 stern get thinner during the optimization iterations. 0.000300 0.000298 0.000296 0.000294 0.000288 0.000286 0.000284 0.000282,, ~ k-m model 0 Rij-c,) model 10 20 30 optimization iterations Figure 7: Evolution of the drag at model scale 8.2.2 Full scale The same exercise is repeated at full scale on a mesh com- prising the same number of points but located in such a way that the "y+ = 0~1) near the wall" criteria is satis- fied. The figure 8 shows the evolution of the drag for the k—w SST and Rij—w turbulence closures, respectively. Here again, the values predicted by both turbulence mod- els are significantly different and one may notice that the

k—~ SST computations encounter difficulties to reduce the drag after the fourth optimization iteration, indicat- ing that a local minimum may have been reached by the optimization tool. On the other hand, the Rij—~ com- putations seem to produce a more regular evolution of the shape, since no plateau is observed prematurately on the drag evolution. The dissimilarities concerning the drag evolution are not confirmed by the study of the x-stations shown in the figures 1 1 to 13 related to full scale. One can observe that the optimized hulls provided by both turbu- lence closures at full scale are similar, although the hull provided by the Rij—w closure is thinner, a fact which may be correlated with the saturation on the drag evolu- tion observed with k—~ computations. However, it seems that the turbulence closure plays a less important role at full scale than at model scale as indicated by the figures 15 related to the full scale configuration. This reassuring and expected result is largely due to the fact that this drag reduction exercise tends to produce, at full scale, similar V-shape hulls, decreasing the overall vorticity. 8.2.3 Comparison of results at model and full scale A cross-comparison of the shapes at model and full scale is also provided by the figures 1 1 to 13 and figures 14 and 15 for the flow fields at the stern. From the author's point of view, these figures seem to indicate that it may be risky to trust an optimization at model scale if full scale ships are considered. It is interesting to compute the full scale drag on the hulls optimized at model scale to have an idea of the error related to the Reynolds number scaling. For the k—~ model, the full scale drag on the shape optimized at model scale is about the same as the drag on the shape optimized at full scale. This result is probably related to the poor improvement obtained during the optimization at full scale. Meanwhile, when the second-order closure is used, the decrease of the drag at full scale performed by the shape optimized at model scale is about the half of the decrease observed for the shape optimized at full scale. 0.000143 _ ( ~ 0.000142 0.000141 Cal 0.000140 0.000139' 0.000138 ~~ A ~ A A A A ~ k-m model 0 Rij-o model ~~s - I 1 ~ ~ 1 , , 1 1 1 , , , , 1 , 1 1 1 1 0 1 0 20 30 40 optimization iterations Figure 8: Evolution of the drag at full scale of the flow field in the disk propeller may result in bet- ter propulsive performances and a reduction of the source of vibrations. Therefore, to achieve this goal, the mean deviation of the longitudinal component of the velocity is computed in the propeller disk and used as cost function. As shown in figure 9 and 10, the mean deviation is sig- nificantly reduced during the design process, the magni- tude of the decrease being 66% at model scale and 72~7o at full scale. One may underline that the mass flux through the propeller disk is slightly increased of some pourcents at the same time. The figures 17 to 19 show the modi- fications of the shapes at model and full scale for several longitudinal stations. One can observe again that the over- all trend is similar from model to full scale, since a bump is generated near the bottom of the stern, the deviation from the initial hull shape being far more pronounced at full scale than at model scale. The explanation of such a behaviour is provided by the figure 16 which shows the isowakes in the propeller plane, the circle representing the propeller disk, where the flow field is homogenized. It is clear that the main consequence of the birth of the bump S.3 Homogenization of the flow at the pro- at the bottom of the hull is the generation of a far more intense longitudinal vortex which achieves satisfactorily the homogenization of the flow field in the propeller disk. Since the initial full scale flow is less vertical than the model scale flow, the deformation of the hull at full scale has to be more important to obtain the same reduction of the mean deviation. At the same time the longitudinal vortex is getting more intense, the drag is increased of 3% at model scale and 13°~7o at full scale. This difference was again predictable, comparing the modifications of the flow fields at model and full scale. One can easily imagine that the same optimizations performed with a Rij—~ model pelter location The goal of this second optimization exercise is by many ways very different from the previous one. Instead of try- ing to optimize a quantity attached to the hull like the drag, this is the control of the velocity distribution in the region where the propeller should be present, which is now considered. Because of the lack of time, all the opti- mizations at model and full scale were performed with the same k—~ SST turbulence closure. This naive but com- plex exercise is based on the idea that the homogenization

would have resulted in a less pronounced deformation of the hull at model and full scales since this deformation ap- pears closely related with the ability of the computations to correlate hull shape and production of longitudinal vor- ticity in the near-wake flow. Lastly, the figures 21 show the wall streamlines for the initial and final hulls at model and full scale. Whereas the pattern of the wall streamlines are very different on the initial hull at model and full scale, it appears that the behaviour of these wall streamlines be- comes more similar when the hull shapes are optimized. In particular, one can notice the development of a long convergence line characteristic of an intense longitudinal vortex. From all these observations, one can state that the in- formation given by the model scale computations is use- ful for the full scale problem, but only in terms of trend. When the fitness of the shape optimized at model scale is evaluated at full scale, one can observe that about 3/4 of the improvement is performed, compared to the optimiza- tion at full scale. Moreover, since the homogenization of the flow is strongly dependent on the vertical character- istics of the near-wake flow, one can deduce that this op- timization will be strongly dependent on the turbulence closure. 0.2 0.18 _ ~ 0.16 o .Ct a' ~5 Cal ~ 0.12 I` o 0.14 0.1 0.08 _ 0.06 - - - \_ -, , . , 1 , , · , 1 , , , , 1 , , , , 1 0 5 10 15 20 optimization iterations Figure 9: Evolution of the mean deviation at model scale 9 Concluding remarks This article was devoted to the study of the influence of turbulence modelisation and Reynolds number scaling on the optimization of hull shape. In order to reinforce the overall robustness of the computational approach, a derivative-free optimizer has been chosen. Two different design problems have been considered in this study, which 0.24 0.22 ~ 0.2 o .'t 0.1 8 ~ 0.16 Ct a) ~ 0.14 - o 0.12 a) 0.1 0.08 0.06 ,: , \ 'I ~ 0 10 20 optimization iterations Figure 10: Evolution of the mean deviation at full scale have both underlined the Reynolds number dependency. For the full scale optimizations, this paper has shown that the role played by turbulence closure is strongly depen- dent on the kind of optimization considered. Whereas a classical SST k—~ closure provides satisfactory results for a drag optimization, one can suspect that it should be risky to trust it when an homogenization of the flow field is aimed at. Therefore, the role played by the turbulence closure appears rather reinforced even for full scale flows as soon as complex optimization objectives are consid- ered. References [1] T. Hino, "Shape optimization of practical ship hull forms using navier-stokes analysis," Proceedings of the 7th International Conference on Numerical Ship Hydrodynamics, 1999. [2] J.I. Madsen, "Response surface techniques for dif- fuser shape optimization," AIAA Journal, Vol. 38, No. 9, 2000, pp. 1512-1518. [3] R. Duvigneau and M. Visonneau, "Single- and multi-objective optimization for high-fidelity cfd using genetic algorithms," EUROGEN 2001 - Evolutionary Methods for Design, Optimisation and Control with Applications to Industrial Problems, September 2001. [4] P. R. Spalart and S. R. Allmaras, "A one-equation turbulence model for aerodynamic flows." AIAA Paper 92-0439, 1991.

[5] F.R. Menter, "Zonal two-equations k—~ turbulence models for aerodynamic flows." AIAA paper, 1993. [6] G. Deng and M. Visonneau, "Comparison of explicit algebraic stress models and second-order turbulence closures for steady flows around the kvlcc2 ship at model and full scales" A Workshop on Numerical ship hydrodynamics, L. Larsson, F. Stern, and V. Bertram, eds., Goteborg, September 2000. [7] I. Demirdzic and S. Muzaferija, "Numerical method for coupled fluid flow, heat transfert and stress anal- ysis using unstructured moving meshes with cells of arbitrary topology," Comput. mesh. Appl. Mech. Eng., Vol. 125, 1995, pp. 235-255. [8] J.H. Ferziger and M. Peric, Computational methods for fluid dynamics. Berlin: Springer-Verlag, 1996. [9] H. Jasak, Error Analysis and Estimation for the Finite Volume Method with Applications to Fluid Jul 1999 Flows. PhD thesis University of London 1996. Y , , [10] B.P. Leonard, "Simple high-accuracy resolution pro- gram for convective modelling of discontinuities," International Journal for Numerical Methods in Fluids, Vol. 8, 1988, pp. 1291-1318. [11] P.H. Gaskell and A.K.C. Lau, "Curvature om- pensated convective transport: SMART, a new boundedness preserving transport algorithme," International Journal for Numerical Methods in Fluids, Vol. 8, 1988, pp. 617-641. [12] C.M. Rhie and W.L. Chow, "A numerical study of the turbulent flow past an isolated airfoil with trail- ing edge separation," AIAA Journal, Vol. 17,1983, pp. 1525-1532. [13] George Karypis and Vipin Kumar, "A fast and high quality multilevel scheme for partitioning ir- regular graphs," Tech. Rep. 95-035, University of Minnesota, Department of Computer Science, Min- neapolis, MN 55455, 1995. Last updated on March 27 1998. [14] MeTiS, level Partitionning http://www-users.cs.umn.edu/ karypis/metis/. "Family of Multi - Algorithms." . . . . t15] MPI, "The Message Passing Interface (MPI) stan- dard." http ://www-unix. mc s. anl . gov/mpi [16] M.M. Gibson and B.E. Launder, "Ground effects on pressure fluctuations in the atmospheric boundar~y layer," Journal of Fluid Mechanics, Vol. 86, 1978, pp. 491-511. [17] B.E. Launder, G.J. Reece, and W. Rodi, "Progress in the development of a reynolds stress turbulent clo- sure," JournalofFluidMechanics,Vol. 183, 1975, pp. 63-73. t18] C.G. Speziale, S. Sarkar, and T.B. Gatski, "Model- ing the pressure-strain correlation of turbulence: an invariant dynamical systems approach," Journal of Fluid Mechanics, Vol. 227, 1991, pp. 245-272. [19] K. Hanjalic and S. Jakirlic, "Contribution towards the second-moment closure modeling of separating turbulent flows," Computers & Fluids, Vol. 37, 1998, pp. 147-156. t20] G. Deng and M. Visonneau, "Comparison of ex- plicit algebraic stress models and second-order tur- bulence closures for steady flows around ships," Proceedings of the Seventh International Conference onNumerical Ship Hydrodynamics, Nantes, France, [21] J.A. Samareh, "A survey of shape parameterization t e c h n i q u e s , " C E A S. / A I A A / I C A S E / N. A S. A L a n . ~ 1 e y International Forum on Aeroelasticity and Structural Dynamics, June 1999. [22] T.W. Sederberg and S.R. Parry, "Free-from defor- mation of solid geometric models," Computer Graphics, Vol. 20, No. 4, 1986, pp. 151-160. t23] W. K. Anderson and V. Venkatakrishnan, "Aero- dynamic design optimization on unstructured grids with a continuous adjoins formulation," Computers and Fluids, Vol. 28, No. 4, 1999, pp. 443-480. [24] A. Jameson, L. Martinelli, and N. A. Pierce, "Op- timum aerodynamic design using the navier-stokes equation," Theorical and Computational Fluid Dynamics, Vol. 10, 1998, pp. 213-237. [25] W. K. Anderson and E. Nielsen, "Aerodynamic de- sign optimization on unstructured meshes using the navier-stockes equations," AIAA Journal, Vol. 37, No. 11, l999,pp. 1411-1419. [26] M. Marazzi and J. Nocedal, "Wedge trust region methods for derivative free optimization." Report Optimization Technolo~y Center, October 2000.

of initial shape ^ k-~o model O Rim model (a) Model scale f initial shape G k-<i) model / O Rij-`,> model (b) Full scale Figure 11: Drag reduction: shapes at x/1 = 0.42 D ; I I~ hap n t s e ', ~ initial shape ~ k-m model O Rij-a) model l ~ k-c,) model .` '. O Rij-`o model ~ ;; (a) Model scale If (a) Model scale (b) Full scale Figure 12: Drag reduction: shapes at x/1 = 0.44 f initial shape G k-m model O Rij-m model Figure 13: Drag reduction: shapes at x/1 = 0.46 (b) Full scale initial shape G k-m model O Rij-m model

~ \ \\ initial final All / ~ (a) k Figure 14: Drag reduction: isowakes at model scale initial final A/// / ~ ~~/ ~ (b) Rij—~ model ~ In initial final (a) k—w model to' ~ (a) Model scale (b) Rij—co model Figure 15: Drag reduction: isowalces at full scale \ As initial final ~ \ Ah/ ~ (b) Full scale Figure 16: Velocity homogenization: isowakes using the k—~ model

// D' initial shape G k-m model - (5 '~ (a) Model scale Figure 17: Velocity homogenization: shapes at x/l initial shape k-a, model (a) Model scale (b) Full scale Figure 18: Velocity homogenization: shapes at x/1 = 0.44 (a) Model scale initial shape G k-m model (b) Full scale Figure 19: Velocity homogenization: shapes at x/1 = 0.46 ~ ^ k-m model (b) Full scale =0.42 OK · initial shape f initial shape G ken model / initial shape G k-`o model

V.;iit...,~__~8=~=:,~ (a) Initial (a) Initial (b) Final Figure 21: Velocity homogenization: streamlines at full scale (b) Final

DISCUSSION C. Yang George Mason University, USA I believe that authors have done excellent research in the shape optimization. We have done the similar work using CFD tools together with gradient-based optimization techniques to minimize wave drag. We have found that the hull form optimized for a single design speed may yield less desirable hull form for other speeds. Have authors tried multi-speed optimization or checked hydrodynamic characteristics of the optimal hull form obtained at other speeds? I would also like to know what is the geometry constraint and how many design variables are used in the current optimization process? AUTHORS' REPLY For this optimization exercise, only 6 control points moved in the cross direction are considered as design variables. Since the shape modification is quite local, the only constraints included concern the range of displacement of the control points. No constraint involving the volume for instance is taken into account. The comparison of the results obtained for different Reynolds numbers shows for these exercises similar tendencies for the shapes, although some differences are noticed which may have a significant influence on the final performance. More global problems involving wave drag minimization were not tested. Multi- objective optimizations using genetic algorithms may provide an interesting answer for these multi-points problems.