src/lagrangian/dieselSpray/spraySubModels/collisionModel/ORourke/sameCell.H

   1 vector v1 = p1().U();
   2 vector v2 = p2().U();
   3
   4 vector vRel = v1 - v2;
   5 scalar magVRel = mag(vRel);
   6
   7 scalar sumD = p1().d() + p2().d();
   8 scalar pc = spray_.p()[p1().cell()];
   9
  10 spray::iterator pMin = p1;
  11 spray::iterator pMax = p2;
  12
  13 scalar dMin = pMin().d();
  14 scalar dMax = pMax().d();
  15
  16 if (dMin > dMax)
  17 {
  18     dMin = pMax().d();
  19     dMax = pMin().d();
  20     pMin = p2;
  21     pMax = p1;
  22 }
  23
  24 scalar rhoMax = spray_.fuels().rho(pc, pMax().T(), pMax().X());
  25 scalar rhoMin = spray_.fuels().rho(pc, pMin().T(), pMin().X());
  26 scalar mMax = pMax().m();
  27 scalar mMin = pMin().m();
  28 scalar mTot = mMax + mMin;
  29
  30 scalar nMax = pMax().N(rhoMax);
  31 scalar nMin = pMin().N(rhoMin);
  32
  33 scalar mdMin = mMin/nMin;
  34
  35 scalar nu0 = 0.25*constant::mathematical::pi*sqr(sumD)*magVRel*dt/vols_[cell1];
  36 scalar nu = nMin*nu0;
  37 scalar collProb = exp(-nu);
  38 scalar xx = rndGen_.sample01<scalar>();
  39
  40 if ((xx > collProb) && (mMin > VSMALL) && (mMax > VSMALL))
  41 {
  42     // collision occurs
  43
  44     scalar gamma = dMax/max(dMin, 1.0e-12);
  45     scalar f = gamma*gamma*gamma + 2.7*gamma - 2.4*gamma*gamma;
  46
  47     vector momMax = mMax*pMax().U();
  48     vector momMin = mMin*pMin().U();
  49
  50     // use mass-averaged temperature to calculate We number
  51     scalar averageTemp = (pMax().T()*mMax + pMin().T()*mMin)/mTot;
  52     // and mass averaged mole fractions ...
  53     scalarField Xav
  54     (
  55         (pMax().m()*pMax().X()+pMin().m()*pMin().X())
  56        /(pMax().m() + pMin().m())
  57     );
  58     scalar sigma = spray_.fuels().sigma(pc, averageTemp, Xav);
  59     sigma = max(1.0e-6, sigma);
  60     scalar rho = spray_.fuels().rho(pc, averageTemp, Xav);
  61
  62     scalar WeColl = max(1.0e-12, 0.5*rho*magVRel*magVRel*dMin/sigma);
  63
  64     scalar coalesceProb = min(1.0, 2.4*f/WeColl);
  65     scalar prob = rndGen_.sample01<scalar>();
  66
  67     // Coalescence
  68     if (prob < coalesceProb && coalescence_)
  69     {
  70         // How 'many' of the droplets coalesce
  71         // This is the kiva way ... which actually works best
  72
  73         scalar zz = collProb;
  74         scalar vnu = nu*collProb;
  75         label n=2;
  76
  77         // xx > collProb=zz
  78         while ((zz < xx) && (n<1000))
  79         {
  80             zz += vnu;
  81             vnu *= nu/n;
  82             n++;
  83         }
  84         scalar nProb = n - 1;
  85
  86         // All droplets coalesce
  87         if (nProb*nMax > nMin)
  88         {
  89             nProb = nMin/nMax;
  90         }
  91
  92         // Conservation of mass, momentum and energy
  93         pMin().m() -= nProb*nMax*mdMin;
  94
  95         scalar newMinMass = pMin().m();
  96         scalar newMaxMass = mMax + (mMin - newMinMass);
  97         pMax().m() = newMaxMass;
  98
  99         pMax().T() = (averageTemp*mTot - newMinMass*pMin().T())/newMaxMass;
 100         rhoMax = spray_.fuels().rho(pc, pMax().T(), pMax().X());
 101         scalar d3 = pow3(dMax) + nProb*pow3(dMin);
 102         pMax().d() = cbrt(d3);
 103         pMax().U() = (momMax + (1.0-newMinMass/mMin)*momMin)/newMaxMass;
 104
 105         // update the liquid molar fractions
 106         scalarField Ymin(spray_.fuels().Y(pMin().X()));
 107         scalarField Ymax(spray_.fuels().Y(pMax().X()));
 108         scalarField Ynew(mMax*Ymax + (mMin - newMinMass)*Ymin);
 109         scalar Wlinv = 0.0;
 110         forAll(Ynew, i)
 111         {
 112             Wlinv += Ynew[i]/spray_.fuels().properties()[i].W();
 113         }
 114         forAll(Ynew, i)
 115         {
 116             pMax().X()[i] = Ynew[i]/(spray_.fuels().properties()[i].W()*Wlinv);
 117         }
 118     }
 119     else
 120     {
 121         // Grazing collision (no coalescence)
 122
 123         scalar gf = sqrt(prob) - sqrt(coalesceProb);
 124         scalar denom = 1.0 - sqrt(coalesceProb);
 125         if (denom < 1.0e-5)
 126         {
 127             denom = 1.0;
 128         }
 129         gf /= denom;
 130
 131         // if gf negative, this means that coalescence is turned off
 132         // and these parcels should have coalesced
 133         gf = max(0.0, gf);
 134
 135         scalar rho1 = spray_.fuels().rho(pc, p1().T(), p1().X());
 136         scalar rho2 = spray_.fuels().rho(pc, p2().T(), p2().X());
 137         scalar m1 = p1().m();
 138         scalar m2 = p2().m();
 139         scalar n1 = p1().N(rho1);
 140         scalar n2 = p2().N(rho2);
 141
 142         // gf -> 1 => v1p -> p1().U() ...
 143         // gf -> 0 => v1p -> momentum/(m1+m2)
 144         vector mr = m1*v1 + m2*v2;
 145         vector v1p = (mr + m2*gf*vRel)/(m1+m2);
 146         vector v2p = (mr - m1*gf*vRel)/(m1+m2);
 147
 148         if (n1 < n2)
 149         {
 150             p1().U() = v1p;
 151             p2().U() = (n1*v2p + (n2-n1)*v2)/n2;
 152         }
 153         else
 154         {
 155             p1().U() = (n2*v1p + (n1-n2)*v1)/n1;
 156             p2().U() = v2p;
 157         }
 158     }
 159 }
 160