applications/solvers/multiphase/compressibleInterDyMFoam/createFields.H

   1     Info<< "Reading field p\n" << endl;
   2     volScalarField p
   3     (
   4         IOobject
   5         (
   6             "p",
   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     volScalarField rho1 = rho10 + psi1*p;
  87     volScalarField rho2 = rho20 + psi2*p;
  88
  89     volScalarField rho
  90     (
  91         IOobject
  92         (
  93             "rho",
  94             runTime.timeName(),
  95             mesh,
  96             IOobject::READ_IF_PRESENT,
  97             IOobject::AUTO_WRITE
  98         ),
  99         alpha1*rho1 + alpha2*rho2
 100     );
 101
 102
 103     // Mass flux
 104     // Initialisation does not matter because rhoPhi is reset after the
 105     // alpha1 solution before it is used in the U equation.
 106     surfaceScalarField rhoPhi
 107     (
 108         IOobject
 109         (
 110             "rho*phi",
 111             runTime.timeName(),
 112             mesh,
 113             IOobject::NO_READ,
 114             IOobject::NO_WRITE
 115         ),
 116         fvc::interpolate(rho)*phi
 117     );
 118
 119     volScalarField dgdt =
 120         pos(alpha2)*fvc::div(phi)/max(alpha2, scalar(0.0001));
 121
 122     // Construct interface from alpha1 distribution
 123     interfaceProperties interface(alpha1, U, twoPhaseProperties);
 124
 125     // Construct incompressible turbulence model
 126     autoPtr<incompressible::turbulenceModel> turbulence
 127     (
 128         incompressible::turbulenceModel::New(U, phi, twoPhaseProperties)
 129     );
 130
 131
 132     wordList pcorrTypes
 133     (
 134         p.boundaryField().size(),
 135         zeroGradientFvPatchScalarField::typeName
 136     );
 137
 138     for (label i=0; i<p.boundaryField().size(); i++)
 139     {
 140         if (p.boundaryField()[i].fixesValue())
 141         {
 142             pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
 143         }
 144     }