SprayParcel.C 12.8 KB
Newer Older
andy's avatar
andy committed
1
2
3
4
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
5
    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
andy's avatar
andy committed
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#include "SprayParcel.H"
27
28
#include "CompositionModel.H"
#include "AtomizationModel.H"
andy's avatar
andy committed
29

30
// * * * * * * * * * * *  Protected Member Functions * * * * * * * * * * * * //
andy's avatar
andy committed
31
32
33
34
35
36
37
38
39
40

template<class ParcelType>
template<class TrackData>
void Foam::SprayParcel<ParcelType>::setCellValues
(
    TrackData& td,
    const scalar dt,
    const label cellI
)
{
41
    ParcelType::setCellValues(td, dt, cellI);
andy's avatar
andy committed
42
43
44
45
46
47
48
49
50
51
52
53
}


template<class ParcelType>
template<class TrackData>
void Foam::SprayParcel<ParcelType>::cellValueSourceCorrection
(
    TrackData& td,
    const scalar dt,
    const label cellI
)
{
54
    ParcelType::cellValueSourceCorrection(td, dt, cellI);
andy's avatar
andy committed
55
56
57
58
59
60
61
62
63
64
65
66
}


template<class ParcelType>
template<class TrackData>
void Foam::SprayParcel<ParcelType>::calc
(
    TrackData& td,
    const scalar dt,
    const label cellI
)
{
67
68
69
70
    typedef typename TrackData::cloudType::reactingCloudType reactingCloudType;
    const CompositionModel<reactingCloudType>& composition =
        td.cloud().composition();

71
    // Check if parcel belongs to liquid core
andy's avatar
andy committed
72
73
    if (liquidCore() > 0.5)
    {
74
        // Liquid core parcels should not experience coupled forces
75
        td.cloud().forces().setCalcCoupled(false);
76
77
    }

78
    // Get old mixture composition
79
80
81
    const scalarField& Y0(this->Y());
    scalarField X0(composition.liquids().X(Y0));

82
    // Check if we have critical or boiling conditions
83
84
85
86
87
    scalar TMax = composition.liquids().Tc(X0);
    const scalar T0 = this->T();
    const scalar pc0 = this->pc_;
    if (composition.liquids().pv(pc0, T0, X0) >= pc0*0.999)
    {
88
        // Set TMax to boiling temperature
89
        TMax = composition.liquids().pvInvert(pc0, X0);
andy's avatar
andy committed
90
91
    }

92
    // Set the maximum temperature limit
93
    td.cloud().constProps().setTMax(TMax);
94

95
    // Store the parcel properties
andy's avatar
andy committed
96
97
    this->Cp() = composition.liquids().Cp(pc0, T0, X0);
    sigma_ = composition.liquids().sigma(pc0, T0, X0);
98
    const scalar rho0 = composition.liquids().rho(pc0, T0, X0);
andy's avatar
andy committed
99
    this->rho() = rho0;
100
    const scalar mass0 = this->mass();
andy's avatar
andy committed
101
    mu_ = composition.liquids().mu(pc0, T0, X0);
andy's avatar
andy committed
102

103
    ParcelType::calc(td, dt, cellI);
andy's avatar
andy committed
104
105
106

    if (td.keepParticle)
    {
107
108
109
110
111
        // Reduce the stripped parcel mass due to evaporation
        // assuming the number of particles remains unchanged
        this->ms() -= this->ms()*(mass0 - this->mass())/mass0;

        // Update Cp, sigma, density and diameter due to change in temperature
112
        // and/or composition
andy's avatar
andy committed
113
114
        scalar T1 = this->T();
        const scalarField& Y1(this->Y());
115
        scalarField X1(composition.liquids().X(Y1));
andy's avatar
andy committed
116

117
        this->Cp() = composition.liquids().Cp(this->pc_, T1, X1);
andy's avatar
andy committed
118

andy's avatar
andy committed
119
        sigma_ = composition.liquids().sigma(this->pc_, T1, X1);
120
121

        scalar rho1 = composition.liquids().rho(this->pc_, T1, X1);
andy's avatar
andy committed
122
        this->rho() = rho1;
123

andy's avatar
andy committed
124
125
        mu_ = composition.liquids().mu(this->pc_, T1, X1);

andy's avatar
andy committed
126
        scalar d1 = this->d()*cbrt(rho0/rho1);
andy's avatar
andy committed
127
128
129
130
131
132
        this->d() = d1;

        if (liquidCore() > 0.5)
        {
            calcAtomization(td, dt, cellI);

133
            // Preserve the total mass/volume by increasing the number of
134
            // particles in parcels due to breakup
andy's avatar
andy committed
135
            scalar d2 = this->d();
136
            this->nParticle() *= pow3(d1/d2);
andy's avatar
andy committed
137
138
139
140
141
142
143
        }
        else
        {
            calcBreakup(td, dt, cellI);
        }
    }

144
    // Restore coupled forces
145
    td.cloud().forces().setCalcCoupled(true);
andy's avatar
andy committed
146
147
148
149
150
151
152
153
154
155
156
157
}


template<class ParcelType>
template<class TrackData>
void Foam::SprayParcel<ParcelType>::calcAtomization
(
    TrackData& td,
    const scalar dt,
    const label cellI
)
{
158
159
160
161
162
163
164
165
    typedef typename TrackData::cloudType::reactingCloudType reactingCloudType;
    const CompositionModel<reactingCloudType>& composition =
        td.cloud().composition();

    typedef typename TrackData::cloudType::sprayCloudType sprayCloudType;
    const AtomizationModel<sprayCloudType>& atomization =
        td.cloud().atomization();

andy's avatar
andy committed
166
    // Average molecular weight of carrier mix - assumes perfect gas
167
    scalar Wc = this->rhoc_*specie::RR*this->Tc()/this->pc();
andy's avatar
andy committed
168
    scalar R = specie::RR/Wc;
169
    scalar Tav = atomization.Taverage(this->T(), this->Tc());
andy's avatar
andy committed
170

171
    // Calculate average gas density based on average temperature
172
    scalar rhoAv = this->pc()/(R*Tav);
andy's avatar
andy committed
173

174
    scalar soi = td.cloud().injectors().timeStart();
175
    scalar currentTime = td.cloud().db().time().value();
andy's avatar
andy committed
176
177
178
    const vector& pos = this->position();
    const vector& injectionPos = this->position0();

179
    // Disregard the continous phase when calculating the relative velocity
andy's avatar
andy committed
180
181
182
    // (in line with the deactivated coupled assumption)
    scalar Urel = mag(this->U());

183
    scalar t0 = max(0.0, currentTime - this->age() - soi);
184
185
    scalar t1 = min(t0 + dt, td.cloud().injectors().timeEnd() - soi);

186
    // This should be the vol flow rate from when the parcel was injected
187
    scalar volFlowRate = td.cloud().injectors().volumeToInject(t0, t1)/dt;
andy's avatar
andy committed
188
189

    scalar chi = 0.0;
190
    if (atomization.calcChi())
andy's avatar
andy committed
191
    {
192
        chi = this->chi(td, composition.liquids().X(this->Y()));
andy's avatar
andy committed
193
194
    }

195
    atomization.update
andy's avatar
andy committed
196
197
198
199
200
    (
        dt,
        this->d(),
        this->liquidCore(),
        this->tc(),
201
202
203
        this->rho(),
        mu_,
        sigma_,
204
        volFlowRate,
andy's avatar
andy committed
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
        rhoAv,
        Urel,
        pos,
        injectionPos,
        td.cloud().pAmbient(),
        chi,
        td.cloud().rndGen()
    );
}


template<class ParcelType>
template<class TrackData>
void Foam::SprayParcel<ParcelType>::calcBreakup
(
    TrackData& td,
    const scalar dt,
    const label cellI
)
{
225
226
227
228
229
230
231
    typedef typename TrackData::cloudType cloudType;
    typedef typename cloudType::parcelType parcelType;
    typedef typename cloudType::forceType forceType;

    const parcelType& p = static_cast<const parcelType&>(*this);
    const forceType& forces = td.cloud().forces();

andy's avatar
andy committed
232
233
234
235
236
237
    if (td.cloud().breakup().solveOscillationEq())
    {
        solveTABEq(td, dt);
    }

    // Average molecular weight of carrier mix - assumes perfect gas
238
    scalar Wc = this->rhoc()*specie::RR*this->Tc()/this->pc();
andy's avatar
andy committed
239
    scalar R = specie::RR/Wc;
240
    scalar Tav = td.cloud().atomization().Taverage(this->T(), this->Tc());
andy's avatar
andy committed
241

242
    // Calculate average gas density based on average temperature
243
244
245
    scalar rhoAv = this->pc()/(R*Tav);
    scalar muAv = this->muc();
    vector Urel = this->U() - this->Uc();
andy's avatar
andy committed
246
    scalar Urmag = mag(Urel);
247
    scalar Re = this->Re(this->U(), this->d(), rhoAv, muAv);
andy's avatar
andy committed
248

249
250
251
    const scalar mass = p.mass();
    const forceSuSp Fcp = forces.calcCoupled(p, dt, mass, Re, muAv);
    const forceSuSp Fncp = forces.calcNonCoupled(p, dt, mass, Re, muAv);
252
    this->tMom() = mass/(Fcp.Sp() + Fncp.Sp());
andy's avatar
andy committed
253
254
255

    const vector g = td.cloud().g().value();

Henry's avatar
Henry committed
256
    scalar parcelMassChild = 0.0;
andy's avatar
andy committed
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
    scalar dChild = 0.0;
    if
    (
        td.cloud().breakup().update
        (
            dt,
            g,
            this->d(),
            this->tc(),
            this->ms(),
            this->nParticle(),
            this->KHindex(),
            this->y(),
            this->yDot(),
            this->d0(),
272
273
274
            this->rho(),
            mu_,
            sigma_,
andy's avatar
andy committed
275
276
277
278
279
            this->U(),
            rhoAv,
            muAv,
            Urel,
            Urmag,
280
            this->tMom(),
andy's avatar
andy committed
281
            dChild,
Henry's avatar
Henry committed
282
            parcelMassChild
andy's avatar
andy committed
283
284
285
286
287
        )
    )
    {
        scalar Re = rhoAv*Urmag*dChild/muAv;

288
289
        // Add child parcel as copy of parent
        SprayParcel<ParcelType>* child = new SprayParcel<ParcelType>(*this);
andy's avatar
andy committed
290
        child->d() = dChild;
Henry's avatar
Henry committed
291
292
293
294
        child->d0() = dChild;
        const scalar massChild = child->mass();
        child->mass0() = massChild;
        child->nParticle() = parcelMassChild/massChild;
295
296
297
298
299
300

        const forceSuSp Fcp =
            forces.calcCoupled(*child, dt, massChild, Re, muAv);
        const forceSuSp Fncp =
            forces.calcNonCoupled(*child, dt, massChild, Re, muAv);

Henry's avatar
Henry committed
301
        child->age() = 0.0;
andy's avatar
andy committed
302
303
304
305
306
307
308
        child->liquidCore() = 0.0;
        child->KHindex() = 1.0;
        child->y() = td.cloud().breakup().y0();
        child->yDot() = td.cloud().breakup().yDot0();
        child->tc() = -GREAT;
        child->ms() = 0.0;
        child->injector() = this->injector();
309
        child->tMom() = massChild/(Fcp.Sp() + Fncp.Sp());
andy's avatar
andy committed
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
        child->user() = 0.0;
        child->setCellValues(td, dt, cellI);

        td.cloud().addParticle(child);
    }
}


template<class ParcelType>
template<class TrackData>
Foam::scalar Foam::SprayParcel<ParcelType>::chi
(
    TrackData& td,
    const scalarField& X
) const
{
326
    // Modifications to take account of the flash boiling on primary break-up
andy's avatar
andy committed
327

328
329
330
    typedef typename TrackData::cloudType::reactingCloudType reactingCloudType;
    const CompositionModel<reactingCloudType>& composition =
        td.cloud().composition();
andy's avatar
andy committed
331
332

    scalar chi = 0.0;
333
334
    scalar T0 = this->T();
    scalar p0 = this->pc();
andy's avatar
andy committed
335
336
    scalar pAmb = td.cloud().pAmbient();

337
    scalar pv = composition.liquids().pv(p0, T0, X);
andy's avatar
andy committed
338

339
340
341
    forAll(composition.liquids(), i)
    {
        if (pv >= 0.999*pAmb)
andy's avatar
andy committed
342
        {
343
            // Liquid is boiling - calc boiling temperature
andy's avatar
andy committed
344

345
            const liquidProperties& liq = composition.liquids().properties()[i];
346
            scalar TBoil = liq.pvInvert(p0);
andy's avatar
andy committed
347

348
            scalar hl = liq.hl(pAmb, TBoil);
Henry's avatar
Henry committed
349
            scalar iTp = liq.h(pAmb, T0) - pAmb/liq.rho(pAmb, T0);
350
            scalar iTb = liq.h(pAmb, TBoil) - pAmb/liq.rho(pAmb, TBoil);
andy's avatar
andy committed
351

352
            chi += X[i]*(iTp - iTb)/hl;
andy's avatar
andy committed
353
354
355
        }
    }

356
    chi = min(1.0, max(chi, 0.0));
andy's avatar
andy committed
357

358
    return chi;
andy's avatar
andy committed
359
360
}

361

andy's avatar
andy committed
362
363
364
365
366
template<class ParcelType>
template<class TrackData>
void Foam::SprayParcel<ParcelType>::solveTABEq
(
    TrackData& td,
367
    const scalar dt
andy's avatar
andy committed
368
369
370
371
372
373
)
{
    const scalar& TABCmu = td.cloud().breakup().TABCmu();
    const scalar& TABWeCrit = td.cloud().breakup().TABWeCrit();
    const scalar& TABComega = td.cloud().breakup().TABComega();

374
    scalar r = 0.5*this->d();
andy's avatar
andy committed
375
376
377
    scalar r2 = r*r;
    scalar r3 = r*r2;

378
    // Inverse of characteristic viscous damping time
andy's avatar
andy committed
379
    scalar rtd = 0.5*TABCmu*mu_/(this->rho()*r2);
andy's avatar
andy committed
380

381
    // Oscillation frequency (squared)
andy's avatar
andy committed
382
    scalar omega2 = TABComega*sigma_/(this->rho()*r3) - rtd*rtd;
andy's avatar
andy committed
383

andy's avatar
andy committed
384
    if (omega2 > 0)
andy's avatar
andy committed
385
386
    {
        scalar omega = sqrt(omega2);
387
        scalar rhoc = this->rhoc();
andy's avatar
andy committed
388
        scalar We = this->We(this->U(), r, rhoc, sigma_)/TABWeCrit;
andy's avatar
andy committed
389

andy's avatar
andy committed
390
        scalar y1 = this->y() - We;
andy's avatar
andy committed
391
392
        scalar y2 = this->yDot()/omega;

393
        // Update distortion parameters
andy's avatar
andy committed
394
395
396
397
398
        scalar c = cos(omega*dt);
        scalar s = sin(omega*dt);
        scalar e = exp(-rtd*dt);
        y2 = (this->yDot() + y1*rtd)/omega;

andy's avatar
andy committed
399
        this->y() = We + e*(y1*c + y2*s);
andy's avatar
andy committed
400
401
402
403
404
405
406
        if (this->y() < 0)
        {
            this->y() = 0.0;
            this->yDot() = 0.0;
        }
        else
        {
andy's avatar
andy committed
407
            this->yDot() = (We - this->y())*rtd + e*omega*(y2*c - y1*s);
andy's avatar
andy committed
408
409
410
411
        }
    }
    else
    {
412
        // Reset distortion parameters
andy's avatar
andy committed
413
414
415
416
417
        this->y() = 0;
        this->yDot() = 0;
    }
}

418
419
420

// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

Henry's avatar
Henry committed
421
template<class ParcelType>
422
423
424
425
426
Foam::SprayParcel<ParcelType>::SprayParcel(const SprayParcel<ParcelType>& p)
:
    ParcelType(p),
    d0_(p.d0_),
    position0_(p.position0_),
427
428
    sigma_(p.sigma_),
    mu_(p.mu_),
429
430
431
432
433
434
435
436
437
438
439
440
    liquidCore_(p.liquidCore_),
    KHindex_(p.KHindex_),
    y_(p.y_),
    yDot_(p.yDot_),
    tc_(p.tc_),
    ms_(p.ms_),
    injector_(p.injector_),
    tMom_(p.tMom_),
    user_(p.user_)
{}


Henry's avatar
Henry committed
441
template<class ParcelType>
442
443
Foam::SprayParcel<ParcelType>::SprayParcel
(
444
445
    const SprayParcel<ParcelType>& p,
    const polyMesh& mesh
446
447
)
:
448
    ParcelType(p, mesh),
449
450
    d0_(p.d0_),
    position0_(p.position0_),
451
452
    sigma_(p.sigma_),
    mu_(p.mu_),
453
454
455
456
457
458
459
460
461
462
463
464
    liquidCore_(p.liquidCore_),
    KHindex_(p.KHindex_),
    y_(p.y_),
    yDot_(p.yDot_),
    tc_(p.tc_),
    ms_(p.ms_),
    injector_(p.injector_),
    tMom_(p.tMom_),
    user_(p.user_)
{}


andy's avatar
andy committed
465
466
467
468
469
470
// * * * * * * * * * * * * * * IOStream operators  * * * * * * * * * * * * * //

#include "SprayParcelIO.C"


// ************************************************************************* //