src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.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 "ThermoParcel.H"
  27 #include "physicoChemicalConstants.H"
  28
  29 using namespace Foam::constant;
  30
  31 // * * * * * * * * * * *  Protected Member Functions * * * * * * * * * * * * //
  32
  33 template<class ParcelType>
  34 template<class TrackData>
  35 void Foam::ThermoParcel<ParcelType>::setCellValues
  36 (
  37     TrackData& td,
  38     const scalar dt,
  39     const label cellI
  40 )
  41 {
  42     ParcelType::setCellValues(td, dt, cellI);
  43
  44     tetIndices tetIs = this->currentTetIndices();
  45
  46     Cpc_ = td.CpInterp().interpolate(this->position(), tetIs);
  47
  48     Tc_ = td.TInterp().interpolate(this->position(), tetIs);
  49
  50     if (Tc_ < td.cloud().constProps().TMin())
  51     {
  52         if (debug)
  53         {
  54             WarningIn
  55             (
  56                 "void Foam::ThermoParcel<ParcelType>::setCellValues"
  57                 "("
  58                     "TrackData&, "
  59                     "const scalar, "
  60                     "const label"
  61                 ")"
  62             )   << "Limiting observed temperature in cell " << cellI << " to "
  63                 << td.cloud().constProps().TMin() <<  nl << endl;
  64         }
  65
  66         Tc_ = td.cloud().constProps().TMin();
  67     }
  68 }
  69
  70
  71 template<class ParcelType>
  72 template<class TrackData>
  73 void Foam::ThermoParcel<ParcelType>::cellValueSourceCorrection
  74 (
  75     TrackData& td,
  76     const scalar dt,
  77     const label cellI
  78 )
  79 {
  80     this->Uc_ += td.cloud().UTrans()[cellI]/this->massCell(cellI);
  81
  82     const scalar CpMean = td.CpInterp().psi()[cellI];
  83     Tc_ += td.cloud().hsTrans()[cellI]/(CpMean*this->massCell(cellI));
  84
  85     if (Tc_ < td.cloud().constProps().TMin())
  86     {
  87         if (debug)
  88         {
  89             WarningIn
  90             (
  91                 "void Foam::ThermoParcel<ParcelType>::cellValueSourceCorrection"
  92                 "("
  93                     "TrackData&, "
  94                     "const scalar, "
  95                     "const label"
  96                 ")"
  97             )   << "Limiting observed temperature in cell " << cellI << " to "
  98                 << td.cloud().constProps().TMin() <<  nl << endl;
  99         }
 100
 101         Tc_ = td.cloud().constProps().TMin();
 102     }
 103 }
 104
 105
 106 template<class ParcelType>
 107 template<class TrackData>
 108 void Foam::ThermoParcel<ParcelType>::calcSurfaceValues
 109 (
 110     TrackData& td,
 111     const label cellI,
 112     const scalar T,
 113     scalar& Ts,
 114     scalar& rhos,
 115     scalar& mus,
 116     scalar& Pr,
 117     scalar& kappas
 118 ) const
 119 {
 120     // Surface temperature using two thirds rule
 121     Ts = (2.0*T + Tc_)/3.0;
 122
 123     if (Ts < td.cloud().constProps().TMin())
 124     {
 125         if (debug)
 126         {
 127             WarningIn
 128             (
 129                 "void Foam::ThermoParcel<ParcelType>::calcSurfaceValues"
 130                 "("
 131                     "TrackData&, "
 132                     "const label, "
 133                     "const scalar, "
 134                     "scalar&, "
 135                     "scalar&, "
 136                     "scalar&, "
 137                     "scalar&, "
 138                     "scalar&"
 139                 ") const"
 140             )   << "Limiting parcel surface temperature to "
 141                 << td.cloud().constProps().TMin() <<  nl << endl;
 142         }
 143
 144         Ts = td.cloud().constProps().TMin();
 145     }
 146
 147     // Assuming thermo props vary linearly with T for small dT
 148     const scalar TRatio = Tc_/Ts;
 149
 150     rhos = this->rhoc_*TRatio;
 151
 152     tetIndices tetIs = this->currentTetIndices();
 153     mus = td.muInterp().interpolate(this->position(), tetIs)/TRatio;
 154
 155     Pr = td.cloud().constProps().Pr();
 156     Pr = max(ROOTVSMALL, Pr);
 157
 158     kappas = Cpc_*mus/Pr;
 159     kappas = max(ROOTVSMALL, kappas);
 160 }
 161
 162
 163 template<class ParcelType>
 164 template<class TrackData>
 165 void Foam::ThermoParcel<ParcelType>::calc
 166 (
 167     TrackData& td,
 168     const scalar dt,
 169     const label cellI
 170 )
 171 {
 172     // Define local properties at beginning of time step
 173     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 174     const scalar np0 = this->nParticle_;
 175     const scalar mass0 = this->mass();
 176
 177
 178     // Calc surface values
 179     // ~~~~~~~~~~~~~~~~~~~
 180     scalar Ts, rhos, mus, Pr, kappas;
 181     calcSurfaceValues(td, cellI, this->T_, Ts, rhos, mus, Pr, kappas);
 182
 183     // Reynolds number
 184     scalar Re = this->Re(this->U_, this->d_, rhos, mus);
 185
 186
 187     // Sources
 188     // ~~~~~~~
 189
 190     // Explicit momentum source for particle
 191     vector Su = vector::zero;
 192
 193     // Linearised momentum source coefficient
 194     scalar Spu = 0.0;
 195
 196     // Momentum transfer from the particle to the carrier phase
 197     vector dUTrans = vector::zero;
 198
 199     // Explicit enthalpy source for particle
 200     scalar Sh = 0.0;
 201
 202     // Linearised enthalpy source coefficient
 203     scalar Sph = 0.0;
 204
 205     // Sensible enthalpy transfer from the particle to the carrier phase
 206     scalar dhsTrans = 0.0;
 207
 208
 209     // Heat transfer
 210     // ~~~~~~~~~~~~~
 211
 212     // Sum Ni*Cpi*Wi of emission species
 213     scalar NCpW = 0.0;
 214
 215     // Calculate new particle temperature
 216     this->T_ =
 217         this->calcHeatTransfer
 218         (
 219             td,
 220             dt,
 221             cellI,
 222             Re,
 223             Pr,
 224             kappas,
 225             NCpW,
 226             Sh,
 227             dhsTrans,
 228             Sph
 229         );
 230
 231
 232     // Motion
 233     // ~~~~~~
 234
 235     // Calculate new particle velocity
 236     this->U_ =
 237         this->calcVelocity(td, dt, cellI, Re, mus, mass0, Su, dUTrans, Spu);
 238
 239
 240     //  Accumulate carrier phase source terms
 241     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 242     if (td.cloud().solution().coupled())
 243     {
 244         // Update momentum transfer
 245         td.cloud().UTrans()[cellI] += np0*dUTrans;
 246
 247         // Update momentum transfer coefficient
 248         td.cloud().UCoeff()[cellI] += np0*Spu;
 249
 250         // Update sensible enthalpy transfer
 251         td.cloud().hsTrans()[cellI] += np0*dhsTrans;
 252
 253         // Update sensible enthalpy coefficient
 254         td.cloud().hsCoeff()[cellI] += np0*Sph;
 255     }
 256 }
 257
 258
 259 template<class ParcelType>
 260 template<class TrackData>
 261 Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
 262 (
 263     TrackData& td,
 264     const scalar dt,
 265     const label cellI,
 266     const scalar Re,
 267     const scalar Pr,
 268     const scalar kappa,
 269     const scalar NCpW,
 270     const scalar Sh,
 271     scalar& dhsTrans,
 272     scalar& Sph
 273 )
 274 {
 275     if (!td.cloud().heatTransfer().active())
 276     {
 277         return T_;
 278     }
 279
 280     const scalar d = this->d();
 281     const scalar rho = this->rho();
 282
 283     // Calc heat transfer coefficient
 284     scalar htc = td.cloud().heatTransfer().htc(d, Re, Pr, kappa, NCpW);
 285
 286     if (mag(htc) < ROOTVSMALL && !td.cloud().radiation())
 287     {
 288         return
 289             max
 290             (
 291                 T_ + dt*Sh/(this->volume(d)*rho*Cp_),
 292                 td.cloud().constProps().TMin()
 293             );
 294     }
 295
 296     htc = max(htc, ROOTVSMALL);
 297     const scalar As = this->areaS(d);
 298
 299     scalar ap = Tc_ + Sh/As/htc;
 300     scalar bp = 6.0*(Sh/As + htc*(Tc_ - T_));
 301     if (td.cloud().radiation())
 302     {
 303         tetIndices tetIs = this->currentTetIndices();
 304         const scalar Gc = td.GInterp().interpolate(this->position(), tetIs);
 305         const scalar sigma = physicoChemical::sigma.value();
 306         const scalar epsilon = td.cloud().constProps().epsilon0();
 307
 308         ap = (ap + epsilon*Gc/(4.0*htc))/(1.0 + epsilon*sigma*pow3(T_)/htc);
 309         bp += 6.0*(epsilon*(Gc/4.0 - sigma*pow4(T_)));
 310     }
 311     bp /= rho*d*Cp_*(ap - T_) + ROOTVSMALL;
 312
 313     // Integrate to find the new parcel temperature
 314     IntegrationScheme<scalar>::integrationResult Tres =
 315         td.cloud().TIntegrator().integrate(T_, dt, ap*bp, bp);
 316
 317     scalar Tnew = max(Tres.value(), td.cloud().constProps().TMin());
 318
 319     Sph = dt*htc*As;
 320
 321     dhsTrans += Sph*(0.5*(T_ + Tnew) - Tc_);
 322
 323     return Tnew;
 324 }
 325
 326
 327 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 328
 329 template<class ParcelType>
 330 Foam::ThermoParcel<ParcelType>::ThermoParcel
 331 (
 332     const ThermoParcel<ParcelType>& p
 333 )
 334 :
 335     ParcelType(p),
 336     T_(p.T_),
 337     Cp_(p.Cp_),
 338     Tc_(p.Tc_),
 339     Cpc_(p.Cpc_)
 340 {}
 341
 342
 343 template<class ParcelType>
 344 Foam::ThermoParcel<ParcelType>::ThermoParcel
 345 (
 346     const ThermoParcel<ParcelType>& p,
 347     const polyMesh& mesh
 348 )
 349 :
 350     ParcelType(p, mesh),
 351     T_(p.T_),
 352     Cp_(p.Cp_),
 353     Tc_(p.Tc_),
 354     Cpc_(p.Cpc_)
 355 {}
 356
 357
 358 // * * * * * * * * * * * * * * IOStream operators  * * * * * * * * * * * * * //
 359
 360 #include "ThermoParcelIO.C"
 361
 362 // ************************************************************************* //