FreeFOAM The Cross-Platform CFD Toolkit
threePhaseInterfaceProperties.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-2010 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  Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
24 
25 Application
26  threePhaseInterfaceProperties
27 
28 Description
29  Properties to aid interFoam :
30  1. Correct the alpha boundary condition for dynamic contact angle.
31  2. Calculate interface curvature.
32 
33 \*---------------------------------------------------------------------------*/
34 
39 #include <finiteVolume/fvcDiv.H>
40 #include <finiteVolume/fvcGrad.H>
41 
42 // * * * * * * * * * * * * * * * Static Member Data * * * * * * * * * * * * //
43 
46 
47 
48 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
49 
50 // Correction for the boundary condition on the unit normal nHat on
51 // walls to produce the correct contact angle.
52 
53 // The dynamic contact angle is calculated from the component of the
54 // velocity on the direction of the interface, parallel to the wall.
55 
56 void Foam::threePhaseInterfaceProperties::correctContactAngle
57 (
58  surfaceVectorField::GeometricBoundaryField& nHatb
59 ) const
60 {
61  const volScalarField::GeometricBoundaryField& alpha1 =
62  mixture_.alpha1().boundaryField();
63  const volScalarField::GeometricBoundaryField& alpha2 =
64  mixture_.alpha2().boundaryField();
65  const volScalarField::GeometricBoundaryField& alpha3 =
66  mixture_.alpha3().boundaryField();
67  const volVectorField::GeometricBoundaryField& U =
68  mixture_.U().boundaryField();
69 
70  const fvMesh& mesh = mixture_.U().mesh();
71  const fvBoundaryMesh& boundary = mesh.boundary();
72 
73  forAll(boundary, patchi)
74  {
75  if (isA<alphaContactAngleFvPatchScalarField>(alpha1[patchi]))
76  {
77  const alphaContactAngleFvPatchScalarField& a2cap =
78  refCast<const alphaContactAngleFvPatchScalarField>
79  (alpha2[patchi]);
80 
81  const alphaContactAngleFvPatchScalarField& a3cap =
82  refCast<const alphaContactAngleFvPatchScalarField>
83  (alpha3[patchi]);
84 
85  scalarField twoPhaseAlpha2 = max(a2cap, scalar(0));
86  scalarField twoPhaseAlpha3 = max(a3cap, scalar(0));
87 
88  scalarField sumTwoPhaseAlpha =
89  twoPhaseAlpha2 + twoPhaseAlpha3 + SMALL;
90 
91  twoPhaseAlpha2 /= sumTwoPhaseAlpha;
92  twoPhaseAlpha3 /= sumTwoPhaseAlpha;
93 
94  fvsPatchVectorField& nHatp = nHatb[patchi];
95 
96  scalarField theta =
98  *(
99  twoPhaseAlpha2*(180 - a2cap.theta(U[patchi], nHatp))
100  + twoPhaseAlpha3*(180 - a3cap.theta(U[patchi], nHatp))
101  );
102 
103  vectorField nf = boundary[patchi].nf();
104 
105  // Reset nHatPatch to correspond to the contact angle
106 
107  scalarField a12 = nHatp & nf;
108 
109  scalarField b1 = cos(theta);
110 
111  scalarField b2(nHatp.size());
112 
113  forAll(b2, facei)
114  {
115  b2[facei] = cos(acos(a12[facei]) - theta[facei]);
116  }
117 
118  scalarField det = 1.0 - a12*a12;
119 
120  scalarField a = (b1 - a12*b2)/det;
121  scalarField b = (b2 - a12*b1)/det;
122 
123  nHatp = a*nf + b*nHatp;
124 
125  nHatp /= (mag(nHatp) + deltaN_.value());
126  }
127  }
128 }
129 
130 
131 void Foam::threePhaseInterfaceProperties::calculateK()
132 {
133  const volScalarField& alpha1 = mixture_.alpha1();
134 
135  const fvMesh& mesh = alpha1.mesh();
136  const surfaceVectorField& Sf = mesh.Sf();
137 
138  // Cell gradient of alpha
139  volVectorField gradAlpha = fvc::grad(alpha1);
140 
141  // Interpolated face-gradient of alpha
142  surfaceVectorField gradAlphaf = fvc::interpolate(gradAlpha);
143 
144  // Face unit interface normal
145  surfaceVectorField nHatfv = gradAlphaf/(mag(gradAlphaf) + deltaN_);
146  correctContactAngle(nHatfv.boundaryField());
147 
148  // Face unit interface normal flux
149  nHatf_ = nHatfv & Sf;
150 
151  // Simple expression for curvature
152  K_ = -fvc::div(nHatf_);
153 
154  // Complex expression for curvature.
155  // Correction is formally zero but numerically non-zero.
156  //volVectorField nHat = gradAlpha/(mag(gradAlpha) + deltaN_);
157  //nHat.boundaryField() = nHatfv.boundaryField();
158  //K_ = -fvc::div(nHatf_) + (nHat & fvc::grad(nHatfv) & nHat);
159 }
160 
161 
162 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
163 
164 Foam::threePhaseInterfaceProperties::threePhaseInterfaceProperties
165 (
166  const threePhaseMixture& mixture
167 )
168 :
169  mixture_(mixture),
170  cAlpha_
171  (
172  readScalar
173  (
174  mixture.U().mesh().solutionDict().subDict("PISO").lookup("cAlpha")
175  )
176  ),
177  sigma12_(mixture.lookup("sigma12")),
178  sigma13_(mixture.lookup("sigma13")),
179 
180  deltaN_
181  (
182  "deltaN",
183  1e-8/pow(average(mixture.U().mesh().V()), 1.0/3.0)
184  ),
185 
186  nHatf_
187  (
188  (
189  fvc::interpolate(fvc::grad(mixture.alpha1()))
190  /(mag(fvc::interpolate(fvc::grad(mixture.alpha1()))) + deltaN_)
191  ) & mixture.alpha1().mesh().Sf()
192  ),
193 
194  K_
195  (
196  IOobject
197  (
198  "Kcond",
199  mixture.alpha1().time().timeName(),
200  mixture.alpha1().mesh()
201  ),
202  -fvc::div(nHatf_)
203  )
204 {
205  calculateK();
206 }
207 
208 
209 // ************************************************************************* //