- Email: [email protected]
A multilevel local mesh refinement projection method for low Mach number flows
Mathematics and Computers in Simulation 61 (2003) 477–488
A multilevel local mesh refinement projection method for low Mach number flows Xavier Coré a , Philippe Angot b , Jean-Claude Latché a,∗ a
b
Département de Recherches en Sécurité, Institut de Protection et de Sˆureté Nucléaire, F-13108 Saint Paul Lez Durance Cedex, France LATP (UMR 6632 CNRS), IMT, Technopˆole de Chˆateau-Gombert, 38 rue F. Joliot Curie, F-13451 Marseille Cedex 20, France
Abstract This paper addresses a sub-problem of low Mach number reactive flows modelling: the solution of the mass and momentum balance equations for a pressure-independent variable density flow. For this objective, we develop a novel numerical scheme: the time discretisation is performed by an incremental projection method, using the original form of the mass balance equation for the projection step; the resulting two semi-discrete problems are then discretised in space by finite volumes on a MAC staggered mesh, with a multilevel local mesh refinement algorithm based on the Flux Interface Correction (FIC) method. The efficiency of this numerical scheme, which offers the opportunity to avoid computing with non-uniform grids or unstructured meshes, is assessed against a practical test case, specifically produced for this purpose. © 2002 IMACS. Published by Elsevier Science B.V. All rights reserved. Keywords: FIC method; Finite volume scheme; Low Mach number flows; Multilevel local mesh refinement; Projection method; Variable density flows
1. Introduction The derivation of the equations for zero Mach number combustion was done in the 1980s by Majda and Sethian [13]. They obtained a system of coupled balance equations, expressing the balance of mass, momentum, energy and chemical species, which, from the point of view of the flow dynamics, shows some similarities with the classical incompressible Navier–Stokes equations: the flow density is independent of the dynamic pressure and the velocity field satisfies a divergence constraint, similar to the usual divergence-free condition, although more complex. So far, using a numerical method originated from the class of incompressible flow solvers appears to be a natural choice for the solution of this problem. ∗
Corresponding author. E-mail address: [email protected] (J.-C. Latch´e). 0378-4754/02/$ – see front matter © 2002 IMACS. Published by Elsevier Science B.V. All rights reserved. PII: S 0 3 7 8 - 4 7 5 4 ( 0 2 ) 0 0 1 4 0 - 4
478
X. Cor´e et al. / Mathematics and Computers in Simulation 61 (2003) 477–488
Indeed, various numerical schemes derived from the projection method initially proposed by Chorin [4] have up to now been implemented; a review can be found in [10]. Most of the proposed numerical schemes use, for the projection step, an explicit expression of the velocity field divergence, obtained by combining the mass balance equation with the energy and chemical species balance equations and the fluid equation of state. The method developed here follows a different path: we use for the projection step the original form of the mass balance equation, i.e. (∂ρ/∂t) + ∇ · (ρ v) = 0. This approach leads to a quite general problem, namely designing a projection method for variable density flows, since the temporal and spatial variations of the density may result from a wide range of physical phenomena. The present paper focuses on this point, the solution of the complete combustion problem being the aim of a further publication. To our knowledge, the only previous attempt to use such a method was in the work of Najm et al. [6]: a reacting flow problem was solved by a projection scheme formally of second-order in time, using a standard finite difference technique for the spatial discretisation. In the present work, we are mainly concerned with the spatial discretisation of the equations: starting from a finite volume approach using a MAC-type staggered mesh, e.g. [14], we develop a local mesh refinement method, using the Flux Interface Correction (FIC) algorithm introduced by Angot and coworkers [1,2,9]. The time semi-discretisation of the equations is performed by a first-order incremental projection method. Special care has been taken to ensure the discrete mass conservation. The paper is organised as follows: first we briefly recall the equations governing the flow dynamics within the low Mach number framework; then the discretisation of the equations is described; the fourth section is devoted to the application of the FIC methodology to the problem in hand; finally, the efficiency of the method is tested against a specifically designed benchmark.
2. Low Mach number balance equations Assuming that the temperature T ( x , t) is given a priori, the so-called low Mach number balance equations for a Newtonian flow such that γ M 2 1 read in a two-dimensional computational domain Ω: Find ( v , P1 ) from Ω × [0; T ] to, respectively, R2 and R such that ∂ρ ∂t + ∇ · (ρ v) = 0,
1 2 ∂(ρ v) T + ∇ · (ρ vv) + ∇P1 = ∇ · µ(∇ v + ∇ v ) − µ(∇ · v)Id , ∂t γ 3 where v is the velocity field, P1 stands for the dynamic part of the pressure, the dynamic viscosity µ = µ(ρ) > 0 and the specific heat ratio γ > 0 are given. In addition, if the physical domain is connected with the atmosphere, the density ρ is known a priori: ρ = ρ( x , t) ≥ ρ0 > 0. This system must be supplemented with initial conditions ( v ( x , 0) = v0 ) and boundary conditions: ∂Ω is supposed to be the union of a finite number of parts, each of them being associated with a boundary condition of an a priori specified type: inflow (ρ v = (ρ v)q ), outflow (P1 = 0, v · t = 0) or symmetry (∇( v · t) · n = 0, v · n = 0), n and t denoting respectively the outward normal and tangential unit vectors on the boundary.
X. Cor´e et al. / Mathematics and Computers in Simulation 61 (2003) 477–488
479
3. Semi-implicit projection method and finite volumes scheme The developed method belongs to the class of so-called fractional step schemes. In a prediction step, the velocity field is evaluated by solving the momentum balance equations, by an Euler time integration scheme. An explicit evaluation of the gradient of the dynamic pressure is kept in the semi-discrete equations, according to the so-called incremental form of the projection method. The convective fluxes are also discretised explicitly in time whereas the viscous fluxes are treated implicitly. Consequently, this step reduces to a system of linear advection–diffusion equations. Moreover, the explicit discretisation of convective fluxes allows an easy implementation of a non-oscillatory strategy. As a counterpart, a CFL-type stability condition must be satisfied. The correction step carries out the projection of the predicted velocity onto the space of velocity fields which satisfy the mass balance equation. It involves the solution of a Poisson equation to obtain the corresponding dynamic pressure. The time semi-discretisation is detailed hereafter from tn = n t. (i) AD: advection–diffusion of velocity field 1 n+1 n+1 1 (ρ v˜ − ρ n vn ) + ∇ · (ρ n vn vn ) + ∇P1n , t γ 2 = ∇ · (µn+1 (∇ v˜ n+1 + (∇ v˜ n+1 )T ) − µn+1 (∇ · v˜ n+1 )Id), in Ω, 3 ˜ q ; outflow : v˜ n+1 · t = 0, inflow : ρ n+1 v˜ n+1 = (ρ v) symmetry : ∇(v˜ n+1 · t) · n = 0, v˜ n+1 · n = 0, on ∂Ω. (ii) PC: dynamic pressure estimation and correction of velocity 1 n+1 n+1 1 (ρ v − ρ n+1 v˜ n+1 ) + (∇P1n+1 − ∇P1n ) = 0, t γ −1 n+1 ∇ · (ρ n+1 vn+1 ) = (ρ − ρ n ), in Ω, t inflow : ρ n+1 vn+1 · n = ρ n+1 v˜ n+1 · n; outflow : P1n+1 = 0, symmetry : ∇P1n+1 · n = 0, on ∂Ω. For the space discretisation, a finite volume MAC-type staggered grid is used. The convective fluxes are discretised using a third-order upwind scheme, the QUICK version proposed by Hayase et al. [7], with the ULTRA slope limiter, developed in [12]. We therefore obtain a positivity and monotonicity preserving method as shown in [8]. Finally, a second-order centred scheme is used to discretise the diffusive fluxes, e.g. [14]. The linear systems are solved by MILU(0) preconditioned conjugate gradients with a stopping criterium on the residuals smaller than = 10−8 .
480
X. Cor´e et al. / Mathematics and Computers in Simulation 61 (2003) 477–488
4. Multilevel local mesh refinement method 4.1. The Flux Interface Correction method: general algorithm The Flux Interface Correction method was initially proposed by Angot and coworkers [1,2,9,11]. It is a multilevel mesh refinement technique, and as such it involves two (or more) levels of discretisation: a mesh GH covering the whole computational domain and one or several finer meshes Gh covering subdomains; for the sake of simplicity, we assume in the sequel that there is only one such refined zone. A differential problem is solved on each computational grid, which is basically the “one-grid discrete problem” described in the preceding section; they are referred to in the following as the coarse and fine grid problems, solved, respectively, in the computational domain Ω H and Ω h to obtain the coarse and fine grid solutions. This multigrid technique is relatively cheap since the iterative solvers have to deal only with a relatively small number of degrees of freedom on each discretisation level. The method relies on the following interactions between the coarse and fine grid problems: • the boundary conditions applied for the fine grid problem at the inner boundaries of Ω h (i.e. ∂Ω h \∂Ω H ) are derived by an interpolation technique from the coarse grid solution via the so-called prolongation operator PHh ; • in the coarse grid problem, some corrective terms appear in the right-hand side, evaluated from the fine grid solution by a correction operator RhH . For each time step, solutions on each discretisation level are iteratively computed an arbitrary number of times: typically, two or three iterations are required when quasi-exact solve of the sub-problems is performed, e.g. [2,11]. At the end of the time step, the coarse grid solution is overloaded on Ω h by a quantity evaluated from the fine grid solution via the so-called restriction operator R˜ hH . The two-level FIC iterative solution procedure with the global coarse grid GH and a local fine grid Gh can be decomposed sequentially as shown in Fig. 1, where AD stands for the first step of the incremental projection method (advection– diffusion of velocity), PC is the second step (dynamic pressure Poisson equation and correction of velocity) and ADC is the advection–diffusion of the coarse grid velocity with additional flux correction residuals derived from the fine grid solution. The grid operators PHh , RhH , R˜ hH are detailed in the next section. 4.2. The Flux Interface Correction method: grid operators A correction is performed only in the prediction step of the coarse grid problem, since the projection step consists in a self-defined mathematical operation. The correction operator is built as proposed in the
Fig. 1. Two-level FIC iterative solution procedure for k = 2 iterations.
X. Cor´e et al. / Mathematics and Computers in Simulation 61 (2003) 477–488
481
original papers, and is recalled here only for the sake of completeness. On the other hand, the prolongation and restriction operators have been developed specifically for the present application. As proposed in [2,3], we choose h = H /3 for the mesh steps to simplify the calculations for the grid operators. 4.2.1. Notations Each domain Ω l for l ∈ {H ; h}, with boundary ∂Ω l is partitioned in a finite number of control volumes Kl . The vector nl stands for the unit exterior normal defined on the boundary ∂Kl of each control volume Kl whereas ∂Kl = ∪γl for all faces γ l of Kl . We denote by ωH the subdomain of Ω h obtained by the union of the coarse control volumes that are included in the interior of the refinement area. The velocity k vH (resp. vhk ) stands for the solution of the coarse (resp. fine) grid problem obtained at the FIC iteration n+1 n and vhn (resp. vH and vhn+1 ) denote these solutions at the beginning (resp. end) of the current time k; vH step. The same convention holds for the representations of ρ used for each discrete problem. Finally, for all these fields, we use the notation (ψ H )i,j to refer to the degree of freedom of any discrete field ψ H at ¯ j¯ i, ¯ zj + j¯h) in Ω h for i, ¯ j¯ ∈ {−1; 0 + 1}. any point (xi , zj ) in Ω H and (ψh )i,j for ψ h at the points (xi + ih, 4.2.2. The flux correction operator RhH If the tensor F stands for the convective and diffusive fluxes of the momentum balance equation, i.e. 2 F( v ) = ρ vv − µ(∇ v + ∇ vT ) + µ(∇ · v)Id, 3 we define the flux correction residual vector on each coarse control volume KH ⊂ ωH by: 1 k k rKk+1 ( v ) = F ( v ) · n − F ( v ) · n , k ≥ 0. h H h H H meas(KH ) γ ⊂∂K γh ∂KH h
H
The correction terms δ fHk added to the right-hand side of the momentum equation of the coarse grid problem are then defined by the following recurrence, where |1|KH denotes the characteristic function of the control volume KH : δ fH0 ( v ) = 0,
δ fk+1 ( v ) = δ fHk ( v ) + KH ⊂ωH |1|KH rKk+1 ( v ), k ≥ 0. H H 4.2.3. The prolongation operator PHh Boundary conditions for the fine grid problem must be prescribed on the inner boundaries of Ω h for both steps of the projection method: for the advection–diffusion step, we impose Dirichlet boundary conditions while homogeneous Neumann boundary conditions are prescribed to the pressure time increment for the Poisson problem of the correction step. As an example, in Fig. 2, we define the prolongation PHh operator j¯ j¯ at the fine grid points (xi+(1/2) , zj ), j¯ ∈ {−1; 0; +1}, where zj = zj + j¯h for the first component of the velocity field by: 0,j¯
(uh )i+(1/2),j = PHh (uH ) :=
(ρH )i+(1/2),j 0,j¯
(ρh )i+(1/2),j
(uH )i+(1/2),j .
482
X. Cor´e et al. / Mathematics and Computers in Simulation 61 (2003) 477–488
Fig. 2. Uniform 2D staggered MAC grids with h = H /3 showing the nodes involved in the computation of the grid operators.
4.2.4. The restriction operator R˜ hH In the restriction step, we overload, for example the first component of the velocity at the coarse grid j¯ point (xi +(1/2) , zj ) located at the center of the control volume vertical face γH = ∪ γh for j¯ ∈ {−1; 0; +1} by:
0,j¯ 0,j¯ j¯ j¯∈{−1;0;+1} meas(γh )(ρh )i+(1/2),j (uh )i+(1/2),j H ˜ (uH )i+(1/2),j = Rh (uh ) := . meas(γH )(ρH )i+(1/2),j 4.2.5. About some coherence constraints for the grid operators Let us suppose that we can choose a representation of the density such that the mass quantity is conserved through the grids for each coarse control volume, i.e. ρH = ρh . KH
Kh ⊂KH Kh
Since the grid operators PHh and R˜ hH are indeed specially designed to satisfy the discrete equation: ρH vH · nH = ρh vh · nh , γH
γh ⊂γH γh
the definition of the grid operators ensures that the mass flux through a face of a coarse control volume is conserved by the prolongation and restriction operators and, consequently, that the following relations are verified: 1 n+1 n (ρ − ρ ) + ρhn+1 vhn+1 · nh = 0, h Ω t h Ωh h 1 n+1 n+1 (ρH − ρHn ) + ρHn+1 vH · nH = 0. KH t ∂KH The first relation is a necessary condition to ensure that the Poisson problem of the correction step on the fine grid level has a solution, whenever Neumann conditions are prescribed on the whole boundary of Ω h (situation encountered, for instance, when the refined zone Ω h is strictly interior to Ω H ). The second relation expresses that the mass balance checked by the coarse grid solution still holds after the
X. Cor´e et al. / Mathematics and Computers in Simulation 61 (2003) 477–488
483
restriction. This feature is highly desirable when, for example, the transport of a passive scalar must be treated in addition with the hydrodynamic problem.
5. Numerical results In this section, we present some numerical tests of the discretisation methods presented in this paper. Our objective is two-fold: firstly, to check the time-accuracy of the projection method; secondly, to assess the spatial accuracy gained when using the multilevel local mesh refinement technique. 5.1. Test case For this purpose, we have derived an unsteady test case from a benchmark steady problem proposed by De Lange and De Goey [5]. The computational domain is the square [0; 0.1 m] × [0; 0.1 m] as shown in Fig. 3. The temperature T ( x , t) is given in time and space, and it schematically mimics the conditions of a combustion problem: 2 r (t) , if r 2 (t) ≤ r02 , 300 1 + 2 exp κ 2 r − 1 0 2 T ( x , t) = 1 − r (t) 300 5 − 2 exp κ , otherwise, r02 with r 2 (t) = x 2 (t) + z2 and x(t) = 5x/2 if 0 ≤ t ≤ 1 s, x(t) = 5x/(3t − 1) if 1s ≤ t ≤ 2s, x(t) = x otherwise, κ = 25 and r0 = 0.05 m. The density is then given by ρ( x , t) = Patm /RT ( x , t) for an ideal gas. Fig. 4 shows its time evolution at point P(0.025 m; 0.025 m). The dynamic viscosity µ is prescribed by µ( x , t) = µref (T ( x , t)/1500)γ , with µref = 4.79 × 10−5 Pa s and γ = 0.77. The inflow boundary condition is given by: v( x , t) = (0.25((1 − z2 )/0.0025), 0)T m/s. The fluid is supposed to be initially at rest.
Fig. 3. Geometry and boundary conditions.
484
X. Cor´e et al. / Mathematics and Computers in Simulation 61 (2003) 477–488
Fig. 4. Time evolution of the density at point P(0.025; 0.025).
5.2. Results of one-grid computations Figs. 5 and 6 show, respectively, at steady state t = 3 s, the fields of velocity for the mesh 21 × 21 and dynamic pressure for the mesh 63 × 63. The computations are carried out for CFL ≈ 0.5. The mesh convergence has been studied, e.g. the horizontal velocity profile at steady state is represented for different meshes in space by two cut planes, inside and outside the density layer, respectively, in Figs. 7(a) and
Fig. 5. Steady-state velocity.
X. Cor´e et al. / Mathematics and Computers in Simulation 61 (2003) 477–488
Fig. 6. Steady-state dynamic pressure.
Fig. 7. Steady-state horizontal velocity profiles for meshes: (- - -) 21 × 21; (· · · ) 63 × 63; (—) 189 × 189 (reference).
485
486
X. Cor´e et al. / Mathematics and Computers in Simulation 61 (2003) 477–488
Fig. 8. Horizontal velocity time evolution at point P(0.025; 0.025) for the mesh 63 × 63 and time steps: (· · · ) t = 10−3 s; (- - -) t = 5 × 10−4 s.
(b). As expected, the mesh 21 × 21 gives rough results in both figures with relative errors of about 10% due to a poor representation of the density gradients, whereas the solutions for the mesh 63 × 63 and the reference mesh 189 × 189 are superposed. As shown in Fig. 8 for the time evolution of the horizontal velocity at point P(0.025 m; 0.025 m), the time step convergence is achieved for the mesh 63 × 63 at least with a time step t = 10−3 s. 5.3. Results of two-grid computations As an example, Fig. 9 shows the coarse grid GH and the fine grid Gh with h = H /3 covering the zoom area which includes the density layer. In Fig. 10, the steady-state horizontal velocity profiles at x = 0.025 m and x = 0.075 m, i.e. inside and outside the zoom area, respectively, are represented for the different meshes. At both locations, the coarse grid solution of the refined mesh 21 × 21 is better than
Fig. 9. Zoom area: coarse grid GH and fine grid Gh where h = H /3.
X. Cor´e et al. / Mathematics and Computers in Simulation 61 (2003) 477–488
487
Fig. 10. Steady-state horizontal velocity profiles for meshes: (- - -) coarse grid solution of refined mesh 21 × 21; (—) 189 × 189; (· · · ) fine grid solution of refined mesh 21 × 21 (only for x = 0.025 m in the zoom area).
the solution with no zoom on the mesh 21 × 21. Inside the refinement zone, at x = 0.025 m, the fine grid solution of the refined mesh 21 × 21 and the reference solution are superimposed. Moreover, outside the refinement zone in Fig. 10, the coarse grid solution of the refined mesh 21 × 21 is superposed on the reference solution. Hence, the multilevel zoom method greatly improves the solution precision not only inside but also outside the refinement area.
6. Conclusion This study proves the efficiency of the designed solver to calculate variable density flows with the multilevel local mesh refinement method FIC associated with the semi-implicit incremental projection. Besides, it shows how the use of such multilevel zoom solvers can be a cheap, efficient and precise enough alternative to more standard computations with non-uniform grids or unstructured meshes.
Acknowledgements This work was supported by the Institut de Protection et de Sˆureté Nucléaire (IPSN), in the framework of its modelling activities on fire consequences in nuclear power plants. Practical implementations were performed using the object oriented framework and software components library PELICANS developed at IPSN.
488
X. Cor´e et al. / Mathematics and Computers in Simulation 61 (2003) 477–488
References [1] P. Angot, J.P. Caltagirone, K. Khadra, Une méthode adaptative de raffinement local: la Correction du Flux à l’Interface, Comptes-Rendus de l’Académie des Sciences de Paris 315 (1) (1992) 739–745. [2] P. Angot, M. Laugier, La méthode FIC de raccordement conservatif de sous-domaines emboˆıtés pour un modèle de circulation océanique, C. R. Acad. Sci. Paris 319 (2) (1994) 993–1000. [3] J.P. Caltagirone, K. Khadra, P. Angot, Sur une méthode de raffinement local multigrille pour la résolution des équations de Navier–Stokes, C. R. Acad. Sci. Paris 320 (2-b) (1995) 295–302. [4] A.J. Chorin, Numerical solution of the Navier–Stokes equations, Math. Comput. 22 (1968) 745–762. [5] H.C. De Lange, P.H. De Goey, Numerical flow modelling in a locally refined grid, Int. J. Numer. Methods Eng. 37 (1994) 497–515. [6] H.N. Najm, P.S. Wyckoff, O.M. Knio, A semi-implicit numerical scheme for reacting flow. I. Stiff chemistry, J. Comput. Phys. 143 (1998) 381–402. [7] T. Hayase, J.A.C. Humphrey, R. Greif, A consistently formulated QUICK scheme for fast and stable convergence using finite-volume iterative calculation procedures, J. Comput. Phys. 98 (1992) 108–118. [8] H.T. Huynh, Accurate monotonic cubic interpolation, SIAM J. Numer. Anal. 30 (1) (1993) 57–100. [9] K. Khadra, P. Angot, J.P. Caltagirone, P. Morel, Concept de zoom adaptatif en architecture multigrille locale; Etude comparative des méthodes LDC, FAC et FIC, Modélisation Mathématique et Analyse Numérique 30 (1) (1996) 39–82. [10] R. Klein, Numerics in combustion, in: L. Vervisch, D. Veynante (Eds.), Introduction to Turbulent Combustion, Lecture Series in Computational Fluid Dynamics, Von Karman Institute, Belgium, March 1999. [11] M. Laugier, P. Angot, L. Mortier, Nested-grid methods for an ocean model: a comparative study, Int. J. Numer. Methods Fluids 23 (11) (1996) 1163–1195. [12] B.P. Leonard, S. Mokthari, Beyond first-order upwinding: the ULTRA-SHARP alternative for non-oscillatory steady-state simulation of convection, Int. J. Numer. Methods Eng. 30 (1990) 729–766. [13] A. Majda, J. Sethian, The derivation and numerical solution of the equations for zero Mach number combustion, Combust. Sci. Tech. 42 (1985) 185–205. [14] S.V. Patankar, Numerical Heat Transfer and Fluid Flow, Computational Methods in Mechanics and Thermal Sciences, McGraw-Hill, New York, 1980.
A multilevel local mesh refinement projection method for low Mach number flows Xavier Coré a , Philippe Angot b , Jean-Claude Latché a,∗ a
b
Département de Recherches en Sécurité, Institut de Protection et de Sˆureté Nucléaire, F-13108 Saint Paul Lez Durance Cedex, France LATP (UMR 6632 CNRS), IMT, Technopˆole de Chˆateau-Gombert, 38 rue F. Joliot Curie, F-13451 Marseille Cedex 20, France
Abstract This paper addresses a sub-problem of low Mach number reactive flows modelling: the solution of the mass and momentum balance equations for a pressure-independent variable density flow. For this objective, we develop a novel numerical scheme: the time discretisation is performed by an incremental projection method, using the original form of the mass balance equation for the projection step; the resulting two semi-discrete problems are then discretised in space by finite volumes on a MAC staggered mesh, with a multilevel local mesh refinement algorithm based on the Flux Interface Correction (FIC) method. The efficiency of this numerical scheme, which offers the opportunity to avoid computing with non-uniform grids or unstructured meshes, is assessed against a practical test case, specifically produced for this purpose. © 2002 IMACS. Published by Elsevier Science B.V. All rights reserved. Keywords: FIC method; Finite volume scheme; Low Mach number flows; Multilevel local mesh refinement; Projection method; Variable density flows
1. Introduction The derivation of the equations for zero Mach number combustion was done in the 1980s by Majda and Sethian [13]. They obtained a system of coupled balance equations, expressing the balance of mass, momentum, energy and chemical species, which, from the point of view of the flow dynamics, shows some similarities with the classical incompressible Navier–Stokes equations: the flow density is independent of the dynamic pressure and the velocity field satisfies a divergence constraint, similar to the usual divergence-free condition, although more complex. So far, using a numerical method originated from the class of incompressible flow solvers appears to be a natural choice for the solution of this problem. ∗
Corresponding author. E-mail address: [email protected] (J.-C. Latch´e). 0378-4754/02/$ – see front matter © 2002 IMACS. Published by Elsevier Science B.V. All rights reserved. PII: S 0 3 7 8 - 4 7 5 4 ( 0 2 ) 0 0 1 4 0 - 4
478
X. Cor´e et al. / Mathematics and Computers in Simulation 61 (2003) 477–488
Indeed, various numerical schemes derived from the projection method initially proposed by Chorin [4] have up to now been implemented; a review can be found in [10]. Most of the proposed numerical schemes use, for the projection step, an explicit expression of the velocity field divergence, obtained by combining the mass balance equation with the energy and chemical species balance equations and the fluid equation of state. The method developed here follows a different path: we use for the projection step the original form of the mass balance equation, i.e. (∂ρ/∂t) + ∇ · (ρ v) = 0. This approach leads to a quite general problem, namely designing a projection method for variable density flows, since the temporal and spatial variations of the density may result from a wide range of physical phenomena. The present paper focuses on this point, the solution of the complete combustion problem being the aim of a further publication. To our knowledge, the only previous attempt to use such a method was in the work of Najm et al. [6]: a reacting flow problem was solved by a projection scheme formally of second-order in time, using a standard finite difference technique for the spatial discretisation. In the present work, we are mainly concerned with the spatial discretisation of the equations: starting from a finite volume approach using a MAC-type staggered mesh, e.g. [14], we develop a local mesh refinement method, using the Flux Interface Correction (FIC) algorithm introduced by Angot and coworkers [1,2,9]. The time semi-discretisation of the equations is performed by a first-order incremental projection method. Special care has been taken to ensure the discrete mass conservation. The paper is organised as follows: first we briefly recall the equations governing the flow dynamics within the low Mach number framework; then the discretisation of the equations is described; the fourth section is devoted to the application of the FIC methodology to the problem in hand; finally, the efficiency of the method is tested against a specifically designed benchmark.
2. Low Mach number balance equations Assuming that the temperature T ( x , t) is given a priori, the so-called low Mach number balance equations for a Newtonian flow such that γ M 2 1 read in a two-dimensional computational domain Ω: Find ( v , P1 ) from Ω × [0; T ] to, respectively, R2 and R such that ∂ρ ∂t + ∇ · (ρ v) = 0,
1 2 ∂(ρ v) T + ∇ · (ρ vv) + ∇P1 = ∇ · µ(∇ v + ∇ v ) − µ(∇ · v)Id , ∂t γ 3 where v is the velocity field, P1 stands for the dynamic part of the pressure, the dynamic viscosity µ = µ(ρ) > 0 and the specific heat ratio γ > 0 are given. In addition, if the physical domain is connected with the atmosphere, the density ρ is known a priori: ρ = ρ( x , t) ≥ ρ0 > 0. This system must be supplemented with initial conditions ( v ( x , 0) = v0 ) and boundary conditions: ∂Ω is supposed to be the union of a finite number of parts, each of them being associated with a boundary condition of an a priori specified type: inflow (ρ v = (ρ v)q ), outflow (P1 = 0, v · t = 0) or symmetry (∇( v · t) · n = 0, v · n = 0), n and t denoting respectively the outward normal and tangential unit vectors on the boundary.
X. Cor´e et al. / Mathematics and Computers in Simulation 61 (2003) 477–488
479
3. Semi-implicit projection method and finite volumes scheme The developed method belongs to the class of so-called fractional step schemes. In a prediction step, the velocity field is evaluated by solving the momentum balance equations, by an Euler time integration scheme. An explicit evaluation of the gradient of the dynamic pressure is kept in the semi-discrete equations, according to the so-called incremental form of the projection method. The convective fluxes are also discretised explicitly in time whereas the viscous fluxes are treated implicitly. Consequently, this step reduces to a system of linear advection–diffusion equations. Moreover, the explicit discretisation of convective fluxes allows an easy implementation of a non-oscillatory strategy. As a counterpart, a CFL-type stability condition must be satisfied. The correction step carries out the projection of the predicted velocity onto the space of velocity fields which satisfy the mass balance equation. It involves the solution of a Poisson equation to obtain the corresponding dynamic pressure. The time semi-discretisation is detailed hereafter from tn = n t. (i) AD: advection–diffusion of velocity field 1 n+1 n+1 1 (ρ v˜ − ρ n vn ) + ∇ · (ρ n vn vn ) + ∇P1n , t γ 2 = ∇ · (µn+1 (∇ v˜ n+1 + (∇ v˜ n+1 )T ) − µn+1 (∇ · v˜ n+1 )Id), in Ω, 3 ˜ q ; outflow : v˜ n+1 · t = 0, inflow : ρ n+1 v˜ n+1 = (ρ v) symmetry : ∇(v˜ n+1 · t) · n = 0, v˜ n+1 · n = 0, on ∂Ω. (ii) PC: dynamic pressure estimation and correction of velocity 1 n+1 n+1 1 (ρ v − ρ n+1 v˜ n+1 ) + (∇P1n+1 − ∇P1n ) = 0, t γ −1 n+1 ∇ · (ρ n+1 vn+1 ) = (ρ − ρ n ), in Ω, t inflow : ρ n+1 vn+1 · n = ρ n+1 v˜ n+1 · n; outflow : P1n+1 = 0, symmetry : ∇P1n+1 · n = 0, on ∂Ω. For the space discretisation, a finite volume MAC-type staggered grid is used. The convective fluxes are discretised using a third-order upwind scheme, the QUICK version proposed by Hayase et al. [7], with the ULTRA slope limiter, developed in [12]. We therefore obtain a positivity and monotonicity preserving method as shown in [8]. Finally, a second-order centred scheme is used to discretise the diffusive fluxes, e.g. [14]. The linear systems are solved by MILU(0) preconditioned conjugate gradients with a stopping criterium on the residuals smaller than = 10−8 .
480
X. Cor´e et al. / Mathematics and Computers in Simulation 61 (2003) 477–488
4. Multilevel local mesh refinement method 4.1. The Flux Interface Correction method: general algorithm The Flux Interface Correction method was initially proposed by Angot and coworkers [1,2,9,11]. It is a multilevel mesh refinement technique, and as such it involves two (or more) levels of discretisation: a mesh GH covering the whole computational domain and one or several finer meshes Gh covering subdomains; for the sake of simplicity, we assume in the sequel that there is only one such refined zone. A differential problem is solved on each computational grid, which is basically the “one-grid discrete problem” described in the preceding section; they are referred to in the following as the coarse and fine grid problems, solved, respectively, in the computational domain Ω H and Ω h to obtain the coarse and fine grid solutions. This multigrid technique is relatively cheap since the iterative solvers have to deal only with a relatively small number of degrees of freedom on each discretisation level. The method relies on the following interactions between the coarse and fine grid problems: • the boundary conditions applied for the fine grid problem at the inner boundaries of Ω h (i.e. ∂Ω h \∂Ω H ) are derived by an interpolation technique from the coarse grid solution via the so-called prolongation operator PHh ; • in the coarse grid problem, some corrective terms appear in the right-hand side, evaluated from the fine grid solution by a correction operator RhH . For each time step, solutions on each discretisation level are iteratively computed an arbitrary number of times: typically, two or three iterations are required when quasi-exact solve of the sub-problems is performed, e.g. [2,11]. At the end of the time step, the coarse grid solution is overloaded on Ω h by a quantity evaluated from the fine grid solution via the so-called restriction operator R˜ hH . The two-level FIC iterative solution procedure with the global coarse grid GH and a local fine grid Gh can be decomposed sequentially as shown in Fig. 1, where AD stands for the first step of the incremental projection method (advection– diffusion of velocity), PC is the second step (dynamic pressure Poisson equation and correction of velocity) and ADC is the advection–diffusion of the coarse grid velocity with additional flux correction residuals derived from the fine grid solution. The grid operators PHh , RhH , R˜ hH are detailed in the next section. 4.2. The Flux Interface Correction method: grid operators A correction is performed only in the prediction step of the coarse grid problem, since the projection step consists in a self-defined mathematical operation. The correction operator is built as proposed in the
Fig. 1. Two-level FIC iterative solution procedure for k = 2 iterations.
X. Cor´e et al. / Mathematics and Computers in Simulation 61 (2003) 477–488
481
original papers, and is recalled here only for the sake of completeness. On the other hand, the prolongation and restriction operators have been developed specifically for the present application. As proposed in [2,3], we choose h = H /3 for the mesh steps to simplify the calculations for the grid operators. 4.2.1. Notations Each domain Ω l for l ∈ {H ; h}, with boundary ∂Ω l is partitioned in a finite number of control volumes Kl . The vector nl stands for the unit exterior normal defined on the boundary ∂Kl of each control volume Kl whereas ∂Kl = ∪γl for all faces γ l of Kl . We denote by ωH the subdomain of Ω h obtained by the union of the coarse control volumes that are included in the interior of the refinement area. The velocity k vH (resp. vhk ) stands for the solution of the coarse (resp. fine) grid problem obtained at the FIC iteration n+1 n and vhn (resp. vH and vhn+1 ) denote these solutions at the beginning (resp. end) of the current time k; vH step. The same convention holds for the representations of ρ used for each discrete problem. Finally, for all these fields, we use the notation (ψ H )i,j to refer to the degree of freedom of any discrete field ψ H at ¯ j¯ i, ¯ zj + j¯h) in Ω h for i, ¯ j¯ ∈ {−1; 0 + 1}. any point (xi , zj ) in Ω H and (ψh )i,j for ψ h at the points (xi + ih, 4.2.2. The flux correction operator RhH If the tensor F stands for the convective and diffusive fluxes of the momentum balance equation, i.e. 2 F( v ) = ρ vv − µ(∇ v + ∇ vT ) + µ(∇ · v)Id, 3 we define the flux correction residual vector on each coarse control volume KH ⊂ ωH by: 1 k k rKk+1 ( v ) = F ( v ) · n − F ( v ) · n , k ≥ 0. h H h H H meas(KH ) γ ⊂∂K γh ∂KH h
H
The correction terms δ fHk added to the right-hand side of the momentum equation of the coarse grid problem are then defined by the following recurrence, where |1|KH denotes the characteristic function of the control volume KH : δ fH0 ( v ) = 0,
δ fk+1 ( v ) = δ fHk ( v ) + KH ⊂ωH |1|KH rKk+1 ( v ), k ≥ 0. H H 4.2.3. The prolongation operator PHh Boundary conditions for the fine grid problem must be prescribed on the inner boundaries of Ω h for both steps of the projection method: for the advection–diffusion step, we impose Dirichlet boundary conditions while homogeneous Neumann boundary conditions are prescribed to the pressure time increment for the Poisson problem of the correction step. As an example, in Fig. 2, we define the prolongation PHh operator j¯ j¯ at the fine grid points (xi+(1/2) , zj ), j¯ ∈ {−1; 0; +1}, where zj = zj + j¯h for the first component of the velocity field by: 0,j¯
(uh )i+(1/2),j = PHh (uH ) :=
(ρH )i+(1/2),j 0,j¯
(ρh )i+(1/2),j
(uH )i+(1/2),j .
482
X. Cor´e et al. / Mathematics and Computers in Simulation 61 (2003) 477–488
Fig. 2. Uniform 2D staggered MAC grids with h = H /3 showing the nodes involved in the computation of the grid operators.
4.2.4. The restriction operator R˜ hH In the restriction step, we overload, for example the first component of the velocity at the coarse grid j¯ point (xi +(1/2) , zj ) located at the center of the control volume vertical face γH = ∪ γh for j¯ ∈ {−1; 0; +1} by:
0,j¯ 0,j¯ j¯ j¯∈{−1;0;+1} meas(γh )(ρh )i+(1/2),j (uh )i+(1/2),j H ˜ (uH )i+(1/2),j = Rh (uh ) := . meas(γH )(ρH )i+(1/2),j 4.2.5. About some coherence constraints for the grid operators Let us suppose that we can choose a representation of the density such that the mass quantity is conserved through the grids for each coarse control volume, i.e. ρH = ρh . KH
Kh ⊂KH Kh
Since the grid operators PHh and R˜ hH are indeed specially designed to satisfy the discrete equation: ρH vH · nH = ρh vh · nh , γH
γh ⊂γH γh
the definition of the grid operators ensures that the mass flux through a face of a coarse control volume is conserved by the prolongation and restriction operators and, consequently, that the following relations are verified: 1 n+1 n (ρ − ρ ) + ρhn+1 vhn+1 · nh = 0, h Ω t h Ωh h 1 n+1 n+1 (ρH − ρHn ) + ρHn+1 vH · nH = 0. KH t ∂KH The first relation is a necessary condition to ensure that the Poisson problem of the correction step on the fine grid level has a solution, whenever Neumann conditions are prescribed on the whole boundary of Ω h (situation encountered, for instance, when the refined zone Ω h is strictly interior to Ω H ). The second relation expresses that the mass balance checked by the coarse grid solution still holds after the
X. Cor´e et al. / Mathematics and Computers in Simulation 61 (2003) 477–488
483
restriction. This feature is highly desirable when, for example, the transport of a passive scalar must be treated in addition with the hydrodynamic problem.
5. Numerical results In this section, we present some numerical tests of the discretisation methods presented in this paper. Our objective is two-fold: firstly, to check the time-accuracy of the projection method; secondly, to assess the spatial accuracy gained when using the multilevel local mesh refinement technique. 5.1. Test case For this purpose, we have derived an unsteady test case from a benchmark steady problem proposed by De Lange and De Goey [5]. The computational domain is the square [0; 0.1 m] × [0; 0.1 m] as shown in Fig. 3. The temperature T ( x , t) is given in time and space, and it schematically mimics the conditions of a combustion problem: 2 r (t) , if r 2 (t) ≤ r02 , 300 1 + 2 exp κ 2 r − 1 0 2 T ( x , t) = 1 − r (t) 300 5 − 2 exp κ , otherwise, r02 with r 2 (t) = x 2 (t) + z2 and x(t) = 5x/2 if 0 ≤ t ≤ 1 s, x(t) = 5x/(3t − 1) if 1s ≤ t ≤ 2s, x(t) = x otherwise, κ = 25 and r0 = 0.05 m. The density is then given by ρ( x , t) = Patm /RT ( x , t) for an ideal gas. Fig. 4 shows its time evolution at point P(0.025 m; 0.025 m). The dynamic viscosity µ is prescribed by µ( x , t) = µref (T ( x , t)/1500)γ , with µref = 4.79 × 10−5 Pa s and γ = 0.77. The inflow boundary condition is given by: v( x , t) = (0.25((1 − z2 )/0.0025), 0)T m/s. The fluid is supposed to be initially at rest.
Fig. 3. Geometry and boundary conditions.
484
X. Cor´e et al. / Mathematics and Computers in Simulation 61 (2003) 477–488
Fig. 4. Time evolution of the density at point P(0.025; 0.025).
5.2. Results of one-grid computations Figs. 5 and 6 show, respectively, at steady state t = 3 s, the fields of velocity for the mesh 21 × 21 and dynamic pressure for the mesh 63 × 63. The computations are carried out for CFL ≈ 0.5. The mesh convergence has been studied, e.g. the horizontal velocity profile at steady state is represented for different meshes in space by two cut planes, inside and outside the density layer, respectively, in Figs. 7(a) and
Fig. 5. Steady-state velocity.
X. Cor´e et al. / Mathematics and Computers in Simulation 61 (2003) 477–488
Fig. 6. Steady-state dynamic pressure.
Fig. 7. Steady-state horizontal velocity profiles for meshes: (- - -) 21 × 21; (· · · ) 63 × 63; (—) 189 × 189 (reference).
485
486
X. Cor´e et al. / Mathematics and Computers in Simulation 61 (2003) 477–488
Fig. 8. Horizontal velocity time evolution at point P(0.025; 0.025) for the mesh 63 × 63 and time steps: (· · · ) t = 10−3 s; (- - -) t = 5 × 10−4 s.
(b). As expected, the mesh 21 × 21 gives rough results in both figures with relative errors of about 10% due to a poor representation of the density gradients, whereas the solutions for the mesh 63 × 63 and the reference mesh 189 × 189 are superposed. As shown in Fig. 8 for the time evolution of the horizontal velocity at point P(0.025 m; 0.025 m), the time step convergence is achieved for the mesh 63 × 63 at least with a time step t = 10−3 s. 5.3. Results of two-grid computations As an example, Fig. 9 shows the coarse grid GH and the fine grid Gh with h = H /3 covering the zoom area which includes the density layer. In Fig. 10, the steady-state horizontal velocity profiles at x = 0.025 m and x = 0.075 m, i.e. inside and outside the zoom area, respectively, are represented for the different meshes. At both locations, the coarse grid solution of the refined mesh 21 × 21 is better than
Fig. 9. Zoom area: coarse grid GH and fine grid Gh where h = H /3.
X. Cor´e et al. / Mathematics and Computers in Simulation 61 (2003) 477–488
487
Fig. 10. Steady-state horizontal velocity profiles for meshes: (- - -) coarse grid solution of refined mesh 21 × 21; (—) 189 × 189; (· · · ) fine grid solution of refined mesh 21 × 21 (only for x = 0.025 m in the zoom area).
the solution with no zoom on the mesh 21 × 21. Inside the refinement zone, at x = 0.025 m, the fine grid solution of the refined mesh 21 × 21 and the reference solution are superimposed. Moreover, outside the refinement zone in Fig. 10, the coarse grid solution of the refined mesh 21 × 21 is superposed on the reference solution. Hence, the multilevel zoom method greatly improves the solution precision not only inside but also outside the refinement area.
6. Conclusion This study proves the efficiency of the designed solver to calculate variable density flows with the multilevel local mesh refinement method FIC associated with the semi-implicit incremental projection. Besides, it shows how the use of such multilevel zoom solvers can be a cheap, efficient and precise enough alternative to more standard computations with non-uniform grids or unstructured meshes.
Acknowledgements This work was supported by the Institut de Protection et de Sˆureté Nucléaire (IPSN), in the framework of its modelling activities on fire consequences in nuclear power plants. Practical implementations were performed using the object oriented framework and software components library PELICANS developed at IPSN.
488
X. Cor´e et al. / Mathematics and Computers in Simulation 61 (2003) 477–488
References [1] P. Angot, J.P. Caltagirone, K. Khadra, Une méthode adaptative de raffinement local: la Correction du Flux à l’Interface, Comptes-Rendus de l’Académie des Sciences de Paris 315 (1) (1992) 739–745. [2] P. Angot, M. Laugier, La méthode FIC de raccordement conservatif de sous-domaines emboˆıtés pour un modèle de circulation océanique, C. R. Acad. Sci. Paris 319 (2) (1994) 993–1000. [3] J.P. Caltagirone, K. Khadra, P. Angot, Sur une méthode de raffinement local multigrille pour la résolution des équations de Navier–Stokes, C. R. Acad. Sci. Paris 320 (2-b) (1995) 295–302. [4] A.J. Chorin, Numerical solution of the Navier–Stokes equations, Math. Comput. 22 (1968) 745–762. [5] H.C. De Lange, P.H. De Goey, Numerical flow modelling in a locally refined grid, Int. J. Numer. Methods Eng. 37 (1994) 497–515. [6] H.N. Najm, P.S. Wyckoff, O.M. Knio, A semi-implicit numerical scheme for reacting flow. I. Stiff chemistry, J. Comput. Phys. 143 (1998) 381–402. [7] T. Hayase, J.A.C. Humphrey, R. Greif, A consistently formulated QUICK scheme for fast and stable convergence using finite-volume iterative calculation procedures, J. Comput. Phys. 98 (1992) 108–118. [8] H.T. Huynh, Accurate monotonic cubic interpolation, SIAM J. Numer. Anal. 30 (1) (1993) 57–100. [9] K. Khadra, P. Angot, J.P. Caltagirone, P. Morel, Concept de zoom adaptatif en architecture multigrille locale; Etude comparative des méthodes LDC, FAC et FIC, Modélisation Mathématique et Analyse Numérique 30 (1) (1996) 39–82. [10] R. Klein, Numerics in combustion, in: L. Vervisch, D. Veynante (Eds.), Introduction to Turbulent Combustion, Lecture Series in Computational Fluid Dynamics, Von Karman Institute, Belgium, March 1999. [11] M. Laugier, P. Angot, L. Mortier, Nested-grid methods for an ocean model: a comparative study, Int. J. Numer. Methods Fluids 23 (11) (1996) 1163–1195. [12] B.P. Leonard, S. Mokthari, Beyond first-order upwinding: the ULTRA-SHARP alternative for non-oscillatory steady-state simulation of convection, Int. J. Numer. Methods Eng. 30 (1990) 729–766. [13] A. Majda, J. Sethian, The derivation and numerical solution of the equations for zero Mach number combustion, Combust. Sci. Tech. 42 (1985) 185–205. [14] S.V. Patankar, Numerical Heat Transfer and Fluid Flow, Computational Methods in Mechanics and Thermal Sciences, McGraw-Hill, New York, 1980.