SprayParcel.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-2016 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 "SprayParcel.H"
27 #include "CompositionModel.H"
28 #include "AtomizationModel.H"
29 
30 // * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * * //
31 
32 template<class ParcelType>
33 template<class TrackData>
35 (
36  TrackData& td,
37  const scalar dt,
38  const label celli
39 )
40 {
41  ParcelType::setCellValues(td, dt, celli);
42 }
43 
44 
45 template<class ParcelType>
46 template<class TrackData>
48 (
49  TrackData& td,
50  const scalar dt,
51  const label celli
52 )
53 {
54  ParcelType::cellValueSourceCorrection(td, dt, celli);
55 }
56 
57 
58 template<class ParcelType>
59 template<class TrackData>
61 (
62  TrackData& td,
63  const scalar dt,
64  const label celli
65 )
66 {
67  typedef typename TrackData::cloudType::reactingCloudType reactingCloudType;
69  td.cloud().composition();
70 
71  // Check if parcel belongs to liquid core
72  if (liquidCore() > 0.5)
73  {
74  // Liquid core parcels should not experience coupled forces
75  td.cloud().forces().setCalcCoupled(false);
76  }
77 
78  // Get old mixture composition
79  scalarField X0(composition.liquids().X(this->Y()));
80 
81  // Check if we have critical or boiling conditions
82  scalar TMax = composition.liquids().Tc(X0);
83  const scalar T0 = this->T();
84  const scalar pc0 = this->pc_;
85  if (composition.liquids().pv(pc0, T0, X0) >= pc0*0.999)
86  {
87  // Set TMax to boiling temperature
88  TMax = composition.liquids().pvInvert(pc0, X0);
89  }
90 
91  // Set the maximum temperature limit
92  td.cloud().constProps().setTMax(TMax);
93 
94  // Store the parcel properties
95  this->Cp() = composition.liquids().Cp(pc0, T0, X0);
96  sigma_ = composition.liquids().sigma(pc0, T0, X0);
97  const scalar rho0 = composition.liquids().rho(pc0, T0, X0);
98  this->rho() = rho0;
99  const scalar mass0 = this->mass();
100  mu_ = composition.liquids().mu(pc0, T0, X0);
101 
102  ParcelType::calc(td, dt, celli);
103 
104  if (td.keepParticle)
105  {
106  // Reduce the stripped parcel mass due to evaporation
107  // assuming the number of particles remains unchanged
108  this->ms() -= this->ms()*(mass0 - this->mass())/mass0;
109 
110  // Update Cp, sigma, density and diameter due to change in temperature
111  // and/or composition
112  scalar T1 = this->T();
113  scalarField X1(composition.liquids().X(this->Y()));
114 
115  this->Cp() = composition.liquids().Cp(this->pc_, T1, X1);
116 
117  sigma_ = composition.liquids().sigma(this->pc_, T1, X1);
118 
119  scalar rho1 = composition.liquids().rho(this->pc_, T1, X1);
120  this->rho() = rho1;
121 
122  mu_ = composition.liquids().mu(this->pc_, T1, X1);
123 
124  scalar d1 = this->d()*cbrt(rho0/rho1);
125  this->d() = d1;
126 
127  if (liquidCore() > 0.5)
128  {
129  calcAtomization(td, dt, celli);
130 
131  // Preserve the total mass/volume by increasing the number of
132  // particles in parcels due to breakup
133  scalar d2 = this->d();
134  this->nParticle() *= pow3(d1/d2);
135  }
136  else
137  {
138  calcBreakup(td, dt, celli);
139  }
140  }
141 
142  // Restore coupled forces
143  td.cloud().forces().setCalcCoupled(true);
144 }
145 
146 
147 template<class ParcelType>
148 template<class TrackData>
150 (
151  TrackData& td,
152  const scalar dt,
153  const label celli
154 )
155 {
156  typedef typename TrackData::cloudType::reactingCloudType reactingCloudType;
158  td.cloud().composition();
159 
160  typedef typename TrackData::cloudType::sprayCloudType sprayCloudType;
161  const AtomizationModel<sprayCloudType>& atomization =
162  td.cloud().atomization();
163 
164  // Average molecular weight of carrier mix - assumes perfect gas
165  scalar Wc = this->rhoc_*RR*this->Tc()/this->pc();
166  scalar R = RR/Wc;
167  scalar Tav = atomization.Taverage(this->T(), this->Tc());
168 
169  // Calculate average gas density based on average temperature
170  scalar rhoAv = this->pc()/(R*Tav);
171 
172  scalar soi = td.cloud().injectors().timeStart();
173  scalar currentTime = td.cloud().db().time().value();
174  const vector& pos = this->position();
175  const vector& injectionPos = this->position0();
176 
177  // Disregard the continous phase when calculating the relative velocity
178  // (in line with the deactivated coupled assumption)
179  scalar Urel = mag(this->U());
180 
181  scalar t0 = max(0.0, currentTime - this->age() - soi);
182  scalar t1 = min(t0 + dt, td.cloud().injectors().timeEnd() - soi);
183 
184  // This should be the vol flow rate from when the parcel was injected
185  scalar volFlowRate = td.cloud().injectors().volumeToInject(t0, t1)/dt;
186 
187  scalar chi = 0.0;
188  if (atomization.calcChi())
189  {
190  chi = this->chi(td, composition.liquids().X(this->Y()));
191  }
192 
193  atomization.update
194  (
195  dt,
196  this->d(),
197  this->liquidCore(),
198  this->tc(),
199  this->rho(),
200  mu_,
201  sigma_,
202  volFlowRate,
203  rhoAv,
204  Urel,
205  pos,
206  injectionPos,
207  td.cloud().pAmbient(),
208  chi,
209  td.cloud().rndGen()
210  );
211 }
212 
213 
214 template<class ParcelType>
215 template<class TrackData>
217 (
218  TrackData& td,
219  const scalar dt,
220  const label celli
221 )
222 {
223  typedef typename TrackData::cloudType cloudType;
224  typedef typename cloudType::parcelType parcelType;
225  typedef typename cloudType::forceType forceType;
226 
227  const parcelType& p = static_cast<const parcelType&>(*this);
228  const forceType& forces = td.cloud().forces();
229 
230  if (td.cloud().breakup().solveOscillationEq())
231  {
232  solveTABEq(td, dt);
233  }
234 
235  // Average molecular weight of carrier mix - assumes perfect gas
236  scalar Wc = this->rhoc()*RR*this->Tc()/this->pc();
237  scalar R = RR/Wc;
238  scalar Tav = td.cloud().atomization().Taverage(this->T(), this->Tc());
239 
240  // Calculate average gas density based on average temperature
241  scalar rhoAv = this->pc()/(R*Tav);
242  scalar muAv = this->muc();
243  vector Urel = this->U() - this->Uc();
244  scalar Urmag = mag(Urel);
245  scalar Re = this->Re(this->U(), this->d(), rhoAv, muAv);
246 
247  const scalar mass = p.mass();
248  const forceSuSp Fcp = forces.calcCoupled(p, dt, mass, Re, muAv);
249  const forceSuSp Fncp = forces.calcNonCoupled(p, dt, mass, Re, muAv);
250  this->tMom() = mass/(Fcp.Sp() + Fncp.Sp());
251 
252  const vector g = td.cloud().g().value();
253 
254  scalar parcelMassChild = 0.0;
255  scalar dChild = 0.0;
256  if
257  (
258  td.cloud().breakup().update
259  (
260  dt,
261  g,
262  this->d(),
263  this->tc(),
264  this->ms(),
265  this->nParticle(),
266  this->KHindex(),
267  this->y(),
268  this->yDot(),
269  this->d0(),
270  this->rho(),
271  mu_,
272  sigma_,
273  this->U(),
274  rhoAv,
275  muAv,
276  Urel,
277  Urmag,
278  this->tMom(),
279  dChild,
280  parcelMassChild
281  )
282  )
283  {
284  scalar Re = rhoAv*Urmag*dChild/muAv;
285 
286  // Add child parcel as copy of parent
288  child->origId() = this->getNewParticleID();
289  child->d() = dChild;
290  child->d0() = dChild;
291  const scalar massChild = child->mass();
292  child->mass0() = massChild;
293  child->nParticle() = parcelMassChild/massChild;
294 
295  const forceSuSp Fcp =
296  forces.calcCoupled(*child, dt, massChild, Re, muAv);
297  const forceSuSp Fncp =
298  forces.calcNonCoupled(*child, dt, massChild, Re, muAv);
299 
300  child->age() = 0.0;
301  child->liquidCore() = 0.0;
302  child->KHindex() = 1.0;
303  child->y() = td.cloud().breakup().y0();
304  child->yDot() = td.cloud().breakup().yDot0();
305  child->tc() = 0.0;
306  child->ms() = -GREAT;
307  child->injector() = this->injector();
308  child->tMom() = massChild/(Fcp.Sp() + Fncp.Sp());
309  child->user() = 0.0;
310  child->setCellValues(td, dt, celli);
311 
312  td.cloud().addParticle(child);
313  }
314 }
315 
316 
317 template<class ParcelType>
318 template<class TrackData>
320 (
321  TrackData& td,
322  const scalarField& X
323 ) const
324 {
325  // Modifications to take account of the flash boiling on primary break-up
326 
327  typedef typename TrackData::cloudType::reactingCloudType reactingCloudType;
329  td.cloud().composition();
330 
331  scalar chi = 0.0;
332  scalar T0 = this->T();
333  scalar p0 = this->pc();
334  scalar pAmb = td.cloud().pAmbient();
335 
336  scalar pv = composition.liquids().pv(p0, T0, X);
337 
338  forAll(composition.liquids(), i)
339  {
340  if (pv >= 0.999*pAmb)
341  {
342  // Liquid is boiling - calc boiling temperature
343 
344  const liquidProperties& liq = composition.liquids().properties()[i];
345  scalar TBoil = liq.pvInvert(p0);
346 
347  scalar hl = liq.hl(pAmb, TBoil);
348  scalar iTp = liq.h(pAmb, T0) - pAmb/liq.rho(pAmb, T0);
349  scalar iTb = liq.h(pAmb, TBoil) - pAmb/liq.rho(pAmb, TBoil);
350 
351  chi += X[i]*(iTp - iTb)/hl;
352  }
353  }
354 
355  chi = min(1.0, max(chi, 0.0));
356 
357  return chi;
358 }
359 
360 
361 template<class ParcelType>
362 template<class TrackData>
364 (
365  TrackData& td,
366  const scalar dt
367 )
368 {
369  const scalar& TABCmu = td.cloud().breakup().TABCmu();
370  const scalar& TABtwoWeCrit = td.cloud().breakup().TABtwoWeCrit();
371  const scalar& TABComega = td.cloud().breakup().TABComega();
372 
373  scalar r = 0.5*this->d();
374  scalar r2 = r*r;
375  scalar r3 = r*r2;
376 
377  // Inverse of characteristic viscous damping time
378  scalar rtd = 0.5*TABCmu*mu_/(this->rho()*r2);
379 
380  // Oscillation frequency (squared)
381  scalar omega2 = TABComega*sigma_/(this->rho()*r3) - rtd*rtd;
382 
383  if (omega2 > 0)
384  {
385  scalar omega = sqrt(omega2);
386  scalar rhoc = this->rhoc();
387  scalar We = this->We(this->U(), r, rhoc, sigma_)/TABtwoWeCrit;
388 
389  // Initial values for y and yDot
390  scalar y0 = this->y() - We;
391  scalar yDot0 = this->yDot() + y0*rtd;
392 
393  // Update distortion parameters
394  scalar c = cos(omega*dt);
395  scalar s = sin(omega*dt);
396  scalar e = exp(-rtd*dt);
397 
398  this->y() = We + e*(y0*c + (yDot0/omega)*s);
399  this->yDot() = (We - this->y())*rtd + e*(yDot0*c - omega*y0*s);
400  }
401  else
402  {
403  // Reset distortion parameters
404  this->y() = 0;
405  this->yDot() = 0;
406  }
407 }
408 
409 
410 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
411 
412 template<class ParcelType>
414 :
415  ParcelType(p),
416  d0_(p.d0_),
417  position0_(p.position0_),
418  sigma_(p.sigma_),
419  mu_(p.mu_),
420  liquidCore_(p.liquidCore_),
421  KHindex_(p.KHindex_),
422  y_(p.y_),
423  yDot_(p.yDot_),
424  tc_(p.tc_),
425  ms_(p.ms_),
426  injector_(p.injector_),
427  tMom_(p.tMom_),
428  user_(p.user_)
429 {}
430 
431 
432 template<class ParcelType>
434 (
435  const SprayParcel<ParcelType>& p,
436  const polyMesh& mesh
437 )
438 :
439  ParcelType(p, mesh),
440  d0_(p.d0_),
441  position0_(p.position0_),
442  sigma_(p.sigma_),
443  mu_(p.mu_),
444  liquidCore_(p.liquidCore_),
445  KHindex_(p.KHindex_),
446  y_(p.y_),
447  yDot_(p.yDot_),
448  tc_(p.tc_),
449  ms_(p.ms_),
450  injector_(p.injector_),
451  tMom_(p.tMom_),
452  user_(p.user_)
453 {}
454 
455 
456 // * * * * * * * * * * * * * * IOStream operators * * * * * * * * * * * * * //
457 
458 #include "SprayParcelIO.C"
459 
460 
461 // ************************************************************************* //
void solveTABEq(TrackData &td, const scalar dt)
Solve the TAB equation.
Definition: SprayParcel.C:364
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
U
Definition: pEqn.H:83
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
virtual scalar rho(scalar p, scalar T) const =0
Liquid density [kg/m^3].
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
scalar tMom_
Momentum relaxation time (needed for calculating parcel acc.)
Definition: SprayParcel.H:170
scalar KHindex() const
Return const access to Kelvin-Helmholtz breakup index.
Definition: SprayParcelI.H:277
scalar tMom() const
Return const access to momentum relaxation time.
Definition: SprayParcelI.H:319
virtual bool calcChi() const =0
Flag to indicate if chi needs to be calculated.
basicMultiComponentMixture & composition
virtual scalar pvInvert(scalar p) const
Invert the vapour pressure relationship to retrieve the.
scalar liquidCore() const
Return const access to liquid core.
Definition: SprayParcelI.H:270
dimensionedScalar sqrt(const dimensionedScalar &ds)
Templated atomization model class.
Definition: SprayCloud.H:47
void calc(TrackData &td, const scalar dt, const label celli)
Update parcel properties over the time interval.
Definition: SprayParcel.C:61
virtual scalar hl(scalar p, scalar T) const =0
Heat of vapourisation [J/kg].
scalar rho0
dimensionedScalar y0(const dimensionedScalar &ds)
scalar chi(TrackData &td, const scalarField &X) const
Calculate the chi-factor for flash-boiling for the.
scalar ms() const
Return const access to stripped parcel mass.
Definition: SprayParcelI.H:305
void calcAtomization(TrackData &td, const scalar dt, const label celli)
Correct parcel properties according to atomization model.
Definition: SprayParcel.C:150
scalar liquidCore_
Part of liquid core ( >0.5=liquid, <0.5=droplet )
Definition: SprayParcel.H:148
SprayParcel(const polyMesh &mesh, const barycentric &coordinates, const label celli, const label tetFacei, const label tetPti)
Construct from mesh, coordinates and topology.
Definition: SprayParcelI.H:108
scalar y() const
Return const access to spherical deviation.
Definition: SprayParcelI.H:284
Helper container for force Su and Sp terms.
Definition: forceSuSp.H:61
scalarList X0(nSpecie, 0.0)
dimensionedScalar pos(const dimensionedScalar &ds)
scalar y
scalar d0() const
Return const access to initial droplet diameter.
Definition: SprayParcelI.H:242
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 dimensionedScalar e
Elementary charge.
Definition: doubleFloat.H:78
dimensionedScalar cos(const dimensionedScalar &ds)
dimensionedScalar exp(const dimensionedScalar &ds)
dimensionedScalar cbrt(const dimensionedScalar &ds)
scalar tc() const
Return const access to atomization characteristic time.
Definition: SprayParcelI.H:298
virtual scalar h(scalar p, scalar T) const =0
Liquid enthalpy [J/kg] - reference to 298.15 K.
The thermophysical properties of a liquid.
scalar injector_
Injected from injector (needed e.g. for calculating distance.
Definition: SprayParcel.H:167
scalar user_
Passive scalar (extra variable to be defined by user)
Definition: SprayParcel.H:173
Urel
Definition: pEqn.H:56
void cellValueSourceCorrection(TrackData &td, const scalar dt, const label celli)
Correct cell values using latest transfer information.
Definition: SprayParcel.C:48
scalar Sp() const
Return const access to the implicit coefficient [kg/s].
Definition: forceSuSpI.H:62
scalar yDot() const
Return const access to rate of change of spherical deviation.
Definition: SprayParcelI.H:291
dimensionedScalar sin(const dimensionedScalar &ds)
const dimensionedVector & g
scalar user() const
Return const access to passive user scalar.
Definition: SprayParcelI.H:326
scalar sigma_
Liquid surface tension [N/m].
Definition: SprayParcel.H:142
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
const volScalarField & T
dimensionedScalar pow3(const dimensionedScalar &ds)
vector position0_
Injection position.
Definition: SprayParcel.H:139
void calcBreakup(TrackData &td, const scalar dt, const label celli)
Correct parcel properties according to breakup model.
Definition: SprayParcel.C:217
#define R(A, B, C, D, E, F, K, M)
scalar injector() const
Return const access to injector id.
Definition: SprayParcelI.H:312
scalar yDot_
Rate of change of spherical deviation.
Definition: SprayParcel.H:157
PtrList< volScalarField > & Y
const dimensionedScalar c
Speed of light in a vacuum.
virtual void update(const scalar dt, scalar &d, scalar &liquidCore, scalar &tc, const scalar rho, const scalar mu, const scalar sigma, const scalar volFlowRate, const scalar rhoAv, const scalar Urel, const vector &pos, const vector &injectionPos, const scalar pAmbient, const scalar chi, cachedRandom &rndGen) const =0
scalar Taverage(const scalar &Tliq, const scalar &Tc) const
Average temperature calculation.
scalar ms_
Stripped parcel mass due to breakup.
Definition: SprayParcel.H:163
dimensioned< scalar > mag(const dimensioned< Type > &)
scalar d0_
Initial droplet diameter.
Definition: SprayParcel.H:136
const scalar RR
Universal gas constant (default in [J/(kmol K)])
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
scalar tc_
Characteristic time (used in atomization and/or breakup model)
Definition: SprayParcel.H:160
scalar mu_
Liquid dynamic viscosity [Pa.s].
Definition: SprayParcel.H:145
const dimensionedScalar & rho1
Definition: createFields.H:39
Templated reacting parcel composition model class Consists of carrier species (via thermo package)...
Definition: ReactingCloud.H:52
scalar y_
Spherical deviation.
Definition: SprayParcel.H:154
scalar KHindex_
Index for KH Breakup.
Definition: SprayParcel.H:151
scalarField Re(const UList< complex > &cf)
Definition: complexFields.C:97
void setCellValues(TrackData &td, const scalar dt, const label celli)
Set cell values.
Definition: SprayParcel.C:35
Reacing spray parcel, with added functionality for atomization and breakup.
Definition: SprayParcel.H:43
const liquidMixtureProperties & liquids() const
Return the global (additional) liquids.