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 | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2011-2017 OpenFOAM Foundation
9-------------------------------------------------------------------------------
10License
11 This file is part of OpenFOAM.
12
13 OpenFOAM is free software: you can redistribute it and/or modify it
14 under the terms of the GNU General Public License as published by
15 the Free Software Foundation, either version 3 of the License, or
16 (at your option) any later version.
17
18 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21 for more details.
22
23 You should have received a copy of the GNU General Public License
24 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25
26\*---------------------------------------------------------------------------*/
27
28#include "SprayParcel.H"
29#include "BreakupModel.H"
30#include "CompositionModel.H"
31#include "AtomizationModel.H"
32
33// * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * * //
34
35template<class ParcelType>
36template<class TrackCloudType>
38(
39 TrackCloudType& cloud,
40 trackingData& td
41)
42{
44}
45
46
47template<class ParcelType>
48template<class TrackCloudType>
50(
51 TrackCloudType& cloud,
52 trackingData& td,
53 const scalar dt
54)
55{
57}
58
59
60template<class ParcelType>
61template<class TrackCloudType>
63(
64 TrackCloudType& cloud,
65 trackingData& td,
66 const scalar dt
67)
68{
69 const auto& composition = cloud.composition();
70 const auto& liquids = composition.liquids();
71
72 // Check if parcel belongs to liquid core
73 if (liquidCore() > 0.5)
74 {
75 // Liquid core parcels should not experience coupled forces
76 cloud.forces().setCalcCoupled(false);
77 }
78
79 // Get old mixture composition
80 scalarField X0(liquids.X(this->Y()));
81
82 // Check if we have critical or boiling conditions
83 scalar TMax = liquids.Tc(X0);
84 const scalar T0 = this->T();
85 const scalar pc0 = td.pc();
86 if (liquids.pv(pc0, T0, X0) >= pc0*0.999)
87 {
88 // Set TMax to boiling temperature
89 TMax = liquids.pvInvert(pc0, X0);
90 }
91
92 // Set the maximum temperature limit
93 cloud.constProps().setTMax(TMax);
94
95 // Store the parcel properties
96 this->Cp() = liquids.Cp(pc0, T0, X0);
97 sigma_ = liquids.sigma(pc0, T0, X0);
98 const scalar rho0 = liquids.rho(pc0, T0, X0);
99 this->rho() = rho0;
100 const scalar mass0 = this->mass();
101 mu_ = liquids.mu(pc0, T0, X0);
102
103 ParcelType::calc(cloud, td, dt);
104
105 if (td.keepParticle)
106 {
107 // Reduce the stripped parcel mass due to evaporation
108 // assuming the number of particles remains unchanged
109 this->ms() -= this->ms()*(mass0 - this->mass())/mass0;
110
111 // Update Cp, sigma, density and diameter due to change in temperature
112 // and/or composition
113 scalar T1 = this->T();
114 scalarField X1(liquids.X(this->Y()));
115
116 this->Cp() = liquids.Cp(td.pc(), T1, X1);
117
118 sigma_ = liquids.sigma(td.pc(), T1, X1);
119
120 scalar rho1 = liquids.rho(td.pc(), T1, X1);
121 this->rho() = rho1;
122
123 mu_ = liquids.mu(td.pc(), T1, X1);
124
125 scalar d1 = this->d()*cbrt(rho0/rho1);
126 this->d() = d1;
127
128 if (liquidCore() > 0.5)
129 {
130 calcAtomization(cloud, td, dt);
131
132 // Preserve the total mass/volume by increasing the number of
133 // particles in parcels due to breakup
134 scalar d2 = this->d();
135 this->nParticle() *= pow3(d1/d2);
136 }
137 else
138 {
139 calcBreakup(cloud, td, dt);
140 }
141 }
142
143 // Restore coupled forces
144 cloud.forces().setCalcCoupled(true);
145}
146
147
148template<class ParcelType>
149template<class TrackCloudType>
151(
152 TrackCloudType& cloud,
153 trackingData& td,
154 const scalar dt
155)
156{
157 const auto& atomization = cloud.atomization();
158
159 if (!atomization.active())
160 {
161 return;
162 }
163
164 const auto& composition = cloud.composition();
165 const auto& liquids = composition.liquids();
166
167 // Average molecular weight of carrier mix - assumes perfect gas
168 scalar Wc = td.rhoc()*RR*td.Tc()/td.pc();
169 scalar R = RR/Wc;
170 scalar Tav = atomization.Taverage(this->T(), td.Tc());
171
172 // Calculate average gas density based on average temperature
173 scalar rhoAv = td.pc()/(R*Tav);
174
175 scalar soi = cloud.injectors().timeStart();
176 scalar currentTime = cloud.db().time().value();
177 const vector pos(this->position());
178 const vector& injectionPos = this->position0();
179
180 // Disregard the continuous phase when calculating the relative velocity
181 // (in line with the deactivated coupled assumption)
182 scalar Urel = mag(this->U());
183
184 scalar t0 = max(0.0, currentTime - this->age() - soi);
185 scalar t1 = min(t0 + dt, cloud.injectors().timeEnd() - soi);
186
187 // This should be the vol flow rate from when the parcel was injected
188 scalar volFlowRate = cloud.injectors().volumeToInject(t0, t1)/dt;
189
190 scalar chi = 0.0;
191 if (atomization.calcChi())
192 {
193 chi = this->chi(cloud, td, liquids.X(this->Y()));
194 }
195
196 atomization.update
197 (
198 dt,
199 this->d(),
200 this->liquidCore(),
201 this->tc(),
202 this->rho(),
203 mu_,
204 sigma_,
205 volFlowRate,
206 rhoAv,
207 Urel,
208 pos,
209 injectionPos,
210 cloud.pAmbient(),
211 chi,
212 cloud.rndGen()
213 );
214}
215
216
217template<class ParcelType>
218template<class TrackCloudType>
220(
221 TrackCloudType& cloud,
222 trackingData& td,
223 const scalar dt
224)
225{
226 auto& breakup = cloud.breakup();
227
228 if (!breakup.active())
229 {
230 return;
231 }
232
233 if (breakup.solveOscillationEq())
234 {
235 solveTABEq(cloud, td, dt);
236 }
237
238 // Average molecular weight of carrier mix - assumes perfect gas
239 scalar Wc = td.rhoc()*RR*td.Tc()/td.pc();
240 scalar R = RR/Wc;
241 scalar Tav = cloud.atomization().Taverage(this->T(), td.Tc());
242
243 // Calculate average gas density based on average temperature
244 scalar rhoAv = td.pc()/(R*Tav);
245 scalar muAv = td.muc();
246 vector Urel = this->U() - td.Uc();
247 scalar Urmag = mag(Urel);
248 scalar Re = this->Re(rhoAv, this->U(), td.Uc(), this->d(), muAv);
249
250 const typename TrackCloudType::parcelType& p =
251 static_cast<const typename TrackCloudType::parcelType&>(*this);
252 typename TrackCloudType::parcelType::trackingData& ttd =
253 static_cast<typename TrackCloudType::parcelType::trackingData&>(td);
254 const scalar mass = p.mass();
255 const typename TrackCloudType::forceType& forces = cloud.forces();
256 const forceSuSp Fcp = forces.calcCoupled(p, ttd, dt, mass, Re, muAv);
257 const forceSuSp Fncp = forces.calcNonCoupled(p, ttd, dt, mass, Re, muAv);
258 this->tMom() = mass/(Fcp.Sp() + Fncp.Sp() + ROOTVSMALL);
259
260 const vector g = cloud.g().value();
261
262 scalar parcelMassChild = 0.0;
263 scalar dChild = 0.0;
264 if
265 (
266 breakup.update
267 (
268 dt,
269 g,
270 this->d(),
271 this->tc(),
272 this->ms(),
273 this->nParticle(),
274 this->KHindex(),
275 this->y(),
276 this->yDot(),
277 this->d0(),
278 this->rho(),
279 mu_,
280 sigma_,
281 this->U(),
282 rhoAv,
283 muAv,
284 Urel,
285 Urmag,
286 this->tMom(),
287 dChild,
288 parcelMassChild
289 )
290 )
291 {
292 scalar Re = rhoAv*Urmag*dChild/muAv;
293
294 // Add child parcel as copy of parent
296 child->origId() = this->getNewParticleID();
297 child->origProc() = Pstream::myProcNo();
298 child->d() = dChild;
299 child->d0() = dChild;
300 const scalar massChild = child->mass();
301 child->mass0() = massChild;
302 child->nParticle() = parcelMassChild/massChild;
303
304 const forceSuSp Fcp =
305 forces.calcCoupled(*child, ttd, dt, massChild, Re, muAv);
306 const forceSuSp Fncp =
307 forces.calcNonCoupled(*child, ttd, dt, massChild, Re, muAv);
308
309 child->age() = 0.0;
310 child->liquidCore() = 0.0;
311 child->KHindex() = 1.0;
312 child->y() = cloud.breakup().y0();
313 child->yDot() = cloud.breakup().yDot0();
314 child->tc() = 0.0;
315 child->ms() = -GREAT;
316 child->injector() = this->injector();
317 child->tMom() = massChild/(Fcp.Sp() + Fncp.Sp());
318 child->user() = 0.0;
319 child->calcDispersion(cloud, td, dt);
320
321 cloud.addParticle(child);
322 }
323}
324
325
326template<class ParcelType>
327template<class TrackCloudType>
329(
330 TrackCloudType& cloud,
331 trackingData& td,
332 const scalarField& X
333) const
334{
335 // Modifications to take account of the flash boiling on primary break-up
336
337 const auto& composition = cloud.composition();
338 const auto& liquids = composition.liquids();
339
340 scalar chi = 0.0;
341 scalar T0 = this->T();
342 scalar p0 = td.pc();
343 scalar pAmb = cloud.pAmbient();
344
345 scalar pv = liquids.pv(p0, T0, X);
346
347 forAll(liquids, i)
348 {
349 if (pv >= 0.999*pAmb)
350 {
351 // Liquid is boiling - calc boiling temperature
352
353 const liquidProperties& liq = liquids.properties()[i];
354 scalar TBoil = liq.pvInvert(p0);
355
356 scalar hl = liq.hl(pAmb, TBoil);
357 scalar iTp = liq.h(pAmb, T0) - pAmb/liq.rho(pAmb, T0);
358 scalar iTb = liq.h(pAmb, TBoil) - pAmb/liq.rho(pAmb, TBoil);
359
360 chi += X[i]*(iTp - iTb)/hl;
361 }
362 }
363
364 chi = min(1.0, max(chi, 0.0));
365
366 return chi;
367}
368
369
370template<class ParcelType>
371template<class TrackCloudType>
373(
374 TrackCloudType& cloud,
375 trackingData& td,
376 const scalar dt
377)
378{
379 const scalar& TABCmu = cloud.breakup().TABCmu();
380 const scalar& TABtwoWeCrit = cloud.breakup().TABtwoWeCrit();
381 const scalar& TABComega = cloud.breakup().TABComega();
382
383 scalar r = 0.5*this->d();
384 scalar r2 = r*r;
385 scalar r3 = r*r2;
386
387 // Inverse of characteristic viscous damping time
388 scalar rtd = 0.5*TABCmu*mu_/(this->rho()*r2);
389
390 // Oscillation frequency (squared)
391 scalar omega2 = TABComega*sigma_/(this->rho()*r3) - rtd*rtd;
392
393 if (omega2 > 0)
394 {
395 scalar omega = sqrt(omega2);
396 scalar We =
397 this->We(td.rhoc(), this->U(), td.Uc(), r, sigma_)/TABtwoWeCrit;
398
399 // Initial values for y and yDot
400 scalar y0 = this->y() - We;
401 scalar yDot0 = this->yDot() + y0*rtd;
402
403 // Update distortion parameters
404 scalar c = cos(omega*dt);
405 scalar s = sin(omega*dt);
406 scalar e = exp(-rtd*dt);
407
408 this->y() = We + e*(y0*c + (yDot0/omega)*s);
409 this->yDot() = (We - this->y())*rtd + e*(yDot0*c - omega*y0*s);
410 }
411 else
412 {
413 // Reset distortion parameters
414 this->y() = 0;
415 this->yDot() = 0;
416 }
417}
418
419
420// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
421
422template<class ParcelType>
424:
425 ParcelType(p),
426 d0_(p.d0_),
427 position0_(p.position0_),
428 sigma_(p.sigma_),
429 mu_(p.mu_),
430 liquidCore_(p.liquidCore_),
431 KHindex_(p.KHindex_),
432 y_(p.y_),
433 yDot_(p.yDot_),
434 tc_(p.tc_),
435 ms_(p.ms_),
436 injector_(p.injector_),
437 tMom_(p.tMom_),
438 user_(p.user_)
439{}
440
441
442template<class ParcelType>
444(
445 const SprayParcel<ParcelType>& p,
446 const polyMesh& mesh
447)
448:
449 ParcelType(p, mesh),
450 d0_(p.d0_),
451 position0_(p.position0_),
452 sigma_(p.sigma_),
453 mu_(p.mu_),
454 liquidCore_(p.liquidCore_),
455 KHindex_(p.KHindex_),
456 y_(p.y_),
457 yDot_(p.yDot_),
458 tc_(p.tc_),
459 ms_(p.ms_),
460 injector_(p.injector_),
461 tMom_(p.tMom_),
462 user_(p.user_)
463{}
464
465
466// * * * * * * * * * * * * * * IOStream operators * * * * * * * * * * * * * //
467
468#include "SprayParcelIO.C"
469
470
471// ************************************************************************* //
scalar y
#define R(A, B, C, D, E, F, K, M)
Y[inertIndex] max(0.0)
const uniformDimensionedVectorField & g
volScalarField & rho1
const objectRegistry & db() const noexcept
Return the local objectRegistry.
Definition: IOobject.C:500
void setCellValues(TrackCloudType &cloud, trackingData &td)
Set cell values.
Reacting spray parcel, with added functionality for atomization and breakup.
Definition: SprayParcel.H:63
scalar injector() const
Return const access to injector id.
Definition: SprayParcelI.H:315
scalar liquidCore() const
Return const access to liquid core.
Definition: SprayParcelI.H:273
scalar tMom() const
Return const access to momentum relaxation time.
Definition: SprayParcelI.H:322
void calcAtomization(TrackCloudType &cloud, trackingData &td, const scalar dt)
Correct parcel properties according to atomization model.
Definition: SprayParcel.C:151
scalar yDot() const
Return const access to rate of change of spherical deviation.
Definition: SprayParcelI.H:294
scalar d0() const
Return const access to initial droplet diameter.
Definition: SprayParcelI.H:245
void solveTABEq(TrackCloudType &cloud, trackingData &td, const scalar dt)
Solve the TAB equation.
Definition: SprayParcel.C:373
scalar KHindex() const
Return const access to Kelvin-Helmholtz breakup index.
Definition: SprayParcelI.H:280
scalar tc() const
Return const access to atomization characteristic time.
Definition: SprayParcelI.H:301
void setCellValues(TrackCloudType &cloud, trackingData &td)
Set cell values.
Definition: SprayParcel.C:38
void calcBreakup(TrackCloudType &cloud, trackingData &td, const scalar dt)
Correct parcel properties according to breakup model.
Definition: SprayParcel.C:220
scalar ms() const
Return const access to stripped parcel mass.
Definition: SprayParcelI.H:308
scalar y() const
Return const access to spherical deviation.
Definition: SprayParcelI.H:287
scalar user() const
Return const access to passive user scalar.
Definition: SprayParcelI.H:329
const Switch cellValueSourceCorrection() const
Return const access to the cell value correction flag.
A cloud is a registry collection of lagrangian particles.
Definition: cloud.H:60
const Type & value() const
Return const reference to value.
Class used to pass data into container.
Helper container for force Su and Sp terms.
Definition: forceSuSp.H:67
scalar Sp() const
Return const access to the implicit coefficient [kg/s].
Definition: forceSuSpI.H:67
The thermophysical properties of a liquid.
virtual scalar hl(scalar p, scalar T) const =0
Heat of vapourisation [J/kg].
virtual scalar h(scalar p, scalar T) const =0
Liquid enthalpy [J/kg] - reference to 298.15 K.
virtual scalar pvInvert(scalar p) const
Invert the vapour pressure relationship to retrieve the.
const Time & time() const noexcept
Return time registry.
virtual scalar rho(scalar p, scalar T) const =0
Liquid density [kg/m^3].
U
Definition: pEqn.H:72
basicSpecieMixture & composition
volScalarField & p
const volScalarField & T
const volScalarField & Cp
Definition: EEqn.H:7
const volScalarField & p0
Definition: EEqn.H:36
dynamicFvMesh & mesh
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.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))
Urel
Definition: pEqn.H:56
const scalar RR
Universal gas constant: default in [J/(kmol K)].
dimensionedScalar pos(const dimensionedScalar &ds)
dimensionedScalar exp(const dimensionedScalar &ds)
scalarField Re(const UList< complex > &cf)
Extract real component.
Definition: complexField.C:159
dimensionedScalar pow3(const dimensionedScalar &ds)
dimensionedScalar y0(const dimensionedScalar &ds)
dimensionedScalar sin(const dimensionedScalar &ds)
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
dimensionedScalar cbrt(const dimensionedScalar &ds)
dimensionedScalar cos(const dimensionedScalar &ds)
scalar rho0
scalarList X0(nSpecie, Zero)
scalar T0
Definition: createFields.H:22
volScalarField & e
Definition: createFields.H:11
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333