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