00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #include "parcel.H"
00027
00028 #include <dieselSpray/spray.H>
00029 #include <dieselSpray/dragModel.H>
00030 #include <dieselSpray/evaporationModel.H>
00031 #include <dieselSpray/heatTransferModel.H>
00032 #include <dieselSpray/wallModel.H>
00033 #include <OpenFOAM/wallPolyPatch.H>
00034 #include <OpenFOAM/wedgePolyPatch.H>
00035 #include <OpenFOAM/processorPolyPatch.H>
00036 #include <reactionThermophysicalModels/basicMultiComponentMixture.H>
00037
00038
00039
00040 namespace Foam
00041 {
00042 defineParticleTypeNameAndDebug(parcel, 0);
00043 defineTemplateTypeNameAndDebug(Cloud<parcel>, 0);
00044 }
00045
00046
00047
00048 Foam::parcel::parcel
00049 (
00050 const Cloud<parcel>& cloud,
00051 const vector& position,
00052 const label cellI,
00053 const vector& n,
00054 const scalar d,
00055 const scalar T,
00056 const scalar m,
00057 const scalar y,
00058 const scalar yDot,
00059 const scalar ct,
00060 const scalar ms,
00061 const scalar tTurb,
00062 const scalar liquidCore,
00063 const scalar injector,
00064 const vector& U,
00065 const vector& Uturb,
00066 const scalarField& X,
00067 const List<word>& liquidNames
00068 )
00069 :
00070 Particle<parcel>(cloud, position, cellI),
00071 liquidComponents_
00072 (
00073 liquidNames
00074 ),
00075 d_(d),
00076 T_(T),
00077 m_(m),
00078 y_(y),
00079 yDot_(yDot),
00080 ct_(ct),
00081 ms_(ms),
00082 tTurb_(tTurb),
00083 liquidCore_(liquidCore),
00084 injector_(injector),
00085 U_(U),
00086 Uturb_(Uturb),
00087 n_(n),
00088 X_(X),
00089 tMom_(GREAT)
00090 {}
00091
00092
00093
00094
00095 bool Foam::parcel::move(spray& sDB)
00096 {
00097 const polyMesh& mesh = cloud().pMesh();
00098 const polyBoundaryMesh& pbMesh = mesh.boundaryMesh();
00099
00100 const liquidMixture& fuels = sDB.fuels();
00101
00102 scalar deltaT = sDB.runTime().deltaT().value();
00103 label Nf = fuels.components().size();
00104 label Ns = sDB.composition().Y().size();
00105
00106
00107 vector Up = sDB.UInterpolator().interpolate(position(), cell())
00108 + Uturb();
00109 scalar rhog = sDB.rhoInterpolator().interpolate(position(), cell());
00110 scalar pg = sDB.pInterpolator().interpolate(position(), cell());
00111 scalar Tg = sDB.TInterpolator().interpolate(position(), cell());
00112
00113 scalarField Yfg(Nf, 0.0);
00114
00115 scalar cpMixture = 0.0;
00116 for(label i=0; i<Ns; i++)
00117 {
00118 const volScalarField& Yi = sDB.composition().Y()[i];
00119 if (sDB.isLiquidFuel()[i])
00120 {
00121 label j = sDB.gasToLiquidIndex()[i];
00122 scalar Yicelli = Yi[cell()];
00123 Yfg[j] = Yicelli;
00124 }
00125 cpMixture += Yi[cell()]*sDB.gasProperties()[i].Cp(Tg);
00126 }
00127
00128
00129
00130 scalar cellV = sDB.mesh().V()[cell()];
00131 scalar cellMass = rhog*cellV;
00132 Tg += sDB.shs()[cell()]/(cpMixture*cellMass);
00133 Tg = max(200.0, Tg);
00134
00135 scalar tauMomentum = GREAT;
00136 scalar tauHeatTransfer = GREAT;
00137 scalarField tauEvaporation(Nf, GREAT);
00138 scalarField tauBoiling(Nf, GREAT);
00139
00140 bool keepParcel = true;
00141
00142 setRelaxationTimes
00143 (
00144 cell(),
00145 tauMomentum,
00146 tauEvaporation,
00147 tauHeatTransfer,
00148 tauBoiling,
00149 sDB,
00150 rhog,
00151 Up,
00152 Tg,
00153 pg,
00154 Yfg,
00155 m()*fuels.Y(X()),
00156 deltaT
00157 );
00158
00159
00160
00161 scalar tEnd = (1.0 - stepFraction())*deltaT;
00162
00163
00164 scalar dtMax = min
00165 (
00166 tEnd,
00167 min
00168 (
00169 tauMomentum,
00170 min
00171 (
00172 1.0e+10*mag(min(tauEvaporation)),
00173 min
00174 (
00175 mag(tauHeatTransfer),
00176 mag(min(tauBoiling))
00177 )
00178 )
00179 )
00180 )/sDB.subCycles();
00181
00182
00183
00184 dtMax = max(dtMax, 1.0e-4*tEnd);
00185
00186 bool switchProcessor = false;
00187 vector planeNormal = vector::zero;
00188 if (sDB.twoD())
00189 {
00190 planeNormal = n() ^ sDB.axisOfSymmetry();
00191 planeNormal /= mag(planeNormal);
00192 }
00193
00194
00195 while (keepParcel && tEnd > SMALL && !switchProcessor)
00196 {
00197
00198 scalar dt = min(dtMax, tEnd);
00199
00200
00201
00202 label celli = cell();
00203 scalar p = sDB.p()[celli];
00204
00205
00206 if (keepParcel)
00207 {
00208
00209 dt *= trackToFace(position() + dt*U_, sDB);
00210
00211
00212 tEnd -= dt;
00213
00214
00215 stepFraction() = 1.0 - tEnd/deltaT;
00216
00217 if (onBoundary())
00218 {
00219 # include <dieselSpray/boundaryTreatment.H>
00220 }
00221 }
00222
00223 if (keepParcel && sDB.twoD())
00224 {
00225 scalar z = position() & sDB.axisOfSymmetry();
00226 vector r = position() - z*sDB.axisOfSymmetry();
00227 if (mag(r) > SMALL)
00228 {
00229 correctNormal(sDB.axisOfSymmetry());
00230 }
00231 }
00232
00233
00234
00235
00236
00237
00238
00239 scalar oRho = fuels.rho(p, T(), X());
00240 scalarField oMass(Nf, 0.0);
00241 scalar oHg = 0.0;
00242 scalar oTotMass = m();
00243 scalarField oYf(fuels.Y(X()));
00244
00245 forAll(oMass, i)
00246 {
00247 oMass[i] = m()*oYf[i];
00248 label j = sDB.liquidToGasIndex()[i];
00249 oHg += oYf[i]*sDB.gasProperties()[j].Hs(T());
00250 }
00251
00252 vector oMom = m()*U();
00253 scalar oHv = fuels.hl(p, T(), X());
00254 scalar oH = oHg - oHv;
00255 scalar oPE = (p - fuels.pv(p, T(), X()))/oRho;
00256
00257
00258 updateParcelProperties
00259 (
00260 dt,
00261 sDB,
00262 celli,
00263 face()
00264 );
00265
00266 scalar nRho = fuels.rho(p, T(), X());
00267 scalar nHg = 0.0;
00268 scalarField nMass(Nf, 0.0);
00269 scalarField nYf(fuels.Y(X()));
00270
00271 forAll(nMass, i)
00272 {
00273 nMass[i] = m()*nYf[i];
00274 label j = sDB.liquidToGasIndex()[i];
00275 nHg += nYf[i]*sDB.gasProperties()[j].Hs(T());
00276 }
00277
00278 vector nMom = m()*U();
00279 scalar nHv = fuels.hl(p, T(), X());
00280 scalar nH = nHg - nHv;
00281 scalar nPE = (p - fuels.pv(p, T(), X()))/nRho;
00282
00283
00284 forAll(nMass, i)
00285 {
00286 sDB.srhos()[i][celli] += oMass[i] - nMass[i];
00287 }
00288 sDB.sms()[celli] += oMom - nMom;
00289
00290 sDB.shs()[celli] +=
00291 oTotMass*(oH + oPE)
00292 - m()*(nH + nPE);
00293
00294
00295 ms() -= ms()*(oTotMass-m())/oTotMass;
00296
00297
00298 if (m() < 1.0e-12)
00299 {
00300 keepParcel = false;
00301
00302
00303 forAll(nMass, i)
00304 {
00305 sDB.srhos()[i][celli] += nMass[i];
00306 }
00307
00308 sDB.sms()[celli] += nMom;
00309 sDB.shs()[celli] += m()*(nH + nPE);
00310 }
00311
00312 if (onBoundary() && keepParcel)
00313 {
00314 if (face() > -1)
00315 {
00316 if (isA<processorPolyPatch>(pbMesh[patch(face())]))
00317 {
00318 switchProcessor = true;
00319 }
00320 }
00321 }
00322 }
00323
00324 return keepParcel;
00325 }
00326
00327
00328 void Foam::parcel::updateParcelProperties
00329 (
00330 const scalar dt,
00331 spray& sDB,
00332 const label celli,
00333 const label facei
00334 )
00335 {
00336 const liquidMixture& fuels = sDB.fuels();
00337
00338 label Nf = sDB.fuels().components().size();
00339 label Ns = sDB.composition().Y().size();
00340
00341
00342 scalar W = 0.0;
00343 for(label i=0; i<Ns; i++)
00344 {
00345 W += sDB.composition().Y()[i][celli]/sDB.gasProperties()[i].W();
00346
00347 }
00348 W = 1.0/W;
00349
00350
00351 vector Up = sDB.UInterpolator().interpolate(position(), celli, facei)
00352 + Uturb();
00353 scalar rhog = sDB.rhoInterpolator().interpolate(position(), celli, facei);
00354 scalar pg = sDB.pInterpolator().interpolate(position(), celli, facei);
00355 scalar Tg0 = sDB.TInterpolator().interpolate(position(), celli, facei);
00356
00357
00358 scalar cpMix = 0.0;
00359 for(label i=0; i<Ns; i++)
00360 {
00361 cpMix += sDB.composition().Y()[i][celli]
00362 *sDB.gasProperties()[i].Cp(T());
00363 }
00364 scalar cellV = sDB.mesh().V()[celli];
00365 scalar rho = sDB.rho()[celli];
00366 scalar cellMass = rho*cellV;
00367 scalar dh = sDB.shs()[celli];
00368 scalarField addedMass(Nf, 0.0);
00369
00370 forAll(addedMass, i)
00371 {
00372 addedMass[i] += sDB.srhos()[i][celli]*cellV;
00373 }
00374
00375 scalar Tg = Tg0 + dh/(cpMix*cellMass);
00376 Tg = max(200.0, Tg);
00377
00378 scalarField Yfg(Nf, 0.0);
00379 forAll(Yfg, i)
00380 {
00381 label j = sDB.liquidToGasIndex()[i];
00382 const volScalarField& Yj = sDB.composition().Y()[j];
00383 scalar Yfg0 = Yj[celli];
00384 Yfg[i] = (Yfg0*cellMass + addedMass[i])/(addedMass[i] + cellMass);
00385 }
00386
00387 scalar tauMomentum = GREAT;
00388 scalar tauHeatTransfer = GREAT;
00389 scalarField tauEvaporation(Nf, GREAT);
00390 scalarField tauBoiling(Nf, GREAT);
00391
00392 setRelaxationTimes
00393 (
00394 celli,
00395 tauMomentum,
00396 tauEvaporation,
00397 tauHeatTransfer,
00398 tauBoiling,
00399 sDB,
00400 rhog,
00401 Up,
00402 Tg,
00403 pg,
00404 Yfg,
00405 m()*fuels.Y(X()),
00406 dt
00407 );
00408
00409 scalar timeRatio = dt/tauMomentum;
00410
00411 vector Ucorr = Up;
00412 vector gcorr = sDB.g();
00413
00414 if (sDB.twoD())
00415 {
00416
00417 scalar v1 = Up & sDB.axisOfSymmetry();
00418 scalar v2 = Up & n();
00419 Ucorr = v1*sDB.axisOfSymmetry() + v2*n();
00420
00421
00422 scalar g1 = gcorr & sDB.axisOfSymmetry();
00423 scalar g2 = gcorr & n();
00424 gcorr = g1*sDB.axisOfSymmetry() + g2*n();
00425 }
00426
00427 U() = (U() + timeRatio*Ucorr + gcorr*dt)/(1.0 + timeRatio);
00428
00429 if (sDB.twoD())
00430 {
00431 vector normal = n() ^ sDB.axisOfSymmetry();
00432 normal /= mag(normal);
00433 scalar dU = U() & normal;
00434 U() -= dU*normal;
00435 }
00436
00437 scalar TDroplet = T();
00438 scalar oldDensity = fuels.rho(pg, T(), X());
00439 scalar oldMass = m();
00440 scalarField Yf0(fuels.Y(X()));
00441 scalarField mi0(m()*Yf0);
00442 scalarField mi(mi0);
00443
00444 scalar oldhg = 0.0;
00445 for (label i=0; i<Nf; i++)
00446 {
00447 label j = sDB.liquidToGasIndex()[i];
00448 oldhg += Yf0[i]*sDB.gasProperties()[j].Hs(T());
00449 }
00450
00451 scalar oldhv = fuels.hl(pg, T(), X());
00452 scalar Np = N(oldDensity);
00453
00454 scalar newDensity = oldDensity;
00455 scalar newMass = oldMass;
00456 scalar newhg = oldhg;
00457 scalar newhv = oldhv;
00458
00459 scalar Tnew = T();
00460
00461
00462
00463
00464
00465
00466
00467 label n = 0;
00468 while ((n < sDB.evaporation().nEvapIter()) && (m() > VSMALL))
00469 {
00470 n++;
00471
00472 if (n > 1)
00473 {
00474 newDensity = fuels.rho(pg, Tnew, X());
00475 newMass = m();
00476 newhg = 0.0;
00477 scalarField Ynew(fuels.Y(X()));
00478 for(label i=0; i<Nf; i++)
00479 {
00480 label j = sDB.liquidToGasIndex()[i];
00481 newhg += Ynew[i]*sDB.gasProperties()[j].Hs(Tnew);
00482 }
00483
00484 newhv = fuels.hl(pg, Tnew, X());
00485
00486 scalar dm = oldMass - newMass;
00487 scalar dhNew = oldMass*(oldhg-oldhv) - newMass*(newhg-newhv);
00488
00489
00490 Tg = Tg0 + (dh+dhNew)/(cpMix*cellMass);
00491 Tg = max(200.0, Tg);
00492
00493 forAll(Yfg, i)
00494 {
00495 label j = sDB.liquidToGasIndex()[i];
00496 const volScalarField& Yj = sDB.composition().Y()[j];
00497 scalar Yfg0 = Yj[celli];
00498 Yfg[i] = (Yfg0*cellMass + addedMass[i] + dm)
00499 /(addedMass[i] + cellMass + dm);
00500 }
00501
00502 setRelaxationTimes
00503 (
00504 celli,
00505 tauMomentum,
00506 tauEvaporation,
00507 tauHeatTransfer,
00508 tauBoiling,
00509 sDB,
00510 rhog,
00511 Up,
00512 Tg,
00513 pg,
00514 Yfg,
00515 m()*fuels.Y(X()),
00516 dt
00517 );
00518
00519 }
00520
00521 scalar Taverage = TDroplet + (Tg - TDroplet)/3.0;
00522
00523 scalar liquidcL = sDB.fuels().cp(pg, TDroplet, X());
00524
00525 cpMix = 0.0;
00526 for(label i=0; i<Ns; i++)
00527 {
00528 if (sDB.isLiquidFuel()[i])
00529 {
00530 label j = sDB.gasToLiquidIndex()[i];
00531 cpMix += Yfg[j]*sDB.gasProperties()[i].Cp(Taverage);
00532 }
00533 else
00534 {
00535 scalar Y = sDB.composition().Y()[i][celli];
00536 cpMix += Y*sDB.gasProperties()[i].Cp(Taverage);
00537 }
00538 }
00539
00540 scalar evaporationSource = 0.0;
00541 scalar z = 0.0;
00542 scalar tauEvap = 0.0;
00543
00544 for(label i=0; i<Nf; i++)
00545 {
00546 tauEvap += X()[i]*fuels.properties()[i].W()/tauEvaporation[i];
00547 }
00548 tauEvap = fuels.W(X())/tauEvap;
00549
00550
00551 if (sDB.evaporation().evaporation())
00552 {
00553 scalar hv = fuels.hl(pg, TDroplet, X());
00554 evaporationSource =
00555 hv/liquidcL/tauEvap;
00556
00557 z = cpMix*tauHeatTransfer/liquidcL/tauEvap;
00558 }
00559
00560 if (sDB.heatTransfer().heatTransfer())
00561 {
00562 scalar fCorrect =
00563 sDB.heatTransfer().fCorrection(z)/tauHeatTransfer;
00564
00565 Tnew =
00566 (TDroplet + dt*(fCorrect * Tg - evaporationSource))
00567 /(1.0 + dt*fCorrect);
00568
00569
00570 Tnew = min(Tnew, fuels.Tc(X()));
00571
00572
00573
00574 Tnew = max(200.0, Tnew);
00575 scalar Td = Tnew;
00576
00577 scalar pAtSurface = fuels.pv(pg, Td, X());
00578 scalar pCompare = 0.999*pg;
00579 scalar boiling = pAtSurface >= pCompare;
00580 if (boiling)
00581 {
00582
00583 scalar Terr = 1.0e-3;
00584 label n=0;
00585 scalar dT = 1.0;
00586 scalar pOld = pAtSurface;
00587 while (dT > Terr)
00588 {
00589 n++;
00590 pAtSurface = fuels.pv(pg, Td, X());
00591 if ((pAtSurface < pCompare) && (pOld < pCompare))
00592 {
00593 Td += dT;
00594 }
00595 else
00596 {
00597 if ((pAtSurface > pCompare) && (pOld > pCompare))
00598 {
00599 Td -= dT;
00600 }
00601 else
00602 {
00603 dT *= 0.5;
00604 if ((pAtSurface > pCompare) && (pOld < pCompare))
00605 {
00606 Td -= dT;
00607 }
00608 else
00609 {
00610 Td += dT;
00611 }
00612 }
00613 }
00614 pOld = pAtSurface;
00615 if (debug)
00616 {
00617 if (n>100)
00618 {
00619 Info << "n = " << n << ", T = " << Td << ", pv = " << pAtSurface << endl;
00620 }
00621 }
00622 }
00623 Tnew = Td;
00624 }
00625 }
00626
00627
00628
00629 if (sDB.evaporation().evaporation())
00630 {
00631 for(label i=0; i<Nf; i++)
00632 {
00633
00634
00635 if (mag(Tnew - fuels.Tc(X())) < SMALL)
00636 {
00637 mi[i] = 0.0;
00638 }
00639 else
00640 {
00641 scalar Td = min(Tnew, 0.999*fuels.properties()[i].Tc());
00642
00643 scalar pAtSurface = fuels.properties()[i].pv(pg, Td);
00644 scalar boiling = pAtSurface >= 0.999*pg;
00645
00646 if (!boiling)
00647 {
00648 scalar fr = dt/tauEvaporation[i];
00649 mi[i] = mi0[i]/(1.0 + fr);
00650 }
00651 else
00652 {
00653 scalar fr = dt/tauBoiling[i];
00654 mi[i] = mi0[i]/(1.0 + fr);
00655 }
00656 }
00657 }
00658
00659 scalar mTot = sum(mi);
00660 if (mTot > VSMALL)
00661 {
00662 scalarField Ynew(mi/mTot);
00663 scalarField Xnew(sDB.fuels().X(Ynew));
00664 forAll(Xnew, i)
00665 {
00666 X()[i] = Xnew[i];
00667 }
00668 m() = mTot;
00669 }
00670 else
00671 {
00672 m() = 0.0;
00673 }
00674 }
00675 T() = Tnew;
00676 scalar rhod = fuels.rho(pg, T(), X());
00677 d() = cbrt(6.0*m()/(Np*rhod*M_PI));
00678 }
00679
00680 T() = Tnew;
00681 scalar rhod = fuels.rho(pg, T(), X());
00682 m() = sum(mi);
00683 d() = cbrt(6.0*m()/(Np*rhod*M_PI));
00684 }
00685
00686
00687 void Foam::parcel::transformProperties(const tensor& T)
00688 {
00689 U_ = transform(T, U_);
00690 }
00691
00692
00693 void Foam::parcel::transformProperties(const vector&)
00694 {}
00695
00696
00697