[PierreAG]

Hello foamers,
I am quite new at using OpenFOAM (I have been using it for 1 month).

I am trying to add temperature and heat equation to the interFoam solver, as well as convection.
I started by following this tutorial : http://www.tfd.chalmers.se/~hani/kur...gLIU-final.pdf
In this tutorial, several constants, variables and methods are implemented, like heat capacity cp (one for each phase, cp1 and cp2) or kappa. The tutorial is a bit old, so I had to adapt a few things, butwithout anything else implemented, it compiles and run without errors.

In addition of what is described in the tutorial, I needed to implement a reference temperature Tref and a coefficient beta for each phase (Tref1, beta1, Tref2, beta2), which are dimensionedScalar constants.

I added them in the incompressibleTwoPhaseMixture model the exact same way cp1 and cp2 had been (Or so I think, I checked the code several times to be sure of that ). The model compiles without problems.

I then changed the createFields in the solver, adding besides the heat capacities the declaration and initialization of Tref1-2 and beta1-2, again in the same way that cp1-2 are.
I also made a few changes in density declaration and definition so it can be modified for convection, and added to TEqn.H a line to change this density.

When compiling the solver, no error occurs but this warning messages appears :

In file included from /opt/openfoam6/src/OpenFOAM/lnInclude/postProcess.H:135:0,
                 from myThermoInterFoam.C:53:
./createFields.H: In function ‘int main(int, char**)’:
./createFields.H:72:26: warning: unused variable ‘beta2’ [-Wunused-variable]
 const dimensionedScalar& beta2 = mixture.beta2();
                          ^~~~~
./createFields.H:75:26: warning: unused variable ‘Tref2’ [-Wunused-variable]
 const dimensionedScalar& Tref2 = mixture.Tref2();
                          ^~~~~
./createFields.H:78:26: warning: unused variable ‘beta1’ [-Wunused-variable]
 const dimensionedScalar& beta1 = mixture.beta1();
                          ^~~~~
./createFields.H:81:26: warning: unused variable ‘Tref1’ [-Wunused-variable]
 const dimensionedScalar& Tref1 = mixture.Tref1();
                          ^~~~~

I did not pay much attention to it, because if I use these scalars later in the solver, they are not used in createFIelds.H itself where I declare and define them, so I thought it was only about that.

Anyway, the solver ready to be used (I named it myThermoInterFoam), I made up a case to see the convection into action : It is just two phases, water and air, in a steady state (water filling a half of the domain, the air filling the other half above water), with water hotter than air. In transportProperties, I added, besides cp, Tref and beta for each phase. For water, beta has the value 0, so there is no convection in it (at the end, the purpose is to see convection in the gas phase only).
Here is the transportProperties dict :

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  6
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

phases (water air);

water
{
    transportModel  Newtonian;
    nu              1e-06;
    rho             1000;
    
    //- Adding thermophysical properties
    
    Pr          Pr [0 0 0 0 0 0 0] 3.366;
    cp          cp [0 2 -2 -1 0 0 0] 1433.4;
    Tref        Tref [0 0 0 1 0 0 0] 300;
    beta        beta [0 0 0 -1 0 0 0] 0;
    //- --------------------------------
}

air
{
    transportModel  Newtonian;
    nu              1.48e-05;
    rho             1;
    
    //- Adding thermophysical properties
    
    Pr          Pr [0 0 0 0 0 0 0] 0.87;
    cp          cp [0 2 -2 -1 0 0 0] 1045.2;
    Tref        Tref [0 0 0 1 0 0 0] 200;
    beta        beta [0 0 0 -1 0 0 0] 2e-3;
    //- --------------------------------
}

sigma            0.07;

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

I create the mesh with blockMesh, then use setFields to set the initial conditions.
I checked the mesh with checkMesh, no problem is detected.

I then run the case with myThermoInterFoam, and before even starting the time loop, this error shows up :

Starting time loop

#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::sigSegv::sigHandler(int) at ??:?
#2  ? in "/lib/x86_64-linux-gnu/libc.so.6"
#3  ? in "/lib/x86_64-linux-gnu/libc.so.6"
 Segmentation fault (core dumped)

I tried to search on the forums something looking like that, but all the error messages I could compare to were more precise than this one.
I also tried to add -DFULLDEBUG -g -O0 in Make/options and to change WM_COMPILE_OPTION to 'Debug', but I could not have a more precise error message.

I searched for what was returning me this and found out it was runTime.run(), used as a condition for the while() of the time loop. I am not the administrator, so I could not add print commands in run() to see what was going wrong. But after some research, I printed the dimensionedScalar that I added or changed directly from the createFields.H.

Here is what I get when trying to execute the solver on the case (again the solver compiled without problem, and without the warning messages, which is probably because now Tref1-2 and beta1-2 are used in the print commands). I ran it several times with different printing orders, since the program stops at the first error. :

//Error on beta1



Reading transportProperties

Selecting incompressible transport model Newtonian
Selecting incompressible transport model Newtonian
Reading field T

mixture rho1c = rho [1 -3 0 0 0 0 0] 1000

mixture rho2c = rho [1 -3 0 0 0 0 0] 1

mixture cp1 = cp [0 2 -2 -1 0 0 0] 1433.4

mixture cp2 = cp [0 2 -2 -1 0 0 0] 1045.2

mixture beta2 =  [2.122e-314 4.94066e-324 4.94066e-324 -1 0 5.73116e-322 0] 0.002

mixture Tref2 =  [6.95254e-310 6.95254e-310 4.64802e-310 6.95254e-310 2.47033e-323 2.17408e-312 0] 6.95254e-310

mixture beta1 = #0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::sigSegv::sigHandler(int) at ??:?
#2  ? in "/lib/x86_64-linux-gnu/libc.so.6"
#3  ? in "/lib/x86_64-linux-gnu/libc.so.6"
#4  _IO_file_xsputn in "/lib/x86_64-linux-gnu/libc.so.6"
#5  fwrite in "/lib/x86_64-linux-gnu/libc.so.6"
#6  std::basic_ostream<char, std::char_traits<char> >& std::__ostream_insert<char, std::char_traits<char> >(std::basic_ostream<char, std::char_traits<char> >&, char const*, long) in "/usr/lib/x86_64-linux-gnu/libstdc++.so.6"
#7  Foam::OSstream::write(Foam::word const&) at ??:?
#8  Foam::operator<<(Foam::Ostream&, Foam::word const&) at ??:?
#9  Foam::Ostream& Foam::operator<< <double>(Foam::Ostream&, Foam::dimensioned<double> const&) in "/home/etudiant/OpenFOAM/etudiant-6/platforms/linux64GccDPInt32Opt/bin/myThermoInterFoam"
#10  ? in "/home/etudiant/OpenFOAM/etudiant-6/platforms/linux64GccDPInt32Opt/bin/myThermoInterFoam"
#11  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#12  ? in "/home/etudiant/OpenFOAM/etudiant-6/platforms/linux64GccDPInt32Opt/bin/myThermoInterFoam"
Erreur de segmentation (core dumped)

//error on Tref1

Reading transportProperties

Selecting incompressible transport model Newtonian
Selecting incompressible transport model Newtonian
Reading field T

mixture rho1c = rho [1 -3 0 0 0 0 0] 1000

mixture rho2c = rho [1 -3 0 0 0 0 0] 1

mixture cp1 = cp [0 2 -2 -1 0 0 0] 1433.4

mixture cp2 = cp [0 2 -2 -1 0 0 0] 1045.2

mixture beta2 =  [2.122e-314 4.94066e-324 4.94066e-324 -1 0 5.73116e-322 0] 0.002

mixture Tref2 =  [6.95335e-310 6.95335e-310 4.64164e-310 6.95335e-310 2.47033e-323 2.17408e-312 0] 6.95335e-310

mixture Tref1 = �?���qU����deltaN�?��rL��Y�h>�����`N���`���qU�����nHatf����IOobjecteta�;Ї��������
0����betap��qU


��tMb`?XM�P_�qUp��qU@�g ;��qU�4��qU`���qU����qU(���IOobjectH����K�h����
0�X��qU�����@��p��qU


|'qU@���qUp��qU��
�P+��qU�4��qU��������T�-�����Z�Z}��Y}�cL�������^L�__-��cL��>-�x����x��}0ř}

--> FOAM FATAL IO ERROR: 
error in IOstream "Sout" for operation Ostream& operator<<(Ostream&, const word&)

file: Sout at line 51.

    From function virtual bool Foam::IOstream::check(const char*) const
    in file db/IOstreams/IOstreams/IOstream.C at line 96.

FOAM exiting

Here are the things I noticed :

• rho1 and rho2, which are already in the interFoam default solver, are read correctly.
• cp1 and cp2, which I implemented following the tutorial, are read correctly.
• Tref1, Tref2, beta1 and beta2, which also are dimensionedScalar, and which I'm pretty sure (again, I have checked the code several times) are implemented the same way cp1 and cp2 are, are completely messed up.
• More precisely, Tref2 and beta2, which are the air phase constants, and are supposed to be in Kelvin and Kelvin^-1, have crazy dimensions, being floats near 0 in all measure units coefficients, including the Kelvin coefficient.
• Tref2 value is also incorrect, but beta2 value is correct.
• beta1 sends back a segmentation fault.
• Tref1 sends back a FATAL IO ERROR.

Again, cp, Tref and beta are supposed to be implemented in the same way.
I tried to recompile the model, the solver, and execute the case with the debug options mentioned above, but I could not have another error messages.
I also searched for these errors on forums, but again I could not find anything to help me.

Does anyone knows anything about my problem ?

I do not have a dropbox account, so I put the files in an archive on Google Drive :

• "Case" is the case
• "myThermoInterFoam" is the solver
• "incompressibleTwoPhaseMixture" is the modified model. The original one can be found in transportModels/incompressible/

https://drive.google.com/file/d/1C-L...ew?usp=sharing
Maybe you will see something I did not.

Thanks in advance for the help,

Pierre

[alexeym]

Hi,

I am not quite sure, OpenFOAM(R)'s immiscibleIncompressibleTwoPhaseMixture uses your incompressibleTwoPhaseMixture class. It would be safer to create your own immiscibleIncompressibleThermalTwoPhaseMixture class with additional properties and use it in the solver.

[ELwardi]

Hi,

Apparently, this is a referenced variable gone out of scope, check how your nuModel class returns beta and Tref.

Also, GDB may be useful in debugging segfaults .

Regards,

[PierreAG]

Hi,

Thank you for your answers.

Apparently, this is a referenced variable gone out of scope, check how your nuModel class returns beta and Tref.

I printed the variables directly from the incompressibleTwoPhaseMixture model, in the constructor :

Foam::incompressibleTwoPhaseMixture::incompressibleTwoPhaseMixture
(
    const volVectorField& U,
    const surfaceScalarField& phi
)
:
    IOdictionary
    (
        IOobject
        (
            "transportProperties",
            U.time().constant(),
            U.db(),
            IOobject::MUST_READ_IF_MODIFIED,
            IOobject::NO_WRITE
        )
    ),
    twoPhaseMixture(U.mesh(), *this),

    nuModel1_
    (
        viscosityModel::New
        (
            "nu1",
            subDict(phase1Name_),
            U,
            phi
        )
    ),
    nuModel2_
    (
        viscosityModel::New
        (
            "nu2",
            subDict(phase2Name_),
            U,
            phi
        )
    ),

    rho1_("rho", dimDensity, nuModel1_->viscosityProperties()),
    rho2_("rho", dimDensity, nuModel2_->viscosityProperties()),
    
    /*//- Adding thermophysical properties ---------------- (old way)
    cp1_(nuModel1_->viscosityProperties().lookup("cp")),
    cp2_(nuModel2_->viscosityProperties().lookup("cp")),
    Pr1_(nuModel1_->viscosityProperties().lookup("Pr")),
    Pr2_(nuModel2_->viscosityProperties().lookup("Pr")),
    Tref1_(nuModel1_->viscosityProperties().lookup("Tref")),
    Tref2_(nuModel2_->viscosityProperties().lookup("Tref")),
    beta1_(nuModel1_->viscosityProperties().lookup("beta")),
    beta2_(nuModel2_->viscosityProperties().lookup("beta")),
    //- -------------------------------------------------*/
    
    //- Adding thermophysical properties ---------------- (new way)
    cp1_("cp", dimSpecificHeatCapacity, nuModel1_->viscosityProperties()),
    cp2_("cp", dimSpecificHeatCapacity, nuModel2_->viscosityProperties()),
    Pr1_("Pr", dimless, nuModel1_->viscosityProperties()),
    Pr2_("Pr", dimless, nuModel2_->viscosityProperties()),
    Tref1_("Tref", dimTemperature, nuModel1_->viscosityProperties()),
    Tref2_("Tref", dimTemperature, nuModel2_->viscosityProperties()),
    beta1_("beta", dimensionSet(0, 0, 0, -1, 0, 0, 0), nuModel1_->viscosityProperties()),
    beta2_("beta", dimensionSet(0, 0, 0, -1, 0, 0, 0), nuModel2_->viscosityProperties()),
    //- -------------------------------------------------
    
    U_(U),
    phi_(phi),

    nu_
    (
        IOobject
        (
            "nu",
            U_.time().timeName(),
            U_.db()
        ),
        U_.mesh(),
        dimensionedScalar("nu", dimViscosity, 0),
        calculatedFvPatchScalarField::typeName
    )
{
    calcNu();
    
    Info<< "Reading from model :\n" << endl ;
    
    Info<< "rho1 = " << rho1_ << endl;
    Info<< "rho2 = " << rho2_ << endl;
    Info<< "cp1 = " << cp1_ << endl;
    Info<< "cp2 = " << cp2_ << endl;
    Info<< "Tref1 = " << Tref1_ << endl;
    Info<< "Tref2 = " << Tref2_ << endl;
    Info<< "beta1 = " << beta1_ << endl;
    Info<< "beta2 = " << beta2_ << "\n" << endl;
    
    Info<< "End model reading\n" << endl;
}

(I also changed the way they are read, but it didn't change anything).

It appears that here the variables are at first read correctly :

Reading from model :

rho1 = rho [1 -3 0 0 0 0 0] 1000
rho2 = rho [1 -3 0 0 0 0 0] 1
cp1 = cp [0 2 -2 -1 0 0 0] 1433.4
cp2 = cp [0 2 -2 -1 0 0 0] 1045.2
Tref1 = Tref [0 0 0 1 0 0 0] 300
Tref2 = Tref [0 0 0 1 0 0 0] 200
beta1 = beta [0 0 0 -1 0 0 0] 0
beta2 = beta [0 0 0 -1 0 0 0] 0.002

End model reading

Apparently, the problem comes from the "get" functions :

//- Return const-access to phase1 density
        const dimensionedScalar& rho1() const
        {
            Info<< "const-access to rho1 : " << rho1_ << endl;
            return rho1_;
        }

        //- Return const-access to phase2 density
        const dimensionedScalar& rho2() const
        {
            Info<< "const-access to rho2 : " << rho2_ << endl;
            return rho2_;
        };

        //- Adding thermodynamical properties
            const dimensionedScalar& Pr1() const
            {
                Info<< "const-access to Pr1 : " << Pr1_ << endl;
                return Pr1_;
            }
            
            const dimensionedScalar& Pr2() const
            {
                Info<< "const-access to Pr2 : " << Pr2_ << endl;
                return Pr2_;
            };
            
            const dimensionedScalar& cp1() const
            {
                Info<< "const-access to cp1 : " << cp1_ << endl;
                return cp1_;
            }
            
            const dimensionedScalar& cp2() const
            {
                Info<< "const-access to cp2 : " << cp2_ << endl;
                return cp2_;
            };
            
            const dimensionedScalar& Tref1() const
            {
                Info<< "const-access to Tref1 : " << Tref1_ << endl;
                return Tref1_;
            }
        
            const dimensionedScalar& Tref2() const
            {
                Info<< "const-access to Tref2 : " << Tref2_ << endl;
                return Tref2_;
            };
            
            const dimensionedScalar& beta1() const
            {
                Info<< "const-access to beta1 : " << beta1_ << endl;
                return beta1_;
            }
            
            const dimensionedScalar& beta2() const
            {
                Info<< "const-access to beta2 : " << beta2_ << endl;
                return beta2_;
            };

Here is what it returns :

const-access to rho1 : rho [1 -3 0 0 0 0 0] 1000
mixture rho1 = rho [1 -3 0 0 0 0 0] 1000

const-access to rho2 : rho [1 -3 0 0 0 0 0] 1
mixture rho2 = rho [1 -3 0 0 0 0 0] 1

const-access to cp1 : cp [0 2 -2 -1 0 0 0] 1433.4
mixture cp1 = cp [0 2 -2 -1 0 0 0] 1433.4

const-access to cp2 : cp [0 2 -2 -1 0 0 0] 1045.2
mixture cp2 = cp [0 2 -2 -1 0 0 0] 1045.2

const-access to beta2 :  [2.122e-314 4.94066e-324 4.94066e-324 -1 0 5.73116e-322 0] 0.002
mixture beta2 =  [2.122e-314 4.94066e-324 4.94066e-324 -1 0 5.73116e-322 0] 0.002

const-access to Tref2 :  [6.95323e-310 6.95323e-310 4.66412e-310 6.95323e-310 2.47033e-323 2.17408e-312 0] 6.95323e-310
mixture Tref2 =  [6.95323e-310 6.95323e-310 4.66412e-310 6.95323e-310 2.47033e-323 2.17408e-312 0] 6.95323e-310

const-access to Tref1 : �?B���Ux�c�deltaN�?��rL��Y�h�c���c�p����U��c�nHatf�cIOobjecteta�;0�c�P�c�
0p�c�beta�����U


��tMb`?XM�#:��U�����U@�g0���U�����Up����U0�%��U��cIOobjectc��c�K��c�
0g@)��U��c�@�������U


|'�U��'��U�����U�����&��U�����UZ�� ��� �S�&�N�&_O�&�S�&�.�&��c�x�]!0�=!�=!��=! �=!

--> FOAM FATAL IO ERROR: 
error in IOstream "Sout" for operation Ostream& operator<<(Ostream&, const word&)

file: Sout at line 74.

    From function virtual bool Foam::IOstream::check(const char*) const
    in file db/IOstreams/IOstreams/IOstream.C at line 96.

 FOAM exiting

The "mixture [variable] = ..." come from the createFields.H.

The errors that Tref1 and beta1 return are the same as before.
It might come from the fact that these "get" functions have a reference type, however I am not skilled enough in C++ nor in OpenFOAM to pinpoint the issue...

Also, GDB may be useful in debugging segfaults .

I am currently trying to use GDB to find something useful, but since OpenFOAM uses a large amount of customed classes, objects, and functions, it makes this work a bit difficult (and I am not even sure of what I am looking for).

If I can't find anything, I will try what Alexei proposed :

OpenFOAM(R)'s immiscibleIncompressibleTwoPhaseMixture uses your incompressibleTwoPhaseMixture class. It would be safer to create your own immiscibleIncompressibleThermalTwoPhaseMixture class with additional properties and use it in the solver.

I did not do it in the first place because the tutorial was referring to the model where rho1 and rho2 were implemented, but since nuModel1 and 2 are protected and not private, and that immiscibleIncompressibleTwoPhaseMixture inherit from incompressibleTwoPhaseMixture, it may be possible.

Feel free to tell if you have any advice.

Thank you in advance,

Pierre

[PierreAG]

Hi again,

I eventually implemented the variables in immiscibleIncompressibleTwoPhaseMixture and it seems to work without issue !

I am not sure to know why though...

Anyway, thank you Alexei !

Regards,

Pierre