[samiam1000]

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:

{
    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:

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

[Chris Lucas]

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

[samiam1000]

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:

{
    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

[Chris Lucas]

Hi,

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

Regards
Christian

[samiam1000]

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?

[Chris Lucas]

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

[samiam1000]

Originally Posted by Chris Lucas

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

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?