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/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/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")); 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(); } 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 24ad89e4e3..75eb401d3d 100644 --- a/applications/solvers/heatTransfer/buoyantFoam/buoyantFoam.C +++ b/applications/solvers/heatTransfer/buoyantFoam/buoyantFoam.C @@ -41,42 +41,41 @@ 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 "readTimeControls.H" + #include "compressibleCourantNo.H" + #include "setInitialDeltaT.H" -# include "setRootCase.H" -# include "createTime.H" -# include "createMesh.H" -# include "readEnvironmentalProperties.H" -# include "createFields.H" -# include "initContinuityErrs.H" -# include "readTimeControls.H" -# include "compressibleCourantNo.H" -# include "setInitialDeltaT.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Info<< "\nStarting time loop\n" << endl; while (runTime.run()) { -# include "readTimeControls.H" -# 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 9535718fbb..a647f0a3ef 100644 --- a/applications/solvers/heatTransfer/buoyantFoam/createFields.H +++ b/applications/solvers/heatTransfer/buoyantFoam/createFields.H @@ -58,6 +58,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 1cdbd0b199..91e190b4cd 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("(rho*(1|A(U)))", 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(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), 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(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), 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; 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..e719d92433 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/UEqn.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/UEqn.H @@ -1,10 +1,25 @@ - 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::reconstruct + ( + ( + fvc::snGrad(pd) + + ghf*fvc::snGrad(rho) + ) * mesh.magSf() + ) + ); + } 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/compressibleContinuityErrors.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleContinuityErrors.H new file mode 100644 index 0000000000..046ca5ec37 --- /dev/null +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleContinuityErrors.H @@ -0,0 +1,21 @@ +{ + 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[i] += globalContErr; + + Info<< "time step continuity errors (" << mesh.name() << ")" + << ": sum local = " << sumLocalContErr + << ", global = " << globalContErr + << ", cumulative = " << cumulativeContErr[i] + << endl; +} diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleCourantNo.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleCourantNo.C index 487fd93825..5f04825a94 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleCourantNo.C +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleCourantNo.C @@ -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/solveContinuityEquation.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleCourantNo.H similarity index 73% rename from applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveContinuityEquation.C rename to applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleCourantNo.H index 5fe2090341..329b8cce65 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveContinuityEquation.C +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleCourantNo.H @@ -23,15 +23,27 @@ License Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA Description - Solve continuity equation + Calculates and outputs the mean and maximum Courant Numbers for the fluid + regions \*---------------------------------------------------------------------------*/ -void solveContinuityEquation -( - volScalarField& rho, - const surfaceScalarField& phi -) +#ifndef compressibleCourantNo_H +#define compressibleCourantNo_H + +#include "fvMesh.H" + +namespace Foam { - solve(fvm::ddt(rho) + fvc::div(phi)); + scalar compressibleCourantNo + ( + const fvMesh& mesh, + const Time& runTime, + const volScalarField& rho, + const surfaceScalarField& phi + ); } + +#endif + +// ************************************************************************* // 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 314b9d028a..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 turb(fluidRegions.size()); - PtrList DpDtf(fluidRegions.size()); - PtrList ghf(fluidRegions.size()); - PtrList pdf(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 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,29 +118,29 @@ IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE ), - linearInterpolate(rhof[i]*Uf[i]) + linearInterpolate(rhoFluid[i]*UFluid[i]) & fluidRegions[i].Sf() ) ); - Info<< " Adding to turb\n" << endl; - turb.set + Info<< " Adding to turbulence\n" << endl; + turbulence.set ( i, autoPtr ( 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/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/initContinuityErrs.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/initContinuityErrs.H new file mode 100644 index 0000000000..1a7f5a3262 --- /dev/null +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/initContinuityErrs.H @@ -0,0 +1 @@ +List cumulativeContErr(fluidRegions.size(), 0.0); diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H index 2b1f5fceb1..e297989809 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H @@ -1,60 +1,75 @@ { - bool closedVolume = false; + bool closedVolume = pd.needReference(); - rhof[i] = thermof[i].rho(); + rho = thermo.rho(); volScalarField rUA = 1.0/UEqn().A(); - Uf[i] = rUA*UEqn().H(); + surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA)); - phif[i] = - fvc::interpolate(rhof[i]) + U = rUA*UEqn().H(); + + surfaceScalarField phiU + ( + 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(); + ); - // Solve pressure difference -# include "pdEqn.H" + phi = phiU - ghf*fvc::snGrad(rho)*rhorUAf*mesh.magSf(); + + 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) + ); + + 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(); + } + } + + // Correct velocity field + U += rUA*fvc::reconstruct((phi - phiU)/rhorUAf); + U.correctBoundaryConditions(); + + // Update pressure field (including bc) + 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]); - - // Correct velocity field - Uf[i] -= rUA*(fvc::grad(pdf[i]) + fvc::grad(rhof[i])*ghf[i]); - Uf[i].correctBoundaryConditions(); +# include "compressibleContinuityErrors.H" // For closed-volume cases adjust the pressure and density levels // to obey overall mass continuity if (closedVolume) { - thermof[i].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(); + p += (massIni - fvc::domainIntegrate(psi*p))/fvc::domainIntegrate(psi); + rho = thermo.rho(); } // Update thermal conductivity - Kf[i] = thermof[i].Cp()*turb[i].alphaEff(); + K = thermoFluid[i].Cp()*turb.alphaEff(); + + // Update pd (including bc) + pd == p - (rho*gh + pRef); } 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/setRegionFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H new file mode 100644 index 0000000000..bc4590a428 --- /dev/null +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H @@ -0,0 +1,18 @@ + const fvMesh& mesh = fluidRegions[i]; + + basicThermo& thermo = thermoFluid[i]; + volScalarField& rho = rhoFluid[i]; + volScalarField& K = KFluid[i]; + volVectorField& U = UFluid[i]; + surfaceScalarField phi = phiFluid[i]; + compressible::turbulenceModel& turb = turbulence[i]; + volScalarField& DpDt = DpDtFluid[i]; + const volScalarField& gh = ghFluid[i]; + const surfaceScalarField& ghf = ghfFluid[i]; + volScalarField& pd = pdFluid[i]; + + volScalarField& p = thermo.p(); + const volScalarField& psi = thermo.psi(); + volScalarField& h = thermo.h(); + + const dimensionedScalar massIni("massIni", dimMass, initialMassFluid[i]); 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/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/setRegionSolidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/setRegionSolidFields.H new file mode 100644 index 0000000000..04c90d2c4c --- /dev/null +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/setRegionSolidFields.H @@ -0,0 +1,6 @@ +// fvMesh& mesh = solidRegions[i]; + + volScalarField& rho = rhos[i]; + volScalarField& cp = cps[i]; + volScalarField& K = Ks[i]; + volScalarField& T = Ts[i]; 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; } 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 phiPtr_; public: @@ -133,12 +133,12 @@ public: const surfaceScalarField& phi() const { - return *phiPtr_; + return phiPtr_(); } surfaceScalarField& phi() { - return *phiPtr_; + return phiPtr_(); } }; diff --git a/applications/utilities/surface/surfaceMeshConvertTesting/Make/options b/applications/utilities/surface/surfaceMeshConvertTesting/Make/options index 5293dabe4c..3b91e457ea 100644 --- a/applications/utilities/surface/surfaceMeshConvertTesting/Make/options +++ b/applications/utilities/surface/surfaceMeshConvertTesting/Make/options @@ -2,4 +2,4 @@ EXE_INC = \ -I$(LIB_SRC)/triSurface/lnInclude \ -I$(LIB_SRC)/surfMesh/lnInclude -EXE_LIBS = -lsurfMesh +EXE_LIBS = -ltriSurface -lsurfMesh diff --git a/etc/cshrc b/etc/cshrc index 5952f6d3c5..bcda5b6318 100644 --- a/etc/cshrc +++ b/etc/cshrc @@ -158,7 +158,7 @@ case Linux: setenv WM_CXXFLAGS '-mabi=64 -fPIC' setenv WM_LDFLAGS '-mabi=64 -G0' setenv WM_MPLIB MPI - ;; + breaksw default: echo Unknown processor type `uname -m` for Linux breaksw diff --git a/src/Pstream/Allwmake b/src/Pstream/Allwmake index f43c2b1a30..cb2c43e843 100755 --- a/src/Pstream/Allwmake +++ b/src/Pstream/Allwmake @@ -6,17 +6,17 @@ wmake libso dummy case "$WM_MPLIB" in GAMMA) - wmake libso gamma - ;; + wmake libso gamma + ;; LAM | *MPI* ) - export WM_OPTIONS=${WM_OPTIONS}$WM_MPLIB - set +x - echo - echo "Note: ignore spurious warnings about missing mpicxx.h headers" - set -x - wmake libso mpi - ;; + WM_OPTIONS=${WM_OPTIONS}$WM_MPLIB + set +x + echo + echo "Note: ignore spurious warnings about missing mpicxx.h headers" + set -x + wmake libso mpi + ;; esac # ----------------------------------------------------------------- end-of-file diff --git a/src/meshTools/Make/files b/src/meshTools/Make/files index 98c7bf0640..a93cd58400 100644 --- a/src/meshTools/Make/files +++ b/src/meshTools/Make/files @@ -55,14 +55,15 @@ indexedOctree/treeDataTriSurface.C searchableSurface = searchableSurface $(searchableSurface)/distributedTriSurfaceMesh.C $(searchableSurface)/searchableBox.C +$(searchableSurface)/searchableCylinder.C $(searchableSurface)/searchablePlane.C $(searchableSurface)/searchablePlate.C $(searchableSurface)/searchableSphere.C $(searchableSurface)/searchableSurface.C $(searchableSurface)/searchableSurfaces.C $(searchableSurface)/searchableSurfacesQueries.C -$(searchableSurface)/triSurfaceMesh.C $(searchableSurface)/searchableSurfaceWithGaps.C +$(searchableSurface)/triSurfaceMesh.C topoSets = sets/topoSets $(topoSets)/cellSet.C diff --git a/src/meshTools/searchableSurface/searchableCylinder.C b/src/meshTools/searchableSurface/searchableCylinder.C new file mode 100644 index 0000000000..e2b8426f1f --- /dev/null +++ b/src/meshTools/searchableSurface/searchableCylinder.C @@ -0,0 +1,468 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 + +\*---------------------------------------------------------------------------*/ + +#include "searchableCylinder.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + +defineTypeNameAndDebug(searchableCylinder, 0); +addToRunTimeSelectionTable(searchableSurface, searchableCylinder, dict); + +} + + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +Foam::pointIndexHit Foam::searchableCylinder::findNearest +( + const point& sample, + const scalar nearestDistSqr +) const +{ + pointIndexHit info(false, sample, -1); + + vector v(sample - point1_); + + // Decompose sample-point1 into normal and parallel component + scalar parallel = (v & unitDir_); + + // Remove the parallel component + v -= parallel*unitDir_; + scalar magV = mag(v); + + if (parallel <= 0) + { + // nearest is at point1 end cap. Clip to radius. + if (magV < ROOTVSMALL) + { + info.setPoint(point1_); + } + else + { + //info.setPoint(point1_ + min(magV, radius_)*v/magV); + info.setPoint(point1_ + radius_*(v/magV)); + } + } + else if (parallel >= magDir_) + { + // nearest is at point2 + if (magV < ROOTVSMALL) + { + info.setPoint(point2_); + } + else + { + info.setPoint(point2_ + min(magV, radius_)*v/magV); + } + } + else + { + if (magV < ROOTVSMALL) + { + info.setPoint(point1_ + parallel*unitDir_); + } + else + { + // Project onto radius + info.setPoint(point1_ + parallel*unitDir_ + radius_*v/magV); + } + } + + if (magSqr(sample - info.rawPoint()) < nearestDistSqr) + { + info.setHit(); + info.setIndex(0); + } + + return info; +} + + +// From http://www.gamedev.net/community/forums/topic.asp?topic_id=467789 - +// intersection of cylinder with ray +void Foam::searchableCylinder::findLineAll +( + const point& start, + const point& end, + pointIndexHit& near, + pointIndexHit& far +) const +{ + near.setMiss(); + far.setMiss(); + + // Line as P = start + t*V + const vector V(end-start); + +//Pout<< "point1:" << point1_ << " point2:" << point2_ +// << " start:" << start << " end:" << end << endl; + + const vector x = (start-point1_) ^ unitDir_; + const vector y = V ^ unitDir_; + const scalar d = sqr(radius_); + + // Second order equation of the form a*t^2 + b*t + c + const scalar a = (y&y); + const scalar b = 2*(x&y); + const scalar c = (x&x)-d; + + const scalar disc = b*b-4*a*c; + +//Pout<< "a:" << a << " b:" << b << " c:" << c << " disc:" << disc +// << endl; + + if (disc < 0) + { + return; + } + else if (disc < ROOTVSMALL) + { + // Single solution + if (mag(a) > ROOTVSMALL) + { + scalar t = -b/(2*a); + if (t >= 0 && t <= 1) + { + near.setPoint(start + t*V); + near.setHit(); + near.setIndex(0); + } + } + } + else + { + if (mag(a) > ROOTVSMALL) + { + scalar sqrtDisc = sqrt(disc); + + scalar t1 = (-b + sqrtDisc)/2*a; + scalar t2 = (-b - sqrtDisc)/2*a; + + if (t1 < t2) + { + if (t1 >= 0 && t1 <= 1) + { + near.setPoint(start + t1*V); + near.setHit(); + near.setIndex(0); + } + if (t2 >= 0 && t2 <= 1) + { + far.setPoint(start + t2*V); + far.setHit(); + far.setIndex(0); + } + } + else + { + if (t2 >= 0 && t2 <= 1) + { + near.setPoint(start + t2*V); + near.setHit(); + near.setIndex(0); + } + if (t1 >= 0 && t1 <= 1) + { + far.setPoint(start + t1*V); + far.setHit(); + far.setIndex(0); + } + } + } + } +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::searchableCylinder::searchableCylinder +( + const IOobject& io, + const point& point1, + const point& point2, + const scalar radius +) +: + searchableSurface(io), + point1_(point1), + point2_(point2), + magDir_(mag(point2_-point1_)), + unitDir_((point2_-point1_)/magDir_), + radius_(radius) +{ + Pout<< "point1_:" << point1_ << endl; + Pout<< "point2_:" << point2_ << endl; + Pout<< "magDir_:" << magDir_ << endl; + Pout<< "unitDir_:" << unitDir_ << endl; + Pout<< "radius_:" << radius_ << endl; +} + + +Foam::searchableCylinder::searchableCylinder +( + const IOobject& io, + const dictionary& dict +) +: + searchableSurface(io), + point1_(dict.lookup("point1")), + point2_(dict.lookup("point2")), + magDir_(mag(point2_-point1_)), + unitDir_((point2_-point1_)/magDir_), + radius_(readScalar(dict.lookup("radius"))) +{ + Pout<< "point1_:" << point1_ << endl; + Pout<< "point2_:" << point2_ << endl; + Pout<< "magDir_:" << magDir_ << endl; + Pout<< "unitDir_:" << unitDir_ << endl; + Pout<< "radius_:" << radius_ << endl; +} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::searchableCylinder::~searchableCylinder() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +const Foam::wordList& Foam::searchableCylinder::regions() const +{ + if (regions_.empty()) + { + regions_.setSize(1); + regions_[0] = "region0"; + } + return regions_; +} + + +void Foam::searchableCylinder::findNearest +( + const pointField& samples, + const scalarField& nearestDistSqr, + List& info +) const +{ + info.setSize(samples.size()); + + forAll(samples, i) + { + info[i] = findNearest(samples[i], nearestDistSqr[i]); + } +} + + +void Foam::searchableCylinder::findLine +( + const pointField& start, + const pointField& end, + List& info +) const +{ + info.setSize(start.size()); + + pointIndexHit b; + + forAll(start, i) + { + // Pick nearest intersection. If none intersected take second one. + findLineAll(start[i], end[i], info[i], b); + if (!info[i].hit() && b.hit()) + { + info[i] = b; + } + } +} + + +void Foam::searchableCylinder::findLineAny +( + const pointField& start, + const pointField& end, + List& info +) const +{ + info.setSize(start.size()); + + pointIndexHit b; + + forAll(start, i) + { + // Discard far intersection + findLineAll(start[i], end[i], info[i], b); + if (!info[i].hit() && b.hit()) + { + info[i] = b; + } + } +} + + +void Foam::searchableCylinder::findLineAll +( + const pointField& start, + const pointField& end, + List >& info +) const +{ + info.setSize(start.size()); + + forAll(start, i) + { + pointIndexHit near, far; + findLineAll(start[i], end[i], near, far); + + if (near.hit()) + { + if (far.hit()) + { + info[i].setSize(2); + info[i][0] = near; + info[i][1] = far; + } + else + { + info[i].setSize(1); + info[i][0] = near; + } + } + else + { + if (far.hit()) + { + info[i].setSize(1); + info[i][0] = far; + } + else + { + info[i].clear(); + } + } + } +} + + +void Foam::searchableCylinder::getRegion +( + const List& info, + labelList& region +) const +{ + region.setSize(info.size()); + region = 0; +} + + +void Foam::searchableCylinder::getNormal +( + const List& info, + vectorField& normal +) const +{ + normal.setSize(info.size()); + normal = vector::zero; + + forAll(info, i) + { + if (info[i].hit()) + { + vector v(info[i].hitPoint() - point1_); + + // Decompose sample-point1 into normal and parallel component + scalar parallel = v & unitDir_; + + if (parallel < 0) + { + normal[i] = -unitDir_; + } + else if (parallel > magDir_) + { + normal[i] = -unitDir_; + } + else + { + // Remove the parallel component + v -= parallel*unitDir_; + normal[i] = v/mag(v); + } + } + } +} + + +void Foam::searchableCylinder::getVolumeType +( + const pointField& points, + List& volType +) const +{ + volType.setSize(points.size()); + volType = INSIDE; + + forAll(points, pointI) + { + const point& pt = points[pointI]; + + vector v(pt - point1_); + + // Decompose sample-point1 into normal and parallel component + scalar parallel = v & unitDir_; + + if (parallel < 0) + { + // left of point1 endcap + volType[pointI] = OUTSIDE; + } + else if (parallel > magDir_) + { + // right of point2 endcap + volType[pointI] = OUTSIDE; + } + else + { + // Remove the parallel component + v -= parallel*unitDir_; + + if (mag(v) > radius_) + { + volType[pointI] = OUTSIDE; + } + else + { + volType[pointI] = INSIDE; + } + } + } +} + + +// ************************************************************************* // diff --git a/src/meshTools/searchableSurface/searchableCylinder.H b/src/meshTools/searchableSurface/searchableCylinder.H new file mode 100644 index 0000000000..55ee0354dd --- /dev/null +++ b/src/meshTools/searchableSurface/searchableCylinder.H @@ -0,0 +1,225 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 + +Class + Foam::searchableCylinder + +Description + Searching on cylinder + +SourceFiles + searchableCylinder.C + +\*---------------------------------------------------------------------------*/ + +#ifndef searchableCylinder_H +#define searchableCylinder_H + +#include "treeBoundBox.H" +#include "searchableSurface.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// Forward declaration of classes + +/*---------------------------------------------------------------------------*\ + Class searchableCylinder Declaration +\*---------------------------------------------------------------------------*/ + +class searchableCylinder +: + public searchableSurface +{ +private: + + // Private Member Data + + //- 'left' point + const point point1_; + + //- 'right' point + const point point2_; + + //- length of vector point2-point1 + const scalar magDir_; + + //- normalised vector point2-point1 + const vector unitDir_; + + //- Radius squared + const scalar radius_; + + //- Names of regions + mutable wordList regions_; + + + // Private Member Functions + + //- Find nearest point on cylinder. + pointIndexHit findNearest + ( + const point& sample, + const scalar nearestDistSqr + ) const; + + //- Find intersection with cylinder + void findLineAll + ( + const point& start, + const point& end, + pointIndexHit& near, + pointIndexHit& far + ) const; + + + //- Disallow default bitwise copy construct + searchableCylinder(const searchableCylinder&); + + //- Disallow default bitwise assignment + void operator=(const searchableCylinder&); + + +public: + + //- Runtime type information + TypeName("searchableCylinder"); + + + // Constructors + + //- Construct from components + searchableCylinder + ( + const IOobject& io, + const point&, + const point&, + const scalar radius + ); + + //- Construct from dictionary (used by searchableSurface) + searchableCylinder + ( + const IOobject& io, + const dictionary& dict + ); + + // Destructor + + virtual ~searchableCylinder(); + + + // Member Functions + + virtual const wordList& regions() const; + + //- Whether supports volume type below + virtual bool hasVolumeType() const + { + return true; + } + + //- Range of local indices that can be returned. + virtual label size() const + { + return 1; + } + + + // Multiple point queries. + + virtual void findNearest + ( + const pointField& sample, + const scalarField& nearestDistSqr, + List& + ) const; + + virtual void findLine + ( + const pointField& start, + const pointField& end, + List& + ) const; + + virtual void findLineAny + ( + const pointField& start, + const pointField& end, + List& + ) const; + + //- Get all intersections in order from start to end. + virtual void findLineAll + ( + const pointField& start, + const pointField& end, + List >& + ) const; + + //- From a set of points and indices get the region + virtual void getRegion + ( + const List&, + labelList& region + ) const; + + //- From a set of points and indices get the normal + virtual void getNormal + ( + const List&, + vectorField& normal + ) const; + + //- Determine type (inside/outside/mixed) for point. unknown if + // cannot be determined (e.g. non-manifold surface) + virtual void getVolumeType + ( + const pointField&, + List& + ) const; + + + // regIOobject implementation + + bool writeData(Ostream&) const + { + notImplemented("searchableCylinder::writeData(Ostream&) const"); + return false; + } + +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/postProcessing/foamCalcFunctions/Make/files b/src/postProcessing/foamCalcFunctions/Make/files index 50dc36be0c..4f4a0ccf6f 100644 --- a/src/postProcessing/foamCalcFunctions/Make/files +++ b/src/postProcessing/foamCalcFunctions/Make/files @@ -7,6 +7,7 @@ field/magSqr/magSqr.C field/magGrad/magGrad.C field/div/div.C field/randomise/randomise.C +field/interpolate/interpolate.C basic/add/add.C diff --git a/src/postProcessing/foamCalcFunctions/field/interpolate/interpolate.C b/src/postProcessing/foamCalcFunctions/field/interpolate/interpolate.C new file mode 100644 index 0000000000..d9126c6c03 --- /dev/null +++ b/src/postProcessing/foamCalcFunctions/field/interpolate/interpolate.C @@ -0,0 +1,136 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 + +\*---------------------------------------------------------------------------*/ + +#include "interpolate.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + namespace calcTypes + { + defineTypeNameAndDebug(interpolate, 0); + addToRunTimeSelectionTable(calcType, interpolate, dictionary); + } +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::calcTypes::interpolate::interpolate() +: + calcType() +{} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::calcTypes::interpolate::~interpolate() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void Foam::calcTypes::interpolate::init() +{ + Foam::argList::validArgs.append("interpolate"); + argList::validArgs.append("fieldName1 .. fieldNameN"); +} + + +void Foam::calcTypes::interpolate::preCalc +( + const argList& args, + const Time& runTime, + const fvMesh& mesh +) +{ + if (args.additionalArgs().size() < 2) + { + Info<< nl << "must specify one or more fields" << nl; + args.printUsage(); + FatalError.exit(); + } +} + + +void Foam::calcTypes::interpolate::calc +( + const argList& args, + const Time& runTime, + const fvMesh& mesh +) +{ + const stringList& params = args.additionalArgs(); + + for (label fieldi=1; fieldi(fieldHeader, mesh, processed); + writeInterpolateField(fieldHeader, mesh, processed); + writeInterpolateField + ( + fieldHeader, + mesh, + processed + ); + writeInterpolateField(fieldHeader, mesh, processed); + writeInterpolateField(fieldHeader, mesh, processed); + + if (!processed) + { + FatalError + << "Unable to process " << fieldName << nl + << "No call to interpolate for fields of type " + << fieldHeader.headerClassName() << nl << nl + << exit(FatalError); + } + } + else + { + Info<< " No " << fieldName << endl; + } + } +} + + +// ************************************************************************* // + diff --git a/src/postProcessing/foamCalcFunctions/field/interpolate/interpolate.H b/src/postProcessing/foamCalcFunctions/field/interpolate/interpolate.H new file mode 100644 index 0000000000..156a2c263f --- /dev/null +++ b/src/postProcessing/foamCalcFunctions/field/interpolate/interpolate.H @@ -0,0 +1,138 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 + +Class + Foam::calcTypes::interpolate + +Description + Interpolates volume fields to surface fields for each time. + +SourceFiles + interpolate.C + +\*---------------------------------------------------------------------------*/ + +#ifndef interpolate_H +#define interpolate_H + +#include "calcType.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +namespace calcTypes +{ + +/*---------------------------------------------------------------------------*\ + Class interpolate Declaration +\*---------------------------------------------------------------------------*/ + +class interpolate +: + public calcType +{ + // Private Member Functions + + //- Disallow default bitwise copy construct + interpolate(const interpolate&); + + //- Disallow default bitwise assignment + void operator=(const interpolate&); + + +protected: + + // Member Functions + + // Calculation routines + + //- Initialise - typically setting static variables, + // e.g. command line arguments + virtual void init(); + + //- Pre-time loop calculations + virtual void preCalc + ( + const argList& args, + const Time& runTime, + const fvMesh& mesh + ); + + //- Time loop calculations + virtual void calc + ( + const argList& args, + const Time& runTime, + const fvMesh& mesh + ); + + + // I-O + + //- Write interpolate fields + template + void writeInterpolateField + ( + const IOobject& header, + const fvMesh& mesh, + bool& processed + ); + + +public: + + //- Runtime type information + TypeName("interpolate"); + + + // Constructors + + //- Construct null + interpolate(); + + + // Destructor + + virtual ~interpolate(); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace calcTypes +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "writeInterpolateField.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveMomentumEquation.C b/src/postProcessing/foamCalcFunctions/field/interpolate/writeInterpolateField.C similarity index 56% rename from applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveMomentumEquation.C rename to src/postProcessing/foamCalcFunctions/field/interpolate/writeInterpolateField.C index 935a02b60a..88b1a295aa 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveMomentumEquation.C +++ b/src/postProcessing/foamCalcFunctions/field/interpolate/writeInterpolateField.C @@ -22,36 +22,41 @@ 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 +template +void Foam::calcTypes::interpolate::writeInterpolateField ( - const bool momentumPredictor, - volVectorField& U, - const volScalarField& rho, - const surfaceScalarField& phi, - const volScalarField& pd, - const volScalarField& gh, - const compressible::turbulenceModel& turb + const IOobject& header, + const fvMesh& mesh, + bool& processed ) { - // Solve the Momentum equation - tmp UEqn - ( - fvm::ddt(rho, U) - + fvm::div(phi, U) - + turb.divDevRhoReff(U) - ); + typedef GeometricField fieldType; + typedef GeometricField surfaceFieldType; - UEqn().relax(); - - if (momentumPredictor) + if (header.headerClassName() == fieldType::typeName) { - solve(UEqn() == -fvc::grad(pd) - fvc::grad(rho)*gh); - } + Info<< " Reading " << header.name() << endl; + fieldType field(header, mesh); - return UEqn; + Info<< " Calculating interpolate" << header.name() << endl; + surfaceFieldType interpolateField + ( + IOobject + ( + "interpolate" + header.name(), + mesh.time().timeName(), + mesh, + IOobject::NO_READ + ), + fvc::interpolate(field) + ); + interpolateField.write(); + + processed = true; + } } + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/postProcessing/functionObjects/field/fieldAverage/fieldAverage/fieldAverage.C b/src/postProcessing/functionObjects/field/fieldAverage/fieldAverage/fieldAverage.C index 92b3f96366..b117112a53 100644 --- a/src/postProcessing/functionObjects/field/fieldAverage/fieldAverage/fieldAverage.C +++ b/src/postProcessing/functionObjects/field/fieldAverage/fieldAverage/fieldAverage.C @@ -82,7 +82,7 @@ void Foam::fieldAverage::initialise() // Add mean fields to the field lists forAll(faItems_, i) { - const word fieldName = faItems_[i].fieldName(); + const word& fieldName = faItems_[i].fieldName(); if (obr_.foundObject(fieldName)) { addMeanField(i, meanScalarFields_); @@ -117,7 +117,7 @@ void Foam::fieldAverage::initialise() { if (faItems_[i].prime2Mean()) { - const word fieldName = faItems_[i].fieldName(); + const word& fieldName = faItems_[i].fieldName(); if (!faItems_[i].mean()) { FatalErrorIn("Foam::fieldAverage::initialise()") diff --git a/src/thermophysicalModels/radiation/radiationModel/radiationModel/newRadiationModel.C b/src/thermophysicalModels/radiation/radiationModel/radiationModel/newRadiationModel.C index 0a213f03d3..9822578114 100644 --- a/src/thermophysicalModels/radiation/radiationModel/radiationModel/newRadiationModel.C +++ b/src/thermophysicalModels/radiation/radiationModel/radiationModel/newRadiationModel.C @@ -86,7 +86,7 @@ autoPtr radiationModel::New // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -} // End radiation +} // End namespace radiation } // End namespace Foam // ************************************************************************* // diff --git a/src/turbulenceModels/compressible/RAS/RASModel/RASModel.C b/src/turbulenceModels/compressible/RAS/RASModel/RASModel.C index 60d3dc3e84..f52646af28 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 ) ), diff --git a/wmake/wmake b/wmake/wmake index 61307c1180..473119fd01 100755 --- a/wmake/wmake +++ b/wmake/wmake @@ -163,8 +163,10 @@ then exit $? elif [ ! -d $MakeDir ] then - $make -k -f $WM_DIR/MakefileApps \ - FOAM_APPS="`find . -maxdepth 1 \( -type d -a ! -name "." -a ! -name "Optional" -a ! -name "Make" \) -printf "%f "`" + # FOAM_APPS=$(find . -maxdepth 1 \( -type d -a ! -name "." -a ! -name Optional -a ! -name Make \) -printf "%f ") + # avoid 'find' with '-printf' ... not entirely portable + FOAM_APPS=$(for d in *; do [ -d "$d" -a "$d" != Optional -a "$d" != Make ] && echo "$d"; done | xargs) + $make -k -f $WM_DIR/MakefileApps FOAM_APPS="$FOAM_APPS" exit $? fi