CFD Online Discussion Forums - Adding transport equation to twoPhaseEulerFoam (OF231)

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

- OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)

- - Adding transport equation to twoPhaseEulerFoam (OF231) (https://www.cfd-online.com/Forums/openfoam-programming-development/148252-adding-transport-equation-twophaseeulerfoam-of231.html)

Adding transport equation to twoPhaseEulerFoam (OF231)

Hello all,
I'm struggling with the implementation of a scalar transport equation in twoPhaseEulerFoam. I've done that in various solvers before but this one seems to be "special".
Also, the explanation in the wiki (http://openfoamwiki.net/index.php/Ho...sport_equation) isn't sufficient for this solver.

At the moment, the psiEqn.H looks like this:

Code:

{

    fvScalarMatrix psiEqn

    (

        fvm::ddt(alpha1, rho1, psi) + fvm::div(alphaRhoPhi1, psi)

      - fvm::Sp(contErr1, psi)



      - fvm::laplacian

        (

            fvc::interpolate(alpha1)

           *fvc::interpolate(thermo1.alphaEff(phase1.turbulence().mut())),

            psi

        )

      ==

    psiSource*rho1*alpha1

    );



    psiEqn.relax();

    psiEqn.solve();



    Info<< "min " <<psi.name()

        << " " << min(psi).value() << endl;

}

In createFields.H, I added the lines

Code:

    volScalarField psi

    (

        IOobject

        (

            "psi",

            runTime.timeName(),

            mesh,

            IOobject::NO_READ,

            IOobject::AUTO_WRITE

        ),

        mesh,

        dimensionedScalar("zero", dimless, 0.0),

    "zeroGradient"

    );

    volScalarField psiSource

    (

        IOobject

        (

            "psiSource",

            runTime.timeName(),

            mesh,

            IOobject::NO_READ,

            IOobject::NO_WRITE

        ),

        mesh,

        dimensionedScalar("zero", dimless/dimTime, 0.01),

    "zeroGradient"

    );

and in twoPhaseEulerFoam.C Ii included the psiEqn.H right after the energy equations.

Also the divScheme for "div\(alphaRhoPhi.*,psi\)" was set to "Gauss limitedLinear 1;" (upwind doesn't help) and an entry in fvSolution was done according to the one of the energy equation (smoothSolver etc.). I also tried different solvers but that didn't make a change at all.

When I run the solver on the injection tutorial case, it crashes in the moment, when it should solve psiEqn:

Code:

Starting time loop



Courant Number mean: 0 max: 0

Max Ur Courant Number = 0

Time = 0.005



PIMPLE: iteration 1

MULES: Solving for alpha.air

MULES: Solving for alpha.air

alpha.air volume fraction = 0.293333  Min(alpha1) = 0  Max(alpha1) = 1

smoothSolver:  Solving for e.air, Initial residual = 1, Final residual = 7.75344e-14, No Iterations 1

smoothSolver:  Solving for e.water, Initial residual = 0.00402093, Final residual = 0.00367231, No Iterations 1000

min T.air 300

min T.water 300

#0  Foam::error::printStack(Foam::Ostream&) at ??:?

#1  Foam::sigFpe::sigHandler(int) at ??:?

#2   in "/lib/x86_64-linux-gnu/libc.so.6"

#3  Foam::symGaussSeidelSmoother::smooth(Foam::word const&, Foam::Field<double>&, Foam::lduMatrix const&, Foam::Field<double> const&, Foam::FieldField<Foam::Field, double> const&, Foam::UPtrList<Foam::lduInterfaceField const> const&, unsigned char, int) at ??:?

#4  Foam::symGaussSeidelSmoother::smooth(Foam::Field<double>&, Foam::Field<double> const&, unsigned char, int) const at ??:?

#5  Foam::smoothSolver::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ??:?

#6  Foam::fvMatrix<double>::solveSegregated(Foam::dictionary const&) at ??:?

#7  Foam::fvMatrix<double>::solve(Foam::dictionary const&) at ??:?

#8  Foam::fvMatrix<double>::solve() at ??:?

#9  

 at ??:?

#10  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"

#11  

 at ??:?

Does anyone have an idea what is wrong here? From what I can read from the error message it complains about signs or possibly division by zero/square root of a negative number,... but I can't figure out the origin of the problem.

Would be very glad to get any help on this! Thanks in advance.

Regards, Alex

Hello Alex,

Perhaps I can help you. I am not sure why your method is wrong, but I can push you in the right direction.

twoPhaseEulerFoam makes use of two phases, so your scalar has to be defined for one or both of the phases, i.e. in phaseModel.C and .H, you have to add your scalar similar to 'alpha'. So for example you add an entry 'phi', which should be and IOobject named "psi", and depending on what you want with regard to BCs/ICs you can make it READ_IF_PRESENT, NO_READ, AUTO_WRITE, etc... Don't forget to declare psi in the phaseModel.H file!

After having done this you can add your scalar to one or both of the phases in twoPhaseEulerFoam, look at all the examples in the beginning of the createFields.H file. For example, you can call your psi for phase 1 with the command phase1.psi. Now you should be able to use this psi in your scalar transport equation.

I hope this makes sense to you, if you have any questions let me know.

Kind regards,
Ramon

Hi Ramon,
thanks for the help.
Does it make a difference if I declare psi like in my first post compared to the implementation through phaseModel.C/H as you mentioned?
I'm trying that approach at the moment but I'm curious what could be the difference.

Thanks, Alex

EDIT: Actually, the results are the same :-/

You could first of all try leaving out any turbulence terms and other complexities, try just a simple convection/diffusion equation, and build it up from there. I know its a little trivial but it that way you can find the part of the equation that is destroying your solver.

K.R.
Ramon

Hi,
may I ask why are you using alphaEff coefficient for diffusion?

fvc::interpolate(thermo1.alphaEff(phase1.turbulenc e().mut()))

This coefficient is the thermal diffusion rate (ratio between kinematic visocsity and Prandtl number). Is this what you really want?

Can you also try do add before the solve the following line?

fvOptions.constrain(psiEqn);

Ciao
Mattia

Try the PBiCG solver with none as preconditioner.

From the solver output I figure that there are zones in your domain with alpha=0. In these zones the terms of your transport equation are multiplied by zero. Solving the linear equation system involves dividing terms in some cases.

Since the floating point exception happens during the smoothing operation, my guess is, that you run into numerical trouble due to the zero-alpha values in parts of your domain.

I found a much easier option using function objects, that doensn't require you to modify the solver.

please find it in this thread

http://www.cfd-online.com/Forums/ope...tml#post574980

Great post vkoppejan, this could be very useful for obtaining some quick results. I will definitely add the post to my ever growing list of OpenFOAM tips & tricks.

Keep in mind that this only works for constant isotropic diffusion type problems. So for simple cases this could be the perfect solution, for more advances problems this seems to lack flexibility.

Kind regards,
Ramon

Hi Ramon,

Your right, it's a very simple solution but the function object can easily adapted to cope with other types of diffusion problems. You just need to make sure the proper files are included and that the makefiles are adapted accordingly.

Glad to have helped you.

Cheers,

Victor