applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H

   1     Info<< "Reading transportProperties\n" << endl;
   2
   3     IOdictionary transportProperties
   4     (
   5         IOobject
   6         (
   7             "transportProperties",
   8             runTime.constant(),
   9             mesh,
  10             IOobject::MUST_READ,
  11             IOobject::NO_WRITE
  12         )
  13     );
  14
  15     autoPtr<phaseModel> phasea = phaseModel::New
  16     (
  17         mesh,
  18         transportProperties,
  19         "a"
  20     );
  21
  22     autoPtr<phaseModel> phaseb = phaseModel::New
  23     (
  24         mesh,
  25         transportProperties,
  26         "b"
  27     );
  28
  29     volVectorField& Ua = phasea->U();
  30     surfaceScalarField& phia = phasea->phi();
  31     const dimensionedScalar& rhoa = phasea->rho();
  32     const dimensionedScalar& nua = phasea->nu();
  33
  34     volVectorField& Ub = phaseb->U();
  35     surfaceScalarField& phib = phaseb->phi();
  36     const dimensionedScalar& rhob = phaseb->rho();
  37     const dimensionedScalar& nub = phaseb->nu();
  38
  39     Info<< "Reading field alpha\n" << endl;
  40     volScalarField alpha
  41     (
  42         IOobject
  43         (
  44             "alpha",
  45             runTime.timeName(),
  46             mesh,
  47             IOobject::MUST_READ,
  48             IOobject::AUTO_WRITE
  49         ),
  50         mesh
  51     );
  52
  53     volScalarField beta
  54     (
  55         IOobject
  56         (
  57             "beta",
  58             runTime.timeName(),
  59             mesh,
  60             IOobject::NO_READ,
  61             IOobject::NO_WRITE
  62         ),
  63         scalar(1) - alpha
  64         //,alpha.boundaryField().types()
  65     );
  66
  67     Info<< "Reading field p\n" << endl;
  68     volScalarField p
  69     (
  70         IOobject
  71         (
  72             "p",
  73             runTime.timeName(),
  74             mesh,
  75             IOobject::MUST_READ,
  76             IOobject::AUTO_WRITE
  77         ),
  78         mesh
  79     );
  80
  81     volVectorField U
  82     (
  83         IOobject
  84         (
  85             "U",
  86             runTime.timeName(),
  87             mesh,
  88             IOobject::NO_READ,
  89             IOobject::AUTO_WRITE
  90         ),
  91         alpha*Ua + beta*Ub
  92     );
  93
  94     Info<< "Reading field k\n" << endl;
  95     volScalarField k
  96     (
  97         IOobject
  98         (
  99             "k",
 100             runTime.timeName(),
 101             mesh,
 102             IOobject::MUST_READ,
 103             IOobject::AUTO_WRITE
 104         ),
 105         mesh
 106     );
 107
 108     Info<< "Reading field epsilon\n" << endl;
 109     volScalarField epsilon
 110     (
 111         IOobject
 112         (
 113             "epsilon",
 114             runTime.timeName(),
 115             mesh,
 116             IOobject::MUST_READ,
 117             IOobject::AUTO_WRITE
 118         ),
 119         mesh
 120     );
 121
 122     dimensionedScalar Cvm
 123     (
 124         transportProperties.lookup("Cvm")
 125     );
 126
 127     dimensionedScalar Cl
 128     (
 129         transportProperties.lookup("Cl")
 130     );
 131
 132     dimensionedScalar Ct
 133     (
 134         transportProperties.lookup("Ct")
 135     );
 136
 137     surfaceScalarField phi
 138     (
 139         IOobject
 140         (
 141             "phi",
 142             runTime.timeName(),
 143             mesh
 144         ),
 145         fvc::interpolate(alpha)*phia + fvc::interpolate(beta)*phib
 146     );
 147
 148
 149     IOdictionary RASProperties
 150     (
 151         IOobject
 152         (
 153             "RASProperties",
 154             runTime.constant(),
 155             mesh,
 156             IOobject::MUST_READ,
 157             IOobject::NO_WRITE
 158         )
 159     );
 160
 161
 162     Switch turbulence
 163     (
 164         RASProperties.lookup("turbulence")
 165     );
 166
 167     dictionary kEpsilonCoeffs
 168     (
 169         RASProperties.subDict("kEpsilonCoeffs")
 170     );
 171
 172
 173     scalar Cmu
 174     (
 175         readScalar(kEpsilonCoeffs.lookup("Cmu"))
 176     );
 177
 178     scalar C1
 179     (
 180         readScalar(kEpsilonCoeffs.lookup("C1"))
 181     );
 182
 183     scalar C2
 184     (
 185         readScalar(kEpsilonCoeffs.lookup("C2"))
 186     );
 187
 188     scalar alphak
 189     (
 190         readScalar(kEpsilonCoeffs.lookup("alphak"))
 191     );
 192
 193     scalar alphaEps
 194     (
 195         readScalar(kEpsilonCoeffs.lookup("alphaEps"))
 196     );
 197
 198     dictionary wallFunctionCoeffs
 199     (
 200         RASProperties.subDict("wallFunctionCoeffs")
 201     );
 202
 203     scalar kappa
 204     (
 205         readScalar(wallFunctionCoeffs.lookup("kappa"))
 206     );
 207
 208     scalar E
 209     (
 210         readScalar(wallFunctionCoeffs.lookup("E"))
 211     );
 212
 213     nearWallDist y(mesh);
 214
 215     Info<< "Calculating field nutb\n" << endl;
 216     volScalarField nutb
 217     (
 218         IOobject
 219         (
 220             "nutb",
 221             runTime.timeName(),
 222             mesh,
 223             IOobject::NO_READ,
 224             IOobject::AUTO_WRITE
 225         ),
 226         Cmu*sqr(k)/epsilon
 227     );
 228
 229     Info<< "Calculating field nuEffa\n" << endl;
 230     volScalarField nuEffa
 231     (
 232         IOobject
 233         (
 234             "nuEffa",
 235             runTime.timeName(),
 236             mesh,
 237             IOobject::NO_READ,
 238             IOobject::NO_WRITE
 239         ),
 240         sqr(Ct)*nutb + nua
 241     );
 242
 243     Info<< "Calculating field nuEffb\n" << endl;
 244     volScalarField nuEffb
 245     (
 246         IOobject
 247         (
 248             "nuEffb",
 249             runTime.timeName(),
 250             mesh,
 251             IOobject::NO_READ,
 252             IOobject::NO_WRITE
 253         ),
 254         nutb + nub
 255     );
 256
 257
 258     Info<< "Calculating field DDtUa and DDtUb\n" << endl;
 259
 260     volVectorField DDtUa =
 261         fvc::ddt(Ua)
 262       + fvc::div(phia, Ua)
 263       - fvc::div(phia)*Ua;
 264
 265     volVectorField DDtUb =
 266         fvc::ddt(Ub)
 267       + fvc::div(phib, Ub)
 268       - fvc::div(phib)*Ub;
 269
 270
 271     Info<< "Calculating field g.h\n" << endl;
 272     volScalarField gh("gh", g & mesh.C());
 273
 274     IOdictionary interfacialProperties
 275     (
 276         IOobject
 277         (
 278             "interfacialProperties",
 279             runTime.constant(),
 280             mesh,
 281             IOobject::MUST_READ,
 282             IOobject::NO_WRITE
 283         )
 284     );
 285
 286     autoPtr<dragModel> draga = dragModel::New
 287     (
 288         interfacialProperties,
 289         alpha,
 290         phasea,
 291         phaseb
 292     );
 293
 294     autoPtr<dragModel> dragb = dragModel::New
 295     (
 296         interfacialProperties,
 297         beta,
 298         phaseb,
 299         phasea
 300     );
 301
 302     word dragPhase("blended");
 303     if (interfacialProperties.found("dragPhase"))
 304     {
 305         dragPhase = word(interfacialProperties.lookup("dragPhase"));
 306
 307         bool validDrag =
 308             dragPhase == "a" || dragPhase == "b" || dragPhase == "blended";
 309
 310         if (!validDrag)
 311         {
 312             FatalErrorIn(args.executable())
 313                 << "invalid dragPhase " << dragPhase
 314                 << exit(FatalError);
 315         }
 316     }
 317
 318     Info << "dragPhase is " << dragPhase << endl;
 319     kineticTheoryModel kineticTheory
 320     (
 321         phasea,
 322         Ub,
 323         alpha,
 324         draga
 325     );
 326
 327     surfaceScalarField rUaAf
 328     (
 329         IOobject
 330         (
 331             "rUaAf",
 332             runTime.timeName(),
 333             mesh,
 334             IOobject::NO_READ,
 335             IOobject::NO_WRITE
 336         ),
 337         mesh,
 338         dimensionedScalar("zero", dimensionSet(0, 0, 1, 0, 0), 0.0)
 339     );
 340
 341     surfaceScalarField ppMagf
 342     (
 343         IOobject
 344         (
 345             "ppMagf",
 346             runTime.timeName(),
 347             mesh,
 348             IOobject::NO_READ,
 349             IOobject::NO_WRITE
 350         ),
 351         mesh,
 352         dimensionedScalar("zero", dimensionSet(0, 2, -1, 0, 0), 0.0)
 353     );
 354
 355
 356     label pRefCell = 0;
 357     scalar pRefValue = 0.0;
 358     setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);