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-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 "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->d() = dChild;
289  child->d0() = dChild;
290  const scalar massChild = child->mass();
291  child->mass0() = massChild;
292  child->nParticle() = parcelMassChild/massChild;
293 
294  const forceSuSp Fcp =
295  forces.calcCoupled(*child, dt, massChild, Re, muAv);
296  const forceSuSp Fncp =
297  forces.calcNonCoupled(*child, dt, massChild, Re, muAv);
298 
299  child->age() = 0.0;
300  child->liquidCore() = 0.0;
301  child->KHindex() = 1.0;
302  child->y() = td.cloud().breakup().y0();
303  child->yDot() = td.cloud().breakup().yDot0();
304  child->tc() = 0.0;
305  child->ms() = -GREAT;
306  child->injector() = this->injector();
307  child->tMom() = massChild/(Fcp.Sp() + Fncp.Sp());
308  child->user() = 0.0;
309  child->setCellValues(td, dt, cellI);
310 
311  td.cloud().addParticle(child);
312  }
313 }
314 
315 
316 template<class ParcelType>
317 template<class TrackData>
319 (
320  TrackData& td,
321  const scalarField& X
322 ) const
323 {
324  // Modifications to take account of the flash boiling on primary break-up
325 
326  typedef typename TrackData::cloudType::reactingCloudType reactingCloudType;
328  td.cloud().composition();
329 
330  scalar chi = 0.0;
331  scalar T0 = this->T();
332  scalar p0 = this->pc();
333  scalar pAmb = td.cloud().pAmbient();
334 
335  scalar pv = composition.liquids().pv(p0, T0, X);
336 
337  forAll(composition.liquids(), i)
338  {
339  if (pv >= 0.999*pAmb)
340  {
341  // Liquid is boiling - calc boiling temperature
342 
343  const liquidProperties& liq = composition.liquids().properties()[i];
344  scalar TBoil = liq.pvInvert(p0);
345 
346  scalar hl = liq.hl(pAmb, TBoil);
347  scalar iTp = liq.h(pAmb, T0) - pAmb/liq.rho(pAmb, T0);
348  scalar iTb = liq.h(pAmb, TBoil) - pAmb/liq.rho(pAmb, TBoil);
349 
350  chi += X[i]*(iTp - iTb)/hl;
351  }
352  }
353 
354  chi = min(1.0, max(chi, 0.0));
355 
356  return chi;
357 }
358 
359 
360 template<class ParcelType>
361 template<class TrackData>
363 (
364  TrackData& td,
365  const scalar dt
366 )
367 {
368  const scalar& TABCmu = td.cloud().breakup().TABCmu();
369  const scalar& TABtwoWeCrit = td.cloud().breakup().TABtwoWeCrit();
370  const scalar& TABComega = td.cloud().breakup().TABComega();
371 
372  scalar r = 0.5*this->d();
373  scalar r2 = r*r;
374  scalar r3 = r*r2;
375 
376  // Inverse of characteristic viscous damping time
377  scalar rtd = 0.5*TABCmu*mu_/(this->rho()*r2);
378 
379  // Oscillation frequency (squared)
380  scalar omega2 = TABComega*sigma_/(this->rho()*r3) - rtd*rtd;
381 
382  if (omega2 > 0)
383  {
384  scalar omega = sqrt(omega2);
385  scalar rhoc = this->rhoc();
386  scalar We = this->We(this->U(), r, rhoc, sigma_)/TABtwoWeCrit;
387 
388  // Initial values for y and yDot
389  scalar y0 = this->y() - We;
390  scalar yDot0 = this->yDot() + y0*rtd;
391 
392  // Update distortion parameters
393  scalar c = cos(omega*dt);
394  scalar s = sin(omega*dt);
395  scalar e = exp(-rtd*dt);
396 
397  this->y() = We + e*(y0*c + (yDot0/omega)*s);
398  this->yDot() = (We - this->y())*rtd + e*(yDot0*c - omega*y0*s);
399  }
400  else
401  {
402  // Reset distortion parameters
403  this->y() = 0;
404  this->yDot() = 0;
405  }
406 }
407 
408 
409 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
410 
411 template<class ParcelType>
413 :
414  ParcelType(p),
415  d0_(p.d0_),
416  position0_(p.position0_),
417  sigma_(p.sigma_),
418  mu_(p.mu_),
419  liquidCore_(p.liquidCore_),
420  KHindex_(p.KHindex_),
421  y_(p.y_),
422  yDot_(p.yDot_),
423  tc_(p.tc_),
424  ms_(p.ms_),
425  injector_(p.injector_),
426  tMom_(p.tMom_),
427  user_(p.user_)
428 {}
429 
430 
431 template<class ParcelType>
433 (
434  const SprayParcel<ParcelType>& p,
435  const polyMesh& mesh
436 )
437 :
438  ParcelType(p, mesh),
439  d0_(p.d0_),
440  position0_(p.position0_),
441  sigma_(p.sigma_),
442  mu_(p.mu_),
443  liquidCore_(p.liquidCore_),
444  KHindex_(p.KHindex_),
445  y_(p.y_),
446  yDot_(p.yDot_),
447  tc_(p.tc_),
448  ms_(p.ms_),
449  injector_(p.injector_),
450  tMom_(p.tMom_),
451  user_(p.user_)
452 {}
453 
454 
455 // * * * * * * * * * * * * * * IOStream operators * * * * * * * * * * * * * //
456 
457 #include "SprayParcelIO.C"
458 
459 
460 // ************************************************************************* //
dimensionedScalar sqrt(const dimensionedScalar &ds)
void calc(const argList &args, const Time &runTime, const fvMesh &mesh)
dimensionedScalar pow3(const dimensionedScalar &ds)
The thermophysical properties of a liquidProperties.
void calcAtomization(TrackData &td, const scalar dt, const label cellI)
Correct parcel properties according to atomization model.
Definition: SprayParcel.C:150
void setCellValues(TrackData &td, const scalar dt, const label cellI)
Set cell values.
Definition: SprayParcel.C:35
scalar yDot() const
Return const access to rate of change of spherical deviation.
Definition: SprayParcelI.H:266
scalar yDot_
Rate of change of spherical deviation.
Definition: SprayParcel.H:157
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 ))
void calcBreakup(TrackData &td, const scalar dt, const label cellI)
Correct parcel properties according to breakup model.
Definition: SprayParcel.C:217
dimensioned< scalar > mag(const dimensioned< Type > &)
void calc(TrackData &td, const scalar dt, const label cellI)
Update parcel properties over the time interval.
Definition: SprayParcel.C:61
scalar d0() const
Return const access to initial droplet diameter.
Definition: SprayParcelI.H:217
This function object calculates the forces and moments by integrating the pressure and skin-friction ...
Definition: forces.H:206
#define R(A, B, C, D, E, F, K, M)
scalarList X0(nSpecie, 0.0)
scalar KHindex_
Index for KH Breakup.
Definition: SprayParcel.H:151
scalar tMom_
Momentum relaxation time (needed for calculating parcel acc.)
Definition: SprayParcel.H:170
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
Templated atomization model class.
Definition: SprayCloud.H:47
scalar Sp() const
Return const access to the implicit coefficient [kg/s].
Definition: forceSuSpI.H:62
scalar mu_
Liquid dynamic viscosity [Pa.s].
Definition: SprayParcel.H:145
Urel
Definition: pEqn.H:53
dimensionedScalar cbrt(const dimensionedScalar &ds)
scalar rho0
dimensionedScalar exp(const dimensionedScalar &ds)
void solveTABEq(TrackData &td, const scalar dt)
Solve the TAB equation.
Definition: SprayParcel.C:363
const liquidMixtureProperties & liquids() const
Return the global (additional) liquids.
dimensionedScalar y0(const dimensionedScalar &ds)
forces(const forces &)
Disallow default bitwise copy construct.
virtual scalar h(scalar p, scalar T) const
Liquid enthalpy [J/kg] - reference to 298.15 K.
const volScalarField & T
Definition: createFields.H:25
scalar tc() const
Return const access to atomization characteristic time.
Definition: SprayParcelI.H:273
scalar sigma_
Liquid surface tension [N/m].
Definition: SprayParcel.H:142
Reacing spray parcel, with added functionality for atomization and breakup.
Definition: SprayParcel.H:43
scalar d0_
Initial droplet diameter.
Definition: SprayParcel.H:136
vector position0_
Injection position.
Definition: SprayParcel.H:139
basicMultiComponentMixture & composition
Definition: createFields.H:35
scalar ms_
Stripped parcel mass due to breakup.
Definition: SprayParcel.H:163
scalar injector_
Injected from injector (needed e.g. for calculating distance.
Definition: SprayParcel.H:167
const dimensionedScalar e
Elementary charge.
Definition: doubleFloat.H:78
scalar chi(TrackData &td, const scalarField &X) const
Calculate the chi-factor for flash-boiling for the.
PtrList< volScalarField > & Y
Definition: createFields.H:36
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
scalar tc_
Characteristic time (used in atomization and/or breakup model)
Definition: SprayParcel.H:160
#define forAll(list, i)
Definition: UList.H:421
rho1
Definition: pEqn.H:124
dimensionedScalar cos(const dimensionedScalar &ds)
void cellValueSourceCorrection(TrackData &td, const scalar dt, const label cellI)
Correct cell values using latest transfer information.
Definition: SprayParcel.C:48
scalar liquidCore() const
Return const access to liquid core.
Definition: SprayParcelI.H:245
scalar user() const
Return const access to passive user scalar.
Definition: SprayParcelI.H:301
const scalar RR
Universal gas constant (default in [J/(kmol K)])
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
const dimensionedVector & g
scalar y
virtual scalar pvInvert(scalar p) const
Invert the vapour pressure relationship to retrieve the.
virtual scalar rho(scalar p, scalar T) const
Liquid rho [kg/m^3].
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 liquidCore_
Part of liquid core ( >0.5=liquid, <0.5=droplet )
Definition: SprayParcel.H:148
Templated reacting parcel composition model class Consists of carrier species (via thermo package)...
Definition: ReactingCloud.H:52
const dimensionedScalar c
Speed of light in a vacuum.
scalar ms() const
Return const access to stripped parcel mass.
Definition: SprayParcelI.H:280
dimensionedScalar pos(const dimensionedScalar &ds)
scalar user_
Passive scalar (extra variable to be defined by user)
Definition: SprayParcel.H:173
SprayParcel(const polyMesh &mesh, const vector &position, const label cellI, const label tetFaceI, const label tetPtI)
Construct from owner, position, and cloud owner.
Definition: SprayParcelI.H:108
scalar injector() const
Return const access to injector id.
Definition: SprayParcelI.H:287
scalar Taverage(const scalar &Tliq, const scalar &Tc) const
Average temperature calculation.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
scalar y_
Spherical deviation.
Definition: SprayParcel.H:154
scalar tMom() const
Return const access to momentum relaxation time.
Definition: SprayParcelI.H:294
scalar KHindex() const
Return const access to Kelvin-Helmholtz breakup index.
Definition: SprayParcelI.H:252
virtual scalar hl(scalar p, scalar T) const
Heat of vapourisation [J/kg].
U
Definition: pEqn.H:82
scalar y() const
Return const access to spherical deviation.
Definition: SprayParcelI.H:259
Helper container for force Su and Sp terms.
Definition: forceSuSp.H:58
scalarField Re(const UList< complex > &cf)
Definition: complexFields.C:97
dimensionedScalar sin(const dimensionedScalar &ds)
virtual bool calcChi() const =0
Flag to indicate if chi needs to be calculated.