Prsentation PowerPoint

1
L.S.E.E.T
University Of Rhode Island
Recent progress in the modeling of non-linear
free surface phenomena in ocean engineering
FRAUNIE, P. (1) GRILLI, S.T. (2)

• L.S.E.E.T, Université de Toulon
• Department of Ocean Engineering, University Of
  Rhode Island

2
WAVE BREAKING IN COASTAL ZONE

• Objectives
• Coastal morphodynamics, sediment transport
• Dammages on coastal areas and structures
• Ocean-atmosphere interactions
• Tools
• - Laboratory/in situ Experiments
• - Numerical Simulation Numerical wave
  tank

Photos www.pvergnaud.free.fr
3
Free surface flows

• Flows containing several fluids/phases
• Several examples
• - wave breaking
• - cavitation
• - slushing of fuel in satellite tanks
• better understanding of the
  physical phenomena
• Tool numerical simulation using interface
  tracking methods

4
Mathematics and numerical modeling
5
Assumptions what is to be modeled ?

• Flows
• Fully 3-D
• Unsteady
• Non hydrostatic
• Laminar/turbulent
• Single or two-phase flows
• Fluids
• Newtonian
• Incompressible or not

6
Mathematical formulation conservation equations
ï
r
(
f
ï
t
in the fluid domain
7
Mathematical formulation Interface boundary
conditions
Interface velocity
Velocity continuity at the interface
Normal to the interface
Viscous fluids only
Interface curvature
Stress balance at the interface
Surface tension coefficient
8
Mathematical formulation boundary conditions

• Solid boundaries slip condition (Euler) or
  no-slip condition (Navier-Stokes), pressure
  extrapolation

• Open boundaries
• - Dirichlet condition (fixed velocity and
  pressure) on inlet boundaries
• - Neumann condition (normal derivative of the
  velocity imposed to zero) for the velocity and
  fixed pressure on outlet boundaries

9
Mathematical formulation summary
Resolve the system composed by
Equation governing the evolution of the interface
Let be C the binary function so that
10
Numerical resolution of the conservation
equations
CFD code EOLE (PRINCIPIA RD)
Navier-Stokes (or Euler) equations in a
curvilinear formulation (?,?,?)
F,G,H flux terms (convective, diffusive,
pressure)
J Jacobian of the transformation
11
Numerical resolution of the conservation
equations
With

• Space discretization Centerred Finite Volume
  scheme (fields computed at the cell center)

• Time discretization second order implicit
  scheme

12
Pseudo-compressibility method (Chorin 1967)
Concept introduction of a time-like variable t,
the pseudo-time and of pseudo-unsteady terms
New unknown introduced in the pseudo-unsteady
terms, the pseudo-density
Additional equation pseudo equation of state
giving the pressure as a function of the
pseudo-density (Viviand)
13
Pseudo-compressibility method

• Adding artificial viscosity terms avoids
  numerical oscillations
• Integration step by step in pseudo-time thanks
  to a five step Runge-Kutta scheme ( dual time
  stepping )

Convergence solution independent on t
corresponding to the numerical solution at time
level n1
Robust method allowing to deal with high density
ratios
14
Algorithm
15
Development and validation of a 3-D Larangian
V.O.F method
16
Interface tracking method aims

• 3-D method allowing to deal with large
  deformations of the interface (large curvatures,
  reconnections, deconnections )
• Accuracy
• Fast compared to classical V.O.F methods

17
Volume Of Fluid (V.O.F) concept

• The interface is tracked thanks to the volumic
  fraction of the denser fluid (fluid 1)
• C 1 in a full cell of fluid 1
• C 0 in a full cell of fluid 2
• 0 lt C lt 1 if a cell an interface occurs in the
  cell

18
V.O.F methods
C function is advected with the fluids and
verifies the transport equation

• Classical discretization schemes (centred,
  upwind, Quick ) are diffusing the interface and
  are not accurate
• Alternative methods with interface
  reconstruction.
  Several possibilities
• SOLAVOF method (Hirt Nichols, 1981)
• CIAM (Li, Zaleski, 1994)
• SL-VOF (Guignard, 2001, Biausser 2003)

19
SL-VOF 3-D method (B. Biausser, 2003)
3 steps allowing to update the interface during a
time step

• Interface reconstruction
• Interface advection
• Computation of the new V.O.F field

20
Step I interface reconstruction
 Piecewise Linear Interface Calculation  (Li 1994)
In each cell, the interface is represented by a
plane portion (intersection of a plane with the
computational cell)
Interface
21
Step I interface reconstruction
Two steps to calculate the interface plane
portion

• Definition of the plane direction
• Translation of the plane in order to verify the
  volume of the cell

Calculation of the plane direction the normal
to the plane (orientated from denser the fluid
towards the less dense fluid) is
22
Step I interface reconstruction
Evaluation of n by finite differences from the
V.O.F of the neighbouring cells a) Computation
of normal vectors at the 8 corners of the cell
b) Normal vector is the mean of the 8 normal
vectors at the corners
23
Step I interface reconstruction
Plane translation the normal to the plane and
the volumic fraction Cijk of the cell determine a
unique plane
Translation of the plane so that the volume
contained under this plane is Cijk
If the equation of the plane of normal nijk (nx,
ny, nz) is nxxnyynzz ? , the problem is
equivalent to calculate ?(Cijk,nijk)
24
Step I interface reconstruction
The calculation of ?(Cijk,nijk) provides a unique
plane portion polygon from 3 to 6 sides whose
corners are known
B
G
A
H
(a)
(b)
(c)
(d)
(e)
(f)
25
Step II interface advection
Calculation of the velocity at the polygon
corners bilinear interpolation from the
velocities computed by the solver at the cell
center
Corners advection first order (in time)
Lagrangian scheme Xn1 Xn U.?t
26
Step II interface advection
After advection, the advected polygon corners are
not necessarily coplanar so that a mean plane to
these corners is defined
Normals to triangular subdivisions
P3
Mean normal of the normals to triangulars subdivis
ions
P4
Pm
Polygon corners after advection
P2
Mean point iso-barycentre of the corners
P1
Mean plane of the corners after advection (nm 
mean normal  , Pm  mean point)
nm
Pm
27
Step III computation of the new V.O.F field

• Two configurations after advection
• Cells containing polygons portions (A type)
• Cells without interface (B type)

type A cells after advection
28
type A cells treatment
Calculation of the mean plane to all polygons
parts in the cell Mean plane defined by
averaging the normals to the plane parts and
their centres (weighted with the portions
surface)
29
type A cells treatment
The new VOF of the cell is the volume generated
by the averaged plane and is calculated by
inversion of the formulae giving ? as a function
of Cijk and nijk
n
New VOF Cijkn1
30
type B cells treatment

• Two configurations are possible
• Cells loosing the interface during the time step
  such cells become full (C 1) or empty (C 0)
  following the stream direction
• Cells without interface before advection the
  value of C remains the same

Cell without interface before and after advection
Cell filled up during advection
31
SL-VOF 3-D summary

• 3-D V.O.F. method with geometrical
  reconstruction of the interface
• PLIC modeling (more precise than HirtNichols)
  allowing to deal with large deformation of the
  interface
• Lagrangian advection (possibly use of larger
  time steps than with classical methods)

32
Evaluation of the methods performances

• Comparison with a classical 3-D V.O.F method
  already developed in the EOLE code FLUVOF (Hirt
  Nichols kind)
• Aability to deal with large changes of the
  interface
• Ability to use large time steps

33
Comparison with FLUVOF 3-D

• Comparison with a HirtNichols method using a
  constant piecewise reconstruction of the
  interface (0 order)

• Pure advection test-case (imposed velocity)
  allows to test the methods performances without
  NS solver
• Advection to a wall the analytic velocity is
  known
• A sphere advected in such a flow is
  progressively changed into ellipsoïds

34
Comparison with FLUVOF 3-D computational domain
In each transverse plane y constant
35
Sphere advection in a distorting velocity field
SL-VOF 3-D simulation
36
Comparison with FLUVOF 3-D
Mesh 100X30X100
37
Comparison with FLUVOF 3-D

• When the curvature is maximum, accuracy is
  better with SL-VOF 3-D than using FLUVOF
  advantage of the P.L.I.C discretization of the
  interface
• The SL-VOF simulation is 4 times faster than
  FLUVOF advantage of the Lagrangian advection
• Volume conservation is quite good 0.13 of
  loss compared to the initial fluid volume after
  70 time steps

The methods approximates (mean plane) are not
penalizing the volume conservation
38
Coupling with the NS solver Rayleigh-Taylor
instability

• Stratified fluids of different densities (the
  denser is above)
• Initial perturbation ? characteristic
  instability involving local vortices
• Overturning of the interface occurs and the flow
  is computed with the full solver good test for
  the method

2-D example (denser fluid in red)
39
Rayleigh-Taylor instability

• Density ratio 2
• Perfect fluids in a cylindrical domain
• Interface initially plane sinusoidal
  perturbation of the velocity

40
Comparisons with 2-D axisymetric results
3-D on a radius
2-D axi
41
Conclusions about the test cases

• Compared to a classical V.O.F method better
  accuracy when the curvature is increased,
  computational time is reduced

• Large accurately deformations are taken into
  account

42
Wave breaking applications
43
First tests of breaking
Evaluation of the methods ability to simulate
wave breaking

• Breaking of an unstable linear wave
• Breaking of a solitary wave on a beach of slope
  1/15

44
Breaking of an unstable linear wave

• Sinusoidal wave of high camber
• Initial velocity field Airy wave
• Periodic boundary conditions over one wavelenght
• L 0.769 m
• T 0.86 s
• D 0.1 m
• H 0.1 m

Fast evolving towards a plunging breaker
45
Propagation direction
Déferlement dune onde linéaire instable
46
Modulus of the velocity
47
Conclusions concerning this test-case

• Results comparable to those of Abadie (1998) on
  the same test-case for 2-D flows (aspect of the
  breaker jet, splash-up, maximal velocity about 2
  times the phase celerity)
• First simulation of breaking conclusive with the
  method (reconnections and deconnections of the
  interface, curvature )
• Artificial breaking, generated by a non-physical
  initial condition

48
Breaking on a beach of slope 1/15

• Solitary wave H0 0.5 m
• Computational domain flat bottom and then
  sloping bottom
• Initialisation with Tanakas algorithm (1986) and
  computation of the initial fields with Boundary
  Integral Equations Method of S. Grilli
  potential code using a Boundary Element Method

49
Boundary Integral Equations Method

• Nonlinear potential flows with a free surface
• Fast and accurate method for wave shoaling and
  overturning applications
• Unable to deal with breaking (no reconnection,
  irrotationnal and inviscid flows )

50
Solitary wave initialization
51
Soliton breaking
Sloping bottom ? kinetic energy is transferred
into potential energy ? camber ? breaking
52
Soliton 2-D results and experimental results
SL-VOF method for 2-D flows has been tested
successfully on wave breaking applications
Result of the simulation of the breaking of a
solitary wave on a bottom of slope 0.0773 - Weak
coupling BEM - SL-VOF
PIV image of the breaking of a solitary wave on a
bottom of slope 0.0773 Experiment made in the
waterl tank in ISITV

53
Soliton comparisons 2-D/3-D

• Test-case runs for 2-D flows compared to
  measurements and BIEM by Guignard (2001)
• Comparison 2-D SL-VOF / 3-D SL-VOF very close
  results
• Few differences (delay for breaking) due to the
  coarser mesh for the 3-D run

54
Conclusions concerning this test-case

• Successful simulation of a physical breaking
• Comparisons with 2-D results ok (similar
  results)
• Pseudo-3-D test-case no variation of the slope
  in the cross direction ? same phenomenon in each
  transverse plane

55
Breaking of a solitary wave on a sloping ridge

• Sloping bottom with a transverse modulation with
  a hyperbolic secant
• Slope 1/15 at the centre of the ridge and 1/36
  on each side
• 350 cells along x, 40 along y and 65 along z
• Solitary wave H0 0.6 m
• Coupled to Grillis BIEM
• Single phase flow in order to reduce
  computational time

350 cells along x, 40 along y, 65 along z
56
Initialisation with BIEM

• First step a part of the shoaling is computed
  using BIEM (accurate and faster than the
  VOF/Navier-Stokes solver but unable to deal with
  breaking)
• The solution of this first simulation is used as
  an initialization of the VOF/Navier-Stokes solver
  (free surface, velocity and pressure)
• The end of the simulation (overturning, breaking
  and post-breaking) is computed with the
  VOF/Navier-Stokes model

Propagation direction
Initial condition for the VOF/Navier-Stokes model
57
Overturning stage

• The breaker jet occurs bottom variations leads
  to the the wave camber and overturning
• Focusing of the energy at the center of the
  ridge because of the steepest slope the breaker
  jet firstly occurs at the center

58
Overturning slices along the x-axis
3-D aspect of overturning the wave is begining
to break at the center while the breaking point
is not reached on the sides
Vertical cross-section along x at y 0 m
Vertical cross-section along x at y 2 m
59
Pressure at breaking point
Due to large vertical accelerations, the pressure
is not hydrostatic in front of the wave
Ratio of the computed pressure to the hydrostatic
pressure in the slice y 0 m
60
Velocity field
Transversal velocity (slice z 0.3 m) focusing
High velocities in the breaker jet high
accelerations (4.9g)
61
Breaking
62
Conclusions concerning this simulation

• Mesh of 900,000 cells (?x 5 cm then 2.5 cm, ?y
  10 cm, ?z ?1.5 cm in the breaking zone) CPU
  time 5 days and 10 h on a Digital Dec alpha
  bi-processor 500 MHz
• Breaking simulation with SL-VOF 3-D consistent
  with the BIEM simulation before breaking
  (focusing, values of the physical fields,
  interface aspect )
• Mass conservation loss is 0.7, Energy
  conservation loss is 10 ? loss of amplitude
  during shoaling and delay to break with respect
  to the time predicted by BIEM
• Errors numerical diffusion (coarse mesh along
  y and x in the shoaling zone), single phase-flow
  run

63
Breaking
3-D aspect of breaking impact occurs firstly at
the center for x19.85 m and progressively on the
sides
64
Post-breaking
The wave continues to collapse, the air tube is
progressively crashed the water jet is
projected with high velocity along the slope
Water jet
65
Conclusions

• Development and validation of a 3-D interface
  tracking method in a CFD code
• PLIC modeling and lagrangian advection ?
  accurate and fast method when compared to
  classical VOF methods
• Efficient method for wave breaking applications
• Loss of energy during the shoaling stage
  numerical diffusion of the CFD code (mesh,
  artificial viscosity, single phase flow )

66
Energy loss numerical diffusion

• Tests sur la discrétisation et les modes
  diphasique/monophasique
• Shoaling dune onde solitaire en fond plat et
  évaluation de la perte dénergie totale

67
Numerical accuracy (pure advection)
Advection of a sphere
Err1
Order 1.65
68
Critical VOF
 
With n(nx,ny,nz)