src/transportModels/interfaceProperties/interfaceProperties.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 \*---------------------------------------------------------------------------*/
  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 =
  37     Foam::mathematicalConstant::pi/180.0;
  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 ) const
  52 {
  53     const fvMesh& mesh = alpha1_.mesh();
  54     const volScalarField::GeometricBoundaryField& gbf = alpha1_.boundaryField();
  55
  56     const fvBoundaryMesh& boundary = mesh.boundary();
  57
  58     forAll(boundary, patchi)
  59     {
  60         if (isA<alphaContactAngleFvPatchScalarField>(gbf[patchi]))
  61         {
  62             const alphaContactAngleFvPatchScalarField& gcap =
  63                 refCast<const alphaContactAngleFvPatchScalarField>
  64                 (gbf[patchi]);
  65
  66             fvsPatchVectorField& nHatp = nHatb[patchi];
  67             scalarField theta =
  68                 convertToRad*gcap.theta(U_.boundaryField()[patchi], nHatp);
  69
  70             vectorField nf = boundary[patchi].nf();
  71
  72             // Reset nHatp to correspond to the contact angle
  73
  74             scalarField a12 = nHatp & nf;
  75
  76             scalarField b1 = cos(theta);
  77
  78             scalarField b2(nHatp.size());
  79
  80             forAll(b2, facei)
  81             {
  82                 b2[facei] = cos(acos(a12[facei]) - theta[facei]);
  83             }
  84
  85             scalarField det = 1.0 - a12*a12;
  86
  87             scalarField a = (b1 - a12*b2)/det;
  88             scalarField b = (b2 - a12*b1)/det;
  89
  90             nHatp = a*nf + b*nHatp;
  91
  92             nHatp /= (mag(nHatp) + deltaN_.value());
  93         }
  94     }
  95 }
  96
  97
  98 void Foam::interfaceProperties::calculateK()
  99 {
 100     const fvMesh& mesh = alpha1_.mesh();
 101     const surfaceVectorField& Sf = mesh.Sf();
 102
 103     // Cell gradient of alpha
 104     volVectorField gradAlpha = fvc::grad(alpha1_);
 105
 106     // Interpolated face-gradient of alpha
 107     surfaceVectorField gradAlphaf = fvc::interpolate(gradAlpha);
 108     //gradAlphaf -=
 109     //    (mesh.Sf()/mesh.magSf())
 110     //   *(fvc::snGrad(alpha1_) - (mesh.Sf() & gradAlphaf)/mesh.magSf());
 111
 112     // Face unit interface normal
 113     surfaceVectorField nHatfv = gradAlphaf/(mag(gradAlphaf) + deltaN_);
 114     correctContactAngle(nHatfv.boundaryField());
 115
 116     // Face unit interface normal flux
 117     nHatf_ = nHatfv & Sf;
 118
 119     // Simple expression for curvature
 120     K_ = -fvc::div(nHatf_);
 121
 122     // Complex expression for curvature.
 123     // Correction is formally zero but numerically non-zero.
 124     /*
 125     volVectorField nHat = gradAlpha/(mag(gradAlpha) + deltaN_);
 126     forAll(nHat.boundaryField(), patchi)
 127     {
 128         nHat.boundaryField()[patchi] = nHatfv.boundaryField()[patchi];
 129     }
 130
 131     K_ = -fvc::div(nHatf_) + (nHat & fvc::grad(nHatfv) & nHat);
 132     */
 133 }
 134
 135
 136 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 137
 138 Foam::interfaceProperties::interfaceProperties
 139 (
 140     const volScalarField& alpha1,
 141     const volVectorField& U,
 142     const IOdictionary& dict
 143 )
 144 :
 145     transportPropertiesDict_(dict),
 146     cAlpha_
 147     (
 148         readScalar
 149         (
 150             alpha1.mesh().solutionDict().subDict("PISO").lookup("cAlpha")
 151         )
 152     ),
 153     sigma_(dict.lookup("sigma")),
 154
 155     deltaN_
 156     (
 157         "deltaN",
 158         1e-8/pow(average(alpha1.mesh().V()), 1.0/3.0)
 159     ),
 160
 161     alpha1_(alpha1),
 162     U_(U),
 163
 164     nHatf_
 165     (
 166         IOobject
 167         (
 168             "nHatf",
 169             alpha1_.time().timeName(),
 170             alpha1_.mesh()
 171         ),
 172         alpha1_.mesh(),
 173         dimensionedScalar("nHatf", dimArea, 0.0)
 174     ),
 175
 176     K_
 177     (
 178         IOobject
 179         (
 180             "K",
 181             alpha1_.time().timeName(),
 182             alpha1_.mesh()
 183         ),
 184         alpha1_.mesh(),
 185         dimensionedScalar("K", dimless/dimLength, 0.0)
 186     )
 187 {
 188     calculateK();
 189 }
 190
 191
 192 // ************************************************************************* //