CFD Online Logo CFD Online URL
Home > Forums > Software User Forums > OpenFOAM

Negative temperature/internal energy in modified rhoCentralFoam

Register Blogs Members List Search Today's Posts Mark Forums Read

LinkBack Thread Tools Search this Thread Display Modes
Old   October 1, 2020, 09:43
Exclamation Negative temperature/internal energy in modified rhoCentralFoam
New Member
Wan Li
Join Date: Jun 2019
Location: Edinburgh
Posts: 3
Rep Power: 7
Dr.Wan_Li is on a distinguished road
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: )

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/"
#3 ? in "/lib/x86_64-linux-gnu/"
#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/"
#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
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",
- (U&fvc::grad(fvc::laplacian(rho)))
+ 0.5 *magSqr(fvc::grad(rho))*((fvc::div(rhoU))/rho)/rho
volScalarField N6

+ fvc::div(phiEp)
- fvc::div(sigmaDotU)
+ fvc::div(Km*tau_v & fvc::grad(lnr))

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


if (!inviscid)
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
rhoE = rho*(e + 0.5*magSqr(U));

Wan LI
Dr.Wan_Li is offline   Reply With Quote


energy equation, negative temperature, openfoam 6, rhocentralfoam, solver modification

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On

Similar Threads
Thread Thread Starter Forum Replies Last Post
Energy equation - where is the turbulent kinetic energy? usv001 OpenFOAM Programming & Development 1 January 25, 2022 15:04
Energy equation in rhoCentralFoam (revisited) usv001 OpenFOAM Programming & Development 3 July 5, 2018 05:52
[blockMesh] Errors during blockMesh meshing Madeleine P. Vincent OpenFOAM Meshing & Mesh Conversion 51 May 30, 2016 10:51
Energy dissipation and Drag coefficient Freeman Main CFD Forum 10 January 27, 2006 07:42
Why FVM for high-Re flows? Zhong Lei Main CFD Forum 23 May 14, 1999 13:22

All times are GMT -4. The time now is 00:58.