openfoam/applications/solvers/multiphase/interDyMFoam/createFields.H

150 lines
3.2 KiB
C

Info<< "Reading field pd\n" << endl;
volScalarField pd
(
IOobject
(
"pd",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading field alpha1\n" << endl;
volScalarField alpha1
(
IOobject
(
"alpha1",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
# include "createPhi.H"
Info<< "Reading transportProperties\n" << endl;
twoPhaseMixture twoPhaseProperties(U, phi);
const dimensionedScalar& rho1 = twoPhaseProperties.rho1();
const dimensionedScalar& rho2 = twoPhaseProperties.rho2();
// Need to store rho for ddt(rho, U)
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT
),
alpha1*rho1 + (scalar(1) - alpha1)*rho2,
alpha1.boundaryField().types()
);
rho.oldTime();
// Mass flux
// Initialisation does not matter because rhoPhi is reset after the
// alpha1 solution before it is used in the U equation.
surfaceScalarField rhoPhi
(
IOobject
(
"rho*phi",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
rho1*phi
);
// Construct interface from alpha1 distribution
interfaceProperties interface(alpha1, U, twoPhaseProperties);
// Construct incompressible turbulence model
autoPtr<incompressible::turbulenceModel> turbulence
(
incompressible::turbulenceModel::New(U, phi, twoPhaseProperties)
);
wordList pcorrTypes(pd.boundaryField().types());
for (label i=0; i<pd.boundaryField().size(); i++)
{
if (pd.boundaryField()[i].fixesValue())
{
pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
}
}
volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
pd + rho*(g & mesh.C())
);
autoPtr<probes> pRefProbe;
label pdRefCell = 0;
scalar pdRefValue = 0.0;
scalar pRefValue = 0.0;
if (pd.needReference())
{
pRefProbe.set
(
new probes
(
"pRefProbe",
mesh,
mesh.solutionDict().subDict("PISO").subDict("pRefProbe")
)
);
pRefValue = readScalar
(
mesh.solutionDict().subDict("PISO").lookup("pRefValue")
);
pdRefCell = pRefProbe->cells()[0];
p += dimensionedScalar
(
"p",
p.dimensions(),
pRefValue - pRefProbe->sample<scalar>("p")()[0]
);
}