ThermoParcel.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2011-2017 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 TrackData>
36 (
37  TrackData& td,
38  const scalar dt,
39  const label celli
40 )
41 {
42  ParcelType::setCellValues(td, dt, celli);
43 
44  tetIndices tetIs = this->currentTetIndices();
45 
46  Cpc_ = td.CpInterp().interpolate(this->coordinates(), tetIs);
47 
48  Tc_ = td.TInterp().interpolate(this->coordinates(), tetIs);
49 
50  if (Tc_ < td.cloud().constProps().TMin())
51  {
52  if (debug)
53  {
55  << "Limiting observed temperature in cell " << celli << " to "
56  << td.cloud().constProps().TMin() << nl << endl;
57  }
58 
59  Tc_ = td.cloud().constProps().TMin();
60  }
61 }
62 
63 
64 template<class ParcelType>
65 template<class TrackData>
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  const scalar CpMean = td.CpInterp().psi()[celli];
76  Tc_ += td.cloud().hsTrans()[celli]/(CpMean*this->massCell(celli));
77 
78  if (Tc_ < td.cloud().constProps().TMin())
79  {
80  if (debug)
81  {
83  << "Limiting observed temperature in cell " << celli << " to "
84  << td.cloud().constProps().TMin() << nl << endl;
85  }
86 
87  Tc_ = td.cloud().constProps().TMin();
88  }
89 }
90 
91 
92 template<class ParcelType>
93 template<class TrackData>
95 (
96  TrackData& td,
97  const label celli,
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 + Tc_)/3.0;
108 
109  if (Ts < td.cloud().constProps().TMin())
110  {
111  if (debug)
112  {
114  << "Limiting parcel surface temperature to "
115  << td.cloud().constProps().TMin() << nl << endl;
116  }
117 
118  Ts = td.cloud().constProps().TMin();
119  }
120 
121  // Assuming thermo props vary linearly with T for small d(T)
122  const scalar TRatio = Tc_/Ts;
123 
124  rhos = this->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 = Cpc_*mus/kappas;
131  Pr = max(ROOTVSMALL, Pr);
132 }
133 
134 
135 template<class ParcelType>
136 template<class TrackData>
138 (
139  TrackData& td,
140  const scalar dt,
141  const label celli
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(td, celli, this->T_, Ts, rhos, mus, Pr, kappas);
157 
158  // Reynolds number
159  scalar Re = this->Re(this->U_, this->d_, rhos, 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  td,
195  dt,
196  celli,
197  Re,
198  Pr,
199  kappas,
200  NCpW,
201  Sh,
202  dhsTrans,
203  Sph
204  );
205 
206 
207  // Motion
208  // ~~~~~~
209 
210  // Calculate new particle velocity
211  this->U_ =
212  this->calcVelocity(td, dt, celli, Re, mus, mass0, Su, dUTrans, Spu);
213 
214 
215  // Accumulate carrier phase source terms
216  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
217  if (td.cloud().solution().coupled())
218  {
219  // Update momentum transfer
220  td.cloud().UTrans()[celli] += np0*dUTrans;
221 
222  // Update momentum transfer coefficient
223  td.cloud().UCoeff()[celli] += np0*Spu;
224 
225  // Update sensible enthalpy transfer
226  td.cloud().hsTrans()[celli] += np0*dhsTrans;
227 
228  // Update sensible enthalpy coefficient
229  td.cloud().hsCoeff()[celli] += np0*Sph;
230 
231  // Update radiation fields
232  if (td.cloud().radiation())
233  {
234  const scalar ap = this->areaP();
235  const scalar T4 = pow4(T0);
236  td.cloud().radAreaP()[celli] += dt*np0*ap;
237  td.cloud().radT4()[celli] += dt*np0*T4;
238  td.cloud().radAreaPT4()[celli] += dt*np0*ap*T4;
239  }
240  }
241 }
242 
243 
244 template<class ParcelType>
245 template<class TrackData>
247 (
248  TrackData& td,
249  const scalar dt,
250  const label celli,
251  const scalar Re,
252  const scalar Pr,
253  const scalar kappa,
254  const scalar NCpW,
255  const scalar Sh,
256  scalar& dhsTrans,
257  scalar& Sph
258 )
259 {
260  if (!td.cloud().heatTransfer().active())
261  {
262  return T_;
263  }
264 
265  const scalar d = this->d();
266  const scalar rho = this->rho();
267 
268  // Calc heat transfer coefficient
269  scalar htc = td.cloud().heatTransfer().htc(d, Re, Pr, kappa, NCpW);
270 
271  if (mag(htc) < ROOTVSMALL && !td.cloud().radiation())
272  {
273  return
274  max
275  (
276  T_ + dt*Sh/(this->volume(d)*rho*Cp_),
277  td.cloud().constProps().TMin()
278  );
279  }
280 
281  htc = max(htc, ROOTVSMALL);
282  const scalar As = this->areaS(d);
283 
284  scalar ap = Tc_ + Sh/(As*htc);
285  scalar bp = 6.0*(Sh/As + htc*(Tc_ - T_));
286  if (td.cloud().radiation())
287  {
288  tetIndices tetIs = this->currentTetIndices();
289  const scalar Gc = td.GInterp().interpolate(this->coordinates(), tetIs);
290  const scalar sigma = physicoChemical::sigma.value();
291  const scalar epsilon = td.cloud().constProps().epsilon0();
292 
293  // Assume constant source
294  scalar s = epsilon*(Gc/4.0 - sigma*pow4(T_));
295 
296  ap += s/htc;
297  bp += 6.0*s;
298  }
299  bp /= rho*d*Cp_*(ap - T_) + ROOTVSMALL;
300 
301  // Integrate to find the new parcel temperature
303  td.cloud().TIntegrator().integrate(T_, dt, ap*bp, bp);
304 
305  scalar Tnew =
306  min
307  (
308  max
309  (
310  Tres.value(),
311  td.cloud().constProps().TMin()
312  ),
313  td.cloud().constProps().TMax()
314  );
315 
316  Sph = dt*htc*As;
317 
318  dhsTrans += Sph*(Tres.average() - Tc_);
319 
320  return Tnew;
321 }
322 
323 
324 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
325 
326 template<class ParcelType>
328 (
329  const ThermoParcel<ParcelType>& p
330 )
331 :
332  ParcelType(p),
333  T_(p.T_),
334  Cp_(p.Cp_),
335  Tc_(p.Tc_),
336  Cpc_(p.Cpc_)
337 {}
338 
339 
340 template<class ParcelType>
342 (
343  const ThermoParcel<ParcelType>& p,
344  const polyMesh& mesh
345 )
346 :
347  ParcelType(p, mesh),
348  T_(p.T_),
349  Cp_(p.Cp_),
350  Tc_(p.Tc_),
351  Cpc_(p.Cpc_)
352 {}
353 
354 
355 // * * * * * * * * * * * * * * IOStream operators * * * * * * * * * * * * * //
356 
357 #include "ThermoParcelIO.C"
358 
359 // ************************************************************************* //
Collection of constants.
dimensionedScalar Pr("Pr", dimless, laminarTransport)
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
void cellValueSourceCorrection(TrackData &td, const scalar dt, const label celli)
Correct cell values using latest transfer information.
Definition: ThermoParcel.C:67
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:233
scalar calcHeatTransfer(TrackData &td, const scalar dt, const label celli, const scalar Re, const scalar Pr, const scalar kappa, const scalar NCpW, const scalar Sh, scalar &dhsTrans, scalar &Sph)
Calculate new particle temperature.
void calc(TrackData &td, const scalar dt, const label celli)
Update parcel properties over the time interval.
Definition: ThermoParcel.C:138
Helper class to supply results of integration.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
const dimensionedScalar sigma
Stefan-Boltzmann constant: default SI units: [W/m2/K4].
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
void setCellValues(TrackData &td, const scalar dt, const label celli)
Set cell values.
Definition: ThermoParcel.C:36
scalar T_
Temperature [K].
Definition: ThermoParcel.H:230
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
const Type & value() const
Return const reference to value.
static const zero Zero
Definition: zero.H:91
Storage and named access for the indices of a tet which is part of the decomposition of a cell...
Definition: tetIndices.H:81
scalar Cpc_
Specific heat capacity [J/(kg.K)].
Definition: ThermoParcel.H:242
void calcSurfaceValues(TrackData &td, const label celli, const scalar T, scalar &Ts, scalar &rhos, scalar &mus, scalar &Pr, scalar &kappas) const
Calculate surface thermo properties.
Definition: ThermoParcel.C:95
static const char nl
Definition: Ostream.H:262
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
PtrList< coordinateSystem > coordinates(solidRegions.size())
#define WarningInFunction
Report a warning using Foam::Warning.
scalar epsilon
dimensionedScalar pow4(const dimensionedScalar &ds)
dimensioned< scalar > mag(const dimensioned< Type > &)
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 Tc_
Temperature [K].
Definition: ThermoParcel.H:239