applications/solvers/multiphase/twoLiquidMixingFoam/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<< "Reading transportProperties\n" << endl;
  46     twoPhaseMixture twoPhaseProperties(U, phi);
  47
  48     const dimensionedScalar& rho1 = twoPhaseProperties.rho1();
  49     const dimensionedScalar& rho2 = twoPhaseProperties.rho2();
  50
  51     dimensionedScalar Dab(twoPhaseProperties.lookup("Dab"));
  52
  53     // Read the reciprocal of the turbulent Schmidt number
  54     dimensionedScalar alphatab(twoPhaseProperties.lookup("alphatab"));
  55
  56     // Need to store rho for ddt(rho, U)
  57     volScalarField rho("rho", alpha1*rho1 + (scalar(1) - alpha1)*rho2);
  58     rho.oldTime();
  59
  60
  61     // Mass flux
  62     // Initialisation does not matter because rhoPhi is reset after the
  63     // alpha1 solution before it is used in the U equation.
  64     surfaceScalarField rhoPhi
  65     (
  66         IOobject
  67         (
  68             "rho*phi",
  69             runTime.timeName(),
  70             mesh,
  71             IOobject::NO_READ,
  72             IOobject::NO_WRITE
  73         ),
  74         rho1*phi
  75     );
  76
  77     Info<< "Calculating field (g.h)f\n" << endl;
  78     volScalarField gh("gh", g & mesh.C());
  79     surfaceScalarField ghf = surfaceScalarField("gh", g & mesh.Cf());
  80
  81     volScalarField p
  82     (
  83         IOobject
  84         (
  85             "p",
  86             runTime.timeName(),
  87             mesh,
  88             IOobject::NO_READ,
  89             IOobject::AUTO_WRITE
  90         ),
  91         p_rgh + rho*gh
  92     );
  93
  94     label p_rghRefCell = 0;
  95     scalar p_rghRefValue = 0.0;
  96     setRefCell
  97     (
  98         p_rgh,
  99         mesh.solutionDict().subDict("PIMPLE"),
 100         p_rghRefCell,
 101         p_rghRefValue
 102     );
 103
 104     scalar pRefValue = 0.0;
 105
 106     if (p_rgh.needReference())
 107     {
 108         pRefValue = readScalar
 109         (
 110             mesh.solutionDict().subDict("PIMPLE").lookup("pRefValue")
 111         );
 112
 113         p += dimensionedScalar
 114         (
 115             "p",
 116             p.dimensions(),
 117             pRefValue - getRefCellValue(p, p_rghRefCell)
 118         );
 119     }
 120
 121
 122     // Construct incompressible turbulence model
 123     autoPtr<incompressible::turbulenceModel> turbulence
 124     (
 125         incompressible::turbulenceModel::New(U, phi, twoPhaseProperties)
 126     );