src/lagrangian/dieselSpray/parcel/parcel.C

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