[joshmccraney]

Hi FOAMers!

I am on OpenFOAM 6. I found an alphaContactAngle model that I'm trying to implement (Kistler contact angle model), the file called dynamicKistlerAlphaContactAngleFvPatchScalarField. H (and .C). It compiles, but issues the warning:

alphaContactAngle/dynamicDavisHockingAlphaContactAngle/dynamicDavisHockingAlphaContactAngleFvPatchScalarField.C:160:33: warning: unused variable ‘sigmap’ [-Wunused-variable]
     const fvPatchField<scalar>& sigmap =
                                 ^~~~~~
alphaContactAngle/dynamicDavisHockingAlphaContactAngle/dynamicDavisHockingAlphaContactAngleFvPatchScalarField.C:163:33: warning: unused variable ‘mup’ [-Wunused-variable]
     const fvPatchField<scalar>& mup =
                                 ^~~
alphaContactAngle/dynamicDavisHockingAlphaContactAngle/dynamicDavisHockingAlphaContactAngleFvPatchScalarField.C:202:12: warning: unused variable ‘InvHoffFuncThetaAroot’ [-Wunused-variable]
     scalar InvHoffFuncThetaAroot = RRInvHoffFuncThetaA.root(0,65);
            ^~~~~~~~~~~~~~~~~~~~~
alphaContactAngle/dynamicDavisHockingAlphaContactAngle/dynamicDavisHockingAlphaContactAngleFvPatchScalarField.C:205:12: warning: unused variable ‘InvHoffFuncThetaRroot’ [-Wunused-variable]
     scalar InvHoffFuncThetaRroot = RRInvHoffFuncThetaR.root(0,65);

In order to use this model I have to use a different solver called kistlerInterFoam (basically interFoam with some small tweaks). This compiles fine with no warnings. When I run the case, I get the error

#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::dynamicKistlerAlphaContactAngleFvPatchScalarField::HoffmanFunction(double const&) const at ??:?
#6  Foam::dynamicKistlerAlphaContactAngleFvPatchScalarField::theta(Foam::fvPatchField<Foam::Vector<double> > const&, Foam::fvsPatchField<Foam::Vector<double> > const&) const at ??:?
#7  Foam::interfaceProperties::correctContactAngle(Foam::GeometricField<Foam::Vector<double>, Foam::fvsPatchField, Foam::surfaceMesh>::Boundary&, Foam::GeometricField<Foam::Vector<double>, Foam::fvsPatchField, Foam::surfaceMesh>::Boundary const&) const at ??:?
#8  Foam::interfaceProperties::calculateK() at ??:?
#9  Foam::interfaceProperties::interfaceProperties(Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&, Foam::IOdictionary const&) at ??:?
#10  Foam::immiscibleIncompressibleTwoPhaseMixture::immiscibleIncompressibleTwoPhaseMixture(Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> const&) at ??:?
#11  ? in "/opt/openfoam6/platforms/linux64GccDPInt32Opt/bin/kistlerInterFoam"
#12  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#13  ? in "/opt/openfoam6/platforms/linux64GccDPInt32Opt/bin/kistlerInterFoam"
./run.sh: line 9:  8374 Floating point exception(core dumped) kistlerInterFoam

Here's what I know

Line #1 contains an Fpe error, which google implies is division by zero.
Line #4 includes pow, which appears only once in the dynamicKistlerAlphaContactAngleFvPatchScalarField. C and .H files
Line #5 references a HoffmanFunction, which is also contained in dynamicKistlerAlphaContactAngleFvPatchScalarField. H and .C files.

Below is the associated .C file, where I highlight in green the Hoffman function and pow components and immediately before that a division step by sigmap, which again threw a warning in compiling.

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright held by original author
     \\/     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 2 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, write to the Free Software Foundation,
    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA

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

#include "dynamicKistlerAlphaContactAngleFvPatchScalarField.H"
#include "addToRunTimeSelectionTable.H"
#include "fvPatchFieldMapper.H"
#include "fvPatchFields.H"
#include "volMesh.H"

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

namespace Foam
{

// * * * * * * * * * * * * * * * Static Member Data  * * * * * * * * * * * * //

const scalar dynamicKistlerAlphaContactAngleFvPatchScalarField::convertToDeg =
    180.0/constant::mathematical::pi;

const scalar dynamicKistlerAlphaContactAngleFvPatchScalarField::convertToRad =
    constant::mathematical::pi/180.0;

const scalar dynamicKistlerAlphaContactAngleFvPatchScalarField::theta0 = 90.0;

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

dynamicKistlerAlphaContactAngleFvPatchScalarField::
dynamicKistlerAlphaContactAngleFvPatchScalarField
(
    const fvPatch& p,
    const DimensionedField<scalar, volMesh>& iF
)
:
    alphaContactAngleFvPatchScalarField(p, iF),
    thetaA_(0.0),
    thetaR_(0.0),
    muName_("undefined"),
    sigmaName_("undefined")
{}

dynamicKistlerAlphaContactAngleFvPatchScalarField::
dynamicKistlerAlphaContactAngleFvPatchScalarField
(
    const dynamicKistlerAlphaContactAngleFvPatchScalarField& acpsf,
    const fvPatch& p,
    const DimensionedField<scalar, volMesh>& iF,
    const fvPatchFieldMapper& mapper
)
:
    alphaContactAngleFvPatchScalarField(acpsf, p, iF, mapper),
    thetaA_(acpsf.thetaA_),
    thetaR_(acpsf.thetaR_),
    muName_(acpsf.muName_),
    sigmaName_(acpsf.sigmaName_)
{}


dynamicKistlerAlphaContactAngleFvPatchScalarField::
dynamicKistlerAlphaContactAngleFvPatchScalarField
(
    const fvPatch& p,
    const DimensionedField<scalar, volMesh>& iF,
    const dictionary& dict
)
:
    alphaContactAngleFvPatchScalarField(p, iF),
    thetaA_(readScalar(dict.lookup("thetaA"))),
    thetaR_(readScalar(dict.lookup("thetaR"))),
    muName_(dict.lookup("muEffKistler")),
    sigmaName_(dict.lookup("sigmaKistler"))
{
    evaluate();
}


dynamicKistlerAlphaContactAngleFvPatchScalarField::
dynamicKistlerAlphaContactAngleFvPatchScalarField
(
    const dynamicKistlerAlphaContactAngleFvPatchScalarField& acpsf
)
:
    alphaContactAngleFvPatchScalarField(acpsf),
    thetaA_(acpsf.thetaA_),
    thetaR_(acpsf.thetaR_),
    muName_(acpsf.muName_),
    sigmaName_(acpsf.sigmaName_)
{}


dynamicKistlerAlphaContactAngleFvPatchScalarField::
dynamicKistlerAlphaContactAngleFvPatchScalarField
(
    const dynamicKistlerAlphaContactAngleFvPatchScalarField& acpsf,
    const DimensionedField<scalar, volMesh>& iF
)
:
    alphaContactAngleFvPatchScalarField(acpsf, iF),
    thetaA_(acpsf.thetaA_),
    thetaR_(acpsf.thetaR_),
    muName_(acpsf.muName_),
    sigmaName_(acpsf.sigmaName_)
{}


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

tmp<scalarField> dynamicKistlerAlphaContactAngleFvPatchScalarField::theta
(
    const fvPatchVectorField& Up,
    const fvsPatchVectorField& nHat
) const
{
    //eb - Lookup and return the patchField of dynamic viscosity of mixture
    //     and surface tension
/*    
    if((muName_ != "muEffKistler") || (sigmaName_ != "sigmaKistler"))
    {
        FatalErrorIn
        (
            "dynamicKistlerAlphaContactAngleFvPatchScalarField"
        )   << " muEffKistler or sigma set inconsitently, muEffKistler = "
            << muName_ << ", sigmaKistler = " << sigmaName_ << '.' << nl
            << "    Set both muEffKistler and sigmaKistler according to the "
            << "definition of dynamicKistlerAlphaContactAngle"
            << "\n    on patch " << this->patch().name()
            << " of field " << this->dimensionedInternalField().name()
            << " in file " << this->dimensionedInternalField().objectPath()
            << exit(FatalError);
    }
*/

    const fvPatchField<scalar>& sigmap =
        patch().lookupPatchField<volScalarField, scalar>(sigmaName_);

    const fvPatchField<scalar>& mup =
        patch().lookupPatchField<volScalarField, scalar>(muName_);

    vectorField nf = patch().nf();

    // Calculate the component of the velocity parallel to the wall
    vectorField Uwall = Up.patchInternalField() - Up;
    Uwall -= (nf & Uwall)*nf;

    // Find the direction of the interface parallel to the wall
    vectorField nWall = nHat - (nf & nHat)*nf;

    // Normalise nWall
    nWall /= (mag(nWall) + SMALL);

    // Calculate Uwall resolved normal to the interface parallel to
    // the wall
    scalarField uwall = nWall & Uwall;

    //eb - Calculate local Capillary number
    scalarField Ca = mup*mag(uwall)/sigmap;
 
    //eb - Instantiate function object InverseHoffmanFunction for thetaA and thetaR
    dynamicKistlerAlphaContactAngleFvPatchScalarField::InverseHoffmanFunction
    InvHoffFuncThetaA
    (
        convertToRad*thetaA_
    );

    dynamicKistlerAlphaContactAngleFvPatchScalarField::InverseHoffmanFunction
    InvHoffFuncThetaR
    (
        convertToRad*thetaR_
    );

   //eb - Calculate InverseHoffmanFunction for thetaA and thetaR using
    // RiddersRoot
    RiddersRoot RRInvHoffFuncThetaA(InvHoffFuncThetaA, 1.e-10);
    scalar InvHoffFuncThetaAroot = RRInvHoffFuncThetaA.root(0,65);

    RiddersRoot RRInvHoffFuncThetaR(InvHoffFuncThetaR, 1.e-10);
    scalar InvHoffFuncThetaRroot = RRInvHoffFuncThetaR.root(0,65);

    //eb - Calculate and return the value of contact angle on patch faces,
    //     a general approach: the product of Uwall and nWall is negative
    //     for advancing and positive for receding motion.
    //     thetaDp is initialized to 90 degrees corresponding to no wall
    //     adhesion
    scalarField thetaDp(patch().size(), convertToRad*theta0);
    forAll(uwall, pfacei)
    {
        if(uwall[pfacei] < 0.0)
        {
            //thetaDp[pfacei] = convertToRad*thetaA_;
            thetaDp[pfacei] = HoffmanFunction(   Ca[pfacei]
                                               + InvHoffFuncThetaAroot);
        }
        else if (uwall[pfacei] > 0.0)
        {
            //thetaDp[pfacei] = convertToRad*thetaR_;
            thetaDp[pfacei] = HoffmanFunction(   Ca[pfacei]
                                               - InvHoffFuncThetaRroot);
        }
    }
    
    return convertToDeg*thetaDp;
}


scalar dynamicKistlerAlphaContactAngleFvPatchScalarField::HoffmanFunction
(
    const scalar& x
) const
{
    return acos(1 - 2*tanh(5.16*pow(x/(1+1.31*pow(x,0.99)),0.706)));
}


void dynamicKistlerAlphaContactAngleFvPatchScalarField::write(Ostream& os) const
{
    fvPatchScalarField::write(os);
    os.writeKeyword("thetaA") << thetaA_ << token::END_STATEMENT << nl;
    os.writeKeyword("thetaR") << thetaR_ << token::END_STATEMENT << nl;
    os.writeKeyword("muEffKistler") << muName_ << token::END_STATEMENT << nl;
    os.writeKeyword("sigmaKistler") << sigmaName_ << token::END_STATEMENT << nl;
    writeEntry("value", os);
}


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

makePatchTypeField
(
    fvPatchScalarField,
    dynamicKistlerAlphaContactAngleFvPatchScalarField
);

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

} // End namespace Foam

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

In case this is relevant, the dynamicKistlerAlphaContactAngleFvPatchScalarField. H file is below, I've highlighted the HoffmanFunction part of the code in green.

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright held by original author
     \\/     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 2 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, write to the Free Software Foundation,
    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA

Class
    dynamicKistlerAlphaContactAngleFvPatchScalarField

Description
    Evaluation of the dynamic contact angle using the Hoffman function,
    based on Kistler (1993)

SourceFiles
    dynamicKistlerAlphaContactAngleFvPatchScalarField.C

Author
    Edin Berberovic 2008 - Updated by Alexander Rattner, 2012

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

#ifndef dynamicKistlerAlphaContactAngleFvPatchScalarField_H
#define dynamicKistlerAlphaContactAngleFvPatchScalarField_H

#include "alphaContactAngleFvPatchScalarField.H"
#include "mathematicalConstants.H"
#include "surfaceFields.H"
#include "RiddersRoot.H"

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

namespace Foam
{

/*---------------------------------------------------------------------------*\
                  Class dynamicAlphaContactAngleFvPatch Declaration
\*---------------------------------------------------------------------------*/

class dynamicKistlerAlphaContactAngleFvPatchScalarField
:
    public alphaContactAngleFvPatchScalarField
{
    // Private data

        //- Limiting advancing contact angle
        scalar thetaA_;

        //- Limiting receeding contact angle
        scalar thetaR_;

        //eb Name of the dynamic viscosity field
        word muName_;

        //eb Name of the surface tension
        word sigmaName_;

public:

    //eb - Function object class, inverse of the Hoffman function
    class InverseHoffmanFunction : public FuncBase
    {

        // Private data

            //- Equilibrium contact angle (advancing or receding), as parameter
            scalar thetaE_;

    public:

        // Constructors

            //- Construct from data
            InverseHoffmanFunction(const scalar& thetaE)
            :
                thetaE_(thetaE)
            {}

        // operator()
            scalar operator()(scalar fHI) const
            {
                return
                  pow(fHI,0.706)
                -   pow((1+1.31*pow(fHI,0.99)),0.706)
                   *(1.0/5.16)*atanh((1-cos(thetaE_))/2);
            }
    };

public:

    //- Runtime type information
    TypeName("dynamicKistlerAlphaContactAngle");

    // Static data members

        //eb - Conversion factor for radians into degrees
        static const scalar convertToDeg;

        //eb - Conversion factor for degrees into radians
        static const scalar convertToRad;

        //eb - Contact angle corresponding to zero contact line velocity
        static const scalar theta0;

    // Constructors

        //- Construct from patch and internal field
        dynamicKistlerAlphaContactAngleFvPatchScalarField
        (
            const fvPatch&,
            const DimensionedField<scalar, volMesh>&
        );

        //- Construct from patch, internal field and dictionary
        dynamicKistlerAlphaContactAngleFvPatchScalarField
        (
            const fvPatch&,
            const DimensionedField<scalar, volMesh>&,
            const dictionary&
        );

        //- Construct by mapping given
        //  dynamicKistlerAlphaContactAngleFvPatchScalarField
        //  onto a new patch
        dynamicKistlerAlphaContactAngleFvPatchScalarField
        (
            const dynamicKistlerAlphaContactAngleFvPatchScalarField&,
            const fvPatch&,
            const DimensionedField<scalar, volMesh>&,
            const fvPatchFieldMapper&
        );

        //- Construct as copy
        dynamicKistlerAlphaContactAngleFvPatchScalarField
        (
            const dynamicKistlerAlphaContactAngleFvPatchScalarField&
        );

        //- Construct and return a clone
        virtual tmp<fvPatchScalarField> clone() const
        {
            return tmp<fvPatchScalarField>
            (
                new dynamicKistlerAlphaContactAngleFvPatchScalarField(*this)
            );
        }

        //- Construct as copy setting internal field reference
        dynamicKistlerAlphaContactAngleFvPatchScalarField
        (
            const dynamicKistlerAlphaContactAngleFvPatchScalarField&,
            const DimensionedField<scalar, volMesh>&
        );

        //- Construct and return a clone setting internal field reference
        virtual tmp<fvPatchScalarField> clone
        (
            const DimensionedField<scalar, volMesh>& iF
        ) const
        {
            return tmp<fvPatchScalarField>
            (
                new dynamicKistlerAlphaContactAngleFvPatchScalarField(*this, iF)
            );
        }


    // Member functions

        //- Evaluate and return dynamic contact-angle
        virtual tmp<scalarField> theta
        (
            const fvPatchVectorField& Up,
            const fvsPatchVectorField& nHat
        ) const;

        //eb - Evaluate the value of the Hoffman function
        scalar HoffmanFunction
        (
            const scalar& x
        ) const;

        //- Write
        virtual void write(Ostream&) const;
};


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

} // End namespace Foam

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

#endif

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

Any help is greatly appreciated!

[frantov]

Hello, have you solved your problem?


Im also trying to implemen Kistler contact angle model in OF7 but I get this error.error:



‘Foam::RiddersRoot’ is not a template
194 | RiddersRoot<InverseHoffmanFunction> RRInvHoffFuncThetaA(InvHoffFuncThetaA, 1.e-10);



Im thinking the RiddlersRoot.H and C. files I have need to be updated. Did you updated this files?
Cheers
Fran

[frantov]

I found some updated .H .C files and did some conversion to OF7.
The library now compiled well..no errors
Now Im trying to link that library with the solver....
Cheers

[joshmccraney]

Originally Posted by frantov

I found some updated .H .C files and did some conversion to OF7.
The library now compiled well..no errors
Now Im trying to link that library with the solver....
Cheers

Sorry, I'm just now seeing these two posts. I ended up moving to OpenFOAM 4.X. If you got it working in OF7, could you post your scripts somewhere? Would be very helpful to see!