diff --git a/ReleaseNotes-dev b/ReleaseNotes-dev new file mode 100644 index 0000000000..d35d1b518e --- /dev/null +++ b/ReleaseNotes-dev @@ -0,0 +1,130 @@ +# -*- mode: org; -*- +# +#+TITLE: OpenFOAM release notes for version dev +#+AUTHOR: OpenCFD Ltd. +#+DATE: TBA +#+LINK: http://www.openfoam.com +#+OPTIONS: author:nil ^:{} +# Copyright (c) 2010 OpenCFD Ltd. + +* Overview + OpenFOAM-dev is the latest major release of OpenFOAM including many new + developments a number of bug-fixes. This release passes our standard tests + and the tutorials have been broadly checked. Please report any bugs by + following the link: http://www.openfoam.com/bugs. + +* GNU/Linux version + This release of OpenFOAM is distributed primarily in 2 ways: (1) as a Debian + pack containing binaries and source; (2) from the SourceForge source code + repository (see [[./README.org][README]]). + + The Ubuntu/Debian pack is available for 32 and 64 bit versions of the 10.04 + LTS operating system using the system compiler and libraries that will be + installed automatically from standard Debian packs. + + To use the source version from the SourceForge repository, we provide a source + pack of third-party packages that can be compiled on the user's system. This + does not include =gcc=, since the system installed version is typically + sufficient, but includes =paraview-3.8.0=, =openmpi-1.4.1=, =scotch_5.1=, + =metis-5.0pre2=, =ParMetis-3.1= and =ParMGridGen-1.0=. + +* Library developments + There have been a number of developments to the libraries to support the + extension of functionality in solver and utility applications. +*** Core library + + Large number of code refinements and consistency improvements to support + other developments. +*** Turbulence modelling +*** Thermo-physical Models +*** DSMC +*** Dynamic Mesh +*** Numerics +*** *Updated* command line help, e.g. `snappyHexMesh -help' now gives: + + Usage: snappyHexMesh [OPTIONS] + options: + -case specify alternate case directory, default is the cwd + -overwrite overwrite existing mesh/results files + -parallel run in parallel + -srcDoc display source code in browser + -doc display application documentation in browser + -help print the usage + +*** *New* Surface film library + + Creation of films by particle addition, or initial film distribution + + Coupled with the lagrangian/intermediate cloud hierarchy library + + Hierarchical design, consisting of + + kinematic film: mass, momentum + + constant thermodynamic properties + + thermodynamic film: mass, momentum and enthalpy + + constant, or temperature dependant thermodynamic properties + + Sub-models: + + detachment/dripping whereby particles (re)enter the originating cloud + + particle sizes set according to PDF + + other properties set to ensure mass, momentum and energy conservation + + heat transfer to/from walls and film surface + + film evaporation and boiling + + Additional wall functions for primary region momentum and temperature + taking film into account + + Parallel aware +*** *Updated* particle tracking algorithm + +* Solvers + A number of new solvers have been developed for a range of engineering + applications. There has been a set of improvements to certain classes of + solver that are introduced in this release. +*** *New* Solvers + + =reactingParcelFilmFoam=: Lagrangian cloud and film transport in a + reacting gas phase system +*** Modifications to multiphase and buoyant solvers + + ... +*** Modifications to solvers for sensible enthalpy + + ... +*** Modifications to steady-state compressible solvers + + ... +*** Other modifications + + ... + +* Boundary conditions + New boundary conditions have been introduced to support new applications in + OpenFOAM. + + *New* wall functions: + + kappatJayatillekeWallFunction: incompressible RAS thermal wall function + +* Utilities + There have been some utilities added and updated in this release. +*** *New* utilities + + =extrudeToRegionMesh=: Extrude faceZones into separate mesh (as a + different region) + + used to e.g. extrude baffles (extrude internal faces) or create + liquid film regions + + if extruding internal faces: + + create baffles in original mesh with directMappedWall patches + + if extruding boundary faces: + + convert boundary faces to directMappedWall patches + + extrude edges of faceZone as a _sidePatch + + extrude edges inbetween different faceZones as a + (nonuniformTransform)cyclic _ + + extrudes into master direction (i.e. away from the owner cell + if flipMap is false) +*** Updated utilities + + ... + +* Post-processing + + =foamToEnsight=: new =-nodeValues= option to generate and output nodal + field data. + + Function objects: + + =residualControl=: new function object to allow users to terminate steady + state calculations when the defined residual levels are achieved + + =abortCalculation=: watches for presence of the named file in the + $FOAM_CASE directory and aborts the calculation if it is present + + =timeActivatedFileUpdate=: performs a file copy/replacement once a + specified time has been reached, e.g. to automagically change fvSchemes and + fvSolution during a calculation + + =streamLine=: generate streamlines; ouputs both trajectory and field data + +* New tutorials + There is a large number of new tutorials to support the new solvers in the + release. + + =reactingParcelFilmFoam= tutorials: + + multipleBoxes, hotBoxes, panel, evaporationTest diff --git a/applications/solvers/combustion/fireFoam/UEqn.H b/applications/solvers/combustion/fireFoam/UEqn.H index 4a2d03490d..a77c262831 100644 --- a/applications/solvers/combustion/fireFoam/UEqn.H +++ b/applications/solvers/combustion/fireFoam/UEqn.H @@ -13,7 +13,9 @@ solve == fvc::reconstruct ( - fvc::interpolate(rho)*(g & mesh.Sf()) - - fvc::snGrad(p)*mesh.magSf() + ( + - ghf*fvc::snGrad(rho) + - fvc::snGrad(p_rgh) + )*mesh.magSf() ) ); diff --git a/applications/solvers/combustion/fireFoam/combustionModels/Make/files b/applications/solvers/combustion/fireFoam/combustionModels/Make/files index 86696b0709..5b1262b70d 100644 --- a/applications/solvers/combustion/fireFoam/combustionModels/Make/files +++ b/applications/solvers/combustion/fireFoam/combustionModels/Make/files @@ -1,5 +1,5 @@ combustionModel/combustionModel.C -combustionModel/combustionModelNew.C +combustionModel/newCombustionModel.C infinitelyFastChemistry/infinitelyFastChemistry.C diff --git a/applications/solvers/combustion/fireFoam/combustionModels/combustionModel/combustionModelNew.C b/applications/solvers/combustion/fireFoam/combustionModels/combustionModel/newCombustionModel.C similarity index 77% rename from applications/solvers/combustion/fireFoam/combustionModels/combustionModel/combustionModelNew.C rename to applications/solvers/combustion/fireFoam/combustionModels/combustionModel/newCombustionModel.C index e23386589c..62781b0fdb 100644 --- a/applications/solvers/combustion/fireFoam/combustionModels/combustionModel/combustionModelNew.C +++ b/applications/solvers/combustion/fireFoam/combustionModels/combustionModel/newCombustionModel.C @@ -29,19 +29,22 @@ License Foam::autoPtr Foam::combustionModel::New ( - const dictionary& propDict, + const dictionary& combustionProperties, const hsCombustionThermo& thermo, const compressible::turbulenceModel& turbulence, const surfaceScalarField& phi, const volScalarField& rho ) { - const word modelType(propDict.lookup("combustionModel")); + word combustionModelTypeName = combustionProperties.lookup + ( + "combustionModel" + ); - Info<< "Selecting combustion model " << modelType << endl; + Info<< "Selecting combustion model " << combustionModelTypeName << endl; dictionaryConstructorTable::iterator cstrIter = - dictionaryConstructorTablePtr_->find(modelType); + dictionaryConstructorTablePtr_->find(combustionModelTypeName); if (cstrIter == dictionaryConstructorTablePtr_->end()) { @@ -49,14 +52,14 @@ Foam::autoPtr Foam::combustionModel::New ( "combustionModel::New" ) << "Unknown combustionModel type " - << modelType << nl << nl - << "Valid combustionModels are : " << endl - << dictionaryConstructorTablePtr_->toc() + << combustionModelTypeName << endl << endl + << "Valid combustionModels are : " << endl + << dictionaryConstructorTablePtr_->sortedToc() << exit(FatalError); } return autoPtr - (cstrIter()(propDict, thermo, turbulence, phi, rho)); + (cstrIter()(combustionProperties, thermo, turbulence, phi, rho)); } diff --git a/applications/solvers/combustion/fireFoam/combustionModels/infinitelyFastChemistry/infinitelyFastChemistry.H b/applications/solvers/combustion/fireFoam/combustionModels/infinitelyFastChemistry/infinitelyFastChemistry.H index 851f119feb..53832f21b9 100644 --- a/applications/solvers/combustion/fireFoam/combustionModels/infinitelyFastChemistry/infinitelyFastChemistry.H +++ b/applications/solvers/combustion/fireFoam/combustionModels/infinitelyFastChemistry/infinitelyFastChemistry.H @@ -88,8 +88,9 @@ public: ); - //- Destructor - virtual ~infinitelyFastChemistry(); + // Destructor + + virtual ~infinitelyFastChemistry(); // Member Functions diff --git a/applications/solvers/combustion/fireFoam/combustionModels/noCombustion/noCombustion.H b/applications/solvers/combustion/fireFoam/combustionModels/noCombustion/noCombustion.H index 9e15f0d300..b3f8bb225e 100644 --- a/applications/solvers/combustion/fireFoam/combustionModels/noCombustion/noCombustion.H +++ b/applications/solvers/combustion/fireFoam/combustionModels/noCombustion/noCombustion.H @@ -83,8 +83,9 @@ public: ); - //- Destructor - virtual ~noCombustion(); + // Destructor + + virtual ~noCombustion(); // Member Functions diff --git a/applications/solvers/combustion/fireFoam/createFields.H b/applications/solvers/combustion/fireFoam/createFields.H index 9a6274bb1a..6a1a40b777 100644 --- a/applications/solvers/combustion/fireFoam/createFields.H +++ b/applications/solvers/combustion/fireFoam/createFields.H @@ -35,6 +35,7 @@ const volScalarField& psi = thermo.psi(); volScalarField& ft = composition.Y("ft"); volScalarField& fu = composition.Y("fu"); + Info<< "Reading field U\n" << endl; volVectorField U @@ -73,7 +74,7 @@ IOdictionary combustionProperties Info<< "Creating combustion model\n" << endl; autoPtr combustion ( - combustionModel::New + combustionModel::combustionModel::New ( combustionProperties, thermo, @@ -83,6 +84,29 @@ autoPtr combustion ) ); + +Info<< "Calculating field g.h\n" << endl; +volScalarField gh("gh", g & mesh.C()); +surfaceScalarField ghf("gh", g & mesh.Cf()); + +Info<< "Reading field p_rgh\n" << endl; +volScalarField p_rgh +( + IOobject + ( + "p_rgh", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh +); + +// Force p_rgh to be consistent with p +p_rgh = p - rho*gh; + + volScalarField dQ ( IOobject @@ -103,15 +127,6 @@ volScalarField DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p); -Info<< "Calculating field g.h\n" << endl; -volScalarField gh("gh", g & mesh.C()); - -surfaceScalarField ghf("gh", g & mesh.Cf()); - -p += rho*gh; - -thermo.correct(); - dimensionedScalar initialMass = fvc::domainIntegrate(rho); diff --git a/applications/solvers/combustion/fireFoam/pEqn.H b/applications/solvers/combustion/fireFoam/pEqn.H index 8d1ed0d9f8..b31a7c6636 100644 --- a/applications/solvers/combustion/fireFoam/pEqn.H +++ b/applications/solvers/combustion/fireFoam/pEqn.H @@ -1,5 +1,3 @@ -bool closedVolume = false; - rho = thermo.rho(); volScalarField rUA = 1.0/UEqn.A(); @@ -15,43 +13,41 @@ surfaceScalarField phiU ) ); -phi = phiU + rhorUAf*fvc::interpolate(rho)*(g & mesh.Sf()); +phi = phiU - rhorUAf*ghf*fvc::snGrad(rho)*mesh.magSf(); for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) { - surfaceScalarField rhorUAf = fvc::interpolate(rho*rUA); + surfaceScalarField rhorUAf = fvc::interpolate(rho*rUA); - fvScalarMatrix pEqn + fvScalarMatrix p_rghEqn + ( + fvm::ddt(psi, p_rgh) + fvc::ddt(psi, rho)*gh + + fvc::div(phi) + - fvm::laplacian(rhorUAf, p_rgh) + ); + + p_rghEqn.solve + ( + mesh.solver ( - fvm::ddt(psi,p) - + fvc::div(phi) - - fvm::laplacian(rhorUAf, p) - ); - - closedVolume = p.needReference(); - - pEqn.solve - ( - mesh.solver + p_rgh.select ( - p.select ( - ( - finalIter - && corr == nCorr-1 - && nonOrth == nNonOrthCorr - ) + finalIter + && corr == nCorr-1 + && nonOrth == nNonOrthCorr ) ) - ); + ) + ); - if (nonOrth == nNonOrthCorr) - { - phi += pEqn.flux(); - } + if (nonOrth == nNonOrthCorr) + { + phi += p_rghEqn.flux(); + } } -DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p); +p = p_rgh + rho*gh; #include "rhoEqn.H" #include "compressibleContinuityErrs.H" @@ -59,12 +55,4 @@ DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p); U += rUA*fvc::reconstruct((phi - phiU)/rhorUAf); 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()); - rho = thermo.rho(); -} +DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p); diff --git a/applications/solvers/compressible/rhoPorousMRFPimpleFoam/createFields.H b/applications/solvers/compressible/rhoPorousMRFPimpleFoam/createFields.H index b9a86ef995..cc8f6436a1 100644 --- a/applications/solvers/compressible/rhoPorousMRFPimpleFoam/createFields.H +++ b/applications/solvers/compressible/rhoPorousMRFPimpleFoam/createFields.H @@ -39,9 +39,14 @@ #include "compressibleCreatePhi.H" - dimensionedScalar pMin + dimensionedScalar rhoMax ( - mesh.solutionDict().subDict("PIMPLE").lookup("pMin") + mesh.solutionDict().subDict("PIMPLE").lookup("rhoMax") + ); + + dimensionedScalar rhoMin + ( + mesh.solutionDict().subDict("PIMPLE").lookup("rhoMin") ); Info<< "Creating turbulence model\n" << endl; diff --git a/applications/solvers/compressible/rhoPorousMRFPimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPorousMRFPimpleFoam/pEqn.H index 56f8e9f3b5..6be6584202 100644 --- a/applications/solvers/compressible/rhoPorousMRFPimpleFoam/pEqn.H +++ b/applications/solvers/compressible/rhoPorousMRFPimpleFoam/pEqn.H @@ -102,6 +102,8 @@ else p.relax(); rho = thermo.rho(); + rho = max(rho, rhoMin); + rho = min(rho, rhoMax); rho.relax(); Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl; @@ -112,8 +114,6 @@ U.correctBoundaryConditions(); DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p); -bound(p, pMin); - // For closed-volume cases adjust the pressure and density levels // to obey overall mass continuity /* diff --git a/applications/solvers/compressible/rhoPorousSimpleFoam/createFields.H b/applications/solvers/compressible/rhoPorousSimpleFoam/createFields.H index d7a6b97ac1..1a242bdfff 100644 --- a/applications/solvers/compressible/rhoPorousSimpleFoam/createFields.H +++ b/applications/solvers/compressible/rhoPorousSimpleFoam/createFields.H @@ -70,16 +70,29 @@ dimensionedScalar initialMass = fvc::domainIntegrate(rho); thermalPorousZones pZones(mesh); + Switch pressureImplicitPorosity(false); // nUCorrectors used for pressureImplicitPorosity int nUCorr = 0; - const bool pressureImplicitPorosity = - ( - pZones.size() - && mesh.solutionDict().subDict("SIMPLE").readIfPresent - ( - "nUCorrectors", - nUCorr - ) - && (nUCorr > 0) - ); + if (pZones.size()) + { + // nUCorrectors for pressureImplicitPorosity + if (mesh.solutionDict().subDict("SIMPLE").found("nUCorrectors")) + { + nUCorr = readInt + ( + mesh.solutionDict().subDict("SIMPLE").lookup("nUCorrectors") + ); + } + + if (nUCorr > 0) + { + pressureImplicitPorosity = true; + Info<< "Using pressure implicit porosity" << endl; + } + else + { + Info<< "Using pressure explicit porosity" << endl; + } + } + diff --git a/applications/solvers/compressible/rhoPorousSimpleFoam/hEqn.H b/applications/solvers/compressible/rhoPorousSimpleFoam/hEqn.H index ff0b91bcb7..9e5138de8a 100644 --- a/applications/solvers/compressible/rhoPorousSimpleFoam/hEqn.H +++ b/applications/solvers/compressible/rhoPorousSimpleFoam/hEqn.H @@ -5,8 +5,8 @@ - fvm::Sp(fvc::div(phi), h) - fvm::laplacian(turbulence->alphaEff(), h) == - fvc::div(phi/fvc::interpolate(rho), rho/psi, "div(U,p)") - - (rho/psi)*fvc::div(phi/fvc::interpolate(rho)) + fvc::div(phi/fvc::interpolate(rho)*fvc::interpolate(p, "div(U,p)")) + - p*fvc::div(phi/fvc::interpolate(rho)) ); pZones.addEnthalpySource(thermo, rho, hEqn); diff --git a/applications/solvers/compressible/rhoPorousSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPorousSimpleFoam/pEqn.H index 0386899612..fe69384c8b 100644 --- a/applications/solvers/compressible/rhoPorousSimpleFoam/pEqn.H +++ b/applications/solvers/compressible/rhoPorousSimpleFoam/pEqn.H @@ -1,8 +1,3 @@ -rho = thermo.rho(); -rho = max(rho, rhoMin); -rho = min(rho, rhoMax); -rho.relax(); - if (pressureImplicitPorosity) { U = trTU()&UEqn().H(); diff --git a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H index f43ec6cf62..96d8280559 100644 --- a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H +++ b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H @@ -30,15 +30,7 @@ if (transonic) pEqn.setReference(pRefCell, pRefValue); - // Retain the residual from the first iteration - if (nonOrth == 0) - { - pEqn.solve(); - } - else - { - pEqn.solve(); - } + pEqn.solve(); if (nonOrth == nNonOrthCorr) { @@ -60,15 +52,7 @@ else pEqn.setReference(pRefCell, pRefValue); - // Retain the residual from the first iteration - if (nonOrth == 0) - { - pEqn.solve(); - } - else - { - pEqn.solve(); - } + pEqn.solve(); if (nonOrth == nNonOrthCorr) { diff --git a/applications/solvers/compressible/rhoSimplecFoam/createFields.H b/applications/solvers/compressible/rhoSimplecFoam/createFields.H index d97ee4705b..fab7b70048 100644 --- a/applications/solvers/compressible/rhoSimplecFoam/createFields.H +++ b/applications/solvers/compressible/rhoSimplecFoam/createFields.H @@ -43,9 +43,14 @@ scalar pRefValue = 0.0; setRefCell(p, mesh.solutionDict().subDict("SIMPLE"), pRefCell, pRefValue); - dimensionedScalar pMin + dimensionedScalar rhoMax ( - mesh.solutionDict().subDict("SIMPLE").lookup("pMin") + mesh.solutionDict().subDict("SIMPLE").lookup("rhoMax") + ); + + dimensionedScalar rhoMin + ( + mesh.solutionDict().subDict("SIMPLE").lookup("rhoMin") ); Info<< "Creating turbulence model\n" << endl; diff --git a/applications/solvers/compressible/rhoSimplecFoam/pEqn.H b/applications/solvers/compressible/rhoSimplecFoam/pEqn.H index a0f17e78bc..43443a507b 100644 --- a/applications/solvers/compressible/rhoSimplecFoam/pEqn.H +++ b/applications/solvers/compressible/rhoSimplecFoam/pEqn.H @@ -101,8 +101,6 @@ U -= (fvc::grad(p0)*(1.0/AU - 1.0/AtU) + fvc::grad(p)/AtU); U.correctBoundaryConditions(); -bound(p, pMin); - // For closed-volume cases adjust the pressure and density levels // to obey overall mass continuity if (closedVolume) @@ -112,6 +110,8 @@ if (closedVolume) } rho = thermo.rho(); +rho = max(rho, rhoMin); +rho = min(rho, rhoMax); if (!transonic) { diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/TEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/TEqn.H index dbfc61739f..9a835792a4 100644 --- a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/TEqn.H +++ b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/TEqn.H @@ -12,7 +12,7 @@ ); TEqn.relax(); - TEqn.solve(); + TEqn.solve(mesh.solver(T.select(finalIter))); rhok = 1.0 - beta*(T - TRef); } diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/UEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/UEqn.H index 35387f4ff4..20a05e5cd4 100644 --- a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/UEqn.H +++ b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/UEqn.H @@ -18,9 +18,10 @@ fvc::reconstruct ( ( - fvc::interpolate(rhok)*(g & mesh.Sf()) - - fvc::snGrad(p)*mesh.magSf() - ) - ) + - ghf*fvc::snGrad(rhok) + - fvc::snGrad(p_rgh) + )*mesh.magSf() + ), + mesh.solver(U.select(finalIter)) ); } diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C index ebf68d409c..54519906a4 100644 --- a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C +++ b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C @@ -87,7 +87,7 @@ int main(int argc, char *argv[]) if (nOuterCorr != 1) { - p.storePrevIter(); + p_rgh.storePrevIter(); } #include "UEqn.H" diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/createFields.H index 23d20cfa96..342af079d8 100644 --- a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/createFields.H +++ b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/createFields.H @@ -14,12 +14,12 @@ mesh ); - Info<< "Reading field p\n" << endl; - volScalarField p + Info<< "Reading field p_rgh\n" << endl; + volScalarField p_rgh ( IOobject ( - "p", + "p_rgh", runTime.timeName(), mesh, IOobject::MUST_READ, @@ -52,6 +52,18 @@ incompressible::RASModel::New(U, phi, laminarTransport) ); + // Kinematic density for buoyancy force + volScalarField rhok + ( + IOobject + ( + "rhok", + runTime.timeName(), + mesh + ), + 1.0 - beta*(T - TRef) + ); + // kinematic turbulent thermal thermal conductivity m2/s Info<< "Reading field kappat\n" << endl; volScalarField kappat @@ -67,25 +79,41 @@ mesh ); + Info<< "Calculating field g.h\n" << endl; + volScalarField gh("gh", g & mesh.C()); + surfaceScalarField ghf("ghf", g & mesh.Cf()); + + volScalarField p + ( + IOobject + ( + "p", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + p_rgh + rhok*gh + ); + label pRefCell = 0; scalar pRefValue = 0.0; setRefCell ( p, + p_rgh, mesh.solutionDict().subDict("PIMPLE"), pRefCell, pRefValue ); - - // Kinematic density for buoyancy force - volScalarField rhok - ( - IOobject + if (p_rgh.needReference()) + { + p += dimensionedScalar ( - "rhok", - runTime.timeName(), - mesh - ), - 1.0 - beta*(T - TRef) - ); + "p", + p.dimensions(), + pRefValue - getRefCellValue(p, pRefCell) + ); + } + diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H index 21be033f9b..60828e18dc 100644 --- a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H +++ b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H @@ -7,22 +7,23 @@ phi = (fvc::interpolate(U) & mesh.Sf()) + fvc::ddtPhiCorr(rUA, U, phi); - surfaceScalarField buoyancyPhi = - rUAf*fvc::interpolate(rhok)*(g & mesh.Sf()); - phi += buoyancyPhi; + surfaceScalarField buoyancyPhi = rUAf*ghf*fvc::snGrad(rhok)*mesh.magSf(); + phi -= buoyancyPhi; for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) { - fvScalarMatrix pEqn + fvScalarMatrix p_rghEqn ( - fvm::laplacian(rUAf, p) == fvc::div(phi) + fvm::laplacian(rUAf, p_rgh) == fvc::div(phi) ); - pEqn.solve + p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell)); + + p_rghEqn.solve ( mesh.solver ( - p.select + p_rgh.select ( ( finalIter @@ -36,17 +37,30 @@ if (nonOrth == nNonOrthCorr) { // Calculate the conservative fluxes - phi -= pEqn.flux(); + phi -= p_rghEqn.flux(); // Explicitly relax pressure for momentum corrector - p.relax(); + p_rgh.relax(); // Correct the momentum source with the pressure gradient flux // calculated from the relaxed pressure - U += rUA*fvc::reconstruct((buoyancyPhi - pEqn.flux())/rUAf); + U -= rUA*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rUAf); U.correctBoundaryConditions(); } } #include "continuityErrs.H" + + p = p_rgh + rhok*gh; + + if (p_rgh.needReference()) + { + p += dimensionedScalar + ( + "p", + p.dimensions(), + pRefValue - getRefCellValue(p, pRefCell) + ); + p_rgh = p - rhok*gh; + } } diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H index 76cc4da8ba..477a322833 100644 --- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H +++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H @@ -96,29 +96,23 @@ p_rgh + rhok*gh ); - label p_rghRefCell = 0; - scalar p_rghRefValue = 0.0; + label pRefCell = 0; + scalar pRefValue = 0.0; setRefCell ( + p, p_rgh, mesh.solutionDict().subDict("SIMPLE"), - p_rghRefCell, - p_rghRefValue + pRefCell, + pRefValue ); - scalar pRefValue = 0.0; - if (p_rgh.needReference()) { - pRefValue = readScalar - ( - mesh.solutionDict().subDict("SIMPLE").lookup("pRefValue") - ); - p += dimensionedScalar ( "p", p.dimensions(), - pRefValue - getRefCellValue(p, p_rghRefCell) + pRefValue - getRefCellValue(p, pRefCell) ); } diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H index 5664bb2154..50149ed360 100644 --- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H +++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H @@ -18,17 +18,9 @@ fvm::laplacian(rUAf, p_rgh) == fvc::div(phi) ); - p_rghEqn.setReference(p_rghRefCell, p_rghRefValue); + p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell)); - // retain the residual from the first iteration - if (nonOrth == 0) - { - p_rghEqn.solve(); - } - else - { - p_rghEqn.solve(); - } + p_rghEqn.solve(); if (nonOrth == nNonOrthCorr) { @@ -55,7 +47,8 @@ ( "p", p.dimensions(), - pRefValue - getRefCellValue(p, p_rghRefCell) + pRefValue - getRefCellValue(p, pRefCell) ); + p_rgh = p - rhok*gh; } } diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/UEqn.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/UEqn.H index d4878d063d..8c6a3f7671 100644 --- a/applications/solvers/heatTransfer/buoyantPimpleFoam/UEqn.H +++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/UEqn.H @@ -17,8 +17,11 @@ == fvc::reconstruct ( - fvc::interpolate(rho)*(g & mesh.Sf()) - - fvc::snGrad(p)*mesh.magSf() - ) + ( + - ghf*fvc::snGrad(rho) + - fvc::snGrad(p_rgh) + )*mesh.magSf() + ), + mesh.solver(U.select(finalIter)) ); } diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C b/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C index 0075c1e3ed..cb4c4d34f6 100644 --- a/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C +++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C @@ -80,7 +80,7 @@ int main(int argc, char *argv[]) if (nOuterCorr != 1) { - p.storePrevIter(); + p_rgh.storePrevIter(); } #include "UEqn.H" diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H index b8ac5595e4..e39d0bb803 100644 --- a/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H +++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H @@ -53,15 +53,30 @@ ) ); + Info<< "Calculating field g.h\n" << endl; + volScalarField gh("gh", g & mesh.C()); + surfaceScalarField ghf("ghf", g & mesh.Cf()); + + Info<< "Reading field p_rgh\n" << endl; + volScalarField p_rgh + ( + IOobject + ( + "p_rgh", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + // Force p_rgh to be consistent with p + p_rgh = p - rho*gh; + Info<< "Creating field DpDt\n" << endl; volScalarField DpDt ( "DpDt", fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p) ); - - thermo.correct(); - - dimensionedScalar initialMass = fvc::domainIntegrate(rho); - - dimensionedScalar totalVolume = sum(mesh.V()); diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/hEqn.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/hEqn.H index 3125cc3ffa..94537508b3 100644 --- a/applications/solvers/heatTransfer/buoyantPimpleFoam/hEqn.H +++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/hEqn.H @@ -9,7 +9,7 @@ ); hEqn.relax(); - hEqn.solve(); + hEqn.solve(mesh.solver(h.select(finalIter))); thermo.correct(); } diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H index a1c3dbfed5..e625f122a3 100644 --- a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H +++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H @@ -1,11 +1,9 @@ { - bool closedVolume = p.needReference(); - rho = thermo.rho(); // Thermodynamic density needs to be updated by psi*d(p) after the // pressure solution - done in 2 parts. Part 1: - thermo.rho() -= psi*p; + thermo.rho() -= psi*p_rgh; volScalarField rUA = 1.0/UEqn.A(); surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA)); @@ -18,24 +16,23 @@ + fvc::ddtPhiCorr(rUA, rho, U, phi) ); - surfaceScalarField buoyancyPhi = - rhorUAf*fvc::interpolate(rho)*(g & mesh.Sf()); + surfaceScalarField buoyancyPhi = -rhorUAf*ghf*fvc::snGrad(rho)*mesh.magSf(); phi += buoyancyPhi; for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) { - fvScalarMatrix pEqn + fvScalarMatrix p_rghEqn ( - fvc::ddt(rho) + psi*correction(fvm::ddt(p)) + fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh)) + fvc::div(phi) - - fvm::laplacian(rhorUAf, p) + - fvm::laplacian(rhorUAf, p_rgh) ); - pEqn.solve + p_rghEqn.solve ( mesh.solver ( - p.select + p_rgh.select ( ( finalIter @@ -49,34 +46,25 @@ if (nonOrth == nNonOrthCorr) { // Calculate the conservative fluxes - phi += pEqn.flux(); + phi += p_rghEqn.flux(); // Explicitly relax pressure for momentum corrector - p.relax(); + p_rgh.relax(); // Correct the momentum source with the pressure gradient flux // calculated from the relaxed pressure - U += rUA*fvc::reconstruct((buoyancyPhi + pEqn.flux())/rhorUAf); + U += rUA*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rhorUAf); U.correctBoundaryConditions(); } } + p = p_rgh + rho*gh; + // Second part of thermodynamic density update - thermo.rho() += psi*p; + thermo.rho() += psi*p_rgh; 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) - { - p += - (initialMass - fvc::domainIntegrate(psi*p)) - /fvc::domainIntegrate(psi); - thermo.rho() = psi*p; - rho += (initialMass - fvc::domainIntegrate(rho))/totalVolume; - } } diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C b/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C index 06fa5d12a0..64cc110cec 100644 --- a/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C +++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C @@ -62,10 +62,7 @@ int main(int argc, char *argv[]) { #include "UEqn.H" #include "hEqn.H" - for (int i=0; i<3; i++) - { #include "pEqn.H" - } } turbulence->correct(); diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H index cf23150332..d6fa9acee9 100644 --- a/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H +++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H @@ -23,20 +23,6 @@ volScalarField& h = thermo.h(); const volScalarField& psi = thermo.psi(); - Info<< "Reading field p_rgh\n" << endl; - volScalarField p_rgh - ( - IOobject - ( - "p_rgh", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh - ); - Info<< "Reading field U\n" << endl; volVectorField U ( @@ -53,7 +39,6 @@ #include "compressibleCreatePhi.H" - Info<< "Creating turbulence model\n" << endl; autoPtr turbulence ( @@ -66,40 +51,39 @@ ) ); + Info<< "Calculating field g.h\n" << endl; volScalarField gh("gh", g & mesh.C()); surfaceScalarField ghf("ghf", g & mesh.Cf()); - p = p_rgh + rho*gh; - thermo.correct(); - rho = thermo.rho(); - p_rgh = p - rho*gh; - - label p_rghRefCell = 0; - scalar p_rghRefValue = 0.0; - setRefCell + Info<< "Reading field p_rgh\n" << endl; + volScalarField p_rgh ( - p_rgh, - mesh.solutionDict().subDict("SIMPLE"), - p_rghRefCell, - p_rghRefValue + IOobject + ( + "p_rgh", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh ); + // Force p_rgh to be consistent with p + p_rgh = p - rho*gh; + + + label pRefCell = 0; scalar pRefValue = 0.0; - - if (p_rgh.needReference()) - { - pRefValue = readScalar - ( - mesh.solutionDict().subDict("SIMPLE").lookup("pRefValue") - ); - - p += dimensionedScalar - ( - "p", - p.dimensions(), - pRefValue - getRefCellValue(p, p_rghRefCell) - ); - } + setRefCell + ( + p, + p_rgh, + mesh.solutionDict().subDict("SIMPLE"), + pRefCell, + pRefValue + ); dimensionedScalar initialMass = fvc::domainIntegrate(rho); + dimensionedScalar totalVolume = sum(mesh.V()); diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H index 3088b42c69..f1a62d4770 100644 --- a/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H +++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H @@ -1,11 +1,12 @@ { rho = thermo.rho(); + rho.relax(); volScalarField rUA = 1.0/UEqn().A(); surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA)); U = rUA*UEqn().H(); - //UEqn.clear(); + UEqn.clear(); phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf()); bool closedVolume = adjustPhi(phi, U, p_rgh); @@ -20,7 +21,7 @@ fvm::laplacian(rhorUAf, p_rgh) == fvc::div(phi) ); - p_rghEqn.setReference(p_rghRefCell, p_rghRefValue); + p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell)); p_rghEqn.solve(); if (nonOrth == nNonOrthCorr) @@ -42,13 +43,13 @@ p = p_rgh + rho*gh; - // For closed-volume cases adjust the pressure and density levels + // For closed-volume cases adjust the pressure level // to obey overall mass continuity if (closedVolume) { p += (initialMass - fvc::domainIntegrate(psi*p)) /fvc::domainIntegrate(psi); - p_rgh == p - rho*gh; + p_rgh = p - rho*gh; } rho = thermo.rho(); diff --git a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C b/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C index 0dbe80c67c..6c35536333 100644 --- a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C +++ b/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C @@ -58,7 +58,7 @@ int main(int argc, char *argv[]) #include "readSIMPLEControls.H" - p.storePrevIter(); + p_rgh.storePrevIter(); rho.storePrevIter(); // Pressure-velocity SIMPLE corrector diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C index a81b0faaf3..36b3f1c3b0 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C @@ -93,6 +93,8 @@ int main(int argc, char *argv[]) // --- PIMPLE loop for (int oCorr=0; oCorr phiFluid(fluidRegions.size()); PtrList gFluid(fluidRegions.size()); PtrList turbulence(fluidRegions.size()); - PtrList DpDtf(fluidRegions.size()); + PtrList p_rghFluid(fluidRegions.size()); + PtrList ghFluid(fluidRegions.size()); + PtrList ghfFluid(fluidRegions.size()); List initialMassFluid(fluidRegions.size()); List