Introduce a drag model
 

I wanted to introduce a drag model call "Tam" into the twoPhaseEulerFoam solver in OF v1.7.1. However, the solver made a complain. I added Tam.C to "files" located at the following folder. Please let me know what mistakes I did in Tam.C.

"run/multiphase/twoPhaseEulerFoamMod/interfacialModels/Make/"

====Tam.C========
...

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

Foam::tmp<Foam::volScalarField> Foam::Tam::K
(
const volScalarField& Ur
) const
{
volScalarField beta = max(scalar(1) - alpha_, scalar(1.0e-6));

forAll(alpha_, celli)
{
if(alpha_[celli] > 0.6666)
{
functAlpha[celli] = 1.0e6;
}
}

return 4.5/(beta*alpha_ + scalar(1.0e-6))*phaseb_.nu()/pow(phasea_.d(),2.0)*functAlpha;


}

================

Below is what I see on my screen.

=====
....

Time = 0.002

DILUPBiCG: Solving for alpha, Initial residual = 0.0021426, Final residual = 5.45235e-12, No Iterations 3
Dispersed phase volume fraction = 0.36753 Min(alpha) = -0.0272179 Max(alpha) = 1.00084
DILUPBiCG: Solving for alpha, Initial residual = 8.90552e-05, Final residual = 1.62296e-12, No Iterations 3
Dispersed phase volume fraction = 0.36753 Min(alpha) = -0.0291739 Max(alpha) = 1.00075
#0 Foam::error::printStack(Foam::Ostream&) in "/home/jeong/OpenFOAM/openfoam171/lib/linuxGccDPOpt/libOpenFOAM.so"
#1 Foam::sigFpe::sigFpeHandler(int) in "/home/jeong/OpenFOAM/openfoam171/lib/linuxGccDPOpt/libOpenFOAM.so"
#2 Uninterpreted:
#3 Foam::sqrt(Foam::Field<double>&, Foam::UList<double> const&) in "/home/jeong/OpenFOAM/openfoam171/lib/linuxGccDPOpt/libOpenFOAM.so"
#4 Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > Foam::sqrt<Foam::fvPatchField, Foam::volMesh>(Foam::tmp<Foam::GeometricField<doub le, Foam::fvPatchField, Foam::volMesh> > const&) in "/home/jeong/OpenFOAM/openfoam171/lib/linuxGccDPOpt/libEulerianInterfacialModels.so"
#5 Foam::Tam::K(Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&) const in "/home/jeong/OpenFOAM/openfoam171/lib/linuxGccDPOpt/libEulerianInterfacialModels.so"
#6
in "/home/jeong/OpenFOAM/openfoam171/applications/bin/linuxGccDPOpt/twoPhaseEulerFoamMod"
#7 __libc_start_main in "/lib/tls/i686/cmov/libc.so.6"
#8
in "/home/jeong/OpenFOAM/openfoam171/applications/bin/linuxGccDPOpt/twoPhaseEulerFoamMod"
Floating point exception

=====

There is somewhere an error in calculating a square root. I suppose the value of K becomes negative in your calculation. Double check your formula or consider bounding the value of K: max(SMALL,4.5/...)

Edit: You can see that your minimum volume fraction is negative (Min(alpha) = -0.0272179), This is not physical but if it works with everything else I would replace all 'alpha' in your code with 'alphaBounded' and add above:

volScalarField alphaBounded = min(1.0,max(0.0,alpha));

Thank you for your comment.