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