applications/solvers/multiphase/settlingFoam/createFields.H

   1     Info<< "Reading field p_rgh\n" << endl;
   2     volScalarField p_rgh
   3     (
   4         IOobject
   5         (
   6             "p_rgh",
   7             runTime.timeName(),
   8             mesh,
   9             IOobject::MUST_READ,
  10             IOobject::AUTO_WRITE
  11         ),
  12         mesh
  13     );
  14
  15     Info<< "Reading field alpha\n" << endl;
  16     volScalarField alpha
  17     (
  18         IOobject
  19         (
  20             "alpha",
  21             runTime.timeName(),
  22             mesh,
  23             IOobject::MUST_READ,
  24             IOobject::AUTO_WRITE
  25         ),
  26         mesh
  27     );
  28
  29     Info<< "Reading field U\n" << endl;
  30     volVectorField U
  31     (
  32         IOobject
  33         (
  34             "U",
  35             runTime.timeName(),
  36             mesh,
  37             IOobject::MUST_READ,
  38             IOobject::AUTO_WRITE
  39         ),
  40         mesh
  41     );
  42
  43
  44     Info<< "Reading transportProperties\n" << endl;
  45
  46     IOdictionary transportProperties
  47     (
  48         IOobject
  49         (
  50             "transportProperties",
  51             runTime.constant(),
  52             mesh,
  53             IOobject::MUST_READ,
  54             IOobject::NO_WRITE
  55         )
  56     );
  57
  58
  59     dimensionedScalar rhoc
  60     (
  61         transportProperties.lookup("rhoc")
  62     );
  63
  64     dimensionedScalar rhod
  65     (
  66         transportProperties.lookup("rhod")
  67     );
  68
  69     dimensionedScalar muc
  70     (
  71         transportProperties.lookup("muc")
  72     );
  73
  74     dimensionedScalar plasticViscosityCoeff
  75     (
  76         transportProperties.lookup("plasticViscosityCoeff")
  77     );
  78
  79     dimensionedScalar plasticViscosityExponent
  80     (
  81         transportProperties.lookup("plasticViscosityExponent")
  82     );
  83
  84     dimensionedScalar yieldStressCoeff
  85     (
  86         transportProperties.lookup("yieldStressCoeff")
  87     );
  88
  89     dimensionedScalar yieldStressExponent
  90     (
  91         transportProperties.lookup("yieldStressExponent")
  92     );
  93
  94     dimensionedScalar yieldStressOffset
  95     (
  96         transportProperties.lookup("yieldStressOffset")
  97     );
  98
  99     Switch BinghamPlastic
 100     (
 101         transportProperties.lookup("BinghamPlastic")
 102     );
 103
 104     volScalarField rho
 105     (
 106         IOobject
 107         (
 108             "rho",
 109             runTime.timeName(),
 110             mesh,
 111             IOobject::NO_READ,
 112             IOobject::NO_WRITE
 113         ),
 114         (scalar(1) - alpha)*rhoc + alpha*rhod
 115     );
 116
 117     volScalarField Alpha
 118     (
 119         IOobject
 120         (
 121             "Alpha",
 122             runTime.timeName(),
 123             mesh,
 124             IOobject::NO_READ,
 125             IOobject::AUTO_WRITE
 126         ),
 127         alpha*rhod/rho,
 128         alpha.boundaryField().types()
 129     );
 130
 131     #include "compressibleCreatePhi.H"
 132
 133
 134     Info<< "Calculating field mul\n" << endl;
 135     volScalarField mul
 136     (
 137         IOobject
 138         (
 139             "mul",
 140             runTime.timeName(),
 141             mesh,
 142             IOobject::NO_READ,
 143             IOobject::AUTO_WRITE
 144         ),
 145         muc +
 146         plasticViscosity
 147         (
 148             plasticViscosityCoeff,
 149             plasticViscosityExponent,
 150             Alpha
 151         )
 152     );
 153
 154
 155     Info<< "Initialising field Vdj\n" << endl;
 156     volVectorField Vdj
 157     (
 158         IOobject
 159         (
 160             "Vdj",
 161             runTime.timeName(),
 162             mesh,
 163             IOobject::NO_READ,
 164             IOobject::AUTO_WRITE
 165         ),
 166         mesh,
 167         dimensionedVector("0.0", U.dimensions(), vector::zero),
 168         U.boundaryField().types()
 169     );
 170
 171
 172     Info<< "Selecting Drift-Flux model " << endl;
 173
 174     word VdjModel
 175     (
 176         transportProperties.lookup("VdjModel")
 177     );
 178
 179     Info<< tab << VdjModel << " selected\n" << endl;
 180
 181     const dictionary& VdjModelCoeffs
 182     (
 183         transportProperties.subDict(VdjModel + "Coeffs")
 184     );
 185
 186     dimensionedVector V0
 187     (
 188         VdjModelCoeffs.lookup("V0")
 189     );
 190
 191     dimensionedScalar a
 192     (
 193         VdjModelCoeffs.lookup("a")
 194     );
 195
 196     dimensionedScalar a1
 197     (
 198         VdjModelCoeffs.lookup("a1")
 199     );
 200
 201     dimensionedScalar alphaMin
 202     (
 203         VdjModelCoeffs.lookup("alphaMin")
 204     );
 205
 206
 207
 208     IOdictionary RASProperties
 209     (
 210         IOobject
 211         (
 212             "RASProperties",
 213             runTime.constant(),
 214             mesh,
 215             IOobject::MUST_READ,
 216             IOobject::NO_WRITE
 217         )
 218     );
 219
 220
 221     Switch turbulence
 222     (
 223         RASProperties.lookup("turbulence")
 224     );
 225
 226     dictionary kEpsilonDict
 227     (
 228         RASProperties.subDictPtr("kEpsilonCoeffs")
 229     );
 230
 231     dimensionedScalar Cmu
 232     (
 233         dimensionedScalar::lookupOrAddToDict
 234         (
 235             "Cmu",
 236             kEpsilonDict,
 237             0.09
 238         )
 239     );
 240
 241     dimensionedScalar C1
 242     (
 243         dimensionedScalar::lookupOrAddToDict
 244         (
 245             "C1",
 246             kEpsilonDict,
 247             1.44
 248         )
 249     );
 250
 251     dimensionedScalar C2
 252     (
 253         dimensionedScalar::lookupOrAddToDict
 254         (
 255             "C2",
 256             kEpsilonDict,
 257             1.92
 258         )
 259     );
 260
 261     dimensionedScalar C3
 262     (
 263         dimensionedScalar::lookupOrAddToDict
 264         (
 265             "C3",
 266             kEpsilonDict,
 267             0.85
 268         )
 269     );
 270
 271     dimensionedScalar alphak
 272     (
 273         dimensionedScalar::lookupOrAddToDict
 274         (
 275             "alphaEps",
 276             kEpsilonDict,
 277             1.0
 278         )
 279     );
 280
 281     dimensionedScalar alphaEps
 282     (
 283         dimensionedScalar::lookupOrAddToDict
 284         (
 285             "alphaEps",
 286             kEpsilonDict,
 287             0.76923
 288         )
 289     );
 290
 291     dictionary wallFunctionDict
 292     (
 293         RASProperties.subDictPtr("wallFunctionCoeffs")
 294     );
 295
 296     dimensionedScalar kappa
 297     (
 298         dimensionedScalar::lookupOrAddToDict
 299         (
 300             "kappa",
 301             wallFunctionDict,
 302             0.41
 303         )
 304     );
 305
 306     dimensionedScalar E
 307     (
 308         dimensionedScalar::lookupOrAddToDict
 309         (
 310             "E",
 311             wallFunctionDict,
 312             9.0
 313         )
 314     );
 315
 316     if (RASProperties.lookupOrDefault<Switch>("printCoeffs", false))
 317     {
 318         Info<< "kEpsilonCoeffs" << kEpsilonDict << nl
 319             << "wallFunctionCoeffs" << wallFunctionDict << endl;
 320     }
 321
 322
 323     nearWallDist y(mesh);
 324
 325     Info<< "Reading field k\n" << endl;
 326     volScalarField k
 327     (
 328         IOobject
 329         (
 330             "k",
 331             runTime.timeName(),
 332             mesh,
 333             IOobject::MUST_READ,
 334             IOobject::AUTO_WRITE
 335         ),
 336         mesh
 337     );
 338
 339     Info<< "Reading field epsilon\n" << endl;
 340     volScalarField epsilon
 341     (
 342         IOobject
 343         (
 344             "epsilon",
 345             runTime.timeName(),
 346             mesh,
 347             IOobject::MUST_READ,
 348             IOobject::AUTO_WRITE
 349         ),
 350         mesh
 351     );
 352
 353     Info<< "Calculating field mut\n" << endl;
 354     volScalarField mut
 355     (
 356         IOobject
 357         (
 358             "mut",
 359             runTime.timeName(),
 360             mesh,
 361             IOobject::NO_READ,
 362             IOobject::AUTO_WRITE
 363         ),
 364         Cmu*rho*sqr(k)/epsilon
 365     );
 366
 367
 368     Info<< "Calculating field mu\n" << endl;
 369     volScalarField mu
 370     (
 371         IOobject
 372         (
 373             "mu",
 374             runTime.timeName(),
 375             mesh,
 376             IOobject::NO_READ,
 377             IOobject::AUTO_WRITE
 378         ),
 379         mut + mul
 380     );
 381
 382
 383     Info<< "Calculating field (g.h)f\n" << endl;
 384     volScalarField gh("gh", g & mesh.C());
 385     surfaceScalarField ghf = surfaceScalarField("gh", g & mesh.Cf());
 386
 387     volScalarField p
 388     (
 389         IOobject
 390         (
 391             "p",
 392             runTime.timeName(),
 393             mesh,
 394             IOobject::NO_READ,
 395             IOobject::AUTO_WRITE
 396         ),
 397         p_rgh + rho*gh
 398     );
 399
 400     label pRefCell = 0;
 401     scalar pRefValue = 0.0;
 402     setRefCell
 403     (
 404         p,
 405         p_rgh,
 406         mesh.solutionDict().subDict("PISO"),
 407         pRefCell,
 408         pRefValue
 409     );
 410
 411     if (p_rgh.needReference())
 412     {
 413         p += dimensionedScalar
 414         (
 415             "p",
 416             p.dimensions(),
 417             pRefValue - getRefCellValue(p, pRefCell)
 418         );
 419         p_rgh = p - rho*gh;
 420     }