src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.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 "SprayParcel.H"
  27 #include "CompositionModel.H"
  28 #include "AtomizationModel.H"
  29
  30 // * * * * * * * * * * *  Protected Member Functions * * * * * * * * * * * * //
  31
  32 template<class ParcelType>
  33 template<class TrackData>
  34 void Foam::SprayParcel<ParcelType>::setCellValues
  35 (
  36     TrackData& td,
  37     const scalar dt,
  38     const label cellI
  39 )
  40 {
  41     ParcelType::setCellValues(td, dt, cellI);
  42 }
  43
  44
  45 template<class ParcelType>
  46 template<class TrackData>
  47 void Foam::SprayParcel<ParcelType>::cellValueSourceCorrection
  48 (
  49     TrackData& td,
  50     const scalar dt,
  51     const label cellI
  52 )
  53 {
  54     ParcelType::cellValueSourceCorrection(td, dt, cellI);
  55 }
  56
  57
  58 template<class ParcelType>
  59 template<class TrackData>
  60 void Foam::SprayParcel<ParcelType>::calc
  61 (
  62     TrackData& td,
  63     const scalar dt,
  64     const label cellI
  65 )
  66 {
  67     bool coupled = td.cloud().solution().coupled();
  68
  69     // check if parcel belongs to liquid core
  70     if (liquidCore() > 0.5)
  71     {
  72         // liquid core parcels should not interact with the gas
  73         if (td.cloud().solution().coupled())
  74         {
  75             td.cloud().solution().coupled() = false;
  76         }
  77     }
  78
  79     // store the parcel properties
  80     const scalarField& Y(this->Y());
  81     scalarField X(td.cloud().composition().liquids().X(Y));
  82
  83     scalar T0 = this->T();
  84     this->Cp() = td.cloud().composition().liquids().Cp(this->pc_, T0, X);
  85     scalar rho0 = td.cloud().composition().liquids().rho(this->pc_, T0, X);
  86     this->rho() = rho0;
  87
  88     ParcelType::calc(td, dt, cellI);
  89
  90     if (td.keepParticle)
  91     {
  92         // update Cp, diameter and density due to change in temperature
  93         // and/or composition
  94         scalar T1 = this->T();
  95         const scalarField& Y1(this->Y());
  96         scalarField X1(td.cloud().composition().liquids().X(Y1));
  97
  98         this->Cp() = td.cloud().composition().liquids().Cp(this->pc_, T1, X1);
  99
 100         scalar rho1 = td.cloud().composition().liquids().rho(this->pc_, T1, X1);
 101         this->rho() = rho1;
 102         scalar d1 = this->d()*cbrt(rho0/rho1);
 103         this->d() = d1;
 104
 105         if (liquidCore() > 0.5)
 106         {
 107             calcAtomization(td, dt, cellI);
 108
 109             // preserve the total mass/volume by increasing the number of
 110             // particles in parcels due to breakup
 111             scalar d2 = this->d();
 112             this->nParticle() *= pow3(d1/d2);
 113         }
 114         else
 115         {
 116             calcBreakup(td, dt, cellI);
 117         }
 118     }
 119
 120     // restore coupled
 121     td.cloud().solution().coupled() = coupled;
 122 }
 123
 124
 125 template<class ParcelType>
 126 template<class TrackData>
 127 void Foam::SprayParcel<ParcelType>::calcAtomization
 128 (
 129     TrackData& td,
 130     const scalar dt,
 131     const label cellI
 132 )
 133 {
 134     typedef typename TrackData::cloudType::reactingCloudType reactingCloudType;
 135     const CompositionModel<reactingCloudType>& composition =
 136         td.cloud().composition();
 137
 138     typedef typename TrackData::cloudType::sprayCloudType sprayCloudType;
 139     const AtomizationModel<sprayCloudType>& atomization =
 140         td.cloud().atomization();
 141
 142
 143     // cell state info is updated in ReactingParcel calc
 144     const scalarField& Y(this->Y());
 145     scalarField X(composition.liquids().X(Y));
 146
 147     scalar rho = composition.liquids().rho(this->pc(), this->T(), X);
 148     scalar mu = composition.liquids().mu(this->pc(), this->T(), X);
 149     scalar sigma = composition.liquids().sigma(this->pc(), this->T(), X);
 150
 151     // Average molecular weight of carrier mix - assumes perfect gas
 152     scalar Wc = this->rhoc_*specie::RR*this->Tc()/this->pc();
 153     scalar R = specie::RR/Wc;
 154     scalar Tav = atomization.Taverage(this->T(), this->Tc());
 155
 156     // calculate average gas density based on average temperature
 157     scalar rhoAv = this->pc()/(R*Tav);
 158
 159     scalar soi = td.cloud().injection().timeStart();
 160     scalar currentTime = td.cloud().db().time().value();
 161     const vector& pos = this->position();
 162     const vector& injectionPos = this->position0();
 163
 164     // disregard the continous phase when calculating the relative velocity
 165     // (in line with the deactivated coupled assumption)
 166     scalar Urel = mag(this->U());
 167
 168     scalar t0 = max(0.0, currentTime - this->age() - soi);
 169     scalar t1 = min(t0 + dt, td.cloud().injection().timeEnd() - soi);
 170     // this should be the vol flow rate from when the parcel was injected
 171     scalar volFlowRate = td.cloud().injection().volumeToInject(t0, t1)/dt;
 172
 173     scalar chi = 0.0;
 174     if (atomization.calcChi())
 175     {
 176         chi = this->chi(td, X);
 177     }
 178
 179     atomization.update
 180     (
 181         dt,
 182         this->d(),
 183         this->liquidCore(),
 184         this->tc(),
 185         rho,
 186         mu,
 187         sigma,
 188         volFlowRate,
 189         rhoAv,
 190         Urel,
 191         pos,
 192         injectionPos,
 193         td.cloud().pAmbient(),
 194         chi,
 195         td.cloud().rndGen()
 196     );
 197 }
 198
 199
 200 template<class ParcelType>
 201 template<class TrackData>
 202 void Foam::SprayParcel<ParcelType>::calcBreakup
 203 (
 204     TrackData& td,
 205     const scalar dt,
 206     const label cellI
 207 )
 208 {
 209     typedef typename TrackData::cloudType::reactingCloudType reactingCloudType;
 210     const CompositionModel<reactingCloudType>& composition =
 211         td.cloud().composition();
 212
 213     typedef typename TrackData::cloudType cloudType;
 214     typedef typename cloudType::parcelType parcelType;
 215     typedef typename cloudType::forceType forceType;
 216
 217     const parcelType& p = static_cast<const parcelType&>(*this);
 218     const forceType& forces = td.cloud().forces();
 219
 220     if (td.cloud().breakup().solveOscillationEq())
 221     {
 222         solveTABEq(td, dt);
 223     }
 224
 225     // cell state info is updated in ReactingParcel calc
 226     const scalarField& Y(this->Y());
 227     scalarField X(composition.liquids().X(Y));
 228
 229     scalar rho = composition.liquids().rho(this->pc(), this->T(), X);
 230     scalar mu = composition.liquids().mu(this->pc(), this->T(), X);
 231     scalar sigma = composition.liquids().sigma(this->pc(), this->T(), X);
 232
 233     // Average molecular weight of carrier mix - assumes perfect gas
 234     scalar Wc = this->rhoc()*specie::RR*this->Tc()/this->pc();
 235     scalar R = specie::RR/Wc;
 236     scalar Tav = td.cloud().atomization().Taverage(this->T(), this->Tc());
 237
 238     // calculate average gas density based on average temperature
 239     scalar rhoAv = this->pc()/(R*Tav);
 240     scalar muAv = this->muc();
 241     vector Urel = this->U() - this->Uc();
 242     scalar Urmag = mag(Urel);
 243     scalar Re = this->Re(this->U(), this->d(), rhoAv, muAv);
 244
 245     const scalar mass = p.mass();
 246     const forceSuSp Fcp = forces.calcCoupled(p, dt, mass, Re, muAv);
 247     const forceSuSp Fncp = forces.calcNonCoupled(p, dt, mass, Re, muAv);
 248     scalar tMom = 1.0/(Fcp.Sp() + Fncp.Sp());
 249
 250     const vector g = td.cloud().g().value();
 251
 252     scalar massChild = 0.0;
 253     scalar dChild = 0.0;
 254     if
 255     (
 256         td.cloud().breakup().update
 257         (
 258             dt,
 259             g,
 260             this->d(),
 261             this->tc(),
 262             this->ms(),
 263             this->nParticle(),
 264             this->KHindex(),
 265             this->y(),
 266             this->yDot(),
 267             this->d0(),
 268             rho,
 269             mu,
 270             sigma,
 271             this->U(),
 272             rhoAv,
 273             muAv,
 274             Urel,
 275             Urmag,
 276             tMom,
 277             td.cloud().averageParcelMass(),
 278             dChild,
 279             massChild,
 280             td.cloud().rndGen()
 281         )
 282     )
 283     {
 284         scalar Re = rhoAv*Urmag*dChild/muAv;
 285         this->mass0() -= massChild;
 286
 287         // Add child parcel as copy of parent
 288         SprayParcel<ParcelType>* child = new SprayParcel<ParcelType>(*this);
 289         child->mass0() = massChild;
 290         child->d() = dChild;
 291         child->nParticle() = massChild/rho*this->volume(dChild);
 292
 293         const forceSuSp Fcp =
 294             forces.calcCoupled(*child, dt, massChild, Re, muAv);
 295         const forceSuSp Fncp =
 296             forces.calcNonCoupled(*child, dt, massChild, Re, muAv);
 297
 298         child->liquidCore() = 0.0;
 299         child->KHindex() = 1.0;
 300         child->y() = td.cloud().breakup().y0();
 301         child->yDot() = td.cloud().breakup().yDot0();
 302         child->tc() = -GREAT;
 303         child->ms() = 0.0;
 304         child->injector() = this->injector();
 305         child->tMom() = 1.0/(Fcp.Sp() + Fncp.Sp());
 306         child->user() = 0.0;
 307         child->setCellValues(td, dt, cellI);
 308
 309         td.cloud().addParticle(child);
 310     }
 311 }
 312
 313
 314 template<class ParcelType>
 315 template<class TrackData>
 316 Foam::scalar Foam::SprayParcel<ParcelType>::chi
 317 (
 318     TrackData& td,
 319     const scalarField& X
 320 ) const
 321 {
 322     // modifications to take account of the flash boiling on primary break-up
 323
 324     static label nIter = 200;
 325
 326     typedef typename TrackData::cloudType::reactingCloudType reactingCloudType;
 327     const CompositionModel<reactingCloudType>& composition =
 328         td.cloud().composition();
 329
 330     scalar chi = 0.0;
 331     scalar T0 = this->T();
 332     scalar Tc0 = this->Tc();
 333     scalar p0 = this->pc();
 334     scalar pAmb = td.cloud().pAmbient();
 335
 336     scalar pv = composition.liquids().pv(p0, T0, X);
 337
 338     forAll(composition.liquids(), i)
 339     {
 340         if (pv >= 0.999*pAmb)
 341         {
 342             // liquid is boiling - calc boiling temperature
 343
 344             const liquidProperties& liq = composition.liquids().properties()[i];
 345             scalar TBoil = T0;
 346
 347             for (label k=0; k<nIter; k++)
 348             {
 349                 scalar pBoil = liq.pv(p0, TBoil);
 350
 351                 if (pBoil > p0)
 352                 {
 353                     TBoil = TBoil - (T0 - Tc0)/nIter;
 354                 }
 355                 else
 356                 {
 357                     break;
 358                 }
 359             }
 360
 361             scalar hl = liq.hl(pAmb, TBoil);
 362             scalar iTp = liq.h(pAmb, T0) - liq.rho(pAmb, T0);
 363             scalar iTb = liq.h(pAmb, TBoil) - pAmb/liq.rho(pAmb, TBoil);
 364
 365             chi += X[i]*(iTp - iTb)/hl;
 366         }
 367     }
 368
 369     chi = min(1.0, max(chi, 0.0));
 370
 371     return chi;
 372 }
 373
 374
 375 template<class ParcelType>
 376 template<class TrackData>
 377 void Foam::SprayParcel<ParcelType>::solveTABEq
 378 (
 379     TrackData& td,
 380     const scalar dt
 381 )
 382 {
 383     typedef typename TrackData::cloudType::reactingCloudType reactingCloudType;
 384     const CompositionModel<reactingCloudType>& composition =
 385         td.cloud().composition();
 386
 387     const scalar& TABCmu = td.cloud().breakup().TABCmu();
 388     const scalar& TABWeCrit = td.cloud().breakup().TABWeCrit();
 389     const scalar& TABComega = td.cloud().breakup().TABComega();
 390
 391     scalar r = 0.5*this->d();
 392     scalar r2 = r*r;
 393     scalar r3 = r*r2;
 394
 395     const scalarField& Y(this->Y());
 396     scalarField X(composition.liquids().X(Y));
 397
 398     scalar rho = composition.liquids().rho(this->pc(), this->T(), X);
 399     scalar mu = composition.liquids().mu(this->pc(), this->T(), X);
 400     scalar sigma = composition.liquids().sigma(this->pc(), this->T(), X);
 401
 402     // inverse of characteristic viscous damping time
 403     scalar rtd = 0.5*TABCmu*mu/(rho*r2);
 404
 405     // oscillation frequency (squared)
 406     scalar omega2 = TABComega*sigma/(rho*r3) - rtd*rtd;
 407
 408     if(omega2 > 0)
 409     {
 410         scalar omega = sqrt(omega2);
 411         scalar rhoc = this->rhoc();
 412         scalar Wetmp = this->We(this->U(), r, rhoc, sigma)/TABWeCrit;
 413
 414         scalar y1 = this->y() - Wetmp;
 415         scalar y2 = this->yDot()/omega;
 416
 417         // update distortion parameters
 418         scalar c = cos(omega*dt);
 419         scalar s = sin(omega*dt);
 420         scalar e = exp(-rtd*dt);
 421         y2 = (this->yDot() + y1*rtd)/omega;
 422
 423         this->y() = Wetmp + e*(y1*c + y2*s);
 424         if (this->y() < 0)
 425         {
 426             this->y() = 0.0;
 427             this->yDot() = 0.0;
 428         }
 429         else
 430         {
 431             this->yDot() = (Wetmp - this->y())*rtd + e*omega*(y2*c - y1*s);
 432         }
 433     }
 434     else
 435     {
 436         // reset distortion parameters
 437         this->y() = 0;
 438         this->yDot() = 0;
 439     }
 440 }
 441
 442
 443 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 444
 445 template <class ParcelType>
 446 Foam::SprayParcel<ParcelType>::SprayParcel(const SprayParcel<ParcelType>& p)
 447 :
 448     ParcelType(p),
 449     d0_(p.d0_),
 450     position0_(p.position0_),
 451     liquidCore_(p.liquidCore_),
 452     KHindex_(p.KHindex_),
 453     y_(p.y_),
 454     yDot_(p.yDot_),
 455     tc_(p.tc_),
 456     ms_(p.ms_),
 457     injector_(p.injector_),
 458     tMom_(p.tMom_),
 459     user_(p.user_)
 460 {}
 461
 462
 463 template <class ParcelType>
 464 Foam::SprayParcel<ParcelType>::SprayParcel
 465 (
 466     const SprayParcel<ParcelType>& p,
 467     const polyMesh& mesh
 468 )
 469 :
 470     ParcelType(p, mesh),
 471     d0_(p.d0_),
 472     position0_(p.position0_),
 473     liquidCore_(p.liquidCore_),
 474     KHindex_(p.KHindex_),
 475     y_(p.y_),
 476     yDot_(p.yDot_),
 477     tc_(p.tc_),
 478     ms_(p.ms_),
 479     injector_(p.injector_),
 480     tMom_(p.tMom_),
 481     user_(p.user_)
 482 {}
 483
 484
 485 // * * * * * * * * * * * * * * IOStream operators  * * * * * * * * * * * * * //
 486
 487 #include "SprayParcelIO.C"
 488
 489
 490 // ************************************************************************* //