[ANacc]

I have run into an issue while trying to apply the Schaeffer frictional stress model to the RAS/fluidisedBed tutorial for twoPhaseEulerFoam.

Error log:

#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::sigFpe::sigHandler(int) at ??:?
#2   in "/lib/x86_64-linux-gnu/libc.so.6"
#3   in "/lib/x86_64-linux-gnu/libm.so.6"
#4  pow in "/lib/x86_64-linux-gnu/libm.so.6"
#5  Foam::pow(Foam::Field<double>&, Foam::UList<double> const&, double const&) at ??:?
#6  Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > Foam::pow<Foam::fvPatchField, Foam::volMesh>(Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&, Foam::dimensioned<double> const&) at ??:?
#7  Foam::dragModels::WenYu::CdRe() const at ??:?
#8  Foam::dragModels::GidaspowErgunWenYu::CdRe() const at ??:?
#9  Foam::dragModel::K() const at ??:?
#10  Foam::BlendedInterfacialModel<Foam::dragModel>::K() const at ??:?
#11  Foam::twoPhaseSystem::dragCoeff() const at ??:?
#12  
 at ??:?
#13  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#14  
 at ??:?

Run log:

Create time

Create mesh for time = 0


Reading g
Creating twoPhaseSystem

Selecting thermodynamics package 
{
    type            heRhoThermo;
    mixture         pureMixture;
    transport       const;
    thermo          hConst;
    equationOfState rhoConst;
    specie          specie;
    energy          sensibleInternalEnergy;
}

Calculating face flux field phi.particles
Selecting diameterModel for phase particles: constant
Selecting turbulence model type RAS
Selecting RAS turbulence model kineticTheory
Selecting viscosityModel Gidaspow
Selecting conductivityModel Gidaspow
Selecting radialModel SinclairJackson
Selecting granularPressureModel Lun
Selecting frictionalStressModel Schaeffer
kineticTheoryCoeffs
{
    equilibrium     off;
    e               0.8;
    alphaMax        0.62;
    alphaMinFriction 0.5;
    residualAlpha   0.0001;
    viscosityModel  Gidaspow;
    conductivityModel Gidaspow;
    granularPressureModel Lun;
    frictionalStressModel Schaeffer;
    radialModel     SinclairJackson;
    SchaefferCoeffs
    {
        Fr              0.05;
        eta             2;
        p               5;
        phi             28.5;
    }
}

Selecting thermodynamics package 
{
    type            heRhoThermo;
    mixture         pureMixture;
    transport       const;
    thermo          hConst;
    equationOfState perfectGas;
    specie          specie;
    energy          sensibleInternalEnergy;
}

Calculating face flux field phi.air
Selecting diameterModel for phase air: constant
Selecting turbulence model type RAS
Selecting RAS turbulence model kEpsilon
kEpsilonCoeffs
{
    Cmu             0.09;
    C1              1.44;
    C2              1.92;
    C3              0;
    sigmak          1;
    sigmaEps        1.3;
}

Selecting default blending method: none
Selecting dragModel for (particles in air): GidaspowErgunWenYu
Selecting swarmCorrection for (particles in air): none
Selecting swarmCorrection for (particles in air): none
Selecting swarmCorrection for (particles in air): none
Selecting heatTransferModel for (particles in air): RanzMarshall
Calculating field DDtU1 and DDtU2

Creating field dpdt

Creating field kinetic energy K

No MRF models present

No finite volume options present

Courant Number mean: 0.00725 max: 0.01
Max Ur Courant Number = 0.01

PIMPLE: no residual control data found. Calculations will employ 3 corrector loops


Starting time loop

fieldAverage fieldAverage1:
    Starting averaging at time 0

Courant Number mean: 0.00725 max: 0.01
Max Ur Courant Number = 0.01
Time = 0.0002

PIMPLE: iteration 1
MULES: Solving for alpha.particles
MULES: Solving for alpha.particles
smoothSolver:  Solving for alpha.particles, Initial residual = 4.69513e-17, Final residual = 1.54587e-17, No Iterations 1
alpha.particles volume fraction = 0.275008  Min(alpha1) = 0  Max(alpha1) = 0.55
smoothSolver:  Solving for e.particles, Initial residual = 1, Final residual = 5.33467e-08, No Iterations 2
smoothSolver:  Solving for e.air, Initial residual = 1, Final residual = 2.54982e-07, No Iterations 3
GAMG:  Solving for p, Initial residual = 1, Final residual = 0.00478867, No Iterations 2
PIMPLE: iteration 2
MULES: Solving for alpha.particles
MULES: Solving for alpha.particles
smoothSolver:  Solving for alpha.particles, Initial residual = 0.999993, Final residual = 8.37276e-05, No Iterations 1000
alpha.particles volume fraction = 2470.01  Min(alpha1) = -66.6511  Max(alpha1) = 465976
smoothSolver:  Solving for e.particles, Initial residual = 1, Final residual = 7.81874e-07, No Iterations 17
smoothSolver:  Solving for e.air, Initial residual = 1, Final residual = 9.81972e-07, No Iterations 53
GAMG:  Solving for p, Initial residual = 1, Final residual = 6.46513e-09, No Iterations 1
PIMPLE: iteration 3
MULES: Solving for alpha.particles
MULES: Solving for alpha.particles
smoothSolver:  Solving for alpha.particles, Initial residual = 4.35647e-39, Final residual = 1.55142e-39, No Iterations 1
alpha.particles volume fraction = 1.27164e+83  Min(alpha1) = -1.66874e+97  Max(alpha1) = 1.66874e+97

All that I have changed from the original tutorial files is the kineticTheoryCoeffs in the turbulenceProperties.particles file:

kineticTheoryCoeffs
    {
        equilibrium             off;

        e                       0.8;
        alphaMax                0.62;
        alphaMinFriction        0.5;
        residualAlpha           1e-4;

        viscosityModel          Gidaspow;
        conductivityModel       Gidaspow;
        granularPressureModel   Lun;
        frictionalStressModel   Schaeffer;
        radialModel             SinclairJackson;

        SchaefferCoeffs
        {
            Fr                      0.05;
            eta                     2;
            p                       5;
            phi                     28.5;
        }
    }

Has anyone had success applying the Schaeffer frictional stress model?

[ANacc]

An update to my debugging question

I am now running into a different issue when trying to apply the Schaeffer boundary condition. Error log posted below. Any guidance is greatly appreciated. Again, version is OF 2.3.0

--> FOAM FATAL ERROR:
Maximum number of iterations exceeded

    From function thermo<Thermo, Type>::T(scalar f, scalar T0, scalar (thermo<Thermo, Type>::*F)(const scalar) const, scalar (thermo<Thermo, Type>::*$
    in file /home/sergio/rpmbuild/BUILD/OpenFOAM-2.3.0/src/thermophysicalModels/specie/lnInclude/thermoI.H at line 76.

FOAM aborting

#0  Foam::error::printStack(Foam::Ostream&) in "/shared/apps/openfoam/openfoam-2.3/install/OpenFOAM-2.3.0/platforms/linux64GccDPOpt/lib/libOpenFOAM.s$
#1  Foam::error::abort() in "/shared/apps/openfoam/openfoam-2.3/install/OpenFOAM-2.3.0/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
#2  Foam::species::thermo<Foam::hConstThermo<Foam::rhoConst<Foam::specie> >, Foam::sensibleInternalEnergy>::T(double, double, double, double (Foam::s$
#3  Foam::heRhoThermo<Foam::rhoThermo, Foam::pureMixture<Foam::constTransport<Foam::species::thermo<Foam::hConstThermo<Foam::rhoConst<Foam::specie> >$
#4  Foam::heRhoThermo<Foam::rhoThermo, Foam::pureMixture<Foam::constTransport<Foam::species::thermo<Foam::hConstThermo<Foam::rhoConst<Foam::specie> >$
#5
 in "/shared/apps/openfoam/openfoam-2.3/install/OpenFOAM-2.3.0/platforms/linux64GccDPOpt/bin/twoPhaseEulerFoam"
#6  __libc_start_main in "/lib64/libc.so.6"
#7
 in "/shared/apps/openfoam/openfoam-2.3/install/OpenFOAM-2.3.0/platforms/linux64GccDPOpt/bin/twoPhaseEulerFoam"

[wyldckat]

Greetings ANacc,

Please provide more information. Because from the 2 posts you've made, I can find at least 3 bug reports that seem to be related, but I can't figure out enough information to diagnose if it's the same problem.

The 3 bug reports I've seen to be related are:

• http://www.openfoam.org/mantisbt/view.php?id=1369
• http://www.openfoam.org/mantisbt/view.php?id=1363
• http://www.openfoam.org/mantisbt/view.php?id=1362

Best regards,
Bruno

[ANacc]

Bruno,

Thank you for the reply. Here is some more information.

To reproduce the first error message, I simply modified the turbulenceProperties.particles file of the RAS/fluidisedBed twoPhaseEulerFoam tutorial case (OF 2.3.0). I changed the frictionalStress model and coefficient dictionary from JohnsonJackson to Schaeffer. No other changes were made.

To reproduce the second error message, I changed the time step in the case I just described from 2e-4 to 1e-7. So Schaeffer frictional stress and smaller time step will reproduce the second error message.

I also have been testing a non-tutorial case in 2.3.x, again trying to apply the Schaeffer frictional stress. This will also replicate the second error message with a time step of 2e-4. I am happy to email it to you if it will help.

If there is any other information I can provide, please let me know.

Anthony

[wyldckat]

Hi Anthony,

With the additional information you've provided, I was able to reproduce the same exact errors, even with an up-to-date build of OpenFOAM 2.3.x (the latest bug fix version).

This is way beyond my expertise, but from I can figure out, when comparing the source code and tutorial cases in 2.2.x vs 2.3.x:

• The solver twoPhaseEulerFoam has changed considerably, where it now works with compressible modelling, instead of incompressible modelling. It uses dedicated alpha field calculations and several more changes were done. This means that it's somewhat hard to compare the two versions, to assess if there was something broken in the evolution.
• The library where the frictional models are has changed quite a bit as well. But when comparing the code side-by-side, at least the part related to friction, I can't spot any flaws in the code's evolution.
  • There is an issue in OpenFOAM 2.3.0, regarding the Schaeffer model, which was fixed in 2.3.x a few weeks after 2.3.0 was released: http://www.openfoam.org/mantisbt/view.php?id=1204 - Since I'm using 2.3.x, this problem shouldn't occur.

Honestly, I think it's best that you report this as a new bug: http://www.openfoam.org/bugs/
I say this because I went to 2.2.x and tried to adapt the similar tutorial named "bed" to use this frictional model and it worked without a problem.


One hypothesis I can think of would be that the Schaeffer friction model only works in an incompressible formulation. But I think that it's most likely something that went wrong with the evolution of the source code. Because the way a couple of the parameters are calculated in this specifci friction model, it makes me wonder how that even works. I'm referring to this part of the code:

Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::frictionalStressModels::Schaeffer::
frictionalPressure
(
    const volScalarField& alpha1,
    const dimensionedScalar& alphaMinFriction,
    const dimensionedScalar& alphaMax
) const
{
    return
        dimensionedScalar("1e24", dimensionSet(1, -1, -2, 0, 0), 1e24)
       *pow(Foam::max(alpha1 - alphaMinFriction, scalar(0)), 10.0);
}


Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::frictionalStressModels::Schaeffer::
frictionalPressurePrime
(
    const volScalarField& alpha1,
    const dimensionedScalar& alphaMinFriction,
    const dimensionedScalar& alphaMax
) const
{
    return
        dimensionedScalar("1e25", dimensionSet(1, -1, -2, 0, 0), 1e25)
       *pow(Foam::max(alpha1 - alphaMinFriction, scalar(0)), 9.0);
}

Multiplying 1e24 and 1e25 by something that is to the power of 10.0 and 9.0, seems to me to be something way too susceptible to numerical problems. Even though it worked fine in 2.2.x, it's possible that there was something triggered when using the compressible modelling, which is why I believe that only the OpenFOAM developers can properly diagnose this issue.

Best regards,
Bruno

[ANacc]

Bruno,

I really appreciate you looking into this.

I will submit a new bug report.

Again, thanks for the help.

Anthony

[dhuckaby]

Hi Anthony,

I agree with Bruno about the instability. The 0.55 initial condition for the volume fraction in the bed leads results in granular pressure (~2e13) for alphaMinFriction = 0.55. This leads to the large granular pressure gradient at the top of the bed and non-physical velocities.

I tried two fixes:
1) set alphaMinFriction to 0.6 - simulation ran to completion without problems
2) initialized the bed volume fraction to 0.45 - simulation crashed at ~0.15s. I tried to restart with lower time-steps, but it continued to crash.

Another issue with the tutorial is the thermal initial conditions. The gas is initialized to 300K while the solid is initialized to 600K. The gas quickly heats and expands creating a large velocity field at the beginning of the simulation.

3) Using a constant temperature field of 300K for both the gas and solid, initializing the volume fraction field to 0.45, and setting alphaMinFriction to 0.5, the simulation ran to completion (2s).

I hope this helps.
Dave

[ANacc]

Dave,

I was also able to run the tutorial case on both OF 2.3.0 and OF2.3.x with Schaeffer frictional stress after changing alphaMinFriction to 0.6. I tried modifying alphaMinFriction in a non-tutorial case that I have been looking at. alphaMax in this case is 0.58, and the only value of alphaMinFriction that would run successfully was 0.579.

I will add this information to the bug report.

Thank you for your help.

Anthony

[wyldckat]

Just for a future reference, the bug report in question is this one: http://www.openfoam.org/mantisbt/view.php?id=1378

[ANacc]

Just thought I would add a little more information to this topic:

I was curious as to why the Schaeffer model runs into the issues discussed in this thread while the JohnsonJackson model does not run into the same problems. I modified the Schaeffer code to calculate frictionalPressure and frictionalPressurePrime using the JohnsonJackson model's formulas:

frictionalPressure
(
const volScalarField& alpha1,
const dimensionedScalar& alphaMinFriction,
const dimensionedScalar& alphaMax
) const
{
return
Fr_*pow(max(alpha1 - alphaMinFriction, scalar(0)), eta_)
/pow(max(alphaMax - alpha1, scalar(5.0e-2)), p_);
}

frictionalPressurePrime
(
const volScalarField& alpha1,
const dimensionedScalar& alphaMinFriction,
const dimensionedScalar& alphaMax
) const
{
return Fr_*
(
eta_*pow(max(alpha1 - alphaMinFriction, scalar(0)), eta_ - 1.0)
*(alphaMax-alpha1)
+ p_*pow(max(alpha1 - alphaMinFriction, scalar(0)), eta_)
)/pow(max(alphaMax - alpha1, scalar(5.0e-2)), p_ + 1.0);
}

The RAS/fluidisedBed tutorial case ran with no issues when using this modified frictional stress model, so I switched the frictionalPressure back to the Schaeffer method. Again, the tutorial case ran with no issues.

So to me it seems like the error lies with the frictionalPressuePrime portion.

[mwaqas]

Hello everyone

The above discussion is very useful.
I am also working with fluidizedBed/RAS (twoPhaseEuler) model. I am using 2.3.0. In my simulation I want to apply partially slip Johanson Jackson boundary condition. my theta.Particles and U.particles boundary conditions (in 0 folder) are as follows

for theeta
walls
{
type zeroGradient;
}

for U.particles
walls
{
type fixedValue;
value uniform (0 0 0);
}

does it mean no slip boundary condition ? If it is so, how can I define partially slip boundary condition. Also where to define "Specularity coefficient"

Thank you

[RjwV]

Hello everyone,

I have also been testing the frictional stress models and playing around with the parameters. I have noticed that at some parameter combinations it seems the model runs into trouble, try for example the following Johnson and Jackson coefficients (taken from Ocone, 1993):

JohnsonJacksonCoeffs
{
    Fr    0.05;
    eta  2;
    p     3;
    phi   28.0;
}

For me this crashes at around 0.71 seconds I believe (if you want more details on this I can recheck...).

Anyone else who noticed this behavior?

Kind regards,
Ramon

[ANacc]

What case are you attempting to run with those settings?

[mwaqas]

Hello everone
I am working with fluidizedBed/RAS (twoPhaseEuler) model. I am using 2.3.0. my simulation is working fine. It does not crash, but I am facing a problem in mass conservation. When I apply following command in terminal window

patchIntegrate phi inlet/outlet

the inlet and outlet mass should remain conserve, but outlet mass is 5-6% higher than the inlet mass. One probable reason may be the solid particle elutriation, but it not the case, as no particle is moving outside the domain because of particle boundary condition. also the following command shows zero particle's flow rate in/out of the domain.

patchIntegrate phi.particles inlet/outlet

Anyone has idea what could be possible error

[RjwV]

Coming back to my post on the J&J frictional pressure Ocone settings... I have merely tried to do a simple fluidized bed simulation with different J&J parameter settings. It seems using Ocone's parameters in twoPhaseEulerFoam causes the solids hold up to go above maximum packing value, I have seen cases with alpha.particles of 0.84, but another case I have found alpha.particles to go up to 0.97-0.98. So depending on which radial distribution function you use, it may or may not crash the simulation (and it it does not crash you have unphysical results).

Not sure why these parameters have such a large impact on the simultion, but it is important to be careful with any changes made to them.

Kind regards,
Ramon

[RL-S]

I know this thread is three years old. However, the question was not really answered here and typically, such questions remain. I encountered the same problem today.

A good answer was in fact provided on the bug report site by will, which I feel warrants repeating here:

The Schaeffer model generates an exponentially increasing frictional pressure for particulate phase fractions above the specified alphaMinFriction of 0.5. The RAS setup also includes a particle pressure which allows the particulate fraction to reach 0.62. At this value, the frictional pressure generated by the Schaeffer model is very large, and presumably not physical. You've already found the answer to the instability; increasing alphaMinFriction.

Models work with some coefficients, and are unstable with others. If you have data which suggests the Schaeffer model is wrong, or you can spot something incorrect in the implementation, then please feel free to re-open this report. As that information is currently lacking, I'm marking this as closed.

NB: Re the CFD Online thread, frictionalPressurePrime is the derivative of the frictional pressure with respect to the volume fraction. The Schaeffer model is implemented consistently in this respect.

With the Schaeffer equation for granular pressure at closest packing being
10^25 * (alpha_max - alpha_minFriction)^9
it is very sensitive to user input. For the difference in parentheses being lower than 0.7, I didn't encounter any problems, though that value might not apply to other cases.