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 | Copyright (C) 2011-2015 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
27 #include "alphaContactAngleFvPatchScalarField.H"
28 #include "mathematicalConstants.H"
29 #include "surfaceInterpolate.H"
30 #include "fvcDiv.H"
31 #include "fvcGrad.H"
32 #include "fvcSnGrad.H"
33 
34 // * * * * * * * * * * * * * * * Static Member Data * * * * * * * * * * * * //
35 
38 
39 
40 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
41 
42 void Foam::threePhaseInterfaceProperties::correctContactAngle
43 (
44  surfaceVectorField::GeometricBoundaryField& nHatb
45 ) const
46 {
47  const volScalarField::GeometricBoundaryField& alpha1 =
48  mixture_.alpha1().boundaryField();
49  const volScalarField::GeometricBoundaryField& alpha2 =
50  mixture_.alpha2().boundaryField();
51  const volScalarField::GeometricBoundaryField& alpha3 =
52  mixture_.alpha3().boundaryField();
53  const volVectorField::GeometricBoundaryField& U =
54  mixture_.U().boundaryField();
55 
56  const fvMesh& mesh = mixture_.U().mesh();
57  const fvBoundaryMesh& boundary = mesh.boundary();
58 
59  forAll(boundary, patchi)
60  {
61  if (isA<alphaContactAngleFvPatchScalarField>(alpha1[patchi]))
62  {
63  const alphaContactAngleFvPatchScalarField& a2cap =
64  refCast<const alphaContactAngleFvPatchScalarField>
65  (alpha2[patchi]);
66 
67  const alphaContactAngleFvPatchScalarField& a3cap =
68  refCast<const alphaContactAngleFvPatchScalarField>
69  (alpha3[patchi]);
70 
71  scalarField twoPhaseAlpha2(max(a2cap, scalar(0)));
72  scalarField twoPhaseAlpha3(max(a3cap, scalar(0)));
73 
74  scalarField sumTwoPhaseAlpha
75  (
76  twoPhaseAlpha2 + twoPhaseAlpha3 + SMALL
77  );
78 
79  twoPhaseAlpha2 /= sumTwoPhaseAlpha;
80  twoPhaseAlpha3 /= sumTwoPhaseAlpha;
81 
82  fvsPatchVectorField& nHatp = nHatb[patchi];
83 
84  scalarField theta
85  (
87  * (
88  twoPhaseAlpha2*(180 - a2cap.theta(U[patchi], nHatp))
89  + twoPhaseAlpha3*(180 - a3cap.theta(U[patchi], nHatp))
90  )
91  );
92 
93  vectorField nf(boundary[patchi].nf());
94 
95  // Reset nHatPatch to correspond to the contact angle
96 
97  scalarField a12(nHatp & nf);
98 
99  scalarField b1(cos(theta));
100 
101  scalarField b2(nHatp.size());
102 
103  forAll(b2, facei)
104  {
105  b2[facei] = cos(acos(a12[facei]) - theta[facei]);
106  }
107 
108  scalarField det(1.0 - a12*a12);
109 
110  scalarField a((b1 - a12*b2)/det);
111  scalarField b((b2 - a12*b1)/det);
112 
113  nHatp = a*nf + b*nHatp;
114 
115  nHatp /= (mag(nHatp) + deltaN_.value());
116  }
117  }
118 }
119 
120 
121 void Foam::threePhaseInterfaceProperties::calculateK()
122 {
123  const volScalarField& alpha1 = mixture_.alpha1();
124 
125  const fvMesh& mesh = alpha1.mesh();
126  const surfaceVectorField& Sf = mesh.Sf();
127 
128  // Cell gradient of alpha
129  volVectorField gradAlpha(fvc::grad(alpha1));
130 
131  // Interpolated face-gradient of alpha
132  surfaceVectorField gradAlphaf(fvc::interpolate(gradAlpha));
133 
134  // Face unit interface normal
135  surfaceVectorField nHatfv(gradAlphaf/(mag(gradAlphaf) + deltaN_));
136 
137  correctContactAngle(nHatfv.boundaryField());
138 
139  // Face unit interface normal flux
140  nHatf_ = nHatfv & Sf;
141 
142  // Simple expression for curvature
143  K_ = -fvc::div(nHatf_);
144 
145  // Complex expression for curvature.
146  // Correction is formally zero but numerically non-zero.
147  // volVectorField nHat = gradAlpha/(mag(gradAlpha) + deltaN_);
148  // nHat.boundaryField() = nHatfv.boundaryField();
149  // K_ = -fvc::div(nHatf_) + (nHat & fvc::grad(nHatfv) & nHat);
150 }
151 
152 
153 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
154 
155 Foam::threePhaseInterfaceProperties::threePhaseInterfaceProperties
156 (
157  const incompressibleThreePhaseMixture& mixture
158 )
159 :
160  mixture_(mixture),
161  cAlpha_
162  (
163  readScalar
164  (
165  mixture.U().mesh().solverDict
166  (
167  mixture_.alpha1().name()
168  ).lookup("cAlpha")
169  )
170  ),
171  sigma12_(mixture.lookup("sigma12")),
172  sigma13_(mixture.lookup("sigma13")),
173 
174  deltaN_
175  (
176  "deltaN",
177  1e-8/pow(average(mixture.U().mesh().V()), 1.0/3.0)
178  ),
179 
180  nHatf_
181  (
182  IOobject
183  (
184  "nHatf",
185  mixture.alpha1().time().timeName(),
186  mixture.alpha1().mesh()
187  ),
188  mixture.alpha1().mesh(),
189  dimensionedScalar("nHatf", dimArea, 0.0)
190  ),
191 
192  K_
193  (
194  IOobject
195  (
196  "interfaceProperties:K",
197  mixture.alpha1().time().timeName(),
198  mixture.alpha1().mesh()
199  ),
200  mixture.alpha1().mesh(),
202  )
203 {
204  calculateK();
205 }
206 
207 
208 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
209 
212 {
213  return fvc::interpolate(sigmaK())*fvc::snGrad(mixture_.alpha1());
214 }
215 
216 
219 {
220  return max
221  (
222  pos(mixture_.alpha1() - 0.01)*pos(0.99 - mixture_.alpha1()),
223  pos(mixture_.alpha2() - 0.01)*pos(0.99 - mixture_.alpha2())
224  );
225 }
226 
227 
228 // ************************************************************************* //
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > interpolate(const GeometricField< Type, fvPatchField, volMesh > &vf, const surfaceScalarField &faceFlux, Istream &schemeData)
dimensioned< Type > average(const DimensionedField< Type, GeoMesh > &df)
Calculate the snGrad of the given volField.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
Definition: createFields.H:28
dimensioned< scalar > mag(const dimensioned< Type > &)
GeometricBoundaryField & boundaryField()
Return reference to GeometricBoundaryField.
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
const dimensionSet dimArea(sqr(dimLength))
Definition: dimensionSets.H:57
Surface Interpolation.
Info<< "Reading field p_rgh\n"<< endl;volScalarField p_rgh(IOobject("p_rgh", runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE), mesh);Info<< "Reading field U\n"<< endl;volVectorField U(IOobject("U", runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE), mesh);Info<< "Creating phaseChangeTwoPhaseMixture\n"<< endl;autoPtr< phaseChangeTwoPhaseMixture > mixture
Definition: createFields.H:33
const Mesh & mesh() const
Return mesh.
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:50
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Calculate the gradient of the given field.
fvsPatchField< vector > fvsPatchVectorField
const double e
Elementary charge.
Definition: doubleFloat.H:78
stressControl lookup("compactNormalStress") >> compactNormalStress
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:52
dimensionedScalar det(const dimensionedSphericalTensor &dt)
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
tmp< surfaceScalarField > surfaceTensionForce() const
Field< vector > vectorField
Specialisation of Field<T> for vector.
#define forAll(list, i)
Definition: UList.H:421
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if succesful.
Definition: doubleScalar.H:63
dimensionedScalar cos(const dimensionedScalar &ds)
label patchi
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > snGrad(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvcSnGrad.C:45
dimensionedScalar acos(const dimensionedScalar &ds)
GeometricField< vector, fvPatchField, volMesh > volVectorField
Definition: volFieldsFwd.H:55
Calculate the divergence of the given field.
const volVectorField & U() const
Return the velocity.
tmp< volScalarField > nearInterface() const
Indicator of the proximity of the interface.
static const scalar convertToRad
Conversion factor for degrees into radians.
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:47
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcDiv.C:47
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcGrad.C:52
dimensionedScalar pos(const dimensionedScalar &ds)
GeometricField< vector, fvsPatchField, surfaceMesh > surfaceVectorField
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
A class for managing temporary objects.
Definition: PtrList.H:118
const Type & value() const
Return const reference to value.
word timeName
Definition: getTimeIndex.H:3