FreeFOAM The Cross-Platform CFD Toolkit
trajectoryCM.H
Go to the documentation of this file.
1 vector v1 = p1().U();
2 vector v2 = p2().U();
3 scalar pc = spray_.p()[p1().cell()];
4 
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 ...
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
204 
205 // ************************ vim: set sw=4 sts=4 et: ************************ //