applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2004-2010 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
  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 "kineticTheoryModel.H"
  27 #include "surfaceInterpolate.H"
  28 #include "mathematicalConstants.H"
  29 #include "fvCFD.H"
  30
  31 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  32
  33 Foam::kineticTheoryModel::kineticTheoryModel
  34 (
  35     const Foam::phaseModel& phasea,
  36     const Foam::volVectorField& Ub,
  37     const Foam::volScalarField& alpha,
  38     const Foam::dragModel& draga
  39 )
  40 :
  41     phasea_(phasea),
  42     Ua_(phasea.U()),
  43     Ub_(Ub),
  44     alpha_(alpha),
  45     phia_(phasea.phi()),
  46     draga_(draga),
  47
  48     rhoa_(phasea.rho()),
  49     da_(phasea.d()),
  50     nua_(phasea.nu()),
  51
  52     kineticTheoryProperties_
  53     (
  54         IOobject
  55         (
  56             "kineticTheoryProperties",
  57             Ua_.time().constant(),
  58             Ua_.mesh(),
  59             IOobject::MUST_READ,
  60             IOobject::NO_WRITE
  61         )
  62     ),
  63     kineticTheory_(kineticTheoryProperties_.lookup("kineticTheory")),
  64     equilibrium_(kineticTheoryProperties_.lookup("equilibrium")),
  65
  66     viscosityModel_
  67     (
  68         kineticTheoryModels::viscosityModel::New
  69         (
  70             kineticTheoryProperties_
  71         )
  72     ),
  73     conductivityModel_
  74     (
  75         conductivityModel::New
  76         (
  77             kineticTheoryProperties_
  78         )
  79     ),
  80     radialModel_
  81     (
  82         radialModel::New
  83         (
  84             kineticTheoryProperties_
  85         )
  86     ),
  87     granularPressureModel_
  88     (
  89         granularPressureModel::New
  90         (
  91             kineticTheoryProperties_
  92         )
  93     ),
  94     frictionalStressModel_
  95     (
  96         frictionalStressModel::New
  97         (
  98             kineticTheoryProperties_
  99         )
 100     ),
 101     e_(kineticTheoryProperties_.lookup("e")),
 102     alphaMax_(kineticTheoryProperties_.lookup("alphaMax")),
 103     alphaMinFriction_(kineticTheoryProperties_.lookup("alphaMinFriction")),
 104     Fr_(kineticTheoryProperties_.lookup("Fr")),
 105     eta_(kineticTheoryProperties_.lookup("eta")),
 106     p_(kineticTheoryProperties_.lookup("p")),
 107     phi_(dimensionedScalar(kineticTheoryProperties_.lookup("phi"))*M_PI/180.0),
 108     Theta_
 109     (
 110         IOobject
 111         (
 112             "Theta",
 113             Ua_.time().timeName(),
 114             Ua_.mesh(),
 115             IOobject::MUST_READ,
 116             IOobject::AUTO_WRITE
 117         ),
 118         Ua_.mesh()
 119     ),
 120     mua_
 121     (
 122         IOobject
 123         (
 124             "mua",
 125             Ua_.time().timeName(),
 126             Ua_.mesh(),
 127             IOobject::NO_READ,
 128             IOobject::AUTO_WRITE
 129         ),
 130         Ua_.mesh(),
 131         dimensionedScalar("zero", dimensionSet(1, -1, -1, 0, 0), 0.0)
 132     ),
 133     lambda_
 134     (
 135         IOobject
 136         (
 137             "lambda",
 138             Ua_.time().timeName(),
 139             Ua_.mesh(),
 140             IOobject::NO_READ,
 141             IOobject::NO_WRITE
 142         ),
 143         Ua_.mesh(),
 144         dimensionedScalar("zero", dimensionSet(1, -1, -1, 0, 0), 0.0)
 145     ),
 146     pa_
 147     (
 148         IOobject
 149         (
 150             "pa",
 151             Ua_.time().timeName(),
 152             Ua_.mesh(),
 153             IOobject::NO_READ,
 154             IOobject::AUTO_WRITE
 155         ),
 156         Ua_.mesh(),
 157         dimensionedScalar("zero", dimensionSet(1, -1, -2, 0, 0), 0.0)
 158     ),
 159     kappa_
 160     (
 161         IOobject
 162         (
 163             "kappa",
 164             Ua_.time().timeName(),
 165             Ua_.mesh(),
 166             IOobject::NO_READ,
 167             IOobject::NO_WRITE
 168         ),
 169         Ua_.mesh(),
 170         dimensionedScalar("zero", dimensionSet(1, -1, -1, 0, 0), 0.0)
 171     ),
 172     gs0_
 173     (
 174         IOobject
 175         (
 176             "gs0",
 177             Ua_.time().timeName(),
 178             Ua_.mesh(),
 179             IOobject::NO_READ,
 180             IOobject::NO_WRITE
 181         ),
 182         Ua_.mesh(),
 183         dimensionedScalar("zero", dimensionSet(0, 0, 0, 0, 0), 1.0)
 184     )
 185 {}
 186
 187
 188 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 189
 190 Foam::kineticTheoryModel::~kineticTheoryModel()
 191 {}
 192
 193
 194 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 195
 196 void Foam::kineticTheoryModel::solve(const volTensorField& gradUat)
 197 {
 198     if (!kineticTheory_)
 199     {
 200         return;
 201     }
 202
 203     const scalar sqrtPi = sqrt(constant::mathematical::pi);
 204
 205     surfaceScalarField phi(1.5*rhoa_*phia_*fvc::interpolate(alpha_));
 206
 207     volTensorField dU(gradUat.T());    //fvc::grad(Ua_);
 208     volSymmTensorField D(symm(dU));
 209
 210     // NB, drag = K*alpha*beta,
 211     // (the alpha and beta has been extracted from the drag function for
 212     // numerical reasons)
 213     volScalarField Ur(mag(Ua_ - Ub_));
 214     volScalarField betaPrim(alpha_*(1.0 - alpha_)*draga_.K(Ur));
 215
 216     // Calculating the radial distribution function (solid volume fraction is
 217     //  limited close to the packing limit, but this needs improvements)
 218     //  The solution is higly unstable close to the packing limit.
 219     gs0_ = radialModel_->g0
 220     (
 221         min(max(alpha_, scalar(1e-6)), alphaMax_ - 0.01),
 222         alphaMax_
 223     );
 224
 225     // particle pressure - coefficient in front of Theta (Eq. 3.22, p. 45)
 226     volScalarField PsCoeff
 227     (
 228         granularPressureModel_->granularPressureCoeff
 229         (
 230             alpha_,
 231             gs0_,
 232             rhoa_,
 233             e_
 234         )
 235     );
 236
 237     // 'thermal' conductivity (Table 3.3, p. 49)
 238     kappa_ = conductivityModel_->kappa(alpha_, Theta_, gs0_, rhoa_, da_, e_);
 239
 240     // particle viscosity (Table 3.2, p.47)
 241     mua_ = viscosityModel_->mua(alpha_, Theta_, gs0_, rhoa_, da_, e_);
 242
 243     dimensionedScalar Tsmall
 244     (
 245         "small",
 246         dimensionSet(0 , 2 ,-2 ,0 , 0, 0, 0),
 247         1.0e-6
 248     );
 249
 250     dimensionedScalar TsmallSqrt = sqrt(Tsmall);
 251     volScalarField ThetaSqrt(sqrt(Theta_));
 252
 253     // dissipation (Eq. 3.24, p.50)
 254     volScalarField gammaCoeff
 255     (
 256         12.0*(1.0 - sqr(e_))*sqr(alpha_)*rhoa_*gs0_*(1.0/da_)*ThetaSqrt/sqrtPi
 257     );
 258
 259     // Eq. 3.25, p. 50 Js = J1 - J2
 260     volScalarField J1(3.0*betaPrim);
 261     volScalarField J2
 262     (
 263         0.25*sqr(betaPrim)*da_*sqr(Ur)
 264        /(max(alpha_, scalar(1e-6))*rhoa_*sqrtPi*(ThetaSqrt + TsmallSqrt))
 265     );
 266
 267     // bulk viscosity  p. 45 (Lun et al. 1984).
 268     lambda_ = (4.0/3.0)*sqr(alpha_)*rhoa_*da_*gs0_*(1.0+e_)*ThetaSqrt/sqrtPi;
 269
 270     // stress tensor, Definitions, Table 3.1, p. 43
 271     volSymmTensorField tau(2.0*mua_*D + (lambda_ - (2.0/3.0)*mua_)*tr(D)*I);
 272
 273     if (!equilibrium_)
 274     {
 275         // construct the granular temperature equation (Eq. 3.20, p. 44)
 276         // NB. note that there are two typos in Eq. 3.20
 277         // no grad infront of Ps
 278         // wrong sign infront of laplacian
 279         fvScalarMatrix ThetaEqn
 280         (
 281             fvm::ddt(1.5*alpha_*rhoa_, Theta_)
 282           + fvm::div(phi, Theta_, "div(phi,Theta)")
 283          ==
 284             fvm::SuSp(-((PsCoeff*I) && dU), Theta_)
 285           + (tau && dU)
 286           + fvm::laplacian(kappa_, Theta_, "laplacian(kappa,Theta)")
 287           + fvm::Sp(-gammaCoeff, Theta_)
 288           + fvm::Sp(-J1, Theta_)
 289           + fvm::Sp(J2/(Theta_ + Tsmall), Theta_)
 290         );
 291
 292         ThetaEqn.relax();
 293         ThetaEqn.solve();
 294     }
 295     else
 296     {
 297         // equilibrium => dissipation == production
 298         // Eq. 4.14, p.82
 299         volScalarField K1(2.0*(1.0 + e_)*rhoa_*gs0_);
 300         volScalarField K3
 301         (
 302             0.5*da_*rhoa_*
 303             (
 304                 (sqrtPi/(3.0*(3.0-e_)))
 305                *(1.0 + 0.4*(1.0 + e_)*(3.0*e_ - 1.0)*alpha_*gs0_)
 306                +1.6*alpha_*gs0_*(1.0 + e_)/sqrtPi
 307             )
 308         );
 309
 310         volScalarField K2
 311         (
 312             4.0*da_*rhoa_*(1.0 + e_)*alpha_*gs0_/(3.0*sqrtPi) - 2.0*K3/3.0
 313         );
 314
 315         volScalarField K4(12.0*(1.0 - sqr(e_))*rhoa_*gs0_/(da_*sqrtPi));
 316
 317         volScalarField trD(tr(D));
 318         volScalarField tr2D(sqr(trD));
 319         volScalarField trD2(tr(D & D));
 320
 321         volScalarField t1(K1*alpha_ + rhoa_);
 322         volScalarField l1(-t1*trD);
 323         volScalarField l2(sqr(t1)*tr2D);
 324         volScalarField l3
 325         (
 326             4.0
 327            *K4
 328            *max(alpha_, scalar(1e-6))
 329            *(2.0*K3*trD2 + K2*tr2D)
 330         );
 331
 332         Theta_ = sqr((l1 + sqrt(l2 + l3))/(2.0*(alpha_ + 1.0e-4)*K4));
 333     }
 334
 335     Theta_.max(1.0e-15);
 336     Theta_.min(1.0e+3);
 337
 338     volScalarField pf
 339     (
 340         frictionalStressModel_->frictionalPressure
 341         (
 342             alpha_,
 343             alphaMinFriction_,
 344             alphaMax_,
 345             Fr_,
 346             eta_,
 347             p_
 348         )
 349     );
 350
 351     PsCoeff += pf/(Theta_+Tsmall);
 352
 353     PsCoeff.min(1.0e+10);
 354     PsCoeff.max(-1.0e+10);
 355
 356     // update particle pressure
 357     pa_ = PsCoeff*Theta_;
 358
 359     // frictional shear stress, Eq. 3.30, p. 52
 360     volScalarField muf
 361     (
 362         frictionalStressModel_->muf
 363         (
 364             alpha_,
 365             alphaMax_,
 366             pf,
 367             D,
 368             phi_
 369         )
 370     );
 371
 372     // add frictional stress
 373     mua_ += muf;
 374     mua_.min(1.0e+2);
 375     mua_.max(0.0);
 376
 377     Info<< "kinTheory: max(Theta) = " << max(Theta_).value() << endl;
 378
 379     volScalarField ktn(mua_/rhoa_);
 380
 381     Info<< "kinTheory: min(nua) = " << min(ktn).value()
 382         << ", max(nua) = " << max(ktn).value() << endl;
 383
 384     Info<< "kinTheory: min(pa) = " << min(pa_).value()
 385         << ", max(pa) = " << max(pa_).value() << endl;
 386 }
 387
 388
 389 // ************************************************************************* //