src/lagrangian/spray/submodels/BreakupModel/PilchErdman/PilchErdman.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 "PilchErdman.H"
  27
  28 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  29
  30 template <class CloudType>
  31 Foam::PilchErdman<CloudType>::PilchErdman
  32 (
  33     const dictionary& dict,
  34     CloudType& owner
  35 )
  36 :
  37     BreakupModel<CloudType>(dict, owner, typeName),
  38     B1_(0.375),
  39     B2_(0.236)
  40 {
  41     if (!this->defaultCoeffs(true))
  42     {
  43         this->coeffDict().lookup("B1") >> B1_;
  44         this->coeffDict().lookup("B2") >> B2_;
  45     }
  46 }
  47
  48
  49 template <class CloudType>
  50 Foam::PilchErdman<CloudType>::PilchErdman(const PilchErdman<CloudType>& bum)
  51 :
  52     BreakupModel<CloudType>(bum),
  53     B1_(bum.B1_),
  54     B2_(bum.B2_)
  55 {}
  56
  57
  58 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
  59
  60 template <class CloudType>
  61 Foam::PilchErdman<CloudType>::~PilchErdman()
  62 {}
  63
  64
  65 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  66
  67 template<class CloudType>
  68 bool Foam::PilchErdman<CloudType>::update
  69 (
  70     const scalar dt,
  71     const vector& g,
  72     scalar& d,
  73     scalar& tc,
  74     scalar& ms,
  75     scalar& nParticle,
  76     scalar& KHindex,
  77     scalar& y,
  78     scalar& yDot,
  79     const scalar d0,
  80     const scalar rho,
  81     const scalar mu,
  82     const scalar sigma,
  83     const vector& U,
  84     const scalar rhoc,
  85     const scalar muc,
  86     const vector& Urel,
  87     const scalar Urmag,
  88     const scalar tMom,
  89     const scalar averageParcelMass,
  90     scalar& dChild,
  91     scalar& massChild,
  92     cachedRandom& rndGen
  93 ) const
  94 {
  95     scalar semiMass = nParticle*pow3(d);
  96     scalar We = 0.5*rhoc*sqr(Urmag)*d/sigma;
  97     scalar Oh = mu/pow(rho*d*sigma, 0.5);
  98
  99     scalar Wec = 6.0*(1.0 + 1.077*pow(Oh, 1.6));
 100
 101     if (We > Wec)
 102     {
 103         // We > 1335, wave crest stripping
 104         scalar taubBar = 5.5;
 105
 106         if (We > 175.0)
 107         {
 108             // sheet stripping
 109             taubBar = 0.766*pow(2.0*We - 12.0, 0.25);
 110         }
 111         else if (We > 22.0)
 112         {
 113             // Bag-and-stamen breakup
 114             taubBar = 14.1*pow(2.0*We - 12.0, -0.25);
 115         }
 116         else if (We > 9.0)
 117         {
 118             // Bag breakup
 119             taubBar = 2.45*pow(2.0*We - 12.0, 0.25);
 120         }
 121         else if (We > 6.0)
 122         {
 123             // Vibrational breakup
 124             taubBar = 6.0*pow(2.0*We - 12.0, -0.25);
 125         }
 126
 127         scalar rho12 = pow(rhoc/rho, 0.5);
 128
 129         scalar Vd = Urmag*rho12*(B1_*taubBar * B2_*taubBar*taubBar);
 130         scalar Vd1 = sqr(1.0 - Vd/Urmag);
 131         Vd1 = max(Vd1, SMALL);
 132         scalar Ds = 2.0*Wec*sigma*Vd1/(Vd1*rhoc*sqr(Urmag));
 133         scalar A = Urmag*rho12/d;
 134
 135         scalar taub = taubBar/A;
 136
 137         scalar frac = dt/taub;
 138
 139         // update the droplet diameter according to the rate eq. (implicitly)
 140         d = (d + frac*Ds)/(1.0 + frac);
 141
 142         // correct the number of particles to conserve mass
 143         nParticle = semiMass/pow3(d);
 144     }
 145
 146     return false;
 147 }
 148
 149
 150 // ************************************************************************* //