diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqFoam/Make/files b/applications/solvers/heatTransfer/buoyantBoussinesqFoam/Make/files
new file mode 100644
index 0000000000000000000000000000000000000000..d6a7819c0e016c560a996d803a3b0e160ebce735
--- /dev/null
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqFoam/Make/files
@@ -0,0 +1,3 @@
+buoyantBoussinesqFoam.C
+
+EXE = $(FOAM_APPBIN)/buoyantBoussinesqFoam
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqFoam/Make/options b/applications/solvers/heatTransfer/buoyantBoussinesqFoam/Make/options
new file mode 100644
index 0000000000000000000000000000000000000000..a53061a93cc503d9368d9d3a900ac6ef801015ba
--- /dev/null
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqFoam/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
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqFoam/TEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqFoam/TEqn.H
new file mode 100644
index 0000000000000000000000000000000000000000..e2a89d1441582292fcb5a33ec7e092eac455a5a5
--- /dev/null
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqFoam/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();
+}
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqFoam/UEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqFoam/UEqn.H
new file mode 100644
index 0000000000000000000000000000000000000000..3b5b5be24923c1b410d47c34fd14da319e117026
--- /dev/null
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqFoam/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()
+ )
+ );
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqFoam/buoyantBoussinesqFoam.C b/applications/solvers/heatTransfer/buoyantBoussinesqFoam/buoyantBoussinesqFoam.C
new file mode 100644
index 0000000000000000000000000000000000000000..72aab39ac044b8fe3ce59501c04204be4ce2808e
--- /dev/null
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqFoam/buoyantBoussinesqFoam.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
+ 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"
+# include "readTimeControls.H"
+# include "CourantNo.H"
+# include "setInitialDeltaT.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ Info<< "\nStarting time loop\n" << endl;
+
+ for (runTime++; !runTime.end(); runTime++)
+ {
+ 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;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqFoam/createFields.H b/applications/solvers/heatTransfer/buoyantBoussinesqFoam/createFields.H
new file mode 100644
index 0000000000000000000000000000000000000000..5f3f13626db0a860d7058e76ce5a90f397558c95
--- /dev/null
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqFoam/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
+ );
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqFoam/pdEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqFoam/pdEqn.H
new file mode 100644
index 0000000000000000000000000000000000000000..a74b4e80221b08df8743866215c7ef8c5311362c
--- /dev/null
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqFoam/pdEqn.H
@@ -0,0 +1,45 @@
+{
+ 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())
+ + fvc::ddtPhiCorr(rUA, U, phi)
+ + betaghf*fvc::snGrad(T)*rUAf*mesh.magSf();
+
+ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+ {
+ fvScalarMatrix pdEqn
+ (
+ fvm::laplacian(rUAf, pd) == fvc::div(phi)
+ );
+
+ // retain the residual from the first iteration
+ if (nonOrth == 0)
+ {
+ pdEqn.solve(mesh.solver(pd.name() + "Final"));
+ }
+ else
+ {
+ pdEqn.solve(mesh.solver(pd.name()));
+ }
+
+ if (nonOrth == nNonOrthCorr)
+ {
+ // Calculate the conservative fluxes
+ phi -= pdEqn.flux();
+
+ // Correct the momentum source with the pressure gradient flux
+ // calculated from the relaxed pressure
+ U -=
+ rUA
+ *fvc::reconstruct((pdEqn.flux() - betaghf*fvc::snGrad(T))/rUAf);
+ U.correctBoundaryConditions();
+ }
+ }
+
+ #include "continuityErrs.H"
+}
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqFoam/readTransportProperties.H b/applications/solvers/heatTransfer/buoyantBoussinesqFoam/readTransportProperties.H
new file mode 100644
index 0000000000000000000000000000000000000000..585128dfdeb7d5e5212f9412d9b415c05c15b752
--- /dev/null
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqFoam/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"));
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqFoam/writeAdditionalFields.H b/applications/solvers/heatTransfer/buoyantBoussinesqFoam/writeAdditionalFields.H
new file mode 100644
index 0000000000000000000000000000000000000000..20f7c6ae1dcff80630dc2ca21253b1e824ea5d77
--- /dev/null
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqFoam/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();
+}