[PositronCascade]

Hello everyone. I use multiphaseEulerFoam solver to solve a cross-diffusion problem, so I just use it as a solver, and my problem is not related to multiphase capabilities of the solver directly. Therefore, assume a flow through a channel. At the left BC, phase1 = 0 and phase2 = 1 while at the right BC, phase1 = 1 and phase2 = 0. If I run the solver with these BCs, it fails. If I select the values such that 0.3 - 0.7 for instance, it works perfectly. How can I solve this issue? What might be the problem?

Anything greater than 0.01 works well with the code. Foam::multiphaseSystem::nearInterface() might be causing this issue. However, I couldn't find any way to go around this problem. I took interfacial compression value zero as I do not need any interface capturing.

Thanks for your help!

[geth03]

my multiphase solver also crashes when the conti phase volume fraction approaches 0. a look what the solver is solving at that point in time shows that the turbulence model causes the crash. i assume that something like a division by 0 takes place and the solver crashes as a result, but i didn't dig deeper into that.

do you use a turbulunce model? if yes, which one?
did you check what the last set of equations are before the crash?

[PositronCascade]

Hello. Well, I use laminar approach. I tried to use Eclipse Debugger to check what might be causing this problem, but I couldn't get any idea. How can I check the last set of equations are before the crash?

Thanks a lot.

Originally Posted by geth03

my multiphase solver also crashes when the conti phase volume fraction approaches 0. a look what the solver is solving at that point in time shows that the turbulence model causes the crash. i assume that something like a division by 0 takes place and the solver crashes as a result, but i didn't dig deeper into that.

do you use a turbulunce model? if yes, which one?
did you check what the last set of equations are before the crash?

[geth03]

if you run your solver you should see how the bash is printing stuff.
there is information like which residuals and how many iteration steps are calculated.
what is the last print before your error message?
can you post a screenshot, or copy the error message and the print before the error, lets say for the last time step?

[PositronCascade]

Originally Posted by geth03

if you run your solver you should see how the bash is printing stuff.
there is information like which residuals and how many iteration steps are calculated.
what is the last print before your error message?
can you post a screenshot, or copy the error message and the print before the error, lets say for the last time step?

This is the terminal window when it crashes. Something turns into zero or NaN, I guess at it crashes.

#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::tmp<Foam::GeometricField<Foam::Vector<double>, Foam::fvsPatchField, Foam::surfaceMesh> > Foam::operator*<Foam::Vector<double>, Foam::fvsPatchField, Foam::surfaceMesh>(Foam::tmp<Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> > const&, Foam::tmp<Foam::GeometricField<Foam::Vector<double>, Foam::fvsPatchField, Foam::surfaceMesh> > const&) at ??:?
#4  Foam::multiphaseSystem::nHatfv(Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&) const at ??:?
#5  Foam::multiphaseSystem::nHatf(Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&) const at ??:?
#6  Foam::multiphaseSystem::solveAlphas() at ??:?
#7  Foam::multiphaseSystem::solve() at ??:?

[geth03]

Originally Posted by PositronCascade

This is the terminal window when it crashes. Something turns into zero or NaN, I guess at it crashes.

#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::tmp<Foam::GeometricField<Foam::Vector<double>, Foam::fvsPatchField, Foam::surfaceMesh> > Foam::operator*<Foam::Vector<double>, Foam::fvsPatchField, Foam::surfaceMesh>(Foam::tmp<Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> > const&, Foam::tmp<Foam::GeometricField<Foam::Vector<double>, Foam::fvsPatchField, Foam::surfaceMesh> > const&) at ??:?
#4  Foam::multiphaseSystem::nHatfv(Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&) const at ??:?
#5  Foam::multiphaseSystem::nHatf(Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&) const at ??:?
#6  Foam::multiphaseSystem::solveAlphas() at ??:?
#7  Foam::multiphaseSystem::solve() at ??:?

can you also post please what is printed before this error message.
the stuff printed in the last time step.

[PositronCascade]

Originally Posted by geth03

can you also post please what is printed before this error message.
the stuff printed in the last time step.

Sure. I attached it. The problem is related to the pressure eqution, I mean the error. When it starts to solve pEqn, the solver crashes. I try to convert the pressure definition in multiphaseEulerFoam into partial pressure. Therefore, instead of just

alpha_i * grad(p)

, I want to define it as

alpha_i * grad(p) + p * grad(alpha_i)

.

So, I added this to the right hand side of UEqn:

    UEqns.set
    (
        phasei,
        new fvVectorMatrix
        (
            fvm::ddt(alpha, U)
          + fvm::div(phase.alphaPhi(), U)

          + (alpha/phase.rho())*fluid.Cvm(phase)*
            (
                fvm::ddt(U)
              + fvm::div(phase.phi(), U)
              - fvm::Sp(fvc::div(phase.phi()), U)
            )

          - fvm::laplacian(alpha*nuEff, U)
          - fvc::div
            (
                alpha*(nuEff*dev(T(fvc::grad(U))) /*- ((2.0/3.0)*I)*k*/),
                "div(Rc)"
            )
         ==
            (alpha/phase.rho())*fluid.Svm(phase)
          - p/phase.rho()*phase.gradAlpha() // ---> added this one
          - fvm::Sp
            (
                slamDampCoeff
               *max
                (
                    mag(U()) - maxSlamVelocity,
                    dimensionedScalar(dimVelocity, 0)
                )
               /pow(mesh.V(), 1.0/3.0),
                U
            )
        )
    );
    MRF.addAcceleration
    (
        alpha*(1 + (1/phase.rho())*fluid.Cvm(phase)),
        UEqns[phasei]
    );

But I am not sure if I it is correct to add it to this side as here I cannot see alpha_i * grad(p). For instance in simpleFoam solver, it defines UEqn so then it calls solve() function and solves

UEqn == -fvc::grad(p)

But I couldn't get how multiphaseEulerFoam handles this solve function.

Here what I want to do,

UEqn[i] == -alpha[i]*fvc::grad(p) - p * fvc::grad(alpha[i])

Maybe it is related to this part or it is related to zero constant BC value. May you explain how multiphaseEulerFoam takes that grad(p) in the solver? I hope you can help.

[geth03]

your simulation output is really strange,
i think your additional equation input screwed up something.

just look at your
Phase-sum volume fraction, min, max = ...
those values are unphysical.
when i run multiphase simulations i always get:
Phase-sum volume fraction, min, max = 1 1 1

[geth03]

Originally Posted by PositronCascade

Sure. I attached it. The problem is related to the pressure eqution, I mean the error. When it starts to solve pEqn, the solver crashes. I try to convert the pressure definition in multiphaseEulerFoam into partial pressure. Therefore, instead of just

alpha_i * grad(p)

, I want to define it as

alpha_i * grad(p) + p * grad(alpha_i)

.

So, I added this to the right hand side of UEqn:

    UEqns.set
    (
        phasei,
        new fvVectorMatrix
        (
            fvm::ddt(alpha, U)
          + fvm::div(phase.alphaPhi(), U)

          + (alpha/phase.rho())*fluid.Cvm(phase)*
            (
                fvm::ddt(U)
              + fvm::div(phase.phi(), U)
              - fvm::Sp(fvc::div(phase.phi()), U)
            )

          - fvm::laplacian(alpha*nuEff, U)
          - fvc::div
            (
                alpha*(nuEff*dev(T(fvc::grad(U))) /*- ((2.0/3.0)*I)*k*/),
                "div(Rc)"
            )
         ==
            (alpha/phase.rho())*fluid.Svm(phase)
          - p/phase.rho()*phase.gradAlpha() // ---> added this one
          - fvm::Sp
            (
                slamDampCoeff
               *max
                (
                    mag(U()) - maxSlamVelocity,
                    dimensionedScalar(dimVelocity, 0)
                )
               /pow(mesh.V(), 1.0/3.0),
                U
            )
        )
    );
    MRF.addAcceleration
    (
        alpha*(1 + (1/phase.rho())*fluid.Cvm(phase)),
        UEqns[phasei]
    );

But I am not sure if I it is correct to add it to this side as here I cannot see alpha_i * grad(p). For instance in simpleFoam solver, it defines UEqn so then it calls solve() function and solves

UEqn == -fvc::grad(p)

But I couldn't get how multiphaseEulerFoam handles this solve function.

Here what I want to do,

UEqn[i] == -alpha[i]*fvc::grad(p) - p * fvc::grad(alpha[i])

Maybe it is related to this part or it is related to zero constant BC value. May you explain how multiphaseEulerFoam takes that grad(p) in the solver? I hope you can help.

hey,
did you find a solution for you problem?

under phaseProperties i use now for blending the type linear and changed minFullyContinuousAlpha = 0.999 and minPartlyContinuousAlpha = 0.
my solver doesn't crash for now.

also in your equation you use the pressure p from the current time, but fvc::grad(p) uses the pressure from the previous time. can you acces the pressure from the previous time step too? will that help?

[PositronCascade]

Originally Posted by geth03

hey,
did you find a solution for you problem?

under phaseProperties i use now for blending the type linear and changed minFullyContinuousAlpha = 0.999 and minPartlyContinuousAlpha = 0.
my solver doesn't crash for now.

also in your equation you use the pressure p from the current time, but fvc::grad(p) uses the pressure from the previous time. can you acces the pressure from the previous time step too? will that help?

Hello, I couldn't find a solution yet. I don't know if there is any option such that phaseProperties in multiphaseEulerFoam, at least, I couldn't find it. May you give a bit more information?

I think fvc::grad(p) might be an issue; however, I don't know if there is any implicit approach for grad in Openfoam. Do you have any idea?

[geth03]

i see. i thought the multiphaseEulerFoam is kind of similar to reactingMultiphaseEulerFoam. in multiphaseEulerFoam it's the transportProperties. Never mind.

do you mind changing the solver to reactingMultiphaseEulerFoam?
there is blending option which you can adjust to your simulation.
i think multiphaseEulerFoam has a fixed blending, so your simulation crashes,
as it is the same when i use fixed blending in reactingMultiphaseEulerFoam.
you can adjust the settings of reactingMultiphaseEulerFoam so it is basically the same as multiphaseEulerFoam. for that you can copy one case from reactingMultiphaseEulerFoam and adjust it, look at the bubbleColumn-tutorial in twoPhaseEulerFoam for that.

regarding the fvc::grad(p): this is fine in the equation. i meant the standalone p in your equation, this one in the UEqn also should be some type of fvc, bc this should be taken from the previous step. so if you use p as stand alone this means the new p, which is not known at that moment, so we take the p from the previous step and therefore need some kind of fvc:: expression.