Merge branch 'master' into molecularDynamics

.gitignore

@@ -55,7 +55,7 @@ doc/[Dd]oxygen/man
 # ignore .timeStamp in the main directory
 /.timeStamp
 
-# ignore .ebrowse in the main directory
-/.ebrowse
+# ignore .tags in the main directory
+/.tags
 
 # end-of-file

applications/solvers/DNS/dnsFoam/dnsFoam.C

@@ -55,7 +55,7 @@ int main(int argc, char *argv[])
 
     Info<< nl << "Starting time loop" << endl;
 
-    for (runTime++; !runTime.end(); runTime++)
+    while (runTime.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
@@ -122,7 +122,7 @@ int main(int argc, char *argv[])
 
     Info<< "End\n" << endl;
 
-    return(0);
+    return 0;
 }
 
 

applications/solvers/Lagrangian/kinematicParcelFoam/kinematicParcelFoam.C

@@ -53,10 +53,8 @@ int main(int argc, char *argv[])
 
     Info<< "\nStarting time loop\n" << endl;
 
-    while (runTime.run())
+    while (runTime.loop())
     {
-        runTime++;
-
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
         Info<< "Evolving " << kinematicCloud.name() << endl;
@@ -72,7 +70,7 @@ int main(int argc, char *argv[])
 
     Info<< "End\n" << endl;
 
-    return(0);
+    return 0;
 }
 
 

applications/solvers/basic/laplacianFoam/laplacianFoam.C

@@ -48,7 +48,7 @@ int main(int argc, char *argv[])
 
     Info<< "\nCalculating temperature distribution\n" << endl;
 
-    for (runTime++; !runTime.end(); runTime++)
+    while (runTime.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
@@ -71,7 +71,7 @@ int main(int argc, char *argv[])
 
     Info<< "End\n" << endl;
 
-    return(0);
+    return 0;
 }
 
 

applications/solvers/basic/potentialFoam/potentialFoam.C

@@ -108,7 +108,7 @@ int main(int argc, char *argv[])
 
     Info<< "End\n" << endl;
 
-    return(0);
+    return 0;
 }
 
 

applications/solvers/basic/scalarTransportFoam/scalarTransportFoam.C

@@ -50,7 +50,7 @@ int main(int argc, char *argv[])
 
 #   include "CourantNo.H"
 
-    for (runTime++; !runTime.end(); runTime++)
+    while (runTime.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
@@ -71,7 +71,7 @@ int main(int argc, char *argv[])
 
     Info<< "End\n" << endl;
 
-    return(0);
+    return 0;
 }
 
 

applications/solvers/combustion/PDRFoam/PDRFoam.C

@@ -130,7 +130,7 @@ scalar StCoNum = 0.0;
 
     Info<< "\n end\n";
 
-    return(0);
+    return 0;
 }
 
 

applications/solvers/combustion/PDRFoam/PDRFoamAutoRefine.C

@@ -195,7 +195,7 @@ scalar StCoNum = 0.0;
 
     Info<< "\n end\n";
 
-    return(0);
+    return 0;
 }
 
 

applications/solvers/combustion/XiFoam/XiFoam.C

@@ -120,7 +120,7 @@ int main(int argc, char *argv[])
 
     Info<< "End\n" << endl;
 
-    return(0);
+    return 0;
 }
 
 

applications/solvers/combustion/coldEngineFoam/coldEngineFoam.C

@@ -94,7 +94,7 @@ int main(int argc, char *argv[])
 
     Info<< "End\n" << endl;
 
-    return(0);
+    return 0;
 }
 
 

applications/solvers/combustion/dieselEngineFoam/createSpray.H

@@ -32,7 +32,5 @@ if (dieselSpray.twoD())
     gasMass0 *= 2.0*mathematicalConstant::pi/dieselSpray.angleOfWedge();
 }
 
-reduce(gasMass0, sumOp<scalar>());
-
 gasMass0 -=
     dieselSpray.injectedMass(runTime.value()) - dieselSpray.liquidMass();

applications/solvers/combustion/dieselEngineFoam/dieselEngineFoam.C

@@ -132,7 +132,7 @@ int main(int argc, char *argv[])
 
     Info<< "End\n" << endl;
 
-    return(0);
+    return 0;
 }
 
 

applications/solvers/combustion/dieselEngineFoam/spraySummary.H

@@ -23,8 +23,6 @@
             gasMass *= 2.0*mathematicalConstant::pi/dieselSpray.angleOfWedge();
         }
 
-        reduce(gasMass, sumOp<scalar>());
-
         scalar addedMass = gasMass - gasMass0;
 
         Info<< "Added gas mass................. | " << 1e6*addedMass << " mg"

applications/solvers/combustion/dieselFoam/dieselFoam.C

@@ -63,7 +63,7 @@ int main(int argc, char *argv[])
 
     Info << "\nStarting time loop\n" << endl;
 
-    while(runTime.run())
+    while (runTime.run())
     {
 #       include "readPISOControls.H"
 #       include "compressibleCourantNo.H"
@@ -124,7 +124,7 @@ int main(int argc, char *argv[])
 
     Info<< "End\n" << endl;
 
-    return(0);
+    return 0;
 }
 
 

applications/solvers/combustion/engineFoam/engineFoam.C

@@ -128,7 +128,7 @@ int main(int argc, char *argv[])
 
     Info<< "End\n" << endl;
 
-    return(0);
+    return 0;
 }
 
 

applications/solvers/combustion/reactingFoam/reactingFoam.C

@@ -97,7 +97,7 @@ int main(int argc, char *argv[])
 
     Info<< "End\n" << endl;
 
-    return(0);
+    return 0;
 }
 
 

applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C

@@ -220,7 +220,7 @@ int main(int argc, char *argv[])
 
     Info<< "End\n" << endl;
 
-    return(0);
+    return 0;
 }
 
 // ************************************************************************* //

applications/solvers/compressible/rhoPisoFoam/rhoPisoFoam.C

@@ -88,7 +88,7 @@ int main(int argc, char *argv[])
 
     Info<< "End\n" << endl;
 
-    return(0);
+    return 0;
 }
 
 

applications/solvers/compressible/rhoPorousSimpleFoam/rhoPorousSimpleFoam.C

@@ -51,7 +51,7 @@ int main(int argc, char *argv[])
 
     Info<< "\nStarting time loop\n" << endl;
 
-    for (runTime++; !runTime.end(); runTime++)
+    while (runTime.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 

applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C

@@ -49,7 +49,7 @@ int main(int argc, char *argv[])
 
     Info<< "\nStarting time loop\n" << endl;
 
-    for (runTime++; !runTime.end(); runTime++)
+    while (runTime.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 

applications/solvers/compressible/rhoSonicFoam/rhoSonicFoam.C

@@ -47,7 +47,7 @@ int main(int argc, char *argv[])
 
     Info<< "\nStarting time loop\n" << endl;
 
-    for (runTime++; !runTime.end(); runTime++)
+    while (runTime.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
@@ -123,7 +123,7 @@ int main(int argc, char *argv[])
 
     Info<< "End\n" << endl;
 
-    return(0);
+    return 0;
 }
 
 

applications/solvers/compressible/rhopSonicFoam/rhopSonicFoam.C

@@ -53,7 +53,7 @@ int main(int argc, char *argv[])
 
     Info<< "\nStarting time loop\n" << endl;
 
-    for (runTime++; !runTime.end(); runTime++)
+    while (runTime.loop())
     {
         Info<< "Time = " << runTime.value() << nl << endl;
 
@@ -186,7 +186,7 @@ int main(int argc, char *argv[])
 
     Info<< "End\n" << endl;
 
-    return(0);
+    return 0;
 }
 
 

applications/solvers/compressible/sonicDyMFoam/sonicDyMFoam.C

@@ -52,7 +52,7 @@ int main(int argc, char *argv[])
 
     autoPtr<Foam::motionSolver> motionPtr = motionSolver::New(mesh);
 
-    for (runTime++; !runTime.end(); runTime++)
+    while (runTime.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
@@ -133,7 +133,7 @@ int main(int argc, char *argv[])
 
     Info<< "End\n" << endl;
 
-    return(0);
+    return 0;
 }
 
 

applications/solvers/compressible/sonicFoam/sonicFoam.C

@@ -49,7 +49,7 @@ int main(int argc, char *argv[])
 
     Info<< "\nStarting time loop\n" << endl;
 
-    for (runTime++; !runTime.end(); runTime++)
+    while (runTime.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
@@ -126,7 +126,7 @@ int main(int argc, char *argv[])
 
     Info<< "End\n" << endl;
 
-    return(0);
+    return 0;
 }
 
 

applications/solvers/compressible/sonicLiquidFoam/sonicLiquidFoam.C

@@ -50,7 +50,7 @@ int main(int argc, char *argv[])
 
     Info<< "\nStarting time loop\n" << endl;
 
-    for (runTime++; !runTime.end(); runTime++)
+    while (runTime.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
@@ -117,7 +117,7 @@ int main(int argc, char *argv[])
 
     Info<< "End\n" << endl;
 
-    return(0);
+    return 0;
 }
 
 

applications/solvers/discreteMethods/dsmc/dsmcFoam/Make/files

@@ -0,0 +1,3 @@
+dsmcFoam.C
+
+EXE = $(FOAM_APPBIN)/dsmcFoam

applications/solvers/discreteMethods/dsmc/dsmcFoam/Make/options

@@ -0,0 +1,12 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/lagrangian/dsmc/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
+
+EXE_LIBS = \
+    -lmeshTools \
+    -lfiniteVolume \
+    -llagrangian \
+    -ldsmc
+

applications/solvers/discreteMethods/dsmc/dsmcFoam/createFields.H

@@ -0,0 +1,162 @@
+
+    Info<< nl << "Reading field boundaryT" << endl;
+    volScalarField boundaryT
+    (
+        IOobject
+        (
+            "boundaryT",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< nl << "Reading field boundaryU" << endl;
+    volVectorField boundaryU
+    (
+        IOobject
+        (
+            "boundaryU",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< nl << "Reading field rhoN (number density)" << endl;
+    volScalarField rhoN
+    (
+        IOobject
+        (
+            "rhoN",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< nl << "Reading field rhoM (mass density)" << endl;
+    volScalarField rhoM
+    (
+        IOobject
+        (
+            "rhoM",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< nl << "Reading field rhoNdsmc (dsmc particle density)" << endl;
+    volScalarField dsmcRhoN
+    (
+        IOobject
+        (
+            "dsmcRhoN",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< nl << "Reading field momentum (momentum density)" << endl;
+    volVectorField momentum
+    (
+        IOobject
+        (
+            "momentum",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< nl << "Reading field linearKE (linear kinetic energy density)"
+        << endl;
+
+    volScalarField linearKE
+    (
+        IOobject
+        (
+            "linearKE",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< nl << "Reading field internalE (internal energy density)" << endl;
+    volScalarField internalE
+    (
+        IOobject
+        (
+            "internalE",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< nl << "Reading field iDof (internal degree of freedom density)"
+        << endl;
+
+    volScalarField iDof
+    (
+        IOobject
+        (
+            "iDof",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< nl << "Reading field q (surface heat transfer)" << endl;
+    volScalarField q
+    (
+        IOobject
+        (
+            "q",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< nl << "Reading field fD (surface force density)" << endl;
+    volVectorField fD
+    (
+        IOobject
+        (
+            "fD",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< nl << "Constructing dsmcCloud " << endl;
+
+    dsmcCloud dsmc("dsmc", boundaryT, boundaryU);

applications/solvers/discreteMethods/dsmc/dsmcFoam/dsmcFoam.C

@@ -0,0 +1,105 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+Application
+    dsmcFoam
+
+Description
+    Direct Simulation Monte Carlo Solver for 3D, transient, multi-species flows
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "dsmcCloud.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createFields.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        // Carry out dsmcCloud timestep
+
+        dsmc.evolve();
+
+        // Retrieve flow field data from dsmcCloud
+
+        rhoN = dsmc.rhoN();
+        rhoN.correctBoundaryConditions();
+
+        rhoM = dsmc.rhoM();
+        rhoM.correctBoundaryConditions();
+
+        dsmcRhoN = dsmc.dsmcRhoN();
+        dsmcRhoN.correctBoundaryConditions();
+
+        momentum = dsmc.momentum();
+        momentum.correctBoundaryConditions();
+
+        linearKE = dsmc.linearKE();
+        linearKE.correctBoundaryConditions();
+
+        internalE = dsmc.internalE();
+        internalE.correctBoundaryConditions();
+
+        iDof = dsmc.iDof();
+        iDof.correctBoundaryConditions();
+
+        // Retrieve surface field data from dsmcCloud
+
+        q = dsmc.q();
+
+        fD = dsmc.fD();
+
+        // Print status of dsmcCloud
+
+        dsmc.info();
+
+        runTime.write();
+
+        Info<< nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return(0);
+}
+
+
+// ************************************************************************* //

applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/mdEquilibrationFoam.C

@@ -52,9 +52,8 @@ int main(int argc, char *argv[])
 
     Info << "\nStarting time loop\n" << endl;
 
-    while (runTime.run())
+    while (runTime.loop())
     {
-        runTime++;
 
         nAveragingSteps++;
 
@@ -82,5 +81,5 @@ int main(int argc, char *argv[])
 
     Info << "End\n" << endl;
 
-    return(0);
+    return 0;
 }

applications/solvers/discreteMethods/molecularDynamics/mdFoam/mdFoam.C

@@ -50,9 +50,8 @@ int main(int argc, char *argv[])
 
     Info << "\nStarting time loop\n" << endl;
 
-    while (runTime.run())
+    while (runTime.loop())
     {
-        runTime++;
 
         nAveragingSteps++;
 
@@ -78,5 +77,5 @@ int main(int argc, char *argv[])
 
     Info << "End\n" << endl;
 
-    return(0);
+    return 0;
 }

applications/solvers/electromagnetics/electrostaticFoam/electrostaticFoam.C

@@ -47,7 +47,7 @@ int main(int argc, char *argv[])
 
     Info<< "\nStarting iteration loop\n" << endl;
 
-    for (runTime++; !runTime.end(); runTime++)
+    while (runTime.loop())
     {
         Info<< "Iteration = " << runTime.timeName() << nl << endl;
 
@@ -72,7 +72,7 @@ int main(int argc, char *argv[])
 
     Info<< "End\n" << endl;
 
-    return(0);
+    return 0;
 }
 
 

applications/solvers/electromagnetics/mhdFoam/mhdFoam.C

@@ -71,7 +71,7 @@ int main(int argc, char *argv[])
 
     Info<< nl << "Starting time loop" << endl;
 
-    for (runTime++; !runTime.end(); runTime++)
+    while (runTime.loop())
     {
 #       include "readPISOControls.H"
 #       include "readBPISOControls.H"
@@ -162,7 +162,7 @@ int main(int argc, char *argv[])
 
     Info<< "End\n" << endl;
 
-    return(0);
+    return 0;
 }
 
 

applications/solvers/financial/financialFoam/financialFoam.C

@@ -55,7 +55,7 @@ int main(int argc, char *argv[])
 
     Info<< "Starting time loop\n" << endl;
 
-    for (runTime++; !runTime.end(); runTime++)
+    while (runTime.loop())
     {
         delta == fvc::grad(V)().component(Foam::vector::X);
 
@@ -84,7 +84,7 @@ int main(int argc, char *argv[])
 
     Info<< "End\n" << endl;
 
-    return(0);
+    return 0;
 }
 
 

applications/solvers/heatTransfer/buoyantBoussinesqPisoFoam/Make/files

@@ -0,0 +1,3 @@
+buoyantBoussinesqPisoFoam.C
+
+EXE = $(FOAM_APPBIN)/buoyantBoussinesqPisoFoam

applications/solvers/heatTransfer/buoyantBoussinesqPisoFoam/Make/options

@@ -0,0 +1,12 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels \
+    -I$(LIB_SRC)/turbulenceModels/incompressible/RAS/lnInclude \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lincompressibleRASModels \
+    -lincompressibleTransportModels

applications/solvers/heatTransfer/buoyantBoussinesqPisoFoam/TEqn.H

@@ -0,0 +1,18 @@
+{
+    volScalarField kappaEff
+    (
+        "kappaEff",
+        turbulence->nu() + turbulence->nut()/Prt
+    );
+
+    fvScalarMatrix TEqn
+    (
+        fvm::ddt(T)
+      + fvm::div(phi, T)
+      - fvm::laplacian(kappaEff, T)
+    );
+
+    TEqn.relax();
+
+    TEqn.solve();
+}

applications/solvers/heatTransfer/buoyantBoussinesqPisoFoam/UEqn.H

@@ -0,0 +1,23 @@
+    // Solve the momentum equation
+
+    tmp<fvVectorMatrix> UEqn
+    (
+        fvm::ddt(U)
+      + fvm::div(phi, U)
+      + turbulence->divDevReff(U)
+    );
+
+    UEqn().relax();
+
+    solve
+    (
+        UEqn()
+      ==
+       -fvc::reconstruct
+        (
+            (
+                fvc::snGrad(pd)
+              - betaghf*fvc::snGrad(T)
+            ) * mesh.magSf()
+        )
+    );

applications/solvers/heatTransfer/buoyantBoussinesqPisoFoam/buoyantBoussinesqPisoFoam.C

@@ -0,0 +1,108 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+Application
+    buoyantBoussinesqPisoFoam
+
+Description
+    Transient solver for buoyant, turbulent flow of incompressible fluids
+
+    Uses the Boussinesq approximation:
+    \f[
+        rho_{eff} = 1 - beta(T - T_{ref})
+    \f]
+
+    where:
+        \f$ rho_{eff} \f$ = the effective (driving) density
+        beta = thermal expansion coefficient [1/K]
+        T = temperature [K]
+        \f$ T_{ref} \f$ = reference temperature [K]
+
+    Valid when:
+    \f[
+        rho_{eff} << 1
+    \f]
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "singlePhaseTransportModel.H"
+#include "RASModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+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 "CourantNo.H"
+#   include "setInitialDeltaT.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.loop())
+    {
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+#       include "readTimeControls.H"
+#       include "readPISOControls.H"
+#       include "CourantNo.H"
+#       include "setDeltaT.H"
+
+#       include "UEqn.H"
+
+        // --- PISO loop
+        for (int corr=0; corr<nCorr; corr++)
+        {
+#           include "TEqn.H"
+#           include "pdEqn.H"
+        }
+
+        turbulence->correct();
+
+        if (runTime.write())
+        {
+#           include "writeAdditionalFields.H"
+        }
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/heatTransfer/buoyantBoussinesqPisoFoam/createFields.H

@@ -0,0 +1,67 @@
+    Info<< "Reading thermophysical properties\n" << endl;
+
+    Info<< "Reading field T\n" << endl;
+    volScalarField T
+    (
+        IOobject
+        (
+            "T",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    // kinematic pd
+    Info<< "Reading field pd\n" << endl;
+    volScalarField pd
+    (
+        IOobject
+        (
+            "pd",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< "Reading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+#   include "createPhi.H"
+
+#   include "readTransportProperties.H"
+
+    Info<< "Creating turbulence model\n" << endl;
+    autoPtr<incompressible::RASModel> turbulence
+    (
+        incompressible::RASModel::New(U, phi, laminarTransport)
+    );
+
+    Info<< "Calculating field beta*(g.h)\n" << endl;
+    surfaceScalarField betaghf("betagh", beta*(g & mesh.Cf()));
+
+    label pdRefCell = 0;
+    scalar pdRefValue = 0.0;
+    setRefCell
+    (
+        pd,
+        mesh.solutionDict().subDict("SIMPLE"),
+        pdRefCell,
+        pdRefValue
+    );

applications/solvers/heatTransfer/buoyantBoussinesqPisoFoam/pdEqn.H

@@ -0,0 +1,42 @@
+{
+    volScalarField rUA("rUA", 1.0/UEqn().A());
+    surfaceScalarField rUAf("(1|A(U))", fvc::interpolate(rUA));
+
+    U = rUA*UEqn().H();
+    UEqn.clear();
+
+    surfaceScalarField phiU
+    (
+        (fvc::interpolate(U) & mesh.Sf())
+      + fvc::ddtPhiCorr(rUA, U, phi)
+    );
+
+    phi = phiU + betaghf*fvc::snGrad(T)*rUAf*mesh.magSf();
+
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    {
+        fvScalarMatrix pdEqn
+        (
+            fvm::laplacian(rUAf, pd) == fvc::div(phi)
+        );
+
+        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();
+        }
+    }
+
+    U -= rUA*fvc::reconstruct((phi - phiU)/rUAf);
+    U.correctBoundaryConditions();
+
+    #include "continuityErrs.H"
+}

applications/solvers/heatTransfer/buoyantBoussinesqPisoFoam/readTransportProperties.H

@@ -0,0 +1,13 @@
+    singlePhaseTransportModel laminarTransport(U, phi);
+
+    // thermal expansion coefficient [1/K]
+    dimensionedScalar beta(laminarTransport.lookup("beta"));
+
+    // reference temperature [K]
+    dimensionedScalar TRef(laminarTransport.lookup("TRef"));
+
+    // reference kinematic pressure [m2/s2]
+    dimensionedScalar pRef(laminarTransport.lookup("pRef"));
+
+    // turbulent Prandtl number
+    dimensionedScalar Prt(laminarTransport.lookup("Prt"));

applications/solvers/heatTransfer/buoyantBoussinesqPisoFoam/writeAdditionalFields.H

@@ -0,0 +1,29 @@
+{
+    volScalarField rhoEff
+    (
+        IOobject
+        (
+            "rhoEff",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        1.0 - beta*(T - TRef)
+    );
+    rhoEff.write();
+
+    volScalarField p
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        pd + rhoEff*(g & mesh.C()) + pRef
+    );
+    p.write();
+}

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/Make/files

@@ -0,0 +1,3 @@
+buoyantBoussinesqSimpleFoam.C
+
+EXE = $(FOAM_APPBIN)/buoyantBoussinesqSimpleFoam

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/Make/options

@@ -0,0 +1,12 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels \
+    -I$(LIB_SRC)/turbulenceModels/incompressible/RAS/lnInclude \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lincompressibleRASModels \
+    -lincompressibleTransportModels

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/TEqn.H

@@ -0,0 +1,19 @@
+{
+    volScalarField kappaEff
+    (
+        "kappaEff",
+        turbulence->nu() + turbulence->nut()/Prt
+    );
+
+    fvScalarMatrix TEqn
+    (
+        fvm::div(phi, T)
+      - fvm::Sp(fvc::div(phi), T)
+      - fvm::laplacian(kappaEff, T)
+    );
+
+    TEqn.relax();
+
+    eqnResidual = TEqn.solve().initialResidual();
+    maxResidual = max(eqnResidual, maxResidual);
+}

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/UEqn.H

@@ -0,0 +1,25 @@
+    // Solve the momentum equation
+
+    tmp<fvVectorMatrix> UEqn
+    (
+        fvm::div(phi, U)
+      - fvm::Sp(fvc::div(phi), U)
+      + turbulence->divDevReff(U)
+    );
+
+    UEqn().relax();
+
+    eqnResidual = solve
+    (
+        UEqn()
+      ==
+       -fvc::reconstruct
+        (
+            (
+                fvc::snGrad(pd)
+              - betaghf*fvc::snGrad(T)
+            ) * mesh.magSf()
+        )
+    ).initialResidual();
+
+    maxResidual = max(eqnResidual, maxResidual);

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/buoyantBoussinesqSimpleFoam.C

@@ -0,0 +1,105 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+Application
+    buoyantBoussinesqSimpleFoam
+
+Description
+    Steady-state solver for buoyant, turbulent flow of incompressible fluids
+
+    Uses the Boussinesq approximation:
+    \f[
+        rho_{eff} = 1 - beta(T - T_{ref})
+    \f]
+
+    where:
+        \f$ rho_{eff} \f$ = the effective (driving) density
+        beta = thermal expansion coefficient [1/K]
+        T = temperature [K]
+        \f$ T_{ref} \f$ = reference temperature [K]
+
+    Valid when:
+    \f[
+        rho_{eff} << 1
+    \f]
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "singlePhaseTransportModel.H"
+#include "RASModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+
+#   include "setRootCase.H"
+#   include "createTime.H"
+#   include "createMesh.H"
+#   include "readEnvironmentalProperties.H"
+#   include "createFields.H"
+#   include "initContinuityErrs.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.loop())
+    {
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+#       include "readSIMPLEControls.H"
+#       include "initConvergenceCheck.H"
+
+        pd.storePrevIter();
+
+        // Pressure-velocity SIMPLE corrector
+        {
+#           include "UEqn.H"
+#           include "TEqn.H"
+#           include "pdEqn.H"
+        }
+
+        turbulence->correct();
+
+        if (runTime.write())
+        {
+#           include "writeAdditionalFields.H"
+        }
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+
+#       include "convergenceCheck.H"
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H

@@ -0,0 +1,67 @@
+    Info<< "Reading thermophysical properties\n" << endl;
+
+    Info<< "Reading field T\n" << endl;
+    volScalarField T
+    (
+        IOobject
+        (
+            "T",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    // kinematic pd
+    Info<< "Reading field pd\n" << endl;
+    volScalarField pd
+    (
+        IOobject
+        (
+            "pd",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< "Reading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+#   include "createPhi.H"
+
+#   include "readTransportProperties.H"
+
+    Info<< "Creating turbulence model\n" << endl;
+    autoPtr<incompressible::RASModel> turbulence
+    (
+        incompressible::RASModel::New(U, phi, laminarTransport)
+    );
+
+    Info<< "Calculating field beta*(g.h)\n" << endl;
+    surfaceScalarField betaghf("betagh", beta*(g & mesh.Cf()));
+
+    label pdRefCell = 0;
+    scalar pdRefValue = 0.0;
+    setRefCell
+    (
+        pd,
+        mesh.solutionDict().subDict("SIMPLE"),
+        pdRefCell,
+        pdRefValue
+    );

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pdEqn.H

@@ -0,0 +1,49 @@
+{
+    volScalarField rUA("rUA", 1.0/UEqn().A());
+    surfaceScalarField rUAf("(1|A(U))", fvc::interpolate(rUA));
+
+    U = rUA*UEqn().H();
+    UEqn.clear();
+
+    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++)
+    {
+        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();
+        }
+    }
+
+    #include "continuityErrs.H"
+}

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/readTransportProperties.H

@@ -0,0 +1,13 @@
+    singlePhaseTransportModel laminarTransport(U, phi);
+
+    // thermal expansion coefficient [1/K]
+    dimensionedScalar beta(laminarTransport.lookup("beta"));
+
+    // reference temperature [K]
+    dimensionedScalar TRef(laminarTransport.lookup("TRef"));
+
+    // reference kinematic pressure [m2/s2]
+    dimensionedScalar pRef(laminarTransport.lookup("pRef"));
+
+    // turbulent Prandtl number
+    dimensionedScalar Prt(laminarTransport.lookup("Prt"));

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/writeAdditionalFields.H

@@ -0,0 +1,29 @@
+{
+    volScalarField rhoEff
+    (
+        IOobject
+        (
+            "rhoEff",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        1.0 - beta*(T - TRef)
+    );
+    rhoEff.write();
+
+    volScalarField p
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        pd + rhoEff*(g & mesh.C()) + pRef
+    );
+    p.write();
+}

applications/solvers/heatTransfer/buoyantFoam/Make/files

@@ -1,3 +0,0 @@
-buoyantFoam.C
-
-EXE = $(FOAM_APPBIN)/buoyantFoam

applications/solvers/heatTransfer/buoyantFoam/UEqn.H

@@ -1,12 +0,0 @@
-    // Solve the Momentum equation
-
-    fvVectorMatrix UEqn
-    (
-        fvm::ddt(rho, U)
-      + fvm::div(phi, U)
-      + turbulence->divDevRhoReff(U)
-    );
-
-    UEqn.relax();
-
-    solve(UEqn == -fvc::grad(pd) - fvc::grad(rho)*gh);

applications/solvers/heatTransfer/buoyantFoam/pEqn.H

@@ -1,60 +0,0 @@
-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
-    (
-        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();
-    }
-}
-
-p == pd + rho*gh + pRef;
-dpdt = fvc::ddt(p);
-
-#include "rhoEqn.H"
-#include "compressibleContinuityErrs.H"
-
-U -= rUA*(fvc::grad(pd) + fvc::grad(rho)*gh);
-U.correctBoundaryConditions();
-
-
-// For closed-volume cases adjust the pressure and density levels
-// to obey overall mass continuity
-if (closedVolume)
-{
-    p += (initialMass - fvc::domainIntegrate(thermo->psi()*p))
-        /fvc::domainIntegrate(thermo->psi());
-    pd == p - (rho*gh + pRef);
-    rho = thermo->rho();
-}

applications/solvers/heatTransfer/buoyantPisoFoam/Make/files

@@ -0,0 +1,3 @@
+buoyantPisoFoam.C
+
+EXE = $(FOAM_APPBIN)/buoyantPisoFoam

applications/solvers/heatTransfer/buoyantPisoFoam/Make/options

@@ -1,7 +1,6 @@
 EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
-    -I../XiFoam \
     -I$(LIB_SRC)/finiteVolume/lnInclude
 
 EXE_LIBS = \

applications/solvers/heatTransfer/buoyantPisoFoam/UEqn.H

@@ -0,0 +1,26 @@
+    // Solve the Momentum equation
+
+    fvVectorMatrix UEqn
+    (
+        fvm::ddt(rho, U)
+      + fvm::div(phi, U)
+      + turbulence->divDevRhoReff(U)
+    );
+
+    UEqn.relax();
+
+    if (momentumPredictor)
+    {
+        solve
+        (
+            UEqn
+         ==
+           -fvc::reconstruct
+            (
+                (
+                    fvc::snGrad(pd)
+                  + ghf*fvc::snGrad(rho)
+                ) * mesh.magSf()
+            )
+        );
+    }

applications/solvers/heatTransfer/buoyantPisoFoam/buoyantPisoFoam.C

@@ -23,7 +23,7 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Application
-    buoyantFoam
+    buoyantPisoFoam
 
 Description
     Transient Solver for buoyant, turbulent flow of compressible fluids for
@@ -41,37 +41,42 @@ 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"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
     while (runTime.run())
     {
-#       include "readPISOControls.H"
-#       include "compressibleCourantNo.H"
-//#       include "setDeltaT.H"
+        #include "readTimeControls.H"
+        #include "readPISOControls.H"
+        #include "compressibleCourantNo.H"
+        #include "setDeltaT.H"
 
         runTime++;
+
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-#       include "rhoEqn.H"
+        #include "rhoEqn.H"
 
-#       include "UEqn.H"
+        #include "UEqn.H"
+
+        #include "hEqn.H"
 
         // --- PISO loop
 
         for (int corr=0; corr<nCorr; corr++)
         {
-#           include "hEqn.H"
-#           include "pEqn.H"
+            #include "pEqn.H"
         }
 
         turbulence->correct();

applications/solvers/heatTransfer/buoyantPisoFoam/createFields.H

@@ -52,11 +52,16 @@
         )
     );
 
-    Info<< "Creating field dpdt\n" << endl;
-    volScalarField dpdt = fvc::ddt(p);
+    Info<< "Creating field DpDt\n" << endl;
+    volScalarField DpDt
+    (
+        "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("ghf", g & mesh.Cf());
 
     dimensionedScalar pRef("pRef", p.dimensions(), thermo->lookup("pRef"));
 

applications/solvers/heatTransfer/buoyantPisoFoam/hEqn.H

@@ -5,9 +5,7 @@
       + fvm::div(phi, h)
       - fvm::laplacian(turbulence->alphaEff(), h)
      ==
-        dpdt
-      + fvc::div(phi/fvc::interpolate(rho)*fvc::interpolate(p))
-      - p*fvc::div(phi/fvc::interpolate(rho))
+        DpDt
     );
 
     hEqn.relax();

applications/solvers/heatTransfer/buoyantPisoFoam/pEqn.H

@@ -0,0 +1,67 @@
+{
+    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
+    (
+        fvc::interpolate(rho)
+       *(
+            (fvc::interpolate(U) & mesh.Sf())
+          + fvc::ddtPhiCorr(rUA, rho, U, phi)
+        )
+    );
+
+    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(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();
+        }
+    }
+
+    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)
+    {
+        p += (initialMass - fvc::domainIntegrate(thermo->psi()*p))
+            /fvc::domainIntegrate(thermo->psi());
+        rho = thermo->rho();
+    }
+
+    pd == p - (rho*gh + pRef);
+}

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);

applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C

@@ -51,7 +51,7 @@ int main(int argc, char *argv[])
 
     Info<< "\nStarting time loop\n" << endl;
 
-    for (runTime++; !runTime.end(); runTime++)
+    while (runTime.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 

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("ghf", g & mesh.Cf());
 
     dimensionedScalar pRef("pRef", p.dimensions(), thermo->lookup("pRef"));
 

applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H

@@ -1,55 +1,66 @@
-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;

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 \

applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/UEqn.H

@@ -1,18 +0,0 @@
-    // Solve the Momentum equation
-
-    tmp<fvVectorMatrix> 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);
-

applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C

@@ -47,13 +47,14 @@ int main(int argc, char *argv[])
 #   include "createMesh.H"
 #   include "readEnvironmentalProperties.H"
 #   include "createFields.H"
+#   include "createRadiationModel.H"
 #   include "initContinuityErrs.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
-    for (runTime++; !runTime.end(); runTime++)
+    while (runTime.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 

applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/createFields.H

@@ -1,96 +0,0 @@
-    Info<< "Reading thermophysical properties\n" << endl;
-
-    autoPtr<basicThermo> thermo
-    (
-        basicThermo::New(mesh)
-    );
-
-    volScalarField rho
-    (
-        IOobject
-        (
-            "rho",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        thermo->rho()
-    );
-
-    volScalarField& p = thermo->p();
-    volScalarField& h = thermo->h();
-    const volScalarField& T = thermo->T();
-
-
-    Info<< "Reading field U\n" << endl;
-    volVectorField U
-    (
-        IOobject
-        (
-            "U",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-#   include "compressibleCreatePhi.H"
-
-
-    Info<< "Creating turbulence model\n" << endl;
-    autoPtr<compressible::RASModel> turbulence
-    (
-        compressible::RASModel::New
-        (
-            rho,
-            U,
-            phi,
-            thermo()
-        )
-    );
-
-    Info<< "Calculating field g.h\n" << endl;
-    volScalarField gh("gh", g & mesh.C());
-
-    dimensionedScalar pRef("pRef", p.dimensions(), thermo->lookup("pRef"));
-
-    Info<< "Creating field pd\n" << endl;
-    volScalarField pd
-    (
-        IOobject
-        (
-            "pd",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    p = pd + rho*gh + pRef;
-    thermo->correct();
-
-
-    label pdRefCell = 0;
-    scalar pdRefValue = 0.0;
-    setRefCell
-    (
-        pd,
-        mesh.solutionDict().subDict("SIMPLE"),
-        pdRefCell,
-        pdRefValue
-    );
-
-
-    Info<< "Creating radiation model\n" << endl;
-    autoPtr<radiation::radiationModel> radiation
-    (
-        radiation::radiationModel::New(T)
-    );
-
-
-    dimensionedScalar initialMass = fvc::domainIntegrate(rho);

applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/createRadiationModel.H

@@ -0,0 +1,5 @@
+    Info<< "Creating radiation model\n" << endl;
+    autoPtr<radiation::radiationModel> radiation
+    (
+        radiation::radiationModel::New(thermo->T())
+    );

applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/pEqn.H

@@ -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;

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

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"
@@ -71,7 +64,7 @@ int main(int argc, char *argv[])
 #       include "setInitialDeltaT.H"
     }
 
-    while(runTime.run())
+    while (runTime.run())
     {
 #       include "readTimeControls.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"
         }
@@ -110,7 +105,7 @@ int main(int argc, char *argv[])
 
     Info << "End\n" << endl;
 
-    return(0);
+    return 0;
 }
 
 

applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.H

@@ -103,8 +103,8 @@ public:
             const dictionary&
         );
 
-        //- Construct by mapping given solidWallMixedTemperatureCoupledFvPatchScalarField
-        //  onto a new patch
+        //- Construct by mapping given
+        //  solidWallMixedTemperatureCoupledFvPatchScalarField onto a new patch
         solidWallMixedTemperatureCoupledFvPatchScalarField
         (
             const solidWallMixedTemperatureCoupledFvPatchScalarField&,

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/UEqn.H

@@ -1,10 +1,25 @@
-    tmp<fvVectorMatrix> UEqn = solveMomentumEquation
+    // Solve the Momentum equation
+    tmp<fvVectorMatrix> 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()
+            )
+        );
+    }

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleContinuityErrors.C

@@ -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;
-}

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;
+}

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,
@@ -39,7 +38,8 @@ scalar compressibleCourantNo
     scalar CoNum = 0.0;
     scalar meanCoNum = 0.0;
 
-    if (mesh.nInternalFaces())
+    //- Can have fluid domains with 0 cells so do not test.
+    //if (mesh.nInternalFaces())
     {
         surfaceScalarField SfUfbyDelta =
             mesh.surfaceInterpolation::deltaCoeffs()

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleCourantNo.H

@@ -0,0 +1,49 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+    Calculates and outputs the mean and maximum Courant Numbers for the fluid
+    regions
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef compressibleCourantNo_H
+#define compressibleCourantNo_H
+
+#include "fvMesh.H"
+
+namespace Foam
+{
+    scalar compressibleCourantNo
+    (
+        const fvMesh& mesh,
+        const Time& runTime,
+        const volScalarField& rho,
+        const surfaceScalarField& phi
+    );
+}
+
+#endif
+
+// ************************************************************************* //

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleMultiRegionCourantNo.H

@@ -7,8 +7,8 @@
             (
                 fluidRegions[regionI],
                 runTime,
-                rhof[regionI],
-                phif[regionI]
+                rhoFluid[regionI],
+                phiFluid[regionI]
             ),
             CoNum
         );

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H

@@ -1,15 +1,16 @@
     // Initialise fluid field pointer lists
-    PtrList<basicThermo> thermof(fluidRegions.size());
-    PtrList<volScalarField> rhof(fluidRegions.size());
-    PtrList<volScalarField> Kf(fluidRegions.size());
-    PtrList<volVectorField> Uf(fluidRegions.size());
-    PtrList<surfaceScalarField> phif(fluidRegions.size());
-    PtrList<compressible::turbulenceModel> turb(fluidRegions.size());
-    PtrList<volScalarField> DpDtf(fluidRegions.size());
-    PtrList<volScalarField> ghf(fluidRegions.size());
-    PtrList<volScalarField> pdf(fluidRegions.size());
+    PtrList<basicThermo> thermoFluid(fluidRegions.size());
+    PtrList<volScalarField> rhoFluid(fluidRegions.size());
+    PtrList<volScalarField> KFluid(fluidRegions.size());
+    PtrList<volVectorField> UFluid(fluidRegions.size());
+    PtrList<surfaceScalarField> phiFluid(fluidRegions.size());
+    PtrList<compressible::turbulenceModel> turbulence(fluidRegions.size());
+    PtrList<volScalarField> DpDtFluid(fluidRegions.size());
+    PtrList<volScalarField> ghFluid(fluidRegions.size());
+    PtrList<surfaceScalarField> ghfFluid(fluidRegions.size());
+    PtrList<volScalarField> pdFluid(fluidRegions.size());
 
-    List<scalar> initialMassf(fluidRegions.size());
+    List<scalar> 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>
             (
                 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();
     }
 
 

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/hEqn.H

@@ -1,9 +1,17 @@
-    solveEnthalpyEquation
+{
+    tmp<fvScalarMatrix> 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()).value() << ' '
+        << max(thermo.T()).value() << endl;
+}

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/initContinuityErrs.H

@@ -0,0 +1 @@
+List<scalar> cumulativeContErr(fluidRegions.size(), 0.0);

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);
 }

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pdEqn.H

@@ -1,14 +0,0 @@
-    solvePressureDifferenceEquation
-    (
-        corr,
-        nCorr,
-        nNonOrthCorr,
-        closedVolume,
-        pdf[i],
-        pRef,
-        rhof[i],
-        thermof[i].psi(),
-        rUA,
-        ghf[i],
-        phif[i]
-    );

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/rhoEqn.H

@@ -1 +1 @@
-    solveContinuityEquation(rhof[i], phif[i]);
+    solve(fvm::ddt(rho) + fvc::div(phi));

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setInitialDeltaT.H

@@ -1,15 +0,0 @@
-    if (adjustTimeStep)
-    {
-        if (CoNum > SMALL)
-        {
-            runTime.setDeltaT
-            (
-                min
-                (
-                    maxCo*runTime.deltaT().value()/CoNum,
-                    maxDeltaT
-                )
-            );
-        }
-    }
-

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]);

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveEnthalpyEquation.C

@@ -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<fvScalarMatrix> 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;
-}

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveFluid.H

@@ -12,4 +12,4 @@
     #       include "pEqn.H"
         }
     }
-    turb[i].correct();
+    turb.correct();

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveMomentumEquation.C

@@ -1,57 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software; you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation; either version 2 of the License, or (at your
-    option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM; if not, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Description
-    Solve momentum equation and return matrix for use in pressure equation
-
-\*---------------------------------------------------------------------------*/
-
-tmp<fvVectorMatrix> solveMomentumEquation
-(
-    const bool momentumPredictor,
-    volVectorField& U,
-    const volScalarField& rho,
-    const surfaceScalarField& phi,
-    const volScalarField& pd,
-    const volScalarField& gh,
-    const compressible::turbulenceModel& turb
-)
-{
-    // Solve the Momentum equation
-    tmp<fvVectorMatrix> UEqn
-    (
-        fvm::ddt(rho, U)
-      + fvm::div(phi, U)
-      + turb.divDevRhoReff(U)
-    );
-
-    UEqn().relax();
-
-    if (momentumPredictor)
-    {
-        solve(UEqn() == -fvc::grad(pd) - fvc::grad(rho)*gh);
-    }
-
-    return UEqn;
-}

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solvePressureDifferenceEquation.C

@@ -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();
-        }
-    }
-}