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