Proceedings of Combustion Institute - Canadian Section Spring Technical Meeting University of Toronto, Ontario May 12-14,2008

Parallel Implicit Adaptive Mesh Refinement Algorithm for Unsteady Laminar Flames Scott A. Northrup∗ , Clinton P. T. Groth† Institute for Aerospace Studies, University of Toronto, 4925 Dufferin Street, Toronto, ON, M3H 5T6, Canada 1

Introduction

Numerical methods have become an essential tool for investigating combustion processes. However despite the significant advances in solution algorithms and computer hardware, obtaining accurate numerical solutions can still place severe demands on available computational resources. Many approaches have been taken to reduce the computational costs of simulating combusting flows. One successful approach is to make use of solution-directed mesh adaptation, such as the adaptive mesh refinement refinement (AMR) algorithms developed for aerospace applications [1–6]. Computational grids that automatically adapt to the solution of the governing equations are very effective in treating problems with disparate length scales, providing the required spatial resolution while minimizing memory and storage requirements. Recent progress in the development and application of AMR algorithms for low-Mach-number reacting flows and premixed turbulent combustion is described by Day and Bell [7]. Another approach for coping with the computational cost of reacting flow prediction is to apply a domain decomposition procedure and solve the problem in a parallel fashion using multiple processors. Large massively parallel distributed-memory computers can provide many fold increases in processing power and memory resources beyond those of conventional single-processor computers and would therefore provide an obvious avenue for greatly reducing the time required to obtain numerical solutions of combusting flows. Douglas et al. [8] describe a parallel algorithm for numerical combustion modelling. Northrup and Groth [9] and Gao and Groth [10] combined these two numerical approaches, producing a parallel AMR method for both non-premixed laminar and turbulent combusting flows. A preconditioned nonlinear multigrid algorithm with multi-stage semi-implicit time marching scheme as a smoother was used to integrate the governing partial differential equations. Although accurate solutions were obtained, the approach was not optimal as in many cases a large number of multigrid cycles and solution residual evaluations were required to obtain steady-state solutions. In order to reduce the time to achieve a solution and deal with the numerical stiffness of the reactive flow problems, a parallel implicit AMR scheme has been developed that combines the block-based AMR approach with Newton’s method [11]. In this approach, Newton’s method is used to solve the system of nonlinear equations arising from a upwind finite-volume spatial discretization procedure and a preconditioned generalized minimal residual (GMRES) method is used to solve the resulting system of linear equations at each step of the Newton algorithm. An additive Schwarz preconditioner is used in combination with block-fill incomplete lower-upper (BFILU) preconditioning to improve performance of the linear iterative solver. The Schwarz preconditioning and block-based AMR readily allow efficient and scalable parallel implementations of the implicit approach on distributed-memory multi-processor architectures. Furthermore, this parallel implicit method is particularly well suited for predicting, in a reliable and efficient fashion, the physically complex flows containing widely disparate spatial and temporal scales encountered in many combustion processes. This paper then considers the extension of the parallel implicit AMR scheme of Northrup and Groth [11,12], originally developed exclusively for steady combustion problems, to the prediction of unsteady laminar diffusion flames. A dualtime-stepping-like procedure [13] is used to maintain time accuracy. Details are provided concerning the mathematical models used to describe the laminar flow combustion processes and a description is given of the parallel implicit AMR finite-volume scheme, including details of the inexact Newton and GMRES iterative methods, dual-time-steppinglike approach, AMR strategy, and parallel implementation. Numerical results are described for methane-air laminar ∗ †

[email protected] [email protected]

diffusion flames. The results are used to demonstrate the validity of the parallel implicit AMR approach and the efficiency of the method for resolving fine-scale features of laminar flame structures. 2

Mathematical Modelling

The Navier-Stokes equations for a thermally perfect reactive mixture, given by ∂ ~) = 0, (ρ) + ∇ · (ρU ∂t ∂ ~ ~U ~ + p~~I) = ∇ · ~~τ − ρ~g , (ρU ) + ∇ · (ρU ∂t ∂ ~ (E + p ) = ∇ · (U ~ · ~~τ − ~q) − ρ~g · U ~, [ρE] + ∇ · ρU ∂t ρ

(1) (2) (3)

∂ ~ ) = ∇ · (ρDs ∇cs ) + ρω˙s , (ρcs ) + ∇ · (ρcs U (4) ∂t are used herein to model the combustion of gaseous fuels and oxidizers. Equations (1)–(3) reflect the conservation of ~ is the mixture velocity, E is mass, momentum, and energy for the reactive mixture, ρ is the mixture mass density, U 1 ~ 2 the total specific energy of the mixture given by E = e + 2 |U | , e is the specific internal energy, p is the mixture pressure, ~~τ is the fluid stress tensor for the mixture, ~q is the heat flux vector, and ~g is the acceleration due to gravitational forces. Equation (4) is the species concentration equation for species s, where cs is the species mass fraction, Ds is the diffusion coefficient, and ω˙s is the time rate of change of the species concentration PNdue to finite-rate chemistry. It follows from the caloric equation of state for a thermally perfect mixture that e = s=1 cs hs − ρp , where hs is the PN species enthalpy, N is the number of species, and the ideal gas law for the mixture is given by p = s=1 ρcs Rs T , where Rs is the species gas constant and T is the mixture temperature. Note that the use of the compressible form of the Navier-Stokes equations readily allows for the often large density variations associated with combusting flows and the prediction of thermoacoustic phenomena. Expressions for individual species and mixture thermodynamic and transport properties are required to complete the system of partial differential equations for the reactive mixture. The empirical expressions complied by Gordon and McBride [14, 15] are used to specify the hs and the species specific heat, cps , entropy, ∆ss , viscosity, µs , and thermal conductivity, κs , as functions of temperature. The Gordon-McBride dataset contains curve fits for over 2000 substances, including 50 reference elements. Perfect mixture rules are used to determine the thermodynamic properties of the reactive mixture and Wilke’s [16] and Mason and Saxena’s [17] mixture rules are used to evaluate the mixture viscosity, µ, and thermal conductivity, κ, respectively. The species diffusion coefficients, Ds are calculated by additionally specifying a Schmidt number, Scs = µ/ρDs for each species. For the cases considered herein, a simplified five-species one-step reduced chemical reaction mechanisms of Westbrook and Dryer [18] is used to model the chemical kinetics of methane-air combustion. Further details and reaction rates for this reduced mechanism are given by Westbrook and Dryer [18]. 3 3.1

Parallel Implicit AMR Scheme Finite-Volume Discretization

The governing equations for the reactive mixture are solved by applying a finite-volume method in which the mixture conservation equations are integrated over quadrilateral cells of a body-fitted multi-block quadrilateral mesh. The finite-volume formulation applied to cell (i, j) can be expressed by 1 X ~ dUi,j =− F · ~n ∆` + Si,j = Ri,j (U) , (5) Γ dt Ai,j i,j,k k

~ is the flux dyad containing contributions from the inviscid and viscous terms, S is the source term associated where F with the axisymmetric geometry, finite rate chemistry, and gravitational forces, Ai,j is the area of cell (i, j), and ∆`

and ~n are the length of the cell face and unit vector normal to the cell face or edge, respectively and Ri,j (U) is the residual vector. The local preconditioning technique proposed by Weiss and Smith [19] is used to alleviate numerical difficulties for low-Mach-number, nearly incompressible flows that generally characterize laminar combusting flows. The preconditioning matrix, Γ, does not affect the steady-state solution but helps control the numerical stiffness and dissipation, making the solution of the governing more tractable. The inviscid component of the numerical fluxes at the faces of each cell are determined using the least-squares piece-wise limited linear solution reconstruction procedure of Barth [20] and Riemann solver based flux functions. The limiter of Venkatakrishnan [21] is used. An extension of the approximate linearized Riemann solver of Roe [22] is used to account for mixture composition. The viscous component of the cell face fluxes are evaluated by employing a diamond-path reconstruction procedure as described by Coirier and Powell [23]. The approach allows for solution-directed block-based AMR and an efficient and highly scalable parallel implementation has been achieved via domain decomposition. Refer to the recent paper by Northrup and Groth [11] for details of the parallel AMR scheme. 3.2

Newton’s Method

The semi-discrete form of the governing equations given in Eq. (5) form a coupled set of non-linear ordinary differential equations. For unsteady flows, time-dependent solutions are obtained by employing a dual-time-stepping-like procedure [13]. In the proposed approach a modified residual is defined by R∗ (U) = dU dt + R(U). Applying an implicit second-order backward discretization of the physical time derivative yields R∗ (U(n+1) ) =

3U(n+1) − 4U(n) + U(n−1) + R(U(n+1) ) = 0, 24t

and using Newton’s method to the solution of Eq. 6 leads to the following linear system of equations: 3 ∂R I+ 4U(n+1,k) = J∆U(n+1,k) = −R∗ (U(n+1,k) ), 24t ∂U

(6)

(7)

for the solution change ∆U(n+1) = U(n+1) − U(n) at time level n, with Newton iteration level k. Using the previous time step as the initial estimate, U(n+1,k=0) = U(n) , successively improved estimates for the solution, U(n+1,k) , are obtained by solving Eq. 7 at each step, k, of the Newton method, where J is the modified residual Jacobian. The iterative procedure is repeated until an appropriate norm of the solution residual is sufficiently small, i.e., ||R∗ (U(n+1,k) )||2 < ||R∗ (U(n) )||2 where is some small parameter (typically, ≈ 10−2 –10−3 ). Each step of Newton’s method requires the solution of a system of linear equations of the form Jx = b. This system is large, sparse, and non-symmetric and a preconditioned GMRES method [24, 25] is used for its solution. In particular, a restarted version of the GMRES algorithm, GMRES(m), is used, where m is the number of steps after which the GMRES algorithm is restarted. Application of this iterative technique leads to an overall solution algorithm with iterations within iterations: the “inner loop” iterations involving the solution of the linear system and the “outer loop” iterations associated with the solution of the nonlinear problem. An inexact Newton method is adopted here in which the inner iterations are not fully converged at each Newton step. The inner iterations are carried out only until ||R∗ + J∆U||2 ≤ ζ||R∗ ||2 , where ζ is typically in the range 0.01–0.5. As discussed by Dembo et al. [26], an exact solution of the linear system is not necessary for rapid convergence of Newton’s method. Preconditioning is required for the linear solver to be effective. Right preconditioning of the form (JM−1 )(Mx) = b is used here where M is the preconditioning matrix. An additive Schwarz global preconditioner with variable overlap [25,27,28] is used in conjunction with local BFILU preconditioners for each sub-domain. The local preconditioner is based on a block-fill ILU(f ) or BFILU(f ) factorization of an approximate Jacobian for each subdomain. Here, f is the level of fill. This combination of preconditioning fits well with the block-based AMR described by Northrup and Groth [11] and is compatible with domain decomposition methods, readily enabling parallel implementation of the overall Newton method. Rather efficient parallel implementations of implicit algorithms via Schwarz preconditioning have been developed by Keyes and co-researchers and successfully applied to the prediction of transonic full potential, low-Mach-number compressible combusting, and three-dimensional inviscid flows [27, 29, 30]. As the GMRES algorithm does not explicitly require the evaluation of the global Jacobian matrix, J, a so-called “matrix-free” or “Jacobian-free” approach can be adopted and is used here. Numerical differentiation based on Fréchet

derivatives is used to approximate the matrix-vector product JM−1 x as follows: JM−1 x ≈

R(U + εM−1 x) − R(U) 3M−1 x + , ε 24t

(8)

where R(U+εM−1 x) is the residual vector evaluated at some perturbed solution state and ε is a small scalar quantity. Although the performance of the Jacobian-free method is sensitive to the choice of ε, Neilsen et al. [31] have found 1/2 that ε = ε◦ /||x||2 seems to work well, with ε◦ ≈ 10−8 –10−7 . 4

Numerical Results

The proposed parallel implicit AMR algorithm is now applied to the solution of a forced, time-dependent, axisymmetric, methane-air coflow laminar diffusion flame. In particular, the solution of the flame studied by Day and Bell [32], and Dworkin et al. [33] is considered. The flame boundary and initial conditions are the same as those used in the previous studies. The computational domain is rectangular in shape with dimensions of 10 cm by 5 cm. The axis of symmetry is aligned with the left boundary of the domain and the right far-field boundary is taken to be a freeslip boundary along which inviscid reflection boundary data is specified. The top or outlet of the flow domain is open to a stagnant reservoir at atmospheric pressure and temperature and Neumann-type boundary conditions are applied to all properties except pressure which is held constant. The bottom or inlet is subdivided into four regions. The innermost region (r ≤ 2 mm) is the fuel inlet or jet, which injects a nitrogen diluted methane fuel mixture (cCH4 = 0.5149, cN2 = 0.4851, cO2 = 0, cCO2 = 0, cCO = 0, and cH2 O = 0) at 298 K with a parabolic axial velocity profile. The time-variation in the flame is produced by imposing a sinusoidal axial velocity fluctuation across the fuel jet, vz = 70.0(1 − r2 /R)(1 + α sin ωt) cm/s where α is the velocity amplitude and ω is the frequency of oscillation. Matching the experiments of Dworkin et al. [33], the velocity amplitude was set at 50% at a frequency of 20Hz. The next region (2 mm < r ≤ 2.38 mm) is a small gap associated with the annular wall separating the fuel and oxidizer. The third region (2.38 mm < r ≤ 2.50 cm) is the co-flowing oxidizer, in this case air at 298K (cO2 = 0.232, cN2 = 0.768, cCH4 = 0, cCO2 = 0, cCO = 0, and cH2 O = 0), with a uniform velocity profile of 0.35 m/s. The final outer region of the lower boundary (2.5 cm < r ≤ 5 cm) is again a far-field boundary along which free-slip boundary conditions are applied. The solution domain is initialized with a uniform solution state corresponding to quiescent air at 298K, except for a thin region across the fuel and oxidizer inlets, which is taken to be air at 1500 K so as to ignite the flame. Additional details concerning the setup for this diffusion flame can be found in the papers by Day and Bell [32], and Dworkin et al. [33] The flame was calculated by first letting the steady-state flame develop with the fuel flow velocity held constant at 70 cm/s, then four-full periods of the flame oscillation were calculated to avoid any non-periodic oscillations created during start-up. Newton’s method was used with a GMRES tolerance of 0.1, ILU(4) preconditioner and at each time step the Newton iterations were converged two orders of magnitude, to a maximum of 10 Newton steps, where 5-7 is typical. The approximate Jacobian preconditioner was only updated for the first Newton step of each time-step, unless the number of GMRES iterations required increased. A fixed time-step of 0.0001s was used during the unsteady calculation. Mesh refinement was carried out every 0.0025s based on the gradient of temperature, and 3 levels of mesh refinement (4 mesh levels) were used as shown in Figure 1. Figure 2 shows the resulting isotherms at five times through the flames periodic fluctuation of 0.05s. The flames shape and structure match well with the experimental results of Dworkin et al. [33]. The most significant difference is in the over prediction of temperature, however this is related to the use of the simplified non-reversible one-step reaction mechanism used to model the methane-air chemistry. 5

Conclusions

A parallel implicit AMR scheme has been developed for solving unsteady laminar combusting flows. The combination of finite-volume discretization procedure, parallel block-based AMR strategy, low-Mach-number preconditioning, dual-time-stepping-like apporach, and Newton method solution procedure has resulted in a powerful computational tool for predicting a wide range of unsteady laminar reactive flows, from compressible to nearly incompressible lowMach-number regimes. The validity and performance of the method has been demonstrated for non-premixed flames. Future work will involve the extension of the the approach to unsteady three-dimensional combusting flows.

a) 0.01s

b) 0.03s

c) 0.05s

0.1

0.08

0.08

0.06

z (m)

z (m)

0.06

0.04

0.04

0.02

0.02

0 -0.02

-0.01

0

0.01

0.02

0 -0.02

-0.01

0

r (m)

0.01

0.02

-0.02

-0.01

0

r (m)

0.01

0.02

r (m)

Fig. 1: Time-varying methane-air co-flow laminar diffusion flame isotherms at three intervals calculated with a) 456, b) 444, and c) 447 (4x8 cell) blocks each with 3 levels of refinement. a) 0.01s

b) 0.02s

c) 0.03s

0.08

0

r (m)

0.005

0.01

z (m)

z (m)

0.06

z (m)

0.06

z (m)

0.06

0.04

0.04

0.04

0.04

0.02

0.02

0.02

0.02

0 -0.01

-0.005

0

r (m)

0.005

0.01

0 -0.01

-0.005

0

r (m)

0.005

0.01

T (K)

0.08

0.08

0.06

0 -0.005

e) 0.05s

d) 0.04s

0.08

2200 2162 2124 2087 2049 2011 1973 1936 1898 1860 1822 1785 1747 1709 1671 1634 1596 1558 1520 1483 1445 1407 1369 1332 1294 1256 1218 1181 1143 1105 1067 1030 992 954 916 879 841 803 765 728 690 652 614 577 539 501 463 426 388 350

0 -0.01

-0.005

0

r (m)

0.005

0.01

-0.01

-0.005

0

0.005

0.01

r (m)

Fig. 2: Time-varying methane-air co-flow laminar diffusion flame isotherms at five 0.01ms, (a-e) intervals.

References [1] [2] [3] [4] [5] [6] [7] [8]

[9] [10] [11]

[12]

[13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33]

M. J. Berger. J. Comput. Phys., 53:484–512, 1984. M. J. Berger and P. Colella. J. Comput. Phys., 82:67–84, 1989. D. De Zeeuw and K. G. Powell. J. Comput. Phys., 104:56–68, 1993. M. J. Aftomis, M. J. Berger, and J. E. Melton. AIAA J., 36(6):952–960, 1998. C. P. T. Groth, D. L. De Zeeuw, K. G. Powell, T. I. Gombosi, and Q. F. Stout. A parallel solution-adaptive scheme for ideal magnetohydrodynamics. Paper 99-3273, AIAA, June 1999. J. S. Sachdev, C. P. T. Groth, and J. J. Gottlieb. Int. J. Comput. Fluid Dyn., 19(2):157–175, 2005. J.B. Bell, M.S. Day, A.S. Almgren, M. J. Lijewski, and C. A. Rendleman. Int. J. Numer. Meth. Fluids, 40:209–216, 2002. C. C. Douglas, A. Ern, and M. D. Smooke. Numerical simulation of flames using multigrid methods. In J. Wang, M. B. Allen, B. M. Chen, and T. Mathew, editors, Iterative Methods in Scientific Computation, volume 4 of IMACS Series in Computational and Applied Mathematics, pages 149–154. New Brunswick, 1998. S. A. Northrup and C. P. T. Groth. Solution of laminar diffusion flames using a parallel adaptive mesh refinement algorithm. Paper 2005–0547, AIAA, January 2005. X. Gao and C. P. T. Groth. Parallel adaptive mesh refinement scheme for turbulent non-premixed combusting flow prediction. Paper 2006–1448, AIAA, January 2006. S. A. Northrup and C. P. T. Groth. Solution of laminar combusting flows using a parallel implicit adaptive mesh refinement algorithm. In Proceedings of the 4th International Conference on Computational Fluid Dynamics, Ghent, Belgium, July 10–14, 2006, 2006. S. A. Northrup and C. P. T. Groth. Parallel implicit adaptive mesh refinement algorithm for solution of laminar combusing flows. In Proceedings of the 14th Annual Conference of the CFD Society of Canada, Kingston, Canada, July 16–18, 2006, pages 221–228. CFD Society of Canada, 2006. A. Jameson. Time dependent calculations using multigrid, with applications to unsteady flows past airfoils and wings. Paper 1991-1596, AIAA, June 1991. S. Gordon and B. J. McBride. Computer program for calculation of complex chemical equilibrium compositions and applications I. analysis. Reference Publication 1311, NASA, 1994. B. J. McBride and S. Gordon. Computer program for calculation of complex chemical equilibrium compositions and applications ii. users manual and program description. Reference Publication 1311, NASA, 1996. C. R. Wilke. J. Chem. Phys., 18:517–519, 1950. G. Dixon-Lewis. Computer modeling of combustion reactions in flowing systems with transport. In W. C. Gardiner, editor, Combustion Chemistry, pages 21–126. Springer-Verlag, New York, 1984. C. K. Westbrook and F. L. Dryer. Combust. Sci. Tech., 27:31–43, 1981. J. M. Weiss and W. A. Smith. AIAA J., 33(11):2050–2057, 1995. T. J. Barth. Recent developments in high order k-exact reconstruction on unstructured meshes. Paper 93-0668, AIAA, January 1993. V. Venkatakrishnan. On the accuracy of limiters and convergence to steady state solutions. Paper 93-0880, AIAA, January 1993. P. L. Roe. J. Comput. Phys., 43:357–372, 1981. W. J. Coirier and K. G. Powell. AIAA J., 34(5):938–945, May 1996. Y. Saad and M. H. Schultz. SIAM J. Sci. Stat. Comput., 7(3):856–869, 1986. Y. Saad. Iterative Methods for Sparse Linear Systems. PWS Publishing Company, Boston, 1996. R. S. Dembo, S. C. Eisenstat, and T. Steihaug. SIAM J. Numer. Anal., 19(2):400–408, 1982. W. D. Gropp, D. K. Kaushik, D. E. Keyes, and B. F. Smith. Parallel Computing, 27:337–362, 2001. D. A. Knoll and D. E. Keyes. J. Comput. Phys., 193:257–397, 2004. D. A. Knoll, P. R. McHugh, and D. E. Keyes. AIAA J., 34(5):961–967, 1996. X.-C. Cai, W. D. Gropp, D. E. Keyes, R. G. Melvin, and D. P. Young. SIAM J. Sci. Comput., 19(1):246–265, 1998. E. J. Nielsen, W. K. Anderson, R. W. Walters, and D. E. Keyes. Application of Newton-Krylov methodology to a threedimensional unstructured Euler code. Paper 95-1733-CP, AIAA, June 1995. M. S. Day and J. B. Bell. Combust. Theory Modelling, 4(4):535–556, 2000. S. B. Dworkin, B. C. Connelly, A. M. Schaffer, B. A. V. Bennett, M. B. Long, M. D. Smooke, M. P. Puccio, B. McAndrews, and J. H. Miller. Computational and experimental study of a forced, time-dependent, methane-air coflow diffusion flame. In Thirty First Symposium (International) on Combustion, pages 971–978, Pittsburgh, 2007. The Combustion Institute.

Parallel Implicit Adaptive Mesh Refinement Algorithm for Unsteady Laminar Flames Scott A. Northrup∗ , Clinton P. T. Groth† Institute for Aerospace Studies, University of Toronto, 4925 Dufferin Street, Toronto, ON, M3H 5T6, Canada 1

Introduction

Numerical methods have become an essential tool for investigating combustion processes. However despite the significant advances in solution algorithms and computer hardware, obtaining accurate numerical solutions can still place severe demands on available computational resources. Many approaches have been taken to reduce the computational costs of simulating combusting flows. One successful approach is to make use of solution-directed mesh adaptation, such as the adaptive mesh refinement refinement (AMR) algorithms developed for aerospace applications [1–6]. Computational grids that automatically adapt to the solution of the governing equations are very effective in treating problems with disparate length scales, providing the required spatial resolution while minimizing memory and storage requirements. Recent progress in the development and application of AMR algorithms for low-Mach-number reacting flows and premixed turbulent combustion is described by Day and Bell [7]. Another approach for coping with the computational cost of reacting flow prediction is to apply a domain decomposition procedure and solve the problem in a parallel fashion using multiple processors. Large massively parallel distributed-memory computers can provide many fold increases in processing power and memory resources beyond those of conventional single-processor computers and would therefore provide an obvious avenue for greatly reducing the time required to obtain numerical solutions of combusting flows. Douglas et al. [8] describe a parallel algorithm for numerical combustion modelling. Northrup and Groth [9] and Gao and Groth [10] combined these two numerical approaches, producing a parallel AMR method for both non-premixed laminar and turbulent combusting flows. A preconditioned nonlinear multigrid algorithm with multi-stage semi-implicit time marching scheme as a smoother was used to integrate the governing partial differential equations. Although accurate solutions were obtained, the approach was not optimal as in many cases a large number of multigrid cycles and solution residual evaluations were required to obtain steady-state solutions. In order to reduce the time to achieve a solution and deal with the numerical stiffness of the reactive flow problems, a parallel implicit AMR scheme has been developed that combines the block-based AMR approach with Newton’s method [11]. In this approach, Newton’s method is used to solve the system of nonlinear equations arising from a upwind finite-volume spatial discretization procedure and a preconditioned generalized minimal residual (GMRES) method is used to solve the resulting system of linear equations at each step of the Newton algorithm. An additive Schwarz preconditioner is used in combination with block-fill incomplete lower-upper (BFILU) preconditioning to improve performance of the linear iterative solver. The Schwarz preconditioning and block-based AMR readily allow efficient and scalable parallel implementations of the implicit approach on distributed-memory multi-processor architectures. Furthermore, this parallel implicit method is particularly well suited for predicting, in a reliable and efficient fashion, the physically complex flows containing widely disparate spatial and temporal scales encountered in many combustion processes. This paper then considers the extension of the parallel implicit AMR scheme of Northrup and Groth [11,12], originally developed exclusively for steady combustion problems, to the prediction of unsteady laminar diffusion flames. A dualtime-stepping-like procedure [13] is used to maintain time accuracy. Details are provided concerning the mathematical models used to describe the laminar flow combustion processes and a description is given of the parallel implicit AMR finite-volume scheme, including details of the inexact Newton and GMRES iterative methods, dual-time-steppinglike approach, AMR strategy, and parallel implementation. Numerical results are described for methane-air laminar ∗ †

[email protected] [email protected]

diffusion flames. The results are used to demonstrate the validity of the parallel implicit AMR approach and the efficiency of the method for resolving fine-scale features of laminar flame structures. 2

Mathematical Modelling

The Navier-Stokes equations for a thermally perfect reactive mixture, given by ∂ ~) = 0, (ρ) + ∇ · (ρU ∂t ∂ ~ ~U ~ + p~~I) = ∇ · ~~τ − ρ~g , (ρU ) + ∇ · (ρU ∂t ∂ ~ (E + p ) = ∇ · (U ~ · ~~τ − ~q) − ρ~g · U ~, [ρE] + ∇ · ρU ∂t ρ

(1) (2) (3)

∂ ~ ) = ∇ · (ρDs ∇cs ) + ρω˙s , (ρcs ) + ∇ · (ρcs U (4) ∂t are used herein to model the combustion of gaseous fuels and oxidizers. Equations (1)–(3) reflect the conservation of ~ is the mixture velocity, E is mass, momentum, and energy for the reactive mixture, ρ is the mixture mass density, U 1 ~ 2 the total specific energy of the mixture given by E = e + 2 |U | , e is the specific internal energy, p is the mixture pressure, ~~τ is the fluid stress tensor for the mixture, ~q is the heat flux vector, and ~g is the acceleration due to gravitational forces. Equation (4) is the species concentration equation for species s, where cs is the species mass fraction, Ds is the diffusion coefficient, and ω˙s is the time rate of change of the species concentration PNdue to finite-rate chemistry. It follows from the caloric equation of state for a thermally perfect mixture that e = s=1 cs hs − ρp , where hs is the PN species enthalpy, N is the number of species, and the ideal gas law for the mixture is given by p = s=1 ρcs Rs T , where Rs is the species gas constant and T is the mixture temperature. Note that the use of the compressible form of the Navier-Stokes equations readily allows for the often large density variations associated with combusting flows and the prediction of thermoacoustic phenomena. Expressions for individual species and mixture thermodynamic and transport properties are required to complete the system of partial differential equations for the reactive mixture. The empirical expressions complied by Gordon and McBride [14, 15] are used to specify the hs and the species specific heat, cps , entropy, ∆ss , viscosity, µs , and thermal conductivity, κs , as functions of temperature. The Gordon-McBride dataset contains curve fits for over 2000 substances, including 50 reference elements. Perfect mixture rules are used to determine the thermodynamic properties of the reactive mixture and Wilke’s [16] and Mason and Saxena’s [17] mixture rules are used to evaluate the mixture viscosity, µ, and thermal conductivity, κ, respectively. The species diffusion coefficients, Ds are calculated by additionally specifying a Schmidt number, Scs = µ/ρDs for each species. For the cases considered herein, a simplified five-species one-step reduced chemical reaction mechanisms of Westbrook and Dryer [18] is used to model the chemical kinetics of methane-air combustion. Further details and reaction rates for this reduced mechanism are given by Westbrook and Dryer [18]. 3 3.1

Parallel Implicit AMR Scheme Finite-Volume Discretization

The governing equations for the reactive mixture are solved by applying a finite-volume method in which the mixture conservation equations are integrated over quadrilateral cells of a body-fitted multi-block quadrilateral mesh. The finite-volume formulation applied to cell (i, j) can be expressed by 1 X ~ dUi,j =− F · ~n ∆` + Si,j = Ri,j (U) , (5) Γ dt Ai,j i,j,k k

~ is the flux dyad containing contributions from the inviscid and viscous terms, S is the source term associated where F with the axisymmetric geometry, finite rate chemistry, and gravitational forces, Ai,j is the area of cell (i, j), and ∆`

and ~n are the length of the cell face and unit vector normal to the cell face or edge, respectively and Ri,j (U) is the residual vector. The local preconditioning technique proposed by Weiss and Smith [19] is used to alleviate numerical difficulties for low-Mach-number, nearly incompressible flows that generally characterize laminar combusting flows. The preconditioning matrix, Γ, does not affect the steady-state solution but helps control the numerical stiffness and dissipation, making the solution of the governing more tractable. The inviscid component of the numerical fluxes at the faces of each cell are determined using the least-squares piece-wise limited linear solution reconstruction procedure of Barth [20] and Riemann solver based flux functions. The limiter of Venkatakrishnan [21] is used. An extension of the approximate linearized Riemann solver of Roe [22] is used to account for mixture composition. The viscous component of the cell face fluxes are evaluated by employing a diamond-path reconstruction procedure as described by Coirier and Powell [23]. The approach allows for solution-directed block-based AMR and an efficient and highly scalable parallel implementation has been achieved via domain decomposition. Refer to the recent paper by Northrup and Groth [11] for details of the parallel AMR scheme. 3.2

Newton’s Method

The semi-discrete form of the governing equations given in Eq. (5) form a coupled set of non-linear ordinary differential equations. For unsteady flows, time-dependent solutions are obtained by employing a dual-time-stepping-like procedure [13]. In the proposed approach a modified residual is defined by R∗ (U) = dU dt + R(U). Applying an implicit second-order backward discretization of the physical time derivative yields R∗ (U(n+1) ) =

3U(n+1) − 4U(n) + U(n−1) + R(U(n+1) ) = 0, 24t

and using Newton’s method to the solution of Eq. 6 leads to the following linear system of equations: 3 ∂R I+ 4U(n+1,k) = J∆U(n+1,k) = −R∗ (U(n+1,k) ), 24t ∂U

(6)

(7)

for the solution change ∆U(n+1) = U(n+1) − U(n) at time level n, with Newton iteration level k. Using the previous time step as the initial estimate, U(n+1,k=0) = U(n) , successively improved estimates for the solution, U(n+1,k) , are obtained by solving Eq. 7 at each step, k, of the Newton method, where J is the modified residual Jacobian. The iterative procedure is repeated until an appropriate norm of the solution residual is sufficiently small, i.e., ||R∗ (U(n+1,k) )||2 < ||R∗ (U(n) )||2 where is some small parameter (typically, ≈ 10−2 –10−3 ). Each step of Newton’s method requires the solution of a system of linear equations of the form Jx = b. This system is large, sparse, and non-symmetric and a preconditioned GMRES method [24, 25] is used for its solution. In particular, a restarted version of the GMRES algorithm, GMRES(m), is used, where m is the number of steps after which the GMRES algorithm is restarted. Application of this iterative technique leads to an overall solution algorithm with iterations within iterations: the “inner loop” iterations involving the solution of the linear system and the “outer loop” iterations associated with the solution of the nonlinear problem. An inexact Newton method is adopted here in which the inner iterations are not fully converged at each Newton step. The inner iterations are carried out only until ||R∗ + J∆U||2 ≤ ζ||R∗ ||2 , where ζ is typically in the range 0.01–0.5. As discussed by Dembo et al. [26], an exact solution of the linear system is not necessary for rapid convergence of Newton’s method. Preconditioning is required for the linear solver to be effective. Right preconditioning of the form (JM−1 )(Mx) = b is used here where M is the preconditioning matrix. An additive Schwarz global preconditioner with variable overlap [25,27,28] is used in conjunction with local BFILU preconditioners for each sub-domain. The local preconditioner is based on a block-fill ILU(f ) or BFILU(f ) factorization of an approximate Jacobian for each subdomain. Here, f is the level of fill. This combination of preconditioning fits well with the block-based AMR described by Northrup and Groth [11] and is compatible with domain decomposition methods, readily enabling parallel implementation of the overall Newton method. Rather efficient parallel implementations of implicit algorithms via Schwarz preconditioning have been developed by Keyes and co-researchers and successfully applied to the prediction of transonic full potential, low-Mach-number compressible combusting, and three-dimensional inviscid flows [27, 29, 30]. As the GMRES algorithm does not explicitly require the evaluation of the global Jacobian matrix, J, a so-called “matrix-free” or “Jacobian-free” approach can be adopted and is used here. Numerical differentiation based on Fréchet

derivatives is used to approximate the matrix-vector product JM−1 x as follows: JM−1 x ≈

R(U + εM−1 x) − R(U) 3M−1 x + , ε 24t

(8)

where R(U+εM−1 x) is the residual vector evaluated at some perturbed solution state and ε is a small scalar quantity. Although the performance of the Jacobian-free method is sensitive to the choice of ε, Neilsen et al. [31] have found 1/2 that ε = ε◦ /||x||2 seems to work well, with ε◦ ≈ 10−8 –10−7 . 4

Numerical Results

The proposed parallel implicit AMR algorithm is now applied to the solution of a forced, time-dependent, axisymmetric, methane-air coflow laminar diffusion flame. In particular, the solution of the flame studied by Day and Bell [32], and Dworkin et al. [33] is considered. The flame boundary and initial conditions are the same as those used in the previous studies. The computational domain is rectangular in shape with dimensions of 10 cm by 5 cm. The axis of symmetry is aligned with the left boundary of the domain and the right far-field boundary is taken to be a freeslip boundary along which inviscid reflection boundary data is specified. The top or outlet of the flow domain is open to a stagnant reservoir at atmospheric pressure and temperature and Neumann-type boundary conditions are applied to all properties except pressure which is held constant. The bottom or inlet is subdivided into four regions. The innermost region (r ≤ 2 mm) is the fuel inlet or jet, which injects a nitrogen diluted methane fuel mixture (cCH4 = 0.5149, cN2 = 0.4851, cO2 = 0, cCO2 = 0, cCO = 0, and cH2 O = 0) at 298 K with a parabolic axial velocity profile. The time-variation in the flame is produced by imposing a sinusoidal axial velocity fluctuation across the fuel jet, vz = 70.0(1 − r2 /R)(1 + α sin ωt) cm/s where α is the velocity amplitude and ω is the frequency of oscillation. Matching the experiments of Dworkin et al. [33], the velocity amplitude was set at 50% at a frequency of 20Hz. The next region (2 mm < r ≤ 2.38 mm) is a small gap associated with the annular wall separating the fuel and oxidizer. The third region (2.38 mm < r ≤ 2.50 cm) is the co-flowing oxidizer, in this case air at 298K (cO2 = 0.232, cN2 = 0.768, cCH4 = 0, cCO2 = 0, cCO = 0, and cH2 O = 0), with a uniform velocity profile of 0.35 m/s. The final outer region of the lower boundary (2.5 cm < r ≤ 5 cm) is again a far-field boundary along which free-slip boundary conditions are applied. The solution domain is initialized with a uniform solution state corresponding to quiescent air at 298K, except for a thin region across the fuel and oxidizer inlets, which is taken to be air at 1500 K so as to ignite the flame. Additional details concerning the setup for this diffusion flame can be found in the papers by Day and Bell [32], and Dworkin et al. [33] The flame was calculated by first letting the steady-state flame develop with the fuel flow velocity held constant at 70 cm/s, then four-full periods of the flame oscillation were calculated to avoid any non-periodic oscillations created during start-up. Newton’s method was used with a GMRES tolerance of 0.1, ILU(4) preconditioner and at each time step the Newton iterations were converged two orders of magnitude, to a maximum of 10 Newton steps, where 5-7 is typical. The approximate Jacobian preconditioner was only updated for the first Newton step of each time-step, unless the number of GMRES iterations required increased. A fixed time-step of 0.0001s was used during the unsteady calculation. Mesh refinement was carried out every 0.0025s based on the gradient of temperature, and 3 levels of mesh refinement (4 mesh levels) were used as shown in Figure 1. Figure 2 shows the resulting isotherms at five times through the flames periodic fluctuation of 0.05s. The flames shape and structure match well with the experimental results of Dworkin et al. [33]. The most significant difference is in the over prediction of temperature, however this is related to the use of the simplified non-reversible one-step reaction mechanism used to model the methane-air chemistry. 5

Conclusions

A parallel implicit AMR scheme has been developed for solving unsteady laminar combusting flows. The combination of finite-volume discretization procedure, parallel block-based AMR strategy, low-Mach-number preconditioning, dual-time-stepping-like apporach, and Newton method solution procedure has resulted in a powerful computational tool for predicting a wide range of unsteady laminar reactive flows, from compressible to nearly incompressible lowMach-number regimes. The validity and performance of the method has been demonstrated for non-premixed flames. Future work will involve the extension of the the approach to unsteady three-dimensional combusting flows.

a) 0.01s

b) 0.03s

c) 0.05s

0.1

0.08

0.08

0.06

z (m)

z (m)

0.06

0.04

0.04

0.02

0.02

0 -0.02

-0.01

0

0.01

0.02

0 -0.02

-0.01

0

r (m)

0.01

0.02

-0.02

-0.01

0

r (m)

0.01

0.02

r (m)

Fig. 1: Time-varying methane-air co-flow laminar diffusion flame isotherms at three intervals calculated with a) 456, b) 444, and c) 447 (4x8 cell) blocks each with 3 levels of refinement. a) 0.01s

b) 0.02s

c) 0.03s

0.08

0

r (m)

0.005

0.01

z (m)

z (m)

0.06

z (m)

0.06

z (m)

0.06

0.04

0.04

0.04

0.04

0.02

0.02

0.02

0.02

0 -0.01

-0.005

0

r (m)

0.005

0.01

0 -0.01

-0.005

0

r (m)

0.005

0.01

T (K)

0.08

0.08

0.06

0 -0.005

e) 0.05s

d) 0.04s

0.08

2200 2162 2124 2087 2049 2011 1973 1936 1898 1860 1822 1785 1747 1709 1671 1634 1596 1558 1520 1483 1445 1407 1369 1332 1294 1256 1218 1181 1143 1105 1067 1030 992 954 916 879 841 803 765 728 690 652 614 577 539 501 463 426 388 350

0 -0.01

-0.005

0

r (m)

0.005

0.01

-0.01

-0.005

0

0.005

0.01

r (m)

Fig. 2: Time-varying methane-air co-flow laminar diffusion flame isotherms at five 0.01ms, (a-e) intervals.

References [1] [2] [3] [4] [5] [6] [7] [8]

[9] [10] [11]

[12]

[13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33]

M. J. Berger. J. Comput. Phys., 53:484–512, 1984. M. J. Berger and P. Colella. J. Comput. Phys., 82:67–84, 1989. D. De Zeeuw and K. G. Powell. J. Comput. Phys., 104:56–68, 1993. M. J. Aftomis, M. J. Berger, and J. E. Melton. AIAA J., 36(6):952–960, 1998. C. P. T. Groth, D. L. De Zeeuw, K. G. Powell, T. I. Gombosi, and Q. F. Stout. A parallel solution-adaptive scheme for ideal magnetohydrodynamics. Paper 99-3273, AIAA, June 1999. J. S. Sachdev, C. P. T. Groth, and J. J. Gottlieb. Int. J. Comput. Fluid Dyn., 19(2):157–175, 2005. J.B. Bell, M.S. Day, A.S. Almgren, M. J. Lijewski, and C. A. Rendleman. Int. J. Numer. Meth. Fluids, 40:209–216, 2002. C. C. Douglas, A. Ern, and M. D. Smooke. Numerical simulation of flames using multigrid methods. In J. Wang, M. B. Allen, B. M. Chen, and T. Mathew, editors, Iterative Methods in Scientific Computation, volume 4 of IMACS Series in Computational and Applied Mathematics, pages 149–154. New Brunswick, 1998. S. A. Northrup and C. P. T. Groth. Solution of laminar diffusion flames using a parallel adaptive mesh refinement algorithm. Paper 2005–0547, AIAA, January 2005. X. Gao and C. P. T. Groth. Parallel adaptive mesh refinement scheme for turbulent non-premixed combusting flow prediction. Paper 2006–1448, AIAA, January 2006. S. A. Northrup and C. P. T. Groth. Solution of laminar combusting flows using a parallel implicit adaptive mesh refinement algorithm. In Proceedings of the 4th International Conference on Computational Fluid Dynamics, Ghent, Belgium, July 10–14, 2006, 2006. S. A. Northrup and C. P. T. Groth. Parallel implicit adaptive mesh refinement algorithm for solution of laminar combusing flows. In Proceedings of the 14th Annual Conference of the CFD Society of Canada, Kingston, Canada, July 16–18, 2006, pages 221–228. CFD Society of Canada, 2006. A. Jameson. Time dependent calculations using multigrid, with applications to unsteady flows past airfoils and wings. Paper 1991-1596, AIAA, June 1991. S. Gordon and B. J. McBride. Computer program for calculation of complex chemical equilibrium compositions and applications I. analysis. Reference Publication 1311, NASA, 1994. B. J. McBride and S. Gordon. Computer program for calculation of complex chemical equilibrium compositions and applications ii. users manual and program description. Reference Publication 1311, NASA, 1996. C. R. Wilke. J. Chem. Phys., 18:517–519, 1950. G. Dixon-Lewis. Computer modeling of combustion reactions in flowing systems with transport. In W. C. Gardiner, editor, Combustion Chemistry, pages 21–126. Springer-Verlag, New York, 1984. C. K. Westbrook and F. L. Dryer. Combust. Sci. Tech., 27:31–43, 1981. J. M. Weiss and W. A. Smith. AIAA J., 33(11):2050–2057, 1995. T. J. Barth. Recent developments in high order k-exact reconstruction on unstructured meshes. Paper 93-0668, AIAA, January 1993. V. Venkatakrishnan. On the accuracy of limiters and convergence to steady state solutions. Paper 93-0880, AIAA, January 1993. P. L. Roe. J. Comput. Phys., 43:357–372, 1981. W. J. Coirier and K. G. Powell. AIAA J., 34(5):938–945, May 1996. Y. Saad and M. H. Schultz. SIAM J. Sci. Stat. Comput., 7(3):856–869, 1986. Y. Saad. Iterative Methods for Sparse Linear Systems. PWS Publishing Company, Boston, 1996. R. S. Dembo, S. C. Eisenstat, and T. Steihaug. SIAM J. Numer. Anal., 19(2):400–408, 1982. W. D. Gropp, D. K. Kaushik, D. E. Keyes, and B. F. Smith. Parallel Computing, 27:337–362, 2001. D. A. Knoll and D. E. Keyes. J. Comput. Phys., 193:257–397, 2004. D. A. Knoll, P. R. McHugh, and D. E. Keyes. AIAA J., 34(5):961–967, 1996. X.-C. Cai, W. D. Gropp, D. E. Keyes, R. G. Melvin, and D. P. Young. SIAM J. Sci. Comput., 19(1):246–265, 1998. E. J. Nielsen, W. K. Anderson, R. W. Walters, and D. E. Keyes. Application of Newton-Krylov methodology to a threedimensional unstructured Euler code. Paper 95-1733-CP, AIAA, June 1995. M. S. Day and J. B. Bell. Combust. Theory Modelling, 4(4):535–556, 2000. S. B. Dworkin, B. C. Connelly, A. M. Schaffer, B. A. V. Bennett, M. B. Long, M. D. Smooke, M. P. Puccio, B. McAndrews, and J. H. Miller. Computational and experimental study of a forced, time-dependent, methane-air coflow diffusion flame. In Thirty First Symposium (International) on Combustion, pages 971–978, Pittsburgh, 2007. The Combustion Institute.