Merge branch 'master' of /home/noisy3/OpenFOAM/OpenFOAM-dev

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/TEqn.H

@@ -2,7 +2,7 @@
     volScalarField kappaEff
     (
         "kappaEff",
-        turbulence->nu() + turbulence->nut()/Pr
+        turbulence->nu() + turbulence->nut()/Prt
     );
 
     fvScalarMatrix TEqn

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/UEqn.H

@@ -1,4 +1,4 @@
-    // Solve the Momentum equation
+    // Solve the momentum equation
 
     tmp<fvVectorMatrix> UEqn
     (
@@ -13,8 +13,13 @@
     (
         UEqn()
       ==
-       -fvc::grad(pd)
+       -fvc::reconstruct
-      + beta*gh*fvc::grad(T)
+        (
+            (
+                fvc::snGrad(pd)
+              - betaghf*fvc::snGrad(T)
+            ) * mesh.magSf()
+        )
     ).initialResidual();
 
     maxResidual = max(eqnResidual, maxResidual);

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H

@@ -53,8 +53,8 @@
         incompressible::RASModel::New(U, phi, laminarTransport)
     );
 
-    Info<< "Calculating field g.h\n" << endl;
+    Info<< "Calculating field beta*(g.h)\n" << endl;
-    volScalarField gh("gh", g & mesh.C());
+    surfaceScalarField betaghf("betagh", beta*(g & mesh.Cf()));
 
     label pdRefCell = 0;
     scalar pdRefValue = 0.0;

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pdEqn.H

@@ -1,16 +1,20 @@
-volScalarField rUA = 1.0/UEqn().A();
+{
+    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);
-phi += fvc::interpolate(beta*gh*rUA)*fvc::snGrad(T)*mesh.magSf();
+    surfaceScalarField buoyancyPhi = -betaghf*fvc::snGrad(T)*rUAf*mesh.magSf();
+    phi -= buoyancyPhi;
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
         fvScalarMatrix pdEqn
         (
-        fvm::laplacian(rUA, pd) == fvc::div(phi)
+            fvm::laplacian(rUAf, pd) == fvc::div(phi)
         );
 
         pdEqn.setReference(pdRefCell, pdRefValue);
@@ -28,14 +32,18 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
 
         if (nonOrth == nNonOrthCorr)
         {
+            // Calculate the conservative fluxes
             phi -= pdEqn.flux();
-    }
-}
-
-#include "continuityErrs.H"
 
             // Explicitly relax pressure for momentum corrector
             pd.relax();
 
-U -= rUA*(fvc::grad(pd) - beta*gh*fvc::grad(T));
+            // 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

@@ -9,5 +9,5 @@
     // reference kinematic pressure [m2/s2]
     dimensionedScalar pRef(laminarTransport.lookup("pRef"));
 
-    // Prandtl number
+    // turbulent Prandtl number
-    dimensionedScalar Pr(laminarTransport.lookup("Pr"));
+    dimensionedScalar Prt(laminarTransport.lookup("Prt"));

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/writeAdditionalFields.H

@@ -23,7 +23,7 @@
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        pd + rhoEff*gh + pRef
+        pd + rhoEff*(g & mesh.C()) + pRef
     );
     p.write();
 }

applications/solvers/heatTransfer/buoyantFoam/UEqn.H

@@ -9,4 +9,18 @@
 
     UEqn.relax();
 
-    solve(UEqn == -fvc::grad(pd) - fvc::grad(rho)*gh);
+    if (momentumPredictor)
+    {
+        solve
+        (
+            UEqn
+         ==
+           -fvc::reconstruct
+            (
+                (
+                    fvc::snGrad(pd)
+                  + ghf*fvc::snGrad(rho)
+                ) * mesh.magSf()
+            )
+        );
+    }

applications/solvers/heatTransfer/buoyantFoam/buoyantFoam.C

@@ -41,7 +41,6 @@ Description
 
 int main(int argc, char *argv[])
 {
-
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
@@ -52,7 +51,7 @@ int main(int argc, char *argv[])
     #include "compressibleCourantNo.H"
     #include "setInitialDeltaT.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 

applications/solvers/heatTransfer/buoyantFoam/createFields.H

@@ -58,6 +58,7 @@
 
     Info<< "Calculating field g.h\n" << endl;
     volScalarField gh("gh", g & mesh.C());
+    surfaceScalarField ghf("gh", g & mesh.Cf());
 
     dimensionedScalar pRef("pRef", p.dimensions(), thermo->lookup("pRef"));
 

applications/solvers/heatTransfer/buoyantFoam/pEqn.H

@@ -1,17 +1,23 @@
+{
     bool closedVolume = pd.needReference();
 
     rho = thermo->rho();
 
     volScalarField rUA = 1.0/UEqn.A();
+    surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA));
+
     U = rUA*UEqn.H();
 
-phi =
+    surfaceScalarField phiU
+    (
         fvc::interpolate(rho)
        *(
            (fvc::interpolate(U) & mesh.Sf())
          + fvc::ddtPhiCorr(rUA, rho, U, phi)
         )
-  - fvc::interpolate(rho*rUA*gh)*fvc::snGrad(rho)*mesh.magSf();
+    );
+
+    phi = phiU - ghf*fvc::snGrad(rho)*rhorUAf*mesh.magSf();
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
@@ -21,7 +27,7 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
           + fvc::ddt(psi)*pRef
           + fvc::ddt(psi, rho)*gh
           + fvc::div(phi)
-      - fvm::laplacian(rho*rUA, pd)
+          - fvm::laplacian(rhorUAf, pd)
         );
 
         if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
@@ -39,22 +45,23 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
         }
     }
 
+    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"
 
-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();
     }
+
+    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/createFields.H

@@ -53,6 +53,7 @@
 
     Info<< "Calculating field g.h\n" << endl;
     volScalarField gh("gh", g & mesh.C());
+    surfaceScalarField ghf("gh", g & mesh.Cf());
 
     dimensionedScalar pRef("pRef", p.dimensions(), thermo->lookup("pRef"));
 

applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H

@@ -1,19 +1,24 @@
+{
     volScalarField rUA = 1.0/UEqn().A();
+    surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA));
+
     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();
+    surfaceScalarField buoyancyPhi = ghf*fvc::snGrad(rho)*rhorUAf*mesh.magSf();
+    phi -= buoyancyPhi;
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
         fvScalarMatrix pdEqn
         (
-        fvm::laplacian(rho*rUA, pd) == fvc::div(phi)
+            fvm::laplacian(rhorUAf, pd) == fvc::div(phi)
         );
 
         pdEqn.setReference(pdRefCell, pdRefValue);
+
         // retain the residual from the first iteration
         if (nonOrth == 0)
         {
@@ -27,19 +32,22 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
 
         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();
         }
     }
 
     #include "continuityErrs.H"
 
-// Explicitly relax pressure for momentum corrector
+    p == pd + rho*gh + pRef;
-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
@@ -47,9 +55,11 @@ 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;
+
+    pd == p - (rho*gh + pRef);
+}

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/createFields.H

@@ -54,6 +54,7 @@
 
     Info<< "Calculating field g.h\n" << endl;
     volScalarField gh("gh", g & mesh.C());
+    surfaceScalarField ghf("gh", g & mesh.Cf());
 
     dimensionedScalar pRef("pRef", p.dimensions(), thermo->lookup("pRef"));
 

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

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

@@ -1,10 +1,25 @@
-    tmp<fvVectorMatrix> UEqn = solveMomentumEquation
+    // Solve the Momentum equation
+    tmp<fvVectorMatrix> UEqn
     (
-        momentumPredictor,
+        fvm::ddt(rho, U)
-        Uf[i],
+      + fvm::div(phi, U)
-        rhof[i],
+      + turb.divDevRhoReff(U)
-        phif[i],
-        pdf[i],
-        ghf[i],
-        turb[i]
     );
+
+    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/compressibleCourantNo.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,13 +22,12 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-Description
-    Calculates and outputs the mean and maximum Courant Numbers for the fluid
-    regions
-
 \*---------------------------------------------------------------------------*/
 
-scalar compressibleCourantNo
+#include "compressibleCourantNo.H"
+#include "fvc.H"
+
+Foam::scalar Foam::compressibleCourantNo
 (
     const fvMesh& mesh,
     const Time& runTime,

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

@@ -7,8 +7,8 @@
             (
                 fluidRegions[regionI],
                 runTime,
-                rhof[regionI],
+                rhoFluid[regionI],
-                phif[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<basicThermo> thermoFluid(fluidRegions.size());
-    PtrList<volScalarField> rhof(fluidRegions.size());
+    PtrList<volScalarField> rhoFluid(fluidRegions.size());
-    PtrList<volScalarField> Kf(fluidRegions.size());
+    PtrList<volScalarField> KFluid(fluidRegions.size());
-    PtrList<volVectorField> Uf(fluidRegions.size());
+    PtrList<volVectorField> UFluid(fluidRegions.size());
-    PtrList<surfaceScalarField> phif(fluidRegions.size());
+    PtrList<surfaceScalarField> phiFluid(fluidRegions.size());
-    PtrList<compressible::turbulenceModel> turb(fluidRegions.size());
+    PtrList<compressible::turbulenceModel> turbulence(fluidRegions.size());
-    PtrList<volScalarField> DpDtf(fluidRegions.size());
+    PtrList<volScalarField> DpDtFluid(fluidRegions.size());
-    PtrList<volScalarField> ghf(fluidRegions.size());
+    PtrList<volScalarField> ghFluid(fluidRegions.size());
-    PtrList<volScalarField> pdf(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;
+        Info<< "    Adding to pdFluid\n" << endl;
-        pdf.set
+        pdFluid.set
         (
             i,
             new volScalarField
@@ -42,16 +43,15 @@
             )
         );
 
-        Info<< "    Adding to thermof\n" << endl;
+        Info<< "    Adding to thermoFluid\n" << endl;
-
+        thermoFluid.set
-        thermof.set
         (
             i,
             basicThermo::New(fluidRegions[i]).ptr()
         );
 
-        Info<< "    Adding to rhof\n" << endl;
+        Info<< "    Adding to rhoFluid\n" << endl;
-        rhof.set
+        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;
+        Info<< "    Adding to KFluid\n" << endl;
-        Kf.set
+        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;
+        Info<< "    Adding to UFluid\n" << endl;
-        Uf.set
+        UFluid.set
         (
             i,
             new volVectorField
@@ -104,8 +104,8 @@
             )
         );
 
-        Info<< "    Adding to phif\n" << endl;
+        Info<< "    Adding to phiFluid\n" << endl;
-        phif.set
+        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;
+        Info<< "    Adding to turbulence\n" << endl;
-        turb.set
+        turbulence.set
         (
             i,
             autoPtr<compressible::turbulenceModel>
             (
                 compressible::turbulenceModel::New
                 (
-                    rhof[i],
+                    rhoFluid[i],
-                    Uf[i],
+                    UFluid[i],
-                    phif[i],
+                    phiFluid[i],
-                    thermof[i]
+                    thermoFluid[i]
                 )
             ).ptr()
         );
 
-        Info<< "    Adding to DpDtf\n" << endl;
+        Info<< "    Adding to DpDtFluid\n" << endl;
-        DpDtf.set
+        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;
+        Info<< "    Adding to ghFluid\n" << endl;
-        ghf.set
+        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;
+        thermoFluid[i].p() == pdFluid[i] + rhoFluid[i]*ghFluid[i] + pRef;
-        thermof[i].correct();
+        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],
+        fvm::ddt(rho, h)
-        DpDtf[i],
+      + fvm::div(phi, h)
-        phif[i],
+      - fvm::laplacian(turb.alphaEff(), h)
-        turb[i],
+     ==
-        thermof[i]
+        DpDt
     );
+    hEqn().relax();
+    hEqn().solve();
 
+    thermo.correct();
+
+    Info<< "Min/max T:" << min(thermo.T()) << ' ' << max(thermo.T())
+        << endl;
+}

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] =
+    U = rUA*UEqn().H();
-        fvc::interpolate(rhof[i])
+
+    surfaceScalarField phiU
+    (
+        fvc::interpolate(rho)
        *(
-            (fvc::interpolate(Uf[i]) & fluidRegions[i].Sf())
+            (fvc::interpolate(U) & mesh.Sf())
-          + fvc::ddtPhiCorr(rUA, rhof[i], Uf[i], phif[i])
+          + fvc::ddtPhiCorr(rUA, rho, U, phi)
         )
-      - fvc::interpolate(rhof[i]*rUA*ghf[i])
+    );
-       *fvc::snGrad(rhof[i])
-       *fluidRegions[i].magSf();
 
-    // Solve pressure difference
+    phi = phiU - ghf*fvc::snGrad(rho)*rhorUAf*mesh.magSf();
-#   include "pdEqn.H"
+
+    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]);
+#   include "compressibleContinuityErrors.H"
-
-    // Correct velocity field
-    Uf[i] -= rUA*(fvc::grad(pdf[i]) + fvc::grad(rhof[i])*ghf[i]);
-    Uf[i].correctBoundaryConditions();
 
     // For closed-volume cases adjust the pressure and density levels
     // to obey overall mass continuity
     if (closedVolume)
     {
-        thermof[i].p() +=
+        p += (massIni - fvc::domainIntegrate(psi*p))/fvc::domainIntegrate(psi);
-            (
+        rho = thermo.rho();
-                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();
     }
 
     // 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/solveContinuityEquation.C

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

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

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H

@@ -3,10 +3,9 @@
     {
         solve
         (
-            fvm::ddt(rhosCps[i], Ts[i]) - fvm::laplacian(Ks[i], Ts[i])
+            fvm::ddt(rho*cp, T) - fvm::laplacian(K, T)
         );
     }
 
-    Info<< "Min/max T:" << min(Ts[i]) << ' ' << max(Ts[i])
+    Info<< "Min/max T:" << min(T) << ' ' << max(T) << endl;
-        << endl;
 }

applications/solvers/multiphase/interFoam/alphaEqn.H

@@ -6,7 +6,7 @@
     phic = min(interface.cAlpha()*phic, max(phic));
     surfaceScalarField phir = phic*interface.nHatf();
 
-    for (int gCorr=0; gCorr<nAlphaCorr; gCorr++)
+    for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
     {
         surfaceScalarField phiAlpha =
             fvc::flux

src/turbulenceModels/compressible/RAS/RASModel/RASModel.C

@@ -88,7 +88,7 @@ RASModel::RASModel
         dimensioned<scalar>::lookupOrAddToDict
         (
             "kappa",
-            subDict("wallFunctionCoeffs"),
+            wallFunctionDict_,
             0.4187
         )
     ),
@@ -97,7 +97,7 @@ RASModel::RASModel
         dimensioned<scalar>::lookupOrAddToDict
         (
             "E",
-            subDict("wallFunctionCoeffs"),
+            wallFunctionDict_
             9.0
         )
     ),