FreeFOAM The Cross-Platform CFD Toolkit
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) 1991-2010 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
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 
28 
29 // * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * * //
30 
31 template<class ParcelType>
32 template<class TrackData>
34 (
35  TrackData& td,
36  const scalar dt,
37  const label cellI
38 )
39 {
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>
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>
83 (
84  TrackData& td,
85  const label cellI,
86  const scalar T,
87  scalar& Ts,
88  scalar& rhos,
89  scalar& mus,
90  scalar& Pr,
91  scalar& kappa
92 ) const
93 {
94  // Surface temperature using two thirds rule
95  Ts = (2.0*T + Tc_)/3.0;
96 
97  // Assuming thermo props vary linearly with T for small dT
98  scalar factor = td.TInterp().interpolate(this->position(), cellI)/Ts;
99  rhos = this->rhoc_*factor;
100  mus = td.muInterp().interpolate(this->position(), cellI)/factor;
101 
102  Pr = td.constProps().Pr();
103  kappa = cpc_*mus/Pr;
104 }
105 
106 
107 template<class ParcelType>
108 template<class TrackData>
110 (
111  TrackData& td,
112  const scalar dt,
113  const label cellI
114 )
115 {
116  // Define local properties at beginning of time step
117  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
118  const scalar np0 = this->nParticle_;
119  const scalar d0 = this->d_;
120  const vector U0 = this->U_;
121  const scalar rho0 = this->rho_;
122  const scalar T0 = this->T_;
123  const scalar cp0 = this->cp_;
124  const scalar mass0 = this->mass();
125 
126 
127  // Calc surface values
128  // ~~~~~~~~~~~~~~~~~~~
129  scalar Ts, rhos, mus, Pr, kappa;
130  calcSurfaceValues(td, cellI, T0, Ts, rhos, mus, Pr, kappa);
131 
132  // Reynolds number
133  scalar Re = this->Re(U0, d0, rhos, mus);
134 
135 
136  // Sources
137  // ~~~~~~~
138 
139  // Explicit momentum source for particle
140  vector Su = vector::zero;
141 
142  // Momentum transfer from the particle to the carrier phase
143  vector dUTrans = vector::zero;
144 
145  // Explicit enthalpy source for particle
146  scalar Sh = 0.0;
147 
148  // Sensible enthalpy transfer from the particle to the carrier phase
149  scalar dhsTrans = 0.0;
150 
151 
152  // Heat transfer
153  // ~~~~~~~~~~~~~
154 
155  // Sum Ni*Cpi*Wi of emission species
156  scalar NCpW = 0.0;
157 
158  // Calculate new particle velocity
159  scalar T1 =
160  calcHeatTransfer
161  (
162  td,
163  dt,
164  cellI,
165  Re,
166  Pr,
167  kappa,
168  d0,
169  rho0,
170  T0,
171  cp0,
172  NCpW,
173  Sh,
174  dhsTrans
175  );
176 
177 
178  // Motion
179  // ~~~~~~
180 
181  // Calculate new particle velocity
182  vector U1 =
183  calcVelocity(td, dt, cellI, Re, mus, d0, U0, rho0, mass0, Su, dUTrans);
184 
185 
186  // Accumulate carrier phase source terms
187  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
188  if (td.cloud().coupled())
189  {
190  // Update momentum transfer
191  td.cloud().UTrans()[cellI] += np0*dUTrans;
192 
193  // Update sensible enthalpy transfer
194  td.cloud().hsTrans()[cellI] += np0*dhsTrans;
195  }
196 
197  // Set new particle properties
198  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
199  this->U_ = U1;
200  T_ = T1;
201 }
202 
203 
204 template<class ParcelType>
205 template <class TrackData>
207 (
208  TrackData& td,
209  const scalar dt,
210  const label cellI,
211  const scalar Re,
212  const scalar Pr,
213  const scalar kappa,
214  const scalar d,
215  const scalar rho,
216  const scalar T,
217  const scalar cp,
218  const scalar NCpW,
219  const scalar Sh,
220  scalar& dhsTrans
221 )
222 {
223  if (!td.cloud().heatTransfer().active())
224  {
225  return T;
226  }
227 
228  // Calc heat transfer coefficient
229  scalar htc = td.cloud().heatTransfer().htc(d, Re, Pr, kappa, NCpW);
230 
231  if (mag(htc) < ROOTVSMALL && !td.cloud().radiation())
232  {
233  return max(T + dt*Sh/(this->volume(d)*rho*cp), td.constProps().TMin());
234  }
235 
236  const scalar As = this->areaS(d);
237  scalar ap = Tc_ + Sh/As/htc;
238  scalar bp = 6.0*(Sh/As + htc*(Tc_ - T));
239  if (td.cloud().radiation())
240  {
241  const scalarField& G =
242  td.cloud().mesh().objectRegistry::lookupObject<volScalarField>("G");
243  const scalar Gc = G[cellI];
244  const scalar sigma = radiation::sigmaSB.value();
245  const scalar epsilon = td.constProps().epsilon0();
246 
247  ap = (ap + epsilon*Gc/(4.0*htc))/(1.0 + epsilon*sigma*pow3(T)/htc);
248  bp += 6.0*(epsilon*(Gc/4.0 - sigma*pow4(T)));
249  }
250  bp /= rho*d*cp*(ap - T);
251 
252  // Integrate to find the new parcel temperature
254  td.cloud().TIntegrator().integrate(T, dt, ap, bp);
255 
256  scalar Tnew = max(Tres.value(), td.constProps().TMin());
257 
258  dhsTrans += dt*htc*As*(0.5*(T + Tnew) - Tc_);
259 
260  return Tnew;
261 }
262 
263 
264 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
265 
266 template <class ParcelType>
268 (
270 )
271 :
272  KinematicParcel<ParcelType>(p),
273  T_(p.T_),
274  cp_(p.cp_),
275  Tc_(p.Tc_),
276  cpc_(p.cpc_)
277 {}
278 
279 
280 // * * * * * * * * * * * * * * IOStream operators * * * * * * * * * * * * * //
281 
282 #include "ThermoParcelIO.C"
283 
284 // ************************ vim: set sw=4 sts=4 et: ************************ //