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