src/transportModels/interfaceProperties/interfaceProperties.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 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 the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     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, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 \*---------------------------------------------------------------------------*/
  26
  27 #include "interfaceProperties.H"
  28 #include "alphaContactAngleFvPatchScalarField.H"
  29 #include "mathematicalConstants.H"
  30 #include "surfaceInterpolate.H"
  31 #include "fvcDiv.H"
  32 #include "fvcGrad.H"
  33 #include "fvcSnGrad.H"
  34
  35 // * * * * * * * * * * * * * * * Static Member Data  * * * * * * * * * * * * //
  36
  37 const Foam::scalar Foam::interfaceProperties::convertToRad =
  38     Foam::mathematicalConstant::pi/180.0;
  39
  40
  41 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  42
  43 // Correction for the boundary condition on the unit normal nHat on
  44 // walls to produce the correct contact angle.
  45
  46 // The dynamic contact angle is calculated from the component of the
  47 // velocity on the direction of the interface, parallel to the wall.
  48
  49 void Foam::interfaceProperties::correctContactAngle
  50 (
  51     surfaceVectorField::GeometricBoundaryField& nHatb
  52 ) const
  53 {
  54     const fvMesh& mesh = alpha1_.mesh();
  55     const volScalarField::GeometricBoundaryField& gbf = alpha1_.boundaryField();
  56
  57     const fvBoundaryMesh& boundary = mesh.boundary();
  58
  59     forAll(boundary, patchi)
  60     {
  61         if (isA<alphaContactAngleFvPatchScalarField>(gbf[patchi]))
  62         {
  63             const alphaContactAngleFvPatchScalarField& gcap =
  64                 refCast<const alphaContactAngleFvPatchScalarField>
  65                 (gbf[patchi]);
  66
  67             fvsPatchVectorField& nHatp = nHatb[patchi];
  68             scalarField theta =
  69                 convertToRad*gcap.theta(U_.boundaryField()[patchi], nHatp);
  70
  71             vectorField nf = boundary[patchi].nf();
  72
  73             // Reset nHatp to correspond to the contact angle
  74
  75             scalarField a12 = nHatp & nf;
  76
  77             scalarField b1 = cos(theta);
  78
  79             scalarField b2(nHatp.size());
  80
  81             forAll(b2, facei)
  82             {
  83                 b2[facei] = cos(acos(a12[facei]) - theta[facei]);
  84             }
  85
  86             scalarField det = 1.0 - a12*a12;
  87
  88             scalarField a = (b1 - a12*b2)/det;
  89             scalarField b = (b2 - a12*b1)/det;
  90
  91             nHatp = a*nf + b*nHatp;
  92
  93             nHatp /= (mag(nHatp) + deltaN_.value());
  94         }
  95     }
  96 }
  97
  98
  99 void Foam::interfaceProperties::calculateK()
 100 {
 101     const fvMesh& mesh = alpha1_.mesh();
 102     const surfaceVectorField& Sf = mesh.Sf();
 103
 104     // Cell gradient of alpha
 105     volVectorField gradAlpha = fvc::grad(alpha1_);
 106
 107     // Interpolated face-gradient of alpha
 108     surfaceVectorField gradAlphaf = fvc::interpolate(gradAlpha);
 109     //gradAlphaf -=
 110     //    (mesh.Sf()/mesh.magSf())
 111     //   *(fvc::snGrad(alpha1_) - (mesh.Sf() & gradAlphaf)/mesh.magSf());
 112
 113     // Face unit interface normal
 114     surfaceVectorField nHatfv = gradAlphaf/(mag(gradAlphaf) + deltaN_);
 115     correctContactAngle(nHatfv.boundaryField());
 116
 117     // Face unit interface normal flux
 118     nHatf_ = nHatfv & Sf;
 119
 120     // Simple expression for curvature
 121     K_ = -fvc::div(nHatf_);
 122
 123     // Complex expression for curvature.
 124     // Correction is formally zero but numerically non-zero.
 125     /*
 126     volVectorField nHat = gradAlpha/(mag(gradAlpha) + deltaN_);
 127     forAll(nHat.boundaryField(), patchi)
 128     {
 129         nHat.boundaryField()[patchi] = nHatfv.boundaryField()[patchi];
 130     }
 131
 132     K_ = -fvc::div(nHatf_) + (nHat & fvc::grad(nHatfv) & nHat);
 133     */
 134 }
 135
 136
 137 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 138
 139 Foam::interfaceProperties::interfaceProperties
 140 (
 141     const volScalarField& alpha1,
 142     const volVectorField& U,
 143     const IOdictionary& dict
 144 )
 145 :
 146     transportPropertiesDict_(dict),
 147     cAlpha_
 148     (
 149         readScalar
 150         (
 151             alpha1.mesh().solutionDict().subDict("PISO").lookup("cAlpha")
 152         )
 153     ),
 154     sigma_(dict.lookup("sigma")),
 155
 156     deltaN_
 157     (
 158         "deltaN",
 159         1e-8/pow(average(alpha1.mesh().V()), 1.0/3.0)
 160     ),
 161
 162     alpha1_(alpha1),
 163     U_(U),
 164
 165     nHatf_
 166     (
 167         IOobject
 168         (
 169             "nHatf",
 170             alpha1_.time().timeName(),
 171             alpha1_.mesh()
 172         ),
 173         alpha1_.mesh(),
 174         dimensionedScalar("nHatf", dimArea, 0.0)
 175     ),
 176
 177     K_
 178     (
 179         IOobject
 180         (
 181             "K",
 182             alpha1_.time().timeName(),
 183             alpha1_.mesh()
 184         ),
 185         alpha1_.mesh(),
 186         dimensionedScalar("K", dimless/dimLength, 0.0)
 187     )
 188 {
 189     calculateK();
 190 }
 191
 192
 193 // ************************************************************************* //