| Copyright © 2009. National Academy of Sciences. All rights reserved. Terms of Use and Privacy Statement |
Below are the first 10 and last 10 pages of uncorrected machine-read text (when available) of this chapter, followed by the top 30 algorithmically extracted key phrases from the chapter as a whole.
Intended to provide our own search engines and external engines with highly rich, chapter-representative searchable text on the opening pages of each chapter.
Because it is UNCORRECTED material, please consider the following text as a useful but insufficient proxy for the authoritative book pages.
Do not use for reproduction, copying, pasting, or reading; exclusively for search engines.
OCR for page 128
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
OCR for page 129
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)
OCR for page 130
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)
OCR for page 131
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
OCR for page 132
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
OCR for page 133
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
OCR for page 134
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)
OCR for page 135
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.
OCR for page 136
\\\\\~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
OCR for page 137
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
OCR for page 138
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.
OCR for page 139
[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.
OCR for page 140
of
initial shape
^ k-~o model
O Rim model
(a) Model scale
f
initial shape
G k- 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
OCR for page 141
~ \ \\ 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
OCR for page 142
//
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
OCR for page 143
V.;iit...,~__~8=~=:,~
(a) Initial
(a) Initial (b) Final
Figure 21: Velocity homogenization: streamlines at full scale
(b) Final
OCR for page 144
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.
Representative terms from entire chapter:
initial shape