CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Running, Solving & CFD

multiphaseEulerFoam fails with zero volume fraction value

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   May 21, 2020, 18:10
Default multiphaseEulerFoam fails with zero volume fraction value
  #1
Member
 
Hasan Celik
Join Date: Sep 2016
Posts: 64
Rep Power: 9
PositronCascade is on a distinguished road
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!

Last edited by PositronCascade; May 21, 2020 at 20:45. Reason: added some info
PositronCascade is offline   Reply With Quote

Old   May 26, 2020, 05:50
Default
  #2
Senior Member
 
Join Date: Dec 2019
Location: Cologne, Germany
Posts: 355
Rep Power: 8
geth03 is on a distinguished road
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 is offline   Reply With Quote

Old   May 26, 2020, 20:28
Default
  #3
Member
 
Hasan Celik
Join Date: Sep 2016
Posts: 64
Rep Power: 9
PositronCascade is on a distinguished road
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.




Quote:
Originally Posted by geth03 View Post
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 is offline   Reply With Quote

Old   May 27, 2020, 03:04
Default
  #4
Senior Member
 
Join Date: Dec 2019
Location: Cologne, Germany
Posts: 355
Rep Power: 8
geth03 is on a distinguished road
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?
geth03 is offline   Reply With Quote

Old   May 27, 2020, 15:44
Default
  #5
Member
 
Hasan Celik
Join Date: Sep 2016
Posts: 64
Rep Power: 9
PositronCascade is on a distinguished road
Quote:
Originally Posted by geth03 View Post
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.


Code:
#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 ??:?
PositronCascade is offline   Reply With Quote

Old   May 28, 2020, 04:04
Default
  #6
Senior Member
 
Join Date: Dec 2019
Location: Cologne, Germany
Posts: 355
Rep Power: 8
geth03 is on a distinguished road
Quote:
Originally Posted by PositronCascade View Post
This is the terminal window when it crashes. Something turns into zero or NaN, I guess at it crashes.


Code:
#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.
geth03 is offline   Reply With Quote

Old   May 28, 2020, 12:17
Default
  #7
Member
 
Hasan Celik
Join Date: Sep 2016
Posts: 64
Rep Power: 9
PositronCascade is on a distinguished road
Quote:
Originally Posted by geth03 View Post
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
Code:
alpha_i * grad(p)
, I want to define it as
Code:
alpha_i * grad(p) + p * grad(alpha_i)
.

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

Code:
    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
Code:
UEqn == -fvc::grad(p)
But I couldn't get how multiphaseEulerFoam handles this solve function.

Here what I want to do,

Code:
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.
Attached Images
File Type: jpg copy1.jpg (121.5 KB, 18 views)
PositronCascade is offline   Reply With Quote

Old   May 29, 2020, 03:52
Default
  #8
Senior Member
 
Join Date: Dec 2019
Location: Cologne, Germany
Posts: 355
Rep Power: 8
geth03 is on a distinguished road
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 is offline   Reply With Quote

Old   June 5, 2020, 03:59
Default
  #9
Senior Member
 
Join Date: Dec 2019
Location: Cologne, Germany
Posts: 355
Rep Power: 8
geth03 is on a distinguished road
Quote:
Originally Posted by PositronCascade View Post
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
Code:
alpha_i * grad(p)
, I want to define it as
Code:
alpha_i * grad(p) + p * grad(alpha_i)
.

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

Code:
    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
Code:
UEqn == -fvc::grad(p)
But I couldn't get how multiphaseEulerFoam handles this solve function.

Here what I want to do,

Code:
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?
geth03 is offline   Reply With Quote

Old   June 5, 2020, 04:26
Default
  #10
Member
 
Hasan Celik
Join Date: Sep 2016
Posts: 64
Rep Power: 9
PositronCascade is on a distinguished road
Quote:
Originally Posted by geth03 View Post
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?
PositronCascade is offline   Reply With Quote

Old   June 5, 2020, 06:53
Default
  #11
Senior Member
 
Join Date: Dec 2019
Location: Cologne, Germany
Posts: 355
Rep Power: 8
geth03 is on a distinguished road
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.
geth03 is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
SU2-7.0.1 on ubuntu 18.04 hyunko SU2 Installation 7 March 16, 2020 05:37
multiphaseEulerFoam (OF2.3.0) : Courant number explodes when running in parallel Mehrez OpenFOAM Running, Solving & CFD 10 May 18, 2016 12:44
Stuck in a Rut- interDyMFoam! xoitx OpenFOAM Running, Solving & CFD 14 March 25, 2016 08:09
mpirun interFoam very stable --> then blows up ghadab OpenFOAM Running, Solving & CFD 3 October 27, 2013 11:34
the Eulerian model and the lower volume fraction hx li FLUENT 2 October 27, 2005 07:17


All times are GMT -4. The time now is 11:59.