[Dr.Wan_Li]

Hello everyone,

I am modifying the rhoCentralFoam by adding additional terms in mass, momentum and energy balance equation, ( Type 1 model of the following paper: https://iopscience.iop.org/article/1...528/ab4b86/pdf )

The solver crashes with a sigFpe error. I have tried debugging the solver and with my limited understanding I can pinpoint the error to a negative temperature at a cellpoint.
Stack trace :

#0 Foam::error::printStack(Foam::Ostream&) at ~/OpenFOAM/OpenFOAM-6-debug/src/OSspecific/POSIX/printStack.C:218
#1 Foam::sigFpe::sigHandler(int) at ~/OpenFOAM/OpenFOAM-6-debug/src/OSspecific/POSIX/signals/sigFpe.C:106
#2 ? in "/lib/x86_64-linux-gnu/libc.so.6"
#3 ? in "/lib/x86_64-linux-gnu/libm.so.6"
#4 Foam::sutherlandTransport<Foam::species::thermo<Fo am::hConstThermo<Foam::perfectGas<Foam::specie> >, Foam::sensibleInternalEnergy> >::mu(double, double) const at ~/OpenFOAM/OpenFOAM-6-debug/src/thermophysicalModels/specie/lnInclude/sutherlandTransportI.H:124
#5 Foam::hePsiThermo<Foam::psiThermo, Foam::pureMixture<Foam::sutherlandTransport<Foam:: species::thermo<Foam::hConstThermo<Foam::perfectGa s<Foam::specie> >, Foam::sensibleInternalEnergy> > > >::calculate() at ~/OpenFOAM/OpenFOAM-6-debug/src/thermophysicalModels/basic/psiThermo/hePsiThermo.C:55 (discriminator 2)
#6 Foam::hePsiThermo<Foam::psiThermo, Foam::pureMixture<Foam::sutherlandTransport<Foam:: species::thermo<Foam::hConstThermo<Foam::perfectGa s<Foam::specie> >, Foam::sensibleInternalEnergy> > > >::correct() at ~/OpenFOAM/OpenFOAM-6-debug/src/thermophysicalModels/basic/psiThermo/hePsiThermo.C:158
#7 ? in "/home/amt4/OpenFOAM/amt4-6-debug/platforms/linux64GccDPInt64Debug/bin/RNSrFoam_debug"
#8 __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#9 ? in "/home/amt4/OpenFOAM/amt4-6-debug/platforms/linux64GccDPInt64Debug/bin/RNSrFoam_debug"

The test case being solved is a 2D lid-driven cavity problem, with dimensions 10^-7 x 10^-7 x 10^-8 .
Top lid moves at 100m/s (BC fixedValue) and other 3 walls have zerogradient BC for p & T and noSlip for U. since its 2D front and back faces of the cuboid are given empty BC.

While debugging with gdb I see that in line 309 of the solver,

308 e = rhoE/rho - 0.5*magSqr(U);
309 e.correctBoundaryConditions();
310 thermo.correct();

the e.correctBoundaryConditions() function returns a negative value at e[39801], which pops up as error in thermo.correct()

(gdb) display *&e[39801]
2: *&e[39801] = 113562.53726693301
(gdb) n
309 e.correctBoundaryConditions();
2: *&e[39801] = -99394.616383268934

I was wondering if someone could suggest any ways to check my implementaion, especially since its my first time modifying a solver. Let me know if additional files are required to check.
I have doubts on my implementation of energy equation too,

Modified energy equation:

// --- Solve energy
surfaceScalarField sigmaDotU
(
"sigmaDotU",
(
fvc::interpolate(muEff)*mesh.magSf()*fvc::snGrad(U )
+ fvc::dotInterpolate(mesh.Sf(), -tau_v2) //tauMCS
)
& (a_pos*U_pos + a_neg*U_neg)
);
// Below -> N1, N2, N3, N4, N5 and N6 based on RNS paper energy equation : Appendix
volVectorField N1("N1", -(U & fvc::grad(rho))*U - 0.5*(U&U)*fvc::grad(rho));
volVectorField N2("N2",(U&fvc::grad(rho))*(fvc::grad(lnr))+0.5*ma gSqr(fvc::grad(rho))*U/rho);
volVectorField N3("N3",-0.5*magSqr(fvc::grad(rho))*((fvc::grad(lnr))/rho)) ;
volScalarField N41("N41",fvc::div(rho*(U*U) + rho*I/psi - tau_RNS )&fvc::grad(lnr));
volScalarField N42("N42",-U&(fvc::grad(lnr)*fvc::div(rhoU) - fvc::grad(fvc::div(rhoU))));
volScalarField N4("N4",N41+N42);
volScalarField N5("N5",
fvc::laplacian(rho)*(U&fvc::grad(lnr))
- (U&fvc::grad(fvc::laplacian(rho)))
+ 0.5 *magSqr(fvc::grad(rho))*((fvc::div(rhoU))/rho)/rho
);
volScalarField N6
(
"N6",
-0.5*magSqr(fvc::grad(rho))*((fvc::laplacian(rho))/rho)/rho
);


solve
(
fvm::ddt(rhoE)
+ fvc::div(phiEp)
- fvc::div(sigmaDotU)
+ fvc::div(Km*tau_v & fvc::grad(lnr))
);

e = rhoE/rho - 0.5*magSqr(U);
e.correctBoundaryConditions();
thermo.correct();
rhoE.boundaryFieldRef() ==
rho.boundaryField()*
(
e.boundaryField() + 0.5*magSqr(U.boundaryField())

);

if (!inviscid)
{
solve
(
fvm::ddt(rho, e) - fvc::ddt(rho, e)
- fvm::laplacian(turbulence->alphaEff(),e)
- fvc::div(Km*rho*e*fvc::grad(lnr))
- fvc::div((Km*rho*I/psi)&fvc::grad(lnr))
+ fvc::div(Km*N1 + Km*Km*N2 + Km*Km*Km*N3)
+ Km*N4
+ Km*Km*N5
+ Km*Km*Km*N6
);
thermo.correct();
rhoE = rho*(e + 0.5*magSqr(U));
}

Regards

Wan LI