[konangsh]

Hil,
I am trying to implement Cox's dynamic contact angle model, and cant figure out what I am doing wrong! I am able to compile transportProperties, interfaceProperties, and interFoam in my local directory. But the solver crashes as soon as I try to run it. Can someone help me pls?

My code is this:

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

const Foam::scalar Foam::dynamicCoxAlphaContactAngleFvPatchScalarFiel d::convertToRad = Foam::constant::mathematical:i/180;
const Foam::scalar Foam::dynamicCoxAlphaContactAngleFvPatchScalarFiel d::convertToDeg = 180/Foam::constant::mathematical:i;
const Foam::scalar Foam::dynamicCoxAlphaContactAngleFvPatchScalarFiel d::theta0 = 90;

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

Foam::dynamicCoxAlphaContactAngleFvPatchScalarFiel d::
dynamicCoxAlphaContactAngleFvPatchScalarField
(
const fvPatch& p,
const DimensionedField<scalar, volMesh>& iF
)
:
alphaContactAngleFvPatchScalarField(p, iF),
thetaA_(0.0),
thetaR_(0.0)
{}

Foam::dynamicCoxAlphaContactAngleFvPatchScalarFiel d::
dynamicCoxAlphaContactAngleFvPatchScalarField
(
const dynamicCoxAlphaContactAngleFvPatchScalarField& gcpsf,
const fvPatch& p,
const DimensionedField<scalar, volMesh>& iF,
const fvPatchFieldMapper& mapper
)
:
alphaContactAngleFvPatchScalarField(gcpsf, p, iF, mapper),
thetaA_(gcpsf.thetaA_),
thetaR_(gcpsf.thetaR_)
{}

Foam::dynamicCoxAlphaContactAngleFvPatchScalarFiel d::
dynamicCoxAlphaContactAngleFvPatchScalarField
(
const fvPatch& p,
const DimensionedField<scalar, volMesh>& iF,
const dictionary& dict
)
:
alphaContactAngleFvPatchScalarField(p, iF, dict),
thetaA_(readScalar(dict.lookup("thetaA"))),
thetaR_(readScalar(dict.lookup("thetaR")))
{
Info<< "Dynamic contact angle (Cox) is being used\n" << endl;
evaluate();
}

Foam::dynamicCoxAlphaContactAngleFvPatchScalarFiel d::
dynamicCoxAlphaContactAngleFvPatchScalarField
(
const dynamicCoxAlphaContactAngleFvPatchScalarField& gcpsf
)
:
alphaContactAngleFvPatchScalarField(gcpsf),
thetaA_(gcpsf.thetaA_),
thetaR_(gcpsf.thetaR_)
{}

Foam::dynamicCoxAlphaContactAngleFvPatchScalarFiel d::
dynamicCoxAlphaContactAngleFvPatchScalarField
(
const dynamicCoxAlphaContactAngleFvPatchScalarField& gcpsf,
const DimensionedField<scalar, volMesh>& iF
)
:
alphaContactAngleFvPatchScalarField(gcpsf, iF),
thetaA_(gcpsf.thetaA_),
thetaR_(gcpsf.thetaR_)
{}


// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::tmp<Foam::scalarField>
Foam::dynamicCoxAlphaContactAngleFvPatchScalarFiel d::theta
(
const fvPatchVectorField& Up,
const fvsPatchVectorField& nHat
) const
{

const fvPatchField<scalar>& nup = patch().lookupPatchField<const::Foam::volScalarFie ld, scalar>("nu");
const fvPatchField<scalar>& rhop = patch().lookupPatchField<const::Foam::volScalarFie ld, scalar>("rho");

const vectorField nf(patch().nf());
// Calculated 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 interface
scalarField uwall(nWall & Uwall);

// Calculate local Capillary number
scalarField Ca = nup*rhop*mag(uwall)/sigmap;
Info<< "Capillary number (Cox) was evaluated\n" << endl;

scalarField thetaDp(patch().size(), (theta0));

/* 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 positiv for receding motion. */

forAll(uwall, pfacei)
{
if(uwall[pfacei] < 0.0)
{
thetaDp[pfacei] = pow((pow(thetaA_,3)+144*Ca[pfacei]),0.3333);
}
else if (uwall[pfacei] > 0.0)
{
thetaDp[pfacei] = pow((pow(thetaR_,3)+144*Ca[pfacei]),0.3333);
}
}

return (thetaDp);
}

[zangthanh]

Hi Konangsh,
I hope you solved your problems. if not, could you please show the error?