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 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 "ThermoParcel.H"
  28 #include "radiationConstants.H"
  29
  30 // * * * * * * * * * * *  Protected Member Functions * * * * * * * * * * * * //
  31
  32 template<class ParcelType>
  33 template<class TrackData>
  34 void Foam::ThermoParcel<ParcelType>::setCellValues
  35 (
  36     TrackData& td,
  37     const scalar dt,
  38     const label cellI
  39 )
  40 {
  41     KinematicParcel<ParcelType>::setCellValues(td, dt, cellI);
  42
  43     cpc_ = td.cpInterp().interpolate(this->position(), cellI);
  44
  45     Tc_ = td.TInterp().interpolate(this->position(), cellI);
  46
  47     if (Tc_ < td.constProps().TMin())
  48     {
  49         WarningIn
  50         (
  51             "void Foam::ThermoParcel<ParcelType>::setCellValues"
  52             "("
  53                 "TrackData&, "
  54                 "const scalar, "
  55                 "const label"
  56             ")"
  57         )   << "Limiting observed temperature in cell " << cellI << " to "
  58             << td.constProps().TMin() <<  nl << endl;
  59
  60         Tc_ = td.constProps().TMin();
  61     }
  62 }
  63
  64
  65 template<class ParcelType>
  66 template<class TrackData>
  67 void Foam::ThermoParcel<ParcelType>::cellValueSourceCorrection
  68 (
  69     TrackData& td,
  70     const scalar dt,
  71     const label cellI
  72 )
  73 {
  74     this->Uc_ += td.cloud().UTrans()[cellI]/this->massCell(cellI);
  75
  76     scalar cpMean = td.cpInterp().psi()[cellI];
  77     Tc_ += td.cloud().hsTrans()[cellI]/(cpMean*this->massCell(cellI));
  78 }
  79
  80
  81 template<class ParcelType>
  82 template<class TrackData>
  83 void Foam::ThermoParcel<ParcelType>::calcSurfaceValues
  84 (
  85     TrackData& td,
  86     const label cellI,
  87     const scalar T,
  88     scalar& Ts,
  89     scalar& rhos,
  90     scalar& mus,
  91     scalar& Pr,
  92     scalar& kappa
  93 ) const
  94 {
  95     // Surface temperature using two thirds rule
  96     Ts = (2.0*T + Tc_)/3.0;
  97
  98     // Assuming thermo props vary linearly with T for small dT
  99     scalar factor = td.TInterp().interpolate(this->position(), cellI)/Ts;
 100     rhos = this->rhoc_*factor;
 101     mus = td.muInterp().interpolate(this->position(), cellI)/factor;
 102
 103     Pr = td.constProps().Pr();
 104     kappa = cpc_*mus/Pr;
 105 }
 106
 107
 108 template<class ParcelType>
 109 template<class TrackData>
 110 void Foam::ThermoParcel<ParcelType>::calc
 111 (
 112     TrackData& td,
 113     const scalar dt,
 114     const label cellI
 115 )
 116 {
 117     // Define local properties at beginning of time step
 118     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 119     const scalar np0 = this->nParticle_;
 120     const scalar d0 = this->d_;
 121     const vector U0 = this->U_;
 122     const scalar rho0 = this->rho_;
 123     const scalar T0 = this->T_;
 124     const scalar cp0 = this->cp_;
 125     const scalar mass0 = this->mass();
 126
 127
 128     // Calc surface values
 129     // ~~~~~~~~~~~~~~~~~~~
 130     scalar Ts, rhos, mus, Pr, kappa;
 131     calcSurfaceValues(td, cellI, T0, Ts, rhos, mus, Pr, kappa);
 132
 133     // Reynolds number
 134     scalar Re = this->Re(U0, d0, rhos, mus);
 135
 136
 137     // Sources
 138     // ~~~~~~~
 139
 140     // Explicit momentum source for particle
 141     vector Su = vector::zero;
 142
 143     // Momentum transfer from the particle to the carrier phase
 144     vector dUTrans = vector::zero;
 145
 146     // Explicit enthalpy source for particle
 147     scalar Sh = 0.0;
 148
 149     // Sensible enthalpy transfer from the particle to the carrier phase
 150     scalar dhsTrans = 0.0;
 151
 152
 153     // Heat transfer
 154     // ~~~~~~~~~~~~~
 155
 156     // Sum Ni*Cpi*Wi of emission species
 157     scalar NCpW = 0.0;
 158
 159     // Calculate new particle velocity
 160     scalar T1 =
 161         calcHeatTransfer
 162         (
 163             td,
 164             dt,
 165             cellI,
 166             Re,
 167             Pr,
 168             kappa,
 169             d0,
 170             rho0,
 171             T0,
 172             cp0,
 173             NCpW,
 174             Sh,
 175             dhsTrans
 176         );
 177
 178
 179     // Motion
 180     // ~~~~~~
 181
 182     // Calculate new particle velocity
 183     vector U1 =
 184         calcVelocity(td, dt, cellI, Re, mus, d0, U0, rho0, mass0, Su, dUTrans);
 185
 186
 187     //  Accumulate carrier phase source terms
 188     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 189     if (td.cloud().coupled())
 190     {
 191         // Update momentum transfer
 192         td.cloud().UTrans()[cellI] += np0*dUTrans;
 193
 194         // Update sensible enthalpy transfer
 195         td.cloud().hsTrans()[cellI] += np0*dhsTrans;
 196     }
 197
 198     // Set new particle properties
 199     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
 200     this->U_ = U1;
 201     T_ = T1;
 202 }
 203
 204
 205 template<class ParcelType>
 206 template <class TrackData>
 207 Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
 208 (
 209     TrackData& td,
 210     const scalar dt,
 211     const label cellI,
 212     const scalar Re,
 213     const scalar Pr,
 214     const scalar kappa,
 215     const scalar d,
 216     const scalar rho,
 217     const scalar T,
 218     const scalar cp,
 219     const scalar NCpW,
 220     const scalar Sh,
 221     scalar& dhsTrans
 222 )
 223 {
 224     if (!td.cloud().heatTransfer().active())
 225     {
 226         return T;
 227     }
 228
 229     // Calc heat transfer coefficient
 230     scalar htc = td.cloud().heatTransfer().htc(d, Re, Pr, kappa, NCpW);
 231
 232     if (mag(htc) < ROOTVSMALL && !td.cloud().radiation())
 233     {
 234         return max(T + dt*Sh/(this->volume(d)*rho*cp), td.constProps().TMin());
 235     }
 236
 237     const scalar As = this->areaS(d);
 238     scalar ap = Tc_ + Sh/As/htc;
 239     scalar bp = 6.0*(Sh/As + htc*(Tc_ - T));
 240     if (td.cloud().radiation())
 241     {
 242         const scalarField& G =
 243             td.cloud().mesh().objectRegistry::lookupObject<volScalarField>("G");
 244         const scalar Gc = G[cellI];
 245         const scalar sigma = radiation::sigmaSB.value();
 246         const scalar epsilon = td.constProps().epsilon0();
 247
 248         ap = (ap + epsilon*Gc/(4.0*htc))/(1.0 + epsilon*sigma*pow3(T)/htc);
 249         bp += 6.0*(epsilon*(Gc/4.0 - sigma*pow4(T)));
 250     }
 251     bp /= rho*d*cp*(ap - T);
 252
 253     // Integrate to find the new parcel temperature
 254     IntegrationScheme<scalar>::integrationResult Tres =
 255         td.cloud().TIntegrator().integrate(T, dt, ap, bp);
 256
 257     scalar Tnew = max(Tres.value(), td.constProps().TMin());
 258
 259     dhsTrans += dt*htc*As*(0.5*(T + Tnew) - Tc_);
 260
 261     return Tnew;
 262 }
 263
 264
 265 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 266
 267 template <class ParcelType>
 268 Foam::ThermoParcel<ParcelType>::ThermoParcel
 269 (
 270     const ThermoParcel<ParcelType>& p
 271 )
 272 :
 273     KinematicParcel<ParcelType>(p),
 274     T_(p.T_),
 275     cp_(p.cp_),
 276     Tc_(p.Tc_),
 277     cpc_(p.cpc_)
 278 {}
 279
 280
 281 // * * * * * * * * * * * * * * IOStream operators  * * * * * * * * * * * * * //
 282
 283 #include "ThermoParcelIO.C"
 284
 285 // ************************************************************************* //