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-2018 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"
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 
36 const Foam::scalar Foam::interfaceProperties::convertToRad =
38 
39 
40 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
41 
42 // Correction for the boundary condition on the unit normal nHat on
43 // walls to produce the correct contact angle.
44 
45 // The dynamic contact angle is calculated from the component of the
46 // velocity on the direction of the interface, parallel to the wall.
47 
48 void Foam::interfaceProperties::correctContactAngle
49 (
50  surfaceVectorField::Boundary& nHatb,
51  const surfaceVectorField::Boundary& gradAlphaf
52 ) const
53 {
54  const fvMesh& mesh = alpha1_.mesh();
55  const volScalarField::Boundary& abf = alpha1_.boundaryField();
56 
57  const fvBoundaryMesh& boundary = mesh.boundary();
58 
59  forAll(boundary, patchi)
60  {
61  if (isA<alphaContactAngleFvPatchScalarField>(abf[patchi]))
62  {
63  alphaContactAngleFvPatchScalarField& acap =
64  const_cast<alphaContactAngleFvPatchScalarField&>
65  (
66  refCast<const alphaContactAngleFvPatchScalarField>
67  (
68  abf[patchi]
69  )
70  );
71 
72  fvsPatchVectorField& nHatp = nHatb[patchi];
73  const scalarField theta
74  (
75  convertToRad*acap.theta(U_.boundaryField()[patchi], nHatp)
76  );
77 
78  const vectorField nf
79  (
80  boundary[patchi].nf()
81  );
82 
83  // Reset nHatp to correspond to the contact angle
84 
85  const scalarField a12(nHatp & nf);
86  const scalarField b1(cos(theta));
87 
88  scalarField b2(nHatp.size());
89  forAll(b2, facei)
90  {
91  b2[facei] = cos(acos(a12[facei]) - theta[facei]);
92  }
93 
94  const scalarField det(1.0 - a12*a12);
95 
96  scalarField a((b1 - a12*b2)/det);
97  scalarField b((b2 - a12*b1)/det);
98 
99  nHatp = a*nf + b*nHatp;
100  nHatp /= (mag(nHatp) + deltaN_.value());
101 
102  acap.gradient() = (nf & nHatp)*mag(gradAlphaf[patchi]);
103  acap.evaluate();
104  }
105  }
106 }
107 
108 
109 void Foam::interfaceProperties::calculateK()
110 {
111  const fvMesh& mesh = alpha1_.mesh();
112  const surfaceVectorField& Sf = mesh.Sf();
113 
114  // Cell gradient of alpha
115  const volVectorField gradAlpha(fvc::grad(alpha1_, "nHat"));
116 
117  // Interpolated face-gradient of alpha
118  surfaceVectorField gradAlphaf(fvc::interpolate(gradAlpha));
119 
120  // gradAlphaf -=
121  // (mesh.Sf()/mesh.magSf())
122  // *(fvc::snGrad(alpha1_) - (mesh.Sf() & gradAlphaf)/mesh.magSf());
123 
124  // Face unit interface normal
125  surfaceVectorField nHatfv(gradAlphaf/(mag(gradAlphaf) + deltaN_));
126  // surfaceVectorField nHatfv
127  // (
128  // (gradAlphaf + deltaN_*vector(0, 0, 1)
129  // *sign(gradAlphaf.component(vector::Z)))/(mag(gradAlphaf) + deltaN_)
130  // );
131  correctContactAngle(nHatfv.boundaryFieldRef(), gradAlphaf.boundaryField());
132 
133  // Face unit interface normal flux
134  nHatf_ = nHatfv & Sf;
135 
136  // Simple expression for curvature
137  K_ = -fvc::div(nHatf_);
138 
139  // Complex expression for curvature.
140  // Correction is formally zero but numerically non-zero.
141  /*
142  volVectorField nHat(gradAlpha/(mag(gradAlpha) + deltaN_));
143  forAll(nHat.boundaryField(), patchi)
144  {
145  nHat.boundaryField()[patchi] = nHatfv.boundaryField()[patchi];
146  }
147 
148  K_ = -fvc::div(nHatf_) + (nHat & fvc::grad(nHatfv) & nHat);
149  */
150 }
151 
152 
153 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
154 
156 (
157  const volScalarField& alpha1,
158  const volVectorField& U,
159  const IOdictionary& dict
160 )
161 :
162  transportPropertiesDict_(dict),
163  cAlpha_
164  (
165  readScalar
166  (
167  alpha1.mesh().solverDict(alpha1.name()).lookup("cAlpha")
168  )
169  ),
170 
171  sigmaPtr_(surfaceTensionModel::New(dict, alpha1.mesh())),
172 
173  deltaN_
174  (
175  "deltaN",
176  1e-8/pow(average(alpha1.mesh().V()), 1.0/3.0)
177  ),
178 
179  alpha1_(alpha1),
180  U_(U),
181 
182  nHatf_
183  (
184  IOobject
185  (
186  "nHatf",
187  alpha1_.time().timeName(),
188  alpha1_.mesh()
189  ),
190  alpha1_.mesh(),
192  ),
193 
194  K_
195  (
196  IOobject
197  (
198  "interfaceProperties:K",
199  alpha1_.time().timeName(),
200  alpha1_.mesh()
201  ),
202  alpha1_.mesh(),
204  )
205 {
206  calculateK();
207 }
208 
209 
210 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
211 
214 {
215  return sigmaPtr_->sigma()*K_;
216 }
217 
218 
221 {
222  return fvc::interpolate(sigmaK())*fvc::snGrad(alpha1_);
223 }
224 
225 
228 {
229  return pos0(alpha1_ - 0.01)*pos0(0.99 - alpha1_);
230 }
231 
232 
234 {
235  calculateK();
236 }
237 
238 
240 {
241  alpha1_.mesh().solverDict(alpha1_.name()).lookup("cAlpha") >> cAlpha_;
242  sigmaPtr_->readDict(transportPropertiesDict_);
243 
244  return true;
245 }
246 
247 
248 // ************************************************************************* //
fvsPatchField< vector > fvsPatchVectorField
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcGrad.C:52
dimensionedScalar acos(const dimensionedScalar &ds)
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
const word & name() const
Return name.
Definition: IOobject.H:295
const Boundary & boundaryField() const
Return const-reference to the boundary field.
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcDiv.C:47
GeometricField< vector, fvsPatchField, surfaceMesh > surfaceVectorField
static autoPtr< surfaceTensionModel > New(const dictionary &dict, const fvMesh &mesh)
bool read()
Read transportProperties dictionary.
Calculate the snGrad of the given volField.
static word timeName(const scalar, const int precision=precision_)
Return time name of given scalar time.
Definition: Time.C:626
dimensionedScalar det(const dimensionedSphericalTensor &dt)
GeometricField< vector, fvPatchField, volMesh > volVectorField
Definition: volFieldsFwd.H:55
tmp< volScalarField > nearInterface() const
Indicator of the proximity of the interface.
tmp< surfaceScalarField > surfaceTensionForce() const
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
Definition: IOdictionary.H:52
stressControl lookup("compactNormalStress") >> compactNormalStress
dimensionedScalar cos(const dimensionedScalar &ds)
const dimensionedScalar b
Wien displacement law constant: default SI units: [m K].
Definition: createFields.H:27
Calculate the gradient of the given field.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
const Type & value() const
Return const reference to value.
dimensioned< Type > average(const DimensionedField< Type, GeoMesh > &df)
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if successful.
Definition: doubleScalar.H:68
static const scalar convertToRad
Conversion factor for degrees into radians.
Calculate the divergence of the given field.
dimensionedScalar pos0(const dimensionedScalar &ds)
const Mesh & mesh() const
Return mesh.
tmp< volScalarField > sigmaK() const
interfaceProperties(const volScalarField &alpha1, const volVectorField &U, const IOdictionary &)
Construct from volume fraction field gamma and IOdictionary.
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.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
label patchi
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:47
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:50
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
const Time & time() const
Return time.
Definition: IOobject.C:360
dimensioned< scalar > mag(const dimensioned< Type > &)
Field< vector > vectorField
Specialisation of Field<T> for vector.
const doubleScalar e
Elementary charge.
Definition: doubleScalar.H:98
mesh Sf()
A class for managing temporary objects.
Definition: PtrList.H:53
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > snGrad(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvcSnGrad.C:45
const dimensionSet dimArea(sqr(dimLength))
Definition: dimensionSets.H:57