Commit 8ec8f595 authored by sergio's avatar sergio
Browse files

ENH: cleaning chtMultiRegionFoam and chtMultiRegionSimpleFoam from

porous fluid and solid now incorporated into fluid and solid regions
parent ee01df6f
......@@ -60,14 +60,9 @@ int main(int argc, char *argv[])
#include "createFluidMeshes.H"
#include "createSolidMeshes.H"
#include "createPorousFluidRegions.H"
#include "createPorousSolidMeshes.H"
#include "createFluidFields.H"
#include "createSolidFields.H"
#include "createPorousFluidFields.H"
#include "createPorousSolidFields.H"
#include "initContinuityErrs.H"
#include "readTimeControls.H"
......@@ -116,24 +111,6 @@ int main(int argc, char *argv[])
#include "solveFluid.H"
}
forAll(porousFluidRegions, i)
{
Info<< "\nSolving for fluid porous region "
<< porousFluidRegions[i].name() << endl;
#include "setPorousFluidFields.H"
#include "readPorousFluidRegionPIMPLEControls.H"
#include "solvePorousFluid.H"
}
forAll(porousSolidRegions, i)
{
Info<< "\nSolving for porous solid region "
<< porousSolidRegions[i].name() << endl;
#include "setPorousRegionSolidFields.H"
#include "readPorousSolidMultiRegionPIMPLEControls.H"
#include "solvePorousSolid.H"
}
forAll(solidRegions, i)
{
Info<< "\nSolving for solid region "
......
EXE_INC = \
-Ifluid \
-Isolid \
-I./porousFluid \
-I./porousSolid \
-I../solid \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/finiteVolume/cfdTools \
......
......@@ -50,13 +50,9 @@ int main(int argc, char *argv[])
#include "createFluidMeshes.H"
#include "createSolidMeshes.H"
#include "createPorousFluidRegions.H"
#include "createPorousSolidMeshes.H"
#include "createFluidFields.H"
#include "createSolidFields.H"
#include "createPorousFluidFields.H"
#include "createPorousSolidFields.H"
#include "initContinuityErrs.H"
......@@ -74,24 +70,6 @@ int main(int argc, char *argv[])
#include "solveFluid.H"
}
forAll(porousFluidRegions, i)
{
Info<< "\nSolving for fluid porous region "
<< porousFluidRegions[i].name() << endl;
#include "setPorousFluidFields.H"
#include "readPorousFluidRegionSIMPLEControls.H"
#include "solvePorousFluid.H"
}
forAll(porousSolidRegions, i)
{
Info<< "\nSolving for porous solid region "
<< porousSolidRegions[i].name() << endl;
#include "setPorousRegionSolidFields.H"
#include "readPorousSolidMultiRegionSIMPLEControls.H"
#include "solvePorousSolid.H"
}
forAll(solidRegions, i)
{
Info<< "\nSolving for solid region "
......
// Solve the Momentum equation
tmp<fvVectorMatrix> porousUEqn
(
fvm::div(porousPhi, porousU)
+ turbPorous.divDevRhoReff(porousU)
+ porousSources(porousRho, porousU)
);
porousUEqn().relax();
solve(porousUEqn() == -fvc::grad(porousP));
// Initialise porous field pointer lists
PtrList<rhoThermo> thermoPorous(porousFluidRegions.size());
PtrList<volScalarField> rhoPorous(porousFluidRegions.size());
PtrList<volScalarField> kappaPorous(porousFluidRegions.size());
PtrList<volVectorField> UPorous(porousFluidRegions.size());
PtrList<surfaceScalarField> phiPorous(porousFluidRegions.size());
PtrList<compressible::turbulenceModel> turbulencePorous
(
porousFluidRegions.size()
);
PtrList<volScalarField> pPorous(porousFluidRegions.size());
List<scalar> initialMassFluidPorous(porousFluidRegions.size());
List<label> pRefCellFluidPorous(porousFluidRegions.size(),0);
List<scalar> pRefValueFluidPorous(porousFluidRegions.size(),0.0);
PtrList<dimensionedScalar> rhoMaxPorous(fluidRegions.size());
PtrList<dimensionedScalar> rhoMinPorous(fluidRegions.size());
PtrList<IObasicSourceList> heatPorousSources
(
porousFluidRegions.size()
);
forAll(porousFluidRegions, i)
{
Info<< "Reading fluid mesh thermophysical properties for porous "
<< porousFluidRegions[i].name() << nl << endl;
Info<< " Adding to thermoFluid porous\n" << endl;
thermoPorous.set
(
i,
rhoThermo::New(porousFluidRegions[i]).ptr()
);
Info<< " Adding to rhoPorous\n" << endl;
rhoPorous.set
(
i,
new volScalarField
(
IOobject
(
"rho",
runTime.timeName(),
porousFluidRegions[i],
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
thermoPorous[i].rho()
)
);
Info<< " Adding to UPorous\n" << endl;
UPorous.set
(
i,
new volVectorField
(
IOobject
(
"U",
runTime.timeName(),
porousFluidRegions[i],
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
porousFluidRegions[i]
)
);
Info<< " Adding to phiPorous\n" << endl;
phiPorous.set
(
i,
new surfaceScalarField
(
IOobject
(
"phi",
runTime.timeName(),
porousFluidRegions[i],
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
linearInterpolate(rhoPorous[i]*UPorous[i])
& porousFluidRegions[i].Sf()
)
);
Info<< " Adding turbulence to porous\n" << endl;
turbulencePorous.set
(
i,
compressible::turbulenceModel::New
(
rhoPorous[i],
UPorous[i],
phiPorous[i],
thermoPorous[i]
).ptr()
);
Info<< " Adding to kappaFluid\n" << endl;
kappaPorous.set
(
i,
new volScalarField
(
IOobject
(
"kappaPorous",
runTime.timeName(),
porousFluidRegions[i],
IOobject::NO_READ,
IOobject::NO_WRITE
),
thermoPorous[i].Cp()*thermoPorous[i].alpha()
)
);
pPorous.set
(
i,
new volScalarField
(
IOobject
(
"p",
runTime.timeName(),
porousFluidRegions[i],
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
porousFluidRegions[i]
)
);
setRefCell
(
thermoPorous[i].p(),
pPorous[i],
porousFluidRegions[i].solutionDict().subDict("SIMPLE"),
pRefCellFluidPorous[i],
pRefValueFluidPorous[i]
);
rhoMaxPorous.set
(
i,
new dimensionedScalar
(
porousFluidRegions[i].solutionDict().subDict("SIMPLE").lookup
(
"rhoMax"
)
)
);
rhoMinPorous.set
(
i,
new dimensionedScalar
(
porousFluidRegions[i].solutionDict().subDict("SIMPLE").lookup
(
"rhoMin"
)
)
);
heatPorousSources.set
(
i,
new IObasicSourceList(porousFluidRegions[i])
);
}
const wordList porousFluidNames(rp["porousFluid"]);
PtrList<fvMesh> porousFluidRegions(porousFluidNames.size());
forAll (porousFluidNames, iPorous)
{
const word porousFluidName = porousFluidNames[iPorous];
Info<< "Create porous fluid region " << porousFluidName
<< nl << endl;
porousFluidRegions.set
(
iPorous,
new fvMesh
(
IOobject
(
porousFluidName,
runTime.timeName(),
runTime,
IOobject::MUST_READ
)
)
);
}
{
fvScalarMatrix hPorousEqn
(
fvm::div(porousPhi, porousH)
- fvm::laplacian(turbPorous.alphaEff(), porousH)
==
- fvc::div(porousPhi, 0.5*magSqr(porousU), "div(phi,K)")
+ porousSources(porousRho, porousH)
);
hPorousEqn.relax();
hPorousEqn.solve();
porousThermo.correct();
Info<< "Min/max in the porous T:"
<< min(porousThermo.T()).value() << ' '
<< max(porousThermo.T()).value() << endl;
}
porousRho = porousThermo.rho();
porousRho = max(porousRho, rhoMin);
porousRho = min(porousRho, rhoMax);
porousRho.relax();
volScalarField rAUPorous(1.0/porousUEqn().A());
porousU = rAUPorous*porousUEqn().H();
porousUEqn.clear();
bool closedVolume = false;
porousPhi =
fvc::interpolate(porousRho)
*(fvc::interpolate(porousU) & porousMesh.Sf());
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqn
(
fvm::laplacian(porousRho*rAUPorous, porousP) == fvc::div(porousPhi)
);
pEqn.setReference(pRefCell, pRefValue);
pEqn.solve();
if (nonOrth == nNonOrthCorr)
{
porousPhi -= pEqn.flux();
}
}
porousP.relax();
porousU -= rAUPorous*fvc::grad(porousP);
porousU.correctBoundaryConditions();
if (closedVolume)
{
porousP += (initialMass - fvc::domainIntegrate(porousPsi*porousP))
/fvc::domainIntegrate(porousPsi);
}
porousRho = porousThermo.rho();
porousRho = max(porousRho, rhoMin);
porousRho = min(porousRho, rhoMax);
porousRho.relax();
Info<< "rho max/min : "
<< max(porousRho).value() << " "
<< min(porousRho).value() << endl;
const dictionary& simple = porousMesh.solutionDict().subDict("SIMPLE");
const int nNonOrthCorr =
simple.lookupOrDefault<int>("nNonOrthogonalCorrectors", 0);
const fvMesh& porousMesh = porousFluidRegions[i];
rhoThermo& porousThermo = thermoPorous[i];
volScalarField& porousRho = rhoPorous[i];
volVectorField& porousU = UPorous[i];
surfaceScalarField& porousPhi = phiPorous[i];
compressible::turbulenceModel& turbPorous = turbulencePorous[i];
volScalarField& porousP = porousThermo.p();
const volScalarField& porousPsi = porousThermo.psi();
volScalarField& porousH = porousThermo.he();
const dimensionedScalar initialMass
(
"initialMass",
dimMass,
initialMassFluidPorous[i]
);
IObasicSourceList& porousSources = heatPorousSources[i];
const label pRefCell = pRefCellFluidPorous[i];
const scalar pRefValue = pRefValueFluidPorous[i];
const scalar rhoMax = rhoMaxPorous[i].value();
const scalar rhoMin = rhoMinPorous[i].value();
// Pressure-velocity SIMPLE corrector
porousP.storePrevIter();
porousRho.storePrevIter();
{
#include "UPorousFluidEqn.H"
#include "hPorousFluidEqn.H"
#include "pPorousFluidEqn.H"
}
turbPorous.correct();
// Initialise solid field pointer lists
PtrList<solidThermo> porousSolidThermos(porousSolidRegions.size());
PtrList<IObasicSourceList> solidHeatSources(porousSolidRegions.size());
PtrList<volScalarField> betavSolid(porousSolidRegions.size());
// Populate solid field pointer lists
forAll(porousSolidRegions, i)
{
Info<< "*** Reading porous solid mesh thermophysical "
<< "properties for region "
<< porousSolidRegions[i].name() << nl << endl;
Info<< " Adding to thermos\n" << endl;
porousSolidThermos.set
(
i,
solidThermo::New(porousSolidRegions[i])
);
Info<< " Adding sources\n" << endl;
solidHeatSources.set
(
i,
new IObasicSourceList(porousSolidRegions[i])
);
Info<< " Adding to betavSolid\n" << endl;
betavSolid.set
(
i,
new volScalarField
(
IOobject
(
"betavSolid",
runTime.timeName(),
porousSolidRegions[i],
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
porousSolidRegions[i]
)
);
}
const wordList porousSolidNames(rp["porousSolid"]);
PtrList<fvMesh> porousSolidRegions(porousSolidNames.size());
forAll(porousSolidNames, i)
{
Info<< "Create solid mesh for region " << porousSolidNames[i]
<< " for time = " << runTime.timeName() << nl << endl;
porousSolidRegions.set
(
i,
new fvMesh
(
IOobject
(
porousSolidNames[i],
runTime.timeName(),
runTime,
IOobject::MUST_READ
)
)
);
}
const dictionary& simple = mesh.solutionDict().subDict("SIMPLE");
const int nNonOrthCorr =
simple.lookupOrDefault<int>("nNonOrthogonalCorrectors", 0);
const fvMesh& mesh = porousSolidRegions[i];
solidThermo& thermo = porousSolidThermos[i];
const volScalarField& betav = betavSolid[i];
tmp<volScalarField> trho = thermo.rho();
const volScalarField& rho = trho();
tmp<volScalarField> tcp = thermo.Cp();
const volScalarField& cp = tcp();
tmp<volScalarField> tkappa = thermo.kappa();
//tmp<volSymmTensorField> tkappa = thermo.directionalKappa()*betav;
const volScalarField& kappa = tkappa();
//const volSymmTensorField& K = tK();
tmp<volScalarField> talpha = thermo.alpha();
const volScalarField& alpha = talpha();
volScalarField& h = thermo.he();
IObasicSourceList& sources = solidHeatSources[i];
scalar DiNum = -GREAT;
forAll(solidRegions, i)
{
# include "setRegionSolidFields.H"
DiNum = max
(
solidRegionDiffNo
(
solidRegions[i],
runTime,
rho*cp,
K
),
DiNum
);
}
{
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
tmp<fvScalarMatrix> hEqn