当前位置:首页 >> >>

IMPLICIT-EXPLICIT RUNGE-KUTTA METHOD FOR COMBUSTION SIMULATION


European Conference on Computational Fluid Dynamics ECCOMAS CFD 2006 P. Wesseling, E. O?ate and J. P?riaux (Eds) n e c TU Delft, The Netherlands, 2006

?

IMPLICIT-EXPLICIT RUNGE-KUTTA METHOD FOR COMBUSTION SIMULATION
u o E. Lindblad1 2 , D.M. Valiev2 , B. M¨ller1 , J. Rantakokko1 3 , P. L¨tstedt1 , M.A. Liberman2
1 Department

of Information Technology, Uppsala University, Box 337, 751 05 Uppsala, Sweden e-mail: erik.lindblad@it.uu.se web page: http://www.it.uu.se/katalog/erikl 2 Department of Physics, Uppsala University, Box 530, 751 21 Uppsala, Sweden e-mail: damir.valiev@fysik.uu.se 3 UPPMAX, Uppsala University, Box 337, 751 05 Uppsala, Sweden web page: http://www.uppmax.uu.se

Key words: Semi-implicit methods, sti? problems, combustion, de?agration-to-detonation transition Abstract. New high order implicit-explicit Runge-Kutta methods have been developed and implemented into a ?nite volume code to solve the Navier-Stokes equations for reacting gas mixtures. The resulting nonlinear systems in each stage are solved by Newton’s method. If only the chemistry is treated implicitly, the linear systems in each Newton iteration are simple and solved directly. If in addition certain convection or di?usion terms are treated implicitly as well, the sparse linear systems in each Newton iteration are solved by preconditioned GMRES. Numerical simulations of de?agration-to-detonation transition (DDT) show the potential of the new time integration for computaional combustion.

1

INTRODUCTION

The distinctive feature of premixed combustion is its ability to propagate as a selfsustained wave of the exothermic chemical reaction spreading through a homogeneous combustible mixture either as a subsonic de?agration (premixed ?ame) or supersonic detonation. Thus, both de?agration and detonation appear to be stable attractors each being linked to its own base of initial data. In uncon?ned obstacle-free systems the concrete realization of the speci?c propagation mode is controlled by the ignition conditions. Normally, de?agrations are initiated by a mild energy discharge, i.e. by a spark, while detonations are provoked by shock waves via localized explosions. It is known, however, that in the presence of obstacles or con?nement (tube walls, wire screens, porous matrix, etc.) the initially formed de?agration undergoes gradual acceleration abruptly converting into detonation [1, 2, 3, 4, 5]. In spite

1

E. Lindblad, D.M. Valiev, B. M¨ ller, J. Rantakokko, P. L¨tstedt and M.A. Liberman u o

of numerous e?orts, the basic mechanisms controlling the spontaneous transition from de?agrative to detonative combustion (DDT) is still remaining a major unsolved challenge of the combustion theory [6, 7, 8, 9]. The classical explanation of DDT [1, 2, 3] was due to ?ame acceleration in a tube caused by the ?ame interaction with the nonuniform ?ow ahead of the ?ame front, which is formed by the ?ame pushing the unburned gas in a tube with adhesive walls. The nonuniform ?ow stretches the ?ame, so that the shape of the ?ame mimics the velocity pro?le in the ?ow ahead. This increases the ?ame surface area, thus accelerating the reaction wave propagation and the ?ow. This ?ame-acceleration phenomenon is presumably enhanced by the ?ame interaction with turbulence generated in the boundary layer formed in the ?ow ahead of the ?ame [1, 2, 10]. Indeed, it has recently been realized [11] that the hydraulic resistance alone is capable of triggering the transition even if the multi-dimensional e?ects, such as the ?ame acceleration due to folding, are completely suppressed and the system is regarded as e?ectively one-dimensional. The basic predictions of the one-dimensional model were recently corroborated by direct numerical simulations of premixed gas combustion in thin channels, where the hydraulic resistance is incorporated through the no-slip boundary condition rather than through the volumetric drag-force [12]. It was also recently shown [13, 14] that a similar e?ect is observed in wide channels and in the channels with rough walls. For both no-slip and rough walls the transition is triggered predominantly by hydraulic resistance inducing formation of an extended preheat zone ahead of the advancing ?ame, and thereby creating conditions pertinent to Zeldovich’s mechanism of soft initiation. The detonation ?rst develops in the near-wall mixture adjacent to the ?ame, corroborating many experimental observations (e.g. [7]). De?agration-to-detonation transition in uncon?ned systems is more problematic. There are reports claiming that in highly sensitive oxygen-based mixtures the transition may be triggered by outwardly propagating ’free-space’ ?ames [15]. In this description, the transition is commonly attributed to the ?ame acceleration induced by the Darrieus-Landau (DL) instability (spontaneous ?ame wrinkling). Yet, the acceleration resulting from the wrinkling seems to be a rather weak e?ect whose ability to cause the transition is not at all obvious. In the foreword to Nettleton’s monograph on gaseous detonation [3], when discussing the problem of transition, Zel’dovich wrote, The role of the internal instability of the plane slow ?ame (Landau, Darrieus) is still not clear. As is now well known, (e.g. [16, 17, 18, 19, 20, 21]) in wide channels the DL instability results in the formation of wrinkled ?ames and the ?ame speed enhancement due to the increase of the ?ame area. The wrinkled ?ame generates a shock with a Mach number of about 1.2 - 1.5, which is too low to trigger detonation. Yet, as has recently been shown [13, 14, 22, 23], there is another previously overlooked aspect of the DL instability. The folded reaction zone creates a low-gradient preheating (preconditioning) of the fresh mixture trapped within the fold interior. This, under favorable conditions, may invoke autoignition triggering the transition. The mechanism of the transition is the temperature increase 2

E. Lindblad, D.M. Valiev, B. M¨ ller, J. Rantakokko, P. L¨tstedt and M.A. Liberman u o

due to the in?ux of heat from the folded reaction zone, followed by autoignition. The transition occurs when the pressure elevation at the accelerating reaction front becomes high enough to produce a shock capable of supporting detonation. This requires the fold to be su?ciently narrow and deep. The e?ect was found to be sensitive to the ?ame’s normal speed and the reaction rate pressure-dependency, favoring fast ?ames and highorder reactions. In the context of colliding elliptic ?ames the folding induced transition was reported also in [24] and for ?ame initiated by corrugated walls in [25]. A central issue in simulations of combustion is to ensure that the computations resolve all the most important ?ow and chemical time and space scales. To simulate combustion processes from ?rst principles, it is necessary to resolve the relevant scales from the size of the system to the ?ame thickness, a range that can cover up to twelve orders of magnitude. This computational challenge in the development of numerical algorithms for solving coupled partial and ordinary di?erential equations resulted in the development of several numerical methods, including adaptive mesh re?nement to deal with multiscale phenomena, domain decomposition, and multiresolution methods using wavelets. The present work is aimed at gaining further improvement of the numerical code performance for modeling complex combustion processes. The present simulations were performed using a parallel version of the code developed by L.-E. Eriksson [26]. This code solves the Navier-Stokes equations for reacting gas mixtures using a third-order upwindbiased ?nite volume method for the inviscid ?uxes and a second-order central discretization of the viscous ?uxes with an explicit second order Runge-Kutta time integration method. The code has been successfully used for solving physics problems of ?ame instability in wide tubes [16, 18, 20] and gaining deeper understanding of DDT in wide tubes with thermoisolated (adiabatic boundary conditions) and rough walls [13, 14, 22, 23]. Further steps in the DDT studies will include investigation of heat losses to the walls, in?uence of complex chemistry and ?ame-turbulence interaction, and simulations in 3D. These additions increase the sti?ness of the governing equations and therefore the time stepping method must be improved. In the present study, new time integration methods have been implemented into the combustion code, namely second and fourth order implicit-explicit Runge-Kutta methods [28] [29], as well as a third order implicit-explicit Runge-Kutta method [36]. Whereas the inviscid and viscous ?uxes are treated explictly, the chemistry is treated implicitly. Since only the species continuity equations have nonzero chemistry terms, the resulting nonlinear systems only involve the species densities and are quickly solved by Newton’s method. Other classi?cations of non-sti? and sti? parts of the equations have been investigated as well. Then, the resulting nonlinear systems are solved using e?cient Jacobian matrix calculations and GMRES. The systems are preconditioned using incomplete LU factorizations. Preliminary results show that the increased stability of the implicit method combined with the e?ciency of the explicit method will be an e?cient solver for the intended combustion problems.

3

E. Lindblad, D.M. Valiev, B. M¨ ller, J. Rantakokko, P. L¨tstedt and M.A. Liberman u o

2 2.1

NAVIER-STOKES EQUATIONS FOR REACTING GAS Thermodynamic, Transport and Chemical Models

We consider a thermally perfect gas mixture of ? species ? with the molecular weights ? and the densities . The mass fraction of species (meaning ? ) is the ratio of the mass of species and the mass of the mixture, i.e. = (1)

where is the density of the mixture. Details of the thermodynamic model we used are given by Kuo [30]. Viscosity, molecular di?usion and heat conduction are described by simpli?ed transport models. For the viscosity , we use the Sutherland law. The Schmidt numbers

? =

(2)

denotes the di?usion coe?cient of all species are assumed to be equal and constant. of species ? and is determined from equation (2) by the speci?ed Schmidt number and the computed density and viscosity. The Prandtl number

?? =

?

(3)

is assumed to be constant. ? and are the speci?c heat at constant pressure and the heat conduction coe?cient of the mixture, respectively. We consider a chemical reaction mechanism with ? reactions [30]. A reaction (here is an index for reactions and not the enthalpy) is of the form
?
=1

? ?

?

?
=1

?? ?

(4)

where ? and ?? are the stoichiometric coe?cients of reaction for species appearing as a reactant and as a product, respectively. The reaction rate of reaction is
?
=1

?

?

(5)

with

? the molecular weight of species and the Arrhenius term
=

? ? ??
4

?

?? ?

(6)

E. Lindblad, D.M. Valiev, B. M¨ ller, J. Rantakokko, P. L¨tstedt and M.A. Liberman u o

where and ? are constants and is the activation energy of reaction . ?? is the universal gas constant, and ? is the temperature. The backward reaction of reaction (5) is
?

??

?
The equilibrium constants
?
=1

=1

?
?? ? ?

(7)

? =

(

)

?

?? ? ? ??

(8)

are functions of temperature, and tables exist for most of them [31]. With the de?nitions above, the production rate of species is given by =? 2.2
?

( ??

=1

?

? )

?
=1

?

?

??

1

?
=1

??

?

(9)

Navier-Stokes equations

In Cartesian coordinates, the 2D compressible Navier-Stokes equations for reacting gas ?ow read in conservative form [30, 17]

 U  (F ? F? )  (G ? G? ) + + = S ? ? ?

(10)

where U is the vector of the conservative variables, F = F1 and G = F2 are the inviscid ?ux vectors for the ?- and ? -directions, and F? = F?1 and Gv = F?2 are the viscous ?ux vectors for the ?- and ? -directions. Let the Cartesian coordinates and velocity components be denoted by (?1 ?2 )? = (? ? )? and (?1 ?2 )? = (? ? )? . The pressure, the total energy per unit mass, the total enthalpy, and the enthalpy of species are denoted by ?, , ? , and , respectively. Then the conservative variables, the inviscid and viscous ?ux vectors are

?

?

?

?1 ?2
U= .  . .
1

F = .  . .

? ?1 ? + ??1 ?2 ? + ??2 ?? 1?
??1 ?

?

?

0
1

?

F? =

?2
. .  .

2

?? ? ?=1 ?? ? + ? + =1  1 1 ?
??1  ??1 ?

 ?

??1

(11)

5

E. Lindblad, D.M. Valiev, B. M¨ ller, J. Rantakokko, P. L¨tstedt and M.A. Liberman u o

and S is the source term S? = (0 0 0 0
1

??1 )

For a Newtonian ?uid, the components of the shear stress tensor are =

? ? + ? ?

?2 3

? =1 ?
2

?

where ? = 1 if = and ? = 0 if = . Integrating equation (10) over a control volume ? (actually a surface in 2D) with the boundary  ? and the outer normal unit vector n = (?? ?? )? and using the Gauss theorem, we obtain the integral form of the 2D Navier-Stokes equations for a reacting gas ?ow

U ? + ? ?
3

?

(F ? F? )??

+

?

(G ? G? )??

=

?

S

?

(12)

FINITE VOLUME METHOD

With the ?nite volume method, equation (12) is discretized for each grid cell by approximating U

?

? ?? +

?

×=1

[(F ? F? )??

+ (G ? G? )?? ]× = S

? ??

(13)

where U is the volume averaged vector of the conservative variables in the cell ? and S the volume averaged source vector. ? ?? is the area of the cell. Since we consider structured grids with quadrilaterals as control volumes, the cells have ? = 4 sides. × is the length of the cell interface ×. The cell averages U are the unknowns in the cell-centered ?nite volume method. Therefore, we have to approximate the ?ux vectors at the cell interfaces. The inviscid ?ux vectors are discretized by a third-order upwind-biased approximation of the characteristic variables and using a total variation diminishing (TVD) limiter [26] [27]. Central discretizations are employed for the viscous ?uxes at the cell interfaces. The volume averaged nonlinear source term is approximated by S S(U ) (14)

After the ?nite volume discretization of equation (13), we have a system of ordinary di?erential equations (ODEs) for the time dependent cell averages U

? ? = (? ) + (? ) (15) where ? denotes the vector of all U , (? ) the vector of all inviscid and viscous ?ux ?? 1 ? ? ? ? contributions, i.e. all ? ? ?? ×=1 [(F ? F? )?? + (G ? G? )?? ]× , and (? ) the vector
6

E. Lindblad, D.M. Valiev, B. M¨ ller, J. Rantakokko, P. L¨tstedt and M.A. Liberman u o

of all source terms S . Thereby, we classify the right hand side of the ODE (15) into a non-sti? part and a sti? part . Other classi?cations are possible and discussed in section 6.2. 4 IMPLICIT-EXPLICIT RUNGE-KUTTA METHODS

We have developed high order implicit-explicit Runge-Kutta (IERK) methods to e?ciently solve separable sti? problems (15) [28, 29]. While ?rst-order IERK methods have frequently been used to treat the sti? source terms from chemistry implicitly in hypersonic ?ow and combustion simulations [32, 33], two of our new IERK methods are second order accurate and one is fourth order accurate. Similar high order IERK methods have only recently been available [34, 35, 36, 37, 38]. An explicit Runge-Kutta (ERK) method is used to solve the non-sti? part and a diagonally implicit Runge-Kutta (DIRK) method is employed to solve the sti? part . A general ×-stage implicit-explicit Runge-Kutta (IERK) method consists of an ×-stage ERK = 1 ×. The and an ×-stage DIRK method with common weighting coe?cients following tableaus de?ne the ERK and DIRK methods of an IERK method [28, 29]: 0
21 11

0

21

22

×1
1

× 1 ?+1 at ? = (? + 1)?? is de?ned by The approximate solution ? × ? ? ?+1 = ? ? + ?? where

? ?

0
× ×?1

?

0

×1

? ? ? ?

(16)
×× ×

= =1

? ? + ?? ×

=1 1 ?

?

+

=1

? ? + ??

?
=1

(17)

The coe?cients of our 4th order accurate 5-stage IERK method denoted as IERK45 method are given in Table 1. 5 PARALLELIZATION

The original code is parallelized using the Message Passing Interface (MPI). Since the solver is designed to handle multi-block grids the processors are ?rst divided into clusters, one cluster for each sub-block. The number of processors in each cluster are chosen to optimize the load balance using the clustering algorithm in [39]. The blocks are then partitioned in two dimensions within the corresponding processor cluster. Almost all calculations can be performed locally within the processors except for the ?ux calculations which require data from its nearest neighbors. The exchange of data is handled with outof-order non-blocking communication using MPI’s persistent communication objects. 7

E. Lindblad, D.M. Valiev, B. M¨ ller, J. Rantakokko, P. L¨tstedt and M.A. Liberman u o

Table 1: Coe?cients for fourth order implicit-explicit Runge-Kutta method IERK45

?21 = 0 39098372452428 ?31 = 1 09436646160460 ?32 = 0 33181504274704 ?41 = 0 14631668003312 ?42 = 0 69488738277516 ?43 = 0 46893381306619 ?51 = ?1 33389883143642 ?52 = 2 90509214801204 ?53 = ?1 06511748457024 ?54 = 0 27210900509137
1 2 3 4 5

11 21 22 31 32 33 41 42 43 44 51 52 53 54 55

= = = = =

51 52 53 54 55

=1 4 = 0 34114705729739 =1 4 = 0 80458720789763 = ?0 07095262154540 =1 4 = ?0 52932607329103 = 1 15137638494253 = ?0 80248263237803 =1 4 = 0 11933093090075 = 0 55125531344927 = ?0 1216872844994 = 0 20110104014943 =1 4

The solver has been tested on a Sun Fire 15k server with UltraSparcIIIcu processors running at 900MHz. The largest partition of the machine has 36 processors and 36GB RAM. Each processor has 64KB L1-cache and 8MB L2-cache. We have run several test cases with di?erent problem sizes (600x200 to 1000x200 grid points) and di?erent numbers of blocks (1 to 5 blocks in the multi-block case). All cases show excellent scaling with superlinear speedup up to full machine size, i.e., speedup close to 40 on 36 processors, cf. Figure 1 [40]. The superlinear speedup is due to better cache utilization as data is divided into smaller pieces when running on multiple processors. 6 6.1 RESULTS De?agration to Detonation Transition (DDT)

One of the typical pictures of the transition due to formation of the appropriate ?ame fold is shown in Fig. 2, where the sequence of zonal images depicting the square of pressure gradient gives a cinematic impression of the dynamics of the ?ame front during the incipient stage of the transition from de?agration to detonation [14]. These images resemble the schlieren photographs of laboratory experiments, though the latter visualize gradients of the density rather than of the pressure. The numerical resolution is 10 cells per ?ame width. The time instants of Fig. 2 are not evenly spaced but rather clustered around the transition point. The earlier images depicting the incipient phase of the ?ame evolution in the vicinity of the tube’s closed end, are not shown. 8

E. Lindblad, D.M. Valiev, B. M¨ ller, J. Rantakokko, P. L¨tstedt and M.A. Liberman u o

Figure 1: Speedup versus number of processors for 3 blocks running on a SunFire 15k, ? ideal speedup, ? 1D decomposition of each block separately, 2D decomposition of each block separately, ?? dividing the processors among the blocks [40]

A zoomed view of the ?ame fold dynamics near the transition point is shown in Fig. 3 [14]. The associated pro?les of temperature, density, ?ow velocity, and pressure along the fold axis are plotted on Fig. 4 [14]. Here one readily observes (i) formation of the large-scale preheat zone (preconditioning) in the unreacted gas trapped within the fold interior, (ii) acceleration of the fold-tip, and (iii) the pressure elevation and formation of a high pressure peak. The transition occurs when the pressure peak becomes high enough to produce a shock capable of supporting detonation. This requires the fold to be narrow and deep enough; otherwise one ends up with a moderately strong pressure wave insu?cient for triggering the transition. This mechanism of transition by an appropriate non-uniformity in the temperature ?eld (preconditioning) may naturally be associated with Zel’dovich’s theory of soft initiation [41, 42]. 6.2 Test of IERK Methods

To show the advantages of using IERK methods for combustion a series of DDT simulations have been performed. The ?ux is calculated using the combustion code of [26]. The original code uses a second order explicit Runge-Kutta method known as Gary’s method [43] for the time stepping. For our tests it has been replaced by implicit-explicit Runge-Kutta methods of order 2, 3 and 4. The third order method is derived in [36], and the second and fourth order methods are derived in [28, 29]. Simulations showing order of accuracy and applications of IERK methods to various model problems are also presented in the cited source papers.

9

E. Lindblad, D.M. Valiev, B. M¨ ller, J. Rantakokko, P. L¨tstedt and M.A. Liberman u o

6.2.1

Static implicit treatment of chemistry

As stated above the code by [26] has been successfully used for combustion calculations. The need for implicit treatment of the chemistry becomes evident when studying the transition from de?agration to detonation for a propagating ?ame. Fig. 5 shows the temperature contours of the DDT simulation obtained using the original explicit code. Soon after the transition the solution becomes unstable. Since the chemical source term (14) only depends on cell data, the implicit iteration can be performed by solving a (3 + ?) ? (3 + ?) system per cell where ? is the number of species. In the present simulations a mixture of two species representing burned and unburned fuel is used. The production rate of the fuel is given by an Arrhenius expression, cf. (9) and (6).
1

=?

??1

??

?

?? ?

(18)

is the mass fraction of the fuel, the reaction rate constant, and the where activation energy. Here ? denotes the reaction order. The temperature ? is dependent on and given by

?=

? ?0 ? 1 (?2 + ?2) 2
?

(19)

where is the total energy per unit mass, is the density and ? and ? are velocities in ? and ? directions, respectively. The enthalpy of formation ?0 and the speci?c heat at constant volume ? for a mixture of ? species are given by

?0 =
?=

?
=1

?0
?

(20) (21)
?1 = ?2

?
=1

?0 and ? are corresponding values for species . For the present case ? = 2, and ?02 = 0 and (19) thus reduces to ?=

? ?0 ? 1 (?2 + ?2 ) 2
1

?1

(22)

This reduces the implicit treatment of the source term (18) to a scalar equation which can be solved using scalar Newton iteration without having to solve any linear systems or compute Jacobians. The overhead for solving the chemistry implicitly is therefore very small. Fig. 6 shows the temperature contours obtained for the same setup as the explicit case above but using the ?rst order implicit Euler method for the chemistry treatment. Here the solution remains stable during the transition. 10

E. Lindblad, D.M. Valiev, B. M¨ ller, J. Rantakokko, P. L¨tstedt and M.A. Liberman u o

The DDT for the implicit-explicit code occurs after roughly twice the time as for the explicit code. This is because of the chaotic behaviour of the perturbations on the ?ame front that cause the cusps and thereby trigger the DDT. Therefore direct comparisons cannot be made between the two cases, but the main observation is that the explicit code becomes unstable shortly after the DDT when choosing the time step based on the stability conditions for the inviscid and viscous parts. 6.2.2 Flexible treatment of ?ux terms using switch technique

Although the implicit treatment of the chemistry term proves e?cient to retain stability during DDT, apart from being more e?cient than fully implicit methods, the real bene?ts come from being able to con?gure the explicit and implicit treatment more ?exibly. Consider an ODE or a semi-discretized PDE consisting of several separable terms:

? ? (?) =

?
=1

(?

?)

(23)

For a given problem some of the terms will be non-sti? and some will be sti?. Every will yield an upper limit on the time step for an explicit method to remain stable, the sti?er the term the smaller the limiting time step. When using an explicit method the most restrictive time step will be used for all terms, resulting in unnecessary evaluations of the other terms. When the sti?ness varies signi?cantly among the terms the use of IERK methods can be more e?cient than fully explicit or implicit methods, as shown above. The sti?ness can vary in the course of a simulation, both in time and space. For combustion simulations the ?ux gradients are steepest near the ?ame, which moves continuously. Depending on the speci?c problem cell shapes, grid con?gurations and boundary layer term will yield the most restrictvie time step. e?ects greatly a?ect which For these situations it is bene?cial to be able to freely choose which terms are treated implicitly. In the Navier-Stokes example the consists of three terms, inviscid ?ux, viscous ?ux and chemical source term. In the next example, a ?ame propagating in a tube similar to the case described above, the limiting time step for the inviscid and viscous ?uxes are calculated in each step, the chemical source term is always treated implicitly. For this problem the viscous terms yield a more than 100 times smaller limiting time step than the inviscid terms. This small step size is not needed to maintain accuracy and therefore we choose to treat also the viscous terms implicitly. This enables us to use a time step limited only by the inviscid terms. To investigate the impact of a larger time step a reference solution is computed explicitly using a time step ??? = 4 ? 10?11 s, calculated from the stability condition of the viscous terms. This is compared with the solutions obtained when using a time step ??? ?? = 4 ? 10?9 and treating the viscous terms and chemistry implicitly.

11

E. Lindblad, D.M. Valiev, B. M¨ ller, J. Rantakokko, P. L¨tstedt and M.A. Liberman u o

The nonlinear equations in each stage of the IERK methods are solved by Newton’s method. The Jacobian matrices are calculated by a variant of TOMS 618 [44]. The linear systems in each Newton iteration are solved by the generalized minimum residual method (GMRES) with a restart vector length of 20. Krylov methods like GMRES are well suited for solving large, sparse linear systems [45, 46, 47]. The implementation used is from [48]. To increase the convergence rate of the GMRES iterations preconditioning is used. Several types of incomplete LU-factorizations (ILU) have been tested and the simulations discussed here use ILU0 preconditioning from SPARSKIT [49, 50]. The ILU0 preconditioning is a very simple technique where no new non-zero elements are introduced during the ILU factorization, thereby allowing for e?cient storage of the sparse matrices. Fig. 7 shows the density contours at ? = 1000??? , ? = 5000??? and ? = 10000??? . The results obtained with IERK45 are in good agreement with those obtained with the explicit second order Runge-Kutta method (ERK2). Although the computational work per time level is larger, IERK45 is more e?cient than ERK2, because the time step for IERK45 can be chosen 100 times larger than for ERK2. The original code is set to recalculate the grid as the wave propagates. The grid update is dependent on the current solution and therefore direct comparisons of the solutions would be misleading. The di?erences in the ?nal solutions would originate both from the solution interpolation and from the longer time step and a direct comparison would not be possible. Therefore a ?xed grid was used to obtain the results shown in ?gure 7. 6.2.3 Dynamical con?guration of explicit and implicit treatment

A general approach to deal with the varying sti?ness in time and space is to use adaptive mesh re?nement, AMR, to concentrate the grid resolution on areas where it is needed. terms are treated explicitly To give even more e?ciency the con?guration of which and which are treated implicitly can be changed dynamically during a simulation. The same mechanisms for deciding when to re?ne or coarsen the grid resolution can be used to switch implicit treatment of di?erent terms on and o?. In this way the expensive implicit solvers are only used where absolutely needed. 7 CONCLUSIONS

A new time discretization has been developed for integration of the equations of combustion. By treating parts of the unknown variables implicitly and the remaining parts explicitly in a Runge-Kutta method, signi?cant savings in computing time are possible. The chemistry equations in each cell are advanced by an implicit time-stepping in one example. The time step restrictions in an explicit method due to the sti?ness in the chemistry model are then avoided. In another model example, the viscous and chemistry terms are integrated implicitly and compared to explicit integration for all components. The resulting CFL number with the implicit-explicit method is up to 100 times larger than with the explicit scheme without hampering the accuracy in the solution. 12

E. Lindblad, D.M. Valiev, B. M¨ ller, J. Rantakokko, P. L¨tstedt and M.A. Liberman u o

ACKNOWLEDGEMENTS Fellowships for E. Lindblad by EU and G¨ran Gustafsson Foundation are gratefully o acknowledged by M. Liberman and B. M¨ ller, respectively. u REFERENCES [1] Schelkin, K. I. and Troshin, Ya. K., Gasdynamics of Combustion, Mono Book Corp. Baltimore, 1965. [2] Zel’dovich, Ya. B., Kompaneets, A.S., Theory of Detonation, Academic press, New York, 1960. [3] Nettleton, M. A., Gaseous Detonation, Chapman and Hall, London, 1987. [4] Urtiew, P., and Oppenheim, A. K., Experimental observation of the transition to detonation in an explosive gas, Proc. Roy. Soc. Lond. Ser. A, 295, 13-28. (1966). [5] Oppenheim, K. and Soloukhin, R. I., Experiments in gasdynamics of explosion, Ann. Rev. Fluid Mech. 5, 31-58 (1973). [6] Lee, J.H.S. and Moen, I., The mechanism of transition from de?agration to detonation in vapor cloud explosion, Prog. Energy Combust. Sci., 6, 359-389 (1980). [7] Shepherd, J. E., and Lee, J. H. S., On the transition from de?agration to detonation, in Major Research Topics in Combustion, (M. Y. Hussaini, A. Kumar, and R. G. Voigt, eds), Springer-Verlag, New York, 439-487, 1992. [8] Buckmaster, J., Clavin, P., Linan, A.,Matalon, M., Peters, N., Sivashinsky, G., Williams, F. A., Some developments in Combustion Theory, Proc. Combust. Inst., 30, 1-14, (2005). [9] Sivashinsky, G. I., Some developments in premixed combustion modeling, Proc. Combust. Inst., 29, 1737 (2002). [10] Kuznetsov, M., Alekseev, V., Matsukov, I., and Dorofeev, S., DDT in a smooth tube ?lled with a hydrogen-oxygen mixture, Shock Waves, 14, 205-215 (2005). [11] Brailovsky, I., and Sivashinsky, G. I., Hydraulic Resistance as a Mechanism for De?agration-to-Detonation Transition, Combust. Flame, 122, 492-499 (2000). [12] Kagan, L., and Sivashinsky, G., The transition from de?agration to detonation in thin channels, Combust. Flame, 134, 389-397 (2003).

13

E. Lindblad, D.M. Valiev, B. M¨ ller, J. Rantakokko, P. L¨tstedt and M.A. Liberman u o

[13] Liberman, M. A., Sivashinsky, G. I., Valiev, D. M. and Eriksson, L-E., Hydrodynamic Instability as a Mechanism for De?agration-to-Detonation Transition, 20th International Colloquium on the Dynamics of Explosions and Reactive Systems (ICDERS), Montreal, 2005. [14] Liberman, M. A., Sivashinsky, G. I., Valiev, D. M. and Eriksson, L-E., Numerical Simulation of De?agration-to-Detonation Transition: The Role of Hydrodynamic Flame Instability, Int. J. Transport Phenomena (2006) to appear. [15] Sokolik, A. S. Self-Ignition, Flame and Detonation in Gases, Israel Program for Scienti?c Translations, NASA TTF-125OTS-63-1179, Jerusalem, 1963. [16] Liberman, M. A., Golberg, S. M., Bychkov, V. V., and Eriksson, L.-E., Numerical studies of hydrodynamically unstable ?ame propagation in 2D channels, Combust. Sci. Tech., 136, 221-242 (1998). [17] Bychkov, V. V. and Liberman, M. A., Dynamics and stability of premixed ?ames, Physics Reports, 325, 115-237 (2000). [18] Travnikov, O.Yu., Bychkov, V.V. and Liberman, M.A., In?uence of compressibility on propagation of curved ?ames, Physics of Fluids, 11, 2657-2666 (1999). [19] Kazakov, K. A. and Liberman, M. A., Nonlinear equation for curved stationary ?ames, Physics of Fluids, 14, 1166-1180 (2002). [20] Liberman, M. A., Ivanov, M. F., Peil, O. E., Valiev, D. M., and Eriksson, L.-E., Numerical modeling of ?ame propagation in wide tube, Combust. Theory Modeling, 7, 653-676 (2003). [21] Kadowaki, S. and Hasegawa, T., Numerical simulation of dynamics of premixed ?ames: ?ame instability and vortex-?ame dynamics, Progr. Energy Combust. Sci., 31, 193-241 (2005). [22] Liberman, M. A., Sivashinsky, G. I., Valiev, D. M., Numerical Simulation of De?agration-to-Detonation Transition: The Role of Hydrodynamic Flame Instability, ECCOMAS Thematic Conference on Computational Combustion, Lisbon, Portugal, (2005). [23] Liberman, M. A., Sivashinsky, G. I., Valiev, D. M. and Eriksson, L-E., The DarrieusLandau Instability as a Mechanism for De?agration-to-Detonation Transition Combust. Theory and Modeling (2006) to appear. [24] Oran, E. S., and Gamezo, V. N., Flame acceleration and detonation transition in narrow tubes, Proc. 20th ICDERS, Montreal, Canada (CD) (2005). 14

E. Lindblad, D.M. Valiev, B. M¨ ller, J. Rantakokko, P. L¨tstedt and M.A. Liberman u o

[25] Kagan L. LibermanM. Sivashinsky G, Detonation initiation by hot corrugated walls, ECCOMAS Thematic Conference on Computational Combustion, Lisbon, Portugal, 21-(2005). [26] Eriksson, L.-E., Development and validation of highly modular ?ow solver versions in G2DFLOW and G3DFLOW series for compressible viscous reacting ?ow, Tech. Report, 9970-1162, Volvo Aero Corporation, Trollh¨ttan, Sweden, 1995. a [27] Andersson, N., Eriksson, L.-E. and Davidson, L., Large-Eddy Simulation of subsonic turbulent jets and their radiated sound. AIAA Journal, 43(9), 1899–1912, (2005). [28] Lindblad, E., High order semi-implicit Runge-Kutta methods for separable sti? problems, Master’s Thesis, Uppsala University, Sweden, Dec. 2004. [29] Lindblad, E. and M¨ ller, B., High order implicit-explicit Runge-Kutta methods for u separable sti? problems. Submitted to BIT, Aug. 2005. [30] Kuo, K.K., Principles of Combustion. John Wiley & Sons, New York, 1986. [31] Kee, R.J., Rupley, F.M. and Miller, J.A., CHEMKIN-II: A Fortran chemical kinetics package for the analysis of gas-phase chemical kinetics. SANDIA Report, SAND898009B, 1991. [32] Bussing, T.R.A. and Murman, E.M., Finite-volume methods for the calculation of compressible chemically reactive ?ows. AIAA J., 26 (9), 1070–1078, (1988). [33] Eriksson, L.-E., Time-marching Euler and Navier-Stokes solution for nozzle ?ows with ?nite rate chemistry. Tech. Report, 9370-794, Volvo Flygmotor AB, Trollh¨ttan, a Sweden, 1993. [34] Ascher, U.M., Ruuth, S.J. and Spiteri, R.J., Implict-explicit Runge-Kutta methods for time-dependent partial di?erential equations, Appl. Numer. Math., 25, 151–167, (1997). [35] Pareschi, L. and Russo, G., Implicit-explicit Runge-Kutta schemes and applications to hyperbolic systems with relaxation, J. Scienti?c Comp., 25 (1/2), 129–155, (2005). [36] Yoh, J.J. and Zhong, X., New hybrid Runge-Kutta methods for unsteady reactive ?ow simulation, AIAA J., 42 (8), 1593–1600, (2004). [37] Yoh, J.J. and Zhong, X., New hybrid Runge-Kutta methods for unsteady reactive ?ow simulation: Applications, AIAA J., 42 (8), 1601–1611, (2004). [38] Kennedy, C.A. and Carpenter, M.H., Additive Runge-Kutta schemes for convectiondi?usion-reaction equations, Appl. Numer. Math., 44, 139–181, (2003). 15

E. Lindblad, D.M. Valiev, B. M¨ ller, J. Rantakokko, P. L¨tstedt and M.A. Liberman u o

[39] Rantakokko, J., Partitioning Strategies for Structured Multiblock Grids, Parallel Computing, 26 (12), 1661-1680, (2000). [40] Myklebust, E., Parallelization of a combustion code, Master’s Thesis, Uppsala University, Sweden, July 2004. [41] Zel’dovich, Ya.B., Librovich, V. B., Makhviladze, G. M., and Sivashinsky, G. I., On the development of detonation in non-uniformly preheated gas Astronautica Acta, 15, 313-321, (1970). [42] Zel’dovich, Ya.B., Regime classi?cation of an exothermic reaction with nonuniform initial conditions, Combust. Flame, 39, 211-226, (1980). [43] Gary, J., On certain ?nite di?erence schemes for hyperbolic systems, Math. Comp., 18, 1–18, (1964). [44] Coleman, T.F., Garbow, B.S. and Mor?, J.J., Algorithm 618: FORTRAN subroue tines for estimating sparse Jacobian Matrices, ACM Transaction on Mathematical Software, 10 (3), 346-347, (1984). [45] Knoll, D.A. and Keyes, D.E., Jacobian-free Newton-Krylov methods: a survey of approaches and applications, J. Comput. Phys., 193, 357-397, (2004). [46] Banas, K., A Newton-Krylov solver with multiplicative Schwarz preconditioning for ?nite element compressible ?ow simulations, Commun. Numer. Meth. Engng., 18, 269-275, (2002). [47] Pernice, M. and Tocci, M.D., A multigrid-preconditioned Newton-Krylov method for the incompressible Navier-Stokes equations, SIAM J. Sci. Comput., 23 (2), 398-418, (2001). [48] Frayss?, V., Giraud, L., Gratton, S. and Langou, J., A set of GMRES roue tines for real and complex arithmetics on high performance computers, CERFACS Technical Report TR/PA/03/0, 2003, public domain software available on www.cerfacs/fr/algor/Softs. [49] Saad, Y., SPARSKIT: a basic tool kit for sparse matrix computations (version 2), Report, 1994, public domain software available on www-users.cs.umn.edu/ saad/software/SPARSKIT/sparskit.html. [50] Saad, Y., Iterative Methods for Sparse Linear Systems, 2nd ed., SIAM, Philadelphia, 2003.

16

E. Lindblad, D.M. Valiev, B. M¨ ller, J. Rantakokko, P. L¨tstedt and M.A. Liberman u o

Figure 2: Time sequence of images for the ?ame/shock dynamics near the transition point. Stronger shading corresponds to higher pressure gradient. The time and distance are referred to ? ? 0 and ? respectively. ? = ? (? ? ? ? 0 ) is the ?ame width, where ? 0 is the incipient velocity of the planar ?ame, ? and ? the density and dynamic viscosity, respectively, of the unburnt fuel. The incipient velocity of the planar ?ame ? 0 corresponds to the Mach number ? 0 = ? 0 ? = 0 05, where ? is the speed of sound in the unburnt fuel. Reaction order ? = 2, tube width = 70? , dimensionless (?? ? ) = 8, expansion ratio Θ = ? ?? = 10. ?? and ? are the temperatures activation energy ? = of the unburnt and burnt gases, respectively.

17

E. Lindblad, D.M. Valiev, B. M¨ ller, J. Rantakokko, P. L¨tstedt and M.A. Liberman u o

Figure 3: Time sequence of zoomed temperature contours for the ?ame dynamics near the transition point for the conditions of Fig. 2. Lighter shading corresponds to higher temperature.

18

E. Lindblad, D.M. Valiev, B. M¨ ller, J. Rantakokko, P. L¨tstedt and M.A. Liberman u o

19
Figure 4: Temperature, density, ?ow velocity and pressure pro?les for the fold dynamics of Fig. 3.

E. Lindblad, D.M. Valiev, B. M¨ ller, J. Rantakokko, P. L¨tstedt and M.A. Liberman u o

Figure 5: Temperature contours with explicit chemistry treatment which yields an earlier DDT, but cannot handle the detonation correctly and becomes unstable shortly after the last state in the series above.

Figure 6: Temperature contours with implicit chemistry treatment. When using the implicit Euler method for chemistry, the DDT evolves with retained stability.

Figure 7: Density contours at ? = 1000??? , ? = 5000??? and ? = 10000??? . The reference explicit calculation is shown on the left and an implicit-explicit solution on the right. Inviscid terms are treated explicitly while viscous terms and chemistry are treated implicitly.

20


赞助商链接
相关文章:
更多相关文章: