当前位置:首页 >> >>

Semi-implicit Runge-Kutta schemes for the non-autonomous differential equations in reactive


Semi-Implicit Runge-Kutta Schemes for Non-Autonomous Di erential Equations in Reactive Flow Computations
Jian Wei Shen and Xiaolin Zhong y University of California, Los Angeles, California 90095

Abstract
This paper is concerned with time-stepping numerical methods for computing sti semi-discrete systems of ordinary di erential equations for transient hypersonic ows with thermo-chemical nonequilibrium. The sti ness of the equations is mainly caused by the viscous ux terms across the boundary layers and by the source terms modeling nite-rate thermo-chemical processes. Implicit methods are needed to treat the sti terms while more e cient explicit methods can still be used for the nonsti terms in the equations. For additively split autonomous di erential equations in the form of u0 = f(u) + g(u), three di erent semiimplicit Runge-Kutta methods have been derived and tested in previous papers, where f is treated by explicit Runge-Kutta methods and g is simultaneously treated by three implicit Runge-Kutta methods. The coe cients of up to third-order accuracy have been derived such that the methods are both high-order accurate and strongly A-stable for the implicit terms. However, these semi-implicit Runge-Kutta methods for the autonomous systems cannot be extended to nonautonomous systems of u0 = f(t; u) + g(t; u) because of the coupling between the f and g terms in the split Runge-Kutta methods. In this paper, we derive and test three di erent semi-implicit Runge-Kutta schemes of up to third-order accuracy for the non-autonomous di erential equations using the A-stability and accuracy conditions with four stages. The new schemes have been tested in computations of unsteady reactive ows with explicit time-dependent terms.

motivated by our research in the direct numerical simulation of on the stability and transition of hypersonic boundary layers involving shock interactions and real gas e ects 1, 2]. In addition to the e ects of viscosity, heat-conduction, and di usion, hypersonic ows often contain nonequilibrium processes of thermal excitations and chemical reactions because of high gas temperature and high speeds. One of the major di culties in computing such ows is the sti ness of the governing equations in temporal integrations. The sti ness is mainly caused by the viscous stress and heat ux terms in the boundary layers and by the source terms modeling nite-rate thermo-chemical processes. The viscous terms across the boundary layer are sti because ne-grid spacing is used in the direction normal to the wall. Finite di erence approximation to the viscous equations with these small-size grids lead to sti systems of ordinary di erential equations. The source terms are sti because the chemical and thermal nonequilibrium processes have a wide range of time scales, some of which are much smaller than the transient ow ones. As a result, if explicit methods are used to integrate the sti governing equations, the computations will become very ine cient because the time-step sizes dictated by the stability requirements are much smaller than those required by the accuracy considerations. In order to remove the stability restriction on the explicit methods, implicit methods need to be used. For computing multi-dimensional reactive ow, global implicit methods are seldom used because it takes a prohibitively large amount of computer time and large memory to convert full implicit equations. Practical implicit methods for multi-dimensional reactive ow calculations include the fractional step method (or time-splitting method) and the additive semi-implicit method. This paper is concerned with the additive semiimplicit methods, which additively split the ordinary di erential equations into sti and nonsti terms. The sti terms are treated implicitly while the nonsti terms are treated explicitly. The semi-implicit methods are more e cient than the full implicit methods 1

Introduction
This paper is concerned with numerical methods for computing sti equations for transient hypersonic ows with thermo-chemical nonequilibrium. This work is
Visiting Scholar, Member AIAA Professor, Mechanical and Aerospace Engineering Department, Member AIAA
y Assistant

Copyright c 1996 by American Institute of Aeronautics and Astronautics, Inc. All rights reserved.

for reactive ow computations because the sti terms can be easily separated from the rest of the equations. The standard semi-implicit method for direct numerical simulation of incompressible turbulence is to use the implicit Crank-Nicolson method for the viscous terms normal to the wall and the explicit AdamsBashford method for the rest of terms 3, 4, 5, 6, 7] (ABCN method). For compressible reactive ow, a semi-implicit MacCormack method 8, 9, 10, 11] has been used to compute the chemical source terms implicitly while the uid terms are computed explicitly. The temporal accuracy of these two methods, however, is usually second-order accurate at most. To obtain simultaneously high-order accuracy and good stability properties, the simultaneous coupling between the explicit and implicit terms need to be considered. This can be accomplished by the additive semi-implicit Runge-Kutta methods. The derivation of an additive semi-implicit method with both high accuracy and good stability is di cult because of the coupling between the explicit and implicit terms. The rst additive Runge-Kutta methods for sti ordinary di erential equations were studied by Cooper and Sayfy 12, 13]. They derived additive Runge-Kutta methods to solve a system of di erential equations in a form of x0 = J(t)x+g(t; x), where the linear term on the right hand side of the equation was sti . In previous papers, Zhong 14, 15] derived and analyzed three di erent sets of additive semi-implicit Runge-Kutta methods for additively split ordinary differential equations in the form of u0 = f(u) + g(u), where the nonsti term f is treated by explicit RungeKutta methods, and the sti term g is simultaneously treated by three implicit Runge-Kutta methods. The three implicit methods for g are a diagonally implicit Runge-Kutta method and two Rosenbrock linearized Runge-Kutta methods 16] with di erent ways of evaluating Jacobian matrices. The semi-implicit RungeKutta methods are derived and analyzed based on the general form of implicit Runge-Kutta formulas. The fully implicitand Rosenbrock semi-implicit RungeKutta methods are both high-order and strongly Astable for the implicit terms. The strongly A-stable methods are needed for numerical results to reach correct asymptotic values for very sti problems. These new schemes have been tested in model equations and in reactive hypersonic ow computations 15]. These additive semi-implicit Runge-Kutta methods, however, cannot be applied to non-autonomous systems relating to unsteady reactive ows with explicit time terms in the equation of u0 = f(t; u) + g(t; u) due to the coupling between the f and g terms. An example of such ow is the unsteady reactive ow with 2

time-dependent perturbations or boundary conditions. The objective of this paper is to extend the previous work of semi-implicit Runge-Kutta methods for the autonomous systems to the non-autonomous systems, where the nonsti f(t; u) term is treated by explicit Runge-Kutta methods, and the sti term g(t; u) is simultaneously treated by three implicit Runge-Kutta methods. The three implicit methods for g(t; u) are a diagonally implicit Runge-Kutta method and two Rosenbrock linearized Runge-Kutta methods 16] with di erent ways of evaluating Jacobian matrices. The coe cients are derived such that the new methods are simultaneously high-order accurate and A-stable for the implicit term. We found that the semi-implicit Runge Kutta method for non-autonomous di erential equations requires four stages to obtain third-order accuracy due to the strong coupling among the coe cients of 18 accuracy and stability conditions. The new schemes are tested in the computations of unsteady reactive ows with explicit time-dependent terms.

Semi-Implicit RK Methods
General Formulas of Semi-Implicit RK Methods for Non-Autonomous Systems
A general partial di erential equation for reactive ows can be written as @ U + @ Fj = W @t @xj (1)

where U is the conservation variables, Fj is the ux terms, and W is the source term due to nite-rate reactions. In the semi-discretization approach, the spatial derivatives in the governing partial di erential equations are rst approximated by spatial discretization methods. The spatial discretization leads to a system of rst-order ordinary di erential equations. For systems of unsteady ow with time-dependent forcing terms or boundary conditions, the semi-implicit di erential equations are non-autonomous, i.e., du = f (t; u) + g(t; u) dt (2)

where u is the vector of discretized ow eld variables. The right hand side of the di erential equation above is additively split into two terms, g and f , where g is the vector resulting from the spatial discretization of the sti terms, and f is the vector resulting from the spatial discretization of the rest of the nonsti ow equations. In general, the splitting of f and g terms is not unique.

The multiple time scales in the ow equations and the di culty in computing fully implicit multidimensional ow equations for reactive ow calculations require that the sti term g to be computed implicitly and the non-sti term f to be computed explicitly. Because of the coupling between f and g, however, a straight-forward combination of an explicit and an implicit time-stepping method will not maintain the original accuracy of the individual methods. Therefore, the derivation of high-order semi-implicit methods needs to simultaneously consider the coupling between the implicit and explicit terms in both stability and accuracy analysis. The stability conditions in the semi-implicit method should be limited by the explicit terms only, i.e., the coupled method should be A-stable for the implicit terms when the explicit terms satis es the CFL condition (or similar conditions). Therefore, we use semi-implicit Runge-Kutta methods to achieve simultaneous high-order accuracy and good stability. The Runge-Kutta methods are one-step methods involving intermediate stages to achieve high-order accuracy 17, 18]. A general r-stage additive semiimplicit Runge-Kutta method integrates Eq. (2) by simultaneously treating f explicitly and g implicitly:

3 2 i?1 X 4I ? hai J(tn + qih; un + dij kj )5 ki
j =1

= hf (tn + ri h; un + +hg(tn + si h; un +

i?1 X j =1 i?1 X j =1

bij kj ) cij kj ) ; r) (6)

(i = 1;

@g Where J = @ u is the Jacobian matrix of the sti term g. The Rosenbrock additive semi-implicit RungeKutta method given by Eqs. (5) and (6) is similar to the implicit methods used in computational uid dynamics and is much more e cient than the diagonally implicit version given by Eqs. (3) and (4). But, for some strongly nonlinear problems, the nonlinear diagonally semi-implicit method given by Eqs. (3) and (4) is necessary because it is more stable than the Rosenbrock additive semi-implicit Runge-Kutta method for nonlinear problems.

un+1 = un +

r X j =1

!j kj
i?1 X j =1 i?1 X n+ cij j j =1

(3) bij kj )

Therefore, three versions of the additive semiimplicit Runge-Kutta methods are derived to be both high-order accurate and strongly A-stable for the implicit terms, i.e., Method A: \Fully implicit" additive semi-implicit Runge-Kutta method given by Eqs. (3) and (4). Method B: Rosenbrock additive semi-implicit RungeKutta method given by Eqs. (5) and (6), and dij = qi = 0. Method C: Rosenbrock additive semi-implicit RungeKutta method given by Eqs. (5) and (6), and qi = si , and dij = cij . The rth-stage additive semi-implicit Runge-Kutta methods are termed ASIRK-rA methods, ASIRK-rB methods, and ASIRK-rC methods for Methods A, B, and C respectively.

ki = hff (tn + rih; un +
+ g(tn + si h; u (i = 1;

k + aiki )g (4)

; r)

where h is the time-step size, and ai , bij , cij , ri, si , and wj are parameters to be determined by accuracy and stability requirements. Because g is treated by a diagonally implicit Runge-Kutta method, Eq. (4) is a nonlinear equation at every stage of the implicit calculations if g is a nonlinear function of u. The computations of this method are relatively ine cient, since nonlinear solvers are required to solve such nonlinear equations. A more computationally e cient additive semiimplicit Runge-Kutta method is a semi-implicit extension of the Rosenbrock Runge-Kutta method 16],

Linear Stability Conditions
The parameters of the additive semi-implicit RungeKutta methods are chosen based on both stability and accuracy requirements with simultaneous coupling between the explicit and implicit terms. The use of an implicit method for the sti term g permits a larger time step than that allowed by a fully explicit method. 3

un+1 = un +

r X j =1

!j kj

(5)

Unlike the explicit Runge-Kutta methods, whose stability conditions are the same for di erent choices of parameters as long as they have the same stages and accuracy, the stability properties of the additive semiimplicit Runge-Kutta methods of the same stages are di erent for di erent choices of parameters because of the coupling between the f and g terms. For simplicity, only a linear stability analysis is conducted in this paper for a special kind of test functions. The stability condition for an additive semi-implicit time-stepping scheme is analyzed by considering a simpli ed linear model equation: du = u + u g dt f (7)

In order to obtain a correct asymptotic decay for sti terms, it is desirable to have a strong A-stability (Lstability) condition for the semi-implicit schemes, i.e.,
hj

lim j fh f ; h g gj = 0 j!1
g

(11)

The strong A-stability for the implicit term assures that the numerical solutions approach the correct solutions as step sizes increase. For the three additive semi-implicit Runge-Kutta methods, the strongly Astable condition can be obtained from Eqs. (8), (9), and (11) as follows: 1+ where
i r X j =1

!j j = 0

(12)

@f where f and g represent the eigenvalues of @ u and @g @ u in Eq. (2). They are complex parameters satisfying Ref f g 0 and Ref g g 0 respectively. In general, j g j is much larger than j f j for sti equations. Though Eq. (2) cannot be reduced to this model equation if the Jacobians of f and g do not commute, Eq. (7) is used as the rst step in analyzing the linear stability properties of the additive semi-implicit RungeKutta methods. Further studies are needed to analyze the general stability properties of the additive RungeKutta methods using the nonlinear stability analysis by Hairer, Bader, and Lubich 19].

= ? 1+

i?1 X j =1

cij j ]=ai

(i = 1;

; r) (13)

Substituting Eq. (7) into any of the three additive semi-implicit Runge-Kutta methods leads to the same equation for the characteristic root as follows:
r n+1 X (8) = uun = 1 + !j kj j =1 Pi P h f (1 + j?1 bij kj ) + h g (1 + ij?1 cij kj ) =1 =1 ki = 1 ? ai h g (i = 1; ; r) (9)

For practical reactive ow problems, it is important that the intermediate variables at each stage of the Runge-Kutta computations maintain their physical meanings, i.e., it is not acceptable to have negative temperatures in an intermediate stage even if the nal results are positive. Therefore, we impose the following additional condition on the additive semi-implicit Runge-Kutta methods: ai > 0 (14)

Accuracy Conditions
Additive semi-implicit Runge-Kutta schemes are derived to be high-order accurate with the simultaneous coupling between the explicit and implicit terms. Taylor series expansions lead to the following accuracy conditions if four steps are kept as a general case:

where = fh f ; h g g is a function of h

f

and h g .

An A( ) stability region of a semi-implicit method in the complex plane of h f is de ned as the region where

1st order:
!1 + !2 + !3 + !4 = 1 (15)

j fh f ; h g gj 1

(10)

for h f within the region and for all h g within a wedge bounded by ? ; + ] in the complex plane. When = =2, the semi-implicit method is A-stable for h g . 4

2nd order:
!2 r2 + !3r3 + !4 r4 = 1 2 (16)

(17) !1s1 + !2 s2 + !3 s3 + !4s4 = 1 2 1 (18) !2b21 + !3 (b31 + b32) + !4 (b41 + b42 + b43) = 2 !1a1 + !2(a2 + c21) + !3(a3 + c31 + c32 ) 1 +!4 (a4 + c41 + c42 + c43) = 2 (19)

!1 a1s1 + !2 s2 (c21 + a2) + !3 s3 (c31 + c32 + a3) (36) 1 +!4 s4 (c41 + c42 + c43 + a4) = 3 For up to second-order accuracy, direct combination of an explicit and an implicit Runge-Kutta methods will result in an additive semi-implicit Runge-Kutta method with the same order of accuracy as long as the two schemes have the same set of !i . However, for accuracy equal to or higher than third order, the direct combination of an explicit and an implicit methods will likely be only second-order accurate because of the coupling between the explicit and implicit terms. We search for the optimal parameters in the additive semi-implicit Runge-Kutta schemes by simultaneously imposing the stability and accuracy conditions discussed above. To be consistent with the explicit Runge-Kutta methods, the coe cients for the t terms are set to be ri =
i?1 X j =1

3rd order:
2 2 2 !2r2 + !3r3 + !4 r4 = 1 3 1 !1s2 + !2 s2 + !3 s2 + !4s2 = 3 1 2 3 4 !1a1 s1 + !2(c21 s1 + a2 s2 ) + !3(c31 s1 + c32 s2 1 +a3 s3 ) + !4(c41s1 + c42s2 + c43s3 + a3 s4 ) = 6 !1a2 + !2 c21a1 + a2 (c21 + a2)] + !3 a1c31 1 +c32(a2 + c21) + a3(c31 + c32 + a3 )] +!4 c41a1 + c42(a2 + c21) + c43(c31 + c32 + a3 ) 1 +a4 (c41 + c42 + c43 + a4)] = 6 1 !3b32r2 + !4 (b42r2 + b43r3) = 6 1 !2b2 + !3 (b31 + b32)2 + !4(b41 + b42 + b43)2 = 3 21 !2r2 b21 + !3 r3(b31 + b32) + !4r4 (b41 1 +b42 + b43) = 3 1 !3b32b21 + !4 b42b21 + b43(b31 + b32)] = 6 !2(b21 a2 + b21a1) + !3 (a1b31 + a2 b32 + c21 b32 +b21c32 + a3 b31 + a3 b32) + !4 b41a1 +b42(a2 + c21) + c42b21 + b43(c31 + c32 + a3 ) 1 +c41(b31 + b32) + a4 (b41 + b42 + b43)] = 3 !2a2 r2 + !3(c32 r2 + a3r3 ) + !4(c42r2 1 +c43r3 + a4 r4) = 6 !2b21s1 + !3 (b31s1 + b32s2 ) + !4(b41s1 1 +b42s2 + b43s3 ) = 6

(20) (21) (22) (23) (24) (25) (26) (27) (28) (29) (30)

bij

(37)

for all the three methods. But for the implicit t terms, method A is di erent from method B and C. For method A, si = ai +
i?1 X j =1

cij

(38)

For method B and method C,

Method A :

!1a2 + !2(c21 + a2 )2 + !3 (c31 + c32 + a3)2 1 1 +!4 (c41 + c42 + c43 + a4 )2 = 3 (31) !1a1 s1 + !2s2 (c21 + a2 ) + !3s3 (c31 + c32 + a3 ) (32) 1 +!4 s4 (c41 + c42 + c43 + a4) = 3

si = ri

(39)

The reason for this is obviously that if equation (38) were available for method C, then combine Equations (21) and (35) we have !1 a2 + !2 a2 + !3a2 + !4a2 = 0 1 2 3 4 (40)

1 !2c2 + !3 (c31 + c32)2 + !4 (c41 + c42 + c43)2 = 3 (33) 21 !2s2 c21 + !3 s3 (c31 + c32) + !4s4 (c41 1 +c42 + c43) = 3 (34)

Method B :

That means there are no solutions that satisfy the conditions of !i > 0 and ai > 0. For method B, s1 =0, but a1 is required to be a1 > 0 for the sti terms. The di erence between the set of accuracy equations derived for conventional explicit Runge-Kutta scheme and for the explicit part of the semi-implicit scheme di ers only in the additional equations generated by 5

!2(c2 + 2a2 c21) + !3 (c31 + c32)2 + 2a3 (c31 + c32)](35) 21 1 +!4 (c41 + c42 + c43)2 + 2a4(c41 + c42 + c43)] = 3

Method C :

the third-order cross terms as a result of the coupling between the explicit and implicit terms. And similarly the accuracy equations for implicit part of the semiimplicit method is same as the conventional implicit method for rst- and second-order schemes only, excluding third-order schemes. Since it is di cult to obtain close form solution of the accuracy and stability equations for a coupled semiimplicit scheme, a numerical searching method is used to nd the optimal parameters. A computer program were developed to conduct extensive search for the optimal parameters numerically. The computational search makes it possible to locate the optimal parameters to satisfy both the stability and accuracy conditions for the methods.

ASIRK-2B Method: 8 > I ? ha1J(tn; un)] k1 = hff (tn; un) + g(tn; un)g > < I ? ha2J(tn; un)] k2 = hff (tn + r2h; un + b21k1) (45) n + > : un+1 = un + !1k1 + !2g(tn + s2 h; u + c21k1 )g k2 ASIRK-2C Method: 8 > I ? ha1J(tn; un)] k1 = hff (tn; un) + g(tn; un)g < I ? ha2J(tn + s2 h; un + c21 k1)] k2 = hff (tn+ > r2 h; un + b21k1 ) + g(tn + s2 h; un + c21k1 )g (46) : un+1 = un + !1k1 + !2k2 Similar to the rst-order case, the coe cients for the second-order methods are the same as the corresponding coe cients for a autonomous systems, and the coe cients are the same for all three methods. Ref. 15] shows two sets of coe cients: Case I:
1 !1 = 2 1 a1 = 4

First-Order Additive Semi-Implicit RungeKutta Methods
The expressions for the rst-order methods are ASIRK-1A Method:

k1 = hff (tn; un) + g(tn + s1 h; un + a1 k1 )g un+1 = un + !1 k1
ASIRK-1B and ASIRK-1C Methods:

(41)

I ? ha1 J(tn; un)] k1 = hff (tn; un) + g(tn; un)g (42) un+1 = un + !1 k1 Case II:
The coe cients are the same as the autonomous systems derived in Ref. 15]: !1 = a1 = s1 = 1 (43)

!2 = 1 2 1 a2 = 3

b21 = 1 5 c21 = 12

!1 = 1 p 2 a1 = 1 ? 22

!2 = 1 p 2 a2 = 1 ? 22

b21 = 1 p c21 = 2 ? 1

The coe cients are the same for Methods A, B, C for the rst-order case. The stability condition for the rst-order additive semi-implicit Runge-Kutta methods is the same as rst-order explicit Runge-Kutta methods for h f and is strongly A-stable for h g .

where ri and si are given by Eqs. (37) and (38) or (39) . The methods are second-order accurate and strongly A-stable for the implicit term h g .

ASIRK-3A Method 8 The expressions for the second-order methods are: > k1 = hff (tn; un) + g(tn + s1 h; un + a1k1 )g > > k2 = hff (tn + r2h; un + b21k1 ) > ASIRK-2A Method: > g(t + s n + c21k1 > < k3 = hff (tn ++3h; n n +2bh; u1 + b32k2)+ a2 k2 )g r u 31k 8 +g(tn + s3 h; un + c31k1 + c32 k2 + a3 k3 )g (47) > > k1 = hff (tn; un) + g(tn + s1 h; un + a1 k1 )g >k < k2 = hff (tn + r2 h; un + b21k1 ) > f (t h; n + + + (44) > +4g= h+f s n + r4+ cu k b41ck1k b42ck2k b43ak3k) )g n + c k + a k )g > (tn 4 h; un 41 1 + 42 2 + 43 3 + 4 4 +g(tn + s2 h; u > 21 1 2 2 > > n+1 n : un+1 = un + !1k1 + !2k2 + !3k3 + !4k4 : u = u + !1 k1 + !2 k2 6

Second-Order Additive Semi-Implicit RungeKutta Methods

Third-Order Additive Semi-Implicit RungeKutta Methods

ASIRK-3B Method

8 > I ? ha1 J(tn; un)] k1 = hff (tn; un) + g(tn ; un)g > > I ? ha2 J(tn; un)] k2 = hff (tn + r2h; un + b21k1 ) > > + (t + s n + > > I ? ha3 J(tn; un)] k3g= n ff (t2nh; ur3h; c21 k1)bg k1 h + un + 31 < n + c k + c k )g ) g(tn s 31 32 2 > I ? +b32k2n;+ n)] k4+= 3hh;fu n + r4h;1un + b41k1 ha4 J(t u f (t > > > +b42k2 + b43k3) + g(tn + s4 h; un + c41k1 > > k 42 > n+1 n : u = u + !1 k1 + !2 k2 +c!3k2 + c43k43 )g + k3 + !4
ASIRK-3C Method

The searching procedure for the numerical solutions is as follows: First, !1 , !2 and !4 are chosen as free parameters to obtain !3 from Equation (15). Then, b32, b42 and b43 are chosen as free parameters to determine r2 , r3 and r4 from Equations (16),(20) and (24). For method B and method C, s1 -s3 are determined by si = ri ; for method A, another parameter s4 is chosen to solve Equations (17),(21) and (30) for s1 , s2 and s3 . Next, Equations (22),(23),(28) and (29) are solved for (48) a1, a2 , a3 and a4 with c32, c42 and c43 as free parameters. Finally, the solutions are checked by stability Equation (10) and (14). If the solutions do not satisfy those conditions, the parameters are changed to repeat the procedure until all ai > 0. The searching range for free parameters of !i and s4 is from 0 to 1. The range for bij and cij is from -2 to 2. The parameters for the explicit and implicit terms are di erent for the three versions of ASIRK3 methods because the explicit/implicit coupling. ASIRK-3A Method:

8 > > > > > > > > < > > > > > > > > :

I ? ha1 J(tn; un)] k1 = hff (tn; un) + g(tn; un)g I ? ha2 J(tn + s2 h; un + c21k1)] k2 = hff (tn + r2h; !1 = 0:13 un + b21k1) + g(tn + s2 h; un + c21k1 )g I ? ha3 J(tn + s3 h; un + c31k1 + c32k2 )] k3 = !3 = 0:52 hff (tn + r3h; un + b31k1 + b32k2 ) (49) b21 = 0:338170 +g(tn + s3 h; un + c31 k1 + c32k2 )g b32 = 0:779584 I ? ha4 J(tn + s4 h; un + c41k1 + c42k2 + c42k3)] k4 b42 = 0:200000 a1 = 1:174810 = hff (tn + r4h; un + b41k1 + b42k2 + b43k3 ) a3 = 0:158717 + g(tn + s4 h; un + c41k1 + c42k2 + c43k3 )g n+1 = un + ! k + ! k + ! k + ! k u a21 = ?0:293999 1 1 2 2 3 3 4 4
c32 = 0:200000 c42 = 1:780818

!2 = 1 41 !4 = 10 b31 = ?0:019084 b41 = ?0:300000 b43 = 0:300000 a2 = 0:526766 a4 = 0:100000 c31 = 0:149135 c41 = ?1:130818 c43 = ?0:500000

We use four stages in ASIRK-3 to get third-order accuracy. In fact, if three-stage ASIRK schemes were used, there would be 18 undetermined parameters for the third-order method. After satisfying the accuracy and strong A-stability conditions, there were only two free parameters for the ASIRK-3 methods. If !1 and !2 were chosen as free parameters, then the third-order accuracy equations and one stability equations would be solved exactly. We searched !1 and !2 from 0 to 1 and found that no solution satis es all ai > 0. For a given set of !1 and !2, the accuracy conditions lead to either a1 or a3 to be zero. This is not an acceptable solution because the resulting ASIRK-3 schemes do not satisfy stability condition (12). That means if a three-stage additive semi-implicit Rungge-Kutta scheme was used to obtain third-order accuracy for non-autonomous differential equations, one stage whould be a full explicit scheme. In order to get the solutions that satisfy the third-order accuracy and stability conditions, meanwhile all ai are positive, four-stage semi-implicit RK schemes are needed to add more parameters to get the third-order accuracy scheme. 7

ASIRK-3B Method: !1 = 1 8 !3 = 0:525 b21 = 0:309921 b32 = 0:591232 b42 = ?0:55000 a1 = 0:130476 a3 = 0:067873 c21 = 0:160000 c32 = 0:400000 c42 = ?0:500000 ASIRK-3C Method: !1 = 1 8 !3 = 0:525 !2 = 1 41 !4 = 10 !2 = 1 41 !4 = 10 b31 = 0:169758 b41 = ?0:37000 b43 = 1:14999 a2 = 0:052913 a4 = 0:424531 c31 = 0:361513 c41 = ?0:974181 c43 = 1:000000

b21 = 0:324692 b32 = 0:767118 b42 = ?1:000000 a1 = 0:170366 a3 = 0:041351 c21 = 0:150000 c32 = 0:706262 c42 = ?1:000000

b31 = ?0:000745 b41 = 0:300000 b43 = 0:890000 a2 = 0:107914 a4 = 0:029692 c31 = 0:033636 c41 = 0:314661 c43 = 0:696340

Table 1: Temporary Accuracy of the ASIRK2A and ASIRK-3A Methods(t = 2:5, and uex = ?0:80114362)
RK-3 e3 h=0.25 1.26D-3 h/2 1.53D-4 h/4 1.88D-5 h/8 2.33D-6 h/16 2.90D-7 h/32 3.62D-8

4t

The parameters with longer signi cant digits can be obtained from the authors. The methods using the coe cients above are third-order accurate and strongly A-stable for the implicit term h g .

R3 8.2 8.1 8.1 8.0 8.0 -

ASIRK-3A e3 R3 1.40D-3 7.1 1.96D-4 7.6 2.58D-5 7.8 3.29D-6 7.9 4.15D-7 8.0 5.20D-8 -

ASIRK-2A e2 R2 1.11D-3 4.2 2.65D-4 4.1 6.50D-5 4.0 1.61D-6 4.0 4.01D-7 4.0 1.00D-7 -

Test Cases
Systems of Ordinary Di erential Equations
In order to check the temporal accuracy of semiimplicit Runge-Kutta method derived in this paper, two cases of systems of ordinary di erential equations are tested. Case I Non-sti Systems of ODE First, we consider a system of ordinary di erential equations as follows du = Au + f (t) dt (50)

Table 1 shows the results of the grid re nement study. The results of 2nd-order and 3rd-order additive semi-implicit Runge-Kutta scheme (ASIRK-2A and ASIRK-3A) using the same time steps are given in this table. As a comparison, the results of conventional third-order explicit Runge-Kutta method(RK3) are also listed in the table. Because this case is a non-sti system of ordinary di erential equations, the right hand part of the equations are treated as the part g in equation (2). In another word, ri and bij in Eq.(47) are set to zero, the equations are calculated with implicit method. Therefore there are some difference between the 3rd-order explicit RK method and semi-implicit RK method. Just as predicted by equation (53), the value of Rp gives out the temporal accuracy of the schemes. Table 1 shows that the ASIRK-2A is second-order accuracy and the ASIRK-3A using four steps is third-order temporal accuracy. Case II Systems of ODE with Sti Terms

Then we tested the scheme of ASIRK-3C for systems where of ordinary di erential equations with sti terms given by 2 3 2 3 0 1 0 0 5(51) A=4 0 0 1 5 f =4 0 du = A(t)u + f (t) + 1 B(t)u + f (t)] (54) ?2 ?5 ?4 ?4 sin(t) ? 2 cos(t) 1 2 dt the initial condition is u(t=0 ) = f1; 0; ?1g. Because this simple case has an exact solution of u1 = cos(t), we can evaluate the temporal order of accuracy of semi-implicit scheme by the re nement of time steps. The errors of numerical solution uh is eh = uex ? uh (52) where 2 3 12:2 5 sin(t) sin(4t) A=4 4 15 2 cos(2t) 5 3 4 cos(3t) 7

2 3 cos(t) f1 = 4 2 sin(t) 5(55) 1

where eh depends on the time step h. For a p-th order numerical method, eh p e h = Rp = 2 2 (53) 8

2 3 2 3 42:2 50:1 ?42:1 sin(2t) sin(t) B = 4 ?66:1 ?58 58:1 cos(2t) 5 f2 = 4 cos(t) 5(56) 26:1 42:1 ?34t0:1 sin(2t)
where is a small parameter chosen to be 0.01. The results of ASIRK-2C and ASIRK-3C are listed in Table 2. The exact solution uex used in equation (52) in

this case is obtained from the Richardson extrapolated solution at the smallest time step. The results show that the ASIRK-3C method is third-order accuracy.

Table 2: Temporary Accuracy of ASIRK-2C and ASIRK-3C Method for Sti Equations(t = 0:3967, and uex = 0:310395394964)
ASIRK-3C 4t e3 R3 h=0.001 0.12 5.6 h/2 2.14D-2 7.5 h/4 2.87D-3 10.6 h/8 2.71D-4 10.7 h/16 2.53D-5 10.0 h/32 2.54D-6 9.2 h/64 2.77D-7 ASIRK-2C e2 R2 0.47 10.4 4.51D-2 3.9 1.14D-2 4.0 2.86D-3 4.0 7.15D-4 4.0 1.78D-4 4.0 4.47D-5 -

distance. Fig.4 shows that the highest pressure region on the wall is near the stagnation line. Figs.5-8 are the species concentrations of N2 and N. Fig.5 is the distributions of mass fractions along the stagnation line and Fig.6 is along the wall. The gure shows that N2 is dissociated behind the bow shock due to high temperature, but it recombines in the boundary layer because of the low wall temperature. Fig. 7 and Fig. 8 are contours of N2 and N mass fractions. The calculations show that the new method is robust and accurate.

Conclusions
The additive semi-implicit Runge-Kutta methods of up to third-order accuracy for non-autonomous di erential equations have been derived and tested in this paper. Unlike the methods for autonomous di erential equations, four steps are required to get third-order temporal accuracy and every step is implicit for the sti terms. The temporal order of accuracy of the four-stage additive semi-implicit Runge-Kutta methods have been validated by the re nement study of time steps. The application of the third-order ASIRK method to two dimensional thermo-chemical reactive ow shows that the semi-implicit methods are e cient and robust for the sti di erential equations.

Two-Dimensional Viscous Reactive Flows
After validating the accuracy of the present ASIRK3 method, we apply these schemes to a two-dimensional viscous ows over a circular cylinder with thermochemical nonequilibrium. There are two features of viscous reactive ows that are sti for time integration, one is the stretched grid in the viscous boundary layer near the wall, and another is the source terms due to the nite rate thermal and chemical processes. The ASIRK-3 methods derived in this paper are applied to these sti equations. The free stream conditions are CN2 = 0:927 u1 = 5590m=s p1 = 2910pa CN = 0:073 T1 = 1833K Tw = 1000K

Acknowledgments
This research was supported by the Air Force O ce of Scienti c Research under grant numbers F4962094-1-0019 and F49620-95-1-0405 monitored by Dr. Len Sakell. More information on this work and related ongoing research can be found on the World Wide Web at http://cfdlab5.seas.ucla.edu.

To resolve the viscous boundary layer on the wall, the grids are stretched in the direction normal to the wall. The total number of grid points used in numerical calculations is 42 142. Fig.1 shows the structured grids used in the calculation. The Navier-Stoke equations are solved by nite volume method. A secondorder TVD scheme is used for the descretization of spatial uxes and the 4-stage semi-implicit Runge-Kutta (ASIRK-3C) method is used to integrate the equations in time. The source terms due to thermo-chemical nonequilibrium, as well as the uxes in the direction normal to the wall, are treated implicitly. Fig. 2 shows the pressure contours of the mixture gas. Fig.3 and Fig.4 are the distributions of pressure along the stagnation line and along the wall, respectively. Fig.3 shows that the pressure in the shock layer increases about two order-of-magnitude as compared with the freestream conditionis. The change of pressure indicates that the thickness of the shock front is about two grids spatial 9

References
1] C.-H. Chiu and X. Zhong. Numerical Simulation of Transient Hypersonic Flow Using the Essentialy Nonoscillatory Schemes. AIAA Journal, Vol. 34, No. 4, pp. 655{661, 1995. (Also AIAA paper 950469, 1995). 2] X. Zhong, X. Joubert, and T. K. Lee. The interaction of freestream disturbances and the bow shock in hypersonic ow past a cylinder. 1995. Presented in 20th International Symp. on Shock Waves, July 23-28, 1995, Pasadena, California.

3] P. Moin and J. Kim. On the numerical solution of the time-dependent viscous incompressible uid ows involving solid boundaries. J. of Computational Physics, 35:381{392, 1980. 4] Thomas A. Zang and M. Y. Hussaini. Numerical experiments on subcritical transition mechanisms. AIAA Paper 81-1227, 1981. 5] R. S. Rogallo and P. Moin. Numerical simulation of turbulent ows. Ann. Rev. Fluid Mech., 16:99{ 137, 1984. 6] C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang. Spectral Methods in Fluid Dynamics. Springer-Verlag, 1988. 7] L. Kleiser and T. A. Zang. Numerical Simulation of Transition in Wall-Bounded Shear Flows. Ann. Rev. Fluid Mech., Vol. 23, pp. 495{535, 1991. 8] T. R. A. Bussing and E. M. Murman. A nite volume method for the calculation of compressible chemically reacting ows. AIAA Paper 85-0296, 1985. 9] J. P. Drummond, R. C. Rogers, and M. Y. Hussaini. Spectral methods for modeling supersonic chemically reacting ow elds. AIAA Paper 850302, 1985. 10] H. C. Yee and J. L. Shinn. Semi-implicit and fully implicit shock capturing methods for hypersonic conservation laws with sti source terms. AIAA Paper 87-1116, 1987. 11] R. J. Leveque and H. C. Yee. A study of numerical methods for hyperbolic conservation laws with sti source terms. J. of Computational Physics, 86:187{210, 1990. 12] G. J. Cooper and A. Sayfy. Additive methods for numerical solutions of ordinary di erential equations. Mathematics of Computation, 35(152):1159{1172, October 1980. 13] G. J. Cooper and A. Sayfy. Additive rungekutta methods for sti ordinary di erential equations. Mathematics of Computation, 40(161):207{ 218, January 1983. 14] X. Zhong. Semi-ImplicitRunge-Kutta Schemes for Direct Numerical Simulation of Transient HighSpeed Reactive Flows. To appear in the Journal of Computational Physics, 1996. 15] X. Zhong. New High-Order Semi-implicit RungeKutta Schemes for Computing Transient Nonequilibrium Hypersonic Flow. AIAA Paper 95-2007, 1995. 10

16] H. H. Rosenbrock. Some general implicit processes for the numerical solution of di erential equations. Comput. J., 5:329, 1963. 17] J. D. Lambert. Computational Methods in Ordinary Di erential Equations. John Wiley & Sons, 1973. 18] J. C. Butcher. Implicit Runge-Kutta processes. Math. Comp., 18:50{64, 1964. 19] E. Hairer, G. Bader, and CH. Lubich. On the stability of semi-implicit methods for ordinary di erential equations. BIT, 22:211{232, 1982.

Pw
1.50E5

1.00E5

5.00E4

-0.040

-0.035

-0.030

X

Figure 1: Computational Grids.

Figure 3: Mixtured pressure on the stagnation line.

Pw
151955 142013 132070 122128 112186 102243 92301.2 82358.9 72416.6 62474.3 52532 42589.7 32647.4 22705.1 12762.8

1.50E5

1.00E5

5.00E4

-0.03

-0.02

-0.01

0.00

0.01

0.02

0.03

S

Figure 2: Pressure contours of the mixtured nitrogen.

Figure 4: Pressure of the mixtured nitrogen on the wall.

11

N2
N2 N
F E D 5 C B A 9 8 7 6 5 B 2 E 4 3 3 2 1 0.92 0.90 0.89 0.88 0.86 0.85 0.84 0.82 0.81 0.79 0.78 0.77 0.75 0.74 0.72

0.90 0.20

F D8 97 C

5

0.85

N2 N

0.15

A 32 41

0.80

0.10

0.75

0.05

0.00 -0.040 -0.035 -0.030

Figure 5: Species concetrations on the stagnation line.

Figure 7: Species contours of N2 .

N2

N
F

N
0.23 0.22 0.20 0.19 0.17 0.16 0.14 0.12 0.11 0.09 0.08 0.06 0.05 0.03 0.02

0.185
A

E D C C B A 9 8 7 6 5 B 4 3 4 2 1

0.790

4 9 2 3 F BE 7 8 5 6 C D ED 1

0.180 0.785

N2 N
0.780

A

0.175

0.170 0.775 -0.03 -0.02 -0.01 0.00 0.01 0.02 0.03

S

Figure 6: Species concentrations on the wall.

Figure 8: Species contours of N.

12


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