src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C

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