All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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-2020 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 "unitConversion.H"
29 #include "surfaceInterpolate.H"
30 #include "fvcDiv.H"
31 #include "fvcGrad.H"
32 #include "fvcSnGrad.H"
33 
34 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
35 
36 // Correction for the boundary condition on the unit normal nHat on
37 // walls to produce the correct contact angle.
38 
39 // The dynamic contact angle is calculated from the component of the
40 // velocity on the direction of the interface, parallel to the wall.
41 
42 void Foam::interfaceProperties::correctContactAngle
43 (
44  surfaceVectorField::Boundary& nHatb,
45  const surfaceVectorField::Boundary& gradAlphaf
46 ) const
47 {
48  const fvMesh& mesh = alpha1_.mesh();
49  const volScalarField::Boundary& abf = alpha1_.boundaryField();
50 
51  const fvBoundaryMesh& boundary = mesh.boundary();
52 
53  forAll(boundary, patchi)
54  {
55  if (isA<alphaContactAngleFvPatchScalarField>(abf[patchi]))
56  {
57  alphaContactAngleFvPatchScalarField& acap =
58  const_cast<alphaContactAngleFvPatchScalarField&>
59  (
60  refCast<const alphaContactAngleFvPatchScalarField>
61  (
62  abf[patchi]
63  )
64  );
65 
66  fvsPatchVectorField& nHatp = nHatb[patchi];
67  const scalarField theta
68  (
69  degToRad(acap.theta(U_.boundaryField()[patchi], nHatp))
70  );
71 
72  const vectorField nf
73  (
74  boundary[patchi].nf()
75  );
76 
77  // Reset nHatp to correspond to the contact angle
78 
79  const scalarField a12(nHatp & nf);
80  const scalarField b1(cos(theta));
81 
82  scalarField b2(nHatp.size());
83  forAll(b2, facei)
84  {
85  b2[facei] = cos(acos(a12[facei]) - theta[facei]);
86  }
87 
88  const scalarField det(1.0 - a12*a12);
89 
90  scalarField a((b1 - a12*b2)/det);
91  scalarField b((b2 - a12*b1)/det);
92 
93  nHatp = a*nf + b*nHatp;
94  nHatp /= (mag(nHatp) + deltaN_.value());
95 
96  acap.gradient() = (nf & nHatp)*mag(gradAlphaf[patchi]);
97  acap.evaluate();
98  }
99  }
100 }
101 
102 
103 void Foam::interfaceProperties::calculateK()
104 {
105  const fvMesh& mesh = alpha1_.mesh();
106  const surfaceVectorField& Sf = mesh.Sf();
107 
108  // Cell gradient of alpha
109  const volVectorField gradAlpha(fvc::grad(alpha1_, "nHat"));
110 
111  // Interpolated face-gradient of alpha
112  surfaceVectorField gradAlphaf(fvc::interpolate(gradAlpha));
113 
114  // gradAlphaf -=
115  // (mesh.Sf()/mesh.magSf())
116  // *(fvc::snGrad(alpha1_) - (mesh.Sf() & gradAlphaf)/mesh.magSf());
117 
118  // Face unit interface normal
119  surfaceVectorField nHatfv(gradAlphaf/(mag(gradAlphaf) + deltaN_));
120  // surfaceVectorField nHatfv
121  // (
122  // (gradAlphaf + deltaN_*vector(0, 0, 1)
123  // *sign(gradAlphaf.component(vector::Z)))/(mag(gradAlphaf) + deltaN_)
124  // );
125  correctContactAngle(nHatfv.boundaryFieldRef(), gradAlphaf.boundaryField());
126 
127  // Face unit interface normal flux
128  nHatf_ = nHatfv & Sf;
129 
130  // Simple expression for curvature
131  K_ = -fvc::div(nHatf_);
132 
133  // Complex expression for curvature.
134  // Correction is formally zero but numerically non-zero.
135  /*
136  volVectorField nHat(gradAlpha/(mag(gradAlpha) + deltaN_));
137  forAll(nHat.boundaryField(), patchi)
138  {
139  nHat.boundaryField()[patchi] = nHatfv.boundaryField()[patchi];
140  }
141 
142  K_ = -fvc::div(nHatf_) + (nHat & fvc::grad(nHatfv) & nHat);
143  */
144 }
145 
146 
147 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
148 
150 (
151  const volScalarField& alpha1,
152  const volVectorField& U,
153  const IOdictionary& dict
154 )
155 :
156  transportPropertiesDict_(dict),
157 
158  sigmaPtr_(surfaceTensionModel::New(dict, alpha1.mesh())),
159 
160  deltaN_
161  (
162  "deltaN",
163  1e-8/pow(average(alpha1.mesh().V()), 1.0/3.0)
164  ),
165 
166  alpha1_(alpha1),
167  U_(U),
168 
169  nHatf_
170  (
171  IOobject
172  (
173  "nHatf",
174  alpha1_.time().timeName(),
175  alpha1_.mesh()
176  ),
177  alpha1_.mesh(),
179  ),
180 
181  K_
182  (
183  IOobject
184  (
185  "interfaceProperties:K",
186  alpha1_.time().timeName(),
187  alpha1_.mesh()
188  ),
189  alpha1_.mesh(),
191  )
192 {
193  calculateK();
194 }
195 
196 
197 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
198 
201 {
202  return sigmaPtr_->sigma()*K_;
203 }
204 
205 
208 {
209  return fvc::interpolate(sigmaK())*fvc::snGrad(alpha1_);
210 }
211 
212 
215 {
216  return pos0(alpha1_ - 0.01)*pos0(0.99 - alpha1_);
217 }
218 
219 
221 {
222  calculateK();
223 }
224 
225 
227 {
228  sigmaPtr_->readDict(transportPropertiesDict_);
229 
230  return true;
231 }
232 
233 
234 // ************************************************************************* //
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 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
Unit conversion functions.
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:622
dimensionedScalar det(const dimensionedSphericalTensor &dt)
GeometricField< vector, fvPatchField, volMesh > volVectorField
Definition: volFieldsFwd.H:58
tmp< volScalarField > nearInterface() const
Indicator of the proximity of the interface.
tmp< surfaceScalarField > surfaceTensionForce() const
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
Definition: IOdictionary.H:52
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)
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:328
dimensioned< scalar > mag(const dimensioned< Type > &)
Field< vector > vectorField
Specialisation of Field<T> for vector.
const doubleScalar e
Elementary charge.
Definition: doubleScalar.H:105
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