ThermoParcel.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2018 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"
28 
29 using namespace Foam::constant;
30 
31 // * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * * //
32 
33 template<class ParcelType>
34 template<class TrackCloudType>
36 (
37  TrackCloudType& cloud,
38  trackingData& td
39 )
40 {
41  ParcelType::setCellValues(cloud, td);
42 
43  tetIndices tetIs = this->currentTetIndices();
44 
45  td.Cpc() = td.CpInterp().interpolate(this->coordinates(), tetIs);
46 
47  td.Tc() = td.TInterp().interpolate(this->coordinates(), tetIs);
48 
49  if (td.Tc() < cloud.constProps().TMin())
50  {
51  if (debug)
52  {
54  << "Limiting observed temperature in cell " << this->cell()
55  << " to " << cloud.constProps().TMin() << nl << endl;
56  }
57 
58  td.Tc() = cloud.constProps().TMin();
59  }
60 }
61 
62 
63 template<class ParcelType>
64 template<class TrackCloudType>
66 (
67  TrackCloudType& cloud,
68  trackingData& td,
69  const scalar dt
70 )
71 {
72  td.Uc() += cloud.UTrans()[this->cell()]/this->massCell(td);
73 
74  const scalar CpMean = td.CpInterp().psi()[this->cell()];
75  td.Tc() += cloud.hsTrans()[this->cell()]/(CpMean*this->massCell(td));
76 
77  if (td.Tc() < cloud.constProps().TMin())
78  {
79  if (debug)
80  {
82  << "Limiting observed temperature in cell " << this->cell()
83  << " to " << cloud.constProps().TMin() << nl << endl;
84  }
85 
86  td.Tc() = cloud.constProps().TMin();
87  }
88 }
89 
90 
91 template<class ParcelType>
92 template<class TrackCloudType>
94 (
95  TrackCloudType& cloud,
96  trackingData& td,
97  const scalar T,
98  scalar& Ts,
99  scalar& rhos,
100  scalar& mus,
101  scalar& Pr,
102  scalar& kappas
103 ) const
104 {
105  // Surface temperature using two thirds rule
106  Ts = (2.0*T + td.Tc())/3.0;
107 
108  if (Ts < cloud.constProps().TMin())
109  {
110  if (debug)
111  {
113  << "Limiting parcel surface temperature to "
114  << cloud.constProps().TMin() << nl << endl;
115  }
116 
117  Ts = cloud.constProps().TMin();
118  }
119 
120  // Assuming thermo props vary linearly with T for small d(T)
121  const scalar TRatio = td.Tc()/Ts;
122 
123  rhos = td.rhoc()*TRatio;
124 
125  tetIndices tetIs = this->currentTetIndices();
126  mus = td.muInterp().interpolate(this->coordinates(), tetIs)/TRatio;
127  kappas = td.kappaInterp().interpolate(this->coordinates(), tetIs)/TRatio;
128 
129  Pr = td.Cpc()*mus/kappas;
130  Pr = max(rootVSmall, Pr);
131 }
132 
133 
134 template<class ParcelType>
135 template<class TrackCloudType>
137 (
138  TrackCloudType& cloud,
139  trackingData& td,
140  const scalar dt
141 )
142 {
143  // Define local properties at beginning of time step
144  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
145  const scalar np0 = this->nParticle_;
146  const scalar mass0 = this->mass();
147 
148  // Store T for consistent radiation source
149  const scalar T0 = this->T_;
150 
151 
152  // Calc surface values
153  // ~~~~~~~~~~~~~~~~~~~
154  scalar Ts, rhos, mus, Pr, kappas;
155  calcSurfaceValues(cloud, td, this->T_, Ts, rhos, mus, Pr, kappas);
156 
157  // Reynolds number
158  scalar Re = this->Re(rhos, this->U_, td.Uc(), this->d_, mus);
159 
160 
161  // Sources
162  // ~~~~~~~
163 
164  // Explicit momentum source for particle
165  vector Su = Zero;
166 
167  // Linearised momentum source coefficient
168  scalar Spu = 0.0;
169 
170  // Momentum transfer from the particle to the carrier phase
171  vector dUTrans = Zero;
172 
173  // Explicit enthalpy source for particle
174  scalar Sh = 0.0;
175 
176  // Linearised enthalpy source coefficient
177  scalar Sph = 0.0;
178 
179  // Sensible enthalpy transfer from the particle to the carrier phase
180  scalar dhsTrans = 0.0;
181 
182 
183  // Heat transfer
184  // ~~~~~~~~~~~~~
185 
186  // Sum Ni*Cpi*Wi of emission species
187  scalar NCpW = 0.0;
188 
189  // Calculate new particle temperature
190  this->T_ =
191  this->calcHeatTransfer
192  (
193  cloud,
194  td,
195  dt,
196  Re,
197  Pr,
198  kappas,
199  NCpW,
200  Sh,
201  dhsTrans,
202  Sph
203  );
204 
205 
206  // Motion
207  // ~~~~~~
208 
209  // Calculate new particle velocity
210  this->U_ =
211  this->calcVelocity(cloud, td, dt, Re, mus, mass0, Su, dUTrans, Spu);
212 
213 
214  // Accumulate carrier phase source terms
215  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
216  if (cloud.solution().coupled())
217  {
218  // Update momentum transfer
219  cloud.UTrans()[this->cell()] += np0*dUTrans;
220 
221  // Update momentum transfer coefficient
222  cloud.UCoeff()[this->cell()] += np0*Spu;
223 
224  // Update sensible enthalpy transfer
225  cloud.hsTrans()[this->cell()] += np0*dhsTrans;
226 
227  // Update sensible enthalpy coefficient
228  cloud.hsCoeff()[this->cell()] += np0*Sph;
229 
230  // Update radiation fields
231  if (cloud.radiation())
232  {
233  const scalar ap = this->areaP();
234  const scalar T4 = pow4(T0);
235  cloud.radAreaP()[this->cell()] += dt*np0*ap;
236  cloud.radT4()[this->cell()] += dt*np0*T4;
237  cloud.radAreaPT4()[this->cell()] += dt*np0*ap*T4;
238  }
239  }
240 }
241 
242 
243 template<class ParcelType>
244 template<class TrackCloudType>
246 (
247  TrackCloudType& cloud,
248  trackingData& td,
249  const scalar dt,
250  const scalar Re,
251  const scalar Pr,
252  const scalar kappa,
253  const scalar NCpW,
254  const scalar Sh,
255  scalar& dhsTrans,
256  scalar& Sph
257 )
258 {
259  if (!cloud.heatTransfer().active())
260  {
261  return T_;
262  }
263 
264  const scalar d = this->d();
265  const scalar rho = this->rho();
266  const scalar As = this->areaS(d);
267  const scalar V = this->volume(d);
268  const scalar m = rho*V;
269 
270  // Calc heat transfer coefficient
271  scalar htc = cloud.heatTransfer().htc(d, Re, Pr, kappa, NCpW);
272 
273  // Calculate the integration coefficients
274  const scalar bcp = htc*As/(m*Cp_);
275  const scalar acp = bcp*td.Tc();
276  scalar ancp = Sh;
277  if (cloud.radiation())
278  {
279  const tetIndices tetIs = this->currentTetIndices();
280  const scalar Gc = td.GInterp().interpolate(this->coordinates(), tetIs);
281  const scalar sigma = physicoChemical::sigma.value();
282  const scalar epsilon = cloud.constProps().epsilon0();
283 
284  ancp += As*epsilon*(Gc/4.0 - sigma*pow4(T_));
285  }
286  ancp /= m*Cp_;
287 
288  // Integrate to find the new parcel temperature
289  const scalar deltaT = cloud.TIntegrator().delta(T_, dt, acp + ancp, bcp);
290  const scalar deltaTncp = ancp*dt;
291  const scalar deltaTcp = deltaT - deltaTncp;
292 
293  // Calculate the new temperature and the enthalpy transfer terms
294  scalar Tnew = T_ + deltaT;
295  Tnew = min(max(Tnew, cloud.constProps().TMin()), cloud.constProps().TMax());
296 
297  dhsTrans -= m*Cp_*deltaTcp;
298 
299  Sph = dt*m*Cp_*bcp;
300 
301  return Tnew;
302 }
303 
304 
305 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
306 
307 template<class ParcelType>
309 (
310  const ThermoParcel<ParcelType>& p
311 )
312 :
313  ParcelType(p),
314  T_(p.T_),
315  Cp_(p.Cp_)
316 {}
317 
318 
319 template<class ParcelType>
321 (
322  const ThermoParcel<ParcelType>& p,
323  const polyMesh& mesh
324 )
325 :
326  ParcelType(p, mesh),
327  T_(p.T_),
328  Cp_(p.Cp_)
329 {}
330 
331 
332 // * * * * * * * * * * * * * * IOStream operators * * * * * * * * * * * * * //
333 
334 #include "ThermoParcelIO.C"
335 
336 // ************************************************************************* //
Collection of constants.
const interpolation< scalar > & kappaInterp() const
Return const access to the interpolator for continuous.
zeroField Su
Definition: alphaSuSp.H:1
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
scalar Cp_
Specific heat capacity [J/kg/K].
Definition: ThermoParcel.H:246
const dimensionedScalar & sigma
Stefan-Boltzmann constant: default SI units: [W/m^2/K^4].
void calc(TrackCloudType &cloud, trackingData &td, const scalar dt)
Update parcel properties over the time interval.
Definition: ThermoParcel.C:137
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
ThermoParcel(const polyMesh &mesh, const barycentric &coordinates, const label celli, const label tetFacei, const label tetPti)
Construct from mesh, coordinates and topology.
Definition: ThermoParcelI.H:75
const interpolation< scalar > & GInterp() const
Return const access to the interpolator for continuous.
const interpolation< scalar > & TInterp() const
Return const access to the interpolator for continuous.
virtual Type interpolate(const vector &position, const label celli, const label facei=-1) const =0
Interpolate field to the given point in the given cell.
const GeometricField< Type, fvPatchField, volMesh > & psi() const
Return the field to be interpolated.
scalar T_
Temperature [K].
Definition: ThermoParcel.H:243
void calcSurfaceValues(TrackCloudType &cloud, trackingData &td, const scalar T, scalar &Ts, scalar &rhos, scalar &mus, scalar &Pr, scalar &kappas) const
Calculate surface thermo properties.
Definition: ThermoParcel.C:94
void cellValueSourceCorrection(TrackCloudType &cloud, trackingData &td, const scalar dt)
Correct cell values using latest transfer information.
Definition: ThermoParcel.C:66
const Type & value() const
Return const reference to value.
static const zero Zero
Definition: zero.H:97
Storage and named access for the indices of a tet which is part of the decomposition of a cell...
Definition: tetIndices.H:81
static const char nl
Definition: Ostream.H:260
scalar Tc() const
Return the continuous phase temperature.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
const interpolation< scalar > & CpInterp() const
Return const access to the interpolator for continuous.
PtrList< coordinateSystem > coordinates(solidRegions.size())
scalar calcHeatTransfer(TrackCloudType &cloud, trackingData &td, const scalar dt, const scalar Re, const scalar Pr, const scalar kappa, const scalar NCpW, const scalar Sh, scalar &dhsTrans, scalar &Sph)
Calculate new particle temperature.
#define WarningInFunction
Report a warning using Foam::Warning.
scalar Cpc() const
Return the continuous phase specific heat capacity.
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:56
scalar epsilon
dimensionedScalar pow4(const dimensionedScalar &ds)
void setCellValues(TrackCloudType &cloud, trackingData &td)
Set cell values.
Definition: ThermoParcel.C:36
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
scalarField Re(const UList< complex > &cf)
Definition: complexFields.C:97
Thermodynamic parcel class with one/two-way coupling with the continuous phase. Includes Kinematic pa...
Definition: ThermoParcel.H:51
scalar T0
Definition: createFields.H:22