src/lagrangian/dieselSpray/spraySubModels/collisionModel/trajectoryModel/trajectoryCM.H

   1 vector v1 = p1().U();
   2 vector v2 = p2().U();
   3 scalar pc = spray_.p()[p1().cell()];
   4
   5 vector vRel = v1 - v2;
   6 scalar magVRel = mag(vRel);
   7
   8 vector p = p2().position() - p1().position();
   9 scalar dist = mag(p);
  10
  11 scalar vAlign = vRel & (p/(dist+SMALL));
  12
  13 if (vAlign > 0)
  14 {
  15     scalar sumD = p1().d() + p2().d();
  16
  17     if (vAlign*dt > dist - 0.5*sumD)
  18     {
  19         scalar v1Mag = mag(v1);
  20         scalar v2Mag = mag(v2);
  21         vector nv1 = v1/v1Mag;
  22         vector nv2 = v2/v2Mag;
  23
  24         scalar v1v2 = nv1 & nv2;
  25         scalar v1p = nv1 & p;
  26         scalar v2p = nv2 & p;
  27
  28         scalar det = 1.0 - v1v2*v1v2;
  29
  30         scalar alpha = 1.0e+20;
  31         scalar beta = 1.0e+20;
  32
  33         if (mag(det) > 1.0e-4)
  34         {
  35             beta = -(v2p - v1v2*v1p)/det;
  36             alpha = v1p + v1v2*beta;
  37         }
  38
  39         alpha /= v1Mag*dt;
  40         beta /= v2Mag*dt;
  41
  42         // is collision possible within this timestep
  43         if ((alpha>0) && (alpha<1.0) && (beta>0) && (beta<1.0))
  44         {
  45             vector p1c = p1().position() + alpha*v1*dt;
  46             vector p2c = p2().position() + beta*v2*dt;
  47
  48             scalar closestDist = mag(p1c-p2c);
  49
  50             scalar collProb =
  51                 pow(0.5*sumD/max(0.5*sumD, closestDist), cSpace_)
  52                *exp(-cTime_*mag(alpha-beta));
  53
  54             scalar xx = rndGen_.scalar01();
  55
  56             spray::iterator pMin = p1;
  57             spray::iterator pMax = p2;
  58
  59             scalar dMin = pMin().d();
  60             scalar dMax = pMax().d();
  61
  62             if (dMin > dMax)
  63             {
  64                 dMin = pMax().d();
  65                 dMax = pMin().d();
  66                 pMin = p2;
  67                 pMax = p1;
  68             }
  69
  70             scalar rhoMax = spray_.fuels().rho(pc, pMax().T(), pMax().X());
  71             scalar rhoMin = spray_.fuels().rho(pc, pMin().T(), pMin().X());
  72             scalar mMax = pMax().m();
  73             scalar mMin = pMin().m();
  74             scalar nMax = pMax().N(rhoMax);
  75             scalar nMin = pMin().N(rhoMin);
  76
  77             scalar mdMin = mMin/nMin;
  78
  79             // collision occur
  80             if ((xx < collProb) && (mMin > VSMALL) && (mMax > VSMALL))
  81             {
  82                 scalar mTot = mMax + mMin;
  83
  84                 scalar gamma = dMax/max(dMin, 1.0e-12);
  85                 scalar f = gamma*gamma*gamma + 2.7*gamma - 2.4*gamma*gamma;
  86
  87                 vector momMax = mMax*pMax().U();
  88                 vector momMin = mMin*pMin().U();
  89
  90                 // use mass-averaged temperature to calculate We number
  91                 scalar averageTemp = (pMax().T()*mMax + pMin().T()*mMin)/mTot;
  92                 // and mass averaged mole fractions ...
  93                 scalarField
  94                     Xav((pMax().m()*pMax().X()+pMin().m()*pMin().X())
  95                    /(pMax().m() + pMin().m()));
  96
  97                 scalar sigma = spray_.fuels().sigma(pc, averageTemp, Xav);
  98                 sigma = max(1.0e-6, sigma);
  99                 scalar rho = spray_.fuels().rho(pc, averageTemp, Xav);
 100
 101                 scalar dMean = sqrt(dMin*dMax);
 102                 scalar WeColl =
 103                     max(1.0e-12, 0.5*rho*magVRel*magVRel*dMean/sigma);
 104
 105                 // coalescence only possible when parcels are close enoug
 106
 107                 scalar coalesceProb = min(1.0, 2.4*f/WeColl);
 108
 109                 scalar prob = rndGen_.scalar01();
 110
 111                 // Coalescence
 112                 if ( prob < coalesceProb && coalescence_)
 113                 {
 114                     // How 'many' of the droplets coalesce
 115                     scalar nProb = prob*nMin/nMax;
 116
 117                     // Conservation of mass, momentum and energy
 118
 119                     pMin().m() -= nMax*nProb*mdMin;
 120
 121                     scalar newMinMass = pMin().m();
 122                     scalar newMaxMass = mMax + (mMin - newMinMass);
 123                     pMax().m() = newMaxMass;
 124
 125                     pMax().T() =
 126                         (averageTemp*mTot - newMinMass*pMin().T())/newMaxMass;
 127                     rhoMax = spray_.fuels().rho(pc, pMax().T(), pMax().X());
 128
 129                     pMax().d() =
 130                         pow
 131                         (
 132                             6.0*newMaxMass/(rhoMax*mathematicalConstant::pi*nMax),
 133                             1.0/3.0
 134                         );
 135
 136                     pMax().U() =
 137                         (momMax + (1.0-newMinMass/mMin)*momMin)/newMaxMass;
 138
 139                     // update the liquid molar fractions
 140                     scalarField Ymin = spray_.fuels().Y(pMin().X());
 141                     scalarField Ymax = spray_.fuels().Y(pMax().X());
 142                     scalarField Ynew = mMax*Ymax + (mMin - newMinMass)*Ymin;
 143                     scalar Wlinv = 0.0;
 144                     forAll(Ynew, i)
 145                     {
 146                         Wlinv += Ynew[i]/spray_.fuels().properties()[i].W();
 147                     }
 148                     forAll(Ynew, i)
 149                     {
 150                         pMax().X()[i] =
 151                             Ynew[i]/(spray_.fuels().properties()[i].W()*Wlinv);
 152                     }
 153
 154
 155                 }
 156                 // Grazing collision (no coalescence)
 157                 else
 158                 {
 159                     scalar gf = sqrt(prob) - sqrt(coalesceProb);
 160                     scalar denom = 1.0 - sqrt(coalesceProb);
 161                     if (denom < 1.0e-5) {
 162                         denom = 1.0;
 163                     }
 164                     gf /= denom;
 165
 166                     // if gf negative, this means that coalescence is turned off
 167                     // and these parcels should have coalesced
 168                     gf = max(0.0, gf);
 169
 170                     scalar rho1 = spray_.fuels().rho(pc, p1().T(), p1().X());
 171                     scalar rho2 = spray_.fuels().rho(0.0, p2().T(), p2().X());
 172                     scalar m1 = p1().m();
 173                     scalar m2 = p2().m();
 174                     scalar n1 = p1().N(rho1);
 175                     scalar n2 = p2().N(rho2);
 176
 177                     // gf -> 1 => v1p -> p1().U() ...
 178                     // gf -> 0 => v1p -> momentum/(m1+m2)
 179
 180                     vector mr = m1*v1 + m2*v2;
 181                     vector v1p = (mr + m2*gf*vRel)/(m1+m2);
 182                     vector v2p = (mr - m1*gf*vRel)/(m1+m2);
 183
 184                     if (n1 < n2)
 185                     {
 186                         p1().U() = v1p;
 187                         p2().U() = (n1*v2p + (n2-n1)*v2)/n2;
 188                     }
 189                     else
 190                     {
 191                         p1().U() = (n2*v1p + (n1-n2)*v1)/n1;
 192                         p2().U() = v2p;
 193                     }
 194
 195                 } // if - coalescence or not
 196
 197             } // if - collision
 198
 199         } // if - possible collision (alpha, beta) in timeinterval
 200
 201     } // if - travelled distance is larger distance between parcels
 202
 203 } // if - parcel travel towards eachother