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) 1991-2009 OpenCFD Ltd.
   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
  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>::calc
  83 (
  84     TrackData& td,
  85     const scalar dt,
  86     const label cellI
  87 )
  88 {
  89     // Define local properties at beginning of time step
  90     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  91     const scalar np0 = this->nParticle_;
  92     const scalar d0 = this->d_;
  93     const vector U0 = this->U_;
  94     const scalar rho0 = this->rho_;
  95     const scalar T0 = this->T_;
  96     const scalar cp0 = this->cp_;
  97     const scalar mass0 = this->mass();
  98
  99     // Explicit momentum source for particle
 100     vector Su = vector::zero;
 101
 102     // Momentum transfer from the particle to the carrier phase
 103     vector dUTrans = vector::zero;
 104
 105     // Explicit enthalpy source for particle
 106     scalar Sh = 0.0;
 107
 108     // Sensible enthalpy transfer from the particle to the carrier phase
 109     scalar dhsTrans = 0.0;
 110
 111     // Heat transfer
 112     // ~~~~~~~~~~~~~
 113
 114     // Calculate new particle velocity
 115     scalar T1 =
 116         calcHeatTransfer(td, dt, cellI, d0, U0, rho0, T0, cp0, Sh, dhsTrans);
 117
 118
 119     // Motion
 120     // ~~~~~~
 121
 122     // Calculate new particle velocity
 123     vector U1 = calcVelocity(td, dt, cellI, d0, U0, rho0, mass0, Su, dUTrans);
 124
 125
 126     //  Accumulate carrier phase source terms
 127     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 128     if (td.cloud().coupled())
 129     {
 130         // Update momentum transfer
 131         td.cloud().UTrans()[cellI] += np0*dUTrans;
 132
 133         // Update sensible enthalpy transfer
 134         td.cloud().hsTrans()[cellI] += np0*dhsTrans;
 135     }
 136
 137     // Set new particle properties
 138     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
 139     this->U_ = U1;
 140     T_ = T1;
 141 }
 142
 143
 144 template<class ParcelType>
 145 template <class TrackData>
 146 Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
 147 (
 148     TrackData& td,
 149     const scalar dt,
 150     const label cellI,
 151     const scalar d,
 152     const vector& U,
 153     const scalar rho,
 154     const scalar T,
 155     const scalar cp,
 156     const scalar Sh,
 157     scalar& dhsTrans
 158 )
 159 {
 160     if (!td.cloud().heatTransfer().active())
 161     {
 162         return T;
 163     }
 164
 165     // Calc heat transfer coefficient
 166     scalar htc = td.cloud().heatTransfer().h
 167     (
 168         d,
 169         U - this->Uc_,
 170         this->rhoc_,
 171         rho,
 172         cpc_,
 173         cp,
 174         this->muc_
 175     );
 176
 177     const scalar As = this->areaS(d);
 178
 179     if (mag(htc) < ROOTVSMALL && !td.cloud().radiation())
 180     {
 181         return  T + dt*Sh/(this->volume(d)*rho*cp);
 182     }
 183
 184     scalar ap;
 185     scalar bp;
 186
 187     if (td.cloud().radiation())
 188     {
 189         const scalarField& G =
 190             td.cloud().mesh().objectRegistry::lookupObject<volScalarField>("G");
 191         const scalar sigma = radiation::sigmaSB.value();
 192         const scalar epsilon = td.constProps().epsilon0();
 193
 194         ap =
 195             (Sh/As + htc*Tc_ + epsilon*G[cellI]/4.0)
 196            /(htc + epsilon*sigma*pow3(T));
 197
 198         bp =
 199             6.0
 200            *(Sh/As + htc*(Tc_ - T) + epsilon*(G[cellI]/4.0 - sigma*pow4(T)))
 201            /(rho*d*cp*(ap - T));
 202     }
 203     else
 204     {
 205         ap = Tc_ + Sh/As/htc;
 206         bp = 6.0*(Sh/As + htc*(Tc_ - T))/(rho*d*cp*(ap - T));
 207     }
 208
 209     // Integrate to find the new parcel temperature
 210     IntegrationScheme<scalar>::integrationResult Tres =
 211         td.cloud().TIntegrator().integrate(T, dt, ap, bp);
 212
 213     dhsTrans += dt*htc*As*(Tres.average() - Tc_);
 214
 215     return Tres.value();
 216 }
 217
 218
 219 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 220
 221 template <class ParcelType>
 222 Foam::ThermoParcel<ParcelType>::ThermoParcel
 223 (
 224     const ThermoParcel<ParcelType>& p
 225 )
 226 :
 227     KinematicParcel<ParcelType>(p),
 228     T_(p.T_),
 229     cp_(p.cp_),
 230     Tc_(p.Tc_),
 231     cpc_(p.cpc_)
 232 {}
 233
 234
 235 // * * * * * * * * * * * * * * IOStream operators  * * * * * * * * * * * * * //
 236
 237 #include "ThermoParcelIO.C"
 238
 239 // ************************************************************************* //