SprayParcel.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2021 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 "forceSuSp.H"
28 #include "CompositionModel.H"
29 #include "AtomisationModel.H"
30 
31 // * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * * //
32 
33 template<class ParcelType>
34 template<class TrackCloudType>
36 (
37  TrackCloudType& cloud,
38  trackingData& td
39 )
40 {
41  ParcelType::setCellValues(cloud, td);
42 }
43 
44 
45 template<class ParcelType>
46 template<class TrackCloudType>
48 (
49  TrackCloudType& cloud,
50  trackingData& td,
51  const scalar dt
52 )
53 {
54  ParcelType::cellValueSourceCorrection(cloud, td, dt);
55 }
56 
57 
58 template<class ParcelType>
59 template<class TrackCloudType>
61 (
62  TrackCloudType& cloud,
63  trackingData& td,
64  const scalar dt
65 )
66 {
67  typedef typename TrackCloudType::thermoCloudType thermoCloudType;
69  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  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 = td.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  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(cloud,td, dt);
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(td.pc(), T1, X1);
116 
117  sigma_ = composition.liquids().sigma(td.pc(), T1, X1);
118 
119  scalar rho1 = composition.liquids().rho(td.pc(), T1, X1);
120  this->rho() = rho1;
121 
122  mu_ = composition.liquids().mu(td.pc(), T1, X1);
123 
124  scalar d1 = this->d()*cbrt(rho0/rho1);
125  this->d() = d1;
126 
127  if (liquidCore() > 0.5)
128  {
129  calcAtomisation(cloud, td, dt);
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(cloud, td, dt);
139  }
140  }
141 
142  // Restore coupled forces
143  cloud.forces().setCalcCoupled(true);
144 }
145 
146 
147 template<class ParcelType>
148 template<class TrackCloudType>
150 (
151  TrackCloudType& cloud,
152  trackingData& td,
153  const scalar dt
154 )
155 {
156  typedef typename TrackCloudType::thermoCloudType thermoCloudType;
158  cloud.composition();
159 
160  typedef typename TrackCloudType::sprayCloudType sprayCloudType;
161  const AtomisationModel<sprayCloudType>& atomisation =
162  cloud.atomisation();
163 
164  // Average molecular weight of carrier mix - assumes perfect gas
165  scalar Wc = td.rhoc()*RR*td.Tc()/td.pc();
166  scalar R = RR/Wc;
167  scalar Tav = atomisation.Taverage(this->T(), td.Tc());
168 
169  // Calculate average gas density based on average temperature
170  scalar rhoAv = td.pc()/(R*Tav);
171 
172  scalar soi = cloud.injectors().timeStart();
173  scalar currentTime = cloud.db().time().value();
174  const vector& pos = this->position();
175  const vector& injectionPos = this->position0();
176 
177  // Disregard the continuous 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, cloud.injectors().timeEnd() - soi);
183 
184  // This should be the vol flow rate from when the parcel was injected
185  scalar volFlowRate = cloud.injectors().volumeToInject(t0, t1)/dt;
186 
187  scalar chi = 0.0;
188  if (atomisation.calcChi())
189  {
190  chi = this->chi(cloud, td, composition.liquids().X(this->Y()));
191  }
192 
193  atomisation.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  cloud.pAmbient(),
208  chi,
209  cloud.rndGen()
210  );
211 }
212 
213 
214 template<class ParcelType>
215 template<class TrackCloudType>
217 (
218  TrackCloudType& cloud,
219  trackingData& td,
220  const scalar dt
221 )
222 {
223  const typename TrackCloudType::parcelType& p =
224  static_cast<const typename TrackCloudType::parcelType&>(*this);
225  typename TrackCloudType::parcelType::trackingData& ttd =
226  static_cast<typename TrackCloudType::parcelType::trackingData&>(td);
227 
228  const typename TrackCloudType::forceType& forces = cloud.forces();
229 
230  if (cloud.breakup().solveOscillationEq())
231  {
232  solveTABEq(cloud, td, dt);
233  }
234 
235  // Average molecular weight of carrier mix - assumes perfect gas
236  scalar Wc = td.rhoc()*RR*td.Tc()/td.pc();
237  scalar R = RR/Wc;
238  scalar Tav = cloud.atomisation().Taverage(this->T(), td.Tc());
239 
240  // Calculate average gas density based on average temperature
241  scalar rhoAv = td.pc()/(R*Tav);
242  scalar muAv = td.muc();
243  vector Urel = this->U() - td.Uc();
244  scalar Urmag = mag(Urel);
245  scalar Re = this->Re(rhoAv, this->U(), td.Uc(), this->d(), muAv);
246 
247  const scalar mass = p.mass();
248  const forceSuSp Fcp = forces.calcCoupled(p, ttd, dt, mass, Re, muAv);
249  const forceSuSp Fncp = forces.calcNonCoupled(p, ttd, dt, mass, Re, muAv);
250  this->tMom() = mass/(Fcp.Sp() + Fncp.Sp());
251 
252  const vector g = cloud.g().value();
253 
254  scalar parcelMassChild = 0.0;
255  scalar dChild = 0.0;
256  if
257  (
258  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, ttd, dt, massChild, Re, muAv);
297  const forceSuSp Fncp =
298  forces.calcNonCoupled(*child, ttd, dt, massChild, Re, muAv);
299 
300  child->age() = 0.0;
301  child->liquidCore() = 0.0;
302  child->KHindex() = 1.0;
303  child->y() = cloud.breakup().y0();
304  child->yDot() = 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->calcDispersion(cloud, td, dt);
311 
312  cloud.addParticle(child);
313  }
314 }
315 
316 
317 template<class ParcelType>
318 template<class TrackCloudType>
320 (
321  TrackCloudType& cloud,
322  trackingData& td,
323  const scalarField& X
324 ) const
325 {
326  // Modifications to take account of the flash boiling on primary break-up
327 
328  typedef typename TrackCloudType::thermoCloudType thermoCloudType;
330  cloud.composition();
331 
332  scalar chi = 0.0;
333  scalar T0 = this->T();
334  scalar p0 = td.pc();
335  scalar pAmb = cloud.pAmbient();
336 
337  scalar pv = composition.liquids().pv(p0, T0, X);
338 
339  forAll(composition.liquids(), i)
340  {
341  if (pv >= 0.999*pAmb)
342  {
343  // Liquid is boiling - calc boiling temperature
344 
345  const liquidProperties& liq = composition.liquids().properties()[i];
346  scalar TBoil = liq.pvInvert(p0);
347 
348  scalar hl = liq.hl(pAmb, TBoil);
349  scalar iTp = liq.Ha(pAmb, T0) - pAmb/liq.rho(pAmb, T0);
350  scalar iTb = liq.Ha(pAmb, TBoil) - pAmb/liq.rho(pAmb, TBoil);
351 
352  chi += X[i]*(iTp - iTb)/hl;
353  }
354  }
355 
356  chi = min(1.0, max(chi, 0.0));
357 
358  return chi;
359 }
360 
361 
362 template<class ParcelType>
363 template<class TrackCloudType>
365 (
366  TrackCloudType& cloud,
367  trackingData& td,
368  const scalar dt
369 )
370 {
371  const scalar& TABCmu = cloud.breakup().TABCmu();
372  const scalar& TABtwoWeCrit = cloud.breakup().TABtwoWeCrit();
373  const scalar& TABComega = cloud.breakup().TABComega();
374 
375  scalar r = 0.5*this->d();
376  scalar r2 = r*r;
377  scalar r3 = r*r2;
378 
379  // Inverse of characteristic viscous damping time
380  scalar rtd = 0.5*TABCmu*mu_/(this->rho()*r2);
381 
382  // Oscillation frequency (squared)
383  scalar omega2 = TABComega*sigma_/(this->rho()*r3) - rtd*rtd;
384 
385  if (omega2 > 0)
386  {
387  scalar omega = sqrt(omega2);
388  scalar We =
389  this->We(td.rhoc(), this->U(), td.Uc(), r, sigma_)/TABtwoWeCrit;
390 
391  // Initial values for y and yDot
392  scalar y0 = this->y() - We;
393  scalar yDot0 = this->yDot() + y0*rtd;
394 
395  // Update distortion parameters
396  scalar c = cos(omega*dt);
397  scalar s = sin(omega*dt);
398  scalar e = exp(-rtd*dt);
399 
400  this->y() = We + e*(y0*c + (yDot0/omega)*s);
401  this->yDot() = (We - this->y())*rtd + e*(yDot0*c - omega*y0*s);
402  }
403  else
404  {
405  // Reset distortion parameters
406  this->y() = 0;
407  this->yDot() = 0;
408  }
409 }
410 
411 
412 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
413 
414 template<class ParcelType>
416 :
417  ParcelType(p),
418  d0_(p.d0_),
419  position0_(p.position0_),
420  sigma_(p.sigma_),
421  mu_(p.mu_),
422  liquidCore_(p.liquidCore_),
423  KHindex_(p.KHindex_),
424  y_(p.y_),
425  yDot_(p.yDot_),
426  tc_(p.tc_),
427  ms_(p.ms_),
428  injector_(p.injector_),
429  tMom_(p.tMom_),
430  user_(p.user_)
431 {}
432 
433 
434 template<class ParcelType>
436 (
437  const SprayParcel<ParcelType>& p,
438  const polyMesh& mesh
439 )
440 :
441  ParcelType(p, mesh),
442  d0_(p.d0_),
443  position0_(p.position0_),
444  sigma_(p.sigma_),
445  mu_(p.mu_),
446  liquidCore_(p.liquidCore_),
447  KHindex_(p.KHindex_),
448  y_(p.y_),
449  yDot_(p.yDot_),
450  tc_(p.tc_),
451  ms_(p.ms_),
452  injector_(p.injector_),
453  tMom_(p.tMom_),
454  user_(p.user_)
455 {}
456 
457 
458 // * * * * * * * * * * * * * * IOStream operators * * * * * * * * * * * * * //
459 
460 #include "SprayParcelIO.C"
461 
462 
463 // ************************************************************************* //
const volScalarField & rho1
virtual bool calcChi() const =0
Flag to indicate if chi needs to be calculated.
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
virtual scalar rho(scalar p, scalar T) const =0
Liquid density [kg/m^3].
U
Definition: pEqn.H:72
scalar tMom_
Momentum relaxation time (needed for calculating parcel acc.)
Definition: SprayParcel.H:183
basicSpecieMixture & composition
scalar KHindex() const
Return const access to Kelvin-Helmholtz breakup index.
Definition: SprayParcelI.H:214
scalar tMom() const
Return const access to momentum relaxation time.
Definition: SprayParcelI.H:256
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:207
dimensionedScalar sqrt(const dimensionedScalar &ds)
void calcAtomisation(TrackCloudType &cloud, trackingData &td, const scalar dt)
Correct parcel properties according to atomisation model.
Definition: SprayParcel.C:150
virtual scalar hl(scalar p, scalar T) const =0
Heat of vapourisation [J/kg].
Templated atomisation model class.
Definition: SprayCloud.H:52
scalar rho0
dimensionedScalar y0(const dimensionedScalar &ds)
void calcBreakup(TrackCloudType &cloud, trackingData &td, const scalar dt)
Correct parcel properties according to breakup model.
Definition: SprayParcel.C:217
scalar ms() const
Return const access to stripped parcel mass.
Definition: SprayParcelI.H:242
const dimensionedScalar c
Speed of light in a vacuum.
scalar liquidCore_
Part of liquid core ( >0.5=liquid, <0.5=droplet )
Definition: SprayParcel.H:161
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:221
Helper container for force Su and Sp terms.
Definition: forceSuSp.H:61
scalarList X0(nSpecie, 0.0)
scalar chi(TrackCloudType &cloud, trackingData &td, const scalarField &X) const
Calculate the chi-factor for flash-boiling for the.
dimensionedScalar pos(const dimensionedScalar &ds)
const dimensionedScalar RR
Universal gas constant: default SI units: [J/kmol/K].
scalar y
scalar d0() const
Return const access to initial droplet diameter.
Definition: SprayParcelI.H:179
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))
dimensionedScalar cos(const dimensionedScalar &ds)
dimensionedScalar exp(const dimensionedScalar &ds)
dimensionedScalar cbrt(const dimensionedScalar &ds)
scalar tc() const
Return const access to atomisation characteristic time.
Definition: SprayParcelI.H:235
virtual scalar Ha(scalar p, scalar T) const =0
Liquid absolute enthalpy [J/kg].
The thermophysical properties of a liquid.
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
scalar injector_
Injected from injector (needed e.g. for calculating distance.
Definition: SprayParcel.H:180
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, Random &rndGen) const =0
scalar user_
Passive scalar (extra variable to be defined by user)
Definition: SprayParcel.H:186
Urel
Definition: pEqn.H:60
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:228
dimensionedScalar sin(const dimensionedScalar &ds)
void calc(TrackCloudType &cloud, trackingData &td, const scalar dt)
Update parcel properties over the time interval.
Definition: SprayParcel.C:61
scalar user() const
Return const access to passive user scalar.
Definition: SprayParcelI.H:263
scalar sigma_
Liquid surface tension [N/m].
Definition: SprayParcel.H:155
scalar Taverage(const scalar &Tliq, const scalar &Tc) const
Average temperature calculation.
dimensionedScalar pow3(const dimensionedScalar &ds)
const volScalarField & T
vector position0_
Injection position.
Definition: SprayParcel.H:152
#define R(A, B, C, D, E, F, K, M)
scalar injector() const
Return const access to injector id.
Definition: SprayParcelI.H:249
scalar yDot_
Rate of change of spherical deviation.
Definition: SprayParcel.H:170
PtrList< volScalarField > & Y
void solveTABEq(TrackCloudType &cloud, trackingData &td, const scalar dt)
Solve the TAB equation.
Definition: SprayParcel.C:365
void setCellValues(TrackCloudType &cloud, trackingData &td)
Set cell values.
Definition: SprayParcel.C:36
scalar ms_
Stripped parcel mass due to breakup.
Definition: SprayParcel.H:176
scalar Cp(const scalar p, const scalar T) const
Definition: EtoHthermo.H:2
void cellValueSourceCorrection(TrackCloudType &cloud, trackingData &td, const scalar dt)
Correct cell values using latest transfer information.
Definition: SprayParcel.C:48
dimensioned< scalar > mag(const dimensioned< Type > &)
scalar d0_
Initial droplet diameter.
Definition: SprayParcel.H:149
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:76
scalar tc_
Characteristic time (used in atomisation and/or breakup model)
Definition: SprayParcel.H:173
scalar mu_
Liquid dynamic viscosity [Pa.s].
Definition: SprayParcel.H:158
const dimensionedVector & g
Templated reacting parcel composition model class Consists of carrier species (via thermo package)...
Definition: ThermoCloud.H:61
scalar y_
Spherical deviation.
Definition: SprayParcel.H:167
ParcelType::trackingData trackingData
Use base tracking data.
Definition: SprayParcel.H:139
scalar KHindex_
Index for KH Breakup.
Definition: SprayParcel.H:164
scalarField Re(const UList< complex > &cf)
Definition: complexFields.C:97
const dimensionedScalar e
Elementary charge.
Definition: doubleScalar.H:105
Reaching spray parcel, with added functionality for atomisation and breakup.
Definition: SprayParcel.H:43
const liquidMixtureProperties & liquids() const
Return the global (additional) liquids.
scalar T0
Definition: createFields.H:22