From 2ab81c3b6998836924806e0ddca1a51a417537c8 Mon Sep 17 00:00:00 2001
From: sergio <s.ferraris@opencfd.co.uk>
Date: Thu, 8 Jul 2021 12:15:41 -0700
Subject: [PATCH] BUG: Correcting surface tension calculation from pair order.
Fixes #2051
---
.../multiphaseMixture/multiphaseMixture.C | 148 ++++++++++--------
.../multiphaseMixture/multiphaseMixture.H | 12 ++
2 files changed, 97 insertions(+), 63 deletions(-)
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
index d3f850db07b..da3547d18a6 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
+++ b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
@@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2017 OpenFOAM Foundation
+ Copyright (C) 2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -26,7 +27,6 @@ License
\*---------------------------------------------------------------------------*/
#include "multiphaseMixture.H"
-#include "alphaContactAngleFvPatchScalarField.H"
#include "Time.H"
#include "subCycle.H"
#include "MULES.H"
@@ -36,6 +36,7 @@ License
#include "fvcDiv.H"
#include "fvcFlux.H"
#include "unitConversion.H"
+#include "alphaContactAngleFvPatchScalarField.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@@ -416,99 +417,120 @@ void Foam::multiphaseMixture::correctContactAngle
surfaceVectorField::Boundary& nHatb
) const
{
- const volScalarField::Boundary& gbf
- = alpha1.boundaryField();
+ const volScalarField::Boundary& gb1f = alpha1.boundaryField();
+ const volScalarField::Boundary& gb2f = alpha2.boundaryField();
const fvBoundaryMesh& boundary = mesh_.boundary();
forAll(boundary, patchi)
{
- if (isA<alphaContactAngleFvPatchScalarField>(gbf[patchi]))
+ if (isA<alphaContactAngleFvPatchScalarField>(gb1f[patchi]))
{
const alphaContactAngleFvPatchScalarField& acap =
- refCast<const alphaContactAngleFvPatchScalarField>(gbf[patchi]);
+ refCast<const alphaContactAngleFvPatchScalarField>(gb1f[patchi]);
- vectorField& nHatPatch = nHatb[patchi];
+ correctBoundaryContactAngle(acap, patchi, alpha1, alpha2, nHatb);
+ }
+ else if (isA<alphaContactAngleFvPatchScalarField>(gb2f[patchi]))
+ {
+ const alphaContactAngleFvPatchScalarField& acap =
+ refCast<const alphaContactAngleFvPatchScalarField>(gb2f[patchi]);
- vectorField AfHatPatch
- (
- mesh_.Sf().boundaryField()[patchi]
- /mesh_.magSf().boundaryField()[patchi]
- );
+ correctBoundaryContactAngle(acap, patchi, alpha2, alpha1, nHatb);
+ }
+ }
+}
- const auto tp =
- acap.thetaProps().cfind(interfacePair(alpha1, alpha2));
- if (!tp.found())
- {
- FatalErrorInFunction
- << "Cannot find interface " << interfacePair(alpha1, alpha2)
- << "\n in table of theta properties for patch "
- << acap.patch().name()
- << exit(FatalError);
- }
+void Foam::multiphaseMixture::correctBoundaryContactAngle
+(
+ const alphaContactAngleFvPatchScalarField& acap,
+ label patchi,
+ const phase& alpha1,
+ const phase& alpha2,
+ surfaceVectorField::Boundary& nHatb
+) const
+{
+ const fvBoundaryMesh& boundary = mesh_.boundary();
- const bool matched = (tp.key().first() == alpha1.name());
+ vectorField& nHatPatch = nHatb[patchi];
- const scalar theta0 = degToRad(tp().theta0(matched));
- scalarField theta(boundary[patchi].size(), theta0);
+ vectorField AfHatPatch
+ (
+ mesh_.Sf().boundaryField()[patchi]
+ /mesh_.magSf().boundaryField()[patchi]
+ );
- const scalar uTheta = tp().uTheta();
+ const auto tp = acap.thetaProps().cfind(interfacePair(alpha1, alpha2));
- // Calculate the dynamic contact angle if required
- if (uTheta > SMALL)
- {
- const scalar thetaA = degToRad(tp().thetaA(matched));
- const scalar thetaR = degToRad(tp().thetaR(matched));
+ if (!tp.found())
+ {
+ FatalErrorInFunction
+ << "Cannot find interface " << interfacePair(alpha1, alpha2)
+ << "\n in table of theta properties for patch "
+ << acap.patch().name()
+ << exit(FatalError);
+ }
- // Calculated the component of the velocity parallel to the wall
- vectorField Uwall
- (
- U_.boundaryField()[patchi].patchInternalField()
- - U_.boundaryField()[patchi]
- );
- Uwall -= (AfHatPatch & Uwall)*AfHatPatch;
+ const bool matched = (tp.key().first() == alpha1.name());
- // Find the direction of the interface parallel to the wall
- vectorField nWall
- (
- nHatPatch - (AfHatPatch & nHatPatch)*AfHatPatch
- );
+ const scalar theta0 = degToRad(tp().theta0(matched));
+ scalarField theta(boundary[patchi].size(), theta0);
- // Normalise nWall
- nWall /= (mag(nWall) + SMALL);
+ const scalar uTheta = tp().uTheta();
- // Calculate Uwall resolved normal to the interface parallel to
- // the interface
- scalarField uwall(nWall & Uwall);
+ // Calculate the dynamic contact angle if required
+ if (uTheta > SMALL)
+ {
+ const scalar thetaA = degToRad(tp().thetaA(matched));
+ const scalar thetaR = degToRad(tp().thetaR(matched));
- theta += (thetaA - thetaR)*tanh(uwall/uTheta);
- }
+ // Calculated the component of the velocity parallel to the wall
+ vectorField Uwall
+ (
+ U_.boundaryField()[patchi].patchInternalField()
+ - U_.boundaryField()[patchi]
+ );
+ Uwall -= (AfHatPatch & Uwall)*AfHatPatch;
+ // Find the direction of the interface parallel to the wall
+ vectorField nWall
+ (
+ nHatPatch - (AfHatPatch & nHatPatch)*AfHatPatch
+ );
- // Reset nHatPatch to correspond to the contact angle
+ // Normalise nWall
+ nWall /= (mag(nWall) + SMALL);
- scalarField a12(nHatPatch & AfHatPatch);
+ // Calculate Uwall resolved normal to the interface parallel to
+ // the interface
+ scalarField uwall(nWall & Uwall);
- scalarField b1(cos(theta));
+ theta += (thetaA - thetaR)*tanh(uwall/uTheta);
+ }
- scalarField b2(nHatPatch.size());
- forAll(b2, facei)
- {
- b2[facei] = cos(acos(a12[facei]) - theta[facei]);
- }
+ // Reset nHatPatch to correspond to the contact angle
- scalarField det(1.0 - a12*a12);
+ scalarField a12(nHatPatch & AfHatPatch);
- scalarField a((b1 - a12*b2)/det);
- scalarField b((b2 - a12*b1)/det);
+ scalarField b1(cos(theta));
- nHatPatch = a*AfHatPatch + b*nHatPatch;
+ scalarField b2(nHatPatch.size());
- nHatPatch /= (mag(nHatPatch) + deltaN_.value());
- }
+ forAll(b2, facei)
+ {
+ b2[facei] = cos(acos(a12[facei]) - theta[facei]);
}
+
+ scalarField det(1.0 - a12*a12);
+
+ scalarField a((b1 - a12*b2)/det);
+ scalarField b((b2 - a12*b1)/det);
+
+ nHatPatch = a*AfHatPatch + b*nHatPatch;
+
+ nHatPatch /= (mag(nHatPatch) + deltaN_.value());
}
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.H b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.H
index ddaa807a931..99bd1c38a96 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.H
+++ b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.H
@@ -57,6 +57,9 @@ SourceFiles
namespace Foam
{
+// Class forward declarations
+class alphaContactAngleFvPatchScalarField;
+
/*---------------------------------------------------------------------------*\
Class multiphaseMixture Declaration
\*---------------------------------------------------------------------------*/
@@ -170,6 +173,15 @@ private:
surfaceVectorField::Boundary& nHatb
) const;
+ void correctBoundaryContactAngle
+ (
+ const alphaContactAngleFvPatchScalarField& acap,
+ label patchi,
+ const phase& alpha1,
+ const phase& alpha2,
+ surfaceVectorField::Boundary& nHatb
+ ) const;
+
tmp<volScalarField> K(const phase& alpha1, const phase& alpha2) const;
--
GitLab