src/thermophysicalModels/liquidMixture/liquidMixture/liquidMixture.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
   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 the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     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, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 \*---------------------------------------------------------------------------*/
  26
  27 #include "liquidMixture.H"
  28 #include "dictionary.H"
  29 #include "specie.H"
  30
  31 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  32
  33 const Foam::scalar Foam::liquidMixture::TrMax = 0.999;
  34
  35 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  36
  37 Foam::liquidMixture::liquidMixture
  38 (
  39     const dictionary& thermophysicalProperties
  40 )
  41 :
  42     components_(thermophysicalProperties.lookup("liquidComponents")),
  43     properties_(components_.size())
  44 {
  45     // use sub-dictionary "liquidProperties" if possible to avoid
  46     // collisions with identically named gas-phase entries
  47     // (eg, H2O liquid vs. gas)
  48     forAll(components_, i)
  49     {
  50         const dictionary* subDictPtr = thermophysicalProperties.subDictPtr
  51         (
  52             "liquidProperties"
  53         );
  54
  55         if (subDictPtr)
  56         {
  57             properties_.set
  58             (
  59                 i,
  60                 liquid::New(subDictPtr->lookup(components_[i]))
  61             );
  62         }
  63         else
  64         {
  65             properties_.set
  66             (
  67                 i,
  68                 liquid::New(thermophysicalProperties.lookup(components_[i]))
  69             );
  70         }
  71     }
  72 }
  73
  74
  75 // * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
  76
  77 Foam::autoPtr<Foam::liquidMixture> Foam::liquidMixture::New
  78 (
  79     const dictionary& thermophysicalProperties
  80 )
  81 {
  82     return autoPtr<liquidMixture>(new liquidMixture(thermophysicalProperties));
  83 }
  84
  85
  86 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  87
  88 Foam::scalar Foam::liquidMixture::Tc
  89 (
  90     const scalarField& x
  91 ) const
  92 {
  93
  94     scalar vTc = 0.0;
  95     scalar vc = 0.0;
  96
  97     forAll(properties_, i)
  98     {
  99         scalar x1 = x[i]*properties_[i].Vc();
 100         vc += x1;
 101         vTc += x1*properties_[i].Tc();
 102     }
 103
 104     return vTc/vc;
 105 }
 106
 107
 108 Foam::scalar Foam::liquidMixture::Tpc
 109 (
 110     const scalarField& x
 111 ) const
 112 {
 113     scalar Tpc = 0.0;
 114     forAll(properties_, i)
 115     {
 116         Tpc += x[i]*properties_[i].Tc();
 117     }
 118
 119     return Tpc;
 120 }
 121
 122
 123 Foam::scalar Foam::liquidMixture::Ppc
 124 (
 125     const scalarField& x
 126 ) const
 127 {
 128     scalar Vc = 0.0;
 129     scalar Zc = 0.0;
 130     forAll(properties_, i)
 131     {
 132         Vc += x[i]*properties_[i].Vc();
 133         Zc += x[i]*properties_[i].Zc();
 134     }
 135
 136     return specie::RR*Zc*Tpc(x)/Vc;
 137 }
 138
 139
 140 Foam::scalar Foam::liquidMixture::omega
 141 (
 142     const scalarField& x
 143 ) const
 144 {
 145     scalar omega = 0.0;
 146     forAll(properties_, i)
 147     {
 148         omega += x[i]*properties_[i].omega();
 149     }
 150
 151     return omega;
 152 }
 153
 154
 155 Foam::scalarField Foam::liquidMixture::Xs
 156 (
 157     const scalar p,
 158     const scalar Tg,
 159     const scalar Tl,
 160     const scalarField& xg,
 161     const scalarField& xl
 162 ) const
 163 {
 164     scalarField xs(xl.size(), 0.0);
 165
 166     // Raoult's Law
 167     forAll(xs, i)
 168     {
 169         scalar Ti = min(TrMax*properties_[i].Tc(), Tl);
 170         xs[i] = properties_[i].pv(p, Ti)*xl[i]/p;
 171     }
 172     return xs;
 173 }
 174
 175
 176 Foam::scalar Foam::liquidMixture::W
 177 (
 178     const scalarField& x
 179 ) const
 180 {
 181     scalar W = 0.0;
 182     forAll(properties_, i)
 183     {
 184         W += x[i]*properties_[i].W();
 185     }
 186
 187     return W;
 188 }
 189
 190
 191 Foam::scalarField Foam::liquidMixture::Y
 192 (
 193     const scalarField& X
 194 ) const
 195 {
 196     scalarField Y(X/W(X));
 197
 198     forAll(Y, i)
 199     {
 200         Y[i] *= properties_[i].W();
 201     }
 202
 203     return Y;
 204 }
 205
 206
 207 Foam::scalarField Foam::liquidMixture::X
 208 (
 209     const scalarField& Y
 210 ) const
 211 {
 212     scalarField X(Y.size());
 213     scalar Winv = 0.0;
 214     forAll(X, i)
 215     {
 216         Winv += Y[i]/properties_[i].W();
 217         X[i] = Y[i]/properties_[i].W();
 218     }
 219
 220     return X/Winv;
 221 }
 222
 223
 224 Foam::scalar Foam::liquidMixture::rho
 225 (
 226     const scalar p,
 227     const scalar T,
 228     const scalarField& x
 229 ) const
 230 {
 231     scalar v = 0.0;
 232
 233     forAll(properties_, i)
 234     {
 235         if (x[i] > SMALL)
 236         {
 237             scalar Ti = min(TrMax*properties_[i].Tc(), T);
 238             scalar rho = SMALL + properties_[i].rho(p, Ti);
 239             v += x[i]*properties_[i].W()/rho;
 240         }
 241     }
 242
 243     return W(x)/v;
 244 }
 245
 246
 247 Foam::scalar Foam::liquidMixture::pv
 248 (
 249     const scalar p,
 250     const scalar T,
 251     const scalarField& x
 252 ) const
 253 {
 254     scalar pv = 0.0;
 255
 256     forAll(properties_, i)
 257     {
 258         if (x[i] > SMALL)
 259         {
 260             scalar Ti = min(TrMax*properties_[i].Tc(), T);
 261             pv += x[i]*properties_[i].pv(p, Ti)*properties_[i].W();
 262         }
 263     }
 264
 265     return pv/W(x);
 266 }
 267
 268
 269 Foam::scalar Foam::liquidMixture::hl
 270 (
 271     const scalar p,
 272     const scalar T,
 273     const scalarField& x
 274 ) const
 275 {
 276     scalar hl = 0.0;
 277
 278     forAll(properties_, i)
 279     {
 280         if (x[i] > SMALL)
 281         {
 282             scalar Ti = min(TrMax*properties_[i].Tc(), T);
 283             hl += x[i]*properties_[i].hl(p, Ti)*properties_[i].W();
 284         }
 285     }
 286
 287     return hl/W(x);
 288 }
 289
 290
 291 Foam::scalar Foam::liquidMixture::cp
 292 (
 293     const scalar p,
 294     const scalar T,
 295     const scalarField& x
 296 ) const
 297 {
 298     scalar cp = 0.0;
 299
 300     forAll(properties_, i)
 301     {
 302         if (x[i] > SMALL)
 303         {
 304             scalar Ti = min(TrMax*properties_[i].Tc(), T);
 305             cp += x[i]*properties_[i].cp(p, Ti)*properties_[i].W();
 306         }
 307     }
 308
 309     return cp/W(x);
 310 }
 311
 312
 313 Foam::scalar Foam::liquidMixture::sigma
 314 (
 315     const scalar p,
 316     const scalar T,
 317     const scalarField& x
 318 ) const
 319 {
 320     // sigma is based on surface mole fractions
 321     // which is estimated from Raoult's Law
 322     scalar sigma = 0.0;
 323     scalarField Xs(x.size(), 0.0);
 324     scalar XsSum = 0.0;
 325     forAll(properties_, i)
 326     {
 327         scalar Ti = min(TrMax*properties_[i].Tc(), T);
 328         scalar Pvs = properties_[i].pv(p, Ti);
 329         scalar xs = x[i]*Pvs/p;
 330         XsSum += xs;
 331         Xs[i] = xs;
 332     }
 333
 334     forAll(properties_, i)
 335     {
 336         if (Xs[i] > SMALL)
 337         {
 338             scalar Ti = min(TrMax*properties_[i].Tc(), T);
 339             sigma += (Xs[i]/XsSum)*properties_[i].sigma(p, Ti);
 340         }
 341     }
 342
 343     return sigma;
 344 }
 345
 346
 347 Foam::scalar Foam::liquidMixture::mu
 348 (
 349     const scalar p,
 350     const scalar T,
 351     const scalarField& x
 352 ) const
 353 {
 354     scalar mu = 0.0;
 355
 356     forAll(properties_, i)
 357     {
 358         if (x[i] > SMALL)
 359         {
 360             scalar Ti = min(TrMax*properties_[i].Tc(), T);
 361             mu += x[i]*log(properties_[i].mu(p, Ti));
 362         }
 363     }
 364
 365     return exp(mu);
 366 }
 367
 368
 369 Foam::scalar Foam::liquidMixture::K
 370 (
 371     const scalar p,
 372     const scalar T,
 373     const scalarField& x
 374 ) const
 375 {
 376     // calculate superficial volume fractions phii
 377     scalarField phii(x.size(), 0.0);
 378     scalar pSum = 0.0;
 379
 380     forAll(properties_, i)
 381     {
 382         scalar Ti = min(TrMax*properties_[i].Tc(), T);
 383
 384         scalar Vi = properties_[i].W()/properties_[i].rho(p, Ti);
 385         phii[i] = x[i]*Vi;
 386         pSum += phii[i];
 387     }
 388
 389     forAll(phii, i)
 390     {
 391         phii[i] /= pSum;
 392     }
 393
 394     scalar K = 0.0;
 395
 396     forAll(properties_, i)
 397     {
 398         scalar Ti = min(TrMax*properties_[i].Tc(), T);
 399
 400         forAll(properties_, j)
 401         {
 402             scalar Tj = min(TrMax*properties_[j].Tc(), T);
 403
 404             scalar Kij =
 405                 2.0
 406                /(
 407                     1.0/properties_[i].K(p, Ti)
 408                   + 1.0/properties_[j].K(p, Tj)
 409                 );
 410             K += phii[i]*phii[j]*Kij;
 411         }
 412     }
 413
 414     return K;
 415 }
 416
 417
 418 Foam::scalar Foam::liquidMixture::D
 419 (
 420     const scalar p,
 421     const scalar T,
 422     const scalarField& x
 423 ) const
 424 {
 425     // Blanc's law
 426     scalar Dinv = 0.0;
 427
 428     forAll(properties_, i)
 429     {
 430         if (x[i] > SMALL)
 431         {
 432             scalar Ti = min(TrMax*properties_[i].Tc(), T);
 433             Dinv += x[i]/properties_[i].D(p, Ti);
 434         }
 435     }
 436
 437     return 1.0/Dinv;
 438 }
 439
 440
 441 // ************************************************************************* //