ReactingParcel.C 11.9 KB
Newer Older
1
2
3
4
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
5
    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
     \\/     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 2 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, write to the Free Software Foundation,
    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA

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

#include "ReactingParcel.H"

// * * * * * * * * * * *  Protected Member Functions * * * * * * * * * * * * //

template<class ParcelType>
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
template<class TrackData>
void Foam::ReactingParcel<ParcelType>::updateCellQuantities
(
    TrackData& td,
    const scalar dt,
    const label celli
)
{
    ThermoParcel<ParcelType>::updateCellQuantities(td, dt, celli);

    pc_ = td.pInterp().interpolate(this->position(), celli);
}


template<class ParcelType>
template<class TrackData>
48
49
void Foam::ReactingParcel<ParcelType>::calcCoupled
(
50
    TrackData& td,
51
    const scalar dt,
52
    const label celli
53
54
55
56
57
)
{
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // Define local properties at beginning of timestep
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
58
    const vector U0 = this->U_;
59
    const scalar mass0 = this->mass();
Andrew Heather's avatar
Andrew Heather committed
60
    const scalar cp0 = this->cp_;
61
62
63
64
65
66
67
68
69
70
71
72
    const scalar np0 = this->nParticle_;
    const scalar T0 = this->T_;

    // ~~~~~~~~~~~~~~~~~~~~~~~~~
    // Initialise transfer terms
    // ~~~~~~~~~~~~~~~~~~~~~~~~~

    // Momentum transfer from the particle to the carrier phase
    vector dUTrans = vector::zero;

    // Enthalpy transfer from the particle to the carrier phase
    scalar dhTrans = 0.0;
73
74
75
76
77
78
79
80
81

    // Mass transfer from particle to carrier phase
    // - components exist in particle already
    scalarList dMassMT(td.cloud().gases().size(), 0.0);

    // Mass transfer due to surface reactions from particle to carrier phase
    // - components do not necessarily exist in particle already
    scalarList dMassSR(td.cloud().gases().size(), 0.0);

82
83
84
85
86
87

    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // Calculate velocity - update U
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    scalar Cud = 0.0;
    const vector U1 = calcVelocity(td, dt, Cud, dUTrans);
88
89
90
91
92
93


    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // Calculate heat transfer - update T
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    scalar htc = 0.0;
94
    scalar T1 = calcHeatTransfer(td, dt, celli, htc, dhTrans);
95
96
97
98
99
100
101
102
103
104
105


    // ~~~~~~~~~~~~~~~~~~~~~~~
    // Calculate mass transfer
    // ~~~~~~~~~~~~~~~~~~~~~~~
    calcMassTransfer(td, dt, T0, T1, dMassMT);


    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // Calculate surface reactions
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
106
107
108
109
110

    // Initialise enthalpy retention to zero
    scalar dhRet = 0.0;

    calcSurfaceReactions(td, dt, celli, T0, T1, dMassMT, dMassSR, dhRet);
111
112

    // New total mass
113
    const scalar mass1 = mass0 - sum(dMassMT) - sum(dMassSR);
114

115
116
117
118
119
    // Correct particle temperature to account for latent heat of
    // devolatilisation
    T1 -=
        td.constProps().Ldevol()
       *sum(dMassMT)
Andrew Heather's avatar
Andrew Heather committed
120
       /(0.5*(mass0 + mass1)*cp0);
121

122
123
    // Add retained enthalpy from surface reaction to particle and remove
    // from gas
Andrew Heather's avatar
Andrew Heather committed
124
    T1 += dhRet/(0.5*(mass0 + mass1)*cp0);
125
126
127
128
129
130
    dhTrans -= dhRet;

    // Correct dhTrans to account for enthalpy of evolved volatiles
    dhTrans +=
        sum(dMassMT)
       *td.cloud().composition().HGas(YGas_, 0.5*(T0 + T1));
131

132
133
134
135
    // Correct dhTrans to account for enthalpy of consumed solids
    dhTrans +=
        sum(dMassSR)
       *td.cloud().composition().HSolid(YSolid_, 0.5*(T0 + T1));
136
137
138
139
140
141
142
143
144


    // ~~~~~~~~~~~~~~~~~~~~~~~
    // Accumulate source terms
    // ~~~~~~~~~~~~~~~~~~~~~~~

    // Transfer mass lost from particle to carrier mass source
    forAll(dMassMT, i)
    {
145
        td.cloud().rhoTrans(i)[celli] += np0*(dMassMT[i] + dMassSR[i]);
146
147
148
    }

    // Update momentum transfer
149
    td.cloud().UTrans()[celli] += np0*dUTrans;
150
151
152
153
154

    // Accumulate coefficient to be applied in carrier phase momentum coupling
    td.cloud().UCoeff()[celli] += np0*mass0*Cud;

    // Update enthalpy transfer
155
156
    // - enthalpy of lost solids already accounted for
    td.cloud().hTrans()[celli] += np0*dhTrans;
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173

    // Accumulate coefficient to be applied in carrier phase enthalpy coupling
    td.cloud().hCoeff()[celli] += np0*htc*this->areaS();


    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // Remove the particle when mass falls below minimum threshold
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    if (mass1 < td.constProps().minParticleMass())
    {
        td.keepParticle = false;

        // Absorb particle(s) into carrier phase
        forAll(dMassMT, i)
        {
            td.cloud().rhoTrans(i)[celli] += np0*dMassMT[i];
        }
174
175
176
177
178
179
        td.cloud().hTrans()[celli] +=
            np0*mass1
           *(
                YMixture_[0]*td.cloud().composition().HGas(YGas_, T1)
              + YMixture_[2]*td.cloud().composition().HSolid(YSolid_, T1)
            );
180
181
182
183
184
185
186
        td.cloud().UTrans()[celli] += np0*mass1*U1;
    }
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // Set new particle properties
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
    else
    {
187
188
        this->U_ = U1;
        this->T_ = T1;
189
190
191
192
        this->cp_ =
            YMixture_[0]*td.cloud().composition().cpGas(YGas_, T1)
          + YMixture_[1]*td.cloud().composition().cpLiquid(YLiquid_, pc_, T1)
          + YMixture_[2]*td.cloud().composition().cpSolid(YSolid_);
193
194

        // Update particle density or diameter
195
        if (td.constProps().constantVolume())
196
        {
197
            this->rho_ = mass1/this->volume();
198
199
200
        }
        else
        {
201
            this->d_ = cbrt(mass1/this->rho_*6.0/mathematicalConstant::pi);
202
203
204
205
206
207
        }
    }
}


template<class ParcelType>
208
template<class TrackData>
209
210
void Foam::ReactingParcel<ParcelType>::calcUncoupled
(
211
    TrackData& td,
212
    const scalar dt,
213
    const label celli
214
215
216
217
218
)
{
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // Define local properties at beginning of timestep
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
219
    const scalar T0 = this->T_;
220
    const scalar mass0 = this->mass();
221
    const scalar cp0 = this->cp_;
222

223
224
225
226
227
228
229
230
231
232
    // ~~~~~~~~~~~~~~~~~~~~~~~~~
    // Initialise transfer terms
    // ~~~~~~~~~~~~~~~~~~~~~~~~~

    // Momentum transfer from the particle to the carrier phase
    vector dUTrans = vector::zero;

    // Enthalpy transfer from the particle to the carrier phase
    scalar dhTrans = 0.0;

233
234
235
236
237
238
239
240
    // Mass transfer from particle to carrier phase
    // - components exist in particle already
    scalarList dMassMT(td.cloud().gases().size(), 0.0);

    // Mass transfer due to surface reactions from particle to carrier phase
    // - components do not necessarily exist in particle already
    scalarList dMassSR(td.cloud().gases().size(), 0.0);

241
242
243
244
245
246

    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // Calculate velocity - update U
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    scalar Cud = 0.0;
    const vector U1 = calcVelocity(td, dt, Cud, dUTrans);
247
248
249
250
251
252


    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // Calculate heat transfer - update T
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    scalar htc = 0.0;
253
    scalar T1 = calcHeatTransfer(td, dt, celli, htc, dhTrans);
254
255
256
257
258
259
260
261
262
263
264
265


    // ~~~~~~~~~~~~~~~~~~~~~~~
    // Calculate mass transfer
    // ~~~~~~~~~~~~~~~~~~~~~~~
    calcMassTransfer(td, dt, T0, T1, dMassMT);


    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // Calculate surface reactions
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~

266
267
    // Initialise enthalpy retention to zero
    scalar dhRet = 0.0;
268

269
    calcSurfaceReactions(td, dt, celli, T0, T1, dMassMT, dMassSR, dhRet);
270

271
272
    // New total mass
    const scalar mass1 = mass0 - sum(dMassMT) - sum(dMassSR);
273

274
275
276
277
278
    // New specific heat capacity
    const scalar cp1 =
        YMixture_[0]*td.cloud().composition().cpGas(YGas_, T1)
      + YMixture_[1]*td.cloud().composition().cpLiquid(YLiquid_, pc_, T1)
      + YMixture_[2]*td.cloud().composition().cpSolid(YSolid_);
279

280
281
    // Add retained enthalpy to particle
    T1 += dhRet/(mass0*0.5*(cp0 + cp1));
282
283
284
285
286
287
288
289
290
291
292
293
294

    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // Remove the particle when mass falls below minimum threshold
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    if (mass1 < td.constProps().minParticleMass())
    {
        td.keepParticle = false;
    }
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // Set new particle properties
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
    else
    {
295
296
297
        this->U_ = U1;
        this->T_ = T1;
        this->cp_ = cp1;
298
299

        // Update particle density or diameter
300
        if (td.constProps().constantVolume())
301
        {
302
            this->rho_ = mass1/this->volume();
303
304
305
        }
        else
        {
306
            this->d_ = cbrt(mass1/this->rho_*6.0/mathematicalConstant::pi);
307
308
309
310
311
312
        }
    }
}


template<class ParcelType>
313
template<class TrackData>
314
315
void Foam::ReactingParcel<ParcelType>::calcMassTransfer
(
316
    TrackData& td,
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
    const scalar dt,
    const scalar T0,
    const scalar T1,
    scalarList& dMassMT
)
{
    if (td.cloud().composition().YMixture0()[1]>SMALL)
    {
        notImplemented
        (
            "void Foam::ReactingParcel<ParcelType>::"
            "calcMassTransfer(...): no treatment currently "
            "available for particles containing liquid species"
        )
    }

    // Check that model is active, and that the parcel temperature is
    // within necessary limits for mass transfer to occur
    if
    (
        !td.cloud().massTransfer().active()
338
339
     || this->T_<td.constProps().Tvap()
     || this->T_<td.constProps().Tbp()
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
    )
    {
        return;
    }

    // Determine mass to add to carrier phase
    const scalar mass = this->mass();
    const scalar dMassTot = td.cloud().massTransfer().calculate
    (
        dt,
        mass0_,
        mass,
        td.cloud().composition().YMixture0(),
        YMixture_,
        T0,
        canCombust_
    );

    // Update (total) mass fractions
    YMixture_[0] = (YMixture_[0]*mass - dMassTot)/(mass - dMassTot);
    YMixture_[1] = YMixture_[1]*mass/(mass - dMassTot);
    YMixture_[2] = 1.0 - YMixture_[0] - YMixture_[1];

    // Add to cummulative mass transfer
    forAll (YGas_, i)
    {
        label id = td.cloud().composition().gasGlobalIds()[i];

        // Mass transfer
        scalar volatileMass = YGas_[i]*dMassTot;
        dMassMT[id] += volatileMass;
    }
}


template<class ParcelType>
376
template<class TrackData>
377
378
void Foam::ReactingParcel<ParcelType>::calcSurfaceReactions
(
379
    TrackData& td,
380
381
382
383
    const scalar dt,
    const label celli,
    const scalar T0,
    const scalar T1,
384
385
386
    const scalarList& dMassMT,
    scalarList& dMassSR,
    scalar& dhRet
387
388
389
390
391
392
393
394
395
396
397
398
399
400
)
{
    // Check that model is active
    if (!td.cloud().surfaceReaction().active() || !canCombust_)
    {
        return;
    }

    // Update mass transfer(s)
    // - Also updates Y()'s
    td.cloud().surfaceReaction().calculate
    (
        dt,
        celli,
401
        this->d_,
402
403
        T0,
        T1,
404
        this->Tc_,
405
        pc_,
406
        this->rhoc_,
407
        this->mass(),
408
        dMassMT,
409
410
411
412
        YGas_,
        YLiquid_,
        YSolid_,
        YMixture_,
413
414
        dMassSR,
        dhRet
415
416
417
418
419
420
421
422
423
424
    );
}


// * * * * * * * * * * * * * * * *  IOStream operators * * * * * * * * * * * //

#include "ReactingParcelIO.C"

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