applications/solvers/multiphase/interPhaseChangeFoam/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     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     #include "createPhi.H"
  44
  45     Info<< "Creating phaseChangeTwoPhaseMixture\n" << endl;
  46     autoPtr<phaseChangeTwoPhaseMixture> twoPhaseProperties =
  47         phaseChangeTwoPhaseMixture::New(U, phi);
  48
  49     const dimensionedScalar& rho1 = twoPhaseProperties->rho1();
  50     const dimensionedScalar& rho2 = twoPhaseProperties->rho2();
  51     const dimensionedScalar& pSat = twoPhaseProperties->pSat();
  52
  53     // Need to store rho for ddt(rho, U)
  54     volScalarField rho
  55     (
  56         IOobject
  57         (
  58             "rho",
  59             runTime.timeName(),
  60             mesh,
  61             IOobject::READ_IF_PRESENT
  62         ),
  63         alpha1*rho1 + (scalar(1) - alpha1)*rho2,
  64         alpha1.boundaryField().types()
  65     );
  66     rho.oldTime();
  67
  68     // Construct interface from alpha1 distribution
  69     interfaceProperties interface(alpha1, U, twoPhaseProperties());
  70
  71     // Construct incompressible turbulence model
  72     autoPtr<incompressible::turbulenceModel> turbulence
  73     (
  74         incompressible::turbulenceModel::New(U, phi, twoPhaseProperties())
  75     );
  76
  77
  78     Info<< "Calculating field g.h\n" << endl;
  79     volScalarField gh("gh", g & mesh.C());
  80     surfaceScalarField ghf("ghf", g & mesh.Cf());
  81
  82     volScalarField p
  83     (
  84         IOobject
  85         (
  86             "p",
  87             runTime.timeName(),
  88             mesh,
  89             IOobject::NO_READ,
  90             IOobject::AUTO_WRITE
  91         ),
  92         p_rgh + rho*gh
  93     );
  94
  95     label pRefCell = 0;
  96     scalar pRefValue = 0.0;
  97     setRefCell
  98     (
  99         p,
 100         p_rgh,
 101         mesh.solutionDict().subDict("PIMPLE"),
 102         pRefCell,
 103         pRefValue
 104     );
 105
 106     if (p_rgh.needReference())
 107     {
 108         p += dimensionedScalar
 109         (
 110             "p",
 111             p.dimensions(),
 112             pRefValue - getRefCellValue(p, pRefCell)
 113         );
 114         p_rgh = p - rho*gh;
 115     }