applications/solvers/multiphase/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation, either version 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Application
  25     threePhaseInterfaceProperties
  26
  27 Description
  28     Properties to aid interFoam :
  29     1. Correct the alpha boundary condition for dynamic contact angle.
  30     2. Calculate interface curvature.
  31
  32 \*---------------------------------------------------------------------------*/
  33
  34 #include "threePhaseInterfaceProperties.H"
  35 #include "alphaContactAngleFvPatchScalarField.H"
  36 #include "mathematicalConstants.H"
  37 #include "surfaceInterpolate.H"
  38 #include "fvcDiv.H"
  39 #include "fvcGrad.H"
  40
  41 // * * * * * * * * * * * * * * * Static Member Data  * * * * * * * * * * * * //
  42
  43 const Foam::scalar Foam::threePhaseInterfaceProperties::convertToRad =
  44     Foam::mathematicalConstant::pi/180.0;
  45
  46
  47 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  48
  49 // Correction for the boundary condition on the unit normal nHat on
  50 // walls to produce the correct contact angle.
  51
  52 // The dynamic contact angle is calculated from the component of the
  53 // velocity on the direction of the interface, parallel to the wall.
  54
  55 void Foam::threePhaseInterfaceProperties::correctContactAngle
  56 (
  57     surfaceVectorField::GeometricBoundaryField& nHatb
  58 ) const
  59 {
  60     const volScalarField::GeometricBoundaryField& alpha1 =
  61         mixture_.alpha1().boundaryField();
  62     const volScalarField::GeometricBoundaryField& alpha2 =
  63         mixture_.alpha2().boundaryField();
  64     const volScalarField::GeometricBoundaryField& alpha3 =
  65         mixture_.alpha3().boundaryField();
  66     const volVectorField::GeometricBoundaryField& U =
  67         mixture_.U().boundaryField();
  68
  69     const fvMesh& mesh = mixture_.U().mesh();
  70     const fvBoundaryMesh& boundary = mesh.boundary();
  71
  72     forAll(boundary, patchi)
  73     {
  74         if (isA<alphaContactAngleFvPatchScalarField>(alpha1[patchi]))
  75         {
  76             const alphaContactAngleFvPatchScalarField& a2cap =
  77                 refCast<const alphaContactAngleFvPatchScalarField>
  78                 (alpha2[patchi]);
  79
  80             const alphaContactAngleFvPatchScalarField& a3cap =
  81                 refCast<const alphaContactAngleFvPatchScalarField>
  82                 (alpha3[patchi]);
  83
  84             scalarField twoPhaseAlpha2 = max(a2cap, scalar(0));
  85             scalarField twoPhaseAlpha3 = max(a3cap, scalar(0));
  86
  87             scalarField sumTwoPhaseAlpha =
  88                 twoPhaseAlpha2 + twoPhaseAlpha3 + SMALL;
  89
  90             twoPhaseAlpha2 /= sumTwoPhaseAlpha;
  91             twoPhaseAlpha3 /= sumTwoPhaseAlpha;
  92
  93             fvsPatchVectorField& nHatp = nHatb[patchi];
  94
  95             scalarField theta =
  96                 convertToRad
  97                *(
  98                    twoPhaseAlpha2*(180 - a2cap.theta(U[patchi], nHatp))
  99                  + twoPhaseAlpha3*(180 - a3cap.theta(U[patchi], nHatp))
 100                );
 101
 102             vectorField nf = boundary[patchi].nf();
 103
 104             // Reset nHatPatch to correspond to the contact angle
 105
 106             scalarField a12 = nHatp & nf;
 107
 108             scalarField b1 = cos(theta);
 109
 110             scalarField b2(nHatp.size());
 111
 112             forAll(b2, facei)
 113             {
 114                 b2[facei] = cos(acos(a12[facei]) - theta[facei]);
 115             }
 116
 117             scalarField det = 1.0 - a12*a12;
 118
 119             scalarField a = (b1 - a12*b2)/det;
 120             scalarField b = (b2 - a12*b1)/det;
 121
 122             nHatp = a*nf + b*nHatp;
 123
 124             nHatp /= (mag(nHatp) + deltaN_.value());
 125         }
 126     }
 127 }
 128
 129
 130 void Foam::threePhaseInterfaceProperties::calculateK()
 131 {
 132     const volScalarField& alpha1 = mixture_.alpha1();
 133
 134     const fvMesh& mesh = alpha1.mesh();
 135     const surfaceVectorField& Sf = mesh.Sf();
 136
 137     // Cell gradient of alpha
 138     volVectorField gradAlpha = fvc::grad(alpha1);
 139
 140     // Interpolated face-gradient of alpha
 141     surfaceVectorField gradAlphaf = fvc::interpolate(gradAlpha);
 142
 143     // Face unit interface normal
 144     surfaceVectorField nHatfv = gradAlphaf/(mag(gradAlphaf) + deltaN_);
 145     correctContactAngle(nHatfv.boundaryField());
 146
 147     // Face unit interface normal flux
 148     nHatf_ = nHatfv & Sf;
 149
 150     // Simple expression for curvature
 151     K_ = -fvc::div(nHatf_);
 152
 153     // Complex expression for curvature.
 154     // Correction is formally zero but numerically non-zero.
 155     //volVectorField nHat = gradAlpha/(mag(gradAlpha) + deltaN_);
 156     //nHat.boundaryField() = nHatfv.boundaryField();
 157     //K_ = -fvc::div(nHatf_) + (nHat & fvc::grad(nHatfv) & nHat);
 158 }
 159
 160
 161 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 162
 163 Foam::threePhaseInterfaceProperties::threePhaseInterfaceProperties
 164 (
 165     const threePhaseMixture& mixture
 166 )
 167 :
 168     mixture_(mixture),
 169     cAlpha_
 170     (
 171         readScalar
 172         (
 173             mixture.U().mesh().solutionDict().subDict("PISO").lookup("cAlpha")
 174         )
 175     ),
 176     sigma12_(mixture.lookup("sigma12")),
 177     sigma13_(mixture.lookup("sigma13")),
 178
 179     deltaN_
 180     (
 181         "deltaN",
 182         1e-8/pow(average(mixture.U().mesh().V()), 1.0/3.0)
 183     ),
 184
 185     nHatf_
 186     (
 187         (
 188             fvc::interpolate(fvc::grad(mixture.alpha1()))
 189            /(mag(fvc::interpolate(fvc::grad(mixture.alpha1()))) + deltaN_)
 190         ) & mixture.alpha1().mesh().Sf()
 191     ),
 192
 193     K_
 194     (
 195         IOobject
 196         (
 197             "K",
 198             mixture.alpha1().time().timeName(),
 199             mixture.alpha1().mesh()
 200         ),
 201         -fvc::div(nHatf_)
 202     )
 203 {
 204     calculateK();
 205 }
 206
 207
 208 // ************************************************************************* //