src/transportModels/interfaceProperties/interfaceProperties.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   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 "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::constant::mathematical::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             const scalarField theta
  74             (
  75                 convertToRad*acap.theta(U_.boundaryField()[patchi], nHatp)
  76             );
  77
  78             const vectorField nf
  79             (
  80                 boundary[patchi].nf()
  81             );
  82
  83             // Reset nHatp to correspond to the contact angle
  84
  85             const scalarField a12(nHatp & nf);
  86             const scalarField b1(cos(theta));
  87
  88             scalarField b2(nHatp.size());
  89             forAll(b2, facei)
  90             {
  91                 b2[facei] = cos(acos(a12[facei]) - theta[facei]);
  92             }
  93
  94             const scalarField det(1.0 - a12*a12);
  95
  96             scalarField a((b1 - a12*b2)/det);
  97             scalarField b((b2 - a12*b1)/det);
  98
  99             nHatp = a*nf + b*nHatp;
 100             nHatp /= (mag(nHatp) + deltaN_.value());
 101
 102             acap.gradient() = (nf & nHatp)*mag(gradAlphaf[patchi]);
 103             acap.evaluate();
 104         }
 105     }
 106 }
 107
 108
 109 void Foam::interfaceProperties::calculateK()
 110 {
 111     const fvMesh& mesh = alpha1_.mesh();
 112     const surfaceVectorField& Sf = mesh.Sf();
 113
 114     // Cell gradient of alpha
 115     const volVectorField gradAlpha(fvc::grad(alpha1_));
 116
 117     // Interpolated face-gradient of alpha
 118     surfaceVectorField gradAlphaf(fvc::interpolate(gradAlpha));
 119
 120     //gradAlphaf -=
 121     //    (mesh.Sf()/mesh.magSf())
 122     //   *(fvc::snGrad(alpha1_) - (mesh.Sf() & gradAlphaf)/mesh.magSf());
 123
 124     // Face unit interface normal
 125     surfaceVectorField nHatfv(gradAlphaf/(mag(gradAlphaf) + deltaN_));
 126     correctContactAngle(nHatfv.boundaryField(), gradAlphaf.boundaryField());
 127
 128     // Face unit interface normal flux
 129     nHatf_ = nHatfv & Sf;
 130
 131     // Simple expression for curvature
 132     K_ = -fvc::div(nHatf_);
 133
 134     // Complex expression for curvature.
 135     // Correction is formally zero but numerically non-zero.
 136     /*
 137     volVectorField nHat(gradAlpha/(mag(gradAlpha) + deltaN_));
 138     forAll(nHat.boundaryField(), patchi)
 139     {
 140         nHat.boundaryField()[patchi] = nHatfv.boundaryField()[patchi];
 141     }
 142
 143     K_ = -fvc::div(nHatf_) + (nHat & fvc::grad(nHatfv) & nHat);
 144     */
 145 }
 146
 147
 148 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 149
 150 Foam::interfaceProperties::interfaceProperties
 151 (
 152     const volScalarField& alpha1,
 153     const volVectorField& U,
 154     const IOdictionary& dict
 155 )
 156 :
 157     transportPropertiesDict_(dict),
 158     cAlpha_
 159     (
 160         readScalar
 161         (
 162             alpha1.mesh().solutionDict().subDict("PIMPLE").lookup("cAlpha")
 163         )
 164     ),
 165     sigma_(dict.lookup("sigma")),
 166
 167     deltaN_
 168     (
 169         "deltaN",
 170         1e-8/pow(average(alpha1.mesh().V()), 1.0/3.0)
 171     ),
 172
 173     alpha1_(alpha1),
 174     U_(U),
 175
 176     nHatf_
 177     (
 178         IOobject
 179         (
 180             "nHatf",
 181             alpha1_.time().timeName(),
 182             alpha1_.mesh()
 183         ),
 184         alpha1_.mesh(),
 185         dimensionedScalar("nHatf", dimArea, 0.0)
 186     ),
 187
 188     K_
 189     (
 190         IOobject
 191         (
 192             "K",
 193             alpha1_.time().timeName(),
 194             alpha1_.mesh()
 195         ),
 196         alpha1_.mesh(),
 197         dimensionedScalar("K", dimless/dimLength, 0.0)
 198     )
 199 {
 200     calculateK();
 201 }
 202
 203
 204 // ************************************************************************* //