M.V Salvetti, F. Beux, M. Bilanceri (University of Pisa)

1
A numerical method for barotropic flow simulation
with applications to cavitation

• M.V Salvetti, F. Beux, M. Bilanceri (University
  of Pisa)
• E. Sinibaldi (now at Scuola Superiore SantAnna,
  Pisa)

Micro-Macro Modelling and Simulation of
Liquid-Vapour Flows Strasbourg, 23-25 January 2008
2
Industrial and engineering motivation
Development of a numerical tool for the
prediction of the performance of axial inducers
typical of turbopumps for liquid propellent
rockets.
rotating inducer
non rotating cylindrical case
The main role of the inducer is to increase the
fluid pressure (velocity decrease) through
rotation.
Fundings from the Italian Space Agency and
European Space Agency.
3
Difficulties

• Complex rotating 3D geometry.
• Severe size limitations which lead to very high
  rotational speed ? cavitation phenomena ? need of
  a model to take into account cavitation.

Choice of the cavitation model
Considering the characteristics of the considered
engineering problem (interest in global
performance predictions, short life-time,
cryogenic propellers, distribution of the active
cavitation nuclei not known)
Homogeneous-flow model, i.e. liquid/vapour
mixture modeled as a single-phase fluid
4
Cavitation model
Model proposed by L. dAgostino et al. (2001),
which takes into account (at least approximately)
for thermal cavitation effects .
Barotropic flow the state equation relates
pressure and density. In particular, for the
considered model the state equation has the
following form
liquid
characteristic fluid physical parameters (known)
liquid-vapour mixture
Starting from this context, the more general
scope of the research activity is to develop a
numerical tool for the simulation of barotropic
flows in complex geometry.
5
Cavitating flow behavior
Difficulty in the homogenous-flow description,
the physical properties of the flow change
dramatically between the zones of pure liquid and
the cavitating regions (fluid/vapour mixture).
speed of sound
high speed of sound (H2O 1400 m/s)
Mltlt1 ? incompressible regime
density
high M ? supersonic-hypersonic flow
6
Modellazione dei flussi cavitanti

• Flows characterized by nearly incompressible
  zones together with highly supersonic flow
  regions.

two possible choices
Numerical solvers for incompressible flows
suitably corrected to take into account
compressibility
Numerical discretization of the compressible flow
equations
? unknown, p from the state equation
p unknown
Examples of applications to cavitating flows van
der Heul et al., ECCOMAS 2000, Senocack and Shyy,
JCP 2002,
7
Mathematical model
Because of the barotropic state law the energy
equation is decoupled ? only mass and momentum
balances are considered

• Equations for 3D laminar viscous compressible
  flows (no turbulence) in conservative variables
• or
• Equations for 3D inviscid compressible flows in
  conservative variables

barotropic state equation (ODE or analytic laws)
8
Numerical discretization outline

• 1D inviscid flows
• Spatial discretization of 1st order of accuracy
  (preconditioning)
• Linearized implicit time advancing
• Extension to 2nd order of accuracy in space
  (MUSCL)
• Time advancing for 2nd order scheme (defect
  correction)
• 3D inviscid flows (in rotating frames)
• extension of the previous ingredients to
  tetrahedral unstructured grids
• 3D laminar viscous flows
• P1 finite-element discretization of viscous flows
  (not shown)

9
1D inviscid flowsspatial discretization
Galerkin projection on the finite-volume basis
functions (piecewise constant)
time discretization
10
Numerical fluxes
Godunov-type flux the exact solution of the
Riemann problem between two neighboring cells is
used.
In the present research activity, a procedure has
been developed for the construction of the
Riemann problem solution for Euler equations and
a generic convex barotropic law.
Reference E. Sinibaldi, Implicit preconditioned
numerical schemes for the simulation of
three-dimensional barotropic flows, Pisa,
Edizioni della Normale, in press, ISBN
978-88-7642-310-9.

• Exact 1D benchmark for generic barotropic state
  laws
• Construction of a Godunov-type scheme

11
Numerical fluxes
Roe scheme approximated solution of the Riemann
problem between two neighboring cells is used. .
centered part
upwind part
Roe matrix
numerical viscosity
Contribution of the present research activity ?
definition of the Roe matrix for a generic
barotropic state law (PhD. Thesis by E. Sinibaldi
or Sinibaldi, Beux Salvetti, INRIA RR4891, 2003
(available on line), Sinibaldi, Beux Salvetti,
Flow Turbulence and Combustion 76(4), 2006).
12
Preconditioning for low-Mach numbers

• Problem the numerical solvers for compressible
  flows suffer in general of accuracy problems if
  applied to low Mach flows. Following Guillard and
  Viozat (1999), an asymptotic analysis for M?0
  (Sinibaldi, Beux Salvetti, 2003 or P.H.D.
  Thesis by E. Sinibaldi) shows that
• the continous solution is characterized by
  pressure variations in space of the order of M2
• the semi-discrete solution is characterized by
  pressure variations in space of the order of M

preconditioning (following Guillard and Viozat,
1999)

• the scheme becomes accurate also for M?0
  (asymptotic analysis)
• preconditioning is applied only to the upwind
  part ? time consistency for unsteady problems

The preconditioning matrix P is of Turkel-type
(for its expression see Sinibaldi, Beux
Salvetti, Flow Turb. Comb. 2006 or P.H.D. Thesis
by E. Sinibaldi).
13
Time discretization for 1st order schemes
Adopted approach implicit linearized scheme

• Backward Euler implicit scheme

• We have shown that for the Roe scheme(Sinibaldi
  et al., 2003 and Sinibaldi P.h.D. Thesis)

• Thus the implicit scheme can be linearized as
  follows

linear system (tridiagonal in 1D)
NB remark that we did not use the homogeneity of
the Eulerian fluxes, which does not hold for
generic barotropic state laws.
14
Space discretizationextension to 2nd order of
accuracy
Adopted approach MUSCL reconstruction
Wij and Wji are the extrapolated values of the
variables at the cell interface
Gradients can be computed in different ways, by
combining different approximations (limited
stencil ad hoc coefficients)
Wi-1,i
Wi,i-1

• different schemes
• 2nd order accurate
• introducing a numerical viscosity proportional to
  2th, 4th or 6th order derivatives (Camarri et
  al., Comp. Fluids 2004).

15
Time advancing for the 2nd order accurate scheme
Adopted approach defect correction

• implicit formulation with a BDF method of order
  q

p-accurate discretization of the spatial
differential operator

• simpler non linear systems are iteratively
  considered (for p2)

second-order
First-order

• s-th DeC iteration after linearization

first-order linearized operator
block tridiagonal linear system
16
Time advancing for the 2nd order accurate scheme
Adopted approach defect correction
Full convergence of the DeC iteration is not
needed to reach the higher order of accuracy in
space and time (Martin and Guillard, Comput.
Fluids, 1996)

• We have shown that, in our case, one DeC
  iteration is sufficient to reach 2-order (space
  and time) accuracy

block tridiagonal linear system

• For comparison, the fully 2-order linearized
  approach gives a block pentadiagonal system in 1D

17
Extension to 2D-3D
Unstructured grids (tetrahedra)

• Easy to build and to refine for 3D complex
  geometry

• With respect to structured grids
• Larger complexity of implementation of numerical
  algorithms
• Larger computational requirements for fixed
  number of nodes.

18
Extension to 2D-3D (methodology developed at
INRIA Sophia-Antipolis)
Finite-volume dual grid (cells) obtained by using
the medians of the tetrahedra faces
normal integrated on the cell boundary
Roe flux
1st order
2nd order
19
Rotating frames
Extension to non-inertial frames rotating with a
constant rotational speed ?

• Incorporation of the non-inertial terms
  (Coriolis and centrifugal effects) in a source
  term (S) in the momentum equation.
• Finitevolume discretization in space of S

with
source term at node i

• Linearized implicit time-discretization (to be
  incorporated in the scheme) through the Jacobian

diagonal term in linear system matrix
known (RHS term)
independent of time
20
Quasi-1D water flow in a nozzle
Water in standard condition ? M? 10-3
1st order of accuracy in space and Roe numerical
flux
21
Quasi-1D water flow in a nozzle
Effects of preconditioning on the solution
accuracy
Pressure distribution along the nozzle axis in
non-cavitating conditions (pure liquid)
non preconditioned

• Preconditioning is actually important and works
  well in improving the numerical solution
  accuracy.

22
Quasi-1D water flow in a nozzle
Effects of preconditioning and of the time
advancing scheme on the numerical efficiency
Test-case (sample) Mach Cav./Non-cav. Time-step Expl. non-prec. Time-step Expl. prec.
TC1 3.5e-3 Non-cav. 1.0e-5 1.0e-6
TC2 3.5e-3 Cav. 1.0e-5 5.0e-7
TC3 7.0e-4 Non-cav. 1.0e-5 5.0e-7
TC4 7.0e-5 Non-cav. 1.0e-5 5.0e-8
Time-step Impl. prec.
8
1.0e-5
8
8

• For explicit time advancing, preconditiong
  significantly decreases the maximum allowable
  time step. This reduction becomes more important
  as the Mach number decreases ? ?precO(M)?noprec
  (see E. Sinibaldi PhD. Thesis)
• The preconditioned linearized implicit scheme
  has practically no time step limitation in non
  cavitating conditions.
• In cavitating conditions, some improvements with
  respect to the explicit scheme are found, but
  severe limitations on the time-step remain.

23
Riemann problems

• 1D numerical experiments for Riemann problems
  characterized by
• different barotropic laws (including the one for
  cavitating flows)
• different characteristic waves
• different regimes (low Mach, transonic)
• Roe and Godunov fluxes (1st order of accuracy)
• implicit linearized scheme

• No differences between the results obtained with
  the Godunov scheme and with the Roe one ? the
  Godunov scheme will not be used in other
  applications (2D, 3D) because much more
  computationally demanding
• For non-cavitating barotropic laws the results
  show
• accuracy consistent with the used 1st order
  accurate approximation,
• satisfactory efficiency of time advancing.

24
Riemann problems

• flow regime low Mach

• flow regime generic Mach

• barotropic law

• barotropic law

• 2 shocks

• shock and rarefaction

p
?
600 cells
4000 cells
u
u
?
stationary contact
blow-up for c(CFL) ?10
blow-up for c(CFL) ?100
25
Riemann problem for the cavitation barotropic law

• flow regime low Mach / high Mach

4000 cells

• barotropic law LdA model for cavitation

2 rarefactions
head
p
u
pressure
tail
pressure (detail)
p
p
barotropic curve, for reference
tail !!!
?
26
Riemann problem for the cavitation barotropic law

• Very fine spatial discretization and small time
  steps are needed to capture pressure and density
  spikes in the cavitating region

4000 cells
27
Water flow around a hydrofoil mounted in a tunnel
(Beux et al.,M2AN, 2005)
Dirichlet
homog. Neumann

• Inviscid flow.
• 1st order of accuracy and Roe scheme
• Linearized implicit time advancing

Free-streams (T 293.16 K)
cavitation number
non-cavitating
cavitating
Grids
GR2
GR1 (det.)
cells
tetrahedra
28
test-section
Water flow around a hydrofoil mounted in a tunnel
(Beux et al.,M2AN, 2005)
Centro Spazio, Pisa
Pressure distribution over the hydrofoil
local preconditioning only in the cavitating
region

• Surprisingly good accuracy.
• Problems of efficiency non-cavitating
  simulations CFL up to 400, with cavitation CFLmax
  10-2

29
Water flow around a hydrofoil mounted in a tunnel
(Beux et al.,M2AN, 2005)
local cavitation number
Mach
sigma
Mach up to 28
? 0.1
less pronounced Mach variation, OK
more extended cavity, OK
Mach up to 11
? 0.01
30
Some Remarks

• First series of test-cases (inviscid flows, 1st
  order of accuracy in space, preconditioning,
  linearized implicit time advancing)
• quasi-1D water flow in a nozzle (non cavitating
  and cavitating conditions)
• Riemann problems with different barotropic state
  laws
• water flow around a hydrofoil (non cavitating
  and cavitating conditions)
• water flow in a turbopump inducer in
  non-cavitating conditions (not shown)

• Satisfactory accuracy (in the limit of the
  assumptions made) in both non-cavitating and
  cavitating conditions.
• Numerical efficiency problems when cavitating
  regions are present.

• Additional series of 1D numerical experiments
• to investigate whether the efficiency problems
  in cavitating conditions are due to the adopted
  linearization technique for time advancing
• to test the efficiency of the defect correction
  approach for 2nd order accuracy simulations in
  non-cavitating conditions

31
Additional series of 1D numerical experiments
linearized implicit vs. fully non-linear implicit
in cavitating conditions
Test-case Riemann problem

• 2 initial liquid states ? 2 rarefactions
• LdA cavitating flow state equation

Solution accuracy
Pressure field at t1s
Discretization
4000 cells
-100
100
2000
-2000
detail
Robustness and computational cost

• No improvement in robustness with the fully
  implicit formulation ? as for non cavitating
  flows, the fully implicit simulations blow up at
  lower CFL than the ones with the linearized
  implicit scheme.
• For the same resolution in space and the same
  time step, the computational costs are much
  larger for the fully implicit scheme.

32
Validation of 2-order formulation 1D numerical
experiments
TEST-CASE 1 Quasi-1D water flow in a
convergent-divergent nozzle

• flow regime steady and supersonic
• barotropic law

FO
Density field
Spatial discretization 400 cells
Comparison of implicit formulations FO
first-order (in space and time) linearized
implicit slope
1 2-order (in space and time) linearized
implicit DeC Defect correction approach
DeC with 1 and 2 inner iterations SO
fully second-order linearized implicit FU
fully implicit second-order
(non linear solver (PETSc
library) based on a gradient-free Newton-GMRES
approach)
Velocity error
log-log scale
slope 2
33
TEST-CASE 2 Riemann problem (shock and
rarefaction)

• flow regime unsteady and subsonic
  barotropic law

velocity field (t1 s)
Spatial refinement
Temporal discretization ?t0.0001
400 cells
40 cells
temporal refinement
Spatial discretization 400 cells
DeC2, DeC3
DeC1
?t0.01
?t0.001
34
Additional series of 1D numerical experiments
Validation of the second-order formulation
Solution accuracy

• No loss of accuracy with the present
  formulation
• neither due to the defect correction ?comparison
  DeC/fully second-order linearized implicit
• nor due to the linearization of the implicit
  time-advancing ?comparison linearized
  implicit/fully implicit formulations
• In accordance with the theoretical appraisal,
  one iteration of defect correction is already
  sufficient to reach 2nd order accuracy
• Nevertheless, for particular cases (large CFL
  number), a second inner iteration can improve the
  solution (stabilization effect)

Computational cost

• Steady regimes the steady solution is obtained
  after very few pseudo-time iterations for all the
  linearized implicit approaches while the fully
  implicit formulation needs a CFL-like condition
  (for fine spatial discretization, i.e. large
  dimension of the non linear system)).
• For the same grid and time step, DeC1 is
  approximately two times cheaper than the fully
  second-order linearized implicit approach. A
  larger ratio is expected for 3D cases due to the
  increase of complexity and stiffness.

35
Concluding Remarks and Developments

• For non-cavitating barotropic flows, the
  proposed numerical methodology shows
  satisfactory
• accuracy (MUSCL reconstruction preconditioning
  for low Mach)
• robustness and efficiency (linearized implicit
  time advancing defect correction)
• For cavitating flows and the homogeneous flow
  model
• severe restriction of the time step are observed
  ? unaffordable CPU requirements for 3D
  simulations
• numerical experiments show that this is not due
  to the adpted linearization of the implicit time
  advancing

• For cavitating flows described through the
  homogenous-flow model
• try more robust numerical fluxes (HLL, HLLC)
  and/or
• relaxion techniques in time
• Change cavitation model (two-phases)

Application of the numerical set-up (as it is) to
the simulation of problems characterized by
barotropic laws less stiff than the cavitating
one (shallow water, atmosphere)
36
Cavitating flow behavior
pressure
liquid incompressible
liquid-vapour mixture highly compressible
density
37
Time advancing for the 2nd order accurate scheme
Full second-order linearized approach

38
Flusso di acqua in un induttore di
turbopompa (Sinibaldi et al., 2006)
inter-blade covering no gap
very complex geomety (detail of hub-blade
intersection)
Free-stream (T 296.16 K)
2.5x106 elements
39
Flusso di acqua in un induttore di
turbopompa (Sinibaldi et al., 2006)
pressure contours
velocity (longitudinal cut plane)
max (red) 177700 Pa min (blue) 79700
Pa spacing 5000 Pa
axial back-flow correctly described!
pretty nice results it seems a promising
scheme! (cavitating simulation not affordable -at
a reasonable cost- due to the efficiency issue)
40
Homogeneous-flow modelsThermal barotropic model
(dAgostino et al., 2001)
pure liquid weakly compressible fluid
mixture
In which ?L, g, pc, ?V are constants dependent
on the considered flow