threePhaseInterfaceProperties.C
Go to the documentation of this file.
1/*---------------------------------------------------------------------------*\
2 ========= |
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4 \\ / O peration |
5 \\ / A nd | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2011-2017 OpenFOAM Foundation
9-------------------------------------------------------------------------------
10License
11 This file is part of OpenFOAM.
12
13 OpenFOAM is free software: you can redistribute it and/or modify it
14 under the terms of the GNU General Public License as published by
15 the Free Software Foundation, either version 3 of the License, or
16 (at your option) any later version.
17
18 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21 for more details.
22
23 You should have received a copy of the GNU General Public License
24 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25
26\*---------------------------------------------------------------------------*/
27
31#include "surfaceInterpolate.H"
32#include "fvcDiv.H"
33#include "fvcGrad.H"
34#include "fvcSnGrad.H"
35#include "unitConversion.H"
36
37// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
38
39void Foam::threePhaseInterfaceProperties::correctContactAngle
40(
41 surfaceVectorField::Boundary& nHatb
42) const
43{
45 mixture_.alpha1().boundaryField();
47 mixture_.alpha2().boundaryField();
48 const volScalarField::Boundary& alpha3 =
49 mixture_.alpha3().boundaryField();
51 mixture_.U().boundaryField();
52
53 const fvMesh& mesh = mixture_.U().mesh();
54 const fvBoundaryMesh& boundary = mesh.boundary();
55
56 forAll(boundary, patchi)
57 {
58 if (isA<alphaContactAngleTwoPhaseFvPatchScalarField>(alpha1[patchi]))
59 {
60 const alphaContactAngleTwoPhaseFvPatchScalarField& a2cap =
61 refCast<const alphaContactAngleTwoPhaseFvPatchScalarField>
62 (alpha2[patchi]);
63
64 const alphaContactAngleTwoPhaseFvPatchScalarField& a3cap =
65 refCast<const alphaContactAngleTwoPhaseFvPatchScalarField>
66 (alpha3[patchi]);
67
68 scalarField twoPhaseAlpha2(max(a2cap, scalar(0)));
69 scalarField twoPhaseAlpha3(max(a3cap, scalar(0)));
70
71 scalarField sumTwoPhaseAlpha
72 (
73 twoPhaseAlpha2 + twoPhaseAlpha3 + SMALL
74 );
75
76 twoPhaseAlpha2 /= sumTwoPhaseAlpha;
77 twoPhaseAlpha3 /= sumTwoPhaseAlpha;
78
79 fvsPatchVectorField& nHatp = nHatb[patchi];
80
81 scalarField theta
82 (
83 degToRad()
84 * (
85 twoPhaseAlpha2*(180 - a2cap.theta(U[patchi], nHatp))
86 + twoPhaseAlpha3*(180 - a3cap.theta(U[patchi], nHatp))
87 )
88 );
89
90 vectorField nf(boundary[patchi].nf());
91
92 // Reset nHatPatch to correspond to the contact angle
93
94 scalarField a12(nHatp & nf);
95
96 scalarField b1(cos(theta));
97
98 scalarField b2(nHatp.size());
99
100 forAll(b2, facei)
101 {
102 b2[facei] = cos(acos(a12[facei]) - theta[facei]);
103 }
104
105 scalarField det(1.0 - a12*a12);
106
107 scalarField a((b1 - a12*b2)/det);
108 scalarField b((b2 - a12*b1)/det);
109
110 nHatp = a*nf + b*nHatp;
111
112 nHatp /= (mag(nHatp) + deltaN_.value());
113 }
114 }
115}
116
117
118void Foam::threePhaseInterfaceProperties::calculateK()
119{
120 const volScalarField& alpha1 = mixture_.alpha1();
121
122 const fvMesh& mesh = alpha1.mesh();
123 const surfaceVectorField& Sf = mesh.Sf();
124
125 // Cell gradient of alpha
126 volVectorField gradAlpha(fvc::grad(alpha1));
127
128 // Interpolated face-gradient of alpha
129 surfaceVectorField gradAlphaf(fvc::interpolate(gradAlpha));
130
131 // Face unit interface normal
132 surfaceVectorField nHatfv(gradAlphaf/(mag(gradAlphaf) + deltaN_));
133
134 correctContactAngle(nHatfv.boundaryFieldRef());
135
136 // Face unit interface normal flux
137 nHatf_ = nHatfv & Sf;
138
139 // Simple expression for curvature
140 K_ = -fvc::div(nHatf_);
141
142 // Complex expression for curvature.
143 // Correction is formally zero but numerically non-zero.
144 // volVectorField nHat = gradAlpha/(mag(gradAlpha) + deltaN_);
145 // nHat.boundaryFieldRef() = nHatfv.boundaryField();
146 // K_ = -fvc::div(nHatf_) + (nHat & fvc::grad(nHatfv) & nHat);
147}
148
149
150// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
151
153(
154 const incompressibleThreePhaseMixture& mixture
155)
156:
157 mixture_(mixture),
158 cAlpha_
159 (
160 mixture.U().mesh().solverDict
161 (
162 mixture_.alpha1().name()
163 ).get<scalar>("cAlpha")
164 ),
165 sigma12_("sigma12", dimMass/sqr(dimTime), mixture),
166 sigma13_("sigma13", dimMass/sqr(dimTime), mixture),
167
168 deltaN_
169 (
170 "deltaN",
171 1e-8/cbrt(average(mixture.U().mesh().V()))
172 ),
173
174 nHatf_
175 (
176 IOobject
177 (
178 "nHatf",
179 mixture.alpha1().time().timeName(),
181 ),
182 mixture.alpha1().mesh(),
184 ),
185
186 K_
187 (
188 IOobject
189 (
190 "interfaceProperties:K",
191 mixture.alpha1().time().timeName(),
193 ),
194 mixture.alpha1().mesh(),
196 )
197{
198 calculateK();
199}
200
201
202// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
203
206{
207 return fvc::interpolate(sigmaK())*fvc::snGrad(mixture_.alpha1());
208}
209
210
213{
214 return max
215 (
216 pos0(mixture_.alpha1() - 0.01)*pos0(0.99 - mixture_.alpha1()),
217 pos0(mixture_.alpha2() - 0.01)*pos0(0.99 - mixture_.alpha2())
218 );
219}
220
221
222// ************************************************************************* //
const volScalarField & alpha1
const volScalarField & alpha2
const Mesh & mesh() const
Return mesh.
GeometricBoundaryField< scalar, fvPatchField, volMesh > Boundary
Type of boundary fields.
const Boundary & boundaryField() const
Return const-reference to the boundary field.
const Type & value() const
Return const reference to value.
const fvBoundaryMesh & boundary() const
Return reference to boundary mesh.
Definition: fvMesh.C:712
const surfaceVectorField & Sf() const
Return cell face area vectors.
const volVectorField & U() const
Return the velocity.
Properties to aid interFoam : 1. Correct the alpha boundary condition for dynamic contact angle....
tmp< volScalarField > nearInterface() const
Indicator of the proximity of the interface.
tmp< surfaceScalarField > surfaceTensionForce() const
A class for managing temporary objects.
Definition: tmp.H:65
U
Definition: pEqn.H:72
faceListList boundary
dynamicFvMesh & mesh
Calculate the divergence of the given field.
Calculate the gradient of the given field.
Calculate the snGrad of the given volField.
word timeName
Definition: getTimeIndex.H:3
List< ReturnType > get(const UPtrList< T > &list, const AccessOp &aop)
static tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > interpolate(const GeometricField< Type, fvPatchField, volMesh > &tvf, const surfaceScalarField &faceFlux, Istream &schemeData)
Interpolate field onto faces using scheme given by Istream.
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh > > grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcGrad.C:54
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > snGrad(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvcSnGrad.C:47
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcDiv.C:49
GeometricField< vector, fvsPatchField, surfaceMesh > surfaceVectorField
dimensionedScalar det(const dimensionedSphericalTensor &dt)
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
GeometricField< vector, fvPatchField, volMesh > volVectorField
Definition: volFieldsFwd.H:83
const dimensionSet dimless
Dimensionless.
dimensionedScalar pos0(const dimensionedScalar &ds)
dimensionedSymmTensor sqr(const dimensionedVector &dv)
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:52
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
Definition: dimensionSets.H:53
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:82
const dimensionSet dimArea(sqr(dimLength))
Definition: dimensionSets.H:59
dimensioned< Type > average(const DimensionedField< Type, GeoMesh > &df)
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
constexpr scalar degToRad() noexcept
Multiplication factor for degrees to radians conversion.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
Field< vector > vectorField
Specialisation of Field<T> for vector.
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
dimensionedScalar cbrt(const dimensionedScalar &ds)
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
fvsPatchField< vector > fvsPatchVectorField
dimensionedScalar cos(const dimensionedScalar &ds)
dimensionedScalar acos(const dimensionedScalar &ds)
const dimensionSet dimMass(1, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:51
volScalarField & b
Definition: createFields.H:27
volScalarField & e
Definition: createFields.H:11
Info<< "Creating temperaturePhaseChangeTwoPhaseMixture\n"<< endl;autoPtr< temperaturePhaseChangeTwoPhaseMixture > mixture
Definition: createFields.H:39
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
Unit conversion functions.