- Email: [email protected]
Singular perturbation techniques in the study of a diatomic gas with reactions of dissociation and recombination
Applied Mathematics and Computation 146 (2003) 509–531 www.elsevier.com/locate/amc
Singular perturbation techniques in the study of a diatomic gas with reactions of dissociation and recombination M. Galli a, M. Groppi a, R. Riganti b, G. Spiga a
a,*
Dipartimento di Matematica, Universit a di Parma, Via D’Azeglio 85, 43100 Parma, Italy b Dipartimento di Matematica, Politecnico di Torino, Corso Duca degli Abruzzi 24, 10129 Torino, Italy
Abstract Model equations for the description of chemical reactions of dissociation and recombination through a transition state constitute necessarily singular perturbation problems because of the very short lifetime of the excited molecules, compared to all other characteristic times. Asymptotic solutions are derived and tested versus an accurate numerical solution by resorting to different expansion algorithms, taking initial layer corrections into account. Since the singular problem turns out to be also singularly perturbed, techniques like Hilbert and Chapman–Enskog expansions, typical of the Boltzmann equation for gas kinetic theory, have been considered. The simple stationary state approximation, very popular in chemistry, is shown to provide very good results out of the initial layer if equipped with the proper initial conditions, as obtained from the present analysis. Ó 2002 Elsevier Inc. All rights reserved. Keywords: Dissociation/recombination reactions; Singular perturbation problems; Asymptotic expansions; Initial layer corrections
*
Corresponding author. E-mail address: [email protected] (G. Spiga).
0096-3003/$ - see front matter Ó 2002 Elsevier Inc. All rights reserved. doi:10.1016/S0096-3003(02)00602-1
510
M. Galli et al. / Appl. Math. Comput. 146 (2003) 509–531
1. Introduction Kinetic theory is being used significantly in the last half century in order to understand evolution problems involving chemical reactions. Among the many contributions related to chemically reacting gas flows we may quote a pioneering work like [1] and some recent papers by the authors [2,3] where an extensive bibliography is given. In this respect, recombination/dissociation reactions are very important in the world of applications, but quite difficult to deal with from a kinetic point of view, and only few approaches are available in the literature, at least to our knowledge. They are based on either three-body collisions [4], or on a sequence of two-body collisions according to the transition-state theory [5]. It is just along the latter line that a Boltzmann-type model has been proposed in [6] for a mesoscopic description of a diatomic gas, considered as a three species mixture of atoms A (mass m), stable molecules A2 , and unstable molecules AH 2 (labeled in the sequel by indices 0, 1, 2, respectively). A key point of the model is that the third species plays the role of a transition state, created by atom-atom collisions only, and does not survive to any interaction. Since the mean collision times of the unstable molecules are asH sumed to be much shorter than all other interaction times, the AH 2 –A2 encounters may be neglected, and the mathematical problem is affected intrinsically by at least one small parameter, which makes the problem itself of singular perturbation type. It was possible to derive in [6], from the kinetic formulation, both conservation and moment equations at the macroscopic level. The essential step of a closure strategy was also addressed, and indeed self-consistent fluid dynamics equations were obtained either in some simplified physical situations, or in general in the asymptotic limit when the small parameter tends to zero. We shall consider here for simplicity the case of an isothermic reaction in a thermal bath in space homogeneous conditions, when the effective collision frequencies maij for an i–j encounter of type a are constant (see [6] for details). In this case the evolution is simply described by a set of three ordinary differential equations for the densities ni . Such a dynamical system retains the essential features of the singular perturbation problem. In particular the problem may be classified as singularly perturbed [7], or as a singular singular perturbation problem [8], since the dominant part of the collision operator has infinite null points, or, in other words, since the reduced problem obtained for vanishing small parameter has an infinity of solutions. This peculiar feature is shared by other important perturbation problems, like derivation of the hydrodynamic limit from the classical Boltzmann equation [9], and in that frame a famous asymptotic expansion, bearing the name of Hilbert, was devised. In this paper some singular perturbation methods suitable for the presented problem are illustrated, aiming at deriving and comparing approximate analytical solutions. The quite usual assumption of quasi-steady-state for the
M. Galli et al. / Appl. Math. Comput. 146 (2003) 509–531
511
excited species [10], considered already in [6], is taken into account. This is called stationary-state approximation, and represents the simplest approach for applications, but ignores completely the initial layer, and thus needs to be equipped with proper (fictitious) initial conditions in order to provide an accurate description outside the initial layer (the so-called bulk region). The singular perturbation analysis presented here has been inspired also by this problem of determining the appropriate initial conditions for such widely used approximation. The paper is organized as follows. The physical problem and the governing equations to be studied are presented and discussed in Section 2, together with the quasi-steady-state approach. The following section deals with the derivation of an analytical solution, first-order approximation with respect to the small parameter, in the spirit of the Hilbert asymptotic expansions. We apply the well know OÕMalley–Hoppensteadt procedure [8,11] in which the solution is sought as a superposition of an outer expansion and an inner expansion, the latter vanishing when time tends to infinity. The inner components involve the stretched time variable and account for initial layer corrections. Section 4 is devoted to the reduction of the same set of equations to a regular singular perturbation problem [8], with appearance of slowly varying dependent variables, and to the exploitation of a different asymptotic procedure. More specifically, the kernel of the dominant part of the collision operator allows to introduce suitable projection operators which split the solution as the sum of a hydrodynamic and of a kinetic part [7], with the former showing a slow evolution. The equivalent singular perturbation problem leads of course, by standard asymptotic expansion, to the same approximate solutions as before, but lends itself to a different asymptotic treatment, which shares the same features of the Chapman–Enskog expansion of kinetic theory [9]. It consists in leaving the outer component of the hydrodynamic part of the solution unexpanded with respect to the small parameter, and it leads to a modified (compressed) asymptotic solution, as proposed in [12] and worked out in detail in [7]. Finally numerical comparison of the various approximate solutions is presented in Section 5, in order to emphasize the different trends of the error in the different approaches.
2. The model equations and the stationary-state approximation According to the adopted physical model, formation of stable molecules A2 from atoms A occurs in two steps according to two irreversible bimolecular chemical reactions. The first is a recombination process A þ A ! AH 2 , labeled by superscript a ¼ r, while the second is an inelastic scattering (a ¼ i) AH 2 þ P ! A2 þ P , where P is either A or A2 , and where the excitation energy is transformed into extra kinetic energy of the products (exothermic process). The
512
M. Galli et al. / Appl. Math. Comput. 146 (2003) 509–531
inverse reaction of dissociation occurs in turn via two possible irreversible bimolecular reactions, both with a ¼ d, namely A2 þ P ! 2A þ P and AH 2 þ P ! 2A þ P . Reaction rates for a binary i–j encounter of type a are expressed in terms of microscopic collision frequencies maij (relative speeds times cross sections [9]). Elastic scattering between any pair of species is also allowed. Without entering the physical details, for which the interested reader is referred to [6], the simplest possible case for an analytical manipulation is the one described in Section 1, in which the evolution of the densities ni is governed by a self-consistent three-dimensional dynamical system, reading as 8 r 2 d d 2 d d > < n_ 0 ¼ 2m00 n0 þ 2m01 n0 n1 þ 2m11 n1 þ 2m02 n0 n2 þ 2m12 n1 n2 ; d d 2 i i ð1Þ n_ 1 ¼ m01 n0 n1 m11 n1 þ m02 n0 n2 þ m12 n1 n2 ; > : n_ ¼ ðmd þ mi Þn n ðmd þ mi Þn n þ mr n2 : 2 02 02 0 2 12 12 1 2 00 0 Indeed the set (1) is of clear phenomenological interpretation; the vector field on the right hand side represents the collision operator, and the quadratic nonlinearities reflect the fact that all collisions are binary. Initial conditions are taken of the form n0 ð0Þ ¼ a n1 ð0Þ ¼ b n2 ð0Þ ¼ c
ð2Þ
and it is easy to see that (1) implies n_ 0 þ 2n_ 1 þ 2n_ 2 ¼ 0, which quantifies the obvious physical requirement that the total mass is constant, since it is conserved in all underlying reactions. More explicitly n0 þ 2n1 þ 2n2 ¼ a þ 2b þ 2c ¼ q:
ð3Þ
a Due to the presence of unstable molecules AH 2 , all collision frequencies mk2 , a k ¼ 0; 1, in (1) are much larger than any of the mkl , k; l ¼ 0; 1. We may introduce a typical density n, typical slow and fast collision frequencies m and mH , with m=mH ¼ 1, and then we can measure densities in units of n (~ ni ¼ ni =n), frequencies makl , k; l 6¼ 2 in units of m (~ makl ¼ makl =m), frequencies mak2 in units of mH (~ mak2 ¼ mak2 =mH ), and time in units of nm inverse (~t ¼ nmt). Such adimensionalization yields a dimensionless dynamical system which differs from (1) only in that the mk2 coefficients are affected by a ‘‘large’’ factor 1=. We are led thus, as anticipated, to a singular perturbation problem [8,13–15]. With the usual abuse of notation we remove now all tildas. Moreover, we take advantage of (3) in order to reduce the dimensionality of the problem by eliminating n0 . In this way, mass conservation is exactly taken care of ‘‘a priori’’, before introducing any approximation techniques. There results ( n_ 1 ¼ A1 ðn1 ; n2 Þ þ 1 B1 ðn1 ; n2 Þ; ð4Þ n_ 2 ¼ A2 ðn1 ; n2 Þ þ 1 B2 ðn1 ; n2 Þ;
M. Galli et al. / Appl. Math. Comput. 146 (2003) 509–531
513
where the vectors A ¼ ðA1 ; A2 Þ and B ¼ ðB1 ; B2 Þ represent respectively the nondominant and the dominant part of the collision operator. They are given by A1 ¼ ð2md01 md11 Þn21 þ 2md01 n1 n2 qmd01 n1 ; A2 ¼ 4mr00 n21 þ 8mr00 n1 n2 þ 4mr00 n22 4qmr00 n1 4qmr00 n2 þ mr00 q2
ð5Þ
and B1 ¼ ð2mi02 mi12 Þn1 n2 2mi02 n22 þ qmi02 n2 ; B2 ¼ ð2md02 þ 2mi02 md12 mi12 Þn1 n2 þ 2ðmd02 þ mi02 Þn22 qðmd02 þ mi02 Þn2 :
ð6Þ
The phase space for the set (4) of ODE is the triangle of the ðn1 ; n2 Þ plane defined by the inequalities 0 6 n1 6 q=2, 0 6 n2 6 q=2, n1 þ n2 6 q=2. Fixed points for (4) have been determined in [6]. There is indeed a unique equilibrium internal to the phase space, determined by a suitable transcendental equation (representing mass action law for the chemical process). Moreover, another equilibrium exists on the border, given by n1 ¼ 0, n2 ¼ q=2, but it will be excluded from any further consideration since it corresponds to the completely unphysical situation of absence of any atom and of any stable molecule. Its formal existence in the model is due to the fact that other de-excitation channels have been disregarded. Anyhow, even from a mathematical point of view, it turns out to be unstable [6]. We will assume then from now on that n0 þ n1 > 0, i.e. n1 þ 2n2 < q. The singular perturbation problem (4) is also singularly perturbed [7]. In fact the limiting equations B1 ¼ 0, B2 ¼ 0 admit evidently at least the 11 solutions n2 ¼ 0, n1 arbitrary. The kernel of the dominant operator B is infinite-dimensional, whereas any positive value of implies a single null point for the full operator B þ A. Asymptotic methods will be applied later for an approximate solution of this singular singular problem. In order to keep all manipulations analytical, all dimensionless collision frequencies will be set in the sequel equal to unity. Then (5) and (6) become A1 ¼ n21 þ 2n1 n2 qn1 ; A2 ¼ 4n21 þ 8n1 n2 þ 4n22 4qn1 4qn2 þ q2
ð7Þ
and B1 ¼ n2 ðq n1 2n2 Þ; B2 ¼ 2n2 ðq n1 2n2 Þ:
ð8Þ
For later use, we report here on the fixed points for this special option. From (1), after excluding the unphysical point ð0; 0; q=2Þ, we get
514
M. Galli et al. / Appl. Math. Comput. 146 (2003) 509–531
n2 ¼
n20 2ðn0 þ n1 Þ
ð9Þ
2 2 plus the algebraic pffiffiffi equation n0 2n0 n1 2n1 ¼ 0, with unique admissible root n0 ¼ n1 ð1 þ 3Þ. Upon taking into account the first integral (3), the unique fixed point is then explicitly given by pffiffiffi 1þ 3 1 p ffiffi ffi pffiffiffi pffiffiffi nH ¼ q; nH q; nH q: ð10Þ 0 1 ¼ 2 ¼ 3 þ 3 þ 2 3 þ 3 þ 2 3 þ 3 þ 2 pffiffiffi pffiffiffi In the limiting case ¼ 0 the point (10) collapses to ðq= 3; q=ð3 þ 3Þ; 0Þ, but fixed points of (1) degenerate, since the only equation ðn0 þ n1 Þn2 ¼ 0 is left. Having assumed n0 þ n1 > 0, this is solved actually by (and only by) the 11 solutions
n0 ¼ q 2n1 ;
n1 arbitrary; n2 ¼ 0:
ð11Þ
The simplest algorithm for a possible approximate analytical solution of problem (1) is provided by the very well known and widely used stationarystate approximation [10], as discussed already in [6]. It amounts to a quasisteady-state assumption on the excited species, which reduces the dimension of the dynamical system by one since the very beginning. It is physically motivated by the fact that, if mk2 collision frequencies are very large, the third species adjusts almost immediately to the collision equilibrium in which the loss rate compensates instantaneously the gain rate due to recombination. From the standpoint of analytical manipulations, if n~2 ¼ n2 = is used instead of n2 in the dimensionless version of (1), disappears from the first two equations, and is left in the third only as a factor in front of n~_ 2 . Therefore, taking the limiting case ¼ 0, it reproduces exactly the stationary-state approximation. In any case, this approximation consists quantitatively in setting the third component of the vector field in (1) equal to zero, which yields at once a linear algebraic equation for n2 , solved by (9). This value of n2 is necessarily OðÞ, and is used in the remaining equations. It is clear that this procedure ignores completely the initial layer and can supply only an outer component hopefully fitting the exact solution in the bulk region. In this respect it is remarkable that it prescribes relaxation to the exact equilibrium, since the approximation does not affect the collision terms. However, the stationary-state algorithm is not self-consistent, in the sense that the proper (fictitious) initial conditions to be applied can not be deduced in terms of the actual initial conditions by the algorithm itself. In order to make this approximation explicit, we start from (4) with n0 expressed by the first integral (3). We get n2 ¼
ðq 2n1 2n2 Þ2 2ðq n1 2n2 Þ
ð12Þ
and the resulting quadratic equation for n2 has only one admissible solution
M. Galli et al. / Appl. Math. Comput. 146 (2003) 509–531
n2 ¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi q n1 þ 2ðq 2n1 Þ ðq n1 Þ2 þ 4n1 ðq 2n1 Þ 4ð1 þ Þ
515
ð13Þ
to be used in the first of (4). In this way the dynamical system collapses to the single first-order autonomous equation for n1 n_ 1 ¼
7 6 42
n21
3q
n 2 1 2ð1 þ Þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi q ð5 þ 2Þn1 q2 ðq n1 Þ2 þ 4n1 ðq 2n1 Þ þ þ 2 4ð1 þ Þ 4ð1 þ Þ2 4ð1 þ Þ
2
ð14Þ
which is amenable to quadrature, and can be solved numerically without any problem. It is easy to check that the right hand side of (14) vanishes for n1 ¼ nH 1 as given by (10). Of course, no information is available on the effective initial conditions for such stationary-state equation, but this will be obtained as a byproduct from the singular perturbation methods proposed in the next sections, which will provide a detailed and mathematically appropriate analysis of the initial layer. 3. Hilbert expansion The first perturbation method we present is based on the Hilbert expansion. We consider the singular and singularly perturbed problem (4), rewritten as 8 dn1 > > ¼ f1 ðn1 ; n2 ; Þ ¼ ðn21 þ 2n1 n2 qn1 Þ n1 n2 2n22 þ qn2 ; > > < dt dn2 ð15Þ ¼ f2 ðn1 ; n2 ; Þ ¼ ð4n21 þ 8n1 n2 þ 4n22 4qn1 4qn2 þ q2 Þ > > > dt > : þ 2n1 n2 þ 4n22 2qn2 with initial conditions n1 ð0Þ ¼ b;
n2 ð0Þ ¼ c
ð16Þ
and we determine an asymptotic solution, uniformly valid with respect to time, by following the classical Hilbert expansion method which is often applied in kinetic theory for the asymptotic treatment of the Boltzmann equation [9]. The solution is sought in the composite form n1 ðt; Þ ¼ no1 ðt; Þ þ ni1 ðs; Þ; n2 ðt; Þ ¼ no2 ðt; Þ þ ni2 ðs; Þ
ð17Þ
as the sum of an outer (or bulk) solution and of an inner initial layer correction, depending on the stretched time variable s ¼ t=, which turns out to be the distinguished limit for this problem [15], and vanishing as s ! þ1. This
516
M. Galli et al. / Appl. Math. Comput. 146 (2003) 509–531
approach is equivalent to resorting to the matching principle [14], which would provide the same results. It is assumed that no1;2 , ni1;2 and the functions f1 , f2 can be expanded in power series of the small parameter X X no1 k n1k ðtÞ; ni1 k N1k ðsÞ; k
no2
X
k
ni2
k
n2k ðtÞ;
k
f1 ðn1 ; n2 ; Þ
X k
X
k
f1k ðn1 ; n2 Þ;
k
f2 ðn1 ; n2 ; Þ
ð18Þ
k N2k ðsÞ;
X
k
f2k ðn1 ; n2 Þ;
k
1 f1k ¼ k! 1 f2k ¼ k!
dk f1 dk dk f2 dk
; ¼0
ð19Þ
¼0
such that the terms in series (18) satisfy the initial conditions n10 ð0Þ þ N10 ð0Þ ¼ b;
ð20Þ
n20 ð0Þ þ N20 ð0Þ ¼ c; n1k ð0Þ þ N1k ð0Þ ¼ 0; n2k ð0Þ þ N2k ð0Þ ¼ 0;
k ¼ 1; 2; . . .
The outer solution is requested to satisfy (15) for every t > 0 8 dno > < 1 ¼ f1 ðno1 ; no2 ; Þ; dt o > dn : 2 ¼ f ðno ; no ; Þ; 2 1 2 dt
ð21Þ
ð22Þ
and this implies that the initial layer correction must be the asymptotically vanishing solution of 8 i dn > > < 1 ¼ f1 ½no1 ðs; Þ þ ni1 ðs; Þ; no2 ðs; Þ þ ni2 ðs; Þ; f1 ½no1 ðs; Þ; no2 ðs; Þ; ; ds i > dn > : 2 ¼ f2 ½no ðs; Þ þ ni ðs; Þ; no ðs; Þ þ ni ðs; Þ; f2 ½no ðs; Þ; no ðs; Þ; : 1 1 2 2 1 2 ds ð23Þ By inserting the expansions (18) and (19) into Eqs. (22) and (23) and by equating the terms with the same power in , we obtain a sequence of equations in the unknown terms of order k ¼ 0; 1; . . . of the expanded solutions, to be solved with the initial conditions (20) and (21).
M. Galli et al. / Appl. Math. Comput. 146 (2003) 509–531
517
3.1. Zero-order solution Let us first consider the reduced problem which occurs when we confine ourselves to the simplest possible approximation for representing the solution in powers of , namely we truncate the series (18) to the order zero, or, in other words, we consider their limiting values for ¼ 0. This leads to the limiting or zero-order solution of problem (15). The terms of order k ¼ 0 of the outer solution which are deduced from (22) must satisfy the nonlinear algebraic system 0 ¼ f1 ðn10 ; n20 ; 0Þ ¼ n20 ðq n10 2n20 Þ;
ð24Þ
0 ¼ f2 ðn10 ; n20 ; 0Þ ¼ n20 ð2q þ 2n10 þ 4n20 Þ:
ð25Þ
Its solutions are represented by the 11 family ðn10 ; 0Þ with arbitrary n10 and by the unstable point ð0; q=2Þ to be discarded as discussed above. Owing to the singularity of the kernel of the limiting vector field, this solution, which gives n20 ¼ 0 and leaves n10 undetermined, is then useless. Nevertheless, n10 ðtÞ can be calculated by pushing the expansion one step further, even though the new terms are not needed in the representation of the solution, and by deriving a compatibility condition for the solution of the outer problem of order k ¼ 1. In fact, by setting n20 ¼ 0, this -order outer problem is defined by the singular system dn10 ¼ f11 ¼ ðq n10 Þn21 qn10 þ n210 ; dt
ð26Þ
0 ¼ f21 ¼ ð2q þ 2n10 Þn21 4qn10 þ 4n210 þ q2 :
ð27Þ
This linear algebraic set of equations for the unknowns n11 and n21 is singular, like (24) and (25), in the sense that now it leaves n11 undetermined and allows a solution for n21 only if the two inhomogeneous terms are properly related. Specifically, the solution is 2
n21 ðtÞ ¼
½2n10 ðtÞ q ; 2½q n10 ðtÞ
ð28Þ
under the condition that n10 satisfies the following nonlinear differential equation: dn10 q ¼ 3n210 3qn10 þ ; 2 dt
ð29Þ
determining the unknown n10 which was left unspecified at the previous step. It is of Bernoulli type, and has the solution
518
M. Galli et al. / Appl. Math. Comput. 146 (2003) 509–531
pffiffiffi q½6n10 ð0Þ ð3 3Þq pffiffi pffiffiffi pffiffiffi pffiffi n10 ðtÞ ¼ pffiffiffi
6 3n10 ð0Þ 1 e 3qt þ 3q 1 3 þ ð1 þ 3Þe 3qt pffiffiffi 3 3 q; þ 6
ð30Þ
where n10 ð0Þ ¼ b N10 ð0Þ is the initial condition that is still unknown and must be determined by matching this outer solution with the limiting inner one. Therefore, let us now consider the limiting inner system, which is related to the terms of order zero in (23). It is dN10 ¼ f1 ðn10 ð0Þ þ N10 ðsÞ; n20 ð0Þ þ N20 ðsÞ; 0Þ f1 ðn10 ð0Þ; n20 ð0Þ; 0Þ ds 2 N10 N20 ; ¼ ½q n10 ð0Þ N20 2N20
ð31Þ
dN20 ¼ f2 ðn10 ð0Þ þ N10 ðsÞ; n20 ð0Þ þ N20 ðsÞ; 0Þ f2 ðn10 ð0Þ; n20 ð0Þ; 0Þ ds 2 þ 2N10 N20 ¼ 2½q n10 ð0Þ N20 þ 4N20
ð32Þ
and, owing again to its singularity, Eq. (31) can be replaced by dN10 1 dN20 ¼ : 2 ds ds
ð33Þ
By using the initial conditions N10 ð0Þ ¼ b n10 ð0Þ, N20 ð0Þ ¼ c and requiring that lim N10 ðsÞ ¼ lim N20 ðsÞ ¼ 0;
s!1
s!1
one obtains from (33) n10 ð0Þ ¼ b þ
c k1 ¼q ; 2 2
with k1 ¼ 2q 2b c > 0;
1 N10 ðsÞ ¼ N20 ðsÞ: 2
ð34Þ ð35Þ
Thus, by inserting this last result into (32) and by integrating, the limiting inner solution is determined as follows: k1 N10 ðsÞ ¼ k1 2 c 3 ek 1 s þ 6 N20 ðsÞ ¼ k1 c
k1 3 ek1 s þ 3
ð36Þ ð37Þ
for c > 0, and N10 ðsÞ ¼ N20 ðsÞ ¼ 0 in the special case c ¼ 0. Moreover, by inserting into (30) the initial condition determined in (34), one also obtains the final result for n10 ðtÞ
M. Galli et al. / Appl. Math. Comput. 146 (2003) 509–531
n10 ðtÞ ¼ e
pffiffi h 3qt
pffiffiffi q 3 3 q pffiffiffii pffiffiffi þ 6 3 þ 3
pffiffi6q ð3þ 3Þq3k1
519
ð38Þ
and the following composite zero-order approximation to the solution of problem (15) ð0Þ
n1 ðt; Þ ¼ n10 ðtÞ þ N10 ðt=Þ; ð0Þ
n2 ðt; Þ ¼ N20 ðt=Þ:
ð39Þ
3.2. Solution of order The -order outer term n21 ðtÞ is given by (28), while n11 ðtÞ will be determined, as in the case of the zero-order solution, through a compatibility condition for the outer system of order 2, and by the matching procedure with the correspondent inner solution N11 ðsÞ. The inner terms N11 , N21 must satisfy the following singular, linear system obtained by equating the terms of order k ¼ 1 in (23) dN11 ¼ a1 ðsÞN11 þ a2 ðsÞN21 þ A0 ðsÞ; ds
ð40Þ
dN21 ¼ 2a1 ðsÞN11 2a2 ðsÞN21 þ B0 ðsÞ ds
ð41Þ
with a1 ðsÞ ¼ N20 ðsÞ;
a2 ðsÞ ¼ k1 =2 4N20 ðsÞ N10 ðsÞ
and where A0 ðsÞ, B0 ðsÞ are known functions of the zero-order solution, also containing the initial data of the -order outer terms h i 2 A0 ðsÞ ¼ sn_ 10 ð0Þ þ n11 ð0Þ N20 n21 ð0ÞðN10 þ 4N20 Þ þ N10 þ ðq k1 ÞN10 þ 2ðq k1 =2ÞN20 þ 2N10 N20 ; h i 2 B0 ðsÞ ¼ 2 sn_ 10 ð0Þ þ n11 ð0Þ N20 þ 2n21 ð0ÞðN10 þ 4N20 Þ þ 4N10 2 þ 4ðq k1 ÞN20 þ 4N20 þ 4ðq k1 ÞN10 þ 8N10 N20 :
ð42Þ
ð43Þ
Eqs. (40) and (41) are related by dN11 1 dN21 B0 ðsÞ þ ¼ A0 ðsÞ þ 2 ds 2 ds 2 2 ¼ 3qN10 þ 4qN20 3k1 ðN10 þ N20 Þ þ 6N10 N20 þ 3N10 þ 2N20 2 ¼ ða 4b cÞN10 N10
ð44Þ
520
M. Galli et al. / Appl. Math. Comput. 146 (2003) 509–531
and a formal integration of this equation, with initial conditions satisfying (21) for k ¼ 1, yields 1 N11 ðsÞ þ n11 ð0Þ þ ½N21 ðsÞ þ n21 ð0Þ ¼ IðsÞ 2 Z s Z s 0 0 2 ¼ ða 4b cÞ N20 ðs Þ ds N20 ðs0 Þ ds0 : 0
ð45Þ
0
Since the initial layer correction terms must vanish as s ! 1, we have n11 ð0Þ ¼ Iðþ1Þ
n21 ð0Þ c q2 ¼ Iðþ1Þ þ b þ 2 2 2k1
ð46Þ
and by using the solution (36) of N10 ðsÞ it is obtained k1 k2 ½k1 s log j ðk1 =c 3Þek1 s þ 3 j 12½ðk1 =c 3Þek1 s þ 3 c k2 log j k1 =c j þ 12
IðsÞ ¼
ð47Þ
with 1 k2 ¼ 6
4 11 1 a b c : 3 3 2
It has the limit lim IðsÞ ¼ Ið1Þ ¼ k2 log
s!þ1
2a þ 2b k1
þ
c 12
and therefore the initial value of n11 is given by
2a þ 2b 7 q2 : n11 ð0Þ ¼ k2 log þbþ c k1 12 2k1
ð48Þ
ð49Þ
We can now derive explicitly the solution of the inner system. From (45) we have N21 ðsÞ ¼ 2½IðsÞ Ið1Þ N11 ðsÞ
ð50Þ
and substitution into (40) provides a differential equation for N11 with the solution Z s
Z s ða1 2a2 Þ ds0 þ f2a2 ½Iðs0 Þ Ið1Þ N11 ðsÞ ¼ n11 ð0Þ exp 0 0 Z s
þ A0 ðs0 Þg exp ða1 2a2 Þ ds00 ds0 : ð51Þ s0
M. Galli et al. / Appl. Math. Comput. 146 (2003) 509–531
521
Finally, let us determine the term n11 ðtÞ by using the singular, nonlinear system of order 2 which results by equating the terms with k ¼ 2 in (22). It is given by dn11 ¼ f12 ¼ ðq n10 Þn22 þ ð2n10 qÞn11 þ 2n10 n21 n11 n21 2n221 ; dt dn21 ¼ f22 ¼ 2ðn10 qÞn22 þ ð8n10 4qÞn11 þ ð8n10 4qÞn21 dt þ 2n11 n21 þ 4n221
ð52Þ
ð53Þ
and has a unique solution for n22 ðtÞ under the condition that n11 satisfies dn11 1 dn21 ¼ 3½2n10 ðtÞ q n11 þ 2½3n10 ðtÞ q n21 ðtÞ : 2 dt dt
ð54Þ
By taking into account that n10 ðtÞ, n21 ðtÞ and their derivatives are known, Eq. (54) may be rewritten as dn11 ¼ aðtÞn11 þ bðtÞ dt with aðtÞ ¼ 3ð2n10 qÞ; 2n10 q 2n10 3q dn10 bðtÞ ¼ ð3n10 qÞð2n10 qÞ þ n10 q 4ðn10 qÞ dt and solved with the initial condition (49), with the result
Z t
2a þ 2b 7 q2 aðt0 Þ dt0 þbþ c exp n11 ðtÞ ¼ k2 log k1 12 2k1 0 Z t
Z t þ bðt0 Þ exp aðnÞ dn dt0 : 0
ð55Þ
t0
The -order approximate solution of problem (15) is thus determined as ð1Þ
n1 ðt; Þ ¼ n10 ðtÞ þ N10 ðt=Þ þ ½n11 ðtÞ þ N11 ðt=Þ ;
ð56Þ
ð1Þ
n2 ðt; Þ ¼ N20 ðt=Þ þ ½n21 ðtÞ þ N21 ðt=Þ : ð0ÞH
ð0ÞH
When tp! ffiffiffi 1 the zero-order solution converges to the point ðn1 ; n2 Þ ¼ ðqð3 3Þ=6; 0Þ, and analogously the first-order solution converges to a point showing OðÞ corrections with respect to this. In fact the final equilibria reached by the previous expansions are necessarily asymptotic representations of the exact one, given by (10). The extension of this solution to higher-order approximations is now a matter of straightforward, though cumbersome, calculations. In fact, the expansions of f1 , f2 in power series of lead to linear equations for the unknown
522
M. Galli et al. / Appl. Math. Comput. 146 (2003) 509–531
terms of order k ¼ 2; 3; . . . having the same coefficient matrix, i.e. the singular Jacobian matrix of the vector field, calculated at ¼ 0
q n10 4n20 n20 : Df ðn10 ; n20 ; 0Þ ¼ 2n20 2ðq n10 4n20 Þ Consequently, the same procedure used for the -order approximation can be applied to determine the outer and inner terms of order k P 2. 4. Fast and slow variables, compressed solution As seen in the previous section, the source of singularity, which forces consideration of the step n þ 1 and a relevant compatibility condition in order to determine the nth order solution, is the operator B, as given by (8). Following the detailed procedure devised in [7], we can easily realize that the null space of B, N ðBÞ, is made up by pairs n ¼ ðn1 ; n2 Þ in the phase space such that n2 ¼ 0. Analogously, the range of B, RðBÞ, is constituted by the pairs n ¼ ðn1 ; n2 Þ with n2 ¼ 2n1 . In order to single out the hydrodynamic variable and to proceed to the compressed (Chapman–Enskog) asymptotic expansion, it is crucial the determination of the projection operator P such that P ðnÞ 2 N ðBÞ and P 2 ðnÞ ¼ P ðnÞ for any n in the phase space, and P ðnÞ ¼ 0 8n 2 RðBÞ. In the present case, it is easily seen by inspection that P is determined as 1 P1 ðn1 ; n2 Þ ¼ n1 þ n2 ; 2 ð57Þ P2 ðn1 ; n2 Þ ¼ 0: The main feature of the operator P is that its application to Eq. (4) eliminates the dominant part of the collision operator, by virtue of the identity P ðBðnÞÞ ¼ 0
ð58Þ
for any n in the phase space. Since P and time derivative obviously commute, this leaves an equation with slow evolution for the quantity P ðnÞ, the hydrodynamic variable of [7]. The complementary projection operator Q is then 1 Q1 ðn1 ; n2 Þ ¼ n2 ; 2 ð59Þ Q2 ðn1 ; n2 Þ ¼ n2 with n ¼ P ðnÞ þ QðnÞ, P ðnÞ 2 N ðBÞ and QðnÞ 2 RðBÞ. Such decomposition introduces spontaneously new dependent variables 1 x 1 ¼ n1 þ n2 ; 2 ð60Þ x 2 ¼ n2 which are respectively the slow and the fast variable in the jargon of singular perturbation theory [8].
M. Galli et al. / Appl. Math. Comput. 146 (2003) 509–531
The set (15), recast in terms of the variables (60), reads as 8 2 1 2 1 1 2 > < x_ 1 ¼ g1 ðx1 ; x2 Þ ¼ 3x1 þ 3x1 x2 4x2 3qx1 2qx2 þ 2q ; _x2 ¼ g2 ðx1 ; x2 ; Þ ¼ 4x21 þ 2ð1 þ 2Þx1 x2 þ ð3 þ Þx22 > : 4qx1 2ð1 þ Þqx2 þ q2 ;
523
ð61Þ
with phase space given by the triangle defined by x2 P 0, x1 P ð1=2Þx2 , 2x1 þ x2 6 q (the unphysical equilibrium corresponds to the vertex ðq=4; q=2Þ). It is only matter of some algebra to show that the unique fixed point in the admissible phase space corresponds to the point (10), and that the problem is now regular singular, since the limiting equations for ¼ 0 have the unique pffiffiffi fixed point of physical meaning x1 ¼ qð3 3Þ=6, x2 ¼ 0, which corresponds to the limiting value of (10) when ! 0. We can apply then the standard composite asymptotic expansion x1 ðt; Þ ¼ xo1 ðt; Þ þ xi1 ðs; Þ; x2 ðt; Þ ¼ xo2 ðt; Þ þ xi2 ðs; Þ
ð62Þ
and the initial conditions x1 ð0Þ ¼ b þ ðc=2Þ, x2 ð0Þ ¼ c. Outer and inner components are expanded in power series of X X xoj ¼ k xjk ðtÞ; xij ¼ k Xjk ðsÞ; j ¼ 1; 2 ð63Þ k
k
with initial values c x10 ð0Þ þ X10 ð0Þ ¼ b þ ; x20 ð0Þ þ X20 ð0Þ ¼ c; 2 x2k ð0Þ þ X2k ð0Þ ¼ 0; k ¼ 1; 2; . . . x1k ð0Þ þ X1k ð0Þ ¼ 0;
ð64Þ
and with initial layer corrections Xij ðsÞ bound to vanish when s ! 1. We will not give all detailed manipulations. The main variation with respect to the previous section is that now each step determines the unknowns relevant to that step, reducing the amount of standard but lengthy algebra. Another advantage is that for the present approach several results about asymptoticity of the series and error estimates are well established in [8,11], providing an adequate mathematical setting to the present formal developments. Of course, after truncation either to the zero order or to the first order, upon going back to the density variables n1 and n2 , the two analytical approximate solutions coincide. It suffices to mention here that in the solution of the set (61) the zero-order outer and inner components have to satisfy respectively dx10 1 1 q2 ¼ 3x210 þ 3x10 x20 x220 3qx20 qx20 þ ; 4 2 dt 2
ð65Þ
0 ¼ 2x10 x20 þ 3x220 2qx20 ;
ð66Þ
524
M. Galli et al. / Appl. Math. Comput. 146 (2003) 509–531
dX10 ¼ 0; ds
ð67Þ
dX20 2 ¼ 2x10 ð0ÞX20 þ 2X10 X20 þ 3X20 2qX20 : ds
ð68Þ
From (66) one has x20 ðtÞ ¼ 0 and therefore also n20 ðtÞ ¼ 0, n10 ¼ x10 . Since X10 must vanish asymptotically, from (67) it is obtained X10 ðsÞ ¼ 0. Thus we have N10 ¼ N20 =2, already obtained in (35), and moreover the initial condition c c ð69Þ x10 ð0Þ ¼ X10 ð0Þ þ b þ ¼ b þ 2 2 for Eq. (65), that after suppression of x20 reads as dx10 q2 ¼ 3x210 3qx10 þ ; dt 2
ð70Þ
and is coincident with (29) since x10 ¼ n10 . Finally, by replacing x10 ð0Þ and suppressing X10 in (68) one obtains dX20 2 ¼ k1 X20 þ 3X20 ; ds
ð71Þ
whose solution satisfying the initial condition X20 ð0Þ ¼ c is the same as N20 ðsÞ given by (37). In analogous fashion, the -order terms in the expansions (63) can be determined by solving the following sets of equations, that result by equating in system (61) the terms of the same order dx11 q ¼ 3ð2x10 qÞx11 þ 3x10 ð72Þ x21 ; 2 dt 0 ¼ 2ðx10 qÞx21 þ 4x210 4qx10 þ q2 ; a dX11 c 1 2 ¼ 2b X20 X20 ; 2 2 4 ds dX21 e 0 ðsÞ ¼ 2X20 X11 þ ð6X20 kÞX21 þ B ds
ð73Þ ð74Þ ð75Þ
with e 0 ðsÞ ¼ ½2_x10 ð0Þs þ 2x11 ð0Þ þ 6x21 ð0Þ 2ða þ cÞ X20 þ X 2 : B 20 Eq. (73) gives x21 ðtÞ ¼ n21 ðtÞ already obtained in (28), and the solution of (74) satisfying the asymptoticity condition requested to the initial layer corrections is X11 ðsÞ ¼ IðsÞ Ið1Þ, where IðÞ are the same integrals shown in (47) and (48). Substitution of these results into (72) and (75) and calculation of the respective initial conditions yield the following two linear first-order ordinary differential equations:
M. Galli et al. / Appl. Math. Comput. 146 (2003) 509–531
525
dx11 ð6x10 qÞð2x10 qÞ2 ¼ 3ð2x10 qÞx11 þ ; dt 4ðq x10 Þ x11 ð0Þ ¼ X11 ð0Þ ¼ Ið1Þ;
ð76Þ
dX21 e 0 ðsÞ; ¼ ðX20 kÞX21 þ 2X20 ½IðsÞ Ið1Þ þ B ds ðq kÞ2 X21 ð0Þ ¼ x21 ð0Þ ¼ k
ð77Þ
whose solution yields the remaining unknown terms at the stage. Thus, by restoring the original density variables, the first-order approximate solution of (15) is determined as t 1 t 1 t 1 ð1Þ n1 ðt; Þ ’ x10 ðtÞ X20 þ x11 ðtÞ x21 ðtÞ þ X11 X21 ; 2 2 2 ð78Þ h t i ð1Þ þ x21 ðtÞ þ X21 : ð79Þ n2 ðt; Þ ’ X20 A rather tedious algebra allows to verify that it is indeed equivalent to the Hilbert approximation of the same order. As mentioned in Section 1, once the problem is reduced to the form (61) and the slow (hydrodynamic) part has been singled out, an alternative algorithm is provided by the compressed asymptotic expansion [7], by leaving the outer component of the slow variable unexpanded. This allows to determine xo1 by solving a unique -dependent differential equation, which hopefully retains more information with respect to a truncation after few powers of . We shall follow this procedure according to the recipe proposed in [12]. We write t
xo1 ðtÞ ¼ x1 ðtÞ; X k uk ðx1 Þ xo2 ðtÞ ¼
ð80Þ
k
and try to determine the unknown functions x1 and uk by requiring the outer solution to satisfy the set (61). We have ! X X dx1 ¼ g1 x1 ; k uk ðx1 Þ ¼ k g1;k ; ð81Þ dt k k
X du dx1 X du X X dxo2 ¼ ¼ k k k k j g1;j ¼ k g2;k dt d x dt d x 1 1 j k k k
with 1 1 1 g1;0 ¼ g1 ðx1 ; u0 ðx1 ÞÞ ¼ 3x21 þ 3x1 u0 u20 3qx1 qu0 þ q2 ; 4 2 2
ð82Þ
526
M. Galli et al. / Appl. Math. Comput. 146 (2003) 509–531
g2;0 ¼ g2 ðx1 u0 ðx1 Þ; 0Þ ¼ 2x1 u0 þ 3u20 2qu0 ;
dg1 u q g1;1 ¼ ¼ u1 ðx1 Þ 3x1 0 ; 2 d ¼0 2
dg2 g2;1 ¼ d ¼0 ¼ u1 ðx1 Þð2x1 þ 6u0 2qÞ þ 4x21 þ 4x1 u0 þ u20 2qu0 4qx1 þ q2 : By equating in (82) the coefficients of 0 and it is obtained the system 8 < g2;0 ¼ ½2x1 þ 3u0 2q u0 ¼ 0; ð83Þ du0 : g2;1 ¼ g1;0 dx1 which has the solution 2
u0 ¼ 0;
u1 ðx1 Þ ¼
ð2x1 qÞ ; 2ðx1 qÞ
ð84Þ
where x1 is determined by solving dx1 ¼ g1 ðx1 ; u1 ðx1 ÞÞ dt " # q2 q2 ð2x1 qÞ4 2 ¼ 3x1 ðx1 qÞ þ þ ðq 6x1 Þ x1 þ : 2 2 4ðx1 qÞ 16ðq x1 Þ ð85Þ The initial condition for (85) and the ones of the inner terms X1k , X2k are now related, up to the -order approximation, as follows: c x1 ð0Þ þ X10 ð0Þ þ X11 ð0Þ ¼ b þ ; 2
ð86Þ
u0 ½x1;0 ð0Þ þ X20 ð0Þ ¼ c;
ð87Þ
u1 ½x1;0 ð0Þ þ X21 ð0Þ ¼ 0;
ð88Þ
where x1;0 ð0Þ ¼ b þ c=2 is the initial value of x1 ðtÞ in the limit ! 0. They are calculated, as usual, by matching the outer solution with the inner one, which must satisfy the following system: 8 i dx > > < 1 ¼ g1 ½x1 þ xi1 ; u1 ðx1 Þ þ xi2 g1 ½x1 ; u1 ðx1 Þ ; ds ð89Þ i > dx > 2 i i : ¼ g2 ½x1 þ x1 ; u1 ðx2 Þ þ x2 ; g2 ½x1 ; u1 ðx1 Þ; : ds
M. Galli et al. / Appl. Math. Comput. 146 (2003) 509–531
527
By equating the terms of order zero and one with respect to , one has the same set provided by the previous approach, e.g. Eqs. (67), (68) and (74), (75) for Xi0 and Xi1 respectively, where the nonhomogeneous term in (75) is turned to e 0 ðsÞ ¼ ½2x_ 1 ð0Þs þ 6u1 ðx1 ð0ÞÞ 2ða þ cÞ X20 þ X 2 : B 0 20 It follows that the initial condition for (85) must be c c x1 ð0Þ ¼ b þ X10 ð0Þ X11 ð0Þ ¼ b þ þ Ið1Þ 2 2 and it can now be used to find the unexpanded outer component of the slow variable, by integration of the differential equation (85). In conclusion, in terms of the original density variables n1 , n2 the -order approximation achieved by using the compressed expansions algorithm is t 1 t 1 t 1 ðcÞ n1 ðt; Þ ’ x1 ðtÞ X20 þ u1 ðx1 Þ þ X11 X21 ; 2 2 2 ð90Þ ðcÞ
n2 ðt; Þ ’ X20
h t i þ u1 ðx1 Þ þ X21 :
t
ð91Þ
We are not reporting on the zero-order approximation since it is not difficult to check that the unexpanded part x1 ðtÞ would coincide with the unknown x10 of the uncompressed procedure, and consequently we would get the same results from the two algorithms.
5. Numerical results The results of the approximate asymptotic solutions derived in the previous Sections are now illustrated and discussed. We consider the dynamical system (15) together with initial conditions (16) and we first compute its ‘‘exact’’ (numerical) solution by means of an accurate ODE solver; the density n0 comes then from the first integral (3). In our numerical examples we start with the following reference numerical values ¼ 0:01;
a ¼ 0;
b ¼ 1;
c ¼ 4:
They have to be considered as dimensionless (for illustrative purposes only), and corresponding to arbitrary scales. In Fig. 1(a) we compare the temporal evolutions of n0 , n1 , n2 obtained by the exact solution (solid line), by the zero-order Hilbert expansion (dotted line), and by the first-order Hilbert expansion (dashed-dotted line). The initial layer clearly shows up, and it is evident that the first-order Hilbert approximation better reproduces the exact behaviour (the very slight differences can be hardly
528
M. Galli et al. / Appl. Math. Comput. 146 (2003) 509–531
6
6 n
n
0
0
5
5
4
4
3
3 n1
n1
2
2
1
1 n
0
n
2
0
0.02
0.04
0.06
0.08
(a)
0.1
0
2
0
0.02
0.04
0.06
0.08
0.1
(b)
Fig. 1. Reference case ¼ 0:01. (a) Exact solution (solid line), zero-order (dotted line) and firstorder (dashed-dotted line) Hilbert expansion; (b) stationary-state approximation (dotted line), outer expansion (dashed-dotted line) and exact solution (solid line).
recognized in the graph). We recall that, as already seen, the zero-order Hilbert expansion coincides analytically with the zero-order compressed expansion and with the zero-order solution of the fast–slow variables method, as well as, at first order, the Hilbert approximation and the fast–slow variables method give exactly the same analytical results. The widely used stationary-state approximation has been compared with the outer solution of the first-order Hilbert expansion in Fig. 1(b) where the temporal evolution of the densities ni is shown. The initial data are the same as for the outer solution in the Hilbert expansion, namely n1 ð0Þ ¼ n10 ð0Þ þ n11 ð0Þ; where n10 ð0Þ and n11 ð0Þ are given by (34) and (49). Similar results could be obtained by choosing initial data in such a way that the least square error between the exact and the stationary-state solutions is minimal. The solid line represents again the reference exact numerical solution, whereas the stationarystate approximation and the outer solution, which overlap almost entirely, are given by the dotted and the dashed-dotted lines respectively. We can notice that, apart from the initial layer of thickness OðÞ, the stationary-state approximation provides an accurate description of the evolution in the bulk region. In Fig. 2(a) we report the results as in Fig. 1(a) obtained by changing the initial data. A different trend with respect to the reference case can be observed in the evolution of the three number densities, and a wide range of behaviours
M. Galli et al. / Appl. Math. Comput. 146 (2003) 509–531 6
529
6 n0
n0
5
5
4
4
3
3
n
n1
1
2
2
1
1 n2
n 0
2
0
0.02
0.04
0.06 (a)
0.08
0.1
0
0
0.02
0.04
0.06
0.08
0.1
(b)
Fig. 2. (a) Exact solution and Hilbert expansions obtained for different initial data (a ¼ 0, b ¼ 2:5, c ¼ 2:5) and the same ¼ 0:01. (b) Exact solution (solid line), zero-order (dotted line) and firstorder (dashed-dotted line) Hilbert expansion obtained for ¼ 0:05 and initial data as in Fig. 1.
could be obtained by playing with initial data. Of course, the initial layer effect tends to disappear when the initial condition approaches the quasi-steady-state condition (9), as already pointed out in [6]. Indeed, it was already observed in (36) and (37) that the zero-order inner solution vanishes for c ¼ 0. The effect of varying is depicted in Fig. 2(b) where we plot the results obtained for ¼ 0:05 when we start with the reference initial data. As expected, a larger initial layer shows up, and the zero- and even first-order Hilbert approximations lead to worse results than in case ¼ 0:01. The first-order compressed expansion gives rise to solutions which differs, in principle, from the ones obtained by the first-order Hilbert approximation, but the discrepancies are not quantitatively appreciable in the graphs ðt; ni Þ. In Fig. 3(a) we then show the temporal trend (in logarithmic scale) of the error vectors, ex defined as the pointwise euclidean norm of the vector (n0 nex 0 , n1 n1 , ex ex n2 n2 ), where ni stands for the exact numerical values. We compare there the errors for the two values of used above versus time. The solid lines are relevant to the error in the first-order Hilbert expansion, whereas the dotted lines represent the error in the first-order compressed expansion. As we can see, the two methods do not show significant discrepancies in the results, although they come from two different asymptotic strategies. The compressed algorithm enhances a bit the accuracy in the bulk region. Fig. 3(b) shows instead the behaviour of the errors (with respect to the exact solution again) in the stationary-state approximation (solid line) and in the first-order Hilbert expansion (dotted line), with exclusion of the very initial region, where the former exhibits
530
10
M. Galli et al. / Appl. Math. Comput. 146 (2003) 509–531 0
4
3
3.5
eps=0.05 10
x 10
2
3 2.5
10
4
2 1.5
eps=0.01 10
6
1 0.5
10
8
0
0.2
0.4
0.6
0.8
0
0.005
0.2
(a)
0.4
0.6
0.8
(b)
Fig. 3. (a) Comparison of the error in the first-order Hilbert expansion (solid line) and in the firstorder compressed expansion (dotted line) obtained for ¼ 0:01 and 0.05 respectively (logarithmic scale). (b) Error in the stationary-state approximation (solid line) and in the Hilbert first-order expansion (dotted line) (case ¼ 0:01).
very high deviations. It is worth noticing that in the steady-state approximation the error tends to vanish for t large enough, since no truncation occurs in the fixed point of the problem with respect to the exact one. To summarize the results, we report finally in Table 1 the values of the density n1 computed at time instant t ¼ 0:01, t ¼ 0:001 by using both the ODE solver and the different asymptotic approximations described above. The values marked by are not significant since they are relevant to the stationarystate approximation out of its region of validity, namely in the initial layer.
Table 1 ð0Þ ð1Þ Computed values of density n1 obtained by zero- (n1 ) and first-order (n1 ) Hilbert expansion, firstðcÞ ðsÞ order compressed method (n1 ) and stationary-state approximation (n1 ) at two different time instants t ¼ 0:01 and t ¼ 0:001 ð0Þ
ð1Þ
ðcÞ
ðsÞ
n1
n1
n1
n1
n1
t ¼ 0:01
0.01 0.005 0.001
2.8912 2.8844 2.8790
3.3924 3.3834 3.3762
2.8961 2.8868 2.8795
2.8962 2.8869 2.8795
2.8915 2.8846 2.8790
t ¼ 0:001
0.01 0.005 0.001
1.6056 2.3716 2.9886
2.4688 3.1602 3.5526
1.6053 2.3720 2.9892
1.6054 2.3721 2.9892
3.0080 2.9950 2.9987
M. Galli et al. / Appl. Math. Comput. 146 (2003) 509–531
531
Acknowledgements This work has been performed in the framework of the activities sponsored by MIUR (Projects ‘‘Problemi matematici non lineari di propagazione e stabilit a nei modelli del continuo’’ for RR and ‘‘Problemi matematici delle teorie cinetiche’’ for the other three authors), by CNR, by GNFM, by University of Parma (Italy), and by the European TMR-Network ‘‘Asymptotic Methods in Kinetic Theory’’.
References [1] I. Prigogine, E. Xhrouet, On the perturbation of Maxwell distribution by chemical reactions in gases, Physica 15 (1949) 913–932. [2] A. Rossani, G. Spiga, A note on the kinetic theory of chemically reacting gases, Physica A 272 (1999) 563–573. [3] M. Groppi, G. Spiga, Kinetic approach to chemical reactions and inelastic transitions in a rarefied gas, J. Math. Chem. 26 (1999) 197–219. [4] J. Mc Lennan, Boltzmann equation for a dissociating gas, J. Stat. Phys. 57 (1989) 887–905. [5] Y. Yoshizawa, Wave structures of a chemically reacting gas by the kinetic theory of gases, in: J.L. Potter (Ed.), Rarefied Gas Dynamics, AIAA, New York, 1975, pp. 501–517. [6] M. Groppi, A. Rossani, G. Spiga, Kinetic theory of a diatomic gas with reactions of dissociation and recombination through a transition state, J. Phys. A: Math. Gen. 33 (2000) 8819–8833. [7] J. Mika, J. Banasiak, Singularly Perturbed Evolution Equations with Application to Kinetic Theory, World Scientific, Singapore, 1995. [8] R.E. OÕMalley Jr., Singular Perturbation Methods for Ordinary Differential Equations, Springer, New York, 1991. [9] C. Cercignani, The Boltzmann Equation and its Applications, Springer, New York, 1988. [10] A.W. Adamson, A Textbook of Physical Chemistry, Academic, New York, 1979. [11] D.R. Smith, Singular Perturbation Theory, Cambridge University Press, Cambridge, 1985. [12] J.R. Mika, A. Palczewski, Asymptotic analysis of singularly perturbed systems of ordinary differential equations, Comput. Math. Appl. 21 (1991) 13–32. [13] J.D. Cole, Perturbation Methods in Applied Mathematics, Blaisdell, Waltham, MA, 1968. [14] W. Eckhaus, Matched Asymptotic Expansions and Singular Perturbations, North Holland, Amsterdam, 1973. [15] A.H. Nayfeh, Perturbation Methods, Wiley, New York, 1993.
Singular perturbation techniques in the study of a diatomic gas with reactions of dissociation and recombination M. Galli a, M. Groppi a, R. Riganti b, G. Spiga a
a,*
Dipartimento di Matematica, Universit a di Parma, Via D’Azeglio 85, 43100 Parma, Italy b Dipartimento di Matematica, Politecnico di Torino, Corso Duca degli Abruzzi 24, 10129 Torino, Italy
Abstract Model equations for the description of chemical reactions of dissociation and recombination through a transition state constitute necessarily singular perturbation problems because of the very short lifetime of the excited molecules, compared to all other characteristic times. Asymptotic solutions are derived and tested versus an accurate numerical solution by resorting to different expansion algorithms, taking initial layer corrections into account. Since the singular problem turns out to be also singularly perturbed, techniques like Hilbert and Chapman–Enskog expansions, typical of the Boltzmann equation for gas kinetic theory, have been considered. The simple stationary state approximation, very popular in chemistry, is shown to provide very good results out of the initial layer if equipped with the proper initial conditions, as obtained from the present analysis. Ó 2002 Elsevier Inc. All rights reserved. Keywords: Dissociation/recombination reactions; Singular perturbation problems; Asymptotic expansions; Initial layer corrections
*
Corresponding author. E-mail address: [email protected] (G. Spiga).
0096-3003/$ - see front matter Ó 2002 Elsevier Inc. All rights reserved. doi:10.1016/S0096-3003(02)00602-1
510
M. Galli et al. / Appl. Math. Comput. 146 (2003) 509–531
1. Introduction Kinetic theory is being used significantly in the last half century in order to understand evolution problems involving chemical reactions. Among the many contributions related to chemically reacting gas flows we may quote a pioneering work like [1] and some recent papers by the authors [2,3] where an extensive bibliography is given. In this respect, recombination/dissociation reactions are very important in the world of applications, but quite difficult to deal with from a kinetic point of view, and only few approaches are available in the literature, at least to our knowledge. They are based on either three-body collisions [4], or on a sequence of two-body collisions according to the transition-state theory [5]. It is just along the latter line that a Boltzmann-type model has been proposed in [6] for a mesoscopic description of a diatomic gas, considered as a three species mixture of atoms A (mass m), stable molecules A2 , and unstable molecules AH 2 (labeled in the sequel by indices 0, 1, 2, respectively). A key point of the model is that the third species plays the role of a transition state, created by atom-atom collisions only, and does not survive to any interaction. Since the mean collision times of the unstable molecules are asH sumed to be much shorter than all other interaction times, the AH 2 –A2 encounters may be neglected, and the mathematical problem is affected intrinsically by at least one small parameter, which makes the problem itself of singular perturbation type. It was possible to derive in [6], from the kinetic formulation, both conservation and moment equations at the macroscopic level. The essential step of a closure strategy was also addressed, and indeed self-consistent fluid dynamics equations were obtained either in some simplified physical situations, or in general in the asymptotic limit when the small parameter tends to zero. We shall consider here for simplicity the case of an isothermic reaction in a thermal bath in space homogeneous conditions, when the effective collision frequencies maij for an i–j encounter of type a are constant (see [6] for details). In this case the evolution is simply described by a set of three ordinary differential equations for the densities ni . Such a dynamical system retains the essential features of the singular perturbation problem. In particular the problem may be classified as singularly perturbed [7], or as a singular singular perturbation problem [8], since the dominant part of the collision operator has infinite null points, or, in other words, since the reduced problem obtained for vanishing small parameter has an infinity of solutions. This peculiar feature is shared by other important perturbation problems, like derivation of the hydrodynamic limit from the classical Boltzmann equation [9], and in that frame a famous asymptotic expansion, bearing the name of Hilbert, was devised. In this paper some singular perturbation methods suitable for the presented problem are illustrated, aiming at deriving and comparing approximate analytical solutions. The quite usual assumption of quasi-steady-state for the
M. Galli et al. / Appl. Math. Comput. 146 (2003) 509–531
511
excited species [10], considered already in [6], is taken into account. This is called stationary-state approximation, and represents the simplest approach for applications, but ignores completely the initial layer, and thus needs to be equipped with proper (fictitious) initial conditions in order to provide an accurate description outside the initial layer (the so-called bulk region). The singular perturbation analysis presented here has been inspired also by this problem of determining the appropriate initial conditions for such widely used approximation. The paper is organized as follows. The physical problem and the governing equations to be studied are presented and discussed in Section 2, together with the quasi-steady-state approach. The following section deals with the derivation of an analytical solution, first-order approximation with respect to the small parameter, in the spirit of the Hilbert asymptotic expansions. We apply the well know OÕMalley–Hoppensteadt procedure [8,11] in which the solution is sought as a superposition of an outer expansion and an inner expansion, the latter vanishing when time tends to infinity. The inner components involve the stretched time variable and account for initial layer corrections. Section 4 is devoted to the reduction of the same set of equations to a regular singular perturbation problem [8], with appearance of slowly varying dependent variables, and to the exploitation of a different asymptotic procedure. More specifically, the kernel of the dominant part of the collision operator allows to introduce suitable projection operators which split the solution as the sum of a hydrodynamic and of a kinetic part [7], with the former showing a slow evolution. The equivalent singular perturbation problem leads of course, by standard asymptotic expansion, to the same approximate solutions as before, but lends itself to a different asymptotic treatment, which shares the same features of the Chapman–Enskog expansion of kinetic theory [9]. It consists in leaving the outer component of the hydrodynamic part of the solution unexpanded with respect to the small parameter, and it leads to a modified (compressed) asymptotic solution, as proposed in [12] and worked out in detail in [7]. Finally numerical comparison of the various approximate solutions is presented in Section 5, in order to emphasize the different trends of the error in the different approaches.
2. The model equations and the stationary-state approximation According to the adopted physical model, formation of stable molecules A2 from atoms A occurs in two steps according to two irreversible bimolecular chemical reactions. The first is a recombination process A þ A ! AH 2 , labeled by superscript a ¼ r, while the second is an inelastic scattering (a ¼ i) AH 2 þ P ! A2 þ P , where P is either A or A2 , and where the excitation energy is transformed into extra kinetic energy of the products (exothermic process). The
512
M. Galli et al. / Appl. Math. Comput. 146 (2003) 509–531
inverse reaction of dissociation occurs in turn via two possible irreversible bimolecular reactions, both with a ¼ d, namely A2 þ P ! 2A þ P and AH 2 þ P ! 2A þ P . Reaction rates for a binary i–j encounter of type a are expressed in terms of microscopic collision frequencies maij (relative speeds times cross sections [9]). Elastic scattering between any pair of species is also allowed. Without entering the physical details, for which the interested reader is referred to [6], the simplest possible case for an analytical manipulation is the one described in Section 1, in which the evolution of the densities ni is governed by a self-consistent three-dimensional dynamical system, reading as 8 r 2 d d 2 d d > < n_ 0 ¼ 2m00 n0 þ 2m01 n0 n1 þ 2m11 n1 þ 2m02 n0 n2 þ 2m12 n1 n2 ; d d 2 i i ð1Þ n_ 1 ¼ m01 n0 n1 m11 n1 þ m02 n0 n2 þ m12 n1 n2 ; > : n_ ¼ ðmd þ mi Þn n ðmd þ mi Þn n þ mr n2 : 2 02 02 0 2 12 12 1 2 00 0 Indeed the set (1) is of clear phenomenological interpretation; the vector field on the right hand side represents the collision operator, and the quadratic nonlinearities reflect the fact that all collisions are binary. Initial conditions are taken of the form n0 ð0Þ ¼ a n1 ð0Þ ¼ b n2 ð0Þ ¼ c
ð2Þ
and it is easy to see that (1) implies n_ 0 þ 2n_ 1 þ 2n_ 2 ¼ 0, which quantifies the obvious physical requirement that the total mass is constant, since it is conserved in all underlying reactions. More explicitly n0 þ 2n1 þ 2n2 ¼ a þ 2b þ 2c ¼ q:
ð3Þ
a Due to the presence of unstable molecules AH 2 , all collision frequencies mk2 , a k ¼ 0; 1, in (1) are much larger than any of the mkl , k; l ¼ 0; 1. We may introduce a typical density n, typical slow and fast collision frequencies m and mH , with m=mH ¼ 1, and then we can measure densities in units of n (~ ni ¼ ni =n), frequencies makl , k; l 6¼ 2 in units of m (~ makl ¼ makl =m), frequencies mak2 in units of mH (~ mak2 ¼ mak2 =mH ), and time in units of nm inverse (~t ¼ nmt). Such adimensionalization yields a dimensionless dynamical system which differs from (1) only in that the mk2 coefficients are affected by a ‘‘large’’ factor 1=. We are led thus, as anticipated, to a singular perturbation problem [8,13–15]. With the usual abuse of notation we remove now all tildas. Moreover, we take advantage of (3) in order to reduce the dimensionality of the problem by eliminating n0 . In this way, mass conservation is exactly taken care of ‘‘a priori’’, before introducing any approximation techniques. There results ( n_ 1 ¼ A1 ðn1 ; n2 Þ þ 1 B1 ðn1 ; n2 Þ; ð4Þ n_ 2 ¼ A2 ðn1 ; n2 Þ þ 1 B2 ðn1 ; n2 Þ;
M. Galli et al. / Appl. Math. Comput. 146 (2003) 509–531
513
where the vectors A ¼ ðA1 ; A2 Þ and B ¼ ðB1 ; B2 Þ represent respectively the nondominant and the dominant part of the collision operator. They are given by A1 ¼ ð2md01 md11 Þn21 þ 2md01 n1 n2 qmd01 n1 ; A2 ¼ 4mr00 n21 þ 8mr00 n1 n2 þ 4mr00 n22 4qmr00 n1 4qmr00 n2 þ mr00 q2
ð5Þ
and B1 ¼ ð2mi02 mi12 Þn1 n2 2mi02 n22 þ qmi02 n2 ; B2 ¼ ð2md02 þ 2mi02 md12 mi12 Þn1 n2 þ 2ðmd02 þ mi02 Þn22 qðmd02 þ mi02 Þn2 :
ð6Þ
The phase space for the set (4) of ODE is the triangle of the ðn1 ; n2 Þ plane defined by the inequalities 0 6 n1 6 q=2, 0 6 n2 6 q=2, n1 þ n2 6 q=2. Fixed points for (4) have been determined in [6]. There is indeed a unique equilibrium internal to the phase space, determined by a suitable transcendental equation (representing mass action law for the chemical process). Moreover, another equilibrium exists on the border, given by n1 ¼ 0, n2 ¼ q=2, but it will be excluded from any further consideration since it corresponds to the completely unphysical situation of absence of any atom and of any stable molecule. Its formal existence in the model is due to the fact that other de-excitation channels have been disregarded. Anyhow, even from a mathematical point of view, it turns out to be unstable [6]. We will assume then from now on that n0 þ n1 > 0, i.e. n1 þ 2n2 < q. The singular perturbation problem (4) is also singularly perturbed [7]. In fact the limiting equations B1 ¼ 0, B2 ¼ 0 admit evidently at least the 11 solutions n2 ¼ 0, n1 arbitrary. The kernel of the dominant operator B is infinite-dimensional, whereas any positive value of implies a single null point for the full operator B þ A. Asymptotic methods will be applied later for an approximate solution of this singular singular problem. In order to keep all manipulations analytical, all dimensionless collision frequencies will be set in the sequel equal to unity. Then (5) and (6) become A1 ¼ n21 þ 2n1 n2 qn1 ; A2 ¼ 4n21 þ 8n1 n2 þ 4n22 4qn1 4qn2 þ q2
ð7Þ
and B1 ¼ n2 ðq n1 2n2 Þ; B2 ¼ 2n2 ðq n1 2n2 Þ:
ð8Þ
For later use, we report here on the fixed points for this special option. From (1), after excluding the unphysical point ð0; 0; q=2Þ, we get
514
M. Galli et al. / Appl. Math. Comput. 146 (2003) 509–531
n2 ¼
n20 2ðn0 þ n1 Þ
ð9Þ
2 2 plus the algebraic pffiffiffi equation n0 2n0 n1 2n1 ¼ 0, with unique admissible root n0 ¼ n1 ð1 þ 3Þ. Upon taking into account the first integral (3), the unique fixed point is then explicitly given by pffiffiffi 1þ 3 1 p ffiffi ffi pffiffiffi pffiffiffi nH ¼ q; nH q; nH q: ð10Þ 0 1 ¼ 2 ¼ 3 þ 3 þ 2 3 þ 3 þ 2 3 þ 3 þ 2 pffiffiffi pffiffiffi In the limiting case ¼ 0 the point (10) collapses to ðq= 3; q=ð3 þ 3Þ; 0Þ, but fixed points of (1) degenerate, since the only equation ðn0 þ n1 Þn2 ¼ 0 is left. Having assumed n0 þ n1 > 0, this is solved actually by (and only by) the 11 solutions
n0 ¼ q 2n1 ;
n1 arbitrary; n2 ¼ 0:
ð11Þ
The simplest algorithm for a possible approximate analytical solution of problem (1) is provided by the very well known and widely used stationarystate approximation [10], as discussed already in [6]. It amounts to a quasisteady-state assumption on the excited species, which reduces the dimension of the dynamical system by one since the very beginning. It is physically motivated by the fact that, if mk2 collision frequencies are very large, the third species adjusts almost immediately to the collision equilibrium in which the loss rate compensates instantaneously the gain rate due to recombination. From the standpoint of analytical manipulations, if n~2 ¼ n2 = is used instead of n2 in the dimensionless version of (1), disappears from the first two equations, and is left in the third only as a factor in front of n~_ 2 . Therefore, taking the limiting case ¼ 0, it reproduces exactly the stationary-state approximation. In any case, this approximation consists quantitatively in setting the third component of the vector field in (1) equal to zero, which yields at once a linear algebraic equation for n2 , solved by (9). This value of n2 is necessarily OðÞ, and is used in the remaining equations. It is clear that this procedure ignores completely the initial layer and can supply only an outer component hopefully fitting the exact solution in the bulk region. In this respect it is remarkable that it prescribes relaxation to the exact equilibrium, since the approximation does not affect the collision terms. However, the stationary-state algorithm is not self-consistent, in the sense that the proper (fictitious) initial conditions to be applied can not be deduced in terms of the actual initial conditions by the algorithm itself. In order to make this approximation explicit, we start from (4) with n0 expressed by the first integral (3). We get n2 ¼
ðq 2n1 2n2 Þ2 2ðq n1 2n2 Þ
ð12Þ
and the resulting quadratic equation for n2 has only one admissible solution
M. Galli et al. / Appl. Math. Comput. 146 (2003) 509–531
n2 ¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi q n1 þ 2ðq 2n1 Þ ðq n1 Þ2 þ 4n1 ðq 2n1 Þ 4ð1 þ Þ
515
ð13Þ
to be used in the first of (4). In this way the dynamical system collapses to the single first-order autonomous equation for n1 n_ 1 ¼
7 6 42
n21
3q
n 2 1 2ð1 þ Þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi q ð5 þ 2Þn1 q2 ðq n1 Þ2 þ 4n1 ðq 2n1 Þ þ þ 2 4ð1 þ Þ 4ð1 þ Þ2 4ð1 þ Þ
2
ð14Þ
which is amenable to quadrature, and can be solved numerically without any problem. It is easy to check that the right hand side of (14) vanishes for n1 ¼ nH 1 as given by (10). Of course, no information is available on the effective initial conditions for such stationary-state equation, but this will be obtained as a byproduct from the singular perturbation methods proposed in the next sections, which will provide a detailed and mathematically appropriate analysis of the initial layer. 3. Hilbert expansion The first perturbation method we present is based on the Hilbert expansion. We consider the singular and singularly perturbed problem (4), rewritten as 8 dn1 > > ¼ f1 ðn1 ; n2 ; Þ ¼ ðn21 þ 2n1 n2 qn1 Þ n1 n2 2n22 þ qn2 ; > > < dt dn2 ð15Þ ¼ f2 ðn1 ; n2 ; Þ ¼ ð4n21 þ 8n1 n2 þ 4n22 4qn1 4qn2 þ q2 Þ > > > dt > : þ 2n1 n2 þ 4n22 2qn2 with initial conditions n1 ð0Þ ¼ b;
n2 ð0Þ ¼ c
ð16Þ
and we determine an asymptotic solution, uniformly valid with respect to time, by following the classical Hilbert expansion method which is often applied in kinetic theory for the asymptotic treatment of the Boltzmann equation [9]. The solution is sought in the composite form n1 ðt; Þ ¼ no1 ðt; Þ þ ni1 ðs; Þ; n2 ðt; Þ ¼ no2 ðt; Þ þ ni2 ðs; Þ
ð17Þ
as the sum of an outer (or bulk) solution and of an inner initial layer correction, depending on the stretched time variable s ¼ t=, which turns out to be the distinguished limit for this problem [15], and vanishing as s ! þ1. This
516
M. Galli et al. / Appl. Math. Comput. 146 (2003) 509–531
approach is equivalent to resorting to the matching principle [14], which would provide the same results. It is assumed that no1;2 , ni1;2 and the functions f1 , f2 can be expanded in power series of the small parameter X X no1 k n1k ðtÞ; ni1 k N1k ðsÞ; k
no2
X
k
ni2
k
n2k ðtÞ;
k
f1 ðn1 ; n2 ; Þ
X k
X
k
f1k ðn1 ; n2 Þ;
k
f2 ðn1 ; n2 ; Þ
ð18Þ
k N2k ðsÞ;
X
k
f2k ðn1 ; n2 Þ;
k
1 f1k ¼ k! 1 f2k ¼ k!
dk f1 dk dk f2 dk
; ¼0
ð19Þ
¼0
such that the terms in series (18) satisfy the initial conditions n10 ð0Þ þ N10 ð0Þ ¼ b;
ð20Þ
n20 ð0Þ þ N20 ð0Þ ¼ c; n1k ð0Þ þ N1k ð0Þ ¼ 0; n2k ð0Þ þ N2k ð0Þ ¼ 0;
k ¼ 1; 2; . . .
The outer solution is requested to satisfy (15) for every t > 0 8 dno > < 1 ¼ f1 ðno1 ; no2 ; Þ; dt o > dn : 2 ¼ f ðno ; no ; Þ; 2 1 2 dt
ð21Þ
ð22Þ
and this implies that the initial layer correction must be the asymptotically vanishing solution of 8 i dn > > < 1 ¼ f1 ½no1 ðs; Þ þ ni1 ðs; Þ; no2 ðs; Þ þ ni2 ðs; Þ; f1 ½no1 ðs; Þ; no2 ðs; Þ; ; ds i > dn > : 2 ¼ f2 ½no ðs; Þ þ ni ðs; Þ; no ðs; Þ þ ni ðs; Þ; f2 ½no ðs; Þ; no ðs; Þ; : 1 1 2 2 1 2 ds ð23Þ By inserting the expansions (18) and (19) into Eqs. (22) and (23) and by equating the terms with the same power in , we obtain a sequence of equations in the unknown terms of order k ¼ 0; 1; . . . of the expanded solutions, to be solved with the initial conditions (20) and (21).
M. Galli et al. / Appl. Math. Comput. 146 (2003) 509–531
517
3.1. Zero-order solution Let us first consider the reduced problem which occurs when we confine ourselves to the simplest possible approximation for representing the solution in powers of , namely we truncate the series (18) to the order zero, or, in other words, we consider their limiting values for ¼ 0. This leads to the limiting or zero-order solution of problem (15). The terms of order k ¼ 0 of the outer solution which are deduced from (22) must satisfy the nonlinear algebraic system 0 ¼ f1 ðn10 ; n20 ; 0Þ ¼ n20 ðq n10 2n20 Þ;
ð24Þ
0 ¼ f2 ðn10 ; n20 ; 0Þ ¼ n20 ð2q þ 2n10 þ 4n20 Þ:
ð25Þ
Its solutions are represented by the 11 family ðn10 ; 0Þ with arbitrary n10 and by the unstable point ð0; q=2Þ to be discarded as discussed above. Owing to the singularity of the kernel of the limiting vector field, this solution, which gives n20 ¼ 0 and leaves n10 undetermined, is then useless. Nevertheless, n10 ðtÞ can be calculated by pushing the expansion one step further, even though the new terms are not needed in the representation of the solution, and by deriving a compatibility condition for the solution of the outer problem of order k ¼ 1. In fact, by setting n20 ¼ 0, this -order outer problem is defined by the singular system dn10 ¼ f11 ¼ ðq n10 Þn21 qn10 þ n210 ; dt
ð26Þ
0 ¼ f21 ¼ ð2q þ 2n10 Þn21 4qn10 þ 4n210 þ q2 :
ð27Þ
This linear algebraic set of equations for the unknowns n11 and n21 is singular, like (24) and (25), in the sense that now it leaves n11 undetermined and allows a solution for n21 only if the two inhomogeneous terms are properly related. Specifically, the solution is 2
n21 ðtÞ ¼
½2n10 ðtÞ q ; 2½q n10 ðtÞ
ð28Þ
under the condition that n10 satisfies the following nonlinear differential equation: dn10 q ¼ 3n210 3qn10 þ ; 2 dt
ð29Þ
determining the unknown n10 which was left unspecified at the previous step. It is of Bernoulli type, and has the solution
518
M. Galli et al. / Appl. Math. Comput. 146 (2003) 509–531
pffiffiffi q½6n10 ð0Þ ð3 3Þq pffiffi pffiffiffi pffiffiffi pffiffi n10 ðtÞ ¼ pffiffiffi
6 3n10 ð0Þ 1 e 3qt þ 3q 1 3 þ ð1 þ 3Þe 3qt pffiffiffi 3 3 q; þ 6
ð30Þ
where n10 ð0Þ ¼ b N10 ð0Þ is the initial condition that is still unknown and must be determined by matching this outer solution with the limiting inner one. Therefore, let us now consider the limiting inner system, which is related to the terms of order zero in (23). It is dN10 ¼ f1 ðn10 ð0Þ þ N10 ðsÞ; n20 ð0Þ þ N20 ðsÞ; 0Þ f1 ðn10 ð0Þ; n20 ð0Þ; 0Þ ds 2 N10 N20 ; ¼ ½q n10 ð0Þ N20 2N20
ð31Þ
dN20 ¼ f2 ðn10 ð0Þ þ N10 ðsÞ; n20 ð0Þ þ N20 ðsÞ; 0Þ f2 ðn10 ð0Þ; n20 ð0Þ; 0Þ ds 2 þ 2N10 N20 ¼ 2½q n10 ð0Þ N20 þ 4N20
ð32Þ
and, owing again to its singularity, Eq. (31) can be replaced by dN10 1 dN20 ¼ : 2 ds ds
ð33Þ
By using the initial conditions N10 ð0Þ ¼ b n10 ð0Þ, N20 ð0Þ ¼ c and requiring that lim N10 ðsÞ ¼ lim N20 ðsÞ ¼ 0;
s!1
s!1
one obtains from (33) n10 ð0Þ ¼ b þ
c k1 ¼q ; 2 2
with k1 ¼ 2q 2b c > 0;
1 N10 ðsÞ ¼ N20 ðsÞ: 2
ð34Þ ð35Þ
Thus, by inserting this last result into (32) and by integrating, the limiting inner solution is determined as follows: k1 N10 ðsÞ ¼ k1 2 c 3 ek 1 s þ 6 N20 ðsÞ ¼ k1 c
k1 3 ek1 s þ 3
ð36Þ ð37Þ
for c > 0, and N10 ðsÞ ¼ N20 ðsÞ ¼ 0 in the special case c ¼ 0. Moreover, by inserting into (30) the initial condition determined in (34), one also obtains the final result for n10 ðtÞ
M. Galli et al. / Appl. Math. Comput. 146 (2003) 509–531
n10 ðtÞ ¼ e
pffiffi h 3qt
pffiffiffi q 3 3 q pffiffiffii pffiffiffi þ 6 3 þ 3
pffiffi6q ð3þ 3Þq3k1
519
ð38Þ
and the following composite zero-order approximation to the solution of problem (15) ð0Þ
n1 ðt; Þ ¼ n10 ðtÞ þ N10 ðt=Þ; ð0Þ
n2 ðt; Þ ¼ N20 ðt=Þ:
ð39Þ
3.2. Solution of order The -order outer term n21 ðtÞ is given by (28), while n11 ðtÞ will be determined, as in the case of the zero-order solution, through a compatibility condition for the outer system of order 2, and by the matching procedure with the correspondent inner solution N11 ðsÞ. The inner terms N11 , N21 must satisfy the following singular, linear system obtained by equating the terms of order k ¼ 1 in (23) dN11 ¼ a1 ðsÞN11 þ a2 ðsÞN21 þ A0 ðsÞ; ds
ð40Þ
dN21 ¼ 2a1 ðsÞN11 2a2 ðsÞN21 þ B0 ðsÞ ds
ð41Þ
with a1 ðsÞ ¼ N20 ðsÞ;
a2 ðsÞ ¼ k1 =2 4N20 ðsÞ N10 ðsÞ
and where A0 ðsÞ, B0 ðsÞ are known functions of the zero-order solution, also containing the initial data of the -order outer terms h i 2 A0 ðsÞ ¼ sn_ 10 ð0Þ þ n11 ð0Þ N20 n21 ð0ÞðN10 þ 4N20 Þ þ N10 þ ðq k1 ÞN10 þ 2ðq k1 =2ÞN20 þ 2N10 N20 ; h i 2 B0 ðsÞ ¼ 2 sn_ 10 ð0Þ þ n11 ð0Þ N20 þ 2n21 ð0ÞðN10 þ 4N20 Þ þ 4N10 2 þ 4ðq k1 ÞN20 þ 4N20 þ 4ðq k1 ÞN10 þ 8N10 N20 :
ð42Þ
ð43Þ
Eqs. (40) and (41) are related by dN11 1 dN21 B0 ðsÞ þ ¼ A0 ðsÞ þ 2 ds 2 ds 2 2 ¼ 3qN10 þ 4qN20 3k1 ðN10 þ N20 Þ þ 6N10 N20 þ 3N10 þ 2N20 2 ¼ ða 4b cÞN10 N10
ð44Þ
520
M. Galli et al. / Appl. Math. Comput. 146 (2003) 509–531
and a formal integration of this equation, with initial conditions satisfying (21) for k ¼ 1, yields 1 N11 ðsÞ þ n11 ð0Þ þ ½N21 ðsÞ þ n21 ð0Þ ¼ IðsÞ 2 Z s Z s 0 0 2 ¼ ða 4b cÞ N20 ðs Þ ds N20 ðs0 Þ ds0 : 0
ð45Þ
0
Since the initial layer correction terms must vanish as s ! 1, we have n11 ð0Þ ¼ Iðþ1Þ
n21 ð0Þ c q2 ¼ Iðþ1Þ þ b þ 2 2 2k1
ð46Þ
and by using the solution (36) of N10 ðsÞ it is obtained k1 k2 ½k1 s log j ðk1 =c 3Þek1 s þ 3 j 12½ðk1 =c 3Þek1 s þ 3 c k2 log j k1 =c j þ 12
IðsÞ ¼
ð47Þ
with 1 k2 ¼ 6
4 11 1 a b c : 3 3 2
It has the limit lim IðsÞ ¼ Ið1Þ ¼ k2 log
s!þ1
2a þ 2b k1
þ
c 12
and therefore the initial value of n11 is given by
2a þ 2b 7 q2 : n11 ð0Þ ¼ k2 log þbþ c k1 12 2k1
ð48Þ
ð49Þ
We can now derive explicitly the solution of the inner system. From (45) we have N21 ðsÞ ¼ 2½IðsÞ Ið1Þ N11 ðsÞ
ð50Þ
and substitution into (40) provides a differential equation for N11 with the solution Z s
Z s ða1 2a2 Þ ds0 þ f2a2 ½Iðs0 Þ Ið1Þ N11 ðsÞ ¼ n11 ð0Þ exp 0 0 Z s
þ A0 ðs0 Þg exp ða1 2a2 Þ ds00 ds0 : ð51Þ s0
M. Galli et al. / Appl. Math. Comput. 146 (2003) 509–531
521
Finally, let us determine the term n11 ðtÞ by using the singular, nonlinear system of order 2 which results by equating the terms with k ¼ 2 in (22). It is given by dn11 ¼ f12 ¼ ðq n10 Þn22 þ ð2n10 qÞn11 þ 2n10 n21 n11 n21 2n221 ; dt dn21 ¼ f22 ¼ 2ðn10 qÞn22 þ ð8n10 4qÞn11 þ ð8n10 4qÞn21 dt þ 2n11 n21 þ 4n221
ð52Þ
ð53Þ
and has a unique solution for n22 ðtÞ under the condition that n11 satisfies dn11 1 dn21 ¼ 3½2n10 ðtÞ q n11 þ 2½3n10 ðtÞ q n21 ðtÞ : 2 dt dt
ð54Þ
By taking into account that n10 ðtÞ, n21 ðtÞ and their derivatives are known, Eq. (54) may be rewritten as dn11 ¼ aðtÞn11 þ bðtÞ dt with aðtÞ ¼ 3ð2n10 qÞ; 2n10 q 2n10 3q dn10 bðtÞ ¼ ð3n10 qÞð2n10 qÞ þ n10 q 4ðn10 qÞ dt and solved with the initial condition (49), with the result
Z t
2a þ 2b 7 q2 aðt0 Þ dt0 þbþ c exp n11 ðtÞ ¼ k2 log k1 12 2k1 0 Z t
Z t þ bðt0 Þ exp aðnÞ dn dt0 : 0
ð55Þ
t0
The -order approximate solution of problem (15) is thus determined as ð1Þ
n1 ðt; Þ ¼ n10 ðtÞ þ N10 ðt=Þ þ ½n11 ðtÞ þ N11 ðt=Þ ;
ð56Þ
ð1Þ
n2 ðt; Þ ¼ N20 ðt=Þ þ ½n21 ðtÞ þ N21 ðt=Þ : ð0ÞH
ð0ÞH
When tp! ffiffiffi 1 the zero-order solution converges to the point ðn1 ; n2 Þ ¼ ðqð3 3Þ=6; 0Þ, and analogously the first-order solution converges to a point showing OðÞ corrections with respect to this. In fact the final equilibria reached by the previous expansions are necessarily asymptotic representations of the exact one, given by (10). The extension of this solution to higher-order approximations is now a matter of straightforward, though cumbersome, calculations. In fact, the expansions of f1 , f2 in power series of lead to linear equations for the unknown
522
M. Galli et al. / Appl. Math. Comput. 146 (2003) 509–531
terms of order k ¼ 2; 3; . . . having the same coefficient matrix, i.e. the singular Jacobian matrix of the vector field, calculated at ¼ 0
q n10 4n20 n20 : Df ðn10 ; n20 ; 0Þ ¼ 2n20 2ðq n10 4n20 Þ Consequently, the same procedure used for the -order approximation can be applied to determine the outer and inner terms of order k P 2. 4. Fast and slow variables, compressed solution As seen in the previous section, the source of singularity, which forces consideration of the step n þ 1 and a relevant compatibility condition in order to determine the nth order solution, is the operator B, as given by (8). Following the detailed procedure devised in [7], we can easily realize that the null space of B, N ðBÞ, is made up by pairs n ¼ ðn1 ; n2 Þ in the phase space such that n2 ¼ 0. Analogously, the range of B, RðBÞ, is constituted by the pairs n ¼ ðn1 ; n2 Þ with n2 ¼ 2n1 . In order to single out the hydrodynamic variable and to proceed to the compressed (Chapman–Enskog) asymptotic expansion, it is crucial the determination of the projection operator P such that P ðnÞ 2 N ðBÞ and P 2 ðnÞ ¼ P ðnÞ for any n in the phase space, and P ðnÞ ¼ 0 8n 2 RðBÞ. In the present case, it is easily seen by inspection that P is determined as 1 P1 ðn1 ; n2 Þ ¼ n1 þ n2 ; 2 ð57Þ P2 ðn1 ; n2 Þ ¼ 0: The main feature of the operator P is that its application to Eq. (4) eliminates the dominant part of the collision operator, by virtue of the identity P ðBðnÞÞ ¼ 0
ð58Þ
for any n in the phase space. Since P and time derivative obviously commute, this leaves an equation with slow evolution for the quantity P ðnÞ, the hydrodynamic variable of [7]. The complementary projection operator Q is then 1 Q1 ðn1 ; n2 Þ ¼ n2 ; 2 ð59Þ Q2 ðn1 ; n2 Þ ¼ n2 with n ¼ P ðnÞ þ QðnÞ, P ðnÞ 2 N ðBÞ and QðnÞ 2 RðBÞ. Such decomposition introduces spontaneously new dependent variables 1 x 1 ¼ n1 þ n2 ; 2 ð60Þ x 2 ¼ n2 which are respectively the slow and the fast variable in the jargon of singular perturbation theory [8].
M. Galli et al. / Appl. Math. Comput. 146 (2003) 509–531
The set (15), recast in terms of the variables (60), reads as 8 2 1 2 1 1 2 > < x_ 1 ¼ g1 ðx1 ; x2 Þ ¼ 3x1 þ 3x1 x2 4x2 3qx1 2qx2 þ 2q ; _x2 ¼ g2 ðx1 ; x2 ; Þ ¼ 4x21 þ 2ð1 þ 2Þx1 x2 þ ð3 þ Þx22 > : 4qx1 2ð1 þ Þqx2 þ q2 ;
523
ð61Þ
with phase space given by the triangle defined by x2 P 0, x1 P ð1=2Þx2 , 2x1 þ x2 6 q (the unphysical equilibrium corresponds to the vertex ðq=4; q=2Þ). It is only matter of some algebra to show that the unique fixed point in the admissible phase space corresponds to the point (10), and that the problem is now regular singular, since the limiting equations for ¼ 0 have the unique pffiffiffi fixed point of physical meaning x1 ¼ qð3 3Þ=6, x2 ¼ 0, which corresponds to the limiting value of (10) when ! 0. We can apply then the standard composite asymptotic expansion x1 ðt; Þ ¼ xo1 ðt; Þ þ xi1 ðs; Þ; x2 ðt; Þ ¼ xo2 ðt; Þ þ xi2 ðs; Þ
ð62Þ
and the initial conditions x1 ð0Þ ¼ b þ ðc=2Þ, x2 ð0Þ ¼ c. Outer and inner components are expanded in power series of X X xoj ¼ k xjk ðtÞ; xij ¼ k Xjk ðsÞ; j ¼ 1; 2 ð63Þ k
k
with initial values c x10 ð0Þ þ X10 ð0Þ ¼ b þ ; x20 ð0Þ þ X20 ð0Þ ¼ c; 2 x2k ð0Þ þ X2k ð0Þ ¼ 0; k ¼ 1; 2; . . . x1k ð0Þ þ X1k ð0Þ ¼ 0;
ð64Þ
and with initial layer corrections Xij ðsÞ bound to vanish when s ! 1. We will not give all detailed manipulations. The main variation with respect to the previous section is that now each step determines the unknowns relevant to that step, reducing the amount of standard but lengthy algebra. Another advantage is that for the present approach several results about asymptoticity of the series and error estimates are well established in [8,11], providing an adequate mathematical setting to the present formal developments. Of course, after truncation either to the zero order or to the first order, upon going back to the density variables n1 and n2 , the two analytical approximate solutions coincide. It suffices to mention here that in the solution of the set (61) the zero-order outer and inner components have to satisfy respectively dx10 1 1 q2 ¼ 3x210 þ 3x10 x20 x220 3qx20 qx20 þ ; 4 2 dt 2
ð65Þ
0 ¼ 2x10 x20 þ 3x220 2qx20 ;
ð66Þ
524
M. Galli et al. / Appl. Math. Comput. 146 (2003) 509–531
dX10 ¼ 0; ds
ð67Þ
dX20 2 ¼ 2x10 ð0ÞX20 þ 2X10 X20 þ 3X20 2qX20 : ds
ð68Þ
From (66) one has x20 ðtÞ ¼ 0 and therefore also n20 ðtÞ ¼ 0, n10 ¼ x10 . Since X10 must vanish asymptotically, from (67) it is obtained X10 ðsÞ ¼ 0. Thus we have N10 ¼ N20 =2, already obtained in (35), and moreover the initial condition c c ð69Þ x10 ð0Þ ¼ X10 ð0Þ þ b þ ¼ b þ 2 2 for Eq. (65), that after suppression of x20 reads as dx10 q2 ¼ 3x210 3qx10 þ ; dt 2
ð70Þ
and is coincident with (29) since x10 ¼ n10 . Finally, by replacing x10 ð0Þ and suppressing X10 in (68) one obtains dX20 2 ¼ k1 X20 þ 3X20 ; ds
ð71Þ
whose solution satisfying the initial condition X20 ð0Þ ¼ c is the same as N20 ðsÞ given by (37). In analogous fashion, the -order terms in the expansions (63) can be determined by solving the following sets of equations, that result by equating in system (61) the terms of the same order dx11 q ¼ 3ð2x10 qÞx11 þ 3x10 ð72Þ x21 ; 2 dt 0 ¼ 2ðx10 qÞx21 þ 4x210 4qx10 þ q2 ; a dX11 c 1 2 ¼ 2b X20 X20 ; 2 2 4 ds dX21 e 0 ðsÞ ¼ 2X20 X11 þ ð6X20 kÞX21 þ B ds
ð73Þ ð74Þ ð75Þ
with e 0 ðsÞ ¼ ½2_x10 ð0Þs þ 2x11 ð0Þ þ 6x21 ð0Þ 2ða þ cÞ X20 þ X 2 : B 20 Eq. (73) gives x21 ðtÞ ¼ n21 ðtÞ already obtained in (28), and the solution of (74) satisfying the asymptoticity condition requested to the initial layer corrections is X11 ðsÞ ¼ IðsÞ Ið1Þ, where IðÞ are the same integrals shown in (47) and (48). Substitution of these results into (72) and (75) and calculation of the respective initial conditions yield the following two linear first-order ordinary differential equations:
M. Galli et al. / Appl. Math. Comput. 146 (2003) 509–531
525
dx11 ð6x10 qÞð2x10 qÞ2 ¼ 3ð2x10 qÞx11 þ ; dt 4ðq x10 Þ x11 ð0Þ ¼ X11 ð0Þ ¼ Ið1Þ;
ð76Þ
dX21 e 0 ðsÞ; ¼ ðX20 kÞX21 þ 2X20 ½IðsÞ Ið1Þ þ B ds ðq kÞ2 X21 ð0Þ ¼ x21 ð0Þ ¼ k
ð77Þ
whose solution yields the remaining unknown terms at the stage. Thus, by restoring the original density variables, the first-order approximate solution of (15) is determined as t 1 t 1 t 1 ð1Þ n1 ðt; Þ ’ x10 ðtÞ X20 þ x11 ðtÞ x21 ðtÞ þ X11 X21 ; 2 2 2 ð78Þ h t i ð1Þ þ x21 ðtÞ þ X21 : ð79Þ n2 ðt; Þ ’ X20 A rather tedious algebra allows to verify that it is indeed equivalent to the Hilbert approximation of the same order. As mentioned in Section 1, once the problem is reduced to the form (61) and the slow (hydrodynamic) part has been singled out, an alternative algorithm is provided by the compressed asymptotic expansion [7], by leaving the outer component of the slow variable unexpanded. This allows to determine xo1 by solving a unique -dependent differential equation, which hopefully retains more information with respect to a truncation after few powers of . We shall follow this procedure according to the recipe proposed in [12]. We write t
xo1 ðtÞ ¼ x1 ðtÞ; X k uk ðx1 Þ xo2 ðtÞ ¼
ð80Þ
k
and try to determine the unknown functions x1 and uk by requiring the outer solution to satisfy the set (61). We have ! X X dx1 ¼ g1 x1 ; k uk ðx1 Þ ¼ k g1;k ; ð81Þ dt k k
X du dx1 X du X X dxo2 ¼ ¼ k k k k j g1;j ¼ k g2;k dt d x dt d x 1 1 j k k k
with 1 1 1 g1;0 ¼ g1 ðx1 ; u0 ðx1 ÞÞ ¼ 3x21 þ 3x1 u0 u20 3qx1 qu0 þ q2 ; 4 2 2
ð82Þ
526
M. Galli et al. / Appl. Math. Comput. 146 (2003) 509–531
g2;0 ¼ g2 ðx1 u0 ðx1 Þ; 0Þ ¼ 2x1 u0 þ 3u20 2qu0 ;
dg1 u q g1;1 ¼ ¼ u1 ðx1 Þ 3x1 0 ; 2 d ¼0 2
dg2 g2;1 ¼ d ¼0 ¼ u1 ðx1 Þð2x1 þ 6u0 2qÞ þ 4x21 þ 4x1 u0 þ u20 2qu0 4qx1 þ q2 : By equating in (82) the coefficients of 0 and it is obtained the system 8 < g2;0 ¼ ½2x1 þ 3u0 2q u0 ¼ 0; ð83Þ du0 : g2;1 ¼ g1;0 dx1 which has the solution 2
u0 ¼ 0;
u1 ðx1 Þ ¼
ð2x1 qÞ ; 2ðx1 qÞ
ð84Þ
where x1 is determined by solving dx1 ¼ g1 ðx1 ; u1 ðx1 ÞÞ dt " # q2 q2 ð2x1 qÞ4 2 ¼ 3x1 ðx1 qÞ þ þ ðq 6x1 Þ x1 þ : 2 2 4ðx1 qÞ 16ðq x1 Þ ð85Þ The initial condition for (85) and the ones of the inner terms X1k , X2k are now related, up to the -order approximation, as follows: c x1 ð0Þ þ X10 ð0Þ þ X11 ð0Þ ¼ b þ ; 2
ð86Þ
u0 ½x1;0 ð0Þ þ X20 ð0Þ ¼ c;
ð87Þ
u1 ½x1;0 ð0Þ þ X21 ð0Þ ¼ 0;
ð88Þ
where x1;0 ð0Þ ¼ b þ c=2 is the initial value of x1 ðtÞ in the limit ! 0. They are calculated, as usual, by matching the outer solution with the inner one, which must satisfy the following system: 8 i dx > > < 1 ¼ g1 ½x1 þ xi1 ; u1 ðx1 Þ þ xi2 g1 ½x1 ; u1 ðx1 Þ ; ds ð89Þ i > dx > 2 i i : ¼ g2 ½x1 þ x1 ; u1 ðx2 Þ þ x2 ; g2 ½x1 ; u1 ðx1 Þ; : ds
M. Galli et al. / Appl. Math. Comput. 146 (2003) 509–531
527
By equating the terms of order zero and one with respect to , one has the same set provided by the previous approach, e.g. Eqs. (67), (68) and (74), (75) for Xi0 and Xi1 respectively, where the nonhomogeneous term in (75) is turned to e 0 ðsÞ ¼ ½2x_ 1 ð0Þs þ 6u1 ðx1 ð0ÞÞ 2ða þ cÞ X20 þ X 2 : B 0 20 It follows that the initial condition for (85) must be c c x1 ð0Þ ¼ b þ X10 ð0Þ X11 ð0Þ ¼ b þ þ Ið1Þ 2 2 and it can now be used to find the unexpanded outer component of the slow variable, by integration of the differential equation (85). In conclusion, in terms of the original density variables n1 , n2 the -order approximation achieved by using the compressed expansions algorithm is t 1 t 1 t 1 ðcÞ n1 ðt; Þ ’ x1 ðtÞ X20 þ u1 ðx1 Þ þ X11 X21 ; 2 2 2 ð90Þ ðcÞ
n2 ðt; Þ ’ X20
h t i þ u1 ðx1 Þ þ X21 :
t
ð91Þ
We are not reporting on the zero-order approximation since it is not difficult to check that the unexpanded part x1 ðtÞ would coincide with the unknown x10 of the uncompressed procedure, and consequently we would get the same results from the two algorithms.
5. Numerical results The results of the approximate asymptotic solutions derived in the previous Sections are now illustrated and discussed. We consider the dynamical system (15) together with initial conditions (16) and we first compute its ‘‘exact’’ (numerical) solution by means of an accurate ODE solver; the density n0 comes then from the first integral (3). In our numerical examples we start with the following reference numerical values ¼ 0:01;
a ¼ 0;
b ¼ 1;
c ¼ 4:
They have to be considered as dimensionless (for illustrative purposes only), and corresponding to arbitrary scales. In Fig. 1(a) we compare the temporal evolutions of n0 , n1 , n2 obtained by the exact solution (solid line), by the zero-order Hilbert expansion (dotted line), and by the first-order Hilbert expansion (dashed-dotted line). The initial layer clearly shows up, and it is evident that the first-order Hilbert approximation better reproduces the exact behaviour (the very slight differences can be hardly
528
M. Galli et al. / Appl. Math. Comput. 146 (2003) 509–531
6
6 n
n
0
0
5
5
4
4
3
3 n1
n1
2
2
1
1 n
0
n
2
0
0.02
0.04
0.06
0.08
(a)
0.1
0
2
0
0.02
0.04
0.06
0.08
0.1
(b)
Fig. 1. Reference case ¼ 0:01. (a) Exact solution (solid line), zero-order (dotted line) and firstorder (dashed-dotted line) Hilbert expansion; (b) stationary-state approximation (dotted line), outer expansion (dashed-dotted line) and exact solution (solid line).
recognized in the graph). We recall that, as already seen, the zero-order Hilbert expansion coincides analytically with the zero-order compressed expansion and with the zero-order solution of the fast–slow variables method, as well as, at first order, the Hilbert approximation and the fast–slow variables method give exactly the same analytical results. The widely used stationary-state approximation has been compared with the outer solution of the first-order Hilbert expansion in Fig. 1(b) where the temporal evolution of the densities ni is shown. The initial data are the same as for the outer solution in the Hilbert expansion, namely n1 ð0Þ ¼ n10 ð0Þ þ n11 ð0Þ; where n10 ð0Þ and n11 ð0Þ are given by (34) and (49). Similar results could be obtained by choosing initial data in such a way that the least square error between the exact and the stationary-state solutions is minimal. The solid line represents again the reference exact numerical solution, whereas the stationarystate approximation and the outer solution, which overlap almost entirely, are given by the dotted and the dashed-dotted lines respectively. We can notice that, apart from the initial layer of thickness OðÞ, the stationary-state approximation provides an accurate description of the evolution in the bulk region. In Fig. 2(a) we report the results as in Fig. 1(a) obtained by changing the initial data. A different trend with respect to the reference case can be observed in the evolution of the three number densities, and a wide range of behaviours
M. Galli et al. / Appl. Math. Comput. 146 (2003) 509–531 6
529
6 n0
n0
5
5
4
4
3
3
n
n1
1
2
2
1
1 n2
n 0
2
0
0.02
0.04
0.06 (a)
0.08
0.1
0
0
0.02
0.04
0.06
0.08
0.1
(b)
Fig. 2. (a) Exact solution and Hilbert expansions obtained for different initial data (a ¼ 0, b ¼ 2:5, c ¼ 2:5) and the same ¼ 0:01. (b) Exact solution (solid line), zero-order (dotted line) and firstorder (dashed-dotted line) Hilbert expansion obtained for ¼ 0:05 and initial data as in Fig. 1.
could be obtained by playing with initial data. Of course, the initial layer effect tends to disappear when the initial condition approaches the quasi-steady-state condition (9), as already pointed out in [6]. Indeed, it was already observed in (36) and (37) that the zero-order inner solution vanishes for c ¼ 0. The effect of varying is depicted in Fig. 2(b) where we plot the results obtained for ¼ 0:05 when we start with the reference initial data. As expected, a larger initial layer shows up, and the zero- and even first-order Hilbert approximations lead to worse results than in case ¼ 0:01. The first-order compressed expansion gives rise to solutions which differs, in principle, from the ones obtained by the first-order Hilbert approximation, but the discrepancies are not quantitatively appreciable in the graphs ðt; ni Þ. In Fig. 3(a) we then show the temporal trend (in logarithmic scale) of the error vectors, ex defined as the pointwise euclidean norm of the vector (n0 nex 0 , n1 n1 , ex ex n2 n2 ), where ni stands for the exact numerical values. We compare there the errors for the two values of used above versus time. The solid lines are relevant to the error in the first-order Hilbert expansion, whereas the dotted lines represent the error in the first-order compressed expansion. As we can see, the two methods do not show significant discrepancies in the results, although they come from two different asymptotic strategies. The compressed algorithm enhances a bit the accuracy in the bulk region. Fig. 3(b) shows instead the behaviour of the errors (with respect to the exact solution again) in the stationary-state approximation (solid line) and in the first-order Hilbert expansion (dotted line), with exclusion of the very initial region, where the former exhibits
530
10
M. Galli et al. / Appl. Math. Comput. 146 (2003) 509–531 0
4
3
3.5
eps=0.05 10
x 10
2
3 2.5
10
4
2 1.5
eps=0.01 10
6
1 0.5
10
8
0
0.2
0.4
0.6
0.8
0
0.005
0.2
(a)
0.4
0.6
0.8
(b)
Fig. 3. (a) Comparison of the error in the first-order Hilbert expansion (solid line) and in the firstorder compressed expansion (dotted line) obtained for ¼ 0:01 and 0.05 respectively (logarithmic scale). (b) Error in the stationary-state approximation (solid line) and in the Hilbert first-order expansion (dotted line) (case ¼ 0:01).
very high deviations. It is worth noticing that in the steady-state approximation the error tends to vanish for t large enough, since no truncation occurs in the fixed point of the problem with respect to the exact one. To summarize the results, we report finally in Table 1 the values of the density n1 computed at time instant t ¼ 0:01, t ¼ 0:001 by using both the ODE solver and the different asymptotic approximations described above. The values marked by are not significant since they are relevant to the stationarystate approximation out of its region of validity, namely in the initial layer.
Table 1 ð0Þ ð1Þ Computed values of density n1 obtained by zero- (n1 ) and first-order (n1 ) Hilbert expansion, firstðcÞ ðsÞ order compressed method (n1 ) and stationary-state approximation (n1 ) at two different time instants t ¼ 0:01 and t ¼ 0:001 ð0Þ
ð1Þ
ðcÞ
ðsÞ
n1
n1
n1
n1
n1
t ¼ 0:01
0.01 0.005 0.001
2.8912 2.8844 2.8790
3.3924 3.3834 3.3762
2.8961 2.8868 2.8795
2.8962 2.8869 2.8795
2.8915 2.8846 2.8790
t ¼ 0:001
0.01 0.005 0.001
1.6056 2.3716 2.9886
2.4688 3.1602 3.5526
1.6053 2.3720 2.9892
1.6054 2.3721 2.9892
3.0080 2.9950 2.9987
M. Galli et al. / Appl. Math. Comput. 146 (2003) 509–531
531
Acknowledgements This work has been performed in the framework of the activities sponsored by MIUR (Projects ‘‘Problemi matematici non lineari di propagazione e stabilit a nei modelli del continuo’’ for RR and ‘‘Problemi matematici delle teorie cinetiche’’ for the other three authors), by CNR, by GNFM, by University of Parma (Italy), and by the European TMR-Network ‘‘Asymptotic Methods in Kinetic Theory’’.
References [1] I. Prigogine, E. Xhrouet, On the perturbation of Maxwell distribution by chemical reactions in gases, Physica 15 (1949) 913–932. [2] A. Rossani, G. Spiga, A note on the kinetic theory of chemically reacting gases, Physica A 272 (1999) 563–573. [3] M. Groppi, G. Spiga, Kinetic approach to chemical reactions and inelastic transitions in a rarefied gas, J. Math. Chem. 26 (1999) 197–219. [4] J. Mc Lennan, Boltzmann equation for a dissociating gas, J. Stat. Phys. 57 (1989) 887–905. [5] Y. Yoshizawa, Wave structures of a chemically reacting gas by the kinetic theory of gases, in: J.L. Potter (Ed.), Rarefied Gas Dynamics, AIAA, New York, 1975, pp. 501–517. [6] M. Groppi, A. Rossani, G. Spiga, Kinetic theory of a diatomic gas with reactions of dissociation and recombination through a transition state, J. Phys. A: Math. Gen. 33 (2000) 8819–8833. [7] J. Mika, J. Banasiak, Singularly Perturbed Evolution Equations with Application to Kinetic Theory, World Scientific, Singapore, 1995. [8] R.E. OÕMalley Jr., Singular Perturbation Methods for Ordinary Differential Equations, Springer, New York, 1991. [9] C. Cercignani, The Boltzmann Equation and its Applications, Springer, New York, 1988. [10] A.W. Adamson, A Textbook of Physical Chemistry, Academic, New York, 1979. [11] D.R. Smith, Singular Perturbation Theory, Cambridge University Press, Cambridge, 1985. [12] J.R. Mika, A. Palczewski, Asymptotic analysis of singularly perturbed systems of ordinary differential equations, Comput. Math. Appl. 21 (1991) 13–32. [13] J.D. Cole, Perturbation Methods in Applied Mathematics, Blaisdell, Waltham, MA, 1968. [14] W. Eckhaus, Matched Asymptotic Expansions and Singular Perturbations, North Holland, Amsterdam, 1973. [15] A.H. Nayfeh, Perturbation Methods, Wiley, New York, 1993.