CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Programming & Development

applying Schaeffer frictional stress

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

Reply
 
LinkBack Thread Tools Display Modes
Old   August 4, 2014, 18:03
Default applying Schaeffer frictional stress
  #1
New Member
 
Join Date: Jun 2014
Posts: 15
Rep Power: 3
ANacc is on a distinguished road
I have run into an issue while trying to apply the Schaeffer frictional stress model to the RAS/fluidisedBed tutorial for twoPhaseEulerFoam.

Error log:
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   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:
Code:
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:

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

Old   August 14, 2014, 10:23
Default
  #2
New Member
 
Join Date: Jun 2014
Posts: 15
Rep Power: 3
ANacc is on a distinguished road
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

Code:
--> 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"
ANacc is offline   Reply With Quote

Old   August 14, 2014, 13:33
Default
  #3
Super Moderator
 
Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 8,526
Blog Entries: 34
Rep Power: 86
wyldckat is just really nicewyldckat is just really nicewyldckat is just really nicewyldckat is just really nice
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:
Best regards,
Bruno
wyldckat is offline   Reply With Quote

Old   August 14, 2014, 14:21
Default
  #4
New Member
 
Join Date: Jun 2014
Posts: 15
Rep Power: 3
ANacc is on a distinguished road
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
ANacc is offline   Reply With Quote

Old   August 14, 2014, 18:52
Default
  #5
Super Moderator
 
Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 8,526
Blog Entries: 34
Rep Power: 86
wyldckat is just really nicewyldckat is just really nicewyldckat is just really nicewyldckat is just really nice
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.
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:
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
wyldckat is offline   Reply With Quote

Old   August 14, 2014, 20:33
Default
  #6
New Member
 
Join Date: Jun 2014
Posts: 15
Rep Power: 3
ANacc is on a distinguished road
Bruno,

I really appreciate you looking into this.

I will submit a new bug report.

Again, thanks for the help.

Anthony
ANacc is offline   Reply With Quote

Old   August 19, 2014, 12:21
Default
  #7
New Member
 
David Huckaby
Join Date: Jul 2009
Posts: 21
Rep Power: 8
dhuckaby is on a distinguished road
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
dhuckaby is offline   Reply With Quote

Old   August 20, 2014, 09:02
Default
  #8
New Member
 
Join Date: Jun 2014
Posts: 15
Rep Power: 3
ANacc is on a distinguished road
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
ANacc is offline   Reply With Quote

Old   August 24, 2014, 08:34
Default
  #9
Super Moderator
 
Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 8,526
Blog Entries: 34
Rep Power: 86
wyldckat is just really nicewyldckat is just really nicewyldckat is just really nicewyldckat is just really nice
Just for a future reference, the bug report in question is this one: http://www.openfoam.org/mantisbt/view.php?id=1378
wyldckat is offline   Reply With Quote

Old   September 3, 2014, 17:18
Default
  #10
New Member
 
Join Date: Jun 2014
Posts: 15
Rep Power: 3
ANacc is on a distinguished road
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:

Code:
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_);
}
Code:
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.
ANacc is offline   Reply With Quote

Old   October 23, 2014, 14:20
Default
  #11
Member
 
Muhammad Waqas
Join Date: Jul 2014
Location: Germany
Posts: 53
Rep Power: 3
mwaqas is on a distinguished road
Send a message via Skype™ to 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
mwaqas is offline   Reply With Quote

Old   November 28, 2014, 11:42
Default
  #12
New Member
 
Ramon
Join Date: Feb 2014
Location: Eindhoven
Posts: 24
Rep Power: 3
RjwV is on a distinguished road
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):

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

Old   December 1, 2014, 17:21
Default
  #13
New Member
 
Join Date: Jun 2014
Posts: 15
Rep Power: 3
ANacc is on a distinguished road
What case are you attempting to run with those settings?

Last edited by ANacc; December 1, 2014 at 17:21. Reason: typo
ANacc is offline   Reply With Quote

Old   December 3, 2014, 14:50
Default
  #14
Member
 
Muhammad Waqas
Join Date: Jul 2014
Location: Germany
Posts: 53
Rep Power: 3
mwaqas is on a distinguished road
Send a message via Skype™ to 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
mwaqas is offline   Reply With Quote

Old   December 8, 2014, 07:15
Default
  #15
New Member
 
Ramon
Join Date: Feb 2014
Location: Eindhoven
Posts: 24
Rep Power: 3
RjwV is on a distinguished road
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
RjwV is offline   Reply With Quote

Reply

Thread Tools
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 On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Average frictional stress on surface weighted by elements area (Ansys Workbench Mech) Giulio89 Structural Mechanics 0 July 8, 2013 09:44
Reynolds Stress Models Jade M Main CFD Forum 0 April 21, 2010 16:38
Applying Shear Stress at a Boundary Condition Stephen CFX 3 March 27, 2007 18:55
shear stress scalar definitions? Novak Elliott CFX 0 April 6, 2003 01:45
What is the detail definition of wall shear stress zjm FLUENT 0 January 2, 2002 08:43


All times are GMT -4. The time now is 07:55.