FreeFOAM The Cross-Platform CFD Toolkit
Hosted by SourceForge:
Get FreeFOAM at SourceForge.net.
            Fast, secure and Free Open Source software downloads

parcel.C

Go to the documentation of this file.
00001 /*---------------------------------------------------------------------------*\
00002   =========                 |
00003   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
00004    \\    /   O peration     |
00005     \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
00006      \\/     M anipulation  |
00007 -------------------------------------------------------------------------------
00008 License
00009     This file is part of OpenFOAM.
00010 
00011     OpenFOAM is free software: you can redistribute it and/or modify it
00012     under the terms of the GNU General Public License as published by
00013     the Free Software Foundation, either version 3 of the License, or
00014     (at your option) any later version.
00015 
00016     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
00017     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00019     for more details.
00020 
00021     You should have received a copy of the GNU General Public License
00022     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
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 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
00039 
00040 namespace Foam
00041 {
00042     defineParticleTypeNameAndDebug(parcel, 0);
00043     defineTemplateTypeNameAndDebug(Cloud<parcel>, 0);
00044 }
00045 
00046 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
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 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
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     // Calculate the interpolated gas properties at the position of the parcel
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     // correct the gaseous temperature for evaporated fuel
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     // set the end-time for the track
00161     scalar tEnd = (1.0 - stepFraction())*deltaT;
00162 
00163     // set the maximum time step for this parcel
00164     scalar dtMax = min
00165     (
00166         tEnd,
00167         min
00168         (
00169             tauMomentum,
00170             min
00171             (
00172                 1.0e+10*mag(min(tauEvaporation)), // evaporation is not an issue
00173                 min
00174                 (
00175                     mag(tauHeatTransfer),
00176                     mag(min(tauBoiling))
00177                 )
00178             )
00179         )
00180     )/sDB.subCycles();
00181 
00182     // prevent the number of subcycles from being too many
00183     // (10 000 seems high enough)
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     // move the parcel until there is no 'timeLeft'
00195     while (keepParcel && tEnd > SMALL && !switchProcessor)
00196     {
00197         // set the lagrangian time-step
00198         scalar dt = min(dtMax, tEnd);
00199 
00200         // remember which cell the parcel is in
00201         // since this will change if a face is hit
00202         label celli = cell();
00203         scalar p = sDB.p()[celli];
00204 
00205         // track parcel to face, or end of trajectory
00206         if (keepParcel)
00207         {
00208             // Track and adjust the time step if the trajectory is not completed
00209             dt *= trackToFace(position() + dt*U_, sDB);
00210 
00211             // Decrement the end-time acording to how much time the track took
00212             tEnd -= dt;
00213 
00214             // Set the current time-step fraction.
00215             stepFraction() = 1.0 - tEnd/deltaT;
00216 
00217             if (onBoundary()) // hit face
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         // **** calculate the lagrangian source terms ****
00234         // First we get the 'old' properties.
00235         // and then 'update' them to get the 'new'
00236         // properties.
00237         // The difference is then added to the source terms.
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         // update the parcel properties (U, T, D)
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         // Update the Spray Source Terms
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         // Remove evaporated mass from stripped mass
00295         ms() -= ms()*(oTotMass-m())/oTotMass;
00296 
00297         // remove parcel if it is 'small'
00298         if (m() < 1.0e-12)
00299         {
00300             keepParcel = false;
00301 
00302             // ... and add the removed 'stuff' to the gas
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     // calculate mean molecular weight
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     // Calculate the interpolated gas properties at the position of the parcel
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     // correct the gaseous temperature for evaporated fuel
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         // remove the tangential velocity component
00417         scalar v1 = Up & sDB.axisOfSymmetry();
00418         scalar v2 = Up & n();
00419         Ucorr     = v1*sDB.axisOfSymmetry() + v2*n();
00420 
00421         // Remove the tangential gravity component
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     // NN.
00462     // first calculate the new temperature and droplet mass,
00463     // then calculate the energy source and correct the
00464     // gaseous temperature, Tg, and mass fraction, Yfg,
00465     // to calculate the new properties for the parcel
00466     // This procedure seems to be more stable
00467     label n = 0;
00468     while ((n < sDB.evaporation().nEvapIter()) && (m() > VSMALL))
00469     {
00470         n++;
00471         // new characteristic times does not need to be calculated the first time
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             // Prediction of new gaseous temperature and fuel mass fraction
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         // for a liquid Cl \approx Cp
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             // Prevent droplet temperature to go above critial value
00570             Tnew = min(Tnew, fuels.Tc(X()));
00571 
00572             // Prevent droplet temperature to go too low
00573             // Mainly a numerical stability issue
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                 // can not go above boiling temperature
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         // Evaporate droplet!
00628         // if the droplet is NOT boiling use implicit scheme.
00629         if (sDB.evaporation().evaporation())
00630         {
00631             for(label i=0; i<Nf; i++)
00632             {
00633                 // immediately evaporate mass that has reached critical
00634                 // condition
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 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines