fireFoam: Updated solver and tutorials from OpenFOAM-1.7.x

applications/solvers/combustion/fireFoam/UEqn.H

@@ -13,7 +13,9 @@ solve
   ==
     fvc::reconstruct
     (
-        fvc::interpolate(rho)*(g & mesh.Sf())
-      - fvc::snGrad(p)*mesh.magSf()
+        (
+          - ghf*fvc::snGrad(rho)
+          - fvc::snGrad(p_rgh)
+        )*mesh.magSf()
     )
 );

applications/solvers/combustion/fireFoam/combustionModels/Make/files

@@ -1,5 +1,5 @@
 combustionModel/combustionModel.C
-combustionModel/combustionModelNew.C
+combustionModel/newCombustionModel.C
 
 infinitelyFastChemistry/infinitelyFastChemistry.C
 

applications/solvers/combustion/fireFoam/combustionModels/combustionModel/newCombustionModel.C

@@ -29,19 +29,22 @@ License
 
 Foam::autoPtr<Foam::combustionModel> Foam::combustionModel::New
 (
-    const dictionary& propDict,
+    const dictionary& combustionProperties,
     const hsCombustionThermo& thermo,
     const compressible::turbulenceModel& turbulence,
     const surfaceScalarField& phi,
     const volScalarField& rho
 )
 {
-    const word modelType(propDict.lookup("combustionModel"));
+    word combustionModelTypeName = combustionProperties.lookup
+    (
+        "combustionModel"
+    );
 
-    Info<< "Selecting combustion model " << modelType << endl;
+    Info<< "Selecting combustion model " << combustionModelTypeName << endl;
 
     dictionaryConstructorTable::iterator cstrIter =
-        dictionaryConstructorTablePtr_->find(modelType);
+        dictionaryConstructorTablePtr_->find(combustionModelTypeName);
 
     if (cstrIter == dictionaryConstructorTablePtr_->end())
     {
@@ -49,14 +52,14 @@ Foam::autoPtr<Foam::combustionModel> Foam::combustionModel::New
         (
             "combustionModel::New"
         )   << "Unknown combustionModel type "
-            << modelType << nl << nl
-            << "Valid combustionModels are : " << endl
-            << dictionaryConstructorTablePtr_->toc()
+            << combustionModelTypeName << endl << endl
+            << "Valid  combustionModels are : " << endl
+            << dictionaryConstructorTablePtr_->sortedToc()
             << exit(FatalError);
     }
 
     return autoPtr<combustionModel>
-        (cstrIter()(propDict, thermo, turbulence, phi, rho));
+        (cstrIter()(combustionProperties, thermo, turbulence, phi, rho));
 }
 
 

applications/solvers/combustion/fireFoam/combustionModels/infinitelyFastChemistry/infinitelyFastChemistry.H

@@ -88,8 +88,9 @@ public:
         );
 
 
-    //- Destructor
-    virtual ~infinitelyFastChemistry();
+    // Destructor
+
+        virtual ~infinitelyFastChemistry();
 
 
     // Member Functions

applications/solvers/combustion/fireFoam/combustionModels/noCombustion/noCombustion.H

@@ -83,8 +83,9 @@ public:
         );
 
 
-    //- Destructor
-    virtual ~noCombustion();
+    // Destructor
+
+        virtual ~noCombustion();
 
     // Member Functions
 

applications/solvers/combustion/fireFoam/createFields.H

@@ -35,6 +35,7 @@ const volScalarField& psi = thermo.psi();
 volScalarField& ft = composition.Y("ft");
 volScalarField& fu = composition.Y("fu");
 
+
 Info<< "Reading field U\n" << endl;
 
 volVectorField U
@@ -73,7 +74,7 @@ IOdictionary combustionProperties
 Info<< "Creating combustion model\n" << endl;
 autoPtr<combustionModel> combustion
 (
-    combustionModel::New
+    combustionModel::combustionModel::New
     (
         combustionProperties,
         thermo,
@@ -83,6 +84,29 @@ autoPtr<combustionModel> combustion
     )
 );
 
+
+Info<< "Calculating field g.h\n" << endl;
+volScalarField gh("gh", g & mesh.C());
+surfaceScalarField ghf("gh", g & mesh.Cf());
+
+Info<< "Reading field p_rgh\n" << endl;
+volScalarField p_rgh
+(
+    IOobject
+    (
+        "p_rgh",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+// Force p_rgh to be consistent with p
+p_rgh = p - rho*gh;
+
+
 volScalarField dQ
 (
     IOobject
@@ -103,15 +127,6 @@ volScalarField DpDt =
     fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
 
 
-Info<< "Calculating field g.h\n" << endl;
-volScalarField gh("gh", g & mesh.C());
-
-surfaceScalarField ghf("gh", g & mesh.Cf());
-
-p += rho*gh;
-
-thermo.correct();
-
 dimensionedScalar initialMass = fvc::domainIntegrate(rho);
 
 

applications/solvers/combustion/fireFoam/pEqn.H

@@ -1,5 +1,3 @@
-bool closedVolume = false;
-
 rho = thermo.rho();
 
 volScalarField rUA = 1.0/UEqn.A();
@@ -15,43 +13,41 @@ surfaceScalarField phiU
     )
 );
 
-phi = phiU + rhorUAf*fvc::interpolate(rho)*(g & mesh.Sf());
+phi = phiU - rhorUAf*ghf*fvc::snGrad(rho)*mesh.magSf();
 
 for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
 {
-       surfaceScalarField rhorUAf = fvc::interpolate(rho*rUA);
+    surfaceScalarField rhorUAf = fvc::interpolate(rho*rUA);
 
-        fvScalarMatrix pEqn
+    fvScalarMatrix p_rghEqn
+    (
+        fvm::ddt(psi, p_rgh) + fvc::ddt(psi, rho)*gh
+      + fvc::div(phi)
+      - fvm::laplacian(rhorUAf, p_rgh)
+    );
+
+    p_rghEqn.solve
+    (
+        mesh.solver
         (
-            fvm::ddt(psi,p)
-          + fvc::div(phi)
-          - fvm::laplacian(rhorUAf, p)
-        );
-
-        closedVolume = p.needReference();
-
-        pEqn.solve
-        (
-            mesh.solver
+            p_rgh.select
             (
-                p.select
                 (
-                    (
-                        finalIter
-                     && corr == nCorr-1
-                     && nonOrth == nNonOrthCorr
-                    )
+                    finalIter
+                 && corr == nCorr-1
+                 && nonOrth == nNonOrthCorr
                 )
             )
-        );
+        )
+    );
 
-        if (nonOrth == nNonOrthCorr)
-        {
-            phi += pEqn.flux();
-        }
+    if (nonOrth == nNonOrthCorr)
+    {
+        phi += p_rghEqn.flux();
+    }
 }
 
-DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
+p = p_rgh + rho*gh;
 
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
@@ -59,12 +55,4 @@ DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
 U += rUA*fvc::reconstruct((phi - phiU)/rhorUAf);
 U.correctBoundaryConditions();
 
-// For closed-volume cases adjust the pressure and density levels
-// to obey overall mass continuity
-if (closedVolume)
-{
-    p +=
-        (initialMass - fvc::domainIntegrate(thermo.psi()*p))
-       /fvc::domainIntegrate(thermo.psi());
-    rho = thermo.rho();
-}
+DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);