- Email: [email protected]
Bifurcation phenomena in a detonation problem
Applied Numerical Mathematics 9 (1992) 427-445 North-Holland
427
Wen Zhang and Ian Gladwell Mathematics Department, Southern Methodist University, Dallas, TX 75275, USA Abstract Zhang, 6~. and 1. Gladwell, Bifurcation phenomena in a detonation problem, Applied Numerical Mathematics 9 (1992) 427-445. In [5,9] a system of paraneterized differential algebraic equations (DAEs), modeling a detonation problem, is posed and its solution computeJ. Due to computational difficulties, these numerical results arc incomplete. A bifurcation curve stops abruptly during continuation in the parameter D, the detonation velocity. In this paper we transform the problem into a parameterized ordinary differential equation (ODE) boundary value problem (BVP), exhibiting a singularity on the free boundary. We present an effective strategy to deal with the singularity. Using this strategy we successfully extend the bifurcation curve. Also we find a new bifurcation phenomenon with a continuous region of so!utions.
I. The detonation problem 1.1. Introduction We investigate a ZND [7] model of an unconfined small divergent detonation along a cylinder. The problem is illustrated in Fig. 1 which is a modified version of the corresponding figure in [4].
Unreacted Explosive
Reaction Products
Fig. 1. Detonation along a cylinder.
016%9274/92/$05.00
0 1992 - Elsevier Science Publishers B.V. All rights reserved
W. Zhang, I. Gladwell / Bifurcation phenomena
428
coordinate system is attached to the shock front which moves at a constan: velocity D into the unreacted explosive. We investigate the reaction zone based on one-dimensional flow theory in Lagrangian coordinates. The position of the CJ plane depends on the solution. Kirby and Leiper [9] gave the physical background and the general formulation of this problem. Braithwaite et al. [4] further described it in mathematical detail. The system governing the reaction zone is an autonomous DAE BVP with time t as the independent variable. It has nine dependent variables, (u, A, A, w, X, p, v, p, E), corresponding to five ODES and four algebraic equatiins (AEs). There are four undetermined constants, (D, d, TcJ, xcJ), and eight boundary conditions (BCs) corresponomg to five first-order ODES. Hence one constant can be considered a free parameter and the remainder ar 9 unknowns. It is of physical interest to determine the relationship between the diameter d and the detonation velocity D or the radius of curvature of the shock front R, and D [7]. The numerical results in [4,5,9] show that d, T,, and xcJ are single-valued functions of D but not vice-versa. Hence D is treated as a free parameter for computational convenience. A further discussion on the roles of these constants is given in Sections 1.2 and 2.1. The physical meanings of the variables, the unknown constants and the parameter are: The
l4
-
A
-
A
-
0
-
x
-
P c
-
P E
-
D d TCJ 'CJ -
particle velocity, area of a stream tube of fluid which had unit area degree of chemical reaction (A = 0 for no reaction, radial velocity, distance from the shock front, pressure, sp\=cific volume, fluid density, specific internal energy, detonation velocity, diameter of the cylinder of explosive, time needed to reach the CJ plane from the shock the (CJ) distance between the shock front and the
at the shock front, A = 1 for complete reaction),
front, CJ plane.
The ODES are puu’
+p’
= 0,
(l.Ia)
A’= 2Ao, A’= (I-
A)F(h,
0’ = -u’/R(x, X4=u.
(l.lb) p), +Jv
(l.lc) d),
(1
.ld)
(Lie)
(1. la) is the momentum conservation equation. (1. lb) describes flow divergence and is derived geometrically [4,7,9]. (1 .lc) describes the chemical reaction rate and (Lid) defines the variation of the radial divergence.
W. Zhang, I. Gladwell / By.furcation phenomena
429
The AEs are E =H(p,
v, A, D),
(1.2a)
E -t pu + ;ti2 = +D’,
(1.2b)
Apu - p,D = 0,
(L2c)
v = l/p.
(1.2d)
(1.2a) defines the energy. (1.2b) is the integral of the enthalpy conservation equation combined with the momentum equation, which is Bernoulli’s equation. (1.2~) states mass conservation. The BCs are
u = u,(D),
(1.3a)
A = 1,
(1.3b)
A =o,
(1.3c)
0 =
[D
- ~,ja)l/~,(+J~
(1.3d)
d),
A:=0
(1.3e)
at t = 0 (the shock front) and
u =c(P,
0, A, D),
CL) = (T(p, v, A, c, D)A’/2, x =‘CJ
(1st CJ condition), (2nd CJ condition),
(1.4a) (1.4b) (1.4c)
at t = (the CJ plane). c is the speed of sound and U- is the thermicity coefficient of the fluid. The value is determined as the time when both (1.4a) and (1.4b) are satisfied. (1.4a,b) are the Chapman-Jouguet (CJ) conditions. (1.4a) determines where particle velocity reaches sound speed and (1.4b) balances the release of energy by the chemical reaction and its withdrawal by the radial flow [4,7]. There are physical constraints which must hold over the whole integration interval: TcJ
TcJ
O.
v, A, D),
(1 Sa)
O
(1.5b)
A 2 1,
(1%)
P, u, p > 0.
(1.5d)
and strict inequality holds The equality in (1.5a) can only be achieved at the CJ plane (t = for t E [O, po, the explosive porous density, is chosen as 1.08 as in [4,5]. R(x, p, V, A, D), c( p, v, A, D), a( p, u, A, c, D), F(A, p), us(D) and Rs(xCJ, d) are modeling functions. R is the radius of curvature of an isobar at a distance x downstream of the shock front. M is the specific internal energy and F describes the pressure dependence. R, has the form bee [2,9,101) TcJ)
TcJ)
[9].
xcJ7
R,(X,,,
d)=(d-&,)/a.
d),
M(
(1 .6)
W Zhang, I. Gladwell / Bifurcation phenomena
430
Kirby, Leiper and Hackett [lo] determined a! = 0.5 and p = 1.2 by experiment. We use A
1-A H(p, v, A, D) =PO
c(p,
v, A, D) =
I
y,‘(D) - 1+
PV(V,
YRW - l Phrl(V)[Y(%
A, D) +
[Y&J)
\
a(p,
v, A, c, D) =
y(4 A, D! - 1 c2
Q +-E,(D)
i
. (18)
- AQ + (1 -A)&(D), 4 -
D) - 11
II2
1’2
. (19)
’ i
1 -PO
Y&J
- 1 - rm
- 1
(1.10) and for the emulsive explosive considered here F(AY P) = Q(A)PP(
(1.1 la)
P)/q-j
where if p>,&,
(P-Pm PP(. P) =
:P
3p 4p,
3 ’
(l.llb) otherwise.
ii I Q, the heat of reaction, and pH, the hot spot critical pressure, are experimentally determined constants. yU(D), the index of the unreacted ingredients; E,(D), the turbulent stored energy; y(v, A, D) and u,(D) depends implicitly on D. As shown in [5] they can be determined by iteration. Y&J) and aH( A) are explicit functions. A detailed description and explanation of these functions is given in [5,14]. In [3,4,13] the problem is investigated with simplified function definitions for R( X, HCp, u, A, DJ,c(p, u, A, D), d p, v, A, c, D), F(A, p) and u&D). Since the relationship between the unknown constants and the parameter is of practical importance, our view is focused on them. represents the solution variable x(t) at t = T,, and can be considered a representative solution variable. The projection of a solution to the DAE BVP (1 .l)-( 1.l 1) in parameter space, such as the (l/d, D) plane, is defined as a solution in that space. Kirby and Leiper show the (l/d, D) diagram in [9]. Their computed curve stops abruptly at D = 0.40: where D, is the ideal detonation velocity. They did not succeed in computing the curve for lower values of D. Braithwaite et al. [5] computed the solution with a similar result. We continue the investigation of this problem below. +J,
d),
.xcJ
1.2. Reformulations of the detonation problem and discussions on the singularity Braithwaite et al. [4,5] computed the solution by solving for the nine unknown variables in the DAE BVP (l.l)--(1 S). The index [6] of this system with respect to the nine unknown variables is shown in [4,5] to be one everywhere except at the boundary t = where it is two. We observe that the size of the system can be reduced as in 1131.Since pO, D > 0, (1.2~) and ( 1.2d) imply &J
V =-
AL1 (I
PqD’
.12)
W Zhang, I. Gladwell / Bifurcation phenomena
431
Substituting (1.2a) into (1.2b), only one AE is needed and the number of unknown variables is reduced from nine to six (u, A, A, w, X, p). p, o, and E, can be obtained directly from them. Therefore the system is reduced to (l.lb)--(l.le), p&4’ +p’A = 0,
(1.13)
and
-
= --
(1.14)
with BCs (1.3) at t = 0, and, at t = TcJ, (1.4c), and (1.15a)
U
Au
------,A
(1.15b)
P, D
The functions H, c and CTare always defined under the rule that the adiabatic equation, y’= -pu2(aA’ - 20)/(1
- u”/c”),
(1.16)
must be satisfied, [9]. In other words, if we differentiate AE (1.14) and combine with the other ODES, we can always obtain (1.16). Thus the index of the reduced DAE system is one everywhere except where the denominator of (1.16) is zero. The denominator only vanishes at t = 7’&, [9], and it can be shown, as in [4,5], that there the index is two. Hence the reduced system has the same index as the unreduced system. One way to solve the reduced DAE BVP is to replace AE (1.14) by (1.16) and solve the enlarged ODE BVP for the six unknown variables (u, A, A, W, x, p) [6]. Solving this way, (1.14) may be satisfied inaccurate!y. Since (1.14) reflects energy conservation, its accuracy is of physical importance. We observe from (1.8) that p in (1.14) can be determined explicitly. Using (1.8), we write (1.14) as
p=P(A,u,A,D,P(,)=
Dpo[ D2/2 - u2/2 + AQ - (I-A)&] . l-h A +1 \ Au --l-----Y” - 1 YR- 1 J
Substituting (1.16) into (1.13) with l/u for ld =
uc2(crA’- 2w)/(c2
- u2).
p
(1.17)
and (1.12) for u, we have (1.18)
Replacing p by (1.17) wherever it appears, equations (1.18) and (I.lb)-(l.le) with BCs (1.3), (1.15), and (1.4~) form a BVP of five ODES with five unknown variables (II, A, A, w, x). (1.14) is satisfied throughout due to the substitution for p. Equation (1.18) exhibits singularity whenever u = c though we do not perceive it in the original form (l.l)-(1.4). We believe this singularity is the cause of the difficulties in the computations reported in [4,5]. There are two ways to deal with it. One way is to bring it out in the ODES, as above. An alternative way, as suggested by Ascher [l], is to keep it inside the AE (1.14) then to form a new set of variables and embed (1.14) in the ODES, [14]. Though we do
W. Zhang, I. Gladweld / Bifmation phenomena
432
not observe the singularity in the new ODE system, this does not mean th-:I it has been removed from the problem. Instead it is converted into a solvability obstacle in the AE. This is demonstrated independently in [3]. The computation is expensive. For instance, if we use shooting methods to solve the new ODE BVP, we need to solve the embedded AE in each call to the derivative evaluation routine. When we cannot solve the AE we cannot proceed with the integration. We choose to work with the form which exposes the singularity in the ODES because it is directly related to the physical phenomena and has an important advantage in investigating the nature of the solution. We observe from (1.7) that R(x, xcJ, d) = R(x, R,(x,, d)) and xcJ, d appear only in (l.ld) and (1.3d) through R, (without taking account of (1.4~)). We introduce a new unknown corzsturzt R, to replace xcJ and d. Since xcJ =x(7’&, it is automatically determined as soon as the solution is found. Then d can be computed directly from (1.6). Thus, the undetermined constants are reduced to R,, D and T,,. The BCs are coupled through xcJ in (1.3d) which is posed at t = 0. After introducing R, (1.3d) simply becomes
40) = [D -
(1.19)
u,(D)l/Rs
and (1.4~) is no longer needed. Therefore the BCs separate. Replacing R(x, xcJ, d) by R(x, R,) in (l.ld), (1.3d) by (1.19) and discarding (1.4~)~ we finally obtain a BVP with five explicit ODES (1.18), (1 .lb-e) corresponding to five unknown variables (u, A, A, o, X) (replacing p by (1.17) wherever necessary); and seven separated BCs (1.3a)-( 1.3c), (1.19), (1.3e) and (1.15) corresponding to five first-order ODES with three unknown constants (R,, D, T,,). As before, one of the constants must be considered as a free parameter and the other two are determined by the solution of the ODE BVP. As shown in Section 2.1, this system is well-posed in a sense that for each (R,, D, T,,) there is at most one solution of the ODE BVP. Further, each solution corresponds to a unique (R,, D, Z’J. Hence the solution of the system can be represented in (R,, D, T,,) space. As there is only one free parameter, the solution expected in CR,, D, T,,) space is a curve. Physically either D or R, can be viewed as the free parameter. Computationally selecting D is advantageous. We restate the reduced system below. The ODES are U'
= uc2(crh’ - 2w)/(c2 - u2),
(1.20a)
A’ = 2Ao,
(1.20b)
A’= (1 - h)F(A, p), 0’ =
(1.2Oc)
-u’;R(x, R,),
(1.20d)
x’=u,
(1.20e)
where P =-
Dp,[ D2/2 - u2/2 + AQ - (1 - A)E,] l-h AU
-+
Yu-
.
A -+1 1
YR-
1
(1.21)
W. Zhang, 1: Gladwell / Bifurcation phenomena
433
The BCs are I/ =
u,(D),
(L22a)
A = 1,
(1.22b)
=o,
(1.22c)
A
0 = [D-u,(D)]&
(1.226)
X=0
(1.22e)
at t = 0, and U
- A)&
p)/2,
(1st CJ condition),
(1.23a)
(2nd CJ condition),
(1.23b)
at t = T,,. The physical constraints (1.5) still hold here. Mathematically there is a singularity in (1.20a) and it occurs at the free boundary t = T,.
2. Numerical simulations and bifurcation phenomena 2.1. Computational strategy for solving the parameterized singular ODE BW There are a variety of numerical methods which might be used to solve the ODE BVP (1.20)-(1.23). We choose a shooting method for following reasons. First there are five BCs (1.22) at the left boundary, t = 0, and five explicit ODES (1.20). These constitute a well-defined initial value problem for given (R,, D). They are already in a standard form for calling an ODE solver. Since the right end of the interval is a free boundary, we may advance the integration as far as needed and T,,, the right boundary position, is determined as the point where the CJ conditions, (1.23) are satisfied. The aim of shooting is to satisfy the two CJ conditions simultaneously. For given R, and D there is at most one 7& because we always terminate the integration the first time (1.23) are satisfied, even if the integration can proceed. It is nonphysical for the first CJ condition (1.23a) to be satisfied in (0, T,,) [9]. So the solution is unique for given R, and D. Using a shooting method we have the opportunity to observe what happens even when a choice R,, D leads to no solution. With this facility we can produce a bifurcation diagram and distinguish different regions in the (R,, D) plane. This information is useful in determining the solution manifold. A simple implementation of a shooting method is to use an ODE solver to perform the integration and use a nonlinear equation solver to match the BCs. There is a choice of the direction of the integration. We choose to integrate from t = 0 to t = T”,. In addition to the reasons above, the decision was also based upon experiments in [4,5]. There a shooting method is also employed but integrating from t = TcJ to t = 0. This choice was made to avoid integrating towards the singularity but there is a difficulty in finding the solution We need to determine ~1’at the singular point t = T’,, in order to start the integration. This requires solving
W Zhang, I. Gladwell / Bifurcation phenomena
434
nonlinear equations involving values of the solution variables at t = T,, which are the initial estimates in the shooting method. Hence the selection of initial estimates must guarantee a solution for u’(T,,), otherwise the integration cannot start. After the shooting method is started,. the nonlinear equation solver adjusts the variables at TcJ internally and there is no guarantee of a solution u’( T,,). Another frequently occurring difficulty in [4,5] when computing the solution for smaller D is the integration stopping before t = 0. This is due to the ODE solver detecting a singularity in the interior of the integration interval. After this we cannot proceed to the matching phase. In order to avoid these interruptions we chose to integrate from t = 0 where the derivatives are well-defined. Since the right boundary is free, it is natural mathematically and physically to end the integration where the singularity occurs. Thus we can always apply the matching procedure after the integration, and then the shooting method continues. In this way we have succeeded in computing the solution on the lower branch in D where previous computations failed. To perform the shooting method we have used standard codes from the NAG Library [12], D02EBF for performing the ODE integration and C05NBF for solving the two nonlinear BCs (1.23). The ODE solver requires as input a fixed integration interval. So, we rescale the original integration intervai to [O,l] by setting t^= t/7& as the new independent variable. The five ODES become u’ = T&C2(CTAf/TCJ- 2w)/(c2
x’ = Tg,
(2.la) (2.lb)
A’ = 2T,,Ao, A’= T,,(I - h)F(h, w’ = -UR(x, R,),
- U2j,
P),
(2.lc) (2.ld) (2Se)
where p is computed from (1.21). There is no change in the written form of the BCs (1.221-t 1.23) except that (1.22) hold at r^= 0 and (1.23) hold at t^= 1. D02EBF uses a variable step size variable order backward differentiation formula to integrate a system of stiff ODES. Here the stiffness property is of no importance. A nonstiff solver will do the job as well. It would become appropriate if the problem was made stiff by the introduction of chemical kinetics terms in a “full” model of the detonation problem. The code has a feature that it outputs the solution at its last successful step when the integration cannot proceed to the end of the interval. This feature makes it possible to use the strategy below for continuing shooting after encountering the singularity. We proceed as follows. Since there are three undetermined constants and two matching equations, we fix one constant as a parameter and search for suitable values for the other two until the two matching equations are satisfied, that is a solution to the ODE BVP is found. We rule out fixing T,, because it makes no sense to fix a free boundary and can only reduce the chance of success in searching for a solution. We choose to fix D because it saves computational work compared to the alternatives. We observe that u,(D) in the initial data, and y,(D) and E,(D) in the function definitions are obtained via solution of a nonlinear equation [5]. If we fix R, we have to solve the extra nonlinear equation each time we change D inside the iterationThere is no such requirement when R, changes. If D is fixed we can compute these quantities before starting shooting and do not need to recompute them,
W. Zhang, I. Gladwell / Bifurcation phenomena
435
We call DKEBF to integrate equations (1.20) with the initial conditions (1.22). There are two possible outcomes. One is that the integration is advanced to t^= 1 without detecting any singularity. In this case we simply evaluate the CJ conditions (1.23). If they are satisfied to the given tolerance, we quit the shooting procedure since the solution is found. Then d is computed from (1.6) where xcJ, the solution x(l) at ? = 1, is available from D02EBF and R, is available from COSNBF. If the CJ conditions are not satisfied, COSNBF provides a new set of estimates for R, and QJ. The other possible outcome is that the integration stops at a point before t^= 1 due to encountering the singularity, that is the first CJ condition (1.23a) is satisfied but the second may not be. The code outputs the result at a point immediately preceding the singularity. We still check the residues of the two CJ conditions using output from D02EBF since it is a reliable result close to the singularity. If the CJ conditions are satisfied within the given tolerance, the solution is considered found. It does not matter that the integration does not proceed to t^= 1 since we are dealing with a free boundary. Let the value of Tcr output from COSNBF be 7’&. Because the integration ends before t^= 1, the length of the integration interval in the original variable t should be T’j = T&i. d can be computed in the same way as before. If the two CJ conditions are not satisfied, we begin the cycle again. For the above strategy the singularity is not a barrier to finding a solution. Instead we use it to prevent unnecessary work. When the estimate of T, from COSNBF is too large, we adjust it to the correct value. So the integration not only finds the solution variables but also finds the correct value of one of the unknown constants, T,,. This is essentially due to integrating towards a free boundary. Clearly the capability of accurately computing the solution “at” the singular point is crucial to the success of this strategy. We have computed the solutions of the problem using the above strategy on an IBM 3081D. All the computations are done in DOUBLE PRECISION with a mixed error tolerance 1.5 x lo-l2 assigned in both D02EBF and COSNBF. We assign a third tolerance 5 X 10m8 on the residues of the two CJ conditions (1.23) which are checked after each call to D02EBF. If they are satisfied within the tolerance, we quit the iteration and claim a solution is found. The reason for using such a stringent tolerance in the call to D02EBF is to limit the error in the integration procedure as much as possible in order to track the solution closely and to prevent violation of the physical constraints. The existence of a solution to the singular BVP depends on whether the two CJ conditions (1.23) can be satisfied simultaneously. Depending on the sensitivity of the CJ conditions, a “large” error in the integration may prevent US in finding a solution as we have seen in a related problem [13]. Also when using a crude integration tolerance, the residues of the CJ conditions become large. Usually we cannot expect the CJ conditions to be satisfied more accurately than the integration result. So we use a stringent tolerance in the ODE solver but a less stringent tolerance in checking the CJ conditions. The value of the tolerance for COSNBF is less important. The recommended value given in the NAG manual 1121is approximately the square root of the machine precision. The value we select is smaller, in order to ensure enough iterations are taken. In the code the difference between the two successive estimated solutions, rather than the residue of the function, is checked against the given error tolerance. So there is a chance that upon exit the residue of the nonlinear equations may not be as small as anticipated. This has happened with a larger tolerance in C05NBF. Using the tolerance selected above, we obtain consistently satisfactory results from the code.
W. Zhang, I. Gladwell / Bifurcation phenomena
436
We have monitored intermediate solution results and the residues of the two CT conditions throughout the integration. We observe a reduction of several orders of magnitude in the ‘_.c;i-icethat both CJ conditions are satisfied without residues as T,, is approached. This is ev!d ambiguity. Monitoring the behavior of the solution variables throughout the integration interval will be useful below. 2.2. Numerical results and analyses We have successfully traced the solution of the original probiem for lower values of i) than in [5]. We refer to this result as the main solution curve. The values of the solution at t = TcJ are referred to as the CJ values. The results are shown in Table 11and Figs. 2 and 3. For D < 1.6 we need a small continuation step in D, a few more iterations and a more accurate guess for I?,, but, in principle, we do not encounter any difficulty computing the solution. However for D < 1.1 the successful continuation step size becomes so small that it requires too much computation to trxe D further. Based on regional properties determined from our numerical results, we empioy the bisection method to replace the nonlinear equation solver Table 1 Numerical results for the main solution curve
D
Rs
TCJ
3.5 3.4 3.3 3.2 3.1 3.0 2.9 2.8 2.7 2.6 2.5 2.4 2.3 2.2 2.1 2.0 1.9 1.8 1.7 1.6 1.5 1.4 1.3 1.2 1.1 1.0 0.5 0.2
0.12+2 0.11+2 0.11+2 0.1@+2 0.99 + 1 0.96 + 1 0.93 + 1 0.91+ 1 0.90 + 1 0.91 -I-1 0.95 + 1 O.E0+2 0.12+2 0.14+2 0.17+2 0.20 + 2 0.25 + 2 0.30 + 2 0.36+2 0.45 + 2 0.56 + 2 0.71+ 2 0.92 + 2 0.12+3 0.17+3 0.25 + 3 0.64 + 4 0.11+8
0.95 + 0.98+ O.lO+ O.ll+ O.ll+ 0.12+ 0.12+ 0.13+ 0.15+ 0.16+ 0.19+ 0.23+ 0.31+ 0.44+ 0.63+ 0.87+ 0.12+ 0.16+ 0.22+ 0.32+ 0.46+ 0.69+ O.ll+ 0.18+ 0.31+ 0.59+ 0.12+ 0.85 +
d 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 3 3 3 3 6 10
0.89+ 0.86+ 0.83+ 0.81+ 0.79+ 0.78+ 0.78+ 0.73+ 0.80+ 0.84+ 0.90+ O.lO+ 0.13+ 0.17+ 0.21+ 0.27+ 0.34+ 0.42+ 0.54+ 0.69+ 0.91+ 0.12+ 0.17+ 0.25+ 0.39+ 0.64+ 0.52+ 0.18+
kJ
1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 5 10
2.9023 2.8386 2.7748 2.7107 2.6462 2.5814 2.5161 2.4501 2.3835 2.3159 2.2471 2.1767 2.1053 2.034 1 1.9637 1.8941 1.8260 1.7583 1.6912 1.6249 1.5593 1.4946 1.4309 1.3684 1.3073 1.2480 0.99459 0.90386
A CJ 1.2177 I.2318 1.2467 1.2626 1.2795 1.2975 1.3167 1.3370 1.3582 1.3803 1.4053 1.4375 1.4827 1.5426 1.6166 1.7839 1.8005 1.9419 2.05 14 2.1641 2.3180 2.4956 2.7033 2.9507 3.2509 3.6230 8.5308 65.423
kJ
OCJ
XC.1
PCJ
0.59481 0.56918 0.54326 0.51709 0.49068 0.46398 0.43690 0.4092 1 0.38048 0.35043 0.3’978 0.29058 0.26732 0.25396 0.24941 0.27326 0.25583 0.26834 0.27076 0.27224 0.27763 0.28304 0.28858 0.29449 0.30113 0.30903 0.41382 0.70323
0.10265 0.10869 0.11473 0.12066 0.12634 0.13154 0.13593 0.13899 0.13996 0.13765 0.13075 0 11866 0.10393 0.09154 0.08306 0.08077 0.073 1I 0.06928 0.0642 1 0.05923 0.05409 0.04877 0.04333 0.03784 0.03238 0.02703 0.00504 0.74 - 4
0.23 + 1 0.23 + 1 0.24 + 1 0.24 + 1 0.24 + 1 0.25 + 1 0.26 + 1 0.27 + 1 0.29 + i 0.32 + 1 0.36 + 1 0.43 + 1 0.57 + 1 0.78 + 1 0.11+2 0.14+2 0.18+2 0.23 + 2 0.30 + 2 0.39 + 2 0.53 + 2 0.73 + 2 0.11+3 0.16+3 0.25 + 3 0.43 + 3 0.41+ 5 0.15+ 10
0.26 + 1 0.24 + 1 0.22 + 1 0.20 + 1 0.18+ 1 0.17+ 1 0.15+ 1 0.14+ 1 0.12+ 1 0.11+1 0.95 0.82 0.71 0.62 0.54 0.48 0.42 0.37 0.33 0.29 0.25 0.22 0.19 0.16 0.13 0.11 0.29 - 1 0.21 - 2
437
W. Zhang, I. Gladweli / Bifurcation phenomena
AIN SOLUTtON
Legend CURVE 1
MAIN CURVE
--0
INVERSE LWMETER
l/d
Fig. 2. Solutions in D-l/d
k
POINT 2 ---
plane.
MAIN SOLWON
UPPERBOUNDARY II.
SINGULARITY
2 p % 0 iii Q
1.5-
LOWER BOUNDARY Legend MAIN CURVE
1 0
100
200
I 300
I 400
I 500
SQLUTICN WEQICN I I I 300 700 300
RADIUS OF CURVATURE Rs
Fig. 3. Solutions in B- R, plane.
CURVE 1 CURVE 2
438
W. Zhang I. Gladwell / Bifurcation phenomena
CO5NBF for D < 1.1. The technical detail is given in Section 2.3. In Table 1 we list the results for 0.2 < D < 3.5. For 0 < D < 0.2 limited machine precision prevents success in the ODE solver because of large solution values. There is no solution to the initial value u,(D) when D = 0 due to violation of physical constraints. For 2.2 < D G 3.5, the range where the solutions were computed successfully in H, our computational results agree to the fifth digit. For values D i 3.5 OUT numerical results (not presented) are similarly accurdte. The point D = 2.2 (referred as the connecting point in Fig. 2) is where the numerical results of [5] stop. As shown in [5] there is a turning point of d. It occurs at D = 2.914 and the corresponding minimum diameter is d, = 7.768. We call the curve above this turning point the upper branch and the one below it the lower branch. It is shown in [“I that the upper branch ends at the ideal detonation velocity D, = 6.1 where d + +m and AcJ, A,, --) 1. By computing the corresponding R, from (1.6), we observe that at D = 6.0, R, = 910.8 and it decreases monotonically as D decays until D = 2.7 where R, = 9.033, which is the minimum. R, increases monotonically as D further decreases. The behavior of R, is similar to that of d but differs from the R, of the simplified problem [3,13]. This difference contributes to the different solution phenomena discussed below. From our intermediate output we observe that for all D 2 2.1, variables u, A, A, m, x and v increase monotonically as t + T,, while p decreases monotonically. On the lower branch, as D tends to zero, R,, d, T,, and xc-- + +*, A,, + 1 and A, increases monotonically as shown in Table 1. Physically xcJ --) +a implies the reaction zone has an infinite length, which is physically unattainable. We illustrate their behavior since the information helps in finding what happens and where the “nonphysical” phenomenon starts. For D < 2.0, variables u and m grow fast near t = T,,. As D gets smaller u becomes flat over almost the whole integration interval and then climbs rapidly to its CJ value. This implies that the second CJ condition (1.23b) is almost satisfied over the whole integration interval and there is a boundary layer at t = T,, where the first CJ condition is also satisfied. The boundary effect becomes more severe as D + 0 revealing the possibility of instability. A similar trend also appears in variables A, A, cc)and p. Note ucJ > D when D G 1.65. In contrast, ucJ 1.65. A new and interesting phenomenon found in our computations is non-isolated solutions of the system We found a solution region corresponding to a surface in (R,, D, 7’J space. The region is shown in Figs. 2 and 3 and is indicated by III. The solution on the boundary of Region III in Fig. 3 is listed in Tables 2 and 3. Though we list results for D >, 1.1, our computations can be continued for D < 1.1. This new region appears at D = 2.4355 and extends to R, + + 00 and D -+ 0. Near the vertex of this region the solution variables and T,, change smoothly. As D decreases from 2.4 a discontinuity gradually appears in T,, at the upper boundary (UB) of Region III as shown in Fig. 4. The corresponding solution variables also exhibit a discontinuity on UB (see Fig 5). This means the solution of the system becomes isolated on UB away from the vertex while non-isolated solutions are inside III and on LB. Detailed results in [14] show that, like in the main solution, there is a boundary effect in the non-isolated solutions as D --) 0. This suggests that the continuum of solutions inside III and LB becomes unstable. Such phenomena are never found in the simplified model [3,4,13]. Instead a limiting point in D exists 1131and the solution is always an isolated curve in (R,, D, T,,) space. Observe from Table 2 that as D -+ 0, A, decreases to 1 and A,. -+ 0 on UB. On LB, A, increases and A,., 4. stays small but nonzero. The trend A cJ + 1 as D -+ 0 on UB is similar to
W. Zhang, I. Gladwell / Bifurcation phenomena
439
Table 2 Numerical results on the upper boundary CUB) of the continuous solution Region ZII D
Rs
TCJ
d-
kJ
ACJ
*CJ
%J
-kJ
PCJ
2.435 2.400 2.300 2.200 2.100 2.000 1.900 1.800 1.700 1.600 1.500 1.400 1.300 1.200 1.100
0.28 + 2 0.38 + 2 0.55 + 2 0.78 + 2 0.11+3 0.17+3 0.27 + 3 0.43 + 3 0.72 + 3 0.12+4 0.21+ 4 0.39+4 0.74 + 4 0.15 +5 0.30 + 5
0.20 + 1 0.20 + 1 0.25 + 1 0.37 + 1 0.55 + 1 0.86 + 1 0.14+2 0.23 + 2 0.40 + 2 0.70 + 2 0.13+3 0.25 + 3 0.50 + 3 0.10+4 0.23 + 4
0.19+2 0.24 + 2 0.33 + 2 0.47 + 2 0.69 + 2 0.10+3 0.16+3 0.26 + 3 0.43 + 3 0.73 + 3 0.13+4 0.23 + 4 0.44 + 4 0.87 + 4 0.18+5
2.1864 2.1588 2.0831 2.0066 1.9288 1.8496 1.7689 1.6867 1.6029 1.5176 1.4308 1.3425 1.2527 1.1616 1.0691
1.0801 1.0532 1.0418 1.0375 1.0341 1.0312 1.0286 1.0261 1.0237 1.0213 1.0131 1.0170 1.0151 1.0133 1.0116
0.10958 0.09495 0.07968 0.06861 0.05917 0.05100 0.04384 0.03751 0.03194 0.02”34 0.02;!75 0.0151;3 0.01574 0.01291 0.01046
0.16- 1 O.lO- 1 0.60 - 2 0.37 - 2 0.23 - 2 0.13 - 2 0.74 - 3 0.40 - 3 0.21- 3 0.11-3 0.51-4 0.24 - 4 0.10-4 0.43 - 5 0.17-5
0.40 + 1 0.38 + 1 0.49 + 1 0.68 + 1 0.99 + 1 0.15+2 0.23 + 2 0.37 + 2 0.60 + 2 0.10-t-3 0.18+3 0.32 + 3 0.60 + 3 0.12+4 0.24 + 4
0.70 0.66 0.56 0.48 0.40 0.34 0.28 0.23 0.18 0.15 0.11 0.89 0.68 0.510.37 -
-
1 1 1 1
what happens on the main solution curve as D + Dr. We observe no boundary layer effect as D + D, on the main solution curve and UB. From Table 2 ucJ
Rs 0.26 + 2
L o-21+ 1
d
2.435
0.18+2
4T.J 2.1869
2.400 2.300 2.200 2.100 2.000 1.900 1.800 1.700 1.600 1.500 1.400 1.300 1.200 1.100
0.27 + 2 0.33 + 2 0.42 + 2 0.54+2 0.76 + 2 0.11+3 0.17+3 0.26 + 3 0.42 + 3 0.71+3 0.12+4 ‘X24+4 0.48 + 4 0.95 + 4
0.31+ 1 0.66 + 1 0.12+2 0.21+ 2 0.38 + 2 0.73 + 2 0.15+3 0.31-t 3 0.69 + 3 o-17+4 0.44 + 4 0.13+5 0.44 + 5 0.14+6
0.21-k 2 0.32 + 2 0.49 + 2 0.74 + 2 0.12+3 0.20 + 3 0.36 + 3 0.69 + 3 0.14+4 0.31+4 0.75 + 4 0.20 + 5 0.62 + 5 0.18+6
2.1622 2.0900 2.0171 1.9436 1.8689 1.7930 1.7161 1.6382 1.5591 1.4790 1.3978 1.3158 1.2332 1.1531
&I 1.0897 1.1270 1.2124 1.2963 1.3854 1.4806 1.5920 1.7295 1.9036 2.1296 2.4305 2.8430 3.4289 4.2985 5.4641
*cJ 0.11303 0.11800 0.12346 0.12603 0.12792 0.12697 0.12497 0.12295 0.12106 0.11939 0.11800 0.11699 0.11648 0.11664 0.11793
*cJ 0.18- 1 0.18- 1 0.17 - 1 0.15 - 1 0.14- 1 0.12- 1 0.98 - 2 0.79 - 2 0.63 - 2 0.50 - 2 0.39 - 2 0.31- 2 0.24 - 2 0.18-2 0.15-2
kJ 0.42 + 1 0.62 + 1 0.13+2 0.23 + 2 0.38 + 2 0.67 + 2 0.12+3 0.23 + 3 0.47 + 3 0.10+4 0.23 + 4 0.57+4 0.16+5 0.50 + 5 0.15+6
PCJ 0.70 0.67 0.58 0.50 0.43 0.36 0.31 0.25 0.21 0.17 0.14 0.11 0.81- 1 0.60 - 1 0.44 - 1
Wt’.Zhang, I. Gladwell / B4jkrcation phenomena
440 10-J
o-
I
I I I I I
8t-
I I
B-
F
Legend
I 5-
0
I I
4-
z
3 I 31 1 25
CURVE 1 CURVE 2
I I I
A
CURVE 3
+
CURVE 4
X
CURVE 5
0
CURVE 6
+
CURVE 7 ---
X
CURVE 8 --c-m
C’ m-e CURVE 9
I
I
I
I
I
I
30
3s
;I:
45
50
55
Fig. 4. Solutions in TcJ -R, plane.
20
1
18 1 ,
I
i
I I I I
I I I
I \
i
I I I I I i
Legend
I
I
0
I I I
A CURVE 3
b
CURVE 2
+
CURVE 4
X
CURVE 5
c3 CURVE 6 -I- CURVE 7 --_ X
2 1 25
I 30
I 35
I ii
I 46
Fig. 5. Solutions in Xc. - Is, plane.
I so
CURVE 1
I 55
CURVE 8 s---w
Cl CURVE 9 ---
441
Legend L!ts D
UB MC
Id0
150
200
2Ijo
400
450
560
W
Fig. 6. Solutions in D- Xcr plane.
from Region III in the (R,, D) and (l/d, D) planes. By connecting the main solution curve to UB, we would obtain a stable isolated solution curve for all D E (0, D,) which exhibi& no boundary layer effect and no violation of ucJ < D. Our limited numerical and graphical experiments seem to confirm this behavior. A schematic ( xCJt D) graph is given in Fig. 6. !n the (xc.D 0) plane, UB is entirely separated from the rest of solutions in III. If we map Region III to the (l/d, D) plane it is much narrower. Since for each pair (R,, E)) there is only one possible solution, just one value of d can be obtained from (1.6). But the converse does not necessarily hold due to the presence of xc3 in (1.6). Our computations show that in Region III solutions with different R, for a fixed D may result in the same diameter d. We use a bold curve on the (l/d, D) plane to indicate this solution region in Fig. 2. On this bold curve more than one solution may exist, I-Iowever only one solution is isolated and “stable”. 2.3. Regiiunal properties and their application Since d does no occur in (1.20)-(1.23), we concentrated on the (R,, D) plane. We performed computations in Regions I and II shown in Fig. 3. The shooting strategy employee (Section 2.1) to solve the singular ODE BVP allows us to perform the integration at any point in the (R,, D) plane. We simply pick a pair (R,, D) and integrate from t - 0. We let the integration proceed as far as it can and monitor the residues of the two CJ conditions (1.23). In ILegion I we find that after continuing the integration for a sufficiently long time the second CJ condition (1.23b) is satisfied and all variables except x tend to constant values. As a consequence the residue of the first CJ condition 41.23a) also tends to a nonvanishing constant. So the c:ondition is, never satisfied in Region I. In other words, as the integration is advanced, the ntunerator in (1.20a)
442
W: Zhang, 1. Gladwell / Bifurcation phenomena
tends to zero before the denominator vanishes. Then the whole system reaches a pseudo steady state and there is no further opportunity for the denominator to decay to zero. Mence there is no solution to the BVP in this region. Since all quantities but x tend to constants for a large t, we call this a steady state region, see Fig. 3. In Region II, we find the integration is always interrupted ivy the singularity and the second CJ condition is never satisfied. Thus the denominator rather than the numerator in (1.20a) reaches zero first. An QDE solver cannot integrate across such a singularity. Again there is no solution to the BVP. We call this a sing&r region (see Fig. 3). When this work was completed the authors learnt that Lynch INI had independently observed the regional properties and used them in numerical experiments. The main solution curve in the (&, D) plane is the boundary between the singular and steady state regions. On the curve both the numerator and the denominator in (1.20al tend to zero simultaneously. The existence of this solution curve seems natural after the identification of Regions I and II. In our computations we observe the immediate switch from one region to another as R, changes. This suggests to us that the value R, which yields a solution to the BVP is isolated. An alternative method for locating the main solution curve is to replace COSNBF by a bisection method. For a fixed D we locate an interval of R, with one end in the steady state region and the other in the singular region. Then we apply bisection based on the two different regional properties. This technique is used in our computations for D Q 1.0 where conventional continuation becomes expensive. Region III is surrounded by the singular Region II. In Region III we find that not only are the two CJ conditions satisfied simultaneously but also the integration can be continued beyond TcJ. Although mathematically the singularity always exists in (1,20a), numerically it may not be detected. As the numerator in (1.20a) tends to zero faster than the denominator (or at the same order), u& is finite. Therefore the integrator does not “see” the singularity and can continue the integration beyond it. When this happens the nonlinear equation solver adjusts the length of the integration interval back to T,, in order to match the two CJ conditions. Computationally this phenomenon is often observed deep in the interior of Region III. Since inside this region there exists a solution for each CR,, D), any R, selected by the nonlinear equation solver is an exact solution. The boundary of III is a fold [8]. In general, more work is required to trace a fold since to define it usually extra equations are added to the original system. In our case we can exploit the regional properties and locate the fold using an amount of work sinGBar to just finding a solution. It turns out that the continuation method succeeds in computing solutions on UB but not on LB. This is due to the solution (R,, T,,) being isolated on UB, which pcevent the nonlinear ell,uation solver searching for a solution inside Region III. Since there is no boundary layer effect on UB, a relatively large continuation step size can be used. When continuation starts on LB, the nonlinear equation solver finds a solution inside Region III rather than traces LB and it finds different solutions when using different continuation steps. So, the bisection method is appropriate for locating LB. 2.4. Stability of the continuum of solution phenomenon We compare the results to those for the simplified problem [4,13]. For large values of D they are qualitatively same. There is a turning point and an upper solution curve in the (l/d, 6))
W. Zhang, I. Gladwell / Bifurcation phenomena
2.8
P
443
MAIN SOLUTION
2.6
z G
0 d
>
II.
2.4
SINGULARITY
I
UPPERBOUNDARY
1.81
Legend CURVE 1
1.6
0
CURVE 2
RADIUS OF CURVATURE Rs
Fig. 7. Solutions in D-R, plane.
plane which ends at the ideal detonation point Dr. For lower values of D the two sets of results are different. For the simplified problem, [!3], the solution curve ends at the lowest critical value of D =D,_( > 0). And as D -+D,, d + +q xcJ + +m and A,, + 1. There is no other solution curve/region in the range (D,, Di). For the current problem there is no value DL( > 0) at which the solution curve/region ends. As Ei + 0, h, + 0 on LJB while A, > 0.7 in [4,13]. This suggests that the rate of the increase in A may be an important factor in determining the bifurcation diagram for lower values of D. 510we tested problem (1.20)-(1.23) with perturbations to (1.20~). First we varied A’ in (1.20~) by inserting a multiplier m, A’ = in@-
A)F(A, p).
(2 .2)
For m = 1.2, 1.5 and 2.0 we computed the solution with the same strategy and error tolerances as in Section 2.1. Our numerical results show that as m increases, Region III in the (R,, D) plane moves closer to the main solution curve. At both m = 1.5 and m = 2, III is connected with the main solution curve as shown in Fig. 7 for m = 2. Region III, instead of being surrounded by II, lies between I and II. When we map to the (l/d, D) plane, Region III connects directly to the main solution curve as shown in Fig. 8. Our second change arises from observing the different dependences on p in F. In [4,13], F is modeled by linear dependence on p and the numerical results always show pcJ 3 1. In the present problem F depends linearly on p for p large but is 0( p4) for p small, Cl.11). The switching point is close to p = 1. The two pieces are smoothly connected. Our computations show that pcJ drops below 1 for lower values of D, so the Q(p4) term yields a much smaller value of F and hence a small A’. We changed F by reducing the power of p to cubic and
CK Zhang, I. Gladwell / Bi&rcation phenomena
444
MAIN SOLUTION CURVE
3 1
Legend CURVE Cl REGION
PNVERSE DIAMETER Fig. 8. Solutions
lid
in D-1 /d plane.
quadratic respectively in the expression for pp, (1.1 lb) whilst preserving the smooth connection to the linear piece [14]. Then we computed the solutions to the perturbed prob!em with the same strategy and error tolerances. We obtained the same solution pattern as in Fig. 7. Region III connects to the main solution curve and lies between I and II. Though the solution values are different, the pattern and the behavior are the same.
We have transformed a parameterized index one DAE BVP with index two on the free boundary into a parameterized singular ODE BVP with a free boundary. By using a shooting technique for solving the singular free boundary value problem, we successfully computed the sob on to the detonation problem studied in [5,9] and presented and justified the complete bifurcation diagram. Because of the difficulty of treating a singularity on the free boundary, we identified regional properties in parameter space. These properties were used in locating the isolated solution and the boundary of a contiwum of solutions, which turned out to be computationally efficient. Bifurcation phenomena of related detonation problems were also presented tikh demonstrates that the phenomena persist.
eferences [l] U. Ascher, Private communication. [2] J.B. Bdzil, Steady s!:?!c. two-dimensional
detonation,
J. Flrrid Me&. 1198(1981) 195-226.
W. Zhang, I. Gladwell / Bifurcation phenomena
445
[3] CC. Beardah and R.M. Thomas, Techniques for the numerical solution of a singular index one differential/
algebraic equation boundary value problem arising in the modelling of unconfined detonations, Numerical Analysis Report, No. 191, Department of Mathematics, UMIST (1990). [4] M. Braithwaite, T. Farran, I. Gladwell, P.M. Lynch, A. Minchinton, I.B. Parker and R.M. Thomas, A detonation problem posed as a differential algebraic boundary value problem, Math. Engrg. Ind. 3 (1) (1?90) 45-57. [S] M. Braithwaite, T. Farran, I. Gladwell, P.M. Lynch, A. Minchinton, I.B. Parker and R.M. Thomas, Private
communication. [6] K.E. Brenan, Sk. Campbell and L.R. Petzold, Numerical Solution of Initial-Value Problems in Differential Algebraic Equations (Elsevier, New York, 1989). [7] W. Fickett and W.C. Davis, Detonation (University of California Press, Berkeley, CA, 1979). [S] H.B. Keller, Numerical Methods in Bifurcation Problems (Springer, Heidelberg, 1987). [9] I. Kirby and G. Leiper, A small divergent detonation theory for intermolecular explosives, in: Proceedings 8th International Symposium on Detonation, Albuquerque, NM (1985). [lo] I. Kirby, G. Leiper and A. Hackett, Determination of reaction rates in intermolecular explosives using the electromagnetic particle velocity gauge, in: Proceedings 8th International Symposium on Detonation, Albuquerque, NM (1985). [ 1l] P.M. Lynch, Private communication. [12] The NAG Fortran Library Manual - Mark 13 (1st ed., 1988). [13] W. Zhang and I. Gladwell, Analysis of a simplified detonation problem, Math. Engrg. Ind. (in press). [14] W. Zhang, and I. Gladwell, Bifurcation phenomena arising in a singular differential algebraic equation detonation problem, Report No. 90-14, Department of Mathematics, Southern Methodist University, Dallas, TX (1990).
427
Wen Zhang and Ian Gladwell Mathematics Department, Southern Methodist University, Dallas, TX 75275, USA Abstract Zhang, 6~. and 1. Gladwell, Bifurcation phenomena in a detonation problem, Applied Numerical Mathematics 9 (1992) 427-445. In [5,9] a system of paraneterized differential algebraic equations (DAEs), modeling a detonation problem, is posed and its solution computeJ. Due to computational difficulties, these numerical results arc incomplete. A bifurcation curve stops abruptly during continuation in the parameter D, the detonation velocity. In this paper we transform the problem into a parameterized ordinary differential equation (ODE) boundary value problem (BVP), exhibiting a singularity on the free boundary. We present an effective strategy to deal with the singularity. Using this strategy we successfully extend the bifurcation curve. Also we find a new bifurcation phenomenon with a continuous region of so!utions.
I. The detonation problem 1.1. Introduction We investigate a ZND [7] model of an unconfined small divergent detonation along a cylinder. The problem is illustrated in Fig. 1 which is a modified version of the corresponding figure in [4].
Unreacted Explosive
Reaction Products
Fig. 1. Detonation along a cylinder.
016%9274/92/$05.00
0 1992 - Elsevier Science Publishers B.V. All rights reserved
W. Zhang, I. Gladwell / Bifurcation phenomena
428
coordinate system is attached to the shock front which moves at a constan: velocity D into the unreacted explosive. We investigate the reaction zone based on one-dimensional flow theory in Lagrangian coordinates. The position of the CJ plane depends on the solution. Kirby and Leiper [9] gave the physical background and the general formulation of this problem. Braithwaite et al. [4] further described it in mathematical detail. The system governing the reaction zone is an autonomous DAE BVP with time t as the independent variable. It has nine dependent variables, (u, A, A, w, X, p, v, p, E), corresponding to five ODES and four algebraic equatiins (AEs). There are four undetermined constants, (D, d, TcJ, xcJ), and eight boundary conditions (BCs) corresponomg to five first-order ODES. Hence one constant can be considered a free parameter and the remainder ar 9 unknowns. It is of physical interest to determine the relationship between the diameter d and the detonation velocity D or the radius of curvature of the shock front R, and D [7]. The numerical results in [4,5,9] show that d, T,, and xcJ are single-valued functions of D but not vice-versa. Hence D is treated as a free parameter for computational convenience. A further discussion on the roles of these constants is given in Sections 1.2 and 2.1. The physical meanings of the variables, the unknown constants and the parameter are: The
l4
-
A
-
A
-
0
-
x
-
P c
-
P E
-
D d TCJ 'CJ -
particle velocity, area of a stream tube of fluid which had unit area degree of chemical reaction (A = 0 for no reaction, radial velocity, distance from the shock front, pressure, sp\=cific volume, fluid density, specific internal energy, detonation velocity, diameter of the cylinder of explosive, time needed to reach the CJ plane from the shock the (CJ) distance between the shock front and the
at the shock front, A = 1 for complete reaction),
front, CJ plane.
The ODES are puu’
+p’
= 0,
(l.Ia)
A’= 2Ao, A’= (I-
A)F(h,
0’ = -u’/R(x, X4=u.
(l.lb) p), +Jv
(l.lc) d),
(1
.ld)
(Lie)
(1. la) is the momentum conservation equation. (1. lb) describes flow divergence and is derived geometrically [4,7,9]. (1 .lc) describes the chemical reaction rate and (Lid) defines the variation of the radial divergence.
W. Zhang, I. Gladwell / By.furcation phenomena
429
The AEs are E =H(p,
v, A, D),
(1.2a)
E -t pu + ;ti2 = +D’,
(1.2b)
Apu - p,D = 0,
(L2c)
v = l/p.
(1.2d)
(1.2a) defines the energy. (1.2b) is the integral of the enthalpy conservation equation combined with the momentum equation, which is Bernoulli’s equation. (1.2~) states mass conservation. The BCs are
u = u,(D),
(1.3a)
A = 1,
(1.3b)
A =o,
(1.3c)
0 =
[D
- ~,ja)l/~,(+J~
(1.3d)
d),
A:=0
(1.3e)
at t = 0 (the shock front) and
u =c(P,
0, A, D),
CL) = (T(p, v, A, c, D)A’/2, x =‘CJ
(1st CJ condition), (2nd CJ condition),
(1.4a) (1.4b) (1.4c)
at t = (the CJ plane). c is the speed of sound and U- is the thermicity coefficient of the fluid. The value is determined as the time when both (1.4a) and (1.4b) are satisfied. (1.4a,b) are the Chapman-Jouguet (CJ) conditions. (1.4a) determines where particle velocity reaches sound speed and (1.4b) balances the release of energy by the chemical reaction and its withdrawal by the radial flow [4,7]. There are physical constraints which must hold over the whole integration interval: TcJ
TcJ
O.
v, A, D),
(1 Sa)
O
(1.5b)
A 2 1,
(1%)
P, u, p > 0.
(1.5d)
and strict inequality holds The equality in (1.5a) can only be achieved at the CJ plane (t = for t E [O, po, the explosive porous density, is chosen as 1.08 as in [4,5]. R(x, p, V, A, D), c( p, v, A, D), a( p, u, A, c, D), F(A, p), us(D) and Rs(xCJ, d) are modeling functions. R is the radius of curvature of an isobar at a distance x downstream of the shock front. M is the specific internal energy and F describes the pressure dependence. R, has the form bee [2,9,101) TcJ)
TcJ)
[9].
xcJ7
R,(X,,,
d)=(d-&,)/a.
d),
M(
(1 .6)
W Zhang, I. Gladwell / Bifurcation phenomena
430
Kirby, Leiper and Hackett [lo] determined a! = 0.5 and p = 1.2 by experiment. We use A
1-A H(p, v, A, D) =PO
c(p,
v, A, D) =
I
y,‘(D) - 1+
PV(V,
YRW - l Phrl(V)[Y(%
A, D) +
[Y&J)
\
a(p,
v, A, c, D) =
y(4 A, D! - 1 c2
Q +-E,(D)
i
. (18)
- AQ + (1 -A)&(D), 4 -
D) - 11
II2
1’2
. (19)
’ i
1 -PO
Y&J
- 1 - rm
- 1
(1.10) and for the emulsive explosive considered here F(AY P) = Q(A)PP(
(1.1 la)
P)/q-j
where if p>,&,
(P-Pm PP(. P) =
:P
3p 4p,
3 ’
(l.llb) otherwise.
ii I Q, the heat of reaction, and pH, the hot spot critical pressure, are experimentally determined constants. yU(D), the index of the unreacted ingredients; E,(D), the turbulent stored energy; y(v, A, D) and u,(D) depends implicitly on D. As shown in [5] they can be determined by iteration. Y&J) and aH( A) are explicit functions. A detailed description and explanation of these functions is given in [5,14]. In [3,4,13] the problem is investigated with simplified function definitions for R( X, HCp, u, A, DJ,c(p, u, A, D), d p, v, A, c, D), F(A, p) and u&D). Since the relationship between the unknown constants and the parameter is of practical importance, our view is focused on them. represents the solution variable x(t) at t = T,, and can be considered a representative solution variable. The projection of a solution to the DAE BVP (1 .l)-( 1.l 1) in parameter space, such as the (l/d, D) plane, is defined as a solution in that space. Kirby and Leiper show the (l/d, D) diagram in [9]. Their computed curve stops abruptly at D = 0.40: where D, is the ideal detonation velocity. They did not succeed in computing the curve for lower values of D. Braithwaite et al. [5] computed the solution with a similar result. We continue the investigation of this problem below. +J,
d),
.xcJ
1.2. Reformulations of the detonation problem and discussions on the singularity Braithwaite et al. [4,5] computed the solution by solving for the nine unknown variables in the DAE BVP (l.l)--(1 S). The index [6] of this system with respect to the nine unknown variables is shown in [4,5] to be one everywhere except at the boundary t = where it is two. We observe that the size of the system can be reduced as in 1131.Since pO, D > 0, (1.2~) and ( 1.2d) imply &J
V =-
AL1 (I
PqD’
.12)
W Zhang, I. Gladwell / Bifurcation phenomena
431
Substituting (1.2a) into (1.2b), only one AE is needed and the number of unknown variables is reduced from nine to six (u, A, A, w, X, p). p, o, and E, can be obtained directly from them. Therefore the system is reduced to (l.lb)--(l.le), p&4’ +p’A = 0,
(1.13)
and
-
= --
(1.14)
with BCs (1.3) at t = 0, and, at t = TcJ, (1.4c), and (1.15a)
U
Au
------,A
(1.15b)
P, D
The functions H, c and CTare always defined under the rule that the adiabatic equation, y’= -pu2(aA’ - 20)/(1
- u”/c”),
(1.16)
must be satisfied, [9]. In other words, if we differentiate AE (1.14) and combine with the other ODES, we can always obtain (1.16). Thus the index of the reduced DAE system is one everywhere except where the denominator of (1.16) is zero. The denominator only vanishes at t = 7’&, [9], and it can be shown, as in [4,5], that there the index is two. Hence the reduced system has the same index as the unreduced system. One way to solve the reduced DAE BVP is to replace AE (1.14) by (1.16) and solve the enlarged ODE BVP for the six unknown variables (u, A, A, W, x, p) [6]. Solving this way, (1.14) may be satisfied inaccurate!y. Since (1.14) reflects energy conservation, its accuracy is of physical importance. We observe from (1.8) that p in (1.14) can be determined explicitly. Using (1.8), we write (1.14) as
p=P(A,u,A,D,P(,)=
Dpo[ D2/2 - u2/2 + AQ - (I-A)&] . l-h A +1 \ Au --l-----Y” - 1 YR- 1 J
Substituting (1.16) into (1.13) with l/u for ld =
uc2(crA’- 2w)/(c2
- u2).
p
(1.17)
and (1.12) for u, we have (1.18)
Replacing p by (1.17) wherever it appears, equations (1.18) and (I.lb)-(l.le) with BCs (1.3), (1.15), and (1.4~) form a BVP of five ODES with five unknown variables (II, A, A, w, x). (1.14) is satisfied throughout due to the substitution for p. Equation (1.18) exhibits singularity whenever u = c though we do not perceive it in the original form (l.l)-(1.4). We believe this singularity is the cause of the difficulties in the computations reported in [4,5]. There are two ways to deal with it. One way is to bring it out in the ODES, as above. An alternative way, as suggested by Ascher [l], is to keep it inside the AE (1.14) then to form a new set of variables and embed (1.14) in the ODES, [14]. Though we do
W. Zhang, I. Gladweld / Bifmation phenomena
432
not observe the singularity in the new ODE system, this does not mean th-:I it has been removed from the problem. Instead it is converted into a solvability obstacle in the AE. This is demonstrated independently in [3]. The computation is expensive. For instance, if we use shooting methods to solve the new ODE BVP, we need to solve the embedded AE in each call to the derivative evaluation routine. When we cannot solve the AE we cannot proceed with the integration. We choose to work with the form which exposes the singularity in the ODES because it is directly related to the physical phenomena and has an important advantage in investigating the nature of the solution. We observe from (1.7) that R(x, xcJ, d) = R(x, R,(x,, d)) and xcJ, d appear only in (l.ld) and (1.3d) through R, (without taking account of (1.4~)). We introduce a new unknown corzsturzt R, to replace xcJ and d. Since xcJ =x(7’&, it is automatically determined as soon as the solution is found. Then d can be computed directly from (1.6). Thus, the undetermined constants are reduced to R,, D and T,,. The BCs are coupled through xcJ in (1.3d) which is posed at t = 0. After introducing R, (1.3d) simply becomes
40) = [D -
(1.19)
u,(D)l/Rs
and (1.4~) is no longer needed. Therefore the BCs separate. Replacing R(x, xcJ, d) by R(x, R,) in (l.ld), (1.3d) by (1.19) and discarding (1.4~)~ we finally obtain a BVP with five explicit ODES (1.18), (1 .lb-e) corresponding to five unknown variables (u, A, A, o, X) (replacing p by (1.17) wherever necessary); and seven separated BCs (1.3a)-( 1.3c), (1.19), (1.3e) and (1.15) corresponding to five first-order ODES with three unknown constants (R,, D, T,,). As before, one of the constants must be considered as a free parameter and the other two are determined by the solution of the ODE BVP. As shown in Section 2.1, this system is well-posed in a sense that for each (R,, D, T,,) there is at most one solution of the ODE BVP. Further, each solution corresponds to a unique (R,, D, Z’J. Hence the solution of the system can be represented in (R,, D, T,,) space. As there is only one free parameter, the solution expected in CR,, D, T,,) space is a curve. Physically either D or R, can be viewed as the free parameter. Computationally selecting D is advantageous. We restate the reduced system below. The ODES are U'
= uc2(crh’ - 2w)/(c2 - u2),
(1.20a)
A’ = 2Ao,
(1.20b)
A’= (1 - h)F(A, p), 0’ =
(1.2Oc)
-u’;R(x, R,),
(1.20d)
x’=u,
(1.20e)
where P =-
Dp,[ D2/2 - u2/2 + AQ - (1 - A)E,] l-h AU
-+
Yu-
.
A -+1 1
YR-
1
(1.21)
W. Zhang, 1: Gladwell / Bifurcation phenomena
433
The BCs are I/ =
u,(D),
(L22a)
A = 1,
(1.22b)
=o,
(1.22c)
A
0 = [D-u,(D)]&
(1.226)
X=0
(1.22e)
at t = 0, and U
- A)&
p)/2,
(1st CJ condition),
(1.23a)
(2nd CJ condition),
(1.23b)
at t = T,,. The physical constraints (1.5) still hold here. Mathematically there is a singularity in (1.20a) and it occurs at the free boundary t = T,.
2. Numerical simulations and bifurcation phenomena 2.1. Computational strategy for solving the parameterized singular ODE BW There are a variety of numerical methods which might be used to solve the ODE BVP (1.20)-(1.23). We choose a shooting method for following reasons. First there are five BCs (1.22) at the left boundary, t = 0, and five explicit ODES (1.20). These constitute a well-defined initial value problem for given (R,, D). They are already in a standard form for calling an ODE solver. Since the right end of the interval is a free boundary, we may advance the integration as far as needed and T,,, the right boundary position, is determined as the point where the CJ conditions, (1.23) are satisfied. The aim of shooting is to satisfy the two CJ conditions simultaneously. For given R, and D there is at most one 7& because we always terminate the integration the first time (1.23) are satisfied, even if the integration can proceed. It is nonphysical for the first CJ condition (1.23a) to be satisfied in (0, T,,) [9]. So the solution is unique for given R, and D. Using a shooting method we have the opportunity to observe what happens even when a choice R,, D leads to no solution. With this facility we can produce a bifurcation diagram and distinguish different regions in the (R,, D) plane. This information is useful in determining the solution manifold. A simple implementation of a shooting method is to use an ODE solver to perform the integration and use a nonlinear equation solver to match the BCs. There is a choice of the direction of the integration. We choose to integrate from t = 0 to t = T”,. In addition to the reasons above, the decision was also based upon experiments in [4,5]. There a shooting method is also employed but integrating from t = TcJ to t = 0. This choice was made to avoid integrating towards the singularity but there is a difficulty in finding the solution We need to determine ~1’at the singular point t = T’,, in order to start the integration. This requires solving
W Zhang, I. Gladwell / Bifurcation phenomena
434
nonlinear equations involving values of the solution variables at t = T,, which are the initial estimates in the shooting method. Hence the selection of initial estimates must guarantee a solution for u’(T,,), otherwise the integration cannot start. After the shooting method is started,. the nonlinear equation solver adjusts the variables at TcJ internally and there is no guarantee of a solution u’( T,,). Another frequently occurring difficulty in [4,5] when computing the solution for smaller D is the integration stopping before t = 0. This is due to the ODE solver detecting a singularity in the interior of the integration interval. After this we cannot proceed to the matching phase. In order to avoid these interruptions we chose to integrate from t = 0 where the derivatives are well-defined. Since the right boundary is free, it is natural mathematically and physically to end the integration where the singularity occurs. Thus we can always apply the matching procedure after the integration, and then the shooting method continues. In this way we have succeeded in computing the solution on the lower branch in D where previous computations failed. To perform the shooting method we have used standard codes from the NAG Library [12], D02EBF for performing the ODE integration and C05NBF for solving the two nonlinear BCs (1.23). The ODE solver requires as input a fixed integration interval. So, we rescale the original integration intervai to [O,l] by setting t^= t/7& as the new independent variable. The five ODES become u’ = T&C2(CTAf/TCJ- 2w)/(c2
x’ = Tg,
(2.la) (2.lb)
A’ = 2T,,Ao, A’= T,,(I - h)F(h, w’ = -UR(x, R,),
- U2j,
P),
(2.lc) (2.ld) (2Se)
where p is computed from (1.21). There is no change in the written form of the BCs (1.221-t 1.23) except that (1.22) hold at r^= 0 and (1.23) hold at t^= 1. D02EBF uses a variable step size variable order backward differentiation formula to integrate a system of stiff ODES. Here the stiffness property is of no importance. A nonstiff solver will do the job as well. It would become appropriate if the problem was made stiff by the introduction of chemical kinetics terms in a “full” model of the detonation problem. The code has a feature that it outputs the solution at its last successful step when the integration cannot proceed to the end of the interval. This feature makes it possible to use the strategy below for continuing shooting after encountering the singularity. We proceed as follows. Since there are three undetermined constants and two matching equations, we fix one constant as a parameter and search for suitable values for the other two until the two matching equations are satisfied, that is a solution to the ODE BVP is found. We rule out fixing T,, because it makes no sense to fix a free boundary and can only reduce the chance of success in searching for a solution. We choose to fix D because it saves computational work compared to the alternatives. We observe that u,(D) in the initial data, and y,(D) and E,(D) in the function definitions are obtained via solution of a nonlinear equation [5]. If we fix R, we have to solve the extra nonlinear equation each time we change D inside the iterationThere is no such requirement when R, changes. If D is fixed we can compute these quantities before starting shooting and do not need to recompute them,
W. Zhang, I. Gladwell / Bifurcation phenomena
435
We call DKEBF to integrate equations (1.20) with the initial conditions (1.22). There are two possible outcomes. One is that the integration is advanced to t^= 1 without detecting any singularity. In this case we simply evaluate the CJ conditions (1.23). If they are satisfied to the given tolerance, we quit the shooting procedure since the solution is found. Then d is computed from (1.6) where xcJ, the solution x(l) at ? = 1, is available from D02EBF and R, is available from COSNBF. If the CJ conditions are not satisfied, COSNBF provides a new set of estimates for R, and QJ. The other possible outcome is that the integration stops at a point before t^= 1 due to encountering the singularity, that is the first CJ condition (1.23a) is satisfied but the second may not be. The code outputs the result at a point immediately preceding the singularity. We still check the residues of the two CJ conditions using output from D02EBF since it is a reliable result close to the singularity. If the CJ conditions are satisfied within the given tolerance, the solution is considered found. It does not matter that the integration does not proceed to t^= 1 since we are dealing with a free boundary. Let the value of Tcr output from COSNBF be 7’&. Because the integration ends before t^= 1, the length of the integration interval in the original variable t should be T’j = T&i. d can be computed in the same way as before. If the two CJ conditions are not satisfied, we begin the cycle again. For the above strategy the singularity is not a barrier to finding a solution. Instead we use it to prevent unnecessary work. When the estimate of T, from COSNBF is too large, we adjust it to the correct value. So the integration not only finds the solution variables but also finds the correct value of one of the unknown constants, T,,. This is essentially due to integrating towards a free boundary. Clearly the capability of accurately computing the solution “at” the singular point is crucial to the success of this strategy. We have computed the solutions of the problem using the above strategy on an IBM 3081D. All the computations are done in DOUBLE PRECISION with a mixed error tolerance 1.5 x lo-l2 assigned in both D02EBF and COSNBF. We assign a third tolerance 5 X 10m8 on the residues of the two CJ conditions (1.23) which are checked after each call to D02EBF. If they are satisfied within the tolerance, we quit the iteration and claim a solution is found. The reason for using such a stringent tolerance in the call to D02EBF is to limit the error in the integration procedure as much as possible in order to track the solution closely and to prevent violation of the physical constraints. The existence of a solution to the singular BVP depends on whether the two CJ conditions (1.23) can be satisfied simultaneously. Depending on the sensitivity of the CJ conditions, a “large” error in the integration may prevent US in finding a solution as we have seen in a related problem [13]. Also when using a crude integration tolerance, the residues of the CJ conditions become large. Usually we cannot expect the CJ conditions to be satisfied more accurately than the integration result. So we use a stringent tolerance in the ODE solver but a less stringent tolerance in checking the CJ conditions. The value of the tolerance for COSNBF is less important. The recommended value given in the NAG manual 1121is approximately the square root of the machine precision. The value we select is smaller, in order to ensure enough iterations are taken. In the code the difference between the two successive estimated solutions, rather than the residue of the function, is checked against the given error tolerance. So there is a chance that upon exit the residue of the nonlinear equations may not be as small as anticipated. This has happened with a larger tolerance in C05NBF. Using the tolerance selected above, we obtain consistently satisfactory results from the code.
W. Zhang, I. Gladwell / Bifurcation phenomena
436
We have monitored intermediate solution results and the residues of the two CT conditions throughout the integration. We observe a reduction of several orders of magnitude in the ‘_.c;i-icethat both CJ conditions are satisfied without residues as T,, is approached. This is ev!d ambiguity. Monitoring the behavior of the solution variables throughout the integration interval will be useful below. 2.2. Numerical results and analyses We have successfully traced the solution of the original probiem for lower values of i) than in [5]. We refer to this result as the main solution curve. The values of the solution at t = TcJ are referred to as the CJ values. The results are shown in Table 11and Figs. 2 and 3. For D < 1.6 we need a small continuation step in D, a few more iterations and a more accurate guess for I?,, but, in principle, we do not encounter any difficulty computing the solution. However for D < 1.1 the successful continuation step size becomes so small that it requires too much computation to trxe D further. Based on regional properties determined from our numerical results, we empioy the bisection method to replace the nonlinear equation solver Table 1 Numerical results for the main solution curve
D
Rs
TCJ
3.5 3.4 3.3 3.2 3.1 3.0 2.9 2.8 2.7 2.6 2.5 2.4 2.3 2.2 2.1 2.0 1.9 1.8 1.7 1.6 1.5 1.4 1.3 1.2 1.1 1.0 0.5 0.2
0.12+2 0.11+2 0.11+2 0.1@+2 0.99 + 1 0.96 + 1 0.93 + 1 0.91+ 1 0.90 + 1 0.91 -I-1 0.95 + 1 O.E0+2 0.12+2 0.14+2 0.17+2 0.20 + 2 0.25 + 2 0.30 + 2 0.36+2 0.45 + 2 0.56 + 2 0.71+ 2 0.92 + 2 0.12+3 0.17+3 0.25 + 3 0.64 + 4 0.11+8
0.95 + 0.98+ O.lO+ O.ll+ O.ll+ 0.12+ 0.12+ 0.13+ 0.15+ 0.16+ 0.19+ 0.23+ 0.31+ 0.44+ 0.63+ 0.87+ 0.12+ 0.16+ 0.22+ 0.32+ 0.46+ 0.69+ O.ll+ 0.18+ 0.31+ 0.59+ 0.12+ 0.85 +
d 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 3 3 3 3 6 10
0.89+ 0.86+ 0.83+ 0.81+ 0.79+ 0.78+ 0.78+ 0.73+ 0.80+ 0.84+ 0.90+ O.lO+ 0.13+ 0.17+ 0.21+ 0.27+ 0.34+ 0.42+ 0.54+ 0.69+ 0.91+ 0.12+ 0.17+ 0.25+ 0.39+ 0.64+ 0.52+ 0.18+
kJ
1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 5 10
2.9023 2.8386 2.7748 2.7107 2.6462 2.5814 2.5161 2.4501 2.3835 2.3159 2.2471 2.1767 2.1053 2.034 1 1.9637 1.8941 1.8260 1.7583 1.6912 1.6249 1.5593 1.4946 1.4309 1.3684 1.3073 1.2480 0.99459 0.90386
A CJ 1.2177 I.2318 1.2467 1.2626 1.2795 1.2975 1.3167 1.3370 1.3582 1.3803 1.4053 1.4375 1.4827 1.5426 1.6166 1.7839 1.8005 1.9419 2.05 14 2.1641 2.3180 2.4956 2.7033 2.9507 3.2509 3.6230 8.5308 65.423
kJ
OCJ
XC.1
PCJ
0.59481 0.56918 0.54326 0.51709 0.49068 0.46398 0.43690 0.4092 1 0.38048 0.35043 0.3’978 0.29058 0.26732 0.25396 0.24941 0.27326 0.25583 0.26834 0.27076 0.27224 0.27763 0.28304 0.28858 0.29449 0.30113 0.30903 0.41382 0.70323
0.10265 0.10869 0.11473 0.12066 0.12634 0.13154 0.13593 0.13899 0.13996 0.13765 0.13075 0 11866 0.10393 0.09154 0.08306 0.08077 0.073 1I 0.06928 0.0642 1 0.05923 0.05409 0.04877 0.04333 0.03784 0.03238 0.02703 0.00504 0.74 - 4
0.23 + 1 0.23 + 1 0.24 + 1 0.24 + 1 0.24 + 1 0.25 + 1 0.26 + 1 0.27 + 1 0.29 + i 0.32 + 1 0.36 + 1 0.43 + 1 0.57 + 1 0.78 + 1 0.11+2 0.14+2 0.18+2 0.23 + 2 0.30 + 2 0.39 + 2 0.53 + 2 0.73 + 2 0.11+3 0.16+3 0.25 + 3 0.43 + 3 0.41+ 5 0.15+ 10
0.26 + 1 0.24 + 1 0.22 + 1 0.20 + 1 0.18+ 1 0.17+ 1 0.15+ 1 0.14+ 1 0.12+ 1 0.11+1 0.95 0.82 0.71 0.62 0.54 0.48 0.42 0.37 0.33 0.29 0.25 0.22 0.19 0.16 0.13 0.11 0.29 - 1 0.21 - 2
437
W. Zhang, I. Gladweli / Bifurcation phenomena
AIN SOLUTtON
Legend CURVE 1
MAIN CURVE
--0
INVERSE LWMETER
l/d
Fig. 2. Solutions in D-l/d
k
POINT 2 ---
plane.
MAIN SOLWON
UPPERBOUNDARY II.
SINGULARITY
2 p % 0 iii Q
1.5-
LOWER BOUNDARY Legend MAIN CURVE
1 0
100
200
I 300
I 400
I 500
SQLUTICN WEQICN I I I 300 700 300
RADIUS OF CURVATURE Rs
Fig. 3. Solutions in B- R, plane.
CURVE 1 CURVE 2
438
W. Zhang I. Gladwell / Bifurcation phenomena
CO5NBF for D < 1.1. The technical detail is given in Section 2.3. In Table 1 we list the results for 0.2 < D < 3.5. For 0 < D < 0.2 limited machine precision prevents success in the ODE solver because of large solution values. There is no solution to the initial value u,(D) when D = 0 due to violation of physical constraints. For 2.2 < D G 3.5, the range where the solutions were computed successfully in H, our computational results agree to the fifth digit. For values D i 3.5 OUT numerical results (not presented) are similarly accurdte. The point D = 2.2 (referred as the connecting point in Fig. 2) is where the numerical results of [5] stop. As shown in [5] there is a turning point of d. It occurs at D = 2.914 and the corresponding minimum diameter is d, = 7.768. We call the curve above this turning point the upper branch and the one below it the lower branch. It is shown in [“I that the upper branch ends at the ideal detonation velocity D, = 6.1 where d + +m and AcJ, A,, --) 1. By computing the corresponding R, from (1.6), we observe that at D = 6.0, R, = 910.8 and it decreases monotonically as D decays until D = 2.7 where R, = 9.033, which is the minimum. R, increases monotonically as D further decreases. The behavior of R, is similar to that of d but differs from the R, of the simplified problem [3,13]. This difference contributes to the different solution phenomena discussed below. From our intermediate output we observe that for all D 2 2.1, variables u, A, A, m, x and v increase monotonically as t + T,, while p decreases monotonically. On the lower branch, as D tends to zero, R,, d, T,, and xc-- + +*, A,, + 1 and A, increases monotonically as shown in Table 1. Physically xcJ --) +a implies the reaction zone has an infinite length, which is physically unattainable. We illustrate their behavior since the information helps in finding what happens and where the “nonphysical” phenomenon starts. For D < 2.0, variables u and m grow fast near t = T,,. As D gets smaller u becomes flat over almost the whole integration interval and then climbs rapidly to its CJ value. This implies that the second CJ condition (1.23b) is almost satisfied over the whole integration interval and there is a boundary layer at t = T,, where the first CJ condition is also satisfied. The boundary effect becomes more severe as D + 0 revealing the possibility of instability. A similar trend also appears in variables A, A, cc)and p. Note ucJ > D when D G 1.65. In contrast, ucJ
W. Zhang, I. Gladwell / Bifurcation phenomena
439
Table 2 Numerical results on the upper boundary CUB) of the continuous solution Region ZII D
Rs
TCJ
d-
kJ
ACJ
*CJ
%J
-kJ
PCJ
2.435 2.400 2.300 2.200 2.100 2.000 1.900 1.800 1.700 1.600 1.500 1.400 1.300 1.200 1.100
0.28 + 2 0.38 + 2 0.55 + 2 0.78 + 2 0.11+3 0.17+3 0.27 + 3 0.43 + 3 0.72 + 3 0.12+4 0.21+ 4 0.39+4 0.74 + 4 0.15 +5 0.30 + 5
0.20 + 1 0.20 + 1 0.25 + 1 0.37 + 1 0.55 + 1 0.86 + 1 0.14+2 0.23 + 2 0.40 + 2 0.70 + 2 0.13+3 0.25 + 3 0.50 + 3 0.10+4 0.23 + 4
0.19+2 0.24 + 2 0.33 + 2 0.47 + 2 0.69 + 2 0.10+3 0.16+3 0.26 + 3 0.43 + 3 0.73 + 3 0.13+4 0.23 + 4 0.44 + 4 0.87 + 4 0.18+5
2.1864 2.1588 2.0831 2.0066 1.9288 1.8496 1.7689 1.6867 1.6029 1.5176 1.4308 1.3425 1.2527 1.1616 1.0691
1.0801 1.0532 1.0418 1.0375 1.0341 1.0312 1.0286 1.0261 1.0237 1.0213 1.0131 1.0170 1.0151 1.0133 1.0116
0.10958 0.09495 0.07968 0.06861 0.05917 0.05100 0.04384 0.03751 0.03194 0.02”34 0.02;!75 0.0151;3 0.01574 0.01291 0.01046
0.16- 1 O.lO- 1 0.60 - 2 0.37 - 2 0.23 - 2 0.13 - 2 0.74 - 3 0.40 - 3 0.21- 3 0.11-3 0.51-4 0.24 - 4 0.10-4 0.43 - 5 0.17-5
0.40 + 1 0.38 + 1 0.49 + 1 0.68 + 1 0.99 + 1 0.15+2 0.23 + 2 0.37 + 2 0.60 + 2 0.10-t-3 0.18+3 0.32 + 3 0.60 + 3 0.12+4 0.24 + 4
0.70 0.66 0.56 0.48 0.40 0.34 0.28 0.23 0.18 0.15 0.11 0.89 0.68 0.510.37 -
-
1 1 1 1
what happens on the main solution curve as D + Dr. We observe no boundary layer effect as D + D, on the main solution curve and UB. From Table 2 ucJ
Rs 0.26 + 2
L o-21+ 1
d
2.435
0.18+2
4T.J 2.1869
2.400 2.300 2.200 2.100 2.000 1.900 1.800 1.700 1.600 1.500 1.400 1.300 1.200 1.100
0.27 + 2 0.33 + 2 0.42 + 2 0.54+2 0.76 + 2 0.11+3 0.17+3 0.26 + 3 0.42 + 3 0.71+3 0.12+4 ‘X24+4 0.48 + 4 0.95 + 4
0.31+ 1 0.66 + 1 0.12+2 0.21+ 2 0.38 + 2 0.73 + 2 0.15+3 0.31-t 3 0.69 + 3 o-17+4 0.44 + 4 0.13+5 0.44 + 5 0.14+6
0.21-k 2 0.32 + 2 0.49 + 2 0.74 + 2 0.12+3 0.20 + 3 0.36 + 3 0.69 + 3 0.14+4 0.31+4 0.75 + 4 0.20 + 5 0.62 + 5 0.18+6
2.1622 2.0900 2.0171 1.9436 1.8689 1.7930 1.7161 1.6382 1.5591 1.4790 1.3978 1.3158 1.2332 1.1531
&I 1.0897 1.1270 1.2124 1.2963 1.3854 1.4806 1.5920 1.7295 1.9036 2.1296 2.4305 2.8430 3.4289 4.2985 5.4641
*cJ 0.11303 0.11800 0.12346 0.12603 0.12792 0.12697 0.12497 0.12295 0.12106 0.11939 0.11800 0.11699 0.11648 0.11664 0.11793
*cJ 0.18- 1 0.18- 1 0.17 - 1 0.15 - 1 0.14- 1 0.12- 1 0.98 - 2 0.79 - 2 0.63 - 2 0.50 - 2 0.39 - 2 0.31- 2 0.24 - 2 0.18-2 0.15-2
kJ 0.42 + 1 0.62 + 1 0.13+2 0.23 + 2 0.38 + 2 0.67 + 2 0.12+3 0.23 + 3 0.47 + 3 0.10+4 0.23 + 4 0.57+4 0.16+5 0.50 + 5 0.15+6
PCJ 0.70 0.67 0.58 0.50 0.43 0.36 0.31 0.25 0.21 0.17 0.14 0.11 0.81- 1 0.60 - 1 0.44 - 1
Wt’.Zhang, I. Gladwell / B4jkrcation phenomena
440 10-J
o-
I
I I I I I
8t-
I I
B-
F
Legend
I 5-
0
I I
4-
z
3 I 31 1 25
CURVE 1 CURVE 2
I I I
A
CURVE 3
+
CURVE 4
X
CURVE 5
0
CURVE 6
+
CURVE 7 ---
X
CURVE 8 --c-m
C’ m-e CURVE 9
I
I
I
I
I
I
30
3s
;I:
45
50
55
Fig. 4. Solutions in TcJ -R, plane.
20
1
18 1 ,
I
i
I I I I
I I I
I \
i
I I I I I i
Legend
I
I
0
I I I
A CURVE 3
b
CURVE 2
+
CURVE 4
X
CURVE 5
c3 CURVE 6 -I- CURVE 7 --_ X
2 1 25
I 30
I 35
I ii
I 46
Fig. 5. Solutions in Xc. - Is, plane.
I so
CURVE 1
I 55
CURVE 8 s---w
Cl CURVE 9 ---
441
Legend L!ts D
UB MC
Id0
150
200
2Ijo
400
450
560
W
Fig. 6. Solutions in D- Xcr plane.
from Region III in the (R,, D) and (l/d, D) planes. By connecting the main solution curve to UB, we would obtain a stable isolated solution curve for all D E (0, D,) which exhibi& no boundary layer effect and no violation of ucJ < D. Our limited numerical and graphical experiments seem to confirm this behavior. A schematic ( xCJt D) graph is given in Fig. 6. !n the (xc.D 0) plane, UB is entirely separated from the rest of solutions in III. If we map Region III to the (l/d, D) plane it is much narrower. Since for each pair (R,, E)) there is only one possible solution, just one value of d can be obtained from (1.6). But the converse does not necessarily hold due to the presence of xc3 in (1.6). Our computations show that in Region III solutions with different R, for a fixed D may result in the same diameter d. We use a bold curve on the (l/d, D) plane to indicate this solution region in Fig. 2. On this bold curve more than one solution may exist, I-Iowever only one solution is isolated and “stable”. 2.3. Regiiunal properties and their application Since d does no occur in (1.20)-(1.23), we concentrated on the (R,, D) plane. We performed computations in Regions I and II shown in Fig. 3. The shooting strategy employee (Section 2.1) to solve the singular ODE BVP allows us to perform the integration at any point in the (R,, D) plane. We simply pick a pair (R,, D) and integrate from t - 0. We let the integration proceed as far as it can and monitor the residues of the two CJ conditions (1.23). In ILegion I we find that after continuing the integration for a sufficiently long time the second CJ condition (1.23b) is satisfied and all variables except x tend to constant values. As a consequence the residue of the first CJ condition 41.23a) also tends to a nonvanishing constant. So the c:ondition is, never satisfied in Region I. In other words, as the integration is advanced, the ntunerator in (1.20a)
442
W: Zhang, 1. Gladwell / Bifurcation phenomena
tends to zero before the denominator vanishes. Then the whole system reaches a pseudo steady state and there is no further opportunity for the denominator to decay to zero. Mence there is no solution to the BVP in this region. Since all quantities but x tend to constants for a large t, we call this a steady state region, see Fig. 3. In Region II, we find the integration is always interrupted ivy the singularity and the second CJ condition is never satisfied. Thus the denominator rather than the numerator in (1.20a) reaches zero first. An QDE solver cannot integrate across such a singularity. Again there is no solution to the BVP. We call this a sing&r region (see Fig. 3). When this work was completed the authors learnt that Lynch INI had independently observed the regional properties and used them in numerical experiments. The main solution curve in the (&, D) plane is the boundary between the singular and steady state regions. On the curve both the numerator and the denominator in (1.20al tend to zero simultaneously. The existence of this solution curve seems natural after the identification of Regions I and II. In our computations we observe the immediate switch from one region to another as R, changes. This suggests to us that the value R, which yields a solution to the BVP is isolated. An alternative method for locating the main solution curve is to replace COSNBF by a bisection method. For a fixed D we locate an interval of R, with one end in the steady state region and the other in the singular region. Then we apply bisection based on the two different regional properties. This technique is used in our computations for D Q 1.0 where conventional continuation becomes expensive. Region III is surrounded by the singular Region II. In Region III we find that not only are the two CJ conditions satisfied simultaneously but also the integration can be continued beyond TcJ. Although mathematically the singularity always exists in (1,20a), numerically it may not be detected. As the numerator in (1.20a) tends to zero faster than the denominator (or at the same order), u& is finite. Therefore the integrator does not “see” the singularity and can continue the integration beyond it. When this happens the nonlinear equation solver adjusts the length of the integration interval back to T,, in order to match the two CJ conditions. Computationally this phenomenon is often observed deep in the interior of Region III. Since inside this region there exists a solution for each CR,, D), any R, selected by the nonlinear equation solver is an exact solution. The boundary of III is a fold [8]. In general, more work is required to trace a fold since to define it usually extra equations are added to the original system. In our case we can exploit the regional properties and locate the fold using an amount of work sinGBar to just finding a solution. It turns out that the continuation method succeeds in computing solutions on UB but not on LB. This is due to the solution (R,, T,,) being isolated on UB, which pcevent the nonlinear ell,uation solver searching for a solution inside Region III. Since there is no boundary layer effect on UB, a relatively large continuation step size can be used. When continuation starts on LB, the nonlinear equation solver finds a solution inside Region III rather than traces LB and it finds different solutions when using different continuation steps. So, the bisection method is appropriate for locating LB. 2.4. Stability of the continuum of solution phenomenon We compare the results to those for the simplified problem [4,13]. For large values of D they are qualitatively same. There is a turning point and an upper solution curve in the (l/d, 6))
W. Zhang, I. Gladwell / Bifurcation phenomena
2.8
P
443
MAIN SOLUTION
2.6
z G
0 d
>
II.
2.4
SINGULARITY
I
UPPERBOUNDARY
1.81
Legend CURVE 1
1.6
0
CURVE 2
RADIUS OF CURVATURE Rs
Fig. 7. Solutions in D-R, plane.
plane which ends at the ideal detonation point Dr. For lower values of D the two sets of results are different. For the simplified problem, [!3], the solution curve ends at the lowest critical value of D =D,_( > 0). And as D -+D,, d + +q xcJ + +m and A,, + 1. There is no other solution curve/region in the range (D,, Di). For the current problem there is no value DL( > 0) at which the solution curve/region ends. As Ei + 0, h, + 0 on LJB while A, > 0.7 in [4,13]. This suggests that the rate of the increase in A may be an important factor in determining the bifurcation diagram for lower values of D. 510we tested problem (1.20)-(1.23) with perturbations to (1.20~). First we varied A’ in (1.20~) by inserting a multiplier m, A’ = in@-
A)F(A, p).
(2 .2)
For m = 1.2, 1.5 and 2.0 we computed the solution with the same strategy and error tolerances as in Section 2.1. Our numerical results show that as m increases, Region III in the (R,, D) plane moves closer to the main solution curve. At both m = 1.5 and m = 2, III is connected with the main solution curve as shown in Fig. 7 for m = 2. Region III, instead of being surrounded by II, lies between I and II. When we map to the (l/d, D) plane, Region III connects directly to the main solution curve as shown in Fig. 8. Our second change arises from observing the different dependences on p in F. In [4,13], F is modeled by linear dependence on p and the numerical results always show pcJ 3 1. In the present problem F depends linearly on p for p large but is 0( p4) for p small, Cl.11). The switching point is close to p = 1. The two pieces are smoothly connected. Our computations show that pcJ drops below 1 for lower values of D, so the Q(p4) term yields a much smaller value of F and hence a small A’. We changed F by reducing the power of p to cubic and
CK Zhang, I. Gladwell / Bi&rcation phenomena
444
MAIN SOLUTION CURVE
3 1
Legend CURVE Cl REGION
PNVERSE DIAMETER Fig. 8. Solutions
lid
in D-1 /d plane.
quadratic respectively in the expression for pp, (1.1 lb) whilst preserving the smooth connection to the linear piece [14]. Then we computed the solutions to the perturbed prob!em with the same strategy and error tolerances. We obtained the same solution pattern as in Fig. 7. Region III connects to the main solution curve and lies between I and II. Though the solution values are different, the pattern and the behavior are the same.
We have transformed a parameterized index one DAE BVP with index two on the free boundary into a parameterized singular ODE BVP with a free boundary. By using a shooting technique for solving the singular free boundary value problem, we successfully computed the sob on to the detonation problem studied in [5,9] and presented and justified the complete bifurcation diagram. Because of the difficulty of treating a singularity on the free boundary, we identified regional properties in parameter space. These properties were used in locating the isolated solution and the boundary of a contiwum of solutions, which turned out to be computationally efficient. Bifurcation phenomena of related detonation problems were also presented tikh demonstrates that the phenomena persist.
eferences [l] U. Ascher, Private communication. [2] J.B. Bdzil, Steady s!:?!c. two-dimensional
detonation,
J. Flrrid Me&. 1198(1981) 195-226.
W. Zhang, I. Gladwell / Bifurcation phenomena
445
[3] CC. Beardah and R.M. Thomas, Techniques for the numerical solution of a singular index one differential/
algebraic equation boundary value problem arising in the modelling of unconfined detonations, Numerical Analysis Report, No. 191, Department of Mathematics, UMIST (1990). [4] M. Braithwaite, T. Farran, I. Gladwell, P.M. Lynch, A. Minchinton, I.B. Parker and R.M. Thomas, A detonation problem posed as a differential algebraic boundary value problem, Math. Engrg. Ind. 3 (1) (1?90) 45-57. [S] M. Braithwaite, T. Farran, I. Gladwell, P.M. Lynch, A. Minchinton, I.B. Parker and R.M. Thomas, Private
communication. [6] K.E. Brenan, Sk. Campbell and L.R. Petzold, Numerical Solution of Initial-Value Problems in Differential Algebraic Equations (Elsevier, New York, 1989). [7] W. Fickett and W.C. Davis, Detonation (University of California Press, Berkeley, CA, 1979). [S] H.B. Keller, Numerical Methods in Bifurcation Problems (Springer, Heidelberg, 1987). [9] I. Kirby and G. Leiper, A small divergent detonation theory for intermolecular explosives, in: Proceedings 8th International Symposium on Detonation, Albuquerque, NM (1985). [lo] I. Kirby, G. Leiper and A. Hackett, Determination of reaction rates in intermolecular explosives using the electromagnetic particle velocity gauge, in: Proceedings 8th International Symposium on Detonation, Albuquerque, NM (1985). [ 1l] P.M. Lynch, Private communication. [12] The NAG Fortran Library Manual - Mark 13 (1st ed., 1988). [13] W. Zhang and I. Gladwell, Analysis of a simplified detonation problem, Math. Engrg. Ind. (in press). [14] W. Zhang, and I. Gladwell, Bifurcation phenomena arising in a singular differential algebraic equation detonation problem, Report No. 90-14, Department of Mathematics, Southern Methodist University, Dallas, TX (1990).