From de76426d86f076e36332c3554753f7bf5cd377c2 Mon Sep 17 00:00:00 2001 From: Henry Weller Date: Mon, 14 Jun 2021 10:42:16 +0100 Subject: [PATCH] multiphaseInterFoam: Added test for contact angle in both phases in interface pair 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 --- .../multiphaseMixture/multiphaseMixture.C | 60 ++++++++++++------- .../multiphaseMixture/multiphaseMixture.H | 5 ++ 2 files changed, 43 insertions(+), 22 deletions(-) diff --git a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C index e003518593..c8927bfcf1 100644 --- a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C +++ b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C @@ -391,12 +391,6 @@ Foam::tmp 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(gbf[patchi])) + if + ( + isA(a1bf[patchi]) + || isA(a2bf[patchi]) + ) { + if + ( + isA(a1bf[patchi]) + && isA(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(gbf[patchi]); + isA(a1bf[patchi]) + ? refCast(a1bf[patchi]) + : refCast(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; diff --git a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.H b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.H index b28040a8c9..e85f7c38a1 100644 --- a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.H +++ b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.H @@ -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,