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-2008 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 Application
  26     interfaceProperties
  27
  28 Description
  29     Properties to aid interFoam :
  30     1. Correct the gamma boundary condition for dynamic contact angle.
  31     2. Calculate interface curvature.
  32
  33 \*---------------------------------------------------------------------------*/
  34
  35 #include "interfaceProperties.H"
  36 #include "gammaContactAngleFvPatchScalarField.H"
  37 #include "mathematicalConstants.H"
  38 #include "surfaceInterpolate.H"
  39 #include "fvcDiv.H"
  40 #include "fvcGrad.H"
  41 #include "fvcSnGrad.H"
  42
  43 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  44
  45 namespace Foam
  46 {
  47
  48 // * * * * * * * * * * * * * * * Static Member Data  * * * * * * * * * * * * //
  49
  50 const scalar interfaceProperties::convertToRad =
  51     mathematicalConstant::pi/180.0;
  52
  53
  54 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  55
  56 // Correction for the boundary condition on the unit normal nHat on
  57 // walls to produce the correct contact angle.
  58
  59 // The dynamic contact angle is calculated from the component of the
  60 // velocity on the direction of the interface, parallel to the wall.
  61
  62 void interfaceProperties::correctContactAngle
  63 (
  64     surfaceVectorField::GeometricBoundaryField& nHatb
  65 ) const
  66 {
  67     const fvMesh& mesh = gamma_.mesh();
  68     const volScalarField::GeometricBoundaryField& gbf = gamma_.boundaryField();
  69
  70     const fvBoundaryMesh& boundary = mesh.boundary();
  71
  72     forAll(boundary, patchi)
  73     {
  74         if (isA<gammaContactAngleFvPatchScalarField>(gbf[patchi]))
  75         {
  76             const gammaContactAngleFvPatchScalarField& gcap =
  77                 refCast<const gammaContactAngleFvPatchScalarField>
  78                 (gbf[patchi]);
  79
  80             fvsPatchVectorField& nHatp = nHatb[patchi];
  81             scalarField theta =
  82                 convertToRad*gcap.theta(U_.boundaryField()[patchi], nHatp);
  83
  84             vectorField nf = boundary[patchi].nf();
  85
  86             // Reset nHatp to correspond to the contact angle
  87
  88             scalarField a12 = nHatp & nf;
  89
  90             scalarField b1 = cos(theta);
  91
  92             scalarField b2(nHatp.size());
  93
  94             forAll(b2, facei)
  95             {
  96                 b2[facei] = cos(acos(a12[facei]) - theta[facei]);
  97             }
  98
  99             scalarField det = 1.0 - a12*a12;
 100
 101             scalarField a = (b1 - a12*b2)/det;
 102             scalarField b = (b2 - a12*b1)/det;
 103
 104             nHatp = a*nf + b*nHatp;
 105
 106             nHatp /= (mag(nHatp) + deltaN_.value());
 107         }
 108     }
 109 }
 110
 111
 112 void interfaceProperties::calculateK()
 113 {
 114     const fvMesh& mesh = gamma_.mesh();
 115     const surfaceVectorField& Sf = mesh.Sf();
 116
 117     // Cell gradient of gamma
 118     volVectorField gradGamma = fvc::grad(gamma_);
 119
 120     // Interpolated face-gradient of gamma
 121     surfaceVectorField gradGammaf = fvc::interpolate(gradGamma);
 122     //gradGammaf -=
 123     //    (mesh.Sf()/mesh.magSf())
 124     //   *(fvc::snGrad(gamma_) - (mesh.Sf() & gradGammaf)/mesh.magSf());
 125
 126     // Face unit interface normal
 127     surfaceVectorField nHatfv = gradGammaf/(mag(gradGammaf) + deltaN_);
 128     correctContactAngle(nHatfv.boundaryField());
 129
 130     // Face unit interface normal flux
 131     nHatf_ = nHatfv & Sf;
 132
 133     // Simple expression for curvature
 134     K_ = -fvc::div(nHatf_);
 135
 136     // Complex expression for curvature.
 137     // Correction is formally zero but numerically non-zero.
 138     /*
 139     volVectorField nHat = gradGamma/(mag(gradGamma) + deltaN_);
 140     forAll(nHat.boundaryField(), patchi)
 141     {
 142         nHat.boundaryField()[patchi] = nHatfv.boundaryField()[patchi];
 143     }
 144
 145     K_ = -fvc::div(nHatf_) + (nHat & fvc::grad(nHatfv) & nHat);
 146     */
 147 }
 148
 149
 150 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 151
 152 interfaceProperties::interfaceProperties
 153 (
 154     const volScalarField& gamma,
 155     const volVectorField& U,
 156     const IOdictionary& dict
 157 )
 158 :
 159     transportPropertiesDict_(dict),
 160     cGamma_
 161     (
 162         readScalar
 163         (
 164             gamma.mesh().solutionDict().subDict("PISO").lookup("cGamma")
 165         )
 166     ),
 167     sigma_(dict.lookup("sigma")),
 168
 169     deltaN_
 170     (
 171         "deltaN",
 172         1e-8/pow(average(gamma.mesh().V()), 1.0/3.0)
 173     ),
 174
 175     gamma_(gamma),
 176     U_(U),
 177
 178     nHatf_
 179     (
 180         IOobject
 181         (
 182             "nHatf",
 183             gamma_.time().timeName(),
 184             gamma_.mesh()
 185         ),
 186         gamma_.mesh(),
 187         dimensionedScalar("nHatf", dimArea, 0.0)
 188     ),
 189
 190     K_
 191     (
 192         IOobject
 193         (
 194             "K",
 195             gamma_.time().timeName(),
 196             gamma_.mesh()
 197         ),
 198         gamma_.mesh(),
 199         dimensionedScalar("K", dimless/dimLength, 0.0)
 200     )
 201 {
 202     calculateK();
 203 }
 204
 205
 206 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 207
 208 } // End namespace Foam
 209
 210 // ************************************************************************* //