applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H

   1     Info<< "Reading field pd\n" << endl;
   2     volScalarField pd
   3     (
   4         IOobject
   5         (
   6             "pd",
   7             runTime.timeName(),
   8             mesh,
   9             IOobject::MUST_READ,
  10             IOobject::AUTO_WRITE
  11         ),
  12         mesh
  13     );
  14
  15     Info<< "Reading field gamma\n" << endl;
  16     volScalarField gamma
  17     (
  18         IOobject
  19         (
  20             "gamma",
  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<< "Reading transportProperties\n" << endl;
  46     twoPhaseMixture twoPhaseProperties(U, phi, "gamma");
  47
  48     const dimensionedScalar& rho1 = twoPhaseProperties.rho1();
  49     const dimensionedScalar& rho2 = twoPhaseProperties.rho2();
  50
  51     dimensionedScalar Dab(twoPhaseProperties.lookup("Dab"));
  52
  53     // Need to store rho for ddt(rho, U)
  54     volScalarField rho("rho", gamma*rho1 + (scalar(1) - gamma)*rho2);
  55     rho.oldTime();
  56
  57
  58     // Mass flux
  59     // Initialisation does not matter because rhoPhi is reset after the
  60     // gamma solution before it is used in the U equation.
  61     surfaceScalarField rhoPhi
  62     (
  63         IOobject
  64         (
  65             "rho*phi",
  66             runTime.timeName(),
  67             mesh,
  68             IOobject::NO_READ,
  69             IOobject::NO_WRITE
  70         ),
  71         rho1*phi
  72     );
  73
  74
  75     Info<< "Calculating field g.h\n" << endl;
  76     surfaceScalarField ghf("gh", g & mesh.Cf());
  77
  78
  79     label pdRefCell = 0;
  80     scalar pdRefValue = 0.0;
  81     setRefCell(pd, mesh.solutionDict().subDict("PISO"), pdRefCell, pdRefValue);