src/lagrangian/dieselSpray/parcel/parcel.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "parcel.H"
  27
  28 #include "spray.H"
  29 #include "dragModel.H"
  30 #include "evaporationModel.H"
  31 #include "heatTransferModel.H"
  32 #include "wallModel.H"
  33 #include "wallPolyPatch.H"
  34 #include "wedgePolyPatch.H"
  35 #include "processorPolyPatch.H"
  36 #include "basicMultiComponentMixture.H"
  37
  38 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  39
  40 Foam::parcel::parcel
  41 (
  42     const polyMesh& mesh,
  43     const vector& position,
  44     const label cellI,
  45     const label tetFaceI,
  46     const label tetPtI,
  47     const vector& n,
  48     const scalar d,
  49     const scalar T,
  50     const scalar m,
  51     const scalar y,
  52     const scalar yDot,
  53     const scalar ct,
  54     const scalar ms,
  55     const scalar tTurb,
  56     const scalar liquidCore,
  57     const scalar injector,
  58     const vector& U,
  59     const vector& Uturb,
  60     const scalarField& X,
  61     const List<word>& liquidNames
  62 )
  63 :
  64     particle(mesh, position, cellI, tetFaceI, tetPtI),
  65     liquidComponents_(liquidNames),
  66     d_(d),
  67     T_(T),
  68     m_(m),
  69     y_(y),
  70     yDot_(yDot),
  71     ct_(ct),
  72     ms_(ms),
  73     tTurb_(tTurb),
  74     liquidCore_(liquidCore),
  75     injector_(injector),
  76     U_(U),
  77     Uturb_(Uturb),
  78     n_(n),
  79     X_(X),
  80     tMom_(GREAT)
  81 {}
  82
  83
  84 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  85
  86 bool Foam::parcel::move(trackingData& td, const scalar trackTime)
  87 {
  88     td.switchProcessor = false;
  89     td.keepParticle = true;
  90
  91     const polyBoundaryMesh& pbMesh = mesh_.boundaryMesh();
  92
  93     spray& sDB = td.cloud();
  94     const liquidMixtureProperties& fuels = sDB.fuels();
  95
  96     label Nf = fuels.components().size();
  97     label Ns = sDB.composition().Y().size();
  98
  99     tetIndices tetIs = this->currentTetIndices();
 100
 101     // Calculate the interpolated gas properties at the position of the parcel
 102     vector Up = sDB.UInterpolator().interpolate(position(), tetIs) + Uturb();
 103     scalar rhog = sDB.rhoInterpolator().interpolate(position(), tetIs);
 104     scalar pg = sDB.pInterpolator().interpolate(position(), tetIs);
 105     scalar Tg = sDB.TInterpolator().interpolate(position(), tetIs);
 106
 107     scalarField Yfg(Nf, 0.0);
 108
 109     scalar cpMixture = 0.0;
 110     for (label i=0; i<Ns; i++)
 111     {
 112         const volScalarField& Yi = sDB.composition().Y()[i];
 113         if (sDB.isLiquidFuel()[i])
 114         {
 115             label j = sDB.gasToLiquidIndex()[i];
 116             scalar YicellI = Yi[cell()];
 117             Yfg[j] = YicellI;
 118         }
 119         cpMixture += Yi[cell()]*sDB.gasProperties()[i].Cp(Tg);
 120     }
 121
 122     // correct the gaseous temperature for evaporated fuel
 123
 124     scalar cellV = sDB.mesh().V()[cell()];
 125     scalar cellMass = rhog*cellV;
 126     Tg += sDB.shs()[cell()]/(cpMixture*cellMass);
 127     Tg = max(200.0, Tg);
 128
 129     scalar tauMomentum = GREAT;
 130     scalar tauHeatTransfer = GREAT;
 131     scalarField tauEvaporation(Nf, GREAT);
 132     scalarField tauBoiling(Nf, GREAT);
 133
 134     setRelaxationTimes
 135     (
 136         cell(),
 137         tauMomentum,
 138         tauEvaporation,
 139         tauHeatTransfer,
 140         tauBoiling,
 141         sDB,
 142         rhog,
 143         Up,
 144         Tg,
 145         pg,
 146         Yfg,
 147         m()*fuels.Y(X()),
 148         trackTime
 149     );
 150
 151
 152     // set the end-time for the track
 153     scalar tEnd = (1.0 - stepFraction())*trackTime;
 154
 155     // set the maximum time step for this parcel
 156     scalar dtMax = min
 157     (
 158         tEnd,
 159         min
 160         (
 161             tauMomentum,
 162             min
 163             (
 164                 1.0e+10*mag(min(tauEvaporation)), // evaporation is not an issue
 165                 min
 166                 (
 167                     mag(tauHeatTransfer),
 168                     mag(min(tauBoiling))
 169                 )
 170             )
 171         )
 172     )/sDB.subCycles();
 173
 174     // prevent the number of subcycles from being too many
 175     // (10 000 seems high enough)
 176     dtMax = max(dtMax, 1.0e-4*tEnd);
 177
 178     vector planeNormal = vector::zero;
 179     if (sDB.twoD())
 180     {
 181         planeNormal = n() ^ sDB.axisOfSymmetry();
 182         planeNormal /= mag(planeNormal);
 183     }
 184
 185     // move the parcel until there is no 'timeLeft'
 186     while (td.keepParticle && !td.switchProcessor && tEnd > SMALL)
 187     {
 188         // set the lagrangian time-step
 189         scalar dt = min(dtMax, tEnd);
 190
 191         // remember which cell the parcel is in
 192         // since this will change if a face is hit
 193         label cellI = cell();
 194         scalar p = sDB.p()[cellI];
 195
 196         // track parcel to face, or end of trajectory
 197         if (td.keepParticle)
 198         {
 199             // Track and adjust the time step if the trajectory is not completed
 200             dt *= trackToFace(position() + dt*U_, td);
 201
 202             // Decrement the end-time acording to how much time the track took
 203             tEnd -= dt;
 204
 205             // Set the current time-step fraction.
 206             stepFraction() = 1.0 - tEnd/trackTime;
 207
 208             if (onBoundary()) // hit face
 209             {
 210                 #include "boundaryTreatment.H"
 211             }
 212         }
 213
 214         if (td.keepParticle && sDB.twoD())
 215         {
 216             scalar z = position() & sDB.axisOfSymmetry();
 217             vector r = position() - z*sDB.axisOfSymmetry();
 218             if (mag(r) > SMALL)
 219             {
 220                 correctNormal(sDB.axisOfSymmetry());
 221             }
 222         }
 223
 224         // **** calculate the lagrangian source terms ****
 225         // First we get the 'old' properties.
 226         // and then 'update' them to get the 'new'
 227         // properties.
 228         // The difference is then added to the source terms.
 229
 230         scalar oRho = fuels.rho(p, T(), X());
 231         scalarField oMass(Nf, 0.0);
 232         scalar oHg = 0.0;
 233         scalar oTotMass = m();
 234         scalarField oYf(fuels.Y(X()));
 235
 236         forAll(oMass, i)
 237         {
 238             oMass[i] = m()*oYf[i];
 239             label j = sDB.liquidToGasIndex()[i];
 240             oHg += oYf[i]*sDB.gasProperties()[j].Hs(T());
 241         }
 242
 243         vector oMom = m()*U();
 244         scalar oHv = fuels.hl(p, T(), X());
 245         scalar oH = oHg - oHv;
 246         scalar oPE = (p - fuels.pv(p, T(), X()))/oRho;
 247
 248         // update the parcel properties (U, T, D)
 249         updateParcelProperties
 250         (
 251             dt,
 252             sDB,
 253             cellI,
 254             face()
 255         );
 256
 257         scalar nRho = fuels.rho(p, T(), X());
 258         scalar nHg = 0.0;
 259         scalarField nMass(Nf, 0.0);
 260         scalarField nYf(fuels.Y(X()));
 261
 262         forAll(nMass, i)
 263         {
 264             nMass[i] = m()*nYf[i];
 265             label j = sDB.liquidToGasIndex()[i];
 266             nHg += nYf[i]*sDB.gasProperties()[j].Hs(T());
 267         }
 268
 269         vector nMom = m()*U();
 270         scalar nHv = fuels.hl(p, T(), X());
 271         scalar nH = nHg - nHv;
 272         scalar nPE = (p - fuels.pv(p, T(), X()))/nRho;
 273
 274         // Update the Spray Source Terms
 275         forAll(nMass, i)
 276         {
 277             sDB.srhos()[i][cellI] += oMass[i] - nMass[i];
 278         }
 279         sDB.sms()[cellI] += oMom - nMom;
 280
 281         sDB.shs()[cellI] += oTotMass*(oH + oPE) - m()*(nH + nPE);
 282
 283         // Remove evaporated mass from stripped mass
 284         ms() -= ms()*(oTotMass - m())/oTotMass;
 285
 286         // remove parcel if it is 'small'
 287         if (m() < 1.0e-12)
 288         {
 289             td.keepParticle = false;
 290
 291             // ... and add the removed 'stuff' to the gas
 292             forAll(nMass, i)
 293             {
 294                 sDB.srhos()[i][cellI] += nMass[i];
 295             }
 296
 297             sDB.sms()[cellI] += nMom;
 298             sDB.shs()[cellI] += m()*(nH + nPE);
 299         }
 300
 301         if (onBoundary() && td.keepParticle)
 302         {
 303             if (face() > -1)
 304             {
 305                 if (isA<processorPolyPatch>(pbMesh[patch(face())]))
 306                 {
 307                     td.switchProcessor = true;
 308                 }
 309             }
 310         }
 311     }
 312
 313     return td.keepParticle;
 314 }
 315
 316
 317 void Foam::parcel::updateParcelProperties
 318 (
 319     const scalar dt,
 320     spray& sDB,
 321     const label cellI,
 322     const label faceI
 323 )
 324 {
 325     const liquidMixtureProperties& fuels = sDB.fuels();
 326
 327     label Nf = sDB.fuels().components().size();
 328     label Ns = sDB.composition().Y().size();
 329
 330     // calculate mean molecular weight
 331     scalar W = 0.0;
 332     for (label i=0; i<Ns; i++)
 333     {
 334         W += sDB.composition().Y()[i][cellI]/sDB.gasProperties()[i].W();
 335
 336     }
 337     W = 1.0/W;
 338
 339     // Calculate the interpolated gas properties at the position of the parcel
 340     vector Up = sDB.UInterpolator().interpolate(position(), cellI, faceI)
 341         + Uturb();
 342     scalar rhog = sDB.rhoInterpolator().interpolate(position(), cellI, faceI);
 343     scalar pg = sDB.pInterpolator().interpolate(position(), cellI, faceI);
 344     scalar Tg0 = sDB.TInterpolator().interpolate(position(), cellI, faceI);
 345
 346     // correct the gaseous temperature for evaporated fuel
 347     scalar cpMix = 0.0;
 348     for (label i=0; i<Ns; i++)
 349     {
 350         cpMix += sDB.composition().Y()[i][cellI]
 351                 *sDB.gasProperties()[i].Cp(T());
 352     }
 353     scalar cellV            = sDB.mesh().V()[cellI];
 354     scalar rho              = sDB.rho()[cellI];
 355     scalar cellMass         = rho*cellV;
 356     scalar dh               = sDB.shs()[cellI];
 357     scalarField addedMass(Nf, 0.0);
 358
 359     forAll(addedMass, i)
 360     {
 361         addedMass[i] += sDB.srhos()[i][cellI]*cellV;
 362     }
 363
 364     scalar Tg  = Tg0 + dh/(cpMix*cellMass);
 365     Tg = max(200.0, Tg);
 366
 367     scalarField Yfg(Nf, 0.0);
 368     forAll(Yfg, i)
 369     {
 370         label j = sDB.liquidToGasIndex()[i];
 371         const volScalarField& Yj = sDB.composition().Y()[j];
 372         scalar Yfg0 = Yj[cellI];
 373         Yfg[i] = (Yfg0*cellMass + addedMass[i])/(addedMass[i] + cellMass);
 374     }
 375
 376     scalar tauMomentum     = GREAT;
 377     scalar tauHeatTransfer = GREAT;
 378     scalarField tauEvaporation(Nf, GREAT);
 379     scalarField tauBoiling(Nf, GREAT);
 380
 381     setRelaxationTimes
 382     (
 383         cellI,
 384         tauMomentum,
 385         tauEvaporation,
 386         tauHeatTransfer,
 387         tauBoiling,
 388         sDB,
 389         rhog,
 390         Up,
 391         Tg,
 392         pg,
 393         Yfg,
 394         m()*fuels.Y(X()),
 395         dt
 396     );
 397
 398     scalar timeRatio = dt/tauMomentum;
 399
 400     vector Ucorr = Up;
 401     vector gcorr = sDB.g();
 402
 403     if (sDB.twoD())
 404     {
 405         // remove the tangential velocity component
 406         scalar v1 = Up & sDB.axisOfSymmetry();
 407         scalar v2 = Up & n();
 408         Ucorr     = v1*sDB.axisOfSymmetry() + v2*n();
 409
 410         // Remove the tangential gravity component
 411         scalar g1 = gcorr & sDB.axisOfSymmetry();
 412         scalar g2 = gcorr & n();
 413         gcorr     = g1*sDB.axisOfSymmetry() + g2*n();
 414     }
 415
 416     U() = (U() + timeRatio*Ucorr + gcorr*dt)/(1.0 + timeRatio);
 417
 418     if (sDB.twoD())
 419     {
 420         vector normal = n() ^ sDB.axisOfSymmetry();
 421         normal /= mag(normal);
 422         scalar dU = U() & normal;
 423         U() -= dU*normal;
 424     }
 425
 426     scalar TDroplet = T();
 427     scalar oldDensity = fuels.rho(pg, T(), X());
 428     scalar oldMass = m();
 429     scalarField Yf0(fuels.Y(X()));
 430     scalarField mi0(m()*Yf0);
 431     scalarField mi(mi0);
 432
 433     scalar oldhg = 0.0;
 434     for (label i=0; i<Nf; i++)
 435     {
 436         label j = sDB.liquidToGasIndex()[i];
 437         oldhg += Yf0[i]*sDB.gasProperties()[j].Hs(T());
 438     }
 439
 440     scalar oldhv = fuels.hl(pg, T(), X());
 441     scalar Np = N(oldDensity);
 442
 443     scalar newMass    = oldMass;
 444     scalar newhg      = oldhg;
 445     scalar newhv      = oldhv;
 446
 447     scalar Tnew = T();
 448
 449     // first calculate the new temperature and droplet mass,
 450     // then calculate the energy source and correct the
 451     // gaseous temperature, Tg, and mass fraction, Yfg,
 452     // to calculate the new properties for the parcel
 453     // This procedure seems to be more stable
 454     label n = 0;
 455     while ((n < sDB.evaporation().nEvapIter()) && (m() > VSMALL))
 456     {
 457         n++;
 458         // new characteristic times does not need to be calculated the
 459         // first time
 460         if (n > 1)
 461         {
 462             newMass = m();
 463             newhg = 0.0;
 464             scalarField Ynew(fuels.Y(X()));
 465             for (label i=0; i<Nf; i++)
 466             {
 467                 label j = sDB.liquidToGasIndex()[i];
 468                 newhg += Ynew[i]*sDB.gasProperties()[j].Hs(Tnew);
 469             }
 470
 471             newhv = fuels.hl(pg, Tnew, X());
 472
 473             scalar dm = oldMass - newMass;
 474             scalar dhNew = oldMass*(oldhg-oldhv) - newMass*(newhg-newhv);
 475
 476             // Prediction of new gaseous temperature and fuel mass fraction
 477             Tg  = Tg0 + (dh+dhNew)/(cpMix*cellMass);
 478             Tg = max(200.0, Tg);
 479
 480             forAll(Yfg, i)
 481             {
 482                 label j = sDB.liquidToGasIndex()[i];
 483                 const volScalarField& Yj = sDB.composition().Y()[j];
 484                 scalar Yfg0 = Yj[cellI];
 485                 Yfg[i] = (Yfg0*cellMass + addedMass[i] + dm)
 486                         /(addedMass[i] + cellMass + dm);
 487             }
 488
 489             setRelaxationTimes
 490             (
 491                 cellI,
 492                 tauMomentum,
 493                 tauEvaporation,
 494                 tauHeatTransfer,
 495                 tauBoiling,
 496                 sDB,
 497                 rhog,
 498                 Up,
 499                 Tg,
 500                 pg,
 501                 Yfg,
 502                 m()*fuels.Y(X()),
 503                 dt
 504             );
 505
 506         }
 507
 508         scalar Taverage = TDroplet + (Tg - TDroplet)/3.0;
 509         // for a liquid Cl \approx Cp
 510         scalar liquidcL = sDB.fuels().Cp(pg, TDroplet, X());
 511
 512         cpMix = 0.0;
 513         for (label i=0; i<Ns; i++)
 514         {
 515             if (sDB.isLiquidFuel()[i])
 516             {
 517                 label j = sDB.gasToLiquidIndex()[i];
 518                 cpMix += Yfg[j]*sDB.gasProperties()[i].Cp(Taverage);
 519             }
 520             else
 521             {
 522                 scalar Y = sDB.composition().Y()[i][cellI];
 523                 cpMix += Y*sDB.gasProperties()[i].Cp(Taverage);
 524             }
 525         }
 526
 527         scalar evaporationSource = 0.0;
 528         scalar z = 0.0;
 529         scalar tauEvap = 0.0;
 530
 531         for (label i=0; i<Nf; i++)
 532         {
 533             tauEvap += X()[i]*fuels.properties()[i].W()/tauEvaporation[i];
 534         }
 535         tauEvap = fuels.W(X())/tauEvap;
 536
 537
 538         if (sDB.evaporation().evaporation())
 539         {
 540             scalar hv = fuels.hl(pg, TDroplet, X());
 541             evaporationSource =
 542                 hv/liquidcL/tauEvap;
 543
 544             z = cpMix*tauHeatTransfer/liquidcL/tauEvap;
 545         }
 546
 547         if (sDB.heatTransfer().heatTransfer())
 548         {
 549             scalar fCorrect =
 550                 sDB.heatTransfer().fCorrection(z)/tauHeatTransfer;
 551
 552             Tnew =
 553                 (TDroplet + dt*(fCorrect * Tg - evaporationSource))
 554                 /(1.0 + dt*fCorrect);
 555
 556             // Prevent droplet temperature to go above critial value
 557             Tnew = min(Tnew, fuels.Tc(X()));
 558
 559             // Prevent droplet temperature to go too low
 560             // Mainly a numerical stability issue
 561             Tnew = max(200.0, Tnew);
 562             scalar Td = Tnew;
 563
 564             scalar pAtSurface = fuels.pv(pg, Td, X());
 565             scalar pCompare = 0.999*pg;
 566             scalar boiling = pAtSurface >= pCompare;
 567             if (boiling)
 568             {
 569                 // can not go above boiling temperature
 570                 scalar Terr = 1.0e-3;
 571                 label n = 0;
 572                 scalar dT = 1.0;
 573                 scalar pOld = pAtSurface;
 574                 while (dT > Terr)
 575                 {
 576                     n++;
 577                     pAtSurface = fuels.pv(pg, Td, X());
 578                     if ((pAtSurface < pCompare) && (pOld < pCompare))
 579                     {
 580                         Td += dT;
 581                     }
 582                     else
 583                     {
 584                         if ((pAtSurface > pCompare) && (pOld > pCompare))
 585                         {
 586                             Td -= dT;
 587                         }
 588                         else
 589                         {
 590                             dT *= 0.5;
 591                             if ((pAtSurface > pCompare) && (pOld < pCompare))
 592                             {
 593                                 Td -= dT;
 594                             }
 595                             else
 596                             {
 597                                 Td += dT;
 598                             }
 599                         }
 600                     }
 601                     pOld = pAtSurface;
 602                     if (debug)
 603                     {
 604                         if (n>100)
 605                         {
 606                             Info<< "n = " << n << ", T = " << Td << ", pv = "
 607                                 << pAtSurface << endl;
 608                         }
 609                     }
 610                 }
 611                 Tnew = Td;
 612             }
 613         }
 614
 615         // Evaporate droplet!
 616         // if the droplet is NOT boiling use implicit scheme.
 617         if (sDB.evaporation().evaporation())
 618         {
 619             for (label i=0; i<Nf; i++)
 620             {
 621                 // immediately evaporate mass that has reached critical
 622                 // condition
 623                 if (mag(Tnew - fuels.Tc(X())) < SMALL)
 624                 {
 625                     mi[i] = 0.0;
 626                 }
 627                 else
 628                 {
 629                     scalar Td = min(Tnew, 0.999*fuels.properties()[i].Tc());
 630
 631                     scalar pAtSurface = fuels.properties()[i].pv(pg, Td);
 632                     scalar boiling = pAtSurface >= 0.999*pg;
 633
 634                     if (!boiling)
 635                     {
 636                         scalar fr = dt/tauEvaporation[i];
 637                         mi[i] = mi0[i]/(1.0 + fr);
 638                     }
 639                     else
 640                     {
 641                         scalar fr = dt/tauBoiling[i];
 642                         mi[i] = mi0[i]/(1.0 + fr);
 643                     }
 644                 }
 645             }
 646
 647             scalar mTot = sum(mi);
 648             if (mTot > VSMALL)
 649             {
 650                 scalarField Ynew(mi/mTot);
 651                 scalarField Xnew(sDB.fuels().X(Ynew));
 652                 forAll(Xnew, i)
 653                 {
 654                     X()[i] = Xnew[i];
 655                 }
 656                 m() = mTot;
 657             }
 658             else
 659             {
 660                 m() = 0.0;
 661             }
 662         }
 663         T() = Tnew;
 664         scalar rhod = fuels.rho(pg, T(), X());
 665         d() = cbrt(6.0*m()/(Np*rhod*M_PI));
 666     }
 667
 668     T() = Tnew;
 669     scalar rhod = fuels.rho(pg, T(), X());
 670     m() = sum(mi);
 671     d() = cbrt(6.0*m()/(Np*rhod*M_PI));
 672 }
 673
 674
 675 void Foam::parcel::transformProperties(const tensor& T)
 676 {
 677     U_ = transform(T, U_);
 678 }
 679
 680
 681 void Foam::parcel::transformProperties(const vector&)
 682 {}
 683
 684
 685 // ************************************************************************* //