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-2015 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->position(), tetIs);
47 
48  Tc_ = td.TInterp().interpolate(this->position(), tetIs);
49 
50  if (Tc_ < td.cloud().constProps().TMin())
51  {
52  if (debug)
53  {
54  WarningIn
55  (
56  "void Foam::ThermoParcel<ParcelType>::setCellValues"
57  "("
58  "TrackData&, "
59  "const scalar, "
60  "const label"
61  ")"
62  ) << "Limiting observed temperature in cell " << cellI << " to "
63  << td.cloud().constProps().TMin() << nl << endl;
64  }
65 
66  Tc_ = td.cloud().constProps().TMin();
67  }
68 }
69 
70 
71 template<class ParcelType>
72 template<class TrackData>
74 (
75  TrackData& td,
76  const scalar dt,
77  const label cellI
78 )
79 {
80  this->Uc_ += td.cloud().UTrans()[cellI]/this->massCell(cellI);
81 
82  const scalar CpMean = td.CpInterp().psi()[cellI];
83  Tc_ += td.cloud().hsTrans()[cellI]/(CpMean*this->massCell(cellI));
84 
85  if (Tc_ < td.cloud().constProps().TMin())
86  {
87  if (debug)
88  {
89  WarningIn
90  (
91  "void Foam::ThermoParcel<ParcelType>::cellValueSourceCorrection"
92  "("
93  "TrackData&, "
94  "const scalar, "
95  "const label"
96  ")"
97  ) << "Limiting observed temperature in cell " << cellI << " to "
98  << td.cloud().constProps().TMin() << nl << endl;
99  }
100 
101  Tc_ = td.cloud().constProps().TMin();
102  }
103 }
104 
105 
106 template<class ParcelType>
107 template<class TrackData>
109 (
110  TrackData& td,
111  const label cellI,
112  const scalar T,
113  scalar& Ts,
114  scalar& rhos,
115  scalar& mus,
116  scalar& Pr,
117  scalar& kappas
118 ) const
119 {
120  // Surface temperature using two thirds rule
121  Ts = (2.0*T + Tc_)/3.0;
122 
123  if (Ts < td.cloud().constProps().TMin())
124  {
125  if (debug)
126  {
127  WarningIn
128  (
129  "void Foam::ThermoParcel<ParcelType>::calcSurfaceValues"
130  "("
131  "TrackData&, "
132  "const label, "
133  "const scalar, "
134  "scalar&, "
135  "scalar&, "
136  "scalar&, "
137  "scalar&, "
138  "scalar&"
139  ") const"
140  ) << "Limiting parcel surface temperature to "
141  << td.cloud().constProps().TMin() << nl << endl;
142  }
143 
144  Ts = td.cloud().constProps().TMin();
145  }
146 
147  // Assuming thermo props vary linearly with T for small d(T)
148  const scalar TRatio = Tc_/Ts;
149 
150  rhos = this->rhoc_*TRatio;
151 
152  tetIndices tetIs = this->currentTetIndices();
153  mus = td.muInterp().interpolate(this->position(), tetIs)/TRatio;
154  kappas = td.kappaInterp().interpolate(this->position(), tetIs)/TRatio;
155 
156  Pr = Cpc_*mus/kappas;
157  Pr = max(ROOTVSMALL, Pr);
158 }
159 
160 
161 template<class ParcelType>
162 template<class TrackData>
164 (
165  TrackData& td,
166  const scalar dt,
167  const label cellI
168 )
169 {
170  // Define local properties at beginning of time step
171  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
172  const scalar np0 = this->nParticle_;
173  const scalar mass0 = this->mass();
174 
175  // Store T for consistent radiation source
176  const scalar T0 = this->T_;
177 
178 
179  // Calc surface values
180  // ~~~~~~~~~~~~~~~~~~~
181  scalar Ts, rhos, mus, Pr, kappas;
182  calcSurfaceValues(td, cellI, this->T_, Ts, rhos, mus, Pr, kappas);
183 
184  // Reynolds number
185  scalar Re = this->Re(this->U_, this->d_, rhos, mus);
186 
187 
188  // Sources
189  // ~~~~~~~
190 
191  // Explicit momentum source for particle
193 
194  // Linearised momentum source coefficient
195  scalar Spu = 0.0;
196 
197  // Momentum transfer from the particle to the carrier phase
198  vector dUTrans = vector::zero;
199 
200  // Explicit enthalpy source for particle
201  scalar Sh = 0.0;
202 
203  // Linearised enthalpy source coefficient
204  scalar Sph = 0.0;
205 
206  // Sensible enthalpy transfer from the particle to the carrier phase
207  scalar dhsTrans = 0.0;
208 
209 
210  // Heat transfer
211  // ~~~~~~~~~~~~~
212 
213  // Sum Ni*Cpi*Wi of emission species
214  scalar NCpW = 0.0;
215 
216  // Calculate new particle temperature
217  this->T_ =
218  this->calcHeatTransfer
219  (
220  td,
221  dt,
222  cellI,
223  Re,
224  Pr,
225  kappas,
226  NCpW,
227  Sh,
228  dhsTrans,
229  Sph
230  );
231 
232 
233  // Motion
234  // ~~~~~~
235 
236  // Calculate new particle velocity
237  this->U_ =
238  this->calcVelocity(td, dt, cellI, Re, mus, mass0, Su, dUTrans, Spu);
239 
240 
241  // Accumulate carrier phase source terms
242  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
243  if (td.cloud().solution().coupled())
244  {
245  // Update momentum transfer
246  td.cloud().UTrans()[cellI] += np0*dUTrans;
247 
248  // Update momentum transfer coefficient
249  td.cloud().UCoeff()[cellI] += np0*Spu;
250 
251  // Update sensible enthalpy transfer
252  td.cloud().hsTrans()[cellI] += np0*dhsTrans;
253 
254  // Update sensible enthalpy coefficient
255  td.cloud().hsCoeff()[cellI] += np0*Sph;
256 
257  // Update radiation fields
258  if (td.cloud().radiation())
259  {
260  const scalar ap = this->areaP();
261  const scalar T4 = pow4(T0);
262  td.cloud().radAreaP()[cellI] += dt*np0*ap;
263  td.cloud().radT4()[cellI] += dt*np0*T4;
264  td.cloud().radAreaPT4()[cellI] += dt*np0*ap*T4;
265  }
266  }
267 }
268 
269 
270 template<class ParcelType>
271 template<class TrackData>
273 (
274  TrackData& td,
275  const scalar dt,
276  const label cellI,
277  const scalar Re,
278  const scalar Pr,
279  const scalar kappa,
280  const scalar NCpW,
281  const scalar Sh,
282  scalar& dhsTrans,
283  scalar& Sph
284 )
285 {
286  if (!td.cloud().heatTransfer().active())
287  {
288  return T_;
289  }
290 
291  const scalar d = this->d();
292  const scalar rho = this->rho();
293 
294  // Calc heat transfer coefficient
295  scalar htc = td.cloud().heatTransfer().htc(d, Re, Pr, kappa, NCpW);
296 
297  if (mag(htc) < ROOTVSMALL && !td.cloud().radiation())
298  {
299  return
300  max
301  (
302  T_ + dt*Sh/(this->volume(d)*rho*Cp_),
303  td.cloud().constProps().TMin()
304  );
305  }
306 
307  htc = max(htc, ROOTVSMALL);
308  const scalar As = this->areaS(d);
309 
310  scalar ap = Tc_ + Sh/(As*htc);
311  scalar bp = 6.0*(Sh/As + htc*(Tc_ - T_));
312  if (td.cloud().radiation())
313  {
314  tetIndices tetIs = this->currentTetIndices();
315  const scalar Gc = td.GInterp().interpolate(this->position(), tetIs);
316  const scalar sigma = physicoChemical::sigma.value();
317  const scalar epsilon = td.cloud().constProps().epsilon0();
318 
319  // Assume constant source
320  scalar s = epsilon*(Gc/4.0 - sigma*pow4(T_));
321 
322  ap += s/htc;
323  bp += 6.0*s;
324  }
325  bp /= rho*d*Cp_*(ap - T_) + ROOTVSMALL;
326 
327  // Integrate to find the new parcel temperature
329  td.cloud().TIntegrator().integrate(T_, dt, ap*bp, bp);
330 
331  scalar Tnew =
332  min
333  (
334  max
335  (
336  Tres.value(),
337  td.cloud().constProps().TMin()
338  ),
339  td.cloud().constProps().TMax()
340  );
341 
342  Sph = dt*htc*As;
343 
344  dhsTrans += Sph*(Tres.average() - Tc_);
345 
346  return Tnew;
347 }
348 
349 
350 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
351 
352 template<class ParcelType>
354 (
355  const ThermoParcel<ParcelType>& p
356 )
357 :
358  ParcelType(p),
359  T_(p.T_),
360  Cp_(p.Cp_),
361  Tc_(p.Tc_),
362  Cpc_(p.Cpc_)
363 {}
364 
365 
366 template<class ParcelType>
368 (
369  const ThermoParcel<ParcelType>& p,
370  const polyMesh& mesh
371 )
372 :
373  ParcelType(p, mesh),
374  T_(p.T_),
375  Cp_(p.Cp_),
376  Tc_(p.Tc_),
377  Cpc_(p.Cpc_)
378 {}
379 
380 
381 // * * * * * * * * * * * * * * IOStream operators * * * * * * * * * * * * * //
382 
383 #include "ThermoParcelIO.C"
384 
385 // ************************************************************************* //
dimensionedScalar Pr("Pr", dimless, laminarTransport)
scalar Tc_
Temperature [K].
Definition: ThermoParcel.H:239
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 ))
dimensioned< scalar > mag(const dimensioned< Type > &)
scalar Cp_
Specific heat capacity [J/(kg.K)].
Definition: ThermoParcel.H:233
Collection of constants.
const dimensionedScalar sigma
Stefan-Boltzmann constant: default SI units: [W/m2/K4].
scalar epsilon
Helper class to supply results of integration.
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 calc(TrackData &td, const scalar dt, const label cellI)
Update parcel properties over the time interval.
Definition: ThermoParcel.C:164
tmp< GeometricField< Type, fvPatchField, volMesh > > Su(const GeometricField< Type, fvPatchField, volMesh > &su, const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: fvcSup.C:44
static const char nl
Definition: Ostream.H:260
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
void cellValueSourceCorrection(TrackData &td, const scalar dt, const label cellI)
Correct cell values using latest transfer information.
Definition: ThermoParcel.C:74
#define WarningIn(functionName)
Report a warning using Foam::Warning.
Thermodynamic parcel class with one/two-way coupling with the continuous phase. Includes Kinematic pa...
Definition: ThermoParcel.H:51
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 setCellValues(TrackData &td, const scalar dt, const label cellI)
Set cell values.
Definition: ThermoParcel.C:36
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
scalar T_
Temperature [K].
Definition: ThermoParcel.H:230
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
static const Vector zero
Definition: Vector.H:80
dimensionedScalar pow4(const dimensionedScalar &ds)
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:109
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
Storage and named access for the indices of a tet which is part of the decomposition of a cell...
Definition: tetIndices.H:73
ThermoParcel(const polyMesh &mesh, const vector &position, const label cellI, const label tetFaceI, const label tetPtI)
Construct from owner, position, and cloud owner.
Definition: ThermoParcelI.H:75
scalar Cpc_
Specific heat capacity [J/(kg.K)].
Definition: ThermoParcel.H:242
const Type & value() const
Return const reference to value.
scalarField Re(const UList< complex > &cf)
Definition: complexFields.C:97