[wang219910611]

Hi everyone, I tried to do some adaptations to the cavitatingFoam. And the compile went well. However, when I tried to apply it to the throttle/ras/cavitatingFoam case, following error appears. I compared it with my succeed log file with cavitatingFoam, and found out the problem might be with the compressibility model. But I can't really find which part of the code leads to this. Can anyone offer any clue?


Create time

Create mesh for time = 0

Reading thermodynamicProperties

Reading field p

Reading field U

Reading/calculating face flux field phiv

Reading/calculating face flux field phi

Reading transportProperties

Selecting incompressible transport model Newtonian
Selecting incompressible transport model Newtonian
#0 Foam::error:rintStack(Foam::Ostream&) at ??:?
#1 Foam::sigSegv::sigHandler(int) at ??:?
#2 in "/lib/x86_64-linux-gnu/libc.so.6"
#3 std::basic_string<char, std::char_traits<char>, std::allocator<char> >::basic_string(std::string const&) in "/usr/lib/x86_64-linux-gnu/libstdc++.so.6"
#4 Foam::regIOobject::regIOobject(Foam::regIOobject const&) at ??:?
#5 Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>::GeometricField(Foam::GeometricFiel d<double, Foam::fvPatchField, Foam::volMesh> const&) at ??:?
#6
at ??:?
#7 __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#8
at ??:?
Segmentation fault (core dumped)

[alexeym]

Hi,

Maybe I missed something but are you proposing people to guess your code modifications?

[wang219910611]

Originally Posted by alexeym

Hi,

Maybe I missed something but are you proposing people to guess your code modifications?

Sorry for this...I have attach the modified code as attachments.

It seems to me that the error happened in the creatField.H. But I only added few lines there as following without other modifications.

volScalarField rhov(twoPhaseProperties.rho1());
volScalarField rhol(twoPhaseProperties.rho2());
volScalarField B=rhol/psil;
double N=11;
double r=1.4;

Regards

[alexeym]

Well,

Don't know what OF version/compiler you're using but this

volScalarField rhov(twoPhaseProperties.rho1());
volScalarField rhol(twoPhaseProperties.rho2());

should not compile (at least OF 2.3.0 and gcc 4.8.2 raise compilation error), as there's no volScalarField constructor with just dimensionedScalar as a parameter.

These lines should be something like:

    rhol
    (
        IOobject
        (
            "rhol",
            mesh.time().timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh,
        twoPhaseProperties.rho1()
    )

if you really need rhol as volScalarField.

[wang219910611]

Originally Posted by alexeym

Well,

Don't know what OF version/compiler you're using but this

volScalarField rhov(twoPhaseProperties.rho1());
volScalarField rhol(twoPhaseProperties.rho2());

should not compile (at least OF 2.3.0 and gcc 4.8.2 raise compilation error), as there's no volScalarField constructor with just dimensionedScalar as a parameter.

These lines should be something like:

    rhol
    (
        IOobject
        (
            "rhol",
            mesh.time().timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh,
        twoPhaseProperties.rho1()
    )

if you really need rhol as volScalarField.

Hi, thanks for the advice. Actually I forced the compile to pass by modifying the declaration of rho1 and rho2 in TwoPhaseMixtuer.H as volScalarField. Is this alright?
And I need the rhol and rhov as volScalarField beacuse there is a function in 0Equn.H:

{
volScalarField pgl=(rhov*rhol*(rhov-rhol))/(magSqr(rhov)*psil-magSqr(rhol)*psiv);
if ( rho>rhol)

p=pSat+(B/N)*(pow((rho/rhol),N)-1.0);

else if (rho<rhov)

p=pSat*pow((rho/rhov),r);

else

p=pSat+pgl*log10(rhov*psil*(rhol+alphav*(rhov-rhol))/(rhol*(rhov*psil-alphav*(rhov*psil-rhol*psiv))));

}
Or is there other method I can do this?

[alexeym]

Hi,

I can't answer your question as I don't see your redefinition. Maybe it's OK, maybe it is the reason for the runtime error.

As I don't know what else you've redefined (and think that one should guess it), I will assume variables in you expression have types as in original cavitatingFoam. This

{
volScalarField pgl=(rhov*rhol*(rhov-rhol))/(magSqr(rhov)*psil-magSqr(rhol)*psiv);
if ( rho>rhol)

p=pSat+(B/N)*(pow((rho/rhol),N)-1.0);

else if (rho<rhov)

p=pSat*pow((rho/rhov),r);

else

p=pSat+pgl*log10(rhov*psil*(rhol+alphav*(rhov-rhol))/(rhol*(rhov*psil-alphav*(rhov*psil-rhol*psiv))));

}

can be rewritten as

1. Iterate over volume field and check your condition is every cell.

dimensionedScalar pgl=(rhov*rhol*(rhov-rhol))/(magSqr(rhov)*psil-magSqr(rhol)*psiv);

forAll(rho, cellI)
{
    if ( rho[cellI] > rhol.value())
        p[cellI] = pSat + (B/N)*(pow((rho/rhol),N)-1.0);
    else if (rho[cellI]<rhov.value())
        p[cellI] = pSat*pow((rho/rhov),r);
    else
        p[cellI] = pSat + pgl*log10(rhov*psil*(rhol+alphav[cellI]*(rhov-rhol))/(rhol*(rhov*psil-alphav[cellI]*(rhov*psil-rhol*psiv))));
}

p.correctBoundaryConditions();

Maybe you'll need to add more .value() things (like in rho[cellI]<rhov.value()) as rho[cellI] is scalar while rhov is dimensionedScalar.

2. These is sign function which return -1 if expression is below zero, 0 is expression is 0, and 1 if expression is above zero. You need 3 conditions: rho - rhol, rho - rhov, everything else. Expr1 is value of p for rhol < rho, expr2 is value of p for rho < rhov, expr3 is everything else. Also assuming rhov < rhol:

const volScalarField s1(sign(rho - rhol));
const volScalarField s2(sign(rho - rhov));

// if rho < rhov: s1 == -1 and s2 == -1 -> expr2
// if rhov < rho < rhol: s1 == 1 and s2 == -1 -> expr3
// if rhol < rho: s1 == 1 and s2 == 1 -> exp1

p = 0.25*(1 - s1)*(1 - s2)*expr2 + 0.25*(1 + s1)*(1 - s2)*exp3 + 0.5*(1 + s2)*exp1

Not quite sure the expression is completely correct but I think you've got the idea.

[wang219910611]

Hi, thanks for the brilliant suggestion, I have been stuck on the expression for a long time. I have learnt a lesson from here. Now I am able to wmake the solver without error.

However, following error appears again when I am trying to apply it to the throttle case.

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create mesh for time = 0

Reading thermodynamicProperties

Reading field p

Reading field U

Reading/calculating face flux field phiv

Reading/calculating face flux field phi

Reading transportProperties

Selecting incompressible transport model Newtonian
Selecting incompressible transport model Newtonian
#0 Foam::error:rintStack(Foam::Ostream&) at ??:?
#1 Foam::sigSegv::sigHandler(int) at ??:?
#2 in "/lib/x86_64-linux-gnu/libc.so.6"
#3 std::basic_string<char, std::char_traits<char>, std::allocator<char> >::basic_string(std::string const&) in "/usr/lib/x86_64-linux-gnu/libstdc++.so.6"
#4 Foam::regIOobject::regIOobject(Foam::regIOobject const&) at ??:?
#5 Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>::GeometricField(Foam::GeometricFiel d<double, Foam::fvPatchField, Foam::volMesh> const&) at ??:?
#6
at ??:?
#7 __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#8
at ??:?
Segmentation fault (core dumped)
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

And this is how my relevant creatField.H looks like now.

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "createPhiv.H"
#include "compressibleCreatePhi.H"

Info<< "Reading transportProperties\n" << endl;

incompressibleTwoPhaseMixture twoPhaseProperties(U, phiv);

volScalarField& alphav(twoPhaseProperties.alpha1());
alphav.oldTime();

volScalarField& alphal(twoPhaseProperties.alpha2());

dimensionedScalar rhov(twoPhaseProperties.rho1());
dimensionedScalar rhol(twoPhaseProperties.rho2());


volScalarField B=rhol/psil;

double N=11;
double r=1.4;
Info<< "Creating compressibilityModel\n" << endl;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

[wang219910611]

Thanks for being so patient.

For the convenient, I have attached all the modified code in the attachment.
For the pEqn.H, it is a immature version now. I planned to add more to it, such as nonorthogonal correction when I am able to run the case.

[alexeym]

Hm, you've changed rhol and rhov into dimensioned scalars and again trying to initialize volScalarField B from ratio of these scalars using = operator. It is OK but you have to create B properly first (via copy constructor from other volume field or via constructor from components). Will B be location dependent or it's another constant like rhol and rhov?

[xianqiejiao]

Originally Posted by alexeym

Hm, you've changed rhol and rhov into dimensioned scalars and again trying to initialize volScalarField B from ratio of these scalars using = operator. It is OK but you have to create B properly first (via copy constructor from other volume field or via constructor from components). Will B be location dependent or it's another constant like rhol and rhov?

Ah, sorry for this. B is not location dependent. So I stated B as dimensionedScalar again as following

Info<< "Reading transportProperties\n" << endl;

    incompressibleTwoPhaseMixture twoPhaseProperties(U, phiv);

    volScalarField& alphav(twoPhaseProperties.alpha1());
    alphav.oldTime();

    volScalarField& alphal(twoPhaseProperties.alpha2());

    dimensionedScalar rhov(twoPhaseProperties.rho1());
    dimensionedScalar rhol(twoPhaseProperties.rho2());

    dimensionedScalar  B=rhol/psil;   
    double  N=11;
    double  r=1.4;
    Info<< "Creating compressibilityModel\n" << endl;

However, after I did this, the terminal told me that the dimensions in 0pEqn didn't match again....

 forAll(rho,cellI)

     {
      if (rho[cellI]>rhol.value())

         p[cellI]=pSat.value()+(B/N)*(pow((rho/rhol),N)-1.0);

      else if (rho[cellI]<rhov.value())
     
         p[cellI]=pSat.value()*pow((rho/rhov),r);
    
      else 
  
         p[cellI]=pSat.value()+pgl*log10(rhov*psil*(rhol+alphav*(rhov-rhol))/(rhol*(rhov*psil-alphav*(rhov*psil-rhol*psiv))));
     }
     p.correctBoundaryConditions();

I noticed that alphav was volScalarField and I can't use .value() to it. So I tried this, which didn't work out.

p[cellI]=pSat+pgl*log10(rhov*psil*(rhol+alphav[cellI]*(rhov-rhol))/(rhol*(rhov*psil-alphav[cellI]*(rhov*psil-rhol*psiv))));

. I also tried to add .value() to every dimensionedScalar, but then something goes wrong with log10 function.......

[xianqiejiao]

Thanks again for being so patient. Or could you please tell me is there any place I can learn how to deal with these?

[alexeym]

Well,

Basically you've got 3 different types here:

1. scalar
2. dimensionedScalar
3. volScalarField

let's assume variables: S is scalar, DS is dimensioned scalar, VSF is volScalarField.

1. DS.value() is scalar.
2. VSF[cellI] (or VSF.internalField()[cellI]) is scalar.

There are certain assignments that are valid, for example, you can assign VSF = DS. For the assignment to be valid, VSF should be correctly constructed. Either as a copy of another volScalarField, or from components.

Returning to your code:

p[cellI]=pSat.value()+(B/N)*(pow((rho/rhol),N)-1.0)

p[cellII] is scalar, pSat.value() is scalar, B/N is dimensionedScalar, rho is volScalarField, so in this line you're trying to assign volScalarField (which is basically as list) to scalar. So the error:

0pEqn.H:23:18: erreur: cannot convert ‘Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> >’ to ‘double’ in assignment

Other errors have the same nature, except this one:

In file included from mycavitatingFoam.C:51:0:
createFields.H:56:53: erreur: no matching function for call to ‘Foam::dimensioned<double>::dimensioned(const volScalarField&)’
     dimensionedScalar rhol(twoPhaseProperties.rho2());

this one is due to your redefinition (funny one btw).

If you'd like to learn about iterating over volScalarField (or volVectorField):

$ cd $FOAM_APP/solvers
$ grep -r forAll *

Your case is similar to multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/turbulentBreakUp/turbulentBreakUp.C or combustion/PDRFoam/PDRModels/XiGModels/basicXiSubG/basicXiSubG.C (just two arbitrary examples).

[xianqiejiao]

Hi, thanks for the advise. I have looked through the example files, redefined my rhov and rhol as scalar and managed to write my 0pEqn like this, no more no match error this time :

dimensionedScalar rhov00(twoPhaseProperties.rho1());
    dimensionedScalar rhol00(twoPhaseProperties.rho2());
    scalar rhov=rhov00.value();
    scalar rhol=rhol00.value();
    scalar B=rhol/psil.value();   
    scalar N=11;
    scalar r=1.4;

forAll(rho,cellI)

     {
      if (rho[cellI]>rhol)

         p[cellI]=pSat.vlaue()+(B/N)*(pow((rho[cellI]/rhol),N)-scalar(1));

      else if (rho[cellI]<rhov)
     
         p[cellI]=pSat.value()*pow((rho[cellI]/rhov),r);
    
      else 
  
         p[cellI]=pSat.value()+pgl.value()*log10(rhov*psil.value()*(rhol+alphav[cellI]*(rhov-rhol))/(rhol*(rhov*psil.value()-alphav[cellI]*(rhov*psil.value()-rhol*psiv.value()))));
     }

However, new errors showed up....."pow" and "log" issue this time:

0pEqn.H:27:55: error: call of overloaded ‘pow(double, Foam::scalar&)’ is ambiguous
          p[cellI]=pSat.value()*pow((rho[cellI]/rhov),r);

0pEqn.H:31:179: error: call of overloaded ‘log10(Foam::scalar)’ is ambiguous
          p[cellI]=pSat.value()+pgl.value()*log10(rhov*psil.value()*(rhol+alphav[cellI]*(rhov-rhol))/(rhol*(rhov*psil.value()-alphav[cellI]*(rhov*psil.value()-rhol*psiv.value()))));

I tried to look for the examples for them from the solver, but I just found one and I didn't see much help: solvers/multiphase/twoPhaseEulerFoam/interfacialModels/aspectRatioModels/VakhrushevEfremov, which is a simple one as following:

{
    volScalarField Ta(pair_.Ta());

    return
        neg(Ta - scalar(1))*scalar(1)
      + pos(Ta - scalar(1))*neg(Ta - scalar(39.8))
       *pow3(0.81 + 0.206*tanh(1.6 - 2*log10(max(Ta, scalar(1)))))
      + pos(Ta - scalar(39.8))*0.24;
}

[alexeym]

There are lots of fancy errors like

In file included from pEqn.H:37:0,
                 from mycavitatingFoam.C:82:
0pEqn.H:23:24: erreur: ‘Foam::dimensionedScalar’ has no member named ‘vlaue’
          p[cellI]=pSat.vlaue()+(B/N)*(std::pow((rho[cellI]/rhol),N)-1.0);

Concerning errors with log10 and pow functions, now you're using simple scalars (aka doubles), so you can use functions from standard library. So just prepend function names with std namespace, i.e.

pow -> std::pow
log10 -> std::log10

Concerning the example you've shown (from solvers/multiphase/twoPhaseEulerFoam/interfacialModels/aspectRatioModels/VakhrushevEfremov), it's more idiomatic version of expression with sign function I've posted in the message #6.

[xianqiejiao]

Hi, I have been able to solve the fancy problem by predefine a scalar pSat

scalar pSat00=pSat.value();

and use pSat00 instead in the following equations.

Now I am able to run the case, although it will explode at last. I think I am going to fix the equations now.

Thanks for all the advise and patience all the way.

Best.