From a434d4f9dec5e16629ee0046e32cc590bb43f494 Mon Sep 17 00:00:00 2001 From: andy Date: Tue, 3 Feb 2009 16:29:39 +0000 Subject: [PATCH 1/8] Prantdl -> turbulent Prandtl number --- .../solvers/heatTransfer/buoyantBoussinesqSimpleFoam/TEqn.H | 2 +- .../buoyantBoussinesqSimpleFoam/readTransportProperties.H | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/TEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/TEqn.H index faddea825b..d00f0877ea 100644 --- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/TEqn.H +++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/TEqn.H @@ -2,7 +2,7 @@ volScalarField kappaEff ( "kappaEff", - turbulence->nu() + turbulence->nut()/Pr + turbulence->nu() + turbulence->nut()/Prt ); fvScalarMatrix TEqn diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/readTransportProperties.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/readTransportProperties.H index eeaa44e155..585128dfde 100644 --- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/readTransportProperties.H +++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/readTransportProperties.H @@ -9,5 +9,5 @@ // reference kinematic pressure [m2/s2] dimensionedScalar pRef(laminarTransport.lookup("pRef")); - // Prandtl number - dimensionedScalar Pr(laminarTransport.lookup("Pr")); + // turbulent Prandtl number + dimensionedScalar Prt(laminarTransport.lookup("Prt")); From bb8bcdab0323f6f3ed02e38b633d882b581323bf Mon Sep 17 00:00:00 2001 From: andy Date: Tue, 3 Feb 2009 16:31:17 +0000 Subject: [PATCH 2/8] using wallFunctionDict instead of using subDict lookup --- src/turbulenceModels/compressible/RAS/RASModel/RASModel.C | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/turbulenceModels/compressible/RAS/RASModel/RASModel.C b/src/turbulenceModels/compressible/RAS/RASModel/RASModel.C index 60d3dc3e84..b349d2fbc8 100644 --- a/src/turbulenceModels/compressible/RAS/RASModel/RASModel.C +++ b/src/turbulenceModels/compressible/RAS/RASModel/RASModel.C @@ -88,7 +88,7 @@ RASModel::RASModel dimensioned::lookupOrAddToDict ( "kappa", - subDict("wallFunctionCoeffs"), + wallFunctionDict_, 0.4187 ) ), @@ -97,7 +97,7 @@ RASModel::RASModel dimensioned::lookupOrAddToDict ( "E", - subDict("wallFunctionCoeffs"), + wallFunctionDict_ 9.0 ) ), From 338b72b1ebbe9b83f81c3d021f9d71c0205c5600 Mon Sep 17 00:00:00 2001 From: henry Date: Tue, 3 Feb 2009 16:51:07 +0000 Subject: [PATCH 3/8] Improved the flux-velocity correspondence for cases where hydrostatic balance is important e.g. in atmospheric flows. --- .../solvers/heatTransfer/buoyantFoam/UEqn.H | 16 ++- .../heatTransfer/buoyantFoam/buoyantFoam.C | 28 ++--- .../heatTransfer/buoyantFoam/createFields.H | 1 + .../solvers/heatTransfer/buoyantFoam/pEqn.H | 103 ++++++++++-------- .../heatTransfer/buoyantSimpleFoam/UEqn.H | 10 +- .../buoyantSimpleFoam/createFields.H | 1 + .../heatTransfer/buoyantSimpleFoam/pEqn.H | 98 +++++++++-------- .../buoyantSimpleRadiationFoam/Make/options | 1 + .../buoyantSimpleRadiationFoam/UEqn.H | 18 --- .../buoyantSimpleRadiationFoam/createFields.H | 1 + .../buoyantSimpleRadiationFoam/pEqn.H | 54 --------- 11 files changed, 151 insertions(+), 180 deletions(-) delete mode 100644 applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/UEqn.H delete mode 100644 applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/pEqn.H diff --git a/applications/solvers/heatTransfer/buoyantFoam/UEqn.H b/applications/solvers/heatTransfer/buoyantFoam/UEqn.H index 988268c4ef..24cfbf4f68 100644 --- a/applications/solvers/heatTransfer/buoyantFoam/UEqn.H +++ b/applications/solvers/heatTransfer/buoyantFoam/UEqn.H @@ -9,4 +9,18 @@ UEqn.relax(); - solve(UEqn == -fvc::grad(pd) - fvc::grad(rho)*gh); + if (momentumPredictor) + { + solve + ( + UEqn + == + -fvc::reconstruct + ( + ( + fvc::snGrad(pd) + + ghf*fvc::snGrad(rho) + ) * mesh.magSf() + ) + ); + } diff --git a/applications/solvers/heatTransfer/buoyantFoam/buoyantFoam.C b/applications/solvers/heatTransfer/buoyantFoam/buoyantFoam.C index cafb677900..8b3cca8d21 100644 --- a/applications/solvers/heatTransfer/buoyantFoam/buoyantFoam.C +++ b/applications/solvers/heatTransfer/buoyantFoam/buoyantFoam.C @@ -41,13 +41,12 @@ Description int main(int argc, char *argv[]) { - -# include "setRootCase.H" -# include "createTime.H" -# include "createMesh.H" -# include "readEnvironmentalProperties.H" -# include "createFields.H" -# include "initContinuityErrs.H" + #include "setRootCase.H" + #include "createTime.H" + #include "createMesh.H" + #include "readEnvironmentalProperties.H" + #include "createFields.H" + #include "initContinuityErrs.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -55,23 +54,24 @@ int main(int argc, char *argv[]) while (runTime.run()) { -# include "readPISOControls.H" -# include "compressibleCourantNo.H" -//# include "setDeltaT.H" + #include "readTimeControls.H" + #include "readPISOControls.H" + #include "compressibleCourantNo.H" + #include "setDeltaT.H" runTime++; Info<< "Time = " << runTime.timeName() << nl << endl; -# include "rhoEqn.H" + #include "rhoEqn.H" -# include "UEqn.H" + #include "UEqn.H" // --- PISO loop for (int corr=0; corrcorrect(); diff --git a/applications/solvers/heatTransfer/buoyantFoam/createFields.H b/applications/solvers/heatTransfer/buoyantFoam/createFields.H index 3d604d4daf..aa32325c34 100644 --- a/applications/solvers/heatTransfer/buoyantFoam/createFields.H +++ b/applications/solvers/heatTransfer/buoyantFoam/createFields.H @@ -57,6 +57,7 @@ Info<< "Calculating field g.h\n" << endl; volScalarField gh("gh", g & mesh.C()); + surfaceScalarField ghf("gh", g & mesh.Cf()); dimensionedScalar pRef("pRef", p.dimensions(), thermo->lookup("pRef")); diff --git a/applications/solvers/heatTransfer/buoyantFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantFoam/pEqn.H index 36378f7f6f..9feededcd3 100644 --- a/applications/solvers/heatTransfer/buoyantFoam/pEqn.H +++ b/applications/solvers/heatTransfer/buoyantFoam/pEqn.H @@ -1,60 +1,67 @@ -bool closedVolume = pd.needReference(); - -rho = thermo->rho(); - -volScalarField rUA = 1.0/UEqn.A(); -U = rUA*UEqn.H(); - -phi = - fvc::interpolate(rho) - *( - (fvc::interpolate(U) & mesh.Sf()) - + fvc::ddtPhiCorr(rUA, rho, U, phi) - ) - - fvc::interpolate(rho*rUA*gh)*fvc::snGrad(rho)*mesh.magSf(); - -for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) { - fvScalarMatrix pdEqn + bool closedVolume = pd.needReference(); + + rho = thermo->rho(); + + volScalarField rUA = 1.0/UEqn.A(); + surfaceScalarField rhorUAf = fvc::interpolate(rho*rUA); + + U = rUA*UEqn.H(); + + surfaceScalarField phiU ( - fvm::ddt(psi, pd) - + fvc::ddt(psi)*pRef - + fvc::ddt(psi, rho)*gh - + fvc::div(phi) - - fvm::laplacian(rho*rUA, pd) + fvc::interpolate(rho) + *( + (fvc::interpolate(U) & mesh.Sf()) + + fvc::ddtPhiCorr(rUA, rho, U, phi) + ) ); - if (corr == nCorr-1 && nonOrth == nNonOrthCorr) + phi = phiU - ghf*fvc::snGrad(rho)*rhorUAf*mesh.magSf(); + + for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) { - pdEqn.solve(mesh.solver(pd.name() + "Final")); - } - else - { - pdEqn.solve(mesh.solver(pd.name())); + fvScalarMatrix pdEqn + ( + fvm::ddt(psi, pd) + + fvc::ddt(psi)*pRef + + fvc::ddt(psi, rho)*gh + + fvc::div(phi) + - fvm::laplacian(rhorUAf, pd) + ); + + if (corr == nCorr-1 && nonOrth == nNonOrthCorr) + { + pdEqn.solve(mesh.solver(pd.name() + "Final")); + } + else + { + pdEqn.solve(mesh.solver(pd.name())); + } + + if (nonOrth == nNonOrthCorr) + { + phi += pdEqn.flux(); + } } - if (nonOrth == nNonOrthCorr) + U += rUA*fvc::reconstruct((phi - phiU)/rhorUAf); + U.correctBoundaryConditions(); + + p == pd + rho*gh + pRef; + dpdt = fvc::ddt(p); + + #include "rhoEqn.H" + #include "compressibleContinuityErrs.H" + + // For closed-volume cases adjust the pressure and density levels + // to obey overall mass continuity + if (closedVolume) { - phi += pdEqn.flux(); + p += (initialMass - fvc::domainIntegrate(thermo->psi()*p)) + /fvc::domainIntegrate(thermo->psi()); + rho = thermo->rho(); } -} -p == pd + rho*gh + pRef; -dpdt = fvc::ddt(p); - -#include "rhoEqn.H" -#include "compressibleContinuityErrs.H" - -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) -{ - p += (initialMass - fvc::domainIntegrate(thermo->psi()*p)) - /fvc::domainIntegrate(thermo->psi()); pd == p - (rho*gh + pRef); - rho = thermo->rho(); } diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/UEqn.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/UEqn.H index 5bc4c4738c..71a57cd2a9 100644 --- a/applications/solvers/heatTransfer/buoyantSimpleFoam/UEqn.H +++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/UEqn.H @@ -11,7 +11,15 @@ eqnResidual = solve ( - UEqn() == -fvc::grad(pd) - fvc::grad(rho)*gh + UEqn() + == + -fvc::reconstruct + ( + ( + fvc::snGrad(pd) + + ghf*fvc::snGrad(rho) + ) * mesh.magSf() + ) ).initialResidual(); maxResidual = max(eqnResidual, maxResidual); diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H index 184242be81..d26da5cb32 100644 --- a/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H +++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H @@ -53,6 +53,7 @@ Info<< "Calculating field g.h\n" << endl; volScalarField gh("gh", g & mesh.C()); + surfaceScalarField ghf("gh", g & mesh.Cf()); dimensionedScalar pRef("pRef", p.dimensions(), thermo->lookup("pRef")); diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H index 8c3621205b..3c257d249d 100644 --- a/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H +++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H @@ -1,55 +1,65 @@ -volScalarField rUA = 1.0/UEqn().A(); -U = rUA*UEqn().H(); -UEqn.clear(); - -phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf()); -bool closedVolume = adjustPhi(phi, U, p); -phi -= fvc::interpolate(rho*gh*rUA)*fvc::snGrad(rho)*mesh.magSf(); - -for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) { - fvScalarMatrix pdEqn - ( - fvm::laplacian(rho*rUA, pd) == fvc::div(phi) - ); + volScalarField rUA = 1.0/UEqn().A(); + surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA)); - pdEqn.setReference(pdRefCell, pdRefValue); - // retain the residual from the first iteration - if (nonOrth == 0) + U = rUA*UEqn().H(); + UEqn.clear(); + + phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf()); + bool closedVolume = adjustPhi(phi, U, p); + surfaceScalarField buoyancyPhi = ghf*fvc::snGrad(rho)*rhorUAf*mesh.magSf(); + phi -= buoyancyPhi; + + for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) { - eqnResidual = pdEqn.solve().initialResidual(); - maxResidual = max(eqnResidual, maxResidual); - } - else - { - pdEqn.solve(); + fvScalarMatrix pdEqn + ( + fvm::laplacian(rhorUAf, pd) == fvc::div(phi) + ); + + pdEqn.setReference(pdRefCell, pdRefValue); + + // retain the residual from the first iteration + if (nonOrth == 0) + { + eqnResidual = pdEqn.solve().initialResidual(); + maxResidual = max(eqnResidual, maxResidual); + } + else + { + pdEqn.solve(); + } + + if (nonOrth == nNonOrthCorr) + { + // Calculate the conservative fluxes + phi -= pdEqn.flux(); + + // Explicitly relax pressure for momentum corrector + pd.relax(); + + // Correct the momentum source with the pressure gradient flux + // calculated from the relaxed pressure + U -= rUA*fvc::reconstruct((buoyancyPhi + pdEqn.flux())/rhorUAf); + U.correctBoundaryConditions(); + } } - if (nonOrth == nNonOrthCorr) + #include "continuityErrs.H" + + p == pd + rho*gh + pRef; + + // For closed-volume cases adjust the pressure and density levels + // to obey overall mass continuity + if (closedVolume) { - phi -= pdEqn.flux(); + p += (initialMass - fvc::domainIntegrate(thermo->psi()*p)) + /fvc::domainIntegrate(thermo->psi()); } -} -#include "continuityErrs.H" + rho = thermo->rho(); + rho.relax(); + Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl; -// Explicitly relax pressure for momentum corrector -pd.relax(); - -p = pd + rho*gh + pRef; - -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) -{ - p += (initialMass - fvc::domainIntegrate(thermo->psi()*p)) - /fvc::domainIntegrate(thermo->psi()); pd == p - (rho*gh + pRef); } - -rho = thermo->rho(); -rho.relax(); -Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl; diff --git a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/Make/options b/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/Make/options index 22d41fdb80..89d7787af6 100644 --- a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/Make/options +++ b/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/Make/options @@ -1,4 +1,5 @@ EXE_INC = \ + -I../buoyantSimpleFoam \ -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \ -I$(LIB_SRC)/turbulenceModels \ diff --git a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/UEqn.H b/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/UEqn.H deleted file mode 100644 index 5bc4c4738c..0000000000 --- a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/UEqn.H +++ /dev/null @@ -1,18 +0,0 @@ - // Solve the Momentum equation - - tmp UEqn - ( - fvm::div(phi, U) - - fvm::Sp(fvc::div(phi), U) - + turbulence->divDevRhoReff(U) - ); - - UEqn().relax(); - - eqnResidual = solve - ( - UEqn() == -fvc::grad(pd) - fvc::grad(rho)*gh - ).initialResidual(); - - maxResidual = max(eqnResidual, maxResidual); - diff --git a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/createFields.H b/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/createFields.H index 76e3805bca..62c06ec38d 100644 --- a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/createFields.H +++ b/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/createFields.H @@ -54,6 +54,7 @@ Info<< "Calculating field g.h\n" << endl; volScalarField gh("gh", g & mesh.C()); + surfaceScalarField ghf("gh", g & mesh.Cf()); dimensionedScalar pRef("pRef", p.dimensions(), thermo->lookup("pRef")); diff --git a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/pEqn.H deleted file mode 100644 index 2a713bac62..0000000000 --- a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/pEqn.H +++ /dev/null @@ -1,54 +0,0 @@ -volScalarField rUA = 1.0/UEqn().A(); -U = rUA*UEqn().H(); -UEqn.clear(); -phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf()); -bool closedVolume = adjustPhi(phi, U, p); -phi -= fvc::interpolate(rho*gh*rUA)*fvc::snGrad(rho)*mesh.magSf(); - -for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) -{ - fvScalarMatrix pdEqn - ( - fvm::laplacian(rho*rUA, pd) == fvc::div(phi) - ); - - pdEqn.setReference(pdRefCell, pdRefValue); - // retain the residual from the first iteration - if (nonOrth == 0) - { - eqnResidual = pdEqn.solve().initialResidual(); - maxResidual = max(eqnResidual, maxResidual); - } - else - { - pdEqn.solve(); - } - - if (nonOrth == nNonOrthCorr) - { - phi -= pdEqn.flux(); - } -} - -#include "continuityErrs.H" - -// Explicitly relax pressure for momentum corrector -pd.relax(); - -p = pd + rho*gh + pRef; - -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) -{ - p += (initialMass - fvc::domainIntegrate(thermo->psi()*p)) - /fvc::domainIntegrate(thermo->psi()); - pd == p - (rho*gh + pRef); -} - -rho = thermo->rho(); -rho.relax(); -Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl; From 6d3a27237230af18f240164916e1fb61f0abeade Mon Sep 17 00:00:00 2001 From: andy Date: Tue, 3 Feb 2009 18:24:23 +0000 Subject: [PATCH 4/8] applying improved flux-velocity correspondence --- .../buoyantBoussinesqSimpleFoam/UEqn.H | 11 ++- .../createFields.H | 4 +- .../buoyantBoussinesqSimpleFoam/pdEqn.H | 76 ++++++++++--------- .../writeAdditionalFields.H | 2 +- 4 files changed, 53 insertions(+), 40 deletions(-) diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/UEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/UEqn.H index 4b6887fff0..cf5c4a1978 100644 --- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/UEqn.H +++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/UEqn.H @@ -1,4 +1,4 @@ - // Solve the Momentum equation + // Solve the momentum equation tmp UEqn ( @@ -13,8 +13,13 @@ ( UEqn() == - -fvc::grad(pd) - + beta*gh*fvc::grad(T) + -fvc::reconstruct + ( + ( + fvc::snGrad(pd) + - betaghf*fvc::snGrad(T) + ) * mesh.magSf() + ) ).initialResidual(); maxResidual = max(eqnResidual, maxResidual); diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H index 03a5be6fd6..5f3f13626d 100644 --- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H +++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H @@ -53,8 +53,8 @@ incompressible::RASModel::New(U, phi, laminarTransport) ); - Info<< "Calculating field g.h\n" << endl; - volScalarField gh("gh", g & mesh.C()); + Info<< "Calculating field beta*(g.h)\n" << endl; + surfaceScalarField betaghf("betagh", beta*(g & mesh.Cf())); label pdRefCell = 0; scalar pdRefValue = 0.0; diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pdEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pdEqn.H index edf3e26d11..e782d26c8d 100644 --- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pdEqn.H +++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pdEqn.H @@ -1,41 +1,49 @@ -volScalarField rUA = 1.0/UEqn().A(); -U = rUA*UEqn().H(); -UEqn.clear(); - -phi = fvc::interpolate(U) & mesh.Sf(); -adjustPhi(phi, U, pd); -phi += fvc::interpolate(beta*gh*rUA)*fvc::snGrad(T)*mesh.magSf(); - -for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) { - fvScalarMatrix pdEqn - ( - fvm::laplacian(rUA, pd) == fvc::div(phi) - ); + volScalarField rUA("rUA", 1.0/UEqn().A()); + surfaceScalarField rUAf("(1|A(U))", fvc::interpolate(rUA)); - pdEqn.setReference(pdRefCell, pdRefValue); + U = rUA*UEqn().H(); + UEqn.clear(); - // retain the residual from the first iteration - if (nonOrth == 0) + phi = fvc::interpolate(U) & mesh.Sf(); + adjustPhi(phi, U, pd); + surfaceScalarField buoyancyPhi = -betaghf*fvc::snGrad(T)*rUAf*mesh.magSf(); + phi -= buoyancyPhi; + + for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) { - eqnResidual = pdEqn.solve().initialResidual(); - maxResidual = max(eqnResidual, maxResidual); - } - else - { - pdEqn.solve(); + fvScalarMatrix pdEqn + ( + fvm::laplacian(rUAf, pd) == fvc::div(phi) + ); + + pdEqn.setReference(pdRefCell, pdRefValue); + + // retain the residual from the first iteration + if (nonOrth == 0) + { + eqnResidual = pdEqn.solve().initialResidual(); + maxResidual = max(eqnResidual, maxResidual); + } + else + { + pdEqn.solve(); + } + + if (nonOrth == nNonOrthCorr) + { + // Calculate the conservative fluxes + phi -= pdEqn.flux(); + + // Explicitly relax pressure for momentum corrector + pd.relax(); + + // Correct the momentum source with the pressure gradient flux + // calculated from the relaxed pressure + U -= rUA*fvc::reconstruct((buoyancyPhi + pdEqn.flux())/rUAf); + U.correctBoundaryConditions(); + } } - if (nonOrth == nNonOrthCorr) - { - phi -= pdEqn.flux(); - } + #include "continuityErrs.H" } - -#include "continuityErrs.H" - -// Explicitly relax pressure for momentum corrector -pd.relax(); - -U -= rUA*(fvc::grad(pd) - beta*gh*fvc::grad(T)); -U.correctBoundaryConditions(); diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/writeAdditionalFields.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/writeAdditionalFields.H index f4b9caf780..20f7c6ae1d 100644 --- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/writeAdditionalFields.H +++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/writeAdditionalFields.H @@ -23,7 +23,7 @@ IOobject::NO_READ, IOobject::AUTO_WRITE ), - pd + rhoEff*gh + pRef + pd + rhoEff*(g & mesh.C()) + pRef ); p.write(); } From 61e1c0c1ea97b600202cd437f859b71cf428f967 Mon Sep 17 00:00:00 2001 From: andy Date: Wed, 4 Feb 2009 09:46:26 +0000 Subject: [PATCH 5/8] major clean-up --- .../chtMultiRegionFoam/Make/files | 2 + .../chtMultiRegionFoam/chtMultiRegionFoam.C | 11 +-- .../chtMultiRegionFoam/fluid/UEqn.H | 20 +++-- .../fluid/compressibleContinuityErrors.C | 58 -------------- .../fluid/compressibleCourantNo.C | 11 ++- .../fluid/createFluidFields.H | 6 +- .../chtMultiRegionFoam/fluid/hEqn.H | 20 +++-- .../chtMultiRegionFoam/fluid/pEqn.H | 78 +++++++++++-------- .../chtMultiRegionFoam/fluid/pdEqn.H | 14 ---- .../chtMultiRegionFoam/fluid/rhoEqn.H | 2 +- .../fluid/setInitialDeltaT.H | 15 ---- .../fluid/solveContinuityEquation.C | 37 --------- .../fluid/solveEnthalpyEquation.C | 56 ------------- .../chtMultiRegionFoam/fluid/solveFluid.H | 2 +- .../fluid/solveMomentumEquation.C | 57 -------------- .../fluid/solvePressureDifferenceEquation.C | 73 ----------------- .../chtMultiRegionFoam/solid/solveSolid.H | 5 +- 17 files changed, 87 insertions(+), 380 deletions(-) delete mode 100644 applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleContinuityErrors.C delete mode 100644 applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pdEqn.H delete mode 100644 applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setInitialDeltaT.H delete mode 100644 applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveContinuityEquation.C delete mode 100644 applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveEnthalpyEquation.C delete mode 100644 applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveMomentumEquation.C delete mode 100644 applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solvePressureDifferenceEquation.C diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/files b/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/files index 3ee71c0ea6..9d9152930a 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/files +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/files @@ -7,6 +7,8 @@ derivedFvPatchFields/solidWallHeatFluxTemperatureCoupled/solidWallHeatFluxTemper derivedFvPatchFields/solidWallTemperatureCoupled/solidWallTemperatureCoupledFvPatchScalarField.C derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.C +fluid/compressibleCourantNo.C + chtMultiRegionFoam.C EXE = $(FOAM_APPBIN)/chtMultiRegionFoam diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C index dad472930c..84c7c11806 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C @@ -36,16 +36,10 @@ Description #include "turbulenceModel.H" #include "fixedGradientFvPatchFields.H" #include "regionProperties.H" +#include "compressibleCourantNo.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -#include "solveContinuityEquation.C" -#include "solveMomentumEquation.C" -#include "compressibleContinuityErrors.C" -#include "solvePressureDifferenceEquation.C" -#include "solveEnthalpyEquation.C" -#include "compressibleCourantNo.C" - int main(int argc, char *argv[]) { @@ -58,7 +52,6 @@ int main(int argc, char *argv[]) # include "createSolidMeshes.H" # include "createFluidFields.H" - # include "createSolidFields.H" # include "initContinuityErrs.H" @@ -89,6 +82,7 @@ int main(int argc, char *argv[]) { Info<< "\nSolving for fluid region " << fluidRegions[i].name() << endl; +# include "setRegionFluidFields.H" # include "readFluidMultiRegionPISOControls.H" # include "solveFluid.H" } @@ -97,6 +91,7 @@ int main(int argc, char *argv[]) { Info<< "\nSolving for solid region " << solidRegions[i].name() << endl; +# include "setRegionSolidFields.H" # include "readSolidMultiRegionPISOControls.H" # include "solveSolid.H" } diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/UEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/UEqn.H index ba3d099a9f..16be7c1ed2 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/UEqn.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/UEqn.H @@ -1,10 +1,14 @@ - tmp UEqn = solveMomentumEquation + // Solve the Momentum equation + tmp UEqn ( - momentumPredictor, - Uf[i], - rhof[i], - phif[i], - pdf[i], - ghf[i], - turb[i] + fvm::ddt(rho, U) + + fvm::div(phi, U) + + turb.divDevRhoReff(U) ); + + UEqn().relax(); + + if (momentumPredictor) + { + solve(UEqn() == -fvc::grad(pd) - fvc::grad(rho)*gh); + } diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleContinuityErrors.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleContinuityErrors.C deleted file mode 100644 index 9cf6446aa0..0000000000 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleContinuityErrors.C +++ /dev/null @@ -1,58 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. - \\/ M anipulation | -------------------------------------------------------------------------------- -License - This file is part of OpenFOAM. - - OpenFOAM is free software; you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation; either version 2 of the License, or (at your - option) any later version. - - OpenFOAM is distributed in the hope that it will be useful, but WITHOUT - ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or - FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License - for more details. - - You should have received a copy of the GNU General Public License - along with OpenFOAM; if not, write to the Free Software Foundation, - Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -Description - Continuity errors for fluid meshes - -\*---------------------------------------------------------------------------*/ - -void compressibleContinuityErrors -( - scalar& cumulativeContErr, - const volScalarField& rho, - const basicThermo& thermo -) -{ - dimensionedScalar totalMass = fvc::domainIntegrate(rho); - - scalar sumLocalContErr = - ( - fvc::domainIntegrate(mag(rho - thermo.rho()))/totalMass - ).value(); - - scalar globalContErr = - ( - fvc::domainIntegrate(rho - thermo.rho())/totalMass - ).value(); - - cumulativeContErr += globalContErr; - - const word& regionName = rho.mesh().name(); - - Info<< "time step continuity errors (" << regionName << ")" - << ": sum local = " << sumLocalContErr - << ", global = " << globalContErr - << ", cumulative = " << cumulativeContErr - << endl; -} diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleCourantNo.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleCourantNo.C index 487fd93825..e671a256c9 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleCourantNo.C +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleCourantNo.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. + \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -22,13 +22,12 @@ License along with OpenFOAM; if not, write to the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA -Description - Calculates and outputs the mean and maximum Courant Numbers for the fluid - regions - \*---------------------------------------------------------------------------*/ -scalar compressibleCourantNo +#include "compressibleCourantNo.H" +#include "fvc.H" + +Foam::scalar Foam::compressibleCourantNo ( const fvMesh& mesh, const Time& runTime, diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H index 314b9d028a..a1fe9e2bc1 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H @@ -4,7 +4,7 @@ PtrList Kf(fluidRegions.size()); PtrList Uf(fluidRegions.size()); PtrList phif(fluidRegions.size()); - PtrList turb(fluidRegions.size()); + PtrList turbulence(fluidRegions.size()); PtrList DpDtf(fluidRegions.size()); PtrList ghf(fluidRegions.size()); PtrList pdf(fluidRegions.size()); @@ -123,8 +123,8 @@ ) ); - Info<< " Adding to turb\n" << endl; - turb.set + Info<< " Adding to turbulence\n" << endl; + turbulence.set ( i, autoPtr diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/hEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/hEqn.H index 301ceddfdb..d421649f13 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/hEqn.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/hEqn.H @@ -1,9 +1,17 @@ - solveEnthalpyEquation +{ + tmp hEqn ( - rhof[i], - DpDtf[i], - phif[i], - turb[i], - thermof[i] + fvm::ddt(rho, h) + + fvm::div(phi, h) + - fvm::laplacian(turb.alphaEff(), h) + == + DpDt ); + hEqn().relax(); + hEqn().solve(); + thermo.correct(); + + Info<< "Min/max T:" << min(thermo.T()) << ' ' << max(thermo.T()) + << endl; +} diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H index 2b1f5fceb1..54507f0d19 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H @@ -1,60 +1,70 @@ { bool closedVolume = false; - rhof[i] = thermof[i].rho(); + rho = thermo.rho(); volScalarField rUA = 1.0/UEqn().A(); - Uf[i] = rUA*UEqn().H(); + U = rUA*UEqn().H(); - phif[i] = - fvc::interpolate(rhof[i]) + phi = + fvc::interpolate(rho) *( - (fvc::interpolate(Uf[i]) & fluidRegions[i].Sf()) - + fvc::ddtPhiCorr(rUA, rhof[i], Uf[i], phif[i]) + (fvc::interpolate(U) & mesh.Sf()) + + fvc::ddtPhiCorr(rUA, rho, U, phi) ) - - fvc::interpolate(rhof[i]*rUA*ghf[i]) - *fvc::snGrad(rhof[i]) - *fluidRegions[i].magSf(); + - fvc::interpolate(rho*rUA*gh)*fvc::snGrad(rho)*mesh.magSf(); - // Solve pressure difference -# include "pdEqn.H" + for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) + { + fvScalarMatrix pdEqn + ( + fvc::ddt(rho) + + fvc::div(phi) + - fvm::laplacian(rho*rUA, pd) + ); + + if (corr == nCorr-1 && nonOrth == nNonOrthCorr) + { + pdEqn.solve(mesh.solver(pd.name() + "Final")); + } + else + { + pdEqn.solve(mesh.solver(pd.name())); + } + + if (nonOrth == nNonOrthCorr) + { + phi += pdEqn.flux(); + } + } + + // Update pressure field (including bc) + thermo.p() == pd + rho*gh + pRef; + DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p); // Solve continuity # include "rhoEqn.H" - // Update pressure field (including bc) - thermof[i].p() == pdf[i] + rhof[i]*ghf[i] + pRef; - DpDtf[i] = fvc::DDt - ( - surfaceScalarField("phiU", phif[i]/fvc::interpolate(rhof[i])), - thermof[i].p() - ); - // Update continuity errors - compressibleContinuityErrors(cumulativeContErr, rhof[i], thermof[i]); +# include "compressibleContinuityErrors.H" // Correct velocity field - Uf[i] -= rUA*(fvc::grad(pdf[i]) + fvc::grad(rhof[i])*ghf[i]); - Uf[i].correctBoundaryConditions(); + 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) { - thermof[i].p() += + thermo.p() += ( - dimensionedScalar - ( - "massIni", - dimMass, - initialMassf[i] - ) - - fvc::domainIntegrate(thermof[i].psi()*thermof[i].p()) - )/fvc::domainIntegrate(thermof[i].psi()); - pdf[i] == thermof[i].p() - (rhof[i]*ghf[i] + pRef); - rhof[i] = thermof[i].rho(); + dimensionedScalar("massIni", dimMass, initialMassf[i]) + - fvc::domainIntegrate(thermo.psi()*thermo.p()) + )/fvc::domainIntegrate(thermo.psi()); + pd == thermo.p() - (rho*gh + pRef); + rho = thermo.rho(); } // Update thermal conductivity - Kf[i] = thermof[i].Cp()*turb[i].alphaEff(); + K = thermo.Cp()*turb.alphaEff(); } diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pdEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pdEqn.H deleted file mode 100644 index a393843256..0000000000 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pdEqn.H +++ /dev/null @@ -1,14 +0,0 @@ - solvePressureDifferenceEquation - ( - corr, - nCorr, - nNonOrthCorr, - closedVolume, - pdf[i], - pRef, - rhof[i], - thermof[i].psi(), - rUA, - ghf[i], - phif[i] - ); diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/rhoEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/rhoEqn.H index a7e2a98c4a..a6b0ac9fe1 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/rhoEqn.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/rhoEqn.H @@ -1 +1 @@ - solveContinuityEquation(rhof[i], phif[i]); + solve(fvm::ddt(rho) + fvc::div(phi)); diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setInitialDeltaT.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setInitialDeltaT.H deleted file mode 100644 index 161eb742e7..0000000000 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setInitialDeltaT.H +++ /dev/null @@ -1,15 +0,0 @@ - if (adjustTimeStep) - { - if (CoNum > SMALL) - { - runTime.setDeltaT - ( - min - ( - maxCo*runTime.deltaT().value()/CoNum, - maxDeltaT - ) - ); - } - } - diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveContinuityEquation.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveContinuityEquation.C deleted file mode 100644 index 5fe2090341..0000000000 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveContinuityEquation.C +++ /dev/null @@ -1,37 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. - \\/ M anipulation | -------------------------------------------------------------------------------- -License - This file is part of OpenFOAM. - - OpenFOAM is free software; you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation; either version 2 of the License, or (at your - option) any later version. - - OpenFOAM is distributed in the hope that it will be useful, but WITHOUT - ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or - FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License - for more details. - - You should have received a copy of the GNU General Public License - along with OpenFOAM; if not, write to the Free Software Foundation, - Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -Description - Solve continuity equation - -\*---------------------------------------------------------------------------*/ - -void solveContinuityEquation -( - volScalarField& rho, - const surfaceScalarField& phi -) -{ - solve(fvm::ddt(rho) + fvc::div(phi)); -} diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveEnthalpyEquation.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveEnthalpyEquation.C deleted file mode 100644 index 2e9cae95ad..0000000000 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveEnthalpyEquation.C +++ /dev/null @@ -1,56 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. - \\/ M anipulation | -------------------------------------------------------------------------------- -License - This file is part of OpenFOAM. - - OpenFOAM is free software; you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation; either version 2 of the License, or (at your - option) any later version. - - OpenFOAM is distributed in the hope that it will be useful, but WITHOUT - ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or - FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License - for more details. - - You should have received a copy of the GNU General Public License - along with OpenFOAM; if not, write to the Free Software Foundation, - Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -Description - Solve enthalpy equation - -\*---------------------------------------------------------------------------*/ - -void solveEnthalpyEquation -( - const volScalarField& rho, - const volScalarField& DpDt, - const surfaceScalarField& phi, - const compressible::turbulenceModel& turb, - basicThermo& thermo -) -{ - volScalarField& h = thermo.h(); - - tmp hEqn - ( - fvm::ddt(rho, h) - + fvm::div(phi, h) - - fvm::laplacian(turb.alphaEff(), h) - == - DpDt - ); - hEqn().relax(); - hEqn().solve(); - - thermo.correct(); - - Info<< "Min/max T:" << min(thermo.T()) << ' ' << max(thermo.T()) - << endl; -} diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveFluid.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveFluid.H index e24771256f..19ec50cac2 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveFluid.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveFluid.H @@ -12,4 +12,4 @@ # include "pEqn.H" } } - turb[i].correct(); + turb.correct(); diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveMomentumEquation.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveMomentumEquation.C deleted file mode 100644 index 935a02b60a..0000000000 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveMomentumEquation.C +++ /dev/null @@ -1,57 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. - \\/ M anipulation | -------------------------------------------------------------------------------- -License - This file is part of OpenFOAM. - - OpenFOAM is free software; you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation; either version 2 of the License, or (at your - option) any later version. - - OpenFOAM is distributed in the hope that it will be useful, but WITHOUT - ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or - FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License - for more details. - - You should have received a copy of the GNU General Public License - along with OpenFOAM; if not, write to the Free Software Foundation, - Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -Description - Solve momentum equation and return matrix for use in pressure equation - -\*---------------------------------------------------------------------------*/ - -tmp solveMomentumEquation -( - const bool momentumPredictor, - volVectorField& U, - const volScalarField& rho, - const surfaceScalarField& phi, - const volScalarField& pd, - const volScalarField& gh, - const compressible::turbulenceModel& turb -) -{ - // Solve the Momentum equation - tmp UEqn - ( - fvm::ddt(rho, U) - + fvm::div(phi, U) - + turb.divDevRhoReff(U) - ); - - UEqn().relax(); - - if (momentumPredictor) - { - solve(UEqn() == -fvc::grad(pd) - fvc::grad(rho)*gh); - } - - return UEqn; -} diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solvePressureDifferenceEquation.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solvePressureDifferenceEquation.C deleted file mode 100644 index 6d8843ab42..0000000000 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solvePressureDifferenceEquation.C +++ /dev/null @@ -1,73 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. - \\/ M anipulation | -------------------------------------------------------------------------------- -License - This file is part of OpenFOAM. - - OpenFOAM is free software; you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation; either version 2 of the License, or (at your - option) any later version. - - OpenFOAM is distributed in the hope that it will be useful, but WITHOUT - ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or - FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License - for more details. - - You should have received a copy of the GNU General Public License - along with OpenFOAM; if not, write to the Free Software Foundation, - Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -Description - Solve pressure difference equation - -\*---------------------------------------------------------------------------*/ - -void solvePressureDifferenceEquation -( - const label corr, - const label nCorr, - const label nNonOrthCorr, - bool& closedVolume, - volScalarField& pd, - const dimensionedScalar& pRef, - const volScalarField& rho, - const volScalarField& psi, - const volScalarField& rUA, - const volScalarField& gh, - surfaceScalarField& phi -) -{ - closedVolume = pd.needReference(); - - for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) - { - fvScalarMatrix pdEqn - ( - fvm::ddt(psi, pd) - + fvc::ddt(psi)*pRef - + fvc::ddt(psi, rho)*gh - + fvc::div(phi) - - fvm::laplacian(rho*rUA, pd) - ); - - //pdEqn.solve(); - if (corr == nCorr-1 && nonOrth == nNonOrthCorr) - { - pdEqn.solve(pd.mesh().solver(pd.name() + "Final")); - } - else - { - pdEqn.solve(pd.mesh().solver(pd.name())); - } - - if (nonOrth == nNonOrthCorr) - { - phi += pdEqn.flux(); - } - } -} diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H index 790d4ec934..5fa731824f 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H @@ -3,10 +3,9 @@ { solve ( - fvm::ddt(rhosCps[i], Ts[i]) - fvm::laplacian(Ks[i], Ts[i]) + fvm::ddt(rho*cp, T) - fvm::laplacian(K, T) ); } - Info<< "Min/max T:" << min(Ts[i]) << ' ' << max(Ts[i]) - << endl; + Info<< "Min/max T:" << min(T) << ' ' << max(T) << endl; } From efdb3fd3fa4d3d3be9a7bad6485ed9773c4f166e Mon Sep 17 00:00:00 2001 From: andy Date: Wed, 4 Feb 2009 10:21:23 +0000 Subject: [PATCH 6/8] improved flux calcs + more clean-up --- .../chtMultiRegionFoam/fluid/UEqn.H | 13 ++- .../fluid/compressibleMultiRegionCourantNo.H | 4 +- .../fluid/createFluidFields.H | 85 ++++++++++--------- .../chtMultiRegionFoam/fluid/pEqn.H | 37 ++++---- 4 files changed, 82 insertions(+), 57 deletions(-) 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); } From a83588ec11c90ca2f1b7be125e276549479f8ad3 Mon Sep 17 00:00:00 2001 From: andy Date: Wed, 4 Feb 2009 10:25:31 +0000 Subject: [PATCH 7/8] using surfaceScalarField constructor to name field --- applications/solvers/heatTransfer/buoyantFoam/pEqn.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/applications/solvers/heatTransfer/buoyantFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantFoam/pEqn.H index e88c69fab3..91e190b4cd 100644 --- a/applications/solvers/heatTransfer/buoyantFoam/pEqn.H +++ b/applications/solvers/heatTransfer/buoyantFoam/pEqn.H @@ -4,7 +4,7 @@ rho = thermo->rho(); volScalarField rUA = 1.0/UEqn.A(); - surfaceScalarField rhorUAf = fvc::interpolate(rho*rUA); + surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA)); U = rUA*UEqn.H(); From 6d1466465b729e5267706c622b7e61c13df8afb2 Mon Sep 17 00:00:00 2001 From: andy Date: Wed, 4 Feb 2009 11:27:35 +0000 Subject: [PATCH 8/8] consistency update --- applications/solvers/multiphase/interFoam/alphaEqn.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/applications/solvers/multiphase/interFoam/alphaEqn.H b/applications/solvers/multiphase/interFoam/alphaEqn.H index 3bf24a845e..0b2fb4ebf8 100644 --- a/applications/solvers/multiphase/interFoam/alphaEqn.H +++ b/applications/solvers/multiphase/interFoam/alphaEqn.H @@ -6,7 +6,7 @@ phic = min(interface.cAlpha()*phic, max(phic)); surfaceScalarField phir = phic*interface.nHatf(); - for (int gCorr=0; gCorr