interfaceProperties.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 
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::GeometricBoundaryField& nHatb,
51  surfaceVectorField::GeometricBoundaryField& gradAlphaf
52 ) const
53 {
54  const fvMesh& mesh = alpha1_.mesh();
55  const volScalarField::GeometricBoundaryField& 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.boundaryField(), 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 
155 Foam::interfaceProperties::interfaceProperties
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  sigma_("sigma", dimensionSet(1, 0, -2, 0, 0), dict),
171 
172  deltaN_
173  (
174  "deltaN",
175  1e-8/pow(average(alpha1.mesh().V()), 1.0/3.0)
176  ),
177 
178  alpha1_(alpha1),
179  U_(U),
180 
181  nHatf_
182  (
183  IOobject
184  (
185  "nHatf",
186  alpha1_.time().timeName(),
187  alpha1_.mesh()
188  ),
189  alpha1_.mesh(),
190  dimensionedScalar("nHatf", dimArea, 0.0)
191  ),
192 
193  K_
194  (
195  IOobject
196  (
197  "interfaceProperties:K",
198  alpha1_.time().timeName(),
199  alpha1_.mesh()
200  ),
201  alpha1_.mesh(),
203  )
204 {
205  calculateK();
206 }
207 
208 
209 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
210 
213 {
214  return fvc::interpolate(sigmaK())*fvc::snGrad(alpha1_);
215 }
216 
217 
220 {
221  return pos(alpha1_ - 0.01)*pos(0.99 - alpha1_);
222 }
223 
224 
226 {
227  alpha1_.mesh().solverDict(alpha1_.name()).lookup("cAlpha") >> cAlpha_;
228  transportPropertiesDict_.lookup("sigma") >> sigma_;
229 
230  return true;
231 }
232 
233 
234 // ************************************************************************* //
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)
static const scalar convertToRad
Conversion factor for degrees into radians.
tmp< volScalarField > nearInterface() const
Indicator of the proximity of the interface.
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
bool read()
Read transportProperties dictionary.
tmp< volScalarField > sigmaK() const
const dimensionSet dimArea(sqr(dimLength))
Definition: dimensionSets.H:57
Surface Interpolation.
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.
dictionary dict
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)
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)
volScalarField & alpha1
Definition: createFields.H:15
const word & name() const
Return name.
Definition: IOobject.H:260
tmp< surfaceScalarField > surfaceTensionForce() const
GeometricField< vector, fvPatchField, volMesh > volVectorField
Definition: volFieldsFwd.H:55
Calculate the divergence of the given field.
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:452
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
U
Definition: pEqn.H:82
const Type & value() const
Return const reference to value.
word timeName
Definition: getTimeIndex.H:3