applications/solvers/multiphase/interMixingFoam/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 alpha1\n" << endl;
  16     volScalarField alpha1
  17     (
  18         IOobject
  19         (
  20             "alpha1",
  21             runTime.timeName(),
  22             mesh,
  23             IOobject::MUST_READ,
  24             IOobject::AUTO_WRITE
  25         ),
  26         mesh
  27     );
  28
  29
  30     Info<< "Reading field alpha2\n" << endl;
  31     volScalarField alpha2
  32     (
  33         IOobject
  34         (
  35             "alpha2",
  36             runTime.timeName(),
  37             mesh,
  38             IOobject::MUST_READ,
  39             IOobject::AUTO_WRITE
  40         ),
  41         mesh
  42     );
  43
  44
  45     Info<< "Reading field alpha3\n" << endl;
  46     volScalarField alpha3
  47     (
  48         IOobject
  49         (
  50             "alpha3",
  51             runTime.timeName(),
  52             mesh,
  53             IOobject::MUST_READ,
  54             IOobject::AUTO_WRITE
  55         ),
  56         mesh
  57     );
  58
  59     alpha3 == 1.0 - alpha1 - alpha2;
  60
  61
  62     Info<< "Reading field U\n" << endl;
  63     volVectorField U
  64     (
  65         IOobject
  66         (
  67             "U",
  68             runTime.timeName(),
  69             mesh,
  70             IOobject::MUST_READ,
  71             IOobject::AUTO_WRITE
  72         ),
  73         mesh
  74     );
  75
  76     #include "createPhi.H"
  77
  78     threePhaseMixture threePhaseProperties(U, phi);
  79
  80     const dimensionedScalar& rho1 = threePhaseProperties.rho1();
  81     const dimensionedScalar& rho2 = threePhaseProperties.rho2();
  82     const dimensionedScalar& rho3 = threePhaseProperties.rho3();
  83
  84     dimensionedScalar D23(threePhaseProperties.lookup("D23"));
  85
  86     // Need to store rho for ddt(rho, U)
  87     volScalarField rho
  88     (
  89         IOobject
  90         (
  91             "rho",
  92             runTime.timeName(),
  93             mesh,
  94             IOobject::READ_IF_PRESENT
  95         ),
  96         alpha1*rho1 + alpha2*rho2 + alpha3*rho3,
  97         alpha1.boundaryField().types()
  98     );
  99     rho.oldTime();
 100
 101
 102     // Mass flux
 103     // Initialisation does not matter because rhoPhi is reset after the
 104     // alpha solution before it is used in the U equation.
 105     surfaceScalarField rhoPhi
 106     (
 107         IOobject
 108         (
 109             "rho*phi",
 110             runTime.timeName(),
 111             mesh,
 112             IOobject::NO_READ,
 113             IOobject::NO_WRITE
 114         ),
 115         rho1*phi
 116     );
 117
 118
 119     // Construct interface from alpha distribution
 120     threePhaseInterfaceProperties interface(threePhaseProperties);
 121
 122
 123     // Construct incompressible turbulence model
 124     autoPtr<incompressible::turbulenceModel> turbulence
 125     (
 126         incompressible::turbulenceModel::New(U, phi, threePhaseProperties)
 127     );
 128
 129
 130     Info<< "Calculating field g.h\n" << endl;
 131     volScalarField gh("gh", g & mesh.C());
 132     surfaceScalarField ghf("ghf", g & mesh.Cf());
 133
 134     volScalarField p
 135     (
 136         IOobject
 137         (
 138             "p",
 139             runTime.timeName(),
 140             mesh,
 141             IOobject::NO_READ,
 142             IOobject::AUTO_WRITE
 143         ),
 144         p_rgh + rho*gh
 145     );
 146
 147     label pRefCell = 0;
 148     scalar pRefValue = 0.0;
 149     setRefCell
 150     (
 151         p,
 152         p_rgh,
 153         mesh.solutionDict().subDict("PIMPLE"),
 154         pRefCell,
 155         pRefValue
 156     );
 157
 158     if (p_rgh.needReference())
 159     {
 160         p += dimensionedScalar
 161         (
 162             "p",
 163             p.dimensions(),
 164             pRefValue - getRefCellValue(p, pRefCell)
 165         );
 166         p_rgh = p - rho*gh;
 167     }