mirror of
https://develop.openfoam.com/Development/openfoam.git
synced 2025-11-28 03:28:01 +00:00
BUG: Correcting surface tension calculation from pair order. Fixes #2051
This commit is contained in:
sergio
@ -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,102 +417,123 @@ 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]
|
||||
);
|
||||
|
||||
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);
|
||||
}
|
||||
|
||||
const bool matched = (tp.key().first() == alpha1.name());
|
||||
|
||||
const scalar theta0 = degToRad(tp().theta0(matched));
|
||||
scalarField theta(boundary[patchi].size(), theta0);
|
||||
|
||||
const scalar uTheta = tp().uTheta();
|
||||
|
||||
// Calculate the dynamic contact angle if required
|
||||
if (uTheta > SMALL)
|
||||
{
|
||||
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
|
||||
(
|
||||
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
|
||||
);
|
||||
|
||||
// Normalise nWall
|
||||
nWall /= (mag(nWall) + SMALL);
|
||||
|
||||
// Calculate Uwall resolved normal to the interface parallel to
|
||||
// the interface
|
||||
scalarField uwall(nWall & Uwall);
|
||||
|
||||
theta += (thetaA - thetaR)*tanh(uwall/uTheta);
|
||||
}
|
||||
|
||||
|
||||
// Reset nHatPatch to correspond to the contact angle
|
||||
|
||||
scalarField a12(nHatPatch & AfHatPatch);
|
||||
|
||||
scalarField b1(cos(theta));
|
||||
|
||||
scalarField b2(nHatPatch.size());
|
||||
|
||||
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());
|
||||
correctBoundaryContactAngle(acap, patchi, alpha2, alpha1, nHatb);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void Foam::multiphaseMixture::correctBoundaryContactAngle
|
||||
(
|
||||
const alphaContactAngleFvPatchScalarField& acap,
|
||||
label patchi,
|
||||
const phase& alpha1,
|
||||
const phase& alpha2,
|
||||
surfaceVectorField::Boundary& nHatb
|
||||
) const
|
||||
{
|
||||
const fvBoundaryMesh& boundary = mesh_.boundary();
|
||||
|
||||
vectorField& nHatPatch = nHatb[patchi];
|
||||
|
||||
vectorField AfHatPatch
|
||||
(
|
||||
mesh_.Sf().boundaryField()[patchi]
|
||||
/mesh_.magSf().boundaryField()[patchi]
|
||||
);
|
||||
|
||||
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);
|
||||
}
|
||||
|
||||
const bool matched = (tp.key().first() == alpha1.name());
|
||||
|
||||
const scalar theta0 = degToRad(tp().theta0(matched));
|
||||
scalarField theta(boundary[patchi].size(), theta0);
|
||||
|
||||
const scalar uTheta = tp().uTheta();
|
||||
|
||||
// Calculate the dynamic contact angle if required
|
||||
if (uTheta > SMALL)
|
||||
{
|
||||
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
|
||||
(
|
||||
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
|
||||
);
|
||||
|
||||
// Normalise nWall
|
||||
nWall /= (mag(nWall) + SMALL);
|
||||
|
||||
// Calculate Uwall resolved normal to the interface parallel to
|
||||
// the interface
|
||||
scalarField uwall(nWall & Uwall);
|
||||
|
||||
theta += (thetaA - thetaR)*tanh(uwall/uTheta);
|
||||
}
|
||||
|
||||
|
||||
// Reset nHatPatch to correspond to the contact angle
|
||||
|
||||
scalarField a12(nHatPatch & AfHatPatch);
|
||||
|
||||
scalarField b1(cos(theta));
|
||||
|
||||
scalarField b2(nHatPatch.size());
|
||||
|
||||
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());
|
||||
}
|
||||
|
||||
|
||||
Foam::tmp<Foam::volScalarField> Foam::multiphaseMixture::K
|
||||
(
|
||||
const phase& alpha1,
|
||||
|
||||
@ -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;
|
||||
|
||||
|
||||
|
||||
Reference in New Issue
Block a user