FreeFOAM The Cross-Platform CFD Toolkit
sameCell.H
Go to the documentation of this file.
1 vector v1 = p1().U();
2 vector v2 = p2().U();
3 
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 
151 
152 // ************************ vim: set sw=4 sts=4 et: ************************ //