FreeFOAM The Cross-Platform CFD Toolkit
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) 1991-2011 OpenCFD Ltd.
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"
30 #include <finiteVolume/fvcDiv.H>
31 #include <finiteVolume/fvcGrad.H>
32 #include <finiteVolume/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  scalarField theta =
74  convertToRad*acap.theta(U_.boundaryField()[patchi], nHatp);
75 
76  vectorField nf = boundary[patchi].nf();
77 
78  // Reset nHatp to correspond to the contact angle
79 
80  scalarField a12 = nHatp & nf;
81 
82  scalarField b1 = cos(theta);
83 
84  scalarField b2(nHatp.size());
85 
86  forAll(b2, facei)
87  {
88  b2[facei] = cos(acos(a12[facei]) - theta[facei]);
89  }
90 
91  scalarField det = 1.0 - a12*a12;
92 
93  scalarField a = (b1 - a12*b2)/det;
94  scalarField b = (b2 - a12*b1)/det;
95 
96  nHatp = a*nf + b*nHatp;
97 
98  nHatp /= (mag(nHatp) + deltaN_.value());
99 
100  acap.gradient() = (nf & nHatp)*mag(gradAlphaf[patchi]);
101  acap.evaluate();
102  }
103  }
104 }
105 
106 
107 void Foam::interfaceProperties::calculateK()
108 {
109  const fvMesh& mesh = alpha1_.mesh();
110  const surfaceVectorField& Sf = mesh.Sf();
111 
112  // Cell gradient of alpha
113  volVectorField gradAlpha = fvc::grad(alpha1_);
114 
115  // Interpolated face-gradient of alpha
116  surfaceVectorField gradAlphaf = fvc::interpolate(gradAlpha);
117  //gradAlphaf -=
118  // (mesh.Sf()/mesh.magSf())
119  // *(fvc::snGrad(alpha1_) - (mesh.Sf() & gradAlphaf)/mesh.magSf());
120 
121  // Face unit interface normal
122  surfaceVectorField nHatfv = gradAlphaf/(mag(gradAlphaf) + deltaN_);
123  correctContactAngle(nHatfv.boundaryField(), 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 
147 Foam::interfaceProperties::interfaceProperties
148 (
149  const volScalarField& alpha1,
150  const volVectorField& U,
151  const IOdictionary& dict
152 )
153 :
154  transportPropertiesDict_(dict),
155  cAlpha_
156  (
157  readScalar
158  (
159  alpha1.mesh().solutionDict().subDict("PISO").lookup("cAlpha")
160  )
161  ),
162  sigma_(dict.lookup("sigma")),
163 
164  deltaN_
165  (
166  "deltaN",
167  1e-8/pow(average(alpha1.mesh().V()), 1.0/3.0)
168  ),
169 
170  alpha1_(alpha1),
171  U_(U),
172 
173  nHatf_
174  (
175  IOobject
176  (
177  "nHatf",
178  alpha1_.time().timeName(),
179  alpha1_.mesh()
180  ),
181  alpha1_.mesh(),
182  dimensionedScalar("nHatf", dimArea, 0.0)
183  ),
184 
185  K_
186  (
187  IOobject
188  (
189  "Kcond",
190  alpha1_.time().timeName(),
191  alpha1_.mesh()
192  ),
193  alpha1_.mesh(),
194  dimensionedScalar("Kcond", dimless/dimLength, 0.0)
195  )
196 {
197  calculateK();
198 }
199 
200 
201 // ************************ vim: set sw=4 sts=4 et: ************************ //