applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H

   1     Info<< "Reading transportProperties\n" << endl;
   2
   3     IOdictionary transportProperties
   4     (
   5         IOobject
   6         (
   7             "transportProperties",
   8             runTime.constant(),
   9             mesh,
  10             IOobject::MUST_READ,
  11             IOobject::NO_WRITE
  12         )
  13     );
  14
  15     autoPtr<phaseModel> phasea = phaseModel::New
  16     (
  17         mesh,
  18         transportProperties,
  19         "a"
  20     );
  21
  22     autoPtr<phaseModel> phaseb = phaseModel::New
  23     (
  24         mesh,
  25         transportProperties,
  26         "b"
  27     );
  28
  29     volVectorField& Ua = phasea->U();
  30     surfaceScalarField& phia = phasea->phi();
  31     const dimensionedScalar& rhoa = phasea->rho();
  32     const dimensionedScalar& nua = phasea->nu();
  33
  34     volVectorField& Ub = phaseb->U();
  35     surfaceScalarField& phib = phaseb->phi();
  36     const dimensionedScalar& rhob = phaseb->rho();
  37     const dimensionedScalar& nub = phaseb->nu();
  38
  39     Info<< "Reading field alpha\n" << endl;
  40     volScalarField alpha
  41     (
  42         IOobject
  43         (
  44             "alpha",
  45             runTime.timeName(),
  46             mesh,
  47             IOobject::MUST_READ,
  48             IOobject::AUTO_WRITE
  49         ),
  50         mesh
  51     );
  52
  53     volScalarField beta
  54     (
  55         IOobject
  56         (
  57             "beta",
  58             runTime.timeName(),
  59             mesh,
  60             IOobject::NO_READ,
  61             IOobject::NO_WRITE
  62         ),
  63         scalar(1) - alpha
  64         //,alpha.boundaryField().types()
  65     );
  66
  67     Info<< "Reading field p\n" << endl;
  68     volScalarField p
  69     (
  70         IOobject
  71         (
  72             "p",
  73             runTime.timeName(),
  74             mesh,
  75             IOobject::MUST_READ,
  76             IOobject::AUTO_WRITE
  77         ),
  78         mesh
  79     );
  80
  81     volVectorField U
  82     (
  83         IOobject
  84         (
  85             "U",
  86             runTime.timeName(),
  87             mesh,
  88             IOobject::NO_READ,
  89             IOobject::AUTO_WRITE
  90         ),
  91         alpha*Ua + beta*Ub
  92     );
  93
  94     dimensionedScalar Cvm
  95     (
  96         transportProperties.lookup("Cvm")
  97     );
  98
  99     dimensionedScalar Cl
 100     (
 101         transportProperties.lookup("Cl")
 102     );
 103
 104     dimensionedScalar Ct
 105     (
 106         transportProperties.lookup("Ct")
 107     );
 108
 109     surfaceScalarField phi
 110     (
 111         IOobject
 112         (
 113             "phi",
 114             runTime.timeName(),
 115             mesh
 116         ),
 117         fvc::interpolate(alpha)*phia + fvc::interpolate(beta)*phib
 118     );
 119
 120     volScalarField rho
 121     (
 122         IOobject
 123         (
 124             "rho",
 125             runTime.timeName(),
 126             mesh
 127         ),
 128         alpha*rhoa + beta*rhob
 129     );
 130
 131     #include "createRASTurbulence.H"
 132
 133     Info<< "Calculating field DDtUa and DDtUb\n" << endl;
 134
 135     volVectorField DDtUa
 136     (
 137         fvc::ddt(Ua)
 138       + fvc::div(phia, Ua)
 139       - fvc::div(phia)*Ua
 140     );
 141
 142     volVectorField DDtUb
 143     (
 144         fvc::ddt(Ub)
 145       + fvc::div(phib, Ub)
 146       - fvc::div(phib)*Ub
 147     );
 148
 149
 150     Info<< "Calculating field g.h\n" << endl;
 151     volScalarField gh("gh", g & mesh.C());
 152
 153     IOdictionary interfacialProperties
 154     (
 155         IOobject
 156         (
 157             "interfacialProperties",
 158             runTime.constant(),
 159             mesh,
 160             IOobject::MUST_READ,
 161             IOobject::NO_WRITE
 162         )
 163     );
 164
 165     autoPtr<dragModel> draga = dragModel::New
 166     (
 167         interfacialProperties,
 168         alpha,
 169         phasea,
 170         phaseb
 171     );
 172
 173     autoPtr<dragModel> dragb = dragModel::New
 174     (
 175         interfacialProperties,
 176         beta,
 177         phaseb,
 178         phasea
 179     );
 180
 181     word dragPhase("blended");
 182     if (interfacialProperties.found("dragPhase"))
 183     {
 184         dragPhase = word(interfacialProperties.lookup("dragPhase"));
 185
 186         bool validDrag =
 187             dragPhase == "a" || dragPhase == "b" || dragPhase == "blended";
 188
 189         if (!validDrag)
 190         {
 191             FatalErrorIn(args.executable())
 192                 << "invalid dragPhase " << dragPhase
 193                 << exit(FatalError);
 194         }
 195     }
 196
 197     Info << "dragPhase is " << dragPhase << endl;
 198     kineticTheoryModel kineticTheory
 199     (
 200         phasea,
 201         Ub,
 202         alpha,
 203         draga
 204     );
 205
 206     surfaceScalarField rUaAf
 207     (
 208         IOobject
 209         (
 210             "rUaAf",
 211             runTime.timeName(),
 212             mesh,
 213             IOobject::NO_READ,
 214             IOobject::NO_WRITE
 215         ),
 216         mesh,
 217         dimensionedScalar("zero", dimensionSet(0, 0, 1, 0, 0), 0.0)
 218     );
 219
 220     surfaceScalarField ppMagf
 221     (
 222         IOobject
 223         (
 224             "ppMagf",
 225             runTime.timeName(),
 226             mesh,
 227             IOobject::NO_READ,
 228             IOobject::NO_WRITE
 229         ),
 230         mesh,
 231         dimensionedScalar("zero", dimensionSet(0, 2, -1, 0, 0), 0.0)
 232     );
 233
 234
 235     label pRefCell = 0;
 236     scalar pRefValue = 0.0;
 237     setRefCell(p, mesh.solutionDict().subDict("PIMPLE"), pRefCell, pRefValue);