Failed Convergence in Detonation Computations

[Andrew Mizener]

Good afternoon, all.

I am attempting to compute the H2-O2 detonation properties for a range of equivalence ratios (0.5 ≤ ϕ ≤ 2.0). For most of these values, I am getting no trouble, but for some of them I'm getting a range of types of failed convergences. Equivalence ratio is the only variable; upstream temperature and pressure are functionally constant*. I am using the CJspeed and PostShock_eq routines and the h2o2_highT mechanism from the Shock & Detonation Toolbox, and Cantera 2.2.1 in MATLAB.

I'm seeing three types of failed convergence behavior. Using sample equivalence ratios:

1.  At ϕ = 0.518, the routine enters what appears to be an infinite regression somewhere inside CJspeed. (At least, close enough to infinite - after letting it sit for a half-hour and seeing no progress, I stopped the code.)

2.  At ϕ = 0.660, the value of U_CJ computed by CJspeed (1588.1677 m/s) is approximately 40% too low, which causes PostShock_eq to fail to converge, and I get an error message stating "shk_eq_calc did not converge for U = 1588.1677", but the script continues to run.

3.  At ϕ = 0.655, at some point in the computation of what I believe is the equation of state for the (upstream?) conditions, the iterated value of density goes negative, and the script crashes. The error message for this is given below**.

It is worth mentioning that ϕ =0.650 and ϕ = 0.665 both work fine. A friend suggested I try changing the solver used in the "equilibrate" routine from 0 to 1, to use the Gibbs-free-energy-minimization method [i.e. changing equilibrate(gas, 'TV') to equilibrate(gas, 'TV', 1)], but this did not solve the problem, and instead just moved around the points where convergence fails.

I have no idea why some conditions work and some do not. Does anyone have any insight into how I can correct this? I have attached a minimal example code, as requested. It is based on my full code, the relevant part of which is in turn based on the "demo_CJstate.m" example from the Shock and Detonation Toolbox.

Thank you very much.

-- Andrew Mizener

* The upstream conditions are actually choked flow through a nozzle, and the insentropic critical pressures and temperatures for which are a function of the ratio of specific heats, which varies slightly with equivalence ratio. However, the difference between - for example - ϕ = 0.65 and ϕ = 0.655 is on the order of 1e-3 percent in both pressure and temperature.

** Here is the error message given for example 3:

Error using ThermoPhase/setDensity (line 12)

The density must be positive.

Error in ThermoPhase/set (line 169)

        setDensity(tp, 1.0/vval);

Error in eq_state (line 22)

    set(gas, 'Density', r1, 'T', T1);

Error in CJ_calc (line 37)

[P, H] = eq_state(gas,r,T);

Error in CJspeed (line 55)

        [gas, w1(i)] = CJ_calc(gas, gas1, ERRFT, ERRFV, x);

Error in Det_Min_Example (line 136)

    [U_CJ] = CJspeed(P_I, T_I, Mixture, Mechanism, fig_num);

[Ray Speth]

Hi Andrew,

Using the code you provided, Cantera 2.2.1 and a stock version of the SD Toolbox (the version updated for Cantera 2.1), I was unable to replicate problems (2) and (3) that you described. The script you provided gives the output:

phi = 0.65

     P_I  = 267576.203 Pa

     T_I  = 208.2677 K

     Mix  = H2:1.3 O2:1

     U_CJ = 2566.1023 m/s

phi = 0.655

     P_I  = 267574.3425 Pa

     T_I  = 208.2658 K

     Mix  = H2:1.31 O2:1

     U_CJ = 2571.9558 m/s

phi = 0.66

     P_I  = 267572.4978 Pa

     T_I  = 208.2639 K

     Mix  = H2:1.32 O2:1

     U_CJ = 2577.7824 m/s

Which seems like it is OK to me. Are there changes that you've made to your copy of the SD Toolbox, perhaps?

For problem (1), the issue seems to be the (somewhat strange) minimization algorithm used in CJspeed, which is sensitive to the numerical properties of the inner iteration (CJ_calc). My first suggestion would be to adjust the tolerances used in the inner function (ERRFT and ERRFV) so that they are tighter than the tolerance specified for the outer iteration (g < 0.99999). This works for this specific case, but I don't know how it affects things more generally. I suspect that there are more robust and efficient minimization algorithms that could be used in CJspeed instead of this apparently homegrown method.

Regards,

Ray

[Andrew Mizener]

Update:

I was getting some assistance from a friend, who suggested I try a different mechanism. He uses the gri30_highT mechanism for H2-O2 calculations to good success, and so I tried it. Unfortunately, it too helps for some points but hurts in others - again, it just appears to just move around the places where the conversion fails. Using gri30_highT, I get modes of failure at these points (among others, these are sample values):

1. ϕ = 1.928 enters the infinite recursion. (But 1.9285 works fine.)

2. ϕ = 0.686 and ϕ = 0.698 give me the "shk_eq_calc did not converge for U = _____" error.

He also suggested that since the thermodynamic curvefits only go down to 200 K, that perhaps I should increase my initial plenum pressure (it was 250 K, which gives a detonation upstream temperature of just over 208 K). Increasing plenum pressure to 300 K results in an upstream temperature of just under 250 K. Again, all this did is to move around the points of non-convergence. Using 300 K, I got the third mode of failure at ϕ = 0.656.

Still stumped.

-- A

[Andrew Mizener]

Ray,

Thank you. I don't think I've made any edits to any of the SD Toolbox files, but just in case, I'll re-download them and try again.

I'll also try altering the convergene tolerance and report back.

-- Andrew

[Andrew Mizener]

Update #2: I reinstalled the SD Toolbox, and it worked for the first pass of results I needed to do (equivalence ratio as the only variable, 0.5 ≤ ϕ ≤ 2.0).

Then, I moved on to the second set of results, in which I generate contour plots of parameters as a function of two variables, one of which is ϕ and... it promptly threw an error (well, threw an error partway through a 15 x 15 matrix). The grid consists of evenly spaced points, and for comparison purposes, I'm using the same grid for all of my cases - hence the numbers not being pretty. I've attached another version of the minimum code example that is set up to run those two cases. 

Isolating the failure points (T0i = 300 K for all): 

• P0i = 10.5 atm, ϕ = 0.607142857142857, mechanism h2o2_highT, I get failure mode #3. The error message is:

Error using ThermoPhase/setDensity (line 12)

The density must be positive.

Error in ThermoPhase/set (line 169)

        setDensity(tp, 1.0/vval);

Error in eq_state (line 20)

    set(gas, 'Density', r1, 'T', T1);

Error in CJ_calc (line 37)

[P, H] = eq_state(gas,r,T);

Error in CJspeed (line 55)

        [gas, w1(i)] = CJ_calc(gas, gas1, ERRFT, ERRFV, x);

Error in Det_Min_Example_2 (line 114)

    [U_CJ] = CJspeed(P_I, T_I, Mixture, Mechanism, fig_num);

I tried running the same grid with the gri30 mechanism (which takes longer, but appears to converge in places where h2o2 does not), and I found a new partial convergence mode (#4!):

• P0i = 9.142857142857142 atm, ϕ = 0.607142857142857, mechanism h2o2_highT, I get the following output:

ans =

j = 500

After this appears a bunch of times (approx. 42) Then it the program moves on to the next data point. However it does not converge. I've attached the contour plot that is generated. The little blip that occurs (highlighted in the center figure) is at the point P0i = 9.142857142857142, ϕ = 0.607142857142857. (Interestingly the repeated "ans = / j = 500" output happens before the previous example as well.)

Is it just a case of "play with the iteration tolerances", or am I doing something more fundamental wrong here?

Thanks again for all your assistance.

-- Andrew

[Andrew Mizener]

... it would appear that I answered my own question. After posting, I realized that I hadn't tried Ray's suggestion of tightening up the error tolerances (ERRFT and ERRFV) in CJspeed for CJ_calc. I did so, and it appears that everything works now.

Sorry for asking a further dumb question without trying out all of the previous suggestions, and thank you again, Ray, for the advice.

-- Andrew