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