CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- OpenFOAM (https://www.cfd-online.com/Forums/openfoam/)

- - thermo.correct (https://www.cfd-online.com/Forums/openfoam/99517-thermo-correct.html)

Dear All,

I have a problem with thermo.correct().

I create my own solver (I needed to add a source term - thanks to the explicitSetValue features - to the buoyanPimpleFoam solver).

I edited the hEqn.H file, since I wanted to impose temperature (i.e. enthalpy) and this is my new file:

Code:

{

    fvScalarMatrix hEqn

    (

        fvm::ddt(rho, h)

      + fvm::div(phi, h)

      - fvm::laplacian(turbulence->alphaEff(), h)

     ==

        dpdt

      - (fvc::ddt(rho, K) + fvc::div(phi, K))

      //- (fvc::ddt(rho, K) + fvc::div(phi, 0.5*magSqr(U), "div(phi,K)"))

      + sources(rho, h)

    );



    sources.constrain(hEqn);

    hEqn.solve();



    thermo.correct();



}

When I compile it everything seems to be ok. But when I try to use it, I get the following error:

Code:

lab@lab-laptop:~/Documenti/cases_OF/OF_case09_steady_vs_unsteady/unsteady_new_solver$ buoyantPimpleFoam_Epta 

/*---------------------------------------------------------------------------*\

| =========                 |                                                 |

| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |

|  \\    /   O peration     | Version:  2.1.0                                 |

|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |

|    \\/     M anipulation  |                                                 |

\*---------------------------------------------------------------------------*/

Build  : 2.1.0-0bc225064152

Exec   : buoyantPimpleFoam_Epta

Date   : Apr 05 2012

Time   : 10:18:35

Host   : "lab-laptop"

PID    : 4307

Case   : /home/lab/Documenti/cases_OF/OF_case09_steady_vs_unsteady/unsteady_new_solver

nProcs : 1

sigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).

fileModificationChecking : Monitoring run-time modified files using timeStampMaster

allowSystemOperations : Disallowing user-supplied system call operations



// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

Create time



Create mesh for time = 0





Reading g

Reading thermophysical properties



Selecting thermodynamics package hRhoThermo<pureMixture<constTransport<specieThermo<hConstThermo<perfectGas>>>>>

Reading field U



Reading/calculating face flux field phi



Creating turbulence model



Selecting turbulence model type RASModel

Selecting RAS turbulence model kEpsilon

--> Upgrading k to employ run-time selectable wall functions

    Backup original k to k.old

    Writing updated k

--> Upgrading epsilon to employ run-time selectable wall functions

    Backup original epsilon to epsilon.old

    Writing updated epsilon

--> Creating mut to employ run-time selectable wall functions

    Writing new mut

--> Creating alphat to employ run-time selectable wall functions

    Writing new alphat

kEpsilonCoeffs

{

    Cmu             0.09;

    C1              1.44;

    C2              1.92;

    C3              -0.33;

    sigmak          1;

    sigmaEps        1.3;

    Prt             1;

}



Calculating field g.h



Reading field p_rgh



Creating field dpdt



Creating field kinetic energy K



Creating field source list from sourcesProperties



Selecting source model type scalarExplicitSetValue

    Source: air_infinite

    - selecting cells using cellSet air_infinite

    - selected 2730 cell(s) with volume 1.2



Courant Number mean: 0.000552417 max: 5.02789



PIMPLE: no residual control data found. Calculations will employ 2 corrector loops





Starting time loop



Courant Number mean: 8.24029e-05 max: 0.75

deltaT = 0.0298336

Time = 0.0298336



diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0

PIMPLE: iteration 1

DILUPBiCG:  Solving for Ux, Initial residual = 1, Final residual = 0.0083111, No Iterations 9

DILUPBiCG:  Solving for Uy, Initial residual = 1, Final residual = 0.0074752, No Iterations 6

DILUPBiCG:  Solving for Uz, Initial residual = 1, Final residual = 0.00669065, No Iterations 4

DILUPBiCG:  Solving for h, Initial residual = 0.00100942, Final residual = 6.31122e-06, No Iterations 5

#0  Foam::error::printStack(Foam::Ostream&) in "/opt/openfoam210/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"

#1  Foam::sigFpe::sigHandler(int) in "/opt/openfoam210/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"

#2   in "/lib/libc.so.6"

#3  Foam::hRhoThermo<Foam::pureMixture<Foam::constTransport<Foam::specieThermo<Foam::hConstThermo<Foam::perfectGas> > > > >::calculate() in "/opt/openfoam210/platforms/linux64GccDPOpt/lib/libbasicThermophysicalModels.so"

#4  Foam::hRhoThermo<Foam::pureMixture<Foam::constTransport<Foam::specieThermo<Foam::hConstThermo<Foam::perfectGas> > > > >::correct() in "/opt/openfoam210/platforms/linux64GccDPOpt/lib/libbasicThermophysicalModels.so"

#5  

 in "/home/lab/OpenFOAM/lab-2.1.0/platforms/linux64GccDPOpt/bin/buoyantPimpleFoam_Epta"

#6  __libc_start_main in "/lib/libc.so.6"

#7  

 in "/home/lab/OpenFOAM/lab-2.1.0/platforms/linux64GccDPOpt/bin/buoyantPimpleFoam_Epta"

Floating point exception

The problem is with the thermo.correct() line.

Could anyone help?

Thanks,
Samuele

Hi,

firstly, could you have a look at the enthalpy field right before you call thermo.correct() (and have already solved the energy equation). I guess there might be something wrong in the enthalpy field. Also, try to find the exact line in the thermo.correct() function where the solver crashes.

One thing that puzzles me is your energy equation. Maybe it is the explicitSetValue features I'm not familiar with, but don't have these two terms in your equation different units?

fvm::ddt(rho, h)
+ sources(rho, h)

Best Regards,
Chritistian

I am not very sure about how the explicitSetValue feature works, too.

But I guess that it is correct since if you look at the /opt/openfoam210/applications/solvers/lagrangian/LTSReactingParcelFoam/hsEqn.H you will find this:

Code:

{

    fvScalarMatrix hsEqn

    (

        fvm::ddt(rho, hs)

      + mvConvection->fvmDiv(phi, hs)

      - fvm::laplacian(turbulence->alphaEff(), hs)

     ==

      - fvc::div(phi, 0.5*magSqr(U), "div(phi,K)")

      + parcels.Sh(hs)

      + radiation->Shs(thermo)

      + combustion->Sh()

      + sources(rho, hs)

    );



    sources.constrain(hsEqn);



    hsEqn.solve();



    thermo.correct();



    radiation->correct();



    Info<< "T gas min/max   = " << min(T).value() << ", "

        << max(T).value() << endl;

}

Also, if you write something different you get an error during compilation. So I think it is correct, isn't it?

Samuele

Hi,

does the solver run with the source term? What about the enthalpy field after solving the energy equation?

Regards
Christian

Hi Christian and thanks for answering.

Do you mean `my' solver?

Of course, it doesn't work since I get the error that I posted above.

It stopped at the thermo.correct() line.

BTW, are you familiar with chtMultiRegionFoam? Could I ask you a think?

Hi,

sorry, there was a typo in my last post.

The question is, does the solver run if you remove he line "sources(rho, h)" from the energy equation in your code.

Additionally, could you write out the enthalpy field after solving the energy equation and have a look at it. Are there values at infinity?

At which line in thermo.correct() (inside the thermo library) does the solver crash?

Best regards,
Christian

Quote:

Originally Posted by Chris Lucas (Post 354198)

Yes, if I remove the source it works since it becomes the `original' buoyantPimpleFoam.

I attach the buoyantSimpleFoam_Epta and the buoyantPimpleFoam_Epta.

The former is steady (and it works), while the ladder is unsteady (and it doesn't).

Any idea?