[pamela]

Hello, everybody.
I am trying to solve for compressible, turbulent flow, steady state, porous media, without assuming local thermal equilibrium. Thus, for LTNE I added an energy equation for the solid matrix with an energy source, linked to the fluid energy equation by the convective heat transfer term.

The enthalpies seem not to converge. I've tried making an internal loop (the residuals tell me that actually no more than 2 or 3 iterations are necessary), and limiting the minimum solid enthalpy.
To no avail: the relaxation factors I need just to make the solver run are terrible (0.00001, for both enthalpies). (I tried everything in the book: starting robustly/finishing accurately, the mesh is extremely orthogonal, probes to monitor enthalpies values, very low residuals to stop Simple: 10⁻8, etc.)

Has anybody tried anything similar?
Should I stop trying to make Simple work, and migrate asap to rhoPimpleFoam (adding porosity explicitly via fvOptions)?
Also, I have just found a coupledFvScalarMatrix class that Klas Jareteg uses to couple regions. I'd like to try that to couple the fields he and hs, in order to solve the PDEs at once, instead of segregatedly. But OpenFOAM 6 does not have it, it seems. Has anybody tried it? The implementation seems kind of cumbersome...

Many, many thanks in advance for any clue...

Here's my EEqn.H file:

{
int valori = 5;//number of calculations of both EEqns in one internal iteration.
int counteri = 1;//counter initialization.

volScalarField& he = thermo.he();

fvScalarMatrix EgEqn
(
fvm::div(phi, he)
+ (
he.name() == "e"
? fvc::div(phi, volScalarField("Ekp", 0.5*magSqr(U) + p/rho))
: fvc::div(phi, volScalarField("K", 0.5*magSqr(U)))
)

- fvm::laplacian(turbulence->alphaEff(), he)
- (hConv*areaEsp)*((hs/Cps)-(he/Cp))
==
fvOptions(rho, he)
);

fvScalarMatrix EsEqn
(
fvm::laplacian(alphas, hs)
+ Q
- hConv*areaEsp*((hs/Cps)-(he/Cp))
);


if (hs < minHs)
hs == minHs;//limiting the solid enthalpy at its minimum value.


while (counteri <= valori)
{

EgEqn.relax();
fvOptions.constrain(EgEqn);
EgEqn.solve();
fvOptions.correct(he);

EsEqn.relax();
EsEqn.solve();

if (hs < minHs)
hs == minHs;//limiting the solid enthalpy at its minimum value.

thermo.correct();

counteri = counteri + 1;
}


}