[hester]

Hello,

I like to implement the Ishii-Zuber drag model for sparsely distributed fluid particles [1], as used in the CFX implementation [2]. I'm using OpenFoam-2.3.0.

Foam::tmp<Foam::volScalarField>Foam::dragModels::IshiiZuberExtended::CdRe() const
{
    volScalarField Eo(pair_.Eo());
    volScalarField d(pair_.dispersed().d());
    volScalarField rho(pair_.continuous().rho());
    volScalarField magUr(pair_.magUr());

    volScalarField CdSphere = (24.0 / pair_.Re() * (1 + 0.15 *
      pow(pair_.Re()+ROOTVSMALL , 0.687)));
    volScalarField CdEllipse = ((4.0/3.0) * pow(2 * Eo+ROOTVSMALL , 0.5));
    volScalarField CdCap = (8.0/3.0) * (pair_.Re()/pair_.Re()+ROOTVSMALL);       
    volScalarField Cd = max(CdSphere, min(CdEllipse, CdCap));

    return Cd * (pair_.Re()+ROOTVSMALL);

}

I added ROOTVSMALL because otherweise I would get a segmentation fault at the beginning of my simulation when U = 0 and threrefor Re = 0. I'm simulating a bubble column after Deen with U(initial) = 0.
The value for CdCap is actually just 8/3 but I didn't know how to create a volScalarField from the scalar value.

I get no compiling errors. But my simulation breaks at the first iteration with a floating point error and the following error message:

#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::divide(Foam::Field<double>&, double const&, Foam::UList<double> const&) at ??:?
#4  Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > Foam::operator/<Foam::fvPatchField, Foam::volMesh>(Foam::dimensioned<double> const&, Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > const&) at ??:?
#5  Foam::dragModels::IshiiZuberExtended::CdRe() const at ??:?
#6  Foam::dragModel::K() const at ??:?
#7  Foam::BlendedInterfacialModel<Foam::dragModel>::K() const at ??:?
#8  Foam::twoPhaseSystem::dragCoeff() const at ??:?
#9  
 at ??:?
#10  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#11  
 at ??:?
Gleitkomma-Ausnahme

I don't know what the problem is. The implementation seems fine to me. Maybe someone could give me hint where to look for a solution to my problem.

thanks
hester


[1] Ishii, M. and Zuber, N., “Drag Coefficient and Relative Velocity in Bubbly, Droplet or Particulate Flows”, AIChE J., 25, 843-855, 1979.
[2] http://www.arc.vt.edu/ansys_help/cfx_thry/i1304903.html, Eq. (5-51)-(5-53)

[GerhardHolzinger]

If you take a close look at the stack-trace you will see that the function Foam::divide was involved in throwing the FPE.
You divide by pair.Re() during your calculation of the drag coefficient.

Most probably, the initial value of pair.Re() evaluates to zero for your initial conditions.

That is why you find code like this max(pair_.Re(), residualRe_) in some of the drag models.

[hester]

Thank you for your answer. I saw terms like max(pair_.Re(), residualRe_) in different drag models but I never understood why it was used. I will try it in my drag model. Maybe it works.

[hester]

Hello Gerhard,

I changed my implementation using residual values for Re and Eo, as suggested. It works now. This is what it looks like:

Foam::dragModels::IshiiZuberExtended::IshiiZuberExtended
(
    const dictionary& dict,
    const phasePair& pair,
    const bool registerObject
)
:
    dragModel(dict, pair, registerObject),
    residualRe_("residualRe", dimless, dict.lookup("residualRe")),
    residualEo_("residualEo", dimless, dict.lookup("residualEo"))
{}


// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

Foam::dragModels::IshiiZuberExtended::~IshiiZuberExtended()
{}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

Foam::tmp<Foam::volScalarField> Foam::dragModels::IshiiZuberExtended::CdRe() const
{
    volScalarField Eo(max(pair_.Eo(), residualEo_));
    volScalarField d(pair_.dispersed().d());
    volScalarField rho(pair_.continuous().rho());
    volScalarField magUr(pair_.magUr());
    volScalarField Re(max(pair_.Re(), residualRe_));

    volScalarField CdSphere = (24.0 / Re * (1. + 0.15 * pow(Re , 0.687)));
    volScalarField CdEllipse = ((4.0/3.0) * pow(2. * Eo , 0.5));
    volScalarField CdCap = (8.0/3.0) * Re/Re;
    volScalarField Cd = max(CdSphere, min(CdEllipse, CdCap));

    return Cd * Re;
}

Regards,
hester