[yang_ab]

Dear foamers,
I implement this drag model in openfoam10 solver multiphaseEulerFoam,
C_D=\frac{4.4}{Re^{1.42}}(1+0.196Ar^{0.381})(1+Eo) ^{0.78}Mo^{-0.15}
And this is my code:

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Copyright (C) 2014-2022 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#include "Hosen.H"
#include "aspectRatioModel.H"
#include "addToRunTimeSelectionTable.H"

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

namespace Foam
{
namespace dragModels
{
    defineTypeNameAndDebug(Hosen, 0);
    addToRunTimeSelectionTable(dragModel, Hosen, dictionary);
}
}


// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

Foam::dragModels::Hosen::Hosen
(
    const dictionary& dict,
    const phaseInterface& interface,
    const bool registerObject
)
:
    dispersedDragModel(dict, interface, registerObject),
    residualRe_("residualRe", dimless, dict),
    residualMo_("residualMo", dimless, dict)

{}


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

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


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

Foam::tmp<Foam::volScalarField>
Foam::dragModels::Hosen::CdRe() const
{
     volScalarField Eo(interface_.Eo());
     volScalarField Mo(max(interface_.Mo(), residualMo_));
     volScalarField Re(max(interface_.Re(), residualRe_));
     volScalarField Ar(interface_.Ar());


    return
        (4.4/pow(Re,1.42)*(1+0.196*pow(Ar,0.381))*pow((1+Eo),0.78)*pow(Mo,-0.15));

}


// ************************************************************************* //

There were no problems during my compilation, and I verified that as long as residualRe is set to 5, the value of Cd will not exceed the boundary. However, when using this model to calculate examples, a floating point overflow error will occur, and the error is as follows:

Create time

Create mesh for time = 0


PIMPLE: No convergence criteria found


PIMPLE: No corrector convergence criteria found
        Calculations will do 3 corrections


PIMPLE: Operating solver in transient mode with 3 outer correctors



Reading g

Reading hRef
Creating phaseSystem

Selecting phaseSystem basicMultiphaseSystem
Selecting phaseModel for air: purePhaseModel
Selecting diameterModel for phase air: constant
Selecting thermodynamics package 
{
    type            heRhoThermo;
    mixture         pureMixture;
    transport       const;
    thermo          hConst;
    equationOfState perfectGas;
    specie          specie;
    energy          sensibleInternalEnergy;
}

Calculating face flux field phi.air
Selecting turbulence model type laminar
Selecting laminar stress model Stokes
Selecting thermophysical transport type laminar
Selecting default laminar thermophysical transport model unityLewisFourier
Selecting phaseModel for water: purePhaseModel
Selecting diameterModel for phase water: constant
Selecting thermodynamics package 
{
    type            heRhoThermo;
    mixture         pureMixture;
    transport       const;
    thermo          eConst;
    equationOfState rPolynomial;
    specie          specie;
    energy          sensibleInternalEnergy;
}

Calculating face flux field phi.water
Selecting turbulence model type laminar
Selecting laminar stress model Stokes
Selecting thermophysical transport type laminar
Selecting default laminar thermophysical transport model unityLewisFourier
No MRF models present

Selecting surfaceTensionModel for air_water: constant
Selecting dragModel blending method for air_water: continuous
Selecting dragModel for air_dispersedIn_water: Hosen
Calculating field g.h

Reading field p_rgh

Creating fvConstraints from "system/fvConstraints"

Selecting finite volume constraint type limitPressure
    Name: limitp
    min 10000

Courant Number mean: 0.000375 max: 0.000375

Starting time loop

Courant Number mean: 0.000375 max: 0.000375
Time = 0.005s

PIMPLE: Iteration 1
MULES: Solving for alpha.air
air fraction, min, max = 1.25e-06 0 0.5
MULES: Solving for alpha.air
air fraction, min, max = 2.5e-06 0 0.5
Constructing momentum equations
#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::symGaussSeidelSmoother::smooth(Foam::word const&, Foam::Field<double>&, Foam::lduMatrix const&, Foam::Field<double> const&, Foam::FieldField<Foam::Field, double> const&, Foam::UPtrList<Foam::lduInterfaceField const> const&, unsigned char, int) at ??:?
#4  Foam::symGaussSeidelSmoother::smooth(Foam::Field<double>&, Foam::Field<double> const&, unsigned char, int) const at ??:?
#5  Foam::smoothSolver::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ??:?
#6  Foam::fvMatrix<double>::solveSegregated(Foam::dictionary const&) at ??:?
#7  Foam::fvMatrix<double>::solve(Foam::dictionary const&) in "/home/dyfluid/OpenFOAM/OpenFOAM-10/platforms/linux64GccDPInt32Opt/bin/multiphaseEulerFoam"
#8  Foam::fvMatrix<double>::solve() in "/home/dyfluid/OpenFOAM/OpenFOAM-10/platforms/linux64GccDPInt32Opt/bin/multiphaseEulerFoam"
#9  ? in "/home/dyfluid/OpenFOAM/OpenFOAM-10/platforms/linux64GccDPInt32Opt/bin/multiphaseEulerFoam"
#10  ? in "/lib/x86_64-linux-gnu/libc.so.6"
#11  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#12  ? in "/home/dyfluid/OpenFOAM/OpenFOAM-10/platforms/linux64GccDPInt32Opt/bin/multiphaseEulerFoam"
Floating point exception (core dumped)

What can I do to avoid this ,many thanks

[amirhosseintaran]

Hello,

There could be a couple of parameters causing this,
To test the worst scenario,

Could you re-compile the code without having Re at the denominator?
If that worked, it means that Re is messing it up. If not, it means that this drag value is causing problems elsewhere in the model.


Bests,
Amir.

[geth03]

1. probably division by 0 in Re, like @amirhosseintaran suggests. you can do something like min(Re, 1E-19).

2. keep in mind that in the code you have Cd*Re! but your formula is for Cd.

[yang_ab]

Thank you so much for your advice ,I missed Re in my formula ,after adding *Re in the formula, it all works well