applications/solvers/multiphase/compressibleInterFoam/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<< "Calculating field alpha1\n" << endl;
  30     volScalarField alpha2("alpha2", scalar(1) - alpha1);
  31
  32     Info<< "Reading field U\n" << endl;
  33     volVectorField U
  34     (
  35         IOobject
  36         (
  37             "U",
  38             runTime.timeName(),
  39             mesh,
  40             IOobject::MUST_READ,
  41             IOobject::AUTO_WRITE
  42         ),
  43         mesh
  44     );
  45
  46     #include "createPhi.H"
  47
  48
  49     Info<< "Reading transportProperties\n" << endl;
  50     twoPhaseMixture twoPhaseProperties(U, phi);
  51
  52     dimensionedScalar rho10
  53     (
  54         twoPhaseProperties.subDict
  55         (
  56             twoPhaseProperties.phase1Name()
  57         ).lookup("rho0")
  58     );
  59
  60     dimensionedScalar rho20
  61     (
  62         twoPhaseProperties.subDict
  63         (
  64             twoPhaseProperties.phase2Name()
  65         ).lookup("rho0")
  66     );
  67
  68     dimensionedScalar psi1
  69     (
  70         twoPhaseProperties.subDict
  71         (
  72             twoPhaseProperties.phase1Name()
  73         ).lookup("psi")
  74     );
  75
  76     dimensionedScalar psi2
  77     (
  78         twoPhaseProperties.subDict
  79         (
  80             twoPhaseProperties.phase2Name()
  81         ).lookup("psi")
  82     );
  83
  84     dimensionedScalar pMin(twoPhaseProperties.lookup("pMin"));
  85
  86     Info<< "Calculating field g.h\n" << endl;
  87     volScalarField gh("gh", g & mesh.C());
  88     surfaceScalarField ghf("ghf", g & mesh.Cf());
  89
  90     volScalarField p
  91     (
  92         IOobject
  93         (
  94             "p",
  95             runTime.timeName(),
  96             mesh,
  97             IOobject::NO_READ,
  98             IOobject::AUTO_WRITE
  99         ),
 100         max
 101         (
 102             (p_rgh + gh*(alpha1*rho10 + alpha2*rho20))
 103            /(1.0 - gh*(alpha1*psi1 + alpha2*psi2)),
 104             pMin
 105         )
 106     );
 107
 108     volScalarField rho1(rho10 + psi1*p);
 109     volScalarField rho2(rho20 + psi2*p);
 110
 111     volScalarField rho
 112     (
 113         IOobject
 114         (
 115             "rho",
 116             runTime.timeName(),
 117             mesh,
 118             IOobject::READ_IF_PRESENT,
 119             IOobject::AUTO_WRITE
 120         ),
 121         alpha1*rho1 + alpha2*rho2
 122     );
 123
 124
 125     // Mass flux
 126     // Initialisation does not matter because rhoPhi is reset after the
 127     // alpha1 solution before it is used in the U equation.
 128     surfaceScalarField rhoPhi
 129     (
 130         IOobject
 131         (
 132             "rho*phi",
 133             runTime.timeName(),
 134             mesh,
 135             IOobject::NO_READ,
 136             IOobject::NO_WRITE
 137         ),
 138         fvc::interpolate(rho)*phi
 139     );
 140
 141     volScalarField dgdt
 142     (
 143         pos(alpha2)*fvc::div(phi)/max(alpha2, scalar(0.0001))
 144     );
 145
 146     // Construct interface from alpha1 distribution
 147     interfaceProperties interface(alpha1, U, twoPhaseProperties);
 148
 149     // Construct incompressible turbulence model
 150     autoPtr<incompressible::turbulenceModel> turbulence
 151     (
 152         incompressible::turbulenceModel::New(U, phi, twoPhaseProperties)
 153     );