interfaceProperties.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2024 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 
26 #include "interfaceProperties.H"
28 #include "surfaceInterpolate.H"
29 #include "fvcDiv.H"
30 #include "fvcGrad.H"
31 #include "fvcSnGrad.H"
32 
33 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
34 
35 // Correction for the boundary condition on the unit normal nHat on
36 // walls to produce the correct contact angle.
37 
38 // The dynamic contact angle is calculated from the component of the
39 // velocity on the direction of the interface, parallel to the wall.
40 
41 void Foam::interfaceProperties::correctContactAngle
42 (
43  surfaceVectorField::Boundary& nHatb,
44  const surfaceVectorField::Boundary& gradAlphaf
45 )
46 {
47  const fvMesh& mesh = alpha1_.mesh();
50 
51  const fvBoundaryMesh& boundary = mesh.boundary();
52 
54  {
55  if (isA<contactAngleFvPatchScalarField>(a1bf[patchi]))
56  {
57  contactAngleFvPatchScalarField& a1cap =
58  refCast<contactAngleFvPatchScalarField>
59  (
60  a1bf[patchi]
61  );
62 
63  fvsPatchVectorField& nHatp = nHatb[patchi];
64  const scalarField cosTheta
65  (
66  a1cap.cosTheta(U_.boundaryField()[patchi], nHatp)
67  );
68 
69  const vectorField nf
70  (
71  boundary[patchi].nf()
72  );
73 
74  // Reset nHatp to correspond to the contact angle
75 
76  const scalarField a12(nHatp & nf);
77  const scalarField b1(cosTheta);
78 
79  scalarField b2(nHatp.size());
80  forAll(b2, facei)
81  {
82  b2[facei] = cos(acos(a12[facei]) - acos(cosTheta[facei]));
83  }
84 
85  const scalarField det(1.0 - a12*a12);
86 
87  scalarField a((b1 - a12*b2)/det);
88  scalarField b((b2 - a12*b1)/det);
89 
90  nHatp = a*nf + b*nHatp;
91  nHatp /= (mag(nHatp) + deltaN_.value());
92 
93  a1cap.gradient() = (nf & nHatp)*mag(gradAlphaf[patchi]);
94  a1cap.evaluate();
95  a2bf[patchi] = 1 - a1cap;
96  }
97  }
98 }
99 
100 
101 void Foam::interfaceProperties::calculateK()
102 {
103  const fvMesh& mesh = alpha1_.mesh();
104  const surfaceVectorField& Sf = mesh.Sf();
105 
106  // Cell gradient of alpha
107  const volVectorField gradAlpha(fvc::grad(alpha1_, "nHat"));
108 
109  // Interpolated face-gradient of alpha
110  surfaceVectorField gradAlphaf(fvc::interpolate(gradAlpha));
111 
112  // gradAlphaf -=
113  // (mesh.Sf()/mesh.magSf())
114  // *(fvc::snGrad(alpha1_) - (mesh.Sf() & gradAlphaf)/mesh.magSf());
115 
116  // Face unit interface normal
117  surfaceVectorField nHatfv(gradAlphaf/(mag(gradAlphaf) + deltaN_));
118  // surfaceVectorField nHatfv
119  // (
120  // (gradAlphaf + deltaN_*vector(0, 0, 1)
121  // *sign(gradAlphaf.component(vector::Z)))/(mag(gradAlphaf) + deltaN_)
122  // );
123  correctContactAngle(nHatfv.boundaryFieldRef(), gradAlphaf.boundaryField());
124 
125  // Face unit interface normal flux
126  nHatf_ = nHatfv & Sf;
127 
128  // Simple expression for curvature
129  K_ = -fvc::div(nHatf_);
130 
131  // Complex expression for curvature.
132  // Correction is formally zero but numerically non-zero.
133  /*
134  volVectorField nHat(gradAlpha/(mag(gradAlpha) + deltaN_));
135  forAll(nHat.boundaryField(), patchi)
136  {
137  nHat.boundaryField()[patchi] = nHatfv.boundaryField()[patchi];
138  }
139 
140  K_ = -fvc::div(nHatf_) + (nHat & fvc::grad(nHatfv) & nHat);
141  */
142 }
143 
144 
145 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
146 
148 (
149  const IOdictionary& dict,
150  volScalarField& alpha1,
151  volScalarField& alpha2,
152  const volVectorField& U
153 )
154 :
155  phasePropertiesDict_(dict),
156 
157  alpha1_(alpha1),
158  alpha2_(alpha2),
159  U_(U),
160 
161  sigmaPtr_(surfaceTensionModel::New(dict, alpha1.mesh())),
162 
163  deltaN_
164  (
165  "deltaN",
166  1e-8/pow(average(alpha1.mesh().V()), 1.0/3.0)
167  ),
168 
169  nHatf_
170  (
171  IOobject
172  (
173  "nHatf",
174  alpha1_.time().name(),
175  alpha1_.mesh()
176  ),
177  alpha1_.mesh(),
179  ),
180 
181  K_
182  (
183  IOobject
184  (
185  "interfaceProperties:K",
186  alpha1_.time().name(),
187  alpha1_.mesh()
188  ),
189  alpha1_.mesh(),
191  )
192 {
193  calculateK();
194 }
195 
196 
197 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
198 
200 {
201  const volVectorField gradAlpha(fvc::grad(alpha1_));
202 
203  return volVectorField::New("n", gradAlpha/(mag(gradAlpha) + deltaN_));
204 }
205 
206 
209 {
210  return sigmaPtr_->sigma()*K_;
211 }
212 
213 
216 {
217  return fvc::interpolate(sigmaK())*fvc::snGrad(alpha1_);
218 }
219 
220 
223 {
224  return pos0(alpha1_ - 0.01)*pos0(0.99 - alpha1_);
225 }
226 
227 
229 {
230  calculateK();
231 }
232 
233 
235 {
236  sigmaPtr_->readDict(phasePropertiesDict_);
237 
238  return true;
239 }
240 
241 
242 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
const Mesh & mesh() const
Return mesh.
Generic GeometricField class.
GeometricBoundaryField< Type, PatchField, GeoMesh > Boundary
Type of the boundary field.
const Boundary & boundaryField() const
Return const-reference to the boundary field.
Boundary & boundaryFieldRef()
Return a reference to the boundary field.
static tmp< GeometricField< Type, PatchField, GeoMesh > > New(const word &name, const Internal &, const PtrList< PatchField< Type >> &, const HashPtrTable< Source > &=HashPtrTable< Source >())
Return a temporary field constructed from name,.
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
Definition: IOdictionary.H:57
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:99
const Type & value() const
Return const reference to value.
tmp< surfaceScalarField > surfaceTensionForce() const
tmp< volScalarField > nearInterface() const
Indicator of the proximity of the interface.
void correct()
Correct the curvature.
interfaceProperties(const IOdictionary &dict, volScalarField &alpha1, volScalarField &alpha2, const volVectorField &U)
Construct from dictionary, volume fraction fields and mixture.
tmp< volVectorField > n() const
tmp< volScalarField > sigmaK() const
bool read()
Read phaseProperties dictionary.
Abstract base-class for surface tension models which return the surface tension coefficient field.
A class for managing temporary objects.
Definition: tmp.H:55
Calculate the divergence of the given field.
Calculate the gradient of the given field.
Calculate the snGrad of the given volField.
label patchi
volScalarField & b
Definition: createFields.H:25
U
Definition: pEqn.H:72
static tmp< SurfaceField< Type > > interpolate(const VolField< Type > &tvf, const surfaceScalarField &faceFlux, Istream &schemeData)
Interpolate field onto faces using scheme given by Istream.
tmp< VolField< typename outerProduct< vector, Type >::type > > grad(const SurfaceField< Type > &ssf)
Definition: fvcGrad.C:46
tmp< VolField< Type > > div(const SurfaceField< Type > &ssf)
Definition: fvcDiv.C:47
tmp< SurfaceField< Type > > snGrad(const VolField< Type > &vf, const word &name)
Definition: fvcSnGrad.C:45
dimensionedScalar det(const dimensionedSphericalTensor &dt)
const doubleScalar e
Definition: doubleScalar.H:106
VolField< vector > volVectorField
Definition: volFieldsFwd.H:65
dimensionedScalar pos0(const dimensionedScalar &ds)
word name(const bool)
Return a word representation of a bool.
Definition: boolIO.C:39
const dimensionSet dimless
const dimensionSet dimLength
dimensioned< Type > average(const DimensionedField< Type, GeoMesh > &df)
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
void correctContactAngle(const volScalarField &alpha1, const volScalarField &alpha2, const volVectorField::Boundary &Ubf, const dimensionedScalar &deltaN, surfaceVectorField::Boundary &nHatbf)
Correct the contact angle for the two volume fraction fields.
dimensioned< scalar > mag(const dimensioned< Type > &)
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh >> &tdf1, const word &name, const dimensionSet &dimensions)
Field< vector > vectorField
Specialisation of Field<T> for vector.
const dimensionSet dimArea
SurfaceField< vector > surfaceVectorField
fvsPatchField< vector > fvsPatchVectorField
dimensionedScalar cos(const dimensionedScalar &ds)
dimensionedScalar acos(const dimensionedScalar &ds)
faceListList boundary(nPatches)
dictionary dict