|
Items in cart [1] |
Questions? Call 888-624-8373 |
| 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 282
A 2D+T VOF Fully Coupled Formulation for
Calculation of Breaking Free Surface Flow
Yann Andrillon, Bertrand A:lessandrini
(Ecole Centrale de Nantes - FRANCE)
~ Abstract
In this paper, we present some results of dimensional stea-
dy how. To obtain this results a 2D+T method have been
used with a two dimensional CFD code. The feature of
our approach, is the used of a 2D Navier-Stokes solver
with an interface capturing method to simulate complex
free surface. This code solves Navier-Stokes equation us-
ing a finite volume method adapted to structured or un-
structured mesh and used a fully coupled method to dis-
cretize the equation. The capturing method employed to
follow the interface is the "Volume Of Fluid" method.
This approach allow to simulate complex flow as break-
. .
1ng or merging wave.
2 Introduction
The main objective of our work is to be able to predict
the flow past a complex hull with breaking and splashing
waves. today, different code exist to predict steady wave
field. But less code could take into account the break-
ing bow splash and they are particularly computationally
expensive. Recently M. Landrini E1] have presented in-
teresting results from a 2D+T decomposition techniques
t2, 31. This approach present the advantage of being able
to give accurate wave field using a simple two dimen-
sional unsteady flow solver. However, it has much rec-
ommended for fine ships like combatants, and especially
for higher froude number where massive breaking of di-
vergent waves is observed. Then initially we restrict sim-
ulation to last kind of hull.
The solver used must solve complex unsteady free sur-
face flow. It solve the Navier-Stokes equation for incom-
pressible fluid and used a finite volume decomposition to
discretise the equations. Then structured as unstructured
mesh is accept to grid the simulation domain. The novelty
in our code is the used at the same time a Volume Of Fluid
capturing method and a verry accurateFully Coupled for-
mulation to solve the Navier Stokes equation. In order
to simulate an interface flow two different methods exist
the front capturing method and the front tracking method.
The last one have been investigated by several authors,
and it have been showed their limitation to predict merg-
ing wave. Harlow and Welch have introduced the first cap-
turing approach: the Marker and Cell (MAC). Later The
Volume Of Fluid method from Hirt and Nichols proved
the versatility and the accuracy of the capturing approach.
It will be described more precisely below and details can
be found in 0. Ubbink t4] and E. Didier ted. The Fully
Coupled Method [6, 73 have been recently investigated.
It has been showed that the method is particularly robust
and accurate. Numerically this technique present the ad-
vantage to accelerate the velocity-pressure coupling and
to converge faster than other weak coupling algorithms as
SIMPLE or PISO.
To check the solver capabilities, different two dimen-
sional flow case have been tested. The first one is the
liquid sloshing in a tank with various amplitude. The sim-
ulation has been compared to theoretical and experimental
data. Second one are the flow made by a collapsing col-
umn water and a Rayleigh Taylor instability. The results
illustrate the robustness of the solver and the ease to pre-
dict merging and rupturing interfaces. Third, the results
from a numerical two dimensional wave tank is presenting
to showed the capabilities of the code to propagate a wave
field. At least, a dimensional steady flow past a Wialey
ship at different froude number with breaking wave has
been simulated. Results prove the ability of the method to
predict breaking bow wave.
. . ~ ,
3 Volume Of Fluid methoc!
The "Volume Of Fluid" method is an interface capturing
method. The computation is performed on a fixed grid,
which extends beyond the surface. The shape of the inter-
face is determined by cell which are partially filled. The
main advantage of this approach is the ability to simu-
late complex free surface geometries like wave breaking.
Also, a new variable is introduced to the system, the vol-
ume fraction, which indicates, whether a cell is filled by
OCR for page 283
fluid 1 (c = 0), filled by fluid 2 (c = 1) or filled by a
mixture, (0 < c < 1~. The physical properties of fluid
is deducted from the volume fraction by the equations (1 )
and (2~.
p = cop + (1 - cope (1)
,~ = car + (1 - c),u2 (2)
This variable is ruled by a advection equation (3) which
are obtained from the second Navier-Stokes equation ( 1 1~.
The integral form of this advection equation is used with
a finite volume method for the spatial discretization. In
order to approximate the time integrals a Euler explicit
scheme is applied (4~. And a cell centered arrangement is
employed.
rt+~t—rt ~ ~ ~
,03~ + ~ (c~) = 0 (3)
dt -S+~;Cft~L=0 (4)
The critical issue in this type of method is the choice
of the differencing schemes for the convective term (4~.
Low order schemes like central differencing scheme are
not suitable because a bounded solution is not ensured.
Other differencing schemes like first-order upwind scheme
are too diffusive, smear the interface and introduce artifi-
cial mixing of the two fluids over a wide region. Then
high order mixing scheme have been designed for this
kind of application by Peric (HRIC scheme) [8, 9, 10] or
by Ubbink (CICSAM scheme) [43. The last one is partic-
ularly accurate to keep a sharp interface and do not create
non-physical deformation. This interpolation is a High
Resolution differencing scheme based on the Normalized
Variable Diagram t11, 12~. This approach is particularly
effective in the case of unstructured mesh.
use
/
Figure 1: normalization variable
cd = cd—cu ~5'
C'd—Cu
The CICSAM scheme use a mix of the compressible
scheme (Hyper-C) and Ultimate-Quickest scheme t11~.
This scheme use the Courant Number and weighting fac-
tor by, based on the angle between the interface and the
direction of motion. The aim of these factor is to avoid
stability problems and to limit overshoot or undershoot.
Another factor, kit introduce the control of the dominance
of the different scheme. More concretely, the UQ operates
where Hyper-C fails to preserve the gradient in the inter-
face and the Hyper-C operates where UQ fails to maintain
the sharpness of the interface.
CfCbc = { min (1, NC ~ if 0 < cd ~ 1 (6)
cot else
- ~ mitt ~ ~ ~ Cfcbc)
CfUq = ~ if 0 < cat < 1 (7)
cot else
Nc = Am, max(O, V ~ (8)
of = mitt (ok C05(~2§f + 1 1\) `9
At Least, the CICSAM have been evaluated on few
convection test with accurate results. The free surface
contour is sharply defined on two or three cells. How-
ever the boundedness criterion is checked only for weak
Courant Number (Nc < 0.1~. An additional subroutine
has been implemented to the last advection routine to treat
the unbounded volume fraction value.
4 Navier-Stokes Equations
The motion of the fluid is driven by the incompressible
Navier-Stokes equations. The dimensionless conserva-
tion mass and the conservation of momentum equations
(10,1 1) are presented below under the conservative form.
/v 0t /. dS+/ dV
is Rep ~ Ad) dS = +/ Fr2 ~dV (10)
/v Rep (~) (( At') dV
/ Aids = 0
s
with
· ~ is the velocity field vector
· p is the pressure
· g is the acceleration of gravity
· p is the mass density
(1 1)
OCR for page 284
. ~ . . .
· p IS dynamic VlSCOSlty
· Re is the Reynolds number
· En is the Froude number
In order to solve the governing equations, a fully im-
plicit finite volume method is choose. It could be apply
on unstructured mesh made by triangle or quadrilateral.
A cell centered arrangement is used to stored all depen-
dent variables. The integrals are calculated with a second
order approximation, and the flux approximation is eval-
uated with a deferred correction to obtain a second order
accuracy. The different term are evaluated at the new time
level, while the old values appear only in the second order
correction and in the approximation of the time derivative.
Uip + U'p + 9aP . Up + ~ 9i7t .Pn = 0 (12)
.0
iterations
Figure 2: Convergence of the system during no-linear
iteration
* an Fi chosen in order to use unstructured mesh. Nevertheless,
UiP + ~ a Uin = a (13) if the convection equation allow to transport an interface
zone wide of two or three elements, the discretisation of
the Navier-Stokes equation require in unstructured mesh
a wider zone. About 10 or more element is necessary to
define the free surface without creating abnormal veloc-
ity. Then in order, to make accurate simulation main of
our grid are structured one.
The pressure equation is derived from the discrete con-
servation momentum equation using a Fully-coupled me-
thod [7, 131. The new variable u*, v* are introduced in
the discrete momentum conservation equation (121. The
Rhie and Chow t14l velocity flux reconstruction applied
to the conservation of momentum Navier-Stokes equation
produce the discrete pressure equation presented below.
Pp + ~ c Pn—˘.4—~ <7 =— (14)
The different discrete equation give the following lin-
ear system.
—~ G ] (U) (Fu) (15)
O D —DO p Fp
In opposition to weak coupling, the fully coupled me-
thod do not require correction step, a single numerical
system is solved to obtain the solution (15~. An iterative
algorithm BiCGSTAB-w is used to solve the system pre-
conditionned by incomplete LU decomposition. One ad-
vantage of this technique is the numerical accelerates of
the coupling between velocity and pressure. So, the con-
vergence is fast, and it allow a reduction of 10 order in a
few number of non linear iteration as the figure 4 shows.
4.! unstructured aspect
The different method used for the discretization of the
Navier-Stokes equation and advection equations have been
5 2D+T method
This method has been introduced by M. Munk into aero-
dynamics for the prediction of loads on inclined slender
bodies of revolution and low aspect ratio wings. This
method has been extended in hydrodynamics domain, in
order to predict wave field due to a hull. More recently, M.
Tulin and M. Landrini [1, 2] have presented result of bow-
wave radiation. The analysis is limited to slender ships,
with a sharp stem. The idea exploited is that longitudinal
gradients of relevant flow quantities are small compared
with vertical and transverse gradients. Then a three di-
mensional steady flow past a ship hull could be generated
by a sequence of unsteady two dimensional problems in
vertical sections. The two dimensional flow are generated
by a deformable body, which correspond to a cross section
of the passing through ship. Then, the section is deducted
from the hull by a cut at the x=u.t coordinate, where u it
the boat speed and t the 2D solver time. It has been show
that this technique give accurate results for flow about a
ship at high Fn (> 1), and for low aspect ratio planning
boats or slender displacement hull((B/L)2 << 1~.
OCR for page 285
S.:l wavemaker
The deformable body described before used to generate
the 2D+T flow could be ranking as a special wavemaker.
To allow for our code to simulate this special wavemaker
a simple moving mesh technic has been implemented. We
define the displacement of the frontier node 16, and the
nodes inside the domain are moved using a translation cal-
culate as a linear function of the x-coordinate, like on the
figure 3.
J (`Ybody ~ =
Figure 3: example of moving boundary t=To and at
t=To+n.dt
The Navier-Stokes equation and the advection equa-
tion of volume fraction equation are modified as described
in Eland 18.
(iv cdv) + is Of (if—Ugrid)) ~dS = 0 `17'
f5t (/v REV) + is ~ (( - ugrid~.~) dS n
iv P + is Rep It\) dS (18)
iv Fr2 9 dV iv Rep (~ ( ~) dV
However, to limit the deformation of the cell near the
bottom of the hull, the ship used in simulation would be
as thin as possible.
6 2D Computation Procedure
The developed numerical method is applied to simulate
several kinds of unsteady free surface flows, from simple
to complex interface topologies. In order, to demonstrate
the versatility of the code, different cases as have been
simulated.
· Oscillating flow in a tank has been simulated t151.
This problem is used to compare our code to an-
other CFD code using front tracking method, and
to experimental data, showing good results.
· Collapse of a dam [16] is a typical problem used
for illustrating the versatility of the current solver
and the capability to predict merging and breaking
interfaces
A two dimensional Rayleigh Taylor instabilities.
(16) · 2D wave tank has been simulated numerically with
success. This 2D unsteady simulation is a particular
2D+T case.
6.l 2D sloshing tank
In this section, we present the flow in a sloshing tank. Tad-
jbakhsh & Keller t17] have developed a theory on slosh-
ing liquid in a tank under the influence of gravity. Ini-
tially the fluid has an interface defined by a half cosine
period with amplitude 0.005 m. The simulation domain is
a structured mesh 1 m long and 0.5 m height. The fluids
is respectively air and water without viscosity. The theo-
retical period of sloshing of the first mode is given by the
equation (19), where h is the depth of the tank and k is the
number of wave.
Ps = 27T>/gk tanh~kh) = 0.3739s (19)
Simulation
...... Theory
0.2'
_ . _
- 0 . 7 5 ~ _
0.5 1
t(s)
1.5 2
Figure 4: height of free surface on the left side
Theoretical result and numerical simulation are com-
pared on the the fig. 4. It show plots of position of the
interface at the left boundary against time for the six peri-
ods. It could be pointed out a good agreement.
Another sloshing flow have been simulated. it have
been compared to experiment made by Corrignan [151.
The tank 6.1 is moved with an horizontal displacement.
OCR for page 286
OCR for page 287
OCR for page 289
OCR for page 290
OCR for page 291
OCR for page 292
OCR for page 293
OCR for page 294
OCR for page 295
Representative terms from entire chapter:
wave field
The movement is ruled by the equ. (20~. The calculations
were performed on a structured and a Restructured mesh
of about 3500 control volume and the time step is about
dt = 0.001.
0.4 m
Figure 5: tank dimension
~ A.(~sin(~27rf~t)—sin(21rf2t))
x(t)= ~ siO
2a
Figure 7: montage colonne d'eau
1.2r
11
o.8t
Cal 0 . 6
0.4 _
0.2 _
n _
_ :
4 _
3.5
3
0.5 _
simulation
· exp. result
I ~ 1 .
, . . . . . . . . . . . . . . .
1 2 3 4
t~75
,~
simulation
· exp. result
1 2
t\/;77;
,
3
Figure 8: The position of the leading edge versus time
6.3 Rayleigh Taylor instability
This simulation showed the versatility of our code to dif-
ferent application domain. The figures 12 give the initial
configuration of the instability and the evolution of the
free surface at different time. The mesh used for the sim-
ulation is a Cartesian grid width 256x64 elements. The
time step is 0.001 s, and the running time is about 30 mn
Figure 9: collapse of a dam with an obstacle
on a HP station with an 500 Mhz Alpha processor for 6000
iterations. The numerical domain is box 1 x4 m, filled by
two fluid. the density of fluid in the bottom are 0.1694
kg/m3 and for the other it is 1.225 kg/m3. they have the
same viscosity 0.00313 kg/m.s . The numerical results
have been compared to other one obtained by Puckett and
al.~191) and it shows a good concordance.
6.4 two dimensional wave tank
The two dimensionnal wave tank allows to test the abil-
ity of the method to propagate wave without amplitude
damping or phase lag that are the two recurrent problems
for CFD hydrodynamics simulations. The first step try
to reproduce a monochromatic wave field, without break-
ing wave and check the accuracy with theoretical solu-
tion. The wave is generated with a piston wave makers as
it have been discussed previously in the wavemaker sec-
tion. The geometry of the numerical tank is 100 m long,
15 m depth and 10 m height. We use a damping zone of
200 m long downstream the outlet boundary to avoid the
wave reflection. The motion of the piston wave makers is
defined by:
x((t) = A.(~cos(
T = 4s ~ = 25m (24)
On the figures 13, the evolution of the free surface dur-
ing the installation of the established mode. On the fig. 14
the numerical solution is compared to the solution of the
linear water potential theory. The wave length don't be
the same one as that predicted by the theory. An largest
length seems to be more appropriate . Moreover, the wave
amplitude is undershoot against the theory too. Along the
wave tank, the amplitude is attenuated. However, it is im-
portant to notice that the simulation used a particularly
inaccurate mesh as it can be viewed on the fig. 14. Only
25 elements to define a wave length. Then, the wave is
correctly propagated for this simulation condition. To test
more qualitatively the cfd code to this kind of application,
it seems to be necessary to refine the grid.
7 3D Compulation Procedure
The well-known Wigley hull is chosen for the first hull
for, since it is a mathematical model made of parabolic
curves as eq. 25.
X2=B2 (1-(E) ) (1 (~) )
where x is the offset of the hull and L,B,d are the length,
breadth and draft of the ship.the length-to-breadth ratio
and length-to-draft ratio is set as below.
L, = 0.1 AT = 0.1
Figure 10: Wigley Hull
the flow past the ship has been calculated for different
Froude number. The computation domain is 2 m long and
1.2 m wide. The number of element is approximately of
18000. The time step adopted is function of Froude num-
ber and betwen 0.001 and 0.00025. The different CPU
time is about 1 6h on a HP DEC 500 MHz station.
Figure IS shows the free surface elevation along the
hull for different Froude Number. For the Froude Number
0.289 the result is compared to the experiment [20, 211.
The difference could be explain by the inaccuracy of the
mesh use for the simulation, and also due to the 2D+T
approximation. Moreover, the oscillation of the elevation
curve arise from the difficulties to locate the free surface
across the interface.
The following figures is the wave field made the wigley
hull. More precisely on the figure for 0.7 and 1.0 Froude
Number, we not wave breaking on the bow.
In fact, even more results with refined mesh is neces-
sary to approve the quality of the simulation, the presented
result give an idea on the power of the globally approach
and the robustness of our two dimensional flow solver.
~ Conclusion and Perspective
We have presented here a new numerical method using
both advantages: The accuracy of the Fully Coupled algo-
rithm first tested on Finite Differences convective formu-
lation for tracking interface method and the robustness of
Volume of Fluid capturing method to simulate very com-
plex interfaces.
The two dimensional test cases show that we achieved
our first objective.: the two dimensional solver is robust
and evade over a wide range of hydrodynamics appli-
cations.
The three dimensional test case present a method to
(25) obtain the flow past a fine hull, with accurate result.
Near future work will give more accurate results using
unstructured grids. Then examine the problem discussed
in the section "unstructured aspect" .To allow the code
to simulate 2D+T flow on more complex hull as serie 60
or DTMB model... Then thereafter it is planned to make
evolve the cfd code to a three dimensional flow solver.
-
- -
—amp ~
-A
-an
aim
- -
Figure 1 1: collapse of a dam at different instant: t = 0.2, 0.2, 0.4, 0.6, 1.0
Figure 12: Rayleigh Taylor instability at instant: O.Os;0.4s;0.8s;1.2s;1.6s;2.0s;2.4s;2.8s
3 ~
2
-
. .. . ... .. . ..
I t: 9.0000
. .. . . . . .. .. .. ... .
so
- ~ 10 20 30 40 5XO
o
-2
-3
3
2
1
o
-1
-2
3
3
2
1
o
-1
-2
-3 1 )
10 ~0 ~
. , .
... .............. .
t~5.0000
t:33.0000
-
10 20 30 40 5XO 60 io 80 90 100
Figure 13: numerical wave tank at different instant.
2 'I ~~ ~11411~
O ~ I _
-1 _ _ _ _ ~ ~ - m _ _
_,9 _ _ L _ ~ _ ~ _ _ ~,~ ~ .~ ~ .~ ~ .~ ~.~ I .~ _ _ _ .
. 111111~i
~l~Llr~lTR~
. ~~
~~ ~ PUTTY i L7
f@; ~
20 30 40
Figure 14: comparison between numerical result and theoretical.
0.04
0.03
0.02
0.01
o
-O.Ot
-0.02
-0.03
11111111111111111111111111 11111111111111111
11~:111111~111~111111111111 1111111111111 11
I 1~1~1 LL - tTI I t~ ~ 1~1 T I IllF
111111111~111111111111~111 · t:45.0000 ~
I.lllll~l.lT.11.7117iT.71,T71TIl',l,l,l,l,!lllll ~
60 70 80 90 100
W~ N~ ~
I, N`~ ~~`
Fn=0.289
exp. results Fn=0.289
.
Figure 15: free surface elevation along the hull
Figure 16: Wave field around the Wigley hull: Fn=0.289
Figure 17: Wave field around the Wigley hull: Fn=0.7
Figure 18: Wave field around the Wigley hull: Fn=l.O
Acknowledgement
The authors express their thanks to the Delegation Generate
pour l'Armement (DGA) which is supporting this work.
References
t1] M. P. Tulin M. Landrini, A. Colagrossi, "Numerical
studies of breaking compared to experimental obser-
vations," 4th Numerical Towing Tank Symposium,
Hamburg, Germany, 2001.
[2] M. Wu M. Tulin, "Divergent bow waves," 2Ist
Symposium on Naval Hydrodynamics, 1997.
[31 M. Wu M. Tulin, "Bow waves on fine ships - nonlin-
ear numerical studies," 9th Int. Workshop on Water
Waves and Floating Structures,Kyushu, Japan, 1994.
[4] O Ubbink, Numerical prediction of two fluid
systems with sharp interface. PhD thesis, University
of London, 1997.
[5] E Didier, "Simulation d'ecoulements a surface libre
sur des maillages Restructures," These de Doctorat,
2001.
t6] X. Vasseur, Etude Numerique de techniques
d'acceleration de convergence lors de la resolution
des equations de Navier-Stokes en Formulation
Decouplee ou Fortement Couplee. PhD thesis, Uni-
versite de Nantes, 1991.
[7] M. Ferry, Resolution Des Equations de
Navier-Stockes Incompressibles En Formulation t19
Vitesse-Pression Fortement Couplee. PhD thesis
Universite de Nantes, 1991.
[8] M Peric and J H Ferziger, Computational Methods
for Fluid Dynamics. Springer,second edition, 1997
[9] R Azcueta, S Muzaferija, and M Peric, "Computa-
tion of water and a~r flow around ships," Euromech
374, Poitiers, 1998, pp. 121-132. [21]
[10] M Peric, "Basics of viscous flow CFD," CFD
for ship and offshore design, 31st Wegemt School,
1999.
[11] B P Leonard, "Bounded higher-order upwind mul-
tidimensional finite-volume convection-diffusion
agorithms," dans W.J. Minkowycz, E.M. Sparrow
(eds), Advances in Numerical Heat Transfer, Chap
1, Taylor and Francis, New York, 1997, pp. 1-57.
[12] H Jasak, H C Weller, R I Issa, and A D Gosman,
"High resolution NVD differencing scheme for ar-
bitrarily unstructured meshes," Web published,
Vol. 0, 1996.
t13l B Alessandrini and G Delhommeau, "A multi-
grid velocity-pressure-free surface elevation fully
coupled solver for turbulent incompressible flow
around a hull calculations," Proc 9th International
Conference on Numerical Methods in Laminar and
Turbulent Flows, Atlanta, 1995, pp. 1173-1184.
rl4] C M Rhie and W L Chow, "A numerical study of
turbulent flow past an isolated airfoil with trailing
edge separation," AIAA Journal, Vol. 21, 1983,
pp. 179-195.
P Corrignan, "Analyse physique des phenomenes as-
socies au ballotement de liquide dans les reservoirs
(sloshing)," These de Doctorat, 1994.
[16] W.J. Moyce J.C. Martin, "An experimental study
of the collapse of liquid columns on a rigid hori-
zontal plane.," Philos. Trans. Roy. Soc. London,
Vol. A244, 1952, pp. 312-324.
t17] J.B. Keller I. Tadjbakhsh, "Standing surface waves
of finite amplitude," J. Fluid Mech., Vol. 442-451,
1960, p. 8.
[18] S Koshizuka, H Tamako, and Y Oka, "A particle
method for incompressible viscous flow with fluid
fragmentation," Computational Fluid Dynamics
JOURNAL, Vol. 4~1), 1995, pp. 29~6.
J.B. Bell E.G. Puckett, A.S. Almgren, "A high or-
der projection method for tracking fluid interfaces
in variable density incompressibl flows," J. Comp.
Phys., Vol. 100, 1997, pp. 269-282.
F. Larrarte, Etude experimentale et theorique des
profies de vagues le long d'une carene. PhD thesis,
Universite de Nantes, 1994.
C.Y. Chen F. Noblesse, "Comparaison beetween
theorical predictions of wave-resistance and exper-
imental data for the wigley hull," J. of Ship
Research, Vol. 27, 1991, pp. 215-226.
DISCUSSION
Hoyte C. Raven
MARIN, Netherlands
In the discussion on Figure 15, you attributed the
imperfect agreement of the hull wave profile at
Fn=0.289 to an insufficient mesh density.
However, at such Froude numbers transverse
waves still play a role; and transverse waves
cannot be modeled by 2D+t. Transverse wave
length here is about 0.5L, and this may well
explain some of the disagreement.
