Transcript Simulazione di flussi turbolenti
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 Sant’Anna, Pisa)
Micro-Macro Modelling and Simulation of Liquid-Vapour Flows Strasbourg, 23-25 January 2008
Industrial and engineering motivation
Development of a numerical tool performance of axial inducers typical of turbopumps for liquid propellent rockets. for the prediction of the
rotating inducer non rotating cylindrical case The main role of the inducer is to increase the fluid pressure decrease) through rotation. (velocity Fundings from the Italian Space Agency and European Space Agency.
Difficulties
Complex rotating 3D geometry.
Severe size limitations which lead to very high rotational speed cavitation phenomena need of a cavitation.
model to take into account INFLOW 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
Cavitation model
Model proposed by L. d’Agostino 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 liquid-vapour mixture
p
p sat p
f p
p sat d
dp
characteristic fluid physical parameters (known) 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.
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).
high M supersonic-hypersonic flow density high speed of sound (H 2 O ~ 1400 m/s) M<<1 regime incompressible
Modellazione dei flussi cavitanti
Flows characterized by nearly incompressible highly supersonic flow regions.
zones together with
two possible choices:
Numerical solvers for incompressible flows suitably corrected to take into account compressibility p unknown Examples of applications to cavitating flows: van der Heul et al., ECCOMAS 2000, Senocack and Shyy, JCP 2002, … Numerical discretization of the compressible flow equations unknown, p from the state equation
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 in conservative variables compressible flows ( no turbulence ) or Equations for 3D inviscid compressible flows in conservative variables + barotropic state equation (ODE or analytic laws)
Numerical discretization: outline
1D inviscid flows Spatial discretization of 1 st order of accuracy (preconditioning) Linearized implicit time advancing Extension to 2 nd order of accuracy in space (MUSCL) Time advancing for 2 nd 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)
1D inviscid flows:spatial discretization
d dt
t W
x
i-1,i
C i
W i-1 C i-1 ) 0
W
i,i+1
F
u
u
2
p
p
p
Galerkin projection on the finite-volume basis functions (piecewise constant) W i W i+1 C i
i
C i+1
C i
C i
1
F
C i
1
C i
0 time discretization x
Numerical fluxes
i dW i dt
i
1,
i
0
Numerical fluxes
Godunov-type flux : the exact solution of the Riemann problem between two neighboring cells is used.
RP
i L i R
1 ,0
i
1,
i
RP
W i L
1 ,
W i R
,0 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
Numerical fluxes
Roe scheme : approximated solution of the Riemann problem between two neighboring cells is used. . (
l
,
r
)
l
) 2 centered part
r
) 1 2
l
,
r
) (
W l
upwind part
W r
) Roe matrix numerical viscosity Contribution of the present research activity Roe matrix for a 2006).
generic barotropic state law definition of the (PhD. Thesis by E. Sinibaldi or Sinibaldi, Beux & Salvetti, INRIA RR4891, 2003 (available on line), Sinibaldi, Beux & Salvetti, Flow Turbulence and Combustion 76(4),
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 M 2 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 unsteady problems time consistency for 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).
Time discretization for 1
st
order schemes
Adopted approach: implicit linearized scheme
Backward Euler implicit scheme:
W i n
1
t
n
1
i i
i n
1 1,
i
W i n
P.h.D. Thesis):
i
We have shown that for the Roe scheme (Sinibaldi et al., 2003 and Sinibaldi
n
1
i n
,
j n
)
W i n
1
W i n
i n
,
j n
)
W j n
1
W j n O
,
A
1 2
A
A
Thus the implicit scheme can be linearized as follows:
n W j B
1
B
1
n W i
1
B
0
n W i
B
1
n W i
1
n i n
1 ,
W i n
)
B
0
x i
n t I
i n
,
W i n
1 ) linear system (tridiagonal in 1D)
B
1 1
n i
1,
i
)
n n i
1 ,
W i i n
) ,
W i n
1 ) NB: remark that we did not use the homogeneity of the Eulerian fluxes, which does not hold for generic barotropic state laws.
Space discretization:extension to 2 nd order of accuracy
Adopted approach: MUSCL reconstruction
W ij and W ji are the extrapolated values of the variables at the cell interface
W i-1,i W i-1 W i,i-1 W i W i+1 Gradients can be computed in different ways, by combining different approximations (limited stencil + ad hoc coefficients) different schemes 2nd order accurate introducing a numerical viscosity proportional to 2 2004).
th , 4 th or 6 th order derivatives (Camarri et al., Comp. Fluids
Time advancing for the 2 nd order accurate scheme
Adopted approach: defect correction
implicit formulation with a BDF method of order q:
p h n
1 ) 0 with
L p
( W ) 1
t
a q
,0 W
i q
1
a W h n i
h
( W )
h
p-accurate discretization of the spatial differential operator simpler non linear systems are iteratively considered (for p=2):
L
1 ( W
s
1 ) W
W h n
0
L
1 1 given ( W
s
W )
M
L
2 ( W
s
)
s
W 0,..., 0 s-th DeC iteration after linearization:
M W h n
) 1 First-order 3 2
t I M h
(1)
M h
(1) ( W
s
)
s
W 3 W
s
4
W h n
2
t
W h n
1 (2)
h
( W
s
) first-order linearized operator block tridiagonal linear system second-order
Time advancing for the 2 nd 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: 3 2
t I
M
(1)
h
(
W h n
)
W h n
1
W h n
W h n
2
t W h n
1 (2)
h
(
W h n
)
O
(
x
block tridiagonal linear system:
x t B
1
n W i
1
B
0
n W i
B
1
n W i
1
S i n t
) 2
B
1
B
0
x i t
I
i n
1 ,
W i n
)
B
1
i n
,
W i n
1 )
i n
,
W i n
1 )
n i
1 ,
W i n
)
S i n
2
x i t
(
W i n
1
W i n
) For comparison, the fully 2-order linearized approach gives a block pentadiagonal system in 1D
n
1
n i
1,
i
)
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.
Extension to 2D-3D
(methodology developed at INRIA Sophia-Antipolis)
Finite-volume dual grid (cells) obtained by using the medians of the tetrahedra faces
nodes in the neighbourhood of
j
node i
C hk
C kh
C k
h i
normal integrated on the cell boundary Roe flux
C h
C h
1 st order 2nd order
k C k
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 . Finite–volume discretization in space of S:
S i
i
0
i
2
u i i
source term at node i with
i
1
i
C i xdV
Linearized implicit time-discretization (to be incorporated in the scheme) through the Jacobian:
S i n
1
S i n
S i
W i
n W i
diagonal term in linear system matrix known (RHS term) independent of time
Quasi-1D water flow in a nozzle
STEADY-STATE IN a C-D SYMMETRICAL NOZZLE:
t
u
x
u
2
u
p
A dA dx
u
u
2
A
cross-sectional area min numerical transient source 2 Symmetrical grid, 360 cells, minimum spacing 0.02 (throat) 1 IN OUT 5.7
5.0
Inlet B.C’s:
p
p
( 0 ),
u
u
0 I.C’s: 0 ,
u
u
0 Outlet B.C’s:
p
/
x
u
/
x
0 Water in standard condition M 10 -3 1 st order of accuracy in space and Roe numerical flux
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.
Quasi-1D water flow in a nozzle
Effects of preconditioning and of the time advancing scheme on the numerical efficiency
Test-case (sample) TC1 TC2 TC3 TC4 Mach 3.5e-3 3.5e-3 7.0e-4 7.0e-5 Cav./Non cav.
Non-cav.
Cav.
Non-cav.
Non-cav.
Time-step Expl. non prec.
1.0e-5 1.0e-5 1.0e-5 1.0e-5 Time-step Expl. prec.
1.0e-6 5.0e-7 5.0e-7 5.0e-8 Time-step Impl. prec.
∞ 1.0e-5 ∞ ∞ For explicit time advancing maximum allowable time step. This reduction becomes more important as the Mach number decreases , preconditiong significantly decreases the prec =O(M)* noprec (see E. Sinibaldi PhD. Thesis) The preconditioned linearized implicit scheme limitation in non cavitating conditions .
has practically no time step In cavitating conditions , some improvements with respect to the explicit scheme are found, but severe limitations on the time-step remain.
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 (1 st 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 1 st order accurate approximation, satisfactory efficiency of time advancing.
flow regime: low Mach barotropic law:
6
2
shock and rarefaction
Riemann problems
flow regime: generic Mach barotropic law:
2 shocks
p
4000 cells
blow-up for c (CFL)
100
600 cells
stationary contact
blow-up for c (CFL)
10
Riemann problem for the cavitation barotropic law
flow regime: low Mach / high Mach barotropic law: LdA model for cavitation
4000 cells 2 rarefactions
head
u p pressure
tail
pressure (detail) p’ p barotropic curve, for reference ’
tail !!!
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
head
4000 cells p pressure, for reference density
tail
p’ barotropic curve, for reference ’ density (detail)
Water flow around a hydrofoil mounted in a tunnel
(Beux et al.,M 2 AN, 2005)
IN Dirichlet Free-streams (T = 293.16 K) : Grids: OUT homog. Neumann cavitation number non-cavitating cavitating Inviscid flow.
1 st order of accuracy and Roe scheme Linearized implicit time advancing GR1 (det.) GR2 cells tetrahedra
Water flow around a hydrofoil mounted in a tunnel
(Beux et al.,M 2 AN, 2005)
Pressure distribution over the hydrofoil test-section
Centro Spazio
, Pisa effect of = 0.1
= 0.01
almost independent of the grid local preconditioning only in the cavitating region Surprisingly good accuracy.
Problems of efficiency cavitating simulations CFL up to 400, with cavitation CFL max : non = 10 -2
Water flow around a hydrofoil mounted in a tunnel
(Beux et al.,M 2 AN, 2005)
local cavitation number Mach Mach up to 28 sigma = 0.1
more extended cavity, OK = 0.01
less pronounced Mach variation, OK Mach up to 11
Some Remarks
First series of test-cases (inviscid flows, 1 st 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 due to the adopted linearization technique in cavitating conditions are for time advancing to test the efficiency accuracy simulations of the defect correction in non-cavitating conditions approach for 2 nd order
-2000 Test-case: Riemann problem 2 initial liquid states 2 rarefactions LdA cavitating flow state equation Solution accuracy Pressure field at t=1s Discretization: 4000 cells
t
5.10
5 2000 -100 Robustness and computational cost detail 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.
100
FO
Validation of 2-order formulation: 1D numerical experiments
TEST-CASE 1: Quasi-1D water flow in a convergent-divergent nozzle Density field flow regime: steady and supersonic
barotropic law: Spatial discretization: 400 cells Velocity error log-log scale 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) slope 2
TEST-CASE 2: Riemann problem (shock and rarefaction)
velocity field (t=1 s) 6 Spatial refinement Temporal discretization : t=0.0001
40 cells 400 cells DeC2, DeC3 temporal refinement Spatial discretization: 400 cells t=0.01 DeC1 t=0.001
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 order linearized implicit comparison DeC/fully second 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 2 nd 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.
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 for 3D simulations unaffordable CPU requirements numerical experiments show that this is not due to the adpted linearization of the implicit time advancing 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…) 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)
Cavitating flow behavior
liquid-vapour mixture: highly compressible liquid: ~ incompressible density
Time advancing for the 2 nd order accurate scheme
Full second-order linearized approach
B
2
n W i
2
B
1
n W i
1
B
0
n W i
B
1
n W i
1
B
2
n W i
2
S i n
B B B B B
2 1 0 1 2 3 2 2 3
h
2 1 2
h
i
h i i h i
1
x
2
h t h i
2
h i
1
A
A
2 1
A i
1,
i
1,
i I i
1
i A i
1,
i
1
A A A n n
i
n
A
1 1 1
n
1
n
1
n n
2 2
i
n
1 2,
i
i
1, 1 1
i n n i
1,
i
1
n
A i
1,
i i
1,
i
2
n
1
n
1
n n
n
2
A
2 2
i
1, 2 1
h i h h
h i
1
i A i
i
1 2
A
2
i
1, 1
h
i i i h
i
i
2,
i
1 1
i
1,
i
i
2
A n
i
1,
i n
1
n
1
n n
A
+
i
1,
i A i
1,
i i
1,
i A
n
1
n n
i
1,
i n
h i h i
1
n
A i
1,
i
A
1
n
1 1
n
n
Flusso di acqua in un induttore di turbopompa
(Sinibaldi et al., 2006)
inducer inter-blade covering: no gap nose afterbody
very complex geomety
(detail of hub-blade intersection) 2.5
x 10 6 elements Free-stream (T = 296.16 K) :
Flusso di acqua in un induttore di turbopompa
(Sinibaldi et al., 2006) velocity (longitudinal cut plane) pressure contours
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)
Homogeneous-flow models Thermal barotropic model
(d’Agostino et al., 2001) pure liquid: weakly compressible fluid
p
p sat
1 ln
sat p
p sat
mixture
p
p sat d
dp
1
p
( 1 ) ( 1
L
)
L p
L g
p c p
V
l
In which
L , g*, p c ,
V
l
;
T R
are constants dependent on the considered flow since
( 1 )
lsat
free model parameter
d
dp
1
a
2
f
( ,
p
( ))