src/lagrangian/spray/submodels/StochasticCollision/ORourkeCollision/ORourkeCollision.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2009-2011 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 "ORourkeCollision.H"
  27
  28 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  29
  30 template <class CloudType>
  31 Foam::ORourkeCollision<CloudType>::ORourkeCollision
  32 (
  33     const dictionary& dict,
  34     CloudType& owner
  35 )
  36 :
  37     StochasticCollisionModel<CloudType>(dict, owner, typeName),
  38     coalescence_(this->coeffDict().lookup("coalescence"))
  39 {}
  40
  41
  42 template <class CloudType>
  43 Foam::ORourkeCollision<CloudType>::ORourkeCollision
  44 (
  45     const ORourkeCollision<CloudType>& cm
  46 )
  47 :
  48     StochasticCollisionModel<CloudType>(cm),
  49     coalescence_(cm.coalescence_)
  50 {}
  51
  52
  53 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
  54
  55 template <class CloudType>
  56 Foam::ORourkeCollision<CloudType>::~ORourkeCollision()
  57 {}
  58
  59
  60 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  61
  62 template<class CloudType>
  63 bool Foam::ORourkeCollision<CloudType>::update
  64 (
  65     const scalar dt,
  66     cachedRandom& rndGen,
  67     vector& pos1,
  68     scalar& m1,
  69     scalar& d1,
  70     scalar& N1,
  71     vector& U1,
  72     scalar& rho1,
  73     scalar& T1,
  74     scalarField& Y1,
  75     const scalar sigma1,
  76     const label celli,
  77     const scalar voli,
  78     vector& pos2,
  79     scalar& m2,
  80     scalar& d2,
  81     scalar& N2,
  82     vector& U2,
  83     scalar& rho2,
  84     scalar& T2,
  85     scalarField& Y2,
  86     const scalar sigma2,
  87     const label cellj,
  88     const scalar volj
  89 ) const
  90 {
  91     // check if parcels belong to same cell
  92     if ((celli != cellj) || (m1 < VSMALL) || (m2 < VSMALL))
  93     {
  94         return false;
  95     }
  96
  97     bool coalescence = false;
  98
  99     scalar magVrel = mag(U1-U2);
 100     scalar sumD = d1 + d2;
 101     scalar nu0 = 0.25*constant::mathematical::pi*sumD*sumD*magVrel*dt/volj;
 102     scalar nMin = min(N1, N2);
 103     scalar nu = nMin*nu0;
 104     scalar collProb = exp(-nu);
 105     scalar xx = rndGen.sample01<scalar>();
 106
 107     // collision occurs
 108     if (xx > collProb)
 109     {
 110         if (d1 > d2)
 111         {
 112             coalescence = collideSorted
 113             (
 114                 dt,
 115                 rndGen,
 116                 pos1,
 117                 m1,
 118                 d1,
 119                 N1,
 120                 U1,
 121                 rho1,
 122                 T1,
 123                 Y1,
 124                 sigma1,
 125                 celli,
 126                 voli,
 127                 pos2,
 128                 m2,
 129                 d2,
 130                 N2,
 131                 U2,
 132                 rho2,
 133                 T2,
 134                 Y2,
 135                 sigma2,
 136                 cellj,
 137                 volj
 138             );
 139         }
 140         else
 141         {
 142             coalescence = collideSorted
 143             (
 144                 dt,
 145                 rndGen,
 146                 pos2,
 147                 m2,
 148                 d2,
 149                 N2,
 150                 U2,
 151                 rho2,
 152                 T2,
 153                 Y2,
 154                 sigma2,
 155                 cellj,
 156                 volj,
 157                 pos1,
 158                 m1,
 159                 d1,
 160                 N1,
 161                 U1,
 162                 rho1,
 163                 T1,
 164                 Y1,
 165                 sigma1,
 166                 celli,
 167                 voli
 168             );
 169         }
 170     }
 171     return coalescence;
 172 }
 173
 174
 175 template<class CloudType>
 176 bool Foam::ORourkeCollision<CloudType>::collideSorted
 177 (
 178     const scalar dt,
 179     cachedRandom& rndGen,
 180     vector& pos1,
 181     scalar& m1,
 182     scalar& d1,
 183     scalar& N1,
 184     vector& U1,
 185     scalar& rho1,
 186     scalar& T1,
 187     scalarField& Y1,
 188     const scalar sigma1,
 189     const label celli,
 190     const scalar voli,
 191     vector& pos2,
 192     scalar& m2,
 193     scalar& d2,
 194     scalar& N2,
 195     vector& U2,
 196     scalar& rho2,
 197     scalar& T2,
 198     scalarField& Y2,
 199     const scalar sigma2,
 200     const label cellj,
 201     const scalar volj
 202 ) const
 203 {
 204     bool coalescence = false;
 205
 206     vector vRel = U1 - U2;
 207     scalar magVRel = mag(vRel);
 208
 209     scalar mdMin = m2/N2;
 210
 211     scalar mTot = m1 + m2;
 212
 213     scalar gamma = d1/max(d2, 1.0e-12);
 214     scalar f = gamma*gamma*gamma + 2.7*gamma - 2.4*gamma*gamma;
 215
 216     vector momMax = m1*U1;
 217     vector momMin = m2*U2;
 218
 219     // use mass-averaged temperature to calculate We number
 220     scalar Tm = (T1*m1 + T2*m2)/mTot;
 221
 222     // interpolate the averaged surface tension
 223     scalar sigma = sigma1 + (sigma2 - sigma1)*(Tm - T1)/(T2 - T1);
 224
 225     sigma = max(1.0e-6, sigma);
 226     scalar Vtot = m1/rho1 + m2/rho2;
 227     scalar rho = mTot/Vtot;
 228
 229     scalar dMean = sqrt(d1*d2);
 230     scalar WeColl = max(1.0e-12, 0.5*rho*magVRel*magVRel*dMean/sigma);
 231
 232     scalar coalesceProb = min(1.0, 2.4*f/WeColl);
 233
 234     scalar prob = rndGen.sample01<scalar>();
 235
 236     // Coalescence
 237     if (prob < coalesceProb && coalescence_)
 238     {
 239         coalescence = true;
 240         // How 'many' of the droplets coalesce
 241         scalar nProb = prob*N2/N1;
 242
 243         // Conservation of mass, momentum and energy
 244         scalar m2Org = m2;
 245         scalar dm = N1*nProb*mdMin;
 246         m2 -= dm;
 247         scalar V2 = constant::mathematical::pi*pow3(d2)/6.0;
 248         N2 = m2/(rho2*V2);
 249
 250         scalar m1Org = m1;
 251         m1 += dm;
 252         T1 = (Tm*mTot - m2*T2)/m1;
 253
 254         U1 =(momMax + (1.0 - m2/m2Org)*momMin)/m1;
 255
 256         // update the liquid mass fractions
 257         Y1 = (m1Org*Y1 + dm*Y2)/m1;
 258     }
 259     // Grazing collision (no coalescence)
 260     else
 261     {
 262         scalar gf = sqrt(prob) - sqrt(coalesceProb);
 263         scalar denom = 1.0 - sqrt(coalesceProb);
 264         if (denom < 1.0e-5)
 265         {
 266             denom = 1.0;
 267         }
 268         gf /= denom;
 269
 270         // if gf negative, this means that coalescence is turned off
 271         // and these parcels should have coalesced
 272         gf = max(0.0, gf);
 273
 274         // gf -> 1 => v1p -> p1().U() ...
 275         // gf -> 0 => v1p -> momentum/(m1+m2)
 276
 277         vector mr = m1*U1 + m2*U2;
 278         vector v1p = (mr + m2*gf*vRel)/(m1+m2);
 279         vector v2p = (mr - m1*gf*vRel)/(m1+m2);
 280
 281         if (N1 < N2)
 282         {
 283             U1 = v1p;
 284             U2 = (N1*v2p + (N2-N1)*U2)/N2;
 285         }
 286         else
 287         {
 288             U1 = (N2*v1p + (N1-N2)*U1)/N1;
 289             U2 = v2p;
 290         }
 291
 292     }
 293
 294     return coalescence;
 295 }
 296
 297
 298 // ************************************************************************* //