applications/solvers/multiphase/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.C

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