diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/UEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/UEqn.H index 16be7c1ed2..e719d92433 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/UEqn.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/UEqn.H @@ -10,5 +10,16 @@ if (momentumPredictor) { - solve(UEqn() == -fvc::grad(pd) - fvc::grad(rho)*gh); + solve + ( + UEqn() + == + -fvc::reconstruct + ( + ( + fvc::snGrad(pd) + + ghf*fvc::snGrad(rho) + ) * mesh.magSf() + ) + ); } diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleMultiRegionCourantNo.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleMultiRegionCourantNo.H index 68c3621ec1..3ca2f68581 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleMultiRegionCourantNo.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleMultiRegionCourantNo.H @@ -7,8 +7,8 @@ ( fluidRegions[regionI], runTime, - rhof[regionI], - phif[regionI] + rhoFluid[regionI], + phiFluid[regionI] ), CoNum ); diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H index a1fe9e2bc1..1f6e50de0a 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H @@ -1,15 +1,16 @@ // Initialise fluid field pointer lists - PtrList thermof(fluidRegions.size()); - PtrList rhof(fluidRegions.size()); - PtrList Kf(fluidRegions.size()); - PtrList Uf(fluidRegions.size()); - PtrList phif(fluidRegions.size()); + PtrList thermoFluid(fluidRegions.size()); + PtrList rhoFluid(fluidRegions.size()); + PtrList KFluid(fluidRegions.size()); + PtrList UFluid(fluidRegions.size()); + PtrList phiFluid(fluidRegions.size()); PtrList turbulence(fluidRegions.size()); - PtrList DpDtf(fluidRegions.size()); - PtrList ghf(fluidRegions.size()); - PtrList pdf(fluidRegions.size()); + PtrList DpDtFluid(fluidRegions.size()); + PtrList ghFluid(fluidRegions.size()); + PtrList ghfFluid(fluidRegions.size()); + PtrList pdFluid(fluidRegions.size()); - List initialMassf(fluidRegions.size()); + List initialMassFluid(fluidRegions.size()); dimensionedScalar pRef ( @@ -24,8 +25,8 @@ Info<< "*** Reading fluid mesh thermophysical properties for region " << fluidRegions[i].name() << nl << endl; - Info<< " Adding to pdf\n" << endl; - pdf.set + Info<< " Adding to pdFluid\n" << endl; + pdFluid.set ( i, new volScalarField @@ -42,16 +43,15 @@ ) ); - Info<< " Adding to thermof\n" << endl; - - thermof.set + Info<< " Adding to thermoFluid\n" << endl; + thermoFluid.set ( i, basicThermo::New(fluidRegions[i]).ptr() ); - Info<< " Adding to rhof\n" << endl; - rhof.set + Info<< " Adding to rhoFluid\n" << endl; + rhoFluid.set ( i, new volScalarField @@ -64,12 +64,12 @@ IOobject::NO_READ, IOobject::AUTO_WRITE ), - thermof[i].rho() + thermoFluid[i].rho() ) ); - Info<< " Adding to Kf\n" << endl; - Kf.set + Info<< " Adding to KFluid\n" << endl; + KFluid.set ( i, new volScalarField @@ -82,12 +82,12 @@ IOobject::NO_READ, IOobject::NO_WRITE ), - thermof[i].Cp()*thermof[i].alpha() + thermoFluid[i].Cp()*thermoFluid[i].alpha() ) ); - Info<< " Adding to Uf\n" << endl; - Uf.set + Info<< " Adding to UFluid\n" << endl; + UFluid.set ( i, new volVectorField @@ -104,8 +104,8 @@ ) ); - Info<< " Adding to phif\n" << endl; - phif.set + Info<< " Adding to phiFluid\n" << endl; + phiFluid.set ( i, new surfaceScalarField @@ -118,7 +118,7 @@ IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE ), - linearInterpolate(rhof[i]*Uf[i]) + linearInterpolate(rhoFluid[i]*UFluid[i]) & fluidRegions[i].Sf() ) ); @@ -131,16 +131,16 @@ ( compressible::turbulenceModel::New ( - rhof[i], - Uf[i], - phif[i], - thermof[i] + rhoFluid[i], + UFluid[i], + phiFluid[i], + thermoFluid[i] ) ).ptr() ); - Info<< " Adding to DpDtf\n" << endl; - DpDtf.set + Info<< " Adding to DpDtFluid\n" << endl; + DpDtFluid.set ( i, new volScalarField @@ -150,9 +150,9 @@ surfaceScalarField ( "phiU", - phif[i]/fvc::interpolate(rhof[i]) + phiFluid[i]/fvc::interpolate(rhoFluid[i]) ), - thermof[i].p() + thermoFluid[i].p() ) ) ); @@ -162,8 +162,8 @@ ("environmentalProperties"); dimensionedVector g(environmentalProperties.lookup("g")); - Info<< " Adding to ghf\n" << endl; - ghf.set + Info<< " Adding to ghFluid\n" << endl; + ghFluid.set ( i, new volScalarField @@ -172,12 +172,21 @@ g & fluidRegions[i].C() ) ); + ghfFluid.set + ( + i, + new surfaceScalarField + ( + "ghf", + g & fluidRegions[i].Cf() + ) + ); Info<< " Updating p from pd\n" << endl; - thermof[i].p() == pdf[i] + rhof[i]*ghf[i] + pRef; - thermof[i].correct(); + thermoFluid[i].p() == pdFluid[i] + rhoFluid[i]*ghFluid[i] + pRef; + thermoFluid[i].correct(); - initialMassf[i] = fvc::domainIntegrate(rhof[i]).value(); + initialMassFluid[i] = fvc::domainIntegrate(rhoFluid[i]).value(); } diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H index 54507f0d19..e297989809 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H @@ -1,24 +1,31 @@ { - bool closedVolume = false; + bool closedVolume = pd.needReference(); rho = thermo.rho(); volScalarField rUA = 1.0/UEqn().A(); + surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA)); + U = rUA*UEqn().H(); - phi = + surfaceScalarField phiU + ( fvc::interpolate(rho) *( (fvc::interpolate(U) & mesh.Sf()) + fvc::ddtPhiCorr(rUA, rho, U, phi) ) - - fvc::interpolate(rho*rUA*gh)*fvc::snGrad(rho)*mesh.magSf(); + ); + + phi = phiU - ghf*fvc::snGrad(rho)*rhorUAf*mesh.magSf(); for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) { fvScalarMatrix pdEqn ( - fvc::ddt(rho) + fvm::ddt(psi, pd) + + fvc::ddt(psi)*pRef + + fvc::ddt(psi, rho)*gh + fvc::div(phi) - fvm::laplacian(rho*rUA, pd) ); @@ -38,8 +45,12 @@ } } + // Correct velocity field + U += rUA*fvc::reconstruct((phi - phiU)/rhorUAf); + U.correctBoundaryConditions(); + // Update pressure field (including bc) - thermo.p() == pd + rho*gh + pRef; + p == pd + rho*gh + pRef; DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p); // Solve continuity @@ -48,23 +59,17 @@ // Update continuity errors # include "compressibleContinuityErrors.H" - // Correct velocity field - U -= rUA*(fvc::grad(pd) + fvc::grad(rho)*gh); - U.correctBoundaryConditions(); - // For closed-volume cases adjust the pressure and density levels // to obey overall mass continuity if (closedVolume) { - thermo.p() += - ( - dimensionedScalar("massIni", dimMass, initialMassf[i]) - - fvc::domainIntegrate(thermo.psi()*thermo.p()) - )/fvc::domainIntegrate(thermo.psi()); - pd == thermo.p() - (rho*gh + pRef); + p += (massIni - fvc::domainIntegrate(psi*p))/fvc::domainIntegrate(psi); rho = thermo.rho(); } // Update thermal conductivity - K = thermo.Cp()*turb.alphaEff(); + K = thermoFluid[i].Cp()*turb.alphaEff(); + + // Update pd (including bc) + pd == p - (rho*gh + pRef); }