zhyang-dev/OpenFOAM-12
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

multiphaseInterFoam: Added test for contact angle in both phases in interface pair

Browse Source 
to ensure that the contact angle specification is used irrespective of which
phase it is specified in.  An error is reported if both phases of the interface
pair have a contact angle specification as the specifications might be
inconsistent.

Resolves bug-report https://bugs.openfoam.org/view.php?id=3688
This commit is contained in:
Henry Weller
2021-06-14 10:42:16 +01:00
parent 28745eca4b
commit de76426d86
2 changed files with 43 additions and 22 deletions
Download Patch File Download Diff File Expand all files Collapse all files

60
applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
View File

@ -391,12 +391,6 @@ Foam::tmp<Foam::surfaceScalarField> Foam::multiphaseMixture::nHatf
}
// Correction for the boundary condition on the unit normal nHat on
// walls to produce the correct contact angle.
// The dynamic contact angle is calculated from the component of the
// velocity on the direction of the interface, parallel to the wall.
void Foam::multiphaseMixture::correctContactAngle
(
const phase& alpha1,
@ -404,21 +398,42 @@ void Foam::multiphaseMixture::correctContactAngle
surfaceVectorField::Boundary& nHatb
) const
{
const volScalarField::Boundary& gbf
= alpha1.boundaryField();
const volScalarField::Boundary& a1bf = alpha1.boundaryField();
const volScalarField::Boundary& a2bf = alpha2.boundaryField();
const fvBoundaryMesh& boundary = mesh_.boundary();
forAll(boundary, patchi)
{
if (isA<alphaContactAngleFvPatchScalarField>(gbf[patchi]))
if
(
isA<alphaContactAngleFvPatchScalarField>(a1bf[patchi])
|| isA<alphaContactAngleFvPatchScalarField>(a2bf[patchi])
)
{
if
(
isA<alphaContactAngleFvPatchScalarField>(a1bf[patchi])
&& isA<alphaContactAngleFvPatchScalarField>(a2bf[patchi])
)
{
FatalErrorInFunction
<< "alphaContactAngle boundary condition "
"specified on patch " << boundary[patchi].name()
<< " for both " << alpha1.name() << " and " << alpha2.name()
<< nl << "which may be inconsistent."
<< exit(FatalError);
}
const alphaContactAngleFvPatchScalarField& acap =
refCast<const alphaContactAngleFvPatchScalarField>(gbf[patchi]);
isA<alphaContactAngleFvPatchScalarField>(a1bf[patchi])
? refCast<const alphaContactAngleFvPatchScalarField>(a1bf[patchi])
: refCast<const alphaContactAngleFvPatchScalarField>(a2bf[patchi])
;
vectorField& nHatPatch = nHatb[patchi];
vectorField AfHatPatch
const vectorField AfHatPatch
(
mesh_.Sf().boundaryField()[patchi]
/mesh_.magSf().boundaryField()[patchi]
@ -437,18 +452,19 @@ void Foam::multiphaseMixture::correctContactAngle
<< exit(FatalError);
}
bool matched = (tp.key().first() == alpha1.name());
const bool matched = (tp.key().first() == alpha1.name());
const scalar theta0 = degToRad(tp().theta0(matched));
scalar theta0 = degToRad(tp().theta0(matched));
scalarField theta(boundary[patchi].size(), theta0);
scalar uTheta = tp().uTheta();
const scalar uTheta = tp().uTheta();
// Calculate the dynamic contact angle if required
if (uTheta > small)
{
scalar thetaA = degToRad(tp().thetaA(matched));
scalar thetaR = degToRad(tp().thetaR(matched));
const scalar thetaA = degToRad(tp().thetaA(matched));
const scalar thetaR = degToRad(tp().thetaR(matched));
// Calculated the component of the velocity parallel to the wall
vectorField Uwall
@ -469,7 +485,7 @@ void Foam::multiphaseMixture::correctContactAngle
// Calculate Uwall resolved normal to the interface parallel to
// the interface
scalarField uwall(nWall & Uwall);
const scalarField uwall(nWall & Uwall);
theta += (thetaA - thetaR)*tanh(uwall/uTheta);
}
@ -477,9 +493,9 @@ void Foam::multiphaseMixture::correctContactAngle
// Reset nHatPatch to correspond to the contact angle
scalarField a12(nHatPatch & AfHatPatch);
const scalarField a12(nHatPatch & AfHatPatch);
scalarField b1(cos(theta));
const scalarField b1(cos(theta));
scalarField b2(nHatPatch.size());
@ -488,10 +504,10 @@ void Foam::multiphaseMixture::correctContactAngle
b2[facei] = cos(acos(a12[facei]) - theta[facei]);
}
scalarField det(1.0 - a12*a12);
const scalarField det(1.0 - a12*a12);
scalarField a((b1 - a12*b2)/det);
scalarField b((b2 - a12*b1)/det);
const scalarField a((b1 - a12*b2)/det);
const scalarField b((b2 - a12*b1)/det);
nHatPatch = a*AfHatPatch + b*nHatPatch;

5
applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.H
View File

@ -173,6 +173,11 @@ private:
const volScalarField& alpha2
) const;
//- Correction for the boundary condition on the unit normal nHat on
// walls to produce the correct contact angle.
//
// The dynamic contact angle is calculated from the component of the
// velocity on the direction of the interface, parallel to the wall.
void correctContactAngle
(
const phase& alpha1,