[KarenRei]

So, I have a mesh and general setup that works just fine with rhoCentralFoam and no chemistry. checkMesh only reports cell concavity errors (which are unavoidable for anything made with snappyHexMesh where the resolution varies - they occur at each resolution transition zone). Non-orthogonality is low, only 36.2 max, 6.08 average.

I'm now trying adding in chemistry (at this point, I'm pretty indifferent to the solver, I just want something that runs; I can work from there). But something curious is happening. It runs fine if I remove aluminum and all aluminum reactions from my reactions list. But if I add them in, it crashes after just 2 or 3 timesteps:

...
min/max(T) = 100, 4000
DICPCG: Solving for p, Initial residual = 0.15503, Final residual = 6.5712e-09, No Iterations 1
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 2.36489e-09, global = -8.20745e-10, cumulative = -3.75218e-08
rho max/min : 8.6484 0.0268814
DICPCG: Solving for p, Initial residual = 4.34903e-06, Final residual = 5.15751e-11, No Iterations 1
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0

time step continuity errors : sum local = 2.29547e-11, global = 8.04112e-13, cumulative = -3.7521e-08
rho max/min : 8.63582 0.0268811
ExecutionTime = 3.15 s ClockTime = 3 s

Courant Number mean: 0.000512311 max: 0.0978574
deltaT = 1.72786e-06
Time = 4.36775e-06

diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
PIMPLE: iteration 1
DILUPBiCGStab: Solving for C30H62_Triacotane, Initial residual = 0.00011374, Final residual = 6.65827e-09, No Iterations 1
DILUPBiCGStab: Solving for O2, Initial residual = 0.000291101, Final residual = 6.37422e-09, No Iterations 1
DILUPBiCGStab: Solving for CO, Initial residual = 0.000536741, Final residual = 1.04301e-07, No Iterations 1
DILUPBiCGStab: Solving for H2O, Initial residual = 0.000255655, Final residual = 2.31429e-09, No Iterations 1
DILUPBiCGStab: Solving for O, Initial residual = 0.00130935, Final residual = 1.39516e-07, No Iterations 1
DILUPBiCGStab: Solving for CO2, Initial residual = 0.0560067, Final residual = 3.14478e-08, No Iterations 1
DILUPBiCGStab: Solving for OH, Initial residual = 0.000779744, Final residual = 1.6847e-08, No Iterations 1
DILUPBiCGStab: Solving for H, Initial residual = 0.000946765, Final residual = 2.53462e-08, No Iterations 1
DILUPBiCGStab: Solving for N, Initial residual = 0.00185759, Final residual = 9.54203e-08, No Iterations 1
DILUPBiCGStab: Solving for NO, Initial residual = 0.00131261, Final residual = 1.39316e-07, No Iterations 1
DILUPBiCGStab: Solving for Al, Initial residual = 0.00011374, Final residual = 6.65827e-09, No Iterations 1
DILUPBiCGStab: Solving for AlH, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCGStab: Solving for AlO, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCGStab: Solving for AlOH, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCGStab: Solving for Al2O2, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCGStab: Solving for AlO2, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCGStab: Solving for Al2O3(l), Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCGStab: Solving for Al2O3, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCGStab: Solving for Al2O3(a), Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCGStab: Solving for AlH3, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCGStab: Solving for h, Initial residual = 0.00333788, Final residual = 5.67018e-09, No Iterations 1


--> FOAM FATAL ERROR:
Maximum number of iterations exceeded

in file /opt/OpenFOAM/OpenFOAM-dev/src/thermophysicalModels/specie/lnInclude/thermoI.H at line 66.

FOAM aborting

(I stripped out some spammy Janaf-bounds errors from that output). You can see that the courant number is quite low. Temperatures are in a reasonable range.

So, since it happens when aluminum is in the list of reactants - even if I don't list any reactions for it - my first thought was, there must be something wrong with my thermodynamics definition:

Al
{
specie
{
nMoles 1;
molWeight 26.9815385111;
}
elements
{
Al 1;
}
thermodynamics
{
Tlow 100.000000;
Thigh 6000.000000;
Tcommon 2790.812000;
highCpCoeffs ( 2.34862 0.000143612 -4.23667e-08 2.72036e-12 3.70085e-16 39104.1 6.37634 );
lowCpCoeffs ( 1.25994 0.00655861 -5.9248e-06 2.24926e-09 -3.04921e-13 242.739 -4.57142 );
}
transport
{
As 1.67212e-06;
Ts 170.672;
}
}

Indeed, if I replace it with my definition for N2, as an example, things run fine. So, that raises the question, what could be wrong with this? I checked out the curves it generates for entropy, enthalpy, and specific heat. They match up to the curves on Citrination quite well. There's no unexpected swings, even outside of the janaf bounds. And since I have no aluminum reactions currently included in my reactions file - simply aluminum listed as a reactant - there couldn't be some sort of other aluminum compound that's causing the problem.

So what could be wrong? What can I do to try to debug this? I've run out of ideas here...

In case it matters: it's a O2/N2 atmosphere (air), N2 listed as an inert gas, and there's a flowRateInletVelocity injecting a mixture of cryogenic aluminum dust, wax, and oxygen (all treated as a gas so that I don't have to track particles - I may migrate to a particle-based system later. I'm using EulerImplicit as the solver, but even if I change it to no solver / chemistry off the crash happens. Thermotype is:

type hePsiThermo;
mixture reactingMixture;
transport sutherland;
thermo janaf;
energy sensibleEnthalpy;
equationOfState perfectGas;
specie specie;

Runtime is adjustable with max courant 0.2.

[KarenRei]

Bout of minor inspiration... it's only the highCpCoeffs in aluminum that are causing the problem; I can replace just them, and it runs, but with them in place the max iter error occurs.
If I swap out the second to last janaf param for a lower one, it works too.

Real, crashes: highCpCoeffs ( 2.348620.000143612 -4.23667e-08 2.72036e-12 3.70085e-16 39104.1 6.37634 );
Altered, works: highCpCoeffs ( 2.348620.000143612 -4.23667e-08 2.72036e-12 3.70085e-16 -958.396 6.37634 );

Need to figure out what's going on here...

[KarenRei]

Update: Well, that's the enthalpy offset. And the reason for the big enthalpy offset is the heat of condensation of aluminum at Tcommon:

https://www.citrination.com/uploads/...tate-reference

I calculate T=2790 (low params) at entropy 84441 and T=2791 (high params) at entropy 382175, which are roughly what show in the graph (T=2791: 87670->381700). I'm calculating as:

R * (a[4]*math.pow(T, 5)/5 + a[3]*math.pow(T, 4)/4 + a[2]*math.pow(T, 3)/3 + a[1]*math.pow(T, 2)/2 + a[0]*T + a[5])

... which should be correct. But my best guess has to be that this problem relates to the large enthalpy difference at that point? The thing is, it's a *real* large jump in enthalpy there, I have to account for it....

[KarenRei]

I think the code that's failing is:

// return the temperature corresponding to the value of the
// thermodynamic property f, given the function f = F(T) and dF(T)/dT
template<class thermo>
inline scalar specieThermo<thermo>::T
(
scalar f,
scalar T0,
scalar (specieThermo<thermo>::*F)(const scalar) const,
scalar (specieThermo<thermo>::*dFdT)(const scalar) const
) const
{
scalar Test = T0;
scalar Tnew = T0;
scalar Ttol = T0*tol_;
int iter = 0;

do
{
Test = Tnew;
Tnew =
(this->*limit)
(Test - ((this->*F)(p, Test) - f)/(this->*dFdT)(p, Test));

if (iter++ > maxIter_)
{
FatalErrorIn
(
"specieThermo<thermo>::T(scalar f, scalar T0, "
"scalar (specieThermo<thermo>::*F)(const scalar) const, "
"scalar (specieThermo<thermo>::*dFdT)(const scalar) const"
") const"
) << "Maximum number of iterations exceeded"
<< abort(FatalError);
}

} while (mag(Tnew - Test) > Ttol);

return Tnew;
}

.... I don't see how to increase "maxIter" (tried setting it in fvSolutions, no matter what value I tried it didn't seem to make a difference). I also tried changing retTol in fvSolutions, no difference. I'm thinking that maybe these aren't the values it's looking at here?

My parsing of what it's doing: it's trying to find what temperature a chemical species is going to end up at... which depends on the specific heat differences en route to the new temperature, aka the net enthalpy change. "Test" is the previous temperature it tried and TNew is the new attempt, which is based on trying to follow the slope of some function (specific heat?) in some manner.... I'm going to instrument the code to try to find out.

[KarenRei]

The results of instrumenting thermophysicalModels/specie/lnInclude/thermoI.H:

T(f=906725, p=62165.7, T0=1035.48)
tol=0.0001, Ttol = 0.103548
Test=1035.48, Tnew=875.878 (1.10861e+06, 201886, 1264.9, 875.878)
Test=875.878, Tnew=1053.82 (684624, -222101, 1248.13, 1053.82)
Test=1053.82, Tnew=876.255 (1.13184e+06, 225111, 1267.73, 876.255)
Test=876.255, Tnew=1053.81 (685095, -221630, 1248.2, 1053.81)
Test=1053.81, Tnew=876.254 (1.13182e+06, 225097, 1267.73, 876.254)
Test=876.254, Tnew=1053.81 (685094, -221631, 1248.2, 1053.81)
Test=1053.81, Tnew=876.254 (1.13182e+06, 225097, 1267.73, 876.254)
Test=876.254, Tnew=1053.81 (685094, -221631, 1248.2, 1053.81)
Test=1053.81, Tnew=876.254 (1.13182e+06, 225097, 1267.73, 876.254)
Test=876.254, Tnew=1053.81 (685094, -221631, 1248.2, 1053.81)
Test=1053.81, Tnew=876.254 (1.13182e+06, 225097, 1267.73, 876.254)
Test=876.254, Tnew=1053.81 (685094, -221631, 1248.2, 1053.81)
Test=1053.81, Tnew=876.254 (1.13182e+06, 225097, 1267.73, 876.254)
Test=876.254, Tnew=1053.81 (685094, -221631, 1248.2, 1053.81)
Test=1053.81, Tnew=876.254 (1.13182e+06, 225097, 1267.73, 876.254)
Test=876.254, Tnew=1053.81 (685094, -221631, 1248.2, 1053.81)
Test=1053.81, Tnew=876.254 (1.13182e+06, 225097, 1267.73, 876.254)
Test=876.254, Tnew=1053.81 (685094, -221631, 1248.2, 1053.81)
Test=1053.81, Tnew=876.254 (1.13182e+06, 225097, 1267.73, 876.254)
Test=876.254, Tnew=1053.81 (685094, -221631, 1248.2, 1053.81)
Test=1053.81, Tnew=876.254 (1.13182e+06, 225097, 1267.73, 876.254)
Test=876.254, Tnew=1053.81 (685094, -221631, 1248.2, 1053.81)
Test=1053.81, Tnew=876.254 (1.13182e+06, 225097, 1267.73, 876.254)
Test=876.254, Tnew=1053.81 (685094, -221631, 1248.2, 1053.81)
Test=1053.81, Tnew=876.254 (1.13182e+06, 225097, 1267.73, 876.254)
Test=876.254, Tnew=1053.81 (685094, -221631, 1248.2, 1053.81)
Test=1053.81, Tnew=876.254 (1.13182e+06, 225097, 1267.73, 876.254)
Test=876.254, Tnew=1053.81 (685094, -221631, 1248.2, 1053.81)
Test=1053.81, Tnew=876.254 (1.13182e+06, 225097, 1267.73, 876.254)
Test=876.254, Tnew=1053.81 (685094, -221631, 1248.2, 1053.81)
Test=1053.81, Tnew=876.254 (1.13182e+06, 225097, 1267.73, 876.254)
Test=876.254, Tnew=1053.81 (685094, -221631, 1248.2, 1053.81)
Test=1053.81, Tnew=876.254 (1.13182e+06, 225097, 1267.73, 876.254)
Test=876.254, Tnew=1053.81 (685094, -221631, 1248.2, 1053.81)
Test=1053.81, Tnew=876.254 (1.13182e+06, 225097, 1267.73, 876.254)
Test=876.254, Tnew=1053.81 (685094, -221631, 1248.2, 1053.81)
Test=1053.81, Tnew=876.254 (1.13182e+06, 225097, 1267.73, 876.254)
Test=876.254, Tnew=1053.81 (685094, -221631, 1248.2, 1053.81)
Test=1053.81, Tnew=876.254 (1.13182e+06, 225097, 1267.73, 876.254)
Test=876.254, Tnew=1053.81 (685094, -221631, 1248.2, 1053.81)
Test=1053.81, Tnew=876.254 (1.13182e+06, 225097, 1267.73, 876.254)
Test=876.254, Tnew=1053.81 (685094, -221631, 1248.2, 1053.81)
Test=1053.81, Tnew=876.254 (1.13182e+06, 225097, 1267.73, 876.254)
Test=876.254, Tnew=1053.81 (685094, -221631, 1248.2, 1053.81)
Test=1053.81, Tnew=876.254 (1.13182e+06, 225097, 1267.73, 876.254)
Test=876.254, Tnew=1053.81 (685094, -221631, 1248.2, 1053.81)
Test=1053.81, Tnew=876.254 (1.13182e+06, 225097, 1267.73, 876.254)
Test=876.254, Tnew=1053.81 (685094, -221631, 1248.2, 1053.81)
Test=1053.81, Tnew=876.254 (1.13182e+06, 225097, 1267.73, 876.254)
Test=876.254, Tnew=1053.81 (685094, -221631, 1248.2, 1053.81)
Test=1053.81, Tnew=876.254 (1.13182e+06, 225097, 1267.73, 876.254)
Test=876.254, Tnew=1053.81 (685094, -221631, 1248.2, 1053.81)
Test=1053.81, Tnew=876.254 (1.13182e+06, 225097, 1267.73, 876.254)
Test=876.254, Tnew=1053.81 (685094, -221631, 1248.2, 1053.81)
Test=1053.81, Tnew=876.254 (1.13182e+06, 225097, 1267.73, 876.254)
Test=876.254, Tnew=1053.81 (685094, -221631, 1248.2, 1053.81)
Test=1053.81, Tnew=876.254 (1.13182e+06, 225097, 1267.73, 876.254)
Test=876.254, Tnew=1053.81 (685094, -221631, 1248.2, 1053.81)
Test=1053.81, Tnew=876.254 (1.13182e+06, 225097, 1267.73, 876.254)
Test=876.254, Tnew=1053.81 (685094, -221631, 1248.2, 1053.81)
Test=1053.81, Tnew=876.254 (1.13182e+06, 225097, 1267.73, 876.254)
Test=876.254, Tnew=1053.81 (685094, -221631, 1248.2, 1053.81)
Test=1053.81, Tnew=876.254 (1.13182e+06, 225097, 1267.73, 876.254)
Test=876.254, Tnew=1053.81 (685094, -221631, 1248.2, 1053.81)
Test=1053.81, Tnew=876.254 (1.13182e+06, 225097, 1267.73, 876.254)
Test=876.254, Tnew=1053.81 (685094, -221631, 1248.2, 1053.81)
Test=1053.81, Tnew=876.254 (1.13182e+06, 225097, 1267.73, 876.254)
Test=876.254, Tnew=1053.81 (685094, -221631, 1248.2, 1053.81)
Test=1053.81, Tnew=876.254 (1.13182e+06, 225097, 1267.73, 876.254)
Test=876.254, Tnew=1053.81 (685094, -221631, 1248.2, 1053.81)
Test=1053.81, Tnew=876.254 (1.13182e+06, 225097, 1267.73, 876.254)
Test=876.254, Tnew=1053.81 (685094, -221631, 1248.2, 1053.81)
Test=1053.81, Tnew=876.254 (1.13182e+06, 225097, 1267.73, 876.254)
Test=876.254, Tnew=1053.81 (685094, -221631, 1248.2, 1053.81)
Test=1053.81, Tnew=876.254 (1.13182e+06, 225097, 1267.73, 876.254)
Test=876.254, Tnew=1053.81 (685094, -221631, 1248.2, 1053.81)
Test=1053.81, Tnew=876.254 (1.13182e+06, 225097, 1267.73, 876.254)
Test=876.254, Tnew=1053.81 (685094, -221631, 1248.2, 1053.81)
Test=1053.81, Tnew=876.254 (1.13182e+06, 225097, 1267.73, 876.254)
Test=876.254, Tnew=1053.81 (685094, -221631, 1248.2, 1053.81)
Test=1053.81, Tnew=876.254 (1.13182e+06, 225097, 1267.73, 876.254)
Test=876.254, Tnew=1053.81 (685094, -221631, 1248.2, 1053.81)
Test=1053.81, Tnew=876.254 (1.13182e+06, 225097, 1267.73, 876.254)
Test=876.254, Tnew=1053.81 (685094, -221631, 1248.2, 1053.81)
Test=1053.81, Tnew=876.254 (1.13182e+06, 225097, 1267.73, 876.254)
Test=876.254, Tnew=1053.81 (685094, -221631, 1248.2, 1053.81)
Test=1053.81, Tnew=876.254 (1.13182e+06, 225097, 1267.73, 876.254)
Test=876.254, Tnew=1053.81 (685094, -221631, 1248.2, 1053.81)
Test=1053.81, Tnew=876.254 (1.13182e+06, 225097, 1267.73, 876.254)
Test=876.254, Tnew=1053.81 (685094, -221631, 1248.2, 1053.81)
Test=1053.81, Tnew=876.254 (1.13182e+06, 225097, 1267.73, 876.254)
Test=876.254, Tnew=1053.81 (685094, -221631, 1248.2, 1053.81)
Test=1053.81, Tnew=876.254 (1.13182e+06, 225097, 1267.73, 876.254)
Test=876.254, Tnew=1053.81 (685094, -221631, 1248.2, 1053.81)
Test=1053.81, Tnew=876.254 (1.13182e+06, 225097, 1267.73, 876.254)
Test=876.254, Tnew=1053.81 (685094, -221631, 1248.2, 1053.81)
Test=1053.81, Tnew=876.254 (1.13182e+06, 225097, 1267.73, 876.254)
Test=876.254, Tnew=1053.81 (685094, -221631, 1248.2, 1053.81)
Test=1053.81, Tnew=876.254 (1.13182e+06, 225097, 1267.73, 876.254)
Test=876.254, Tnew=1053.81 (685094, -221631, 1248.2, 1053.81)
Test=1053.81, Tnew=876.254 (1.13182e+06, 225097, 1267.73, 876.254)
Test=876.254, Tnew=1053.81 (685094, -221631, 1248.2, 1053.81)



--> FOAM FATAL ERROR:
Maximum number of iterations exceeded

in file /opt/OpenFOAM/OpenFOAM-dev/src/thermophysicalModels/specie/lnInclude/thermoI.H at line 76.

FOAM aborting

#0 Foam::error:rintStack(Foam::Ostream&) at ??:?
#1 Foam::error::abort() at ??:?
#5 ? at ??:?
#6 __libc_start_main in "/lib64/libc.so.6"
#7 ? at ??:?


So, to narrow this down. Tnew is:

Tnew =
(this->*limit)
(Test - ((this->*F)(p, Test) - f)/(this->*dFdT)(p, Test));

So, when we see:

Test=1035.48, Tnew=875.878 (1.10861e+06, 201886, 1264.9, 875.878)

it's saying:

875.878 = (this->*limit)(1035.48 - ((this->*F)(62165.7, 1035.48) - 906725) / (this->*dFdT)(62165.7, 1035.48))
875.878 = (this->*limit)(1035.48 - (1108611 - 906725) / 1264.9)
875.878 = (this->*limit)(1035.48 - 201886 / 1264.9)
875.878 = (this->*limit)(1035.48 - 159,59)
875.878 = (this->*limit)(875.878)

Then this runs it in reverse:

Test=875.878, Tnew=1053.82 (684624, -222101, 1248.13, 1053.82)

1053.82 = (this->*limit)(875.878 - ((this->*F)(62165.7, 875.878) - 906725) / (this->*dFdT)(62165.7, 875.878))
1053.82 = (this->*limit)(875.878 - (684624 - 906725) / 1264.9)
1053.82 = (this->*limit)(875.878 - -222101 / 1264.9)
1053.82 = (this->*limit)(875.878 + 175,59)
1053.82 = (this->*limit)(1053.82)

So it just cycles between these two points (~875 and ~1053). But why? Also: strange that this is *below Tcommon*. Tcommon is 2790.812.... Yet when I had previously tried swapping out the high Cp coeffs for aluminum with those from N2, it had worked...

Still not sure what function "*F" is supposed to be here... looks too high to be enthalpy J/K-mol (1.1e6 at 1035K?), and way too high for specific heat or entropy (J/mol).... Also, it's something that's pressure-dependent... hmm.... maybe it's not per-mole, maybe it's.... per... something else?

ED: If I swap out the low Cp coeffs on aluminum, the output is exactly the same. Huh??? Okay, so maybe we're not dealing with aluminum here.... even though when I change the thermo data for aluminum, or remove aluminum from the run altogether, it works? *mind boggles*

[KarenRei]

Okay, so some results of further tests.

*F is the function Hs. So, enthalpy.
*dFdT is the function Cp. So, specific heat. Which makes sense, that is indeed the derivative of enthalpy.

For p=62165.7, T=1035.48

Hs returns hs(p, T) / W(). hs is 26182800 and W is 23.6177, hence the 1108610
Cp returns cp(p, T) / W(). cp is 29873.9 and W is 23.6177, hence the 1264.9

Again, these don't look like H or Cp per mole, they must be per mass or per volume. cp() appears to be:

const coeffArray& a = coeffs(T);
return RR*
(
((((a[4]/5.0*T + a[3]/4.0)*T + a[2]/3.0)*T + a[1]/2.0)*T + a[0])*T
+ a[5]
)
+ EquationOfState::h(p, T);

... while hs is:

const coeffArray& a = coeffs(T);
return RR*
(
((((a[4]/5.0*T + a[3]/4.0)*T + a[2]/3.0)*T + a[1]/2.0)*T + a[0])*T
+ a[5]
)
+ EquationOfState::h(p, T) - hc();

In practice:

cp: 29940.7 + 0
hs: 2.61367e+07 + 0 - -594257

Sooooo the numbers actually are this large, straight from the janaf params....? Okay, need to figure out how *that's* happening.... enthalpy should only be on the order of 30k J/K-mol and specific heat on the order of 30 J/mol... Let+s instrument that better:

29940.7 = 8314.47*((((-5.02086e-15*1053.81 + 7.87828e-11)*1053.81 + -4.41075e-07)*1053.81 + 0.00112498)*1053.81 + 2.81934))
2.61367e+07 = 8314.47*((((-5.02086e-15/5.0*1053.81 + 7.87828e-11/4.0)*1053.81 + -4.41075e-07/3.0)*1053.81 + 0.00112498/2.0)*1053.81 + 2.81934)*1053.81 + -303.117)) + 0 - -594257
25542443 = 8314.47*((((-5.02086e-15/5.0*1053.81 + 7.87828e-11/4.0)*1053.81 + -4.41075e-07/3.0)*1053.81 + 0.00112498/2.0)*1053.81 + 2.81934)*1053.81 + -303.117))

Huh, wait a minute.... R is 8314.47? What sort of units are those? I've been using R=8.314, aka J/mol-K. That R has to be, what, cm3 kPa/K-mol or something weird? Then that raises the question.... should I be accepting that OpenFOAM is wanting cp values in these sort of unexpected units... or are they expecting the janaf coefficients to be expressed in some sort of weird units, and then converting them to J/mol-K? My parameters *must* be at least approximately right, certainly not off by an order of 1000, because they match pretty closely to janaf params in thermo files that come with the distributions (where matching species can be found). The coefficients should all be off by about three orders of magnitude if I was doing it wrong.

Another oddity: I don't recognize those janaf coefficients it lists! They're not anywhere in my thermo.compressibleGas file! Yet strangely the results aren't that far off for aluminum, if you scale the R by 1000.....

God, I am so baffled here....

Ed: I think trying to understand how OpenFOAM is representing the data internally is leading me astray. The real problem is, I think, that we have some sort of local discontinuity between these ends - each side thinks it has to go a long way because of a low gradient, but in-between the two there's a high gradient. Something like that, right?

[KarenRei]

So, I'm going to basically class this as, "Hey OpenFOAM - it's not me, it's you". It's my view that their convergence algorithm is fundamentally flawed - that it'll work most of the time, but not all of the time. I have gone and patched it.

For others experiencing this problem, you should do the following.

1) Open src/thermophysicalModels/specie/lnInclude/thermoI.H and find inline Foam::scalar Foam::species::thermo<Thermo, Type>::T

2) Replace the current maxIter check that aborts with:

if (iter++ > maxIter_ / 2)
{
iter = 0;
do
{
const scalar scaledown = 1.0 - iter / (maxIter_ / 2.0);
Test = Tnew;
Tnew =
(this->*limit)
(Test - scaledown * ((this->*F)(p, Test) - f)/(this->*dFdT)(

if (iter++ > maxIter_ / 2)
{
FatalErrorInFunction
<< "Maximum number of iterations exceeded"
<< abort(FatalError);
}

} while (mag(Tnew - Test) > Ttol);

break;
}


This basically forces convergence by cutting the distance that it may travel per hop, more and more with each cycle, if it fails the first half of the tolerance checks. Though it might be better to just do the whole outer loop like that rather than having the forced convergence be a "in case of failure" scenario - I'm not sure which would be best.

But anyway, to anyone else who encounters this (the main reason I've basically been debugging this out loud).... enjoy!

[arjun]

In case of openfoam, lot of magic with density scaling too happens so take care of that too.

I wish i could help but unfotunately I can not in this case. Good luck.