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