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-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
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "interfaceProperties.H"
  27 #include <interfaceProperties/alphaContactAngleFvPatchScalarField.H>
  28 #include <OpenFOAM/mathematicalConstants.H>
  29 #include <finiteVolume/surfaceInterpolate.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 =
  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     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         }
 102     }
 103 }
 104
 105
 106 void Foam::interfaceProperties::calculateK()
 107 {
 108     const fvMesh& mesh = alpha1_.mesh();
 109     const surfaceVectorField& Sf = mesh.Sf();
 110
 111     // Cell gradient of alpha
 112     volVectorField gradAlpha = fvc::grad(alpha1_);
 113
 114     // Interpolated face-gradient of alpha
 115     surfaceVectorField gradAlphaf = fvc::interpolate(gradAlpha);
 116     //gradAlphaf -=
 117     //    (mesh.Sf()/mesh.magSf())
 118     //   *(fvc::snGrad(alpha1_) - (mesh.Sf() & gradAlphaf)/mesh.magSf());
 119
 120     // Face unit interface normal
 121     surfaceVectorField nHatfv = gradAlphaf/(mag(gradAlphaf) + deltaN_);
 122     correctContactAngle(nHatfv.boundaryField(), gradAlphaf.boundaryField());
 123
 124     // Face unit interface normal flux
 125     nHatf_ = nHatfv & Sf;
 126
 127     // Simple expression for curvature
 128     K_ = -fvc::div(nHatf_);
 129
 130     // Complex expression for curvature.
 131     // Correction is formally zero but numerically non-zero.
 132     /*
 133     volVectorField nHat = gradAlpha/(mag(gradAlpha) + deltaN_);
 134     forAll(nHat.boundaryField(), patchi)
 135     {
 136         nHat.boundaryField()[patchi] = nHatfv.boundaryField()[patchi];
 137     }
 138
 139     K_ = -fvc::div(nHatf_) + (nHat & fvc::grad(nHatfv) & nHat);
 140     */
 141 }
 142
 143
 144 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 145
 146 Foam::interfaceProperties::interfaceProperties
 147 (
 148     const volScalarField& alpha1,
 149     const volVectorField& U,
 150     const IOdictionary& dict
 151 )
 152 :
 153     transportPropertiesDict_(dict),
 154     cAlpha_
 155     (
 156         readScalar
 157         (
 158             alpha1.mesh().solutionDict().subDict("PISO").lookup("cAlpha")
 159         )
 160     ),
 161     sigma_(dict.lookup("sigma")),
 162
 163     deltaN_
 164     (
 165         "deltaN",
 166         1e-8/pow(average(alpha1.mesh().V()), 1.0/3.0)
 167     ),
 168
 169     alpha1_(alpha1),
 170     U_(U),
 171
 172     nHatf_
 173     (
 174         IOobject
 175         (
 176             "nHatf",
 177             alpha1_.time().timeName(),
 178             alpha1_.mesh()
 179         ),
 180         alpha1_.mesh(),
 181         dimensionedScalar("nHatf", dimArea, 0.0)
 182     ),
 183
 184     K_
 185     (
 186         IOobject
 187         (
 188             "Kcond",
 189             alpha1_.time().timeName(),
 190             alpha1_.mesh()
 191         ),
 192         alpha1_.mesh(),
 193         dimensionedScalar("Kcond", dimless/dimLength, 0.0)
 194     )
 195 {
 196     calculateK();
 197 }
 198
 199
 200 // ************************ vim: set sw=4 sts=4 et: ************************ //