diff --git a/applications/solvers/combustion/PDRFoam/EaEqn.H b/applications/solvers/combustion/PDRFoam/EaEqn.H new file mode 100644 index 0000000000..03d13f2693 --- /dev/null +++ b/applications/solvers/combustion/PDRFoam/EaEqn.H @@ -0,0 +1,17 @@ +{ + volScalarField& hea = thermo.he(); + + solve + ( + betav*fvm::ddt(rho, hea) + mvConvection->fvmDiv(phi, hea) + + betav*fvc::ddt(rho, K) + fvc::div(phi, K) + + ( + hea.name() == "ea" + ? fvc::div(phi, volScalarField("Ep", p/rho)) + : -betav*dpdt + ) + - fvm::laplacian(Db, hea) + ); + + thermo.correct(); +} diff --git a/applications/solvers/combustion/PDRFoam/EauEqn.H b/applications/solvers/combustion/PDRFoam/EauEqn.H new file mode 100644 index 0000000000..96412cc3c3 --- /dev/null +++ b/applications/solvers/combustion/PDRFoam/EauEqn.H @@ -0,0 +1,22 @@ +if (ign.ignited()) +{ + volScalarField& heau = thermo.heu(); + + solve + ( + betav*fvm::ddt(rho, heau) + mvConvection->fvmDiv(phi, heau) + + (betav*fvc::ddt(rho, K) + fvc::div(phi, K))*rho/thermo.rhou() + + ( + heau.name() == "eau" + ? fvc::div(phi, volScalarField("Ep", p/rho))*rho/thermo.rhou() + : -betav*dpdt*rho/thermo.rhou() + ) + - fvm::laplacian(Db, heau) + + // These terms cannot be used in partially-premixed combustion due to + // the resultant inconsistency between ft and heau transport. + // A possible solution would be to solve for ftu as well as ft. + //- fvm::div(muEff*fvc::grad(b)/(b + 0.001), heau) + //+ fvm::Sp(fvc::div(muEff*fvc::grad(b)/(b + 0.001)), heau) + ); +} diff --git a/applications/solvers/combustion/PDRFoam/PDRFoam.C b/applications/solvers/combustion/PDRFoam/PDRFoam.C index ae63ee80c8..ad6e2d0298 100644 --- a/applications/solvers/combustion/PDRFoam/PDRFoam.C +++ b/applications/solvers/combustion/PDRFoam/PDRFoam.C @@ -123,12 +123,12 @@ int main(int argc, char *argv[]) { #include "bEqn.H" #include "ftEqn.H" - #include "hauEqn.H" - #include "haEqn.H" + #include "EauEqn.H" + #include "EaEqn.H" if (!ign.ignited()) { - hau == ha; + thermo.heu() == thermo.he(); } #include "pEqn.H" diff --git a/applications/solvers/combustion/PDRFoam/createFields.H b/applications/solvers/combustion/PDRFoam/createFields.H index 64d3058c0e..c94bc4fe81 100644 --- a/applications/solvers/combustion/PDRFoam/createFields.H +++ b/applications/solvers/combustion/PDRFoam/createFields.H @@ -5,6 +5,8 @@ psiuReactionThermo::New(mesh) ); psiuReactionThermo& thermo = pThermo(); + thermo.validate(args.executable(), "ha", "ea"); + basicMultiComponentMixture& composition = thermo.composition(); volScalarField rho @@ -22,14 +24,10 @@ volScalarField& p = thermo.p(); const volScalarField& psi = thermo.psi(); - volScalarField& ha = thermo.he(); - volScalarField& hau = thermo.heu(); volScalarField& b = composition.Y("b"); Info<< "min(b) = " << min(b).value() << endl; - //const volScalarField& T = thermo->T(); - Info<< "\nReading field U\n" << endl; volVectorField U ( @@ -196,6 +194,6 @@ } fields.add(b); - fields.add(ha); - fields.add(hau); + fields.add(thermo.he()); + fields.add(thermo.heu()); flameWrinkling->addXi(fields); diff --git a/applications/solvers/combustion/PDRFoam/haEqn.H b/applications/solvers/combustion/PDRFoam/haEqn.H deleted file mode 100644 index 20f605f79c..0000000000 --- a/applications/solvers/combustion/PDRFoam/haEqn.H +++ /dev/null @@ -1,13 +0,0 @@ -{ - solve - ( - betav*fvm::ddt(rho, ha) - + mvConvection->fvmDiv(phi, ha) - - fvm::laplacian(Db, ha) - == - betav*dpdt - - betav*(fvc::ddt(rho, K) + fvc::div(phi, K)) - ); - - thermo.correct(); -} diff --git a/applications/solvers/combustion/PDRFoam/hauEqn.H b/applications/solvers/combustion/PDRFoam/hauEqn.H deleted file mode 100644 index 1cc1a1fd35..0000000000 --- a/applications/solvers/combustion/PDRFoam/hauEqn.H +++ /dev/null @@ -1,18 +0,0 @@ -if (ign.ignited()) -{ - solve - ( - betav*fvm::ddt(rho, hau) - + mvConvection->fvmDiv(phi, hau) - - fvm::laplacian(Db, hau) - - // These terms cannot be used in partially-premixed combustion due to - // the resultant inconsistency between ft and hau transport. - // A possible solution would be to solve for ftu as well as ft. - //- fvm::div(muEff*fvc::grad(b)/(b + 0.001), hau) - //+ fvm::Sp(fvc::div(muEff*fvc::grad(b)/(b + 0.001)), hau) - - == - betav*(dpdt - (fvc::ddt(rho, K) + fvc::div(phi, K)))*rho/thermo.rhou() - ); -} diff --git a/applications/solvers/combustion/XiFoam/EaEqn.H b/applications/solvers/combustion/XiFoam/EaEqn.H new file mode 100644 index 0000000000..a9e72b4627 --- /dev/null +++ b/applications/solvers/combustion/XiFoam/EaEqn.H @@ -0,0 +1,20 @@ +{ + volScalarField& hea = thermo.he(); + + fvScalarMatrix EaEqn + ( + fvm::ddt(rho, hea) + mvConvection->fvmDiv(phi, hea) + + fvc::ddt(rho, K) + fvc::div(phi, K) + + ( + hea.name() == "ea" + ? fvc::div(phi, volScalarField("Ep", p/rho)) + : -dpdt + ) + - fvm::laplacian(turbulence->alphaEff(), hea) + ); + + EaEqn.relax(); + EaEqn.solve(); + + thermo.correct(); +} diff --git a/applications/solvers/combustion/XiFoam/EauEqn.H b/applications/solvers/combustion/XiFoam/EauEqn.H new file mode 100644 index 0000000000..0d81bd6e60 --- /dev/null +++ b/applications/solvers/combustion/XiFoam/EauEqn.H @@ -0,0 +1,22 @@ +if (ign.ignited()) +{ + volScalarField& heau = thermo.heu(); + + solve + ( + fvm::ddt(rho, heau) + mvConvection->fvmDiv(phi, heau) + + (fvc::ddt(rho, K) + fvc::div(phi, K))*rho/thermo.rhou() + + ( + heau.name() == "eau" + ? fvc::div(phi, volScalarField("Ep", p/rho))*rho/thermo.rhou() + : -dpdt*rho/thermo.rhou() + ) + - fvm::laplacian(turbulence->alphaEff(), heau) + + // These terms cannot be used in partially-premixed combustion due to + // the resultant inconsistency between ft and heau transport. + // A possible solution would be to solve for ftu as well as ft. + //- fvm::div(muEff*fvc::grad(b)/(b + 0.001), heau) + //+ fvm::Sp(fvc::div(muEff*fvc::grad(b)/(b + 0.001)), heau) + ); +} diff --git a/applications/solvers/combustion/XiFoam/XiFoam.C b/applications/solvers/combustion/XiFoam/XiFoam.C index 874af9a1fd..1573ede331 100644 --- a/applications/solvers/combustion/XiFoam/XiFoam.C +++ b/applications/solvers/combustion/XiFoam/XiFoam.C @@ -97,12 +97,12 @@ int main(int argc, char *argv[]) #include "ftEqn.H" #include "bEqn.H" - #include "hauEqn.H" - #include "haEqn.H" + #include "EauEqn.H" + #include "EaEqn.H" if (!ign.ignited()) { - hau == ha; + thermo.heu() == thermo.he(); } // --- Pressure corrector loop diff --git a/applications/solvers/combustion/XiFoam/createFields.H b/applications/solvers/combustion/XiFoam/createFields.H index c266010adb..17103885e5 100644 --- a/applications/solvers/combustion/XiFoam/createFields.H +++ b/applications/solvers/combustion/XiFoam/createFields.H @@ -5,6 +5,8 @@ psiuReactionThermo::New(mesh) ); psiuReactionThermo& thermo = pThermo(); + thermo.validate(args.executable(), "ha", "ea"); + basicMultiComponentMixture& composition = thermo.composition(); volScalarField rho @@ -22,14 +24,10 @@ volScalarField& p = thermo.p(); const volScalarField& psi = thermo.psi(); - volScalarField& ha = thermo.he(); - volScalarField& hau = thermo.heu(); volScalarField& b = composition.Y("b"); Info<< "min(b) = " << min(b).value() << endl; - const volScalarField& T = thermo.T(); - Info<< "\nReading field U\n" << endl; volVectorField U @@ -138,5 +136,5 @@ } fields.add(b); - fields.add(ha); - fields.add(hau); + fields.add(thermo.he()); + fields.add(thermo.heu()); diff --git a/applications/solvers/combustion/XiFoam/haEqn.H b/applications/solvers/combustion/XiFoam/haEqn.H deleted file mode 100644 index ffbac948ae..0000000000 --- a/applications/solvers/combustion/XiFoam/haEqn.H +++ /dev/null @@ -1,16 +0,0 @@ -{ - fvScalarMatrix haEqn - ( - fvm::ddt(rho, ha) - + mvConvection->fvmDiv(phi, ha) - - fvm::laplacian(turbulence->alphaEff(), ha) - == - dpdt - - (fvc::ddt(rho, K) + fvc::div(phi, K)) - ); - - haEqn.relax(); - haEqn.solve(); - - thermo.correct(); -} diff --git a/applications/solvers/combustion/XiFoam/hauEqn.H b/applications/solvers/combustion/XiFoam/hauEqn.H deleted file mode 100644 index d9f5767643..0000000000 --- a/applications/solvers/combustion/XiFoam/hauEqn.H +++ /dev/null @@ -1,18 +0,0 @@ -if (ign.ignited()) -{ - solve - ( - fvm::ddt(rho, hau) - + mvConvection->fvmDiv(phi, hau) - - fvm::laplacian(turbulence->alphaEff(), hau) - - // These terms cannot be used in partially-premixed combustion due to - // the resultant inconsistency between ft and hau transport. - // A possible solution would be to solve for ftu as well as ft. - //- fvm::div(muEff*fvc::grad(b)/(b + 0.001), hau) - //+ fvm::Sp(fvc::div(muEff*fvc::grad(b)/(b + 0.001)), hau) - - == - (dpdt - (fvc::ddt(rho, K) + fvc::div(phi, K)))*rho/thermo.rhou() - ); -} diff --git a/applications/solvers/combustion/chemFoam/createFields.H b/applications/solvers/combustion/chemFoam/createFields.H index 3beca02bec..8855e4501a 100644 --- a/applications/solvers/combustion/chemFoam/createFields.H +++ b/applications/solvers/combustion/chemFoam/createFields.H @@ -30,6 +30,8 @@ scalar dtChem = refCast(chemistry).deltaTChem()[0]; psiReactionThermo& thermo = chemistry.thermo(); + thermo.validate(args.executable(), "h"); + basicMultiComponentMixture& composition = thermo.composition(); PtrList& Y = composition.Y(); @@ -47,7 +49,6 @@ ); volScalarField& p = thermo.p(); - volScalarField& hs = thermo.he(); volScalarField Rspecific ( diff --git a/applications/solvers/combustion/chemFoam/hEqn.H b/applications/solvers/combustion/chemFoam/hEqn.H index 855ca66a42..2478e6017c 100644 --- a/applications/solvers/combustion/chemFoam/hEqn.H +++ b/applications/solvers/combustion/chemFoam/hEqn.H @@ -1,10 +1,12 @@ { + volScalarField& h = thermo.he(); + if (constProp == "volume") { - hs[0] = u0 + p[0]/rho[0] + integratedHeat; + h[0] = u0 + p[0]/rho[0] + integratedHeat; } else { - hs[0] = hs0 + integratedHeat; + h[0] = h0 + integratedHeat; } -} \ No newline at end of file +} diff --git a/applications/solvers/combustion/chemFoam/readInitialConditions.H b/applications/solvers/combustion/chemFoam/readInitialConditions.H index 9a03c53b1f..d929b14128 100644 --- a/applications/solvers/combustion/chemFoam/readInitialConditions.H +++ b/applications/solvers/combustion/chemFoam/readInitialConditions.H @@ -83,20 +83,19 @@ } } - scalar hs0 = 0.0; + scalar h0 = 0.0; forAll(Y, i) { Y[i] = Y0[i]; - hs0 += Y0[i]*specieData[i].Hs(p[i], T0); + h0 += Y0[i]*specieData[i].Hs(p[i], T0); } - hs = dimensionedScalar("h", dimEnergy/dimMass, hs0); - + thermo.he() = dimensionedScalar("h", dimEnergy/dimMass, h0); thermo.correct(); rho = thermo.rho(); scalar rho0 = rho[0]; - scalar u0 = hs0 - p0/rho0; + scalar u0 = h0 - p0/rho0; scalar R0 = p0/(rho0*T0); Rspecific[0] = R0; diff --git a/applications/solvers/combustion/coldEngineFoam/Make/options b/applications/solvers/combustion/coldEngineFoam/Make/options index 5af1f2fd61..a59f3b8053 100644 --- a/applications/solvers/combustion/coldEngineFoam/Make/options +++ b/applications/solvers/combustion/coldEngineFoam/Make/options @@ -1,6 +1,7 @@ EXE_INC = \ -I../engineFoam \ -I../XiFoam \ + -I../../compressible/rhoPimpleFoam \ -I$(LIB_SRC)/engine/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ diff --git a/applications/solvers/combustion/coldEngineFoam/coldEngineFoam.C b/applications/solvers/combustion/coldEngineFoam/coldEngineFoam.C index b517f93c9b..87474e778f 100644 --- a/applications/solvers/combustion/coldEngineFoam/coldEngineFoam.C +++ b/applications/solvers/combustion/coldEngineFoam/coldEngineFoam.C @@ -81,7 +81,7 @@ int main(int argc, char *argv[]) // --- Pressure corrector loop while (pimple.correct()) { - #include "hEqn.H" + #include "EEqn.H" #include "pEqn.H" } diff --git a/applications/solvers/combustion/coldEngineFoam/createFields.H b/applications/solvers/combustion/coldEngineFoam/createFields.H index f97c40eab9..f78e5bde29 100644 --- a/applications/solvers/combustion/coldEngineFoam/createFields.H +++ b/applications/solvers/combustion/coldEngineFoam/createFields.H @@ -5,6 +5,7 @@ psiThermo::New(mesh) ); psiThermo& thermo = pThermo(); + thermo.validate(args.executable(), "h", "e"); volScalarField rho ( @@ -21,7 +22,6 @@ volScalarField& p = thermo.p(); const volScalarField& psi = thermo.psi(); - volScalarField& h = thermo.he(); const volScalarField& T = thermo.T(); diff --git a/applications/solvers/combustion/coldEngineFoam/hEqn.H b/applications/solvers/combustion/coldEngineFoam/hEqn.H deleted file mode 100644 index 2e1a2703a8..0000000000 --- a/applications/solvers/combustion/coldEngineFoam/hEqn.H +++ /dev/null @@ -1,12 +0,0 @@ -{ - solve - ( - fvm::ddt(rho, h) - + fvm::div(phi, h) - - fvm::laplacian(turbulence->alphaEff(), h) - == - - fvc::div(phi, 0.5*magSqr(U)) - ); - - thermo.correct(); -} diff --git a/applications/solvers/combustion/engineFoam/engineFoam.C b/applications/solvers/combustion/engineFoam/engineFoam.C index 267dfb1119..0ee3df6b7b 100644 --- a/applications/solvers/combustion/engineFoam/engineFoam.C +++ b/applications/solvers/combustion/engineFoam/engineFoam.C @@ -103,12 +103,12 @@ int main(int argc, char *argv[]) #include "ftEqn.H" #include "bEqn.H" - #include "hauEqn.H" - #include "haEqn.H" + #include "EauEqn.H" + #include "EaEqn.H" if (!ign.ignited()) { - hau == ha; + thermo.heu() == thermo.he(); } // --- Pressure corrector loop diff --git a/applications/solvers/combustion/engineFoam/logSummary.H b/applications/solvers/combustion/engineFoam/logSummary.H index 24d8812008..c243b889ae 100644 --- a/applications/solvers/combustion/engineFoam/logSummary.H +++ b/applications/solvers/combustion/engineFoam/logSummary.H @@ -1,5 +1,6 @@ Info<< "Mean pressure:" << p.weightedAverage(mesh.V()).value() << endl; -Info<< "Mean temperature:" << T.weightedAverage(mesh.V()).value() << endl; +Info<< "Mean temperature:" << thermo.T().weightedAverage(mesh.V()).value() + << endl; Info<< "Mean u':" << (sqrt((2.0/3.0)*turbulence->k()))().weightedAverage(mesh.V()).value() << endl; @@ -7,7 +8,7 @@ Info<< "Mean u':" logSummaryFile << runTime.theta() << tab << p.weightedAverage(mesh.V()).value() << tab - << T.weightedAverage(mesh.V()).value() << tab + << thermo.T().weightedAverage(mesh.V()).value() << tab << (sqrt((2.0/3.0)*turbulence->k()))().weightedAverage(mesh.V()).value() << tab << 1 - b.weightedAverage(mesh.V()).value() diff --git a/applications/solvers/combustion/fireFoam/YhsEqn.H b/applications/solvers/combustion/fireFoam/YEEqn.H similarity index 74% rename from applications/solvers/combustion/fireFoam/YhsEqn.H rename to applications/solvers/combustion/fireFoam/YEEqn.H index 1319a4d5df..340cfc63b1 100644 --- a/applications/solvers/combustion/fireFoam/YhsEqn.H +++ b/applications/solvers/combustion/fireFoam/YEEqn.H @@ -47,22 +47,27 @@ tmp > mvConvection Y[inertIndex] = scalar(1) - Yt; Y[inertIndex].max(0.0); - fvScalarMatrix hsEqn + volScalarField& he = thermo.he(); + + fvScalarMatrix EEqn ( - fvm::ddt(rho, hs) - + mvConvection->fvmDiv(phi, hs) - - fvm::laplacian(turbulence->alphaEff(), hs) + fvm::ddt(rho, he) + mvConvection->fvmDiv(phi, he) + + fvc::ddt(rho, K) + fvc::div(phi, K) + + ( + he.name() == "e" + ? fvc::div(phi, volScalarField("Ep", p/rho)) + : -dpdt + ) + - fvm::laplacian(turbulence->alphaEff(), he) == - dpdt - - (fvc::ddt(rho, K) + fvc::div(phi, K)) - + combustion->Sh() + combustion->Sh() + radiation->Sh(thermo) - + parcels.Sh(hs) + + parcels.Sh(he) + surfaceFilm.Sh() ); - hsEqn.relax(); - hsEqn.solve(); + EEqn.relax(); + EEqn.solve(); thermo.correct(); diff --git a/applications/solvers/combustion/fireFoam/createFields.H b/applications/solvers/combustion/fireFoam/createFields.H index 873ae28a07..8159ee73aa 100644 --- a/applications/solvers/combustion/fireFoam/createFields.H +++ b/applications/solvers/combustion/fireFoam/createFields.H @@ -11,6 +11,7 @@ Info<< "Reading thermophysical properties\n" << endl; psiReactionThermo& thermo = combustion->thermo(); + thermo.validate(args.executable(), "h", "e"); SLGThermo slgThermo(mesh, thermo); @@ -34,7 +35,6 @@ ); volScalarField& p = thermo.p(); - volScalarField& hs = thermo.he(); const volScalarField& T = thermo.T(); const volScalarField& psi = thermo.psi(); @@ -128,7 +128,7 @@ { fields.add(Y[i]); } - fields.add(hs); + fields.add(thermo.he()); IOdictionary additionalControlsDict ( diff --git a/applications/solvers/combustion/fireFoam/fireFoam.C b/applications/solvers/combustion/fireFoam/fireFoam.C index 1d58470cbf..9f97843a96 100644 --- a/applications/solvers/combustion/fireFoam/fireFoam.C +++ b/applications/solvers/combustion/fireFoam/fireFoam.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -94,7 +94,7 @@ int main(int argc, char *argv[]) while (pimple.loop()) { #include "UEqn.H" - #include "YhsEqn.H" + #include "YEEqn.H" // --- Pressure corrector loop while (pimple.correct()) diff --git a/applications/solvers/combustion/reactingFoam/EEqn.H b/applications/solvers/combustion/reactingFoam/EEqn.H new file mode 100644 index 0000000000..33dbc46287 --- /dev/null +++ b/applications/solvers/combustion/reactingFoam/EEqn.H @@ -0,0 +1,26 @@ +{ + volScalarField& he = thermo.he(); + + fvScalarMatrix EEqn + ( + fvm::ddt(rho, he) + mvConvection->fvmDiv(phi, he) + + fvc::ddt(rho, K) + fvc::div(phi, K) + + ( + he.name() == "e" + ? fvc::div(phi, volScalarField("Ep", p/rho)) + : -dpdt + ) + - fvm::laplacian(turbulence->alphaEff(), he) +// - fvm::laplacian(turbulence->muEff(), he) // unit lewis no. + == + reaction->Sh() + ); + + EEqn.relax(); + EEqn.solve(); + + thermo.correct(); + + Info<< "min/max(T) = " + << min(T).value() << ", " << max(T).value() << endl; +} diff --git a/applications/solvers/combustion/reactingFoam/createFields.H b/applications/solvers/combustion/reactingFoam/createFields.H index 24f4135d7d..32bdd372ee 100644 --- a/applications/solvers/combustion/reactingFoam/createFields.H +++ b/applications/solvers/combustion/reactingFoam/createFields.H @@ -6,6 +6,7 @@ autoPtr reaction ); psiReactionThermo& thermo = reaction->thermo(); +thermo.validate(args.executable(), "h", "e"); basicMultiComponentMixture& composition = thermo.composition(); PtrList& Y = composition.Y(); @@ -40,7 +41,6 @@ volVectorField U volScalarField& p = thermo.p(); const volScalarField& psi = thermo.psi(); -volScalarField& hs = thermo.he(); const volScalarField& T = thermo.T(); #include "compressibleCreatePhi.H" @@ -84,7 +84,7 @@ forAll(Y, i) { fields.add(Y[i]); } -fields.add(hs); +fields.add(thermo.he()); volScalarField dQ ( diff --git a/applications/solvers/combustion/reactingFoam/hsEqn.H b/applications/solvers/combustion/reactingFoam/hsEqn.H deleted file mode 100644 index 5bc7c5b169..0000000000 --- a/applications/solvers/combustion/reactingFoam/hsEqn.H +++ /dev/null @@ -1,21 +0,0 @@ -{ - fvScalarMatrix hsEqn - ( - fvm::ddt(rho, hs) - + mvConvection->fvmDiv(phi, hs) - - fvm::laplacian(turbulence->alphaEff(), hs) -// - fvm::laplacian(turbulence->muEff(), hs) // unit lewis no. - == - dpdt - - (fvc::ddt(rho, K) + fvc::div(phi, K)) - + reaction->Sh() - ); - - hsEqn.relax(); - hsEqn.solve(); - - thermo.correct(); - - Info<< "T gas min/max = " << min(T).value() << ", " - << max(T).value() << endl; -} diff --git a/applications/solvers/combustion/reactingFoam/reactingFoam.C b/applications/solvers/combustion/reactingFoam/reactingFoam.C index f2517529e2..cd7d371651 100644 --- a/applications/solvers/combustion/reactingFoam/reactingFoam.C +++ b/applications/solvers/combustion/reactingFoam/reactingFoam.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -70,7 +70,7 @@ int main(int argc, char *argv[]) { #include "UEqn.H" #include "YEqn.H" - #include "hsEqn.H" + #include "EEqn.H" // --- Pressure corrector loop while (pimple.correct()) diff --git a/applications/solvers/combustion/rhoReactingFoam/Make/options b/applications/solvers/combustion/rhoReactingFoam/Make/options index 8d18c08338..6dbd401b63 100644 --- a/applications/solvers/combustion/rhoReactingFoam/Make/options +++ b/applications/solvers/combustion/rhoReactingFoam/Make/options @@ -1,4 +1,5 @@ EXE_INC = \ + -I../reactingFoam \ -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \ -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \ diff --git a/applications/solvers/combustion/rhoReactingFoam/UEqn.H b/applications/solvers/combustion/rhoReactingFoam/UEqn.H deleted file mode 100644 index b9bc567aae..0000000000 --- a/applications/solvers/combustion/rhoReactingFoam/UEqn.H +++ /dev/null @@ -1,16 +0,0 @@ - fvVectorMatrix UEqn - ( - fvm::ddt(rho, U) - + fvm::div(phi, U) - + turbulence->divDevRhoReff(U) - == - rho*g - ); - - UEqn.relax(); - - if (pimple.momentumPredictor()) - { - solve(UEqn == -fvc::grad(p)); - K = 0.5*magSqr(U); - } diff --git a/applications/solvers/combustion/rhoReactingFoam/YEqn.H b/applications/solvers/combustion/rhoReactingFoam/YEqn.H deleted file mode 100644 index ccc4b135a5..0000000000 --- a/applications/solvers/combustion/rhoReactingFoam/YEqn.H +++ /dev/null @@ -1,47 +0,0 @@ -tmp > mvConvection -( - fv::convectionScheme::New - ( - mesh, - fields, - phi, - mesh.divScheme("div(phi,Yi_h)") - ) -); - -{ - reaction->correct(); - dQ = reaction->dQ(); - label inertIndex = -1; - volScalarField Yt(0.0*Y[0]); - - forAll(Y, i) - { - if (Y[i].name() != inertSpecie) - { - volScalarField& Yi = Y[i]; - - fvScalarMatrix YiEqn - ( - fvm::ddt(rho, Yi) - + mvConvection->fvmDiv(phi, Yi) - - fvm::laplacian(turbulence->muEff(), Yi) - == - reaction->R(Yi) - ); - - YiEqn.relax(); - YiEqn.solve(mesh.solver("Yi")); - - Yi.max(0.0); - Yt += Yi; - } - else - { - inertIndex = i; - } - } - - Y[inertIndex] = scalar(1) - Yt; - Y[inertIndex].max(0.0); -} diff --git a/applications/solvers/combustion/rhoReactingFoam/createFields.H b/applications/solvers/combustion/rhoReactingFoam/createFields.H index c0b1a76381..e99639e189 100644 --- a/applications/solvers/combustion/rhoReactingFoam/createFields.H +++ b/applications/solvers/combustion/rhoReactingFoam/createFields.H @@ -6,6 +6,7 @@ autoPtr reaction ); rhoReactionThermo& thermo = reaction->thermo(); +thermo.validate(args.executable(), "h", "e"); basicMultiComponentMixture& composition = thermo.composition(); PtrList& Y = composition.Y(); @@ -40,7 +41,6 @@ volVectorField U volScalarField& p = thermo.p(); const volScalarField& psi = thermo.psi(); -volScalarField& hs = thermo.he(); const volScalarField& T = thermo.T(); @@ -86,7 +86,7 @@ forAll(Y, i) { fields.add(Y[i]); } -fields.add(hs); +fields.add(thermo.he()); volScalarField dQ ( diff --git a/applications/solvers/combustion/rhoReactingFoam/hsEqn.H b/applications/solvers/combustion/rhoReactingFoam/hsEqn.H deleted file mode 100644 index f92cfd140b..0000000000 --- a/applications/solvers/combustion/rhoReactingFoam/hsEqn.H +++ /dev/null @@ -1,21 +0,0 @@ -{ - fvScalarMatrix hsEqn - ( - fvm::ddt(rho, hs) - + mvConvection->fvmDiv(phi, hs) - - fvm::laplacian(turbulence->alphaEff(), hs) -// - fvm::laplacian(turbulence->muEff(), hs) // unit lewis no. - == - dpdt - - (fvc::ddt(rho, K) + fvc::div(phi, K)) - + reaction->Sh() - ); - - hsEqn.relax(); - hsEqn.solve(); - - thermo.correct(); - - Info<< "min/max(T) = " - << min(T).value() << ", " << max(T).value() << endl; -} diff --git a/applications/solvers/combustion/rhoReactingFoam/rhoReactingFoam.C b/applications/solvers/combustion/rhoReactingFoam/rhoReactingFoam.C index f3bab92f2f..3ac060612d 100644 --- a/applications/solvers/combustion/rhoReactingFoam/rhoReactingFoam.C +++ b/applications/solvers/combustion/rhoReactingFoam/rhoReactingFoam.C @@ -72,7 +72,7 @@ int main(int argc, char *argv[]) { #include "UEqn.H" #include "YEqn.H" - #include "hsEqn.H" + #include "EEqn.H" // --- Pressure corrector loop while (pimple.correct()) diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/EEqn.H b/applications/solvers/compressible/rhoPimpleFoam/EEqn.H similarity index 100% rename from applications/solvers/heatTransfer/buoyantPimpleFoam/EEqn.H rename to applications/solvers/compressible/rhoPimpleFoam/EEqn.H diff --git a/applications/solvers/compressible/rhoPimpleFoam/UEqn.H b/applications/solvers/compressible/rhoPimpleFoam/UEqn.H index b7b4b04db5..0f5ec2487b 100644 --- a/applications/solvers/compressible/rhoPimpleFoam/UEqn.H +++ b/applications/solvers/compressible/rhoPimpleFoam/UEqn.H @@ -4,7 +4,6 @@ tmp UEqn ( fvm::ddt(rho, U) + fvm::div(phi, U) - - fvm::Sp(fvc::ddt(rho) + fvc::div(phi), U) + turbulence->divDevRhoReff(U) ); diff --git a/applications/solvers/compressible/rhoPimpleFoam/createFields.H b/applications/solvers/compressible/rhoPimpleFoam/createFields.H index 8d3a6ad21d..67cc0c3e45 100644 --- a/applications/solvers/compressible/rhoPimpleFoam/createFields.H +++ b/applications/solvers/compressible/rhoPimpleFoam/createFields.H @@ -5,9 +5,9 @@ psiThermo::New(mesh) ); psiThermo& thermo = pThermo(); + thermo.validate(args.executable(), "h", "e"); volScalarField& p = thermo.p(); - volScalarField& h = thermo.he(); const volScalarField& psi = thermo.psi(); volScalarField rho diff --git a/applications/solvers/compressible/rhoPimpleFoam/hEqn.H b/applications/solvers/compressible/rhoPimpleFoam/hEqn.H deleted file mode 100644 index 57d72887e6..0000000000 --- a/applications/solvers/compressible/rhoPimpleFoam/hEqn.H +++ /dev/null @@ -1,20 +0,0 @@ -{ - fvScalarMatrix hEqn - ( - fvm::ddt(rho, h) - + fvm::div(phi, h) - - fvm::Sp(fvc::ddt(rho) + fvc::div(phi), h) - - fvm::laplacian(turbulence->alphaEff(), h) - == - dpdt - - ( - fvc::ddt(rho, K) + fvc::div(phi, K) - - (fvc::ddt(rho) + fvc::div(phi))*K - ) - ); - - hEqn.relax(); - hEqn.solve(); - - thermo.correct(); -} diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C index 4a4a82db7c..adb6882268 100644 --- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C +++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C @@ -75,7 +75,7 @@ int main(int argc, char *argv[]) while (pimple.loop()) { #include "UEqn.H" - #include "hEqn.H" + #include "EEqn.H" // --- Pressure corrector loop while (pimple.correct()) diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/rhoPimplecFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/rhoPimplecFoam.C index 73a6f5bd2f..cf8aed0bca 100644 --- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/rhoPimplecFoam.C +++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/rhoPimplecFoam.C @@ -75,7 +75,7 @@ int main(int argc, char *argv[]) while (pimple.loop()) { #include "UEqn.H" - #include "hEqn.H" + #include "EEqn.H" // --- Pressure corrector loop while (pimple.correct()) diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/rhoPorousMRFLTSPimpleFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/rhoPorousMRFLTSPimpleFoam.C index ae57bda417..3a8eb884ba 100644 --- a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/rhoPorousMRFLTSPimpleFoam.C +++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/rhoPorousMRFLTSPimpleFoam.C @@ -85,7 +85,7 @@ int main(int argc, char *argv[]) turbulence->correct(); #include "UEqn.H" - #include "hEqn.H" + #include "EEqn.H" // --- Pressure corrector loop while (pimple.correct()) diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C index a5adc28054..16e6c0ad74 100644 --- a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C +++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C @@ -78,7 +78,7 @@ int main(int argc, char *argv[]) while (pimple.loop()) { #include "UEqn.H" - #include "hEqn.H" + #include "EEqn.H" // --- Pressure corrector loop while (pimple.correct()) diff --git a/applications/solvers/compressible/rhoSimpleFoam/UEqn.H b/applications/solvers/compressible/rhoSimpleFoam/UEqn.H index f1bed4d071..21ec2646be 100644 --- a/applications/solvers/compressible/rhoSimpleFoam/UEqn.H +++ b/applications/solvers/compressible/rhoSimpleFoam/UEqn.H @@ -3,7 +3,6 @@ tmp UEqn ( fvm::div(phi, U) - - fvm::Sp(fvc::div(phi), U) + turbulence->divDevRhoReff(U) ); diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/TEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/TEqn.H index 180a33d918..76adf004cb 100644 --- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/TEqn.H +++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/TEqn.H @@ -7,7 +7,6 @@ fvScalarMatrix TEqn ( fvm::div(phi, T) - - fvm::Sp(fvc::div(phi), T) - fvm::laplacian(kappaEff, T) ); diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/Make/options b/applications/solvers/heatTransfer/buoyantPimpleFoam/Make/options index db07a71c86..f6e12a3b7c 100644 --- a/applications/solvers/heatTransfer/buoyantPimpleFoam/Make/options +++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/Make/options @@ -1,4 +1,5 @@ EXE_INC = \ + -I../../compressible/rhoPimpleFoam \ -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \ -I$(LIB_SRC)/finiteVolume/lnInclude diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/EEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/EEqn.H new file mode 100644 index 0000000000..be84cfce0e --- /dev/null +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/EEqn.H @@ -0,0 +1,26 @@ +{ + volScalarField& he = thermo.he(); + + fvScalarMatrix EEqn + ( + fvm::div(phi, he) + + ( + he.name() == "e" + ? fvc::div(phi, volScalarField("Ekp", 0.5*magSqr(U) + p/rho)) + : fvc::div(phi, volScalarField("K", 0.5*magSqr(U))) + ) + - fvm::laplacian(turb.alphaEff(), he) + == + rad.Sh(thermo) + + sources(rho, he) + ); + + EEqn.relax(); + EEqn.solve(); + + thermo.correct(); + rad.correct(); + + Info<< "Min/max T:" << min(thermo.T()).value() << ' ' + << max(thermo.T()).value() << endl; +} diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/hEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/hEqn.H deleted file mode 100644 index e55b6d141d..0000000000 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/hEqn.H +++ /dev/null @@ -1,23 +0,0 @@ -{ - fvScalarMatrix hEqn - ( - fvm::div(phi, h) - - fvm::Sp(fvc::div(phi), h) - - fvm::laplacian(turb.alphaEff(), h) - == - - fvc::div(phi, 0.5*magSqr(U), "div(phi,K)") - + rad.Sh(thermo) - + sources(rho, h) - ); - - hEqn.relax(); - - hEqn.solve(); - - thermo.correct(); - - rad.correct(); - - Info<< "Min/max T:" << min(thermo.T()).value() << ' ' - << max(thermo.T()).value() << endl; -} diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H index e6fc3d6f9a..d2ea510ced 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H @@ -7,23 +7,31 @@ volScalarField rAU(1.0/UEqn().A()); surfaceScalarField rhorAUf("Dp", fvc::interpolate(rho*rAU)); - U = rAU*UEqn().H(); + volVectorField HbyA("HbyA", U); + HbyA = rAU*UEqn().H(); UEqn.clear(); - phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf()); - bool closedVolume = adjustPhi(phi, U, p_rgh); + surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf()); + + surfaceScalarField phiHbyA + ( + "phiHbyA", + fvc::interpolate(rho)*(fvc::interpolate(HbyA) & mesh.Sf()) + ); + + bool closedVolume = adjustPhi(phiHbyA, U, p_rgh); + + phiHbyA += phig; + dimensionedScalar compressibility = fvc::domainIntegrate(psi); bool compressible = (compressibility.value() > SMALL); - surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf()); - phi += phig; - // Solve pressure for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) { fvScalarMatrix p_rghEqn ( - fvm::laplacian(rhorAUf, p_rgh) == fvc::div(phi) + fvm::laplacian(rhorAUf, p_rgh) == fvc::div(phiHbyA) ); p_rghEqn.setReference @@ -37,14 +45,14 @@ if (nonOrth == nNonOrthCorr) { // Calculate the conservative fluxes - phi -= p_rghEqn.flux(); + phi = phiHbyA - p_rghEqn.flux(); // Explicitly relax pressure for momentum corrector p_rgh.relax(); // Correct the momentum source with the pressure gradient flux // calculated from the relaxed pressure - U += rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rhorAUf); + U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rhorAUf); U.correctBoundaryConditions(); } } diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H index d7e8d6e5c5..3abad39417 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H @@ -1,6 +1,8 @@ const fvMesh& mesh = fluidRegions[i]; rhoThermo& thermo = thermoFluid[i]; + thermo.validate(args.executable(), "h", "e"); + volScalarField& rho = rhoFluid[i]; volVectorField& U = UFluid[i]; surfaceScalarField& phi = phiFluid[i]; @@ -9,7 +11,6 @@ volScalarField& p = thermo.p(); const volScalarField& psi = thermo.psi(); - volScalarField& h = thermo.he(); IObasicSourceList& sources = heatSources[i]; diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/solveFluid.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/solveFluid.H index 2b6de83ca3..a17893fb83 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/solveFluid.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/solveFluid.H @@ -4,7 +4,7 @@ rho.storePrevIter(); { #include "UEqn.H" - #include "hEqn.H" + #include "EEqn.H" #include "pEqn.H" } diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousFluid/hPorousFluidEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousFluid/hPorousFluidEqn.H index 653389b221..fcc8b054a2 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousFluid/hPorousFluidEqn.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousFluid/hPorousFluidEqn.H @@ -2,7 +2,6 @@ fvScalarMatrix hPorousEqn ( fvm::div(porousPhi, porousH) - - fvm::Sp(fvc::div(porousPhi), porousH) - fvm::laplacian(turbPorous.alphaEff(), porousH) == - fvc::div(porousPhi, 0.5*magSqr(porousU), "div(phi,K)") diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/EEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/EEqn.H new file mode 100644 index 0000000000..c5680e7777 --- /dev/null +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/EEqn.H @@ -0,0 +1,27 @@ +{ + volScalarField& he = thermo.he(); + + fvScalarMatrix EEqn + ( + fvm::ddt(rho, he) + fvm::div(phi, he) + + fvc::ddt(rho, K) + fvc::div(phi, K) + + ( + he.name() == "e" + ? fvc::div(phi, volScalarField("Ep", p/rho)) + : -dpdt + ) + - fvm::laplacian(turb.alphaEff(), he) + == + rad.Sh(thermo) + + sources(rho, he) + ); + + EEqn.relax(); + EEqn.solve(mesh.solver(he.select(finalIter))); + + thermo.correct(); + rad.correct(); + + Info<< "Min/max T:" << min(thermo.T()).value() << ' ' + << max(thermo.T()).value() << endl; +} diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/hEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/hEqn.H deleted file mode 100644 index 9a94f52b75..0000000000 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/hEqn.H +++ /dev/null @@ -1,24 +0,0 @@ -{ - fvScalarMatrix hEqn - ( - fvm::ddt(rho, h) - + fvm::div(phi, h) - - fvm::laplacian(turb.alphaEff(), h) - == - dpdt - - (fvc::ddt(rho, K) + fvc::div(phi, K)) - + rad.Sh(thermo) - + sources(rho, h) - ); - - hEqn.relax(); - hEqn.solve(mesh.solver(h.select(finalIter))); - - thermo.correct(); - - rad.correct(); - - Info<< "Min/max T:" << min(thermo.T()).value() << ' ' - << max(thermo.T()).value() << endl; - -} diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H index 126b8f0961..04986d1289 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H @@ -1,6 +1,8 @@ fvMesh& mesh = fluidRegions[i]; rhoThermo& thermo = thermoFluid[i]; + thermo.validate(args.executable(), "h", "e"); + volScalarField& rho = rhoFluid[i]; volVectorField& U = UFluid[i]; surfaceScalarField& phi = phiFluid[i]; @@ -11,7 +13,6 @@ volScalarField& p = thermo.p(); const volScalarField& psi = thermo.psi(); - volScalarField& h = thermo.he(); volScalarField& p_rgh = p_rghFluid[i]; const volScalarField& gh = ghFluid[i]; diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveFluid.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveFluid.H index b36cf89d34..354c4a137c 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveFluid.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveFluid.H @@ -9,8 +9,7 @@ if (oCorr == 0) } #include "UEqn.H" - -#include "hEqn.H" +#include "EEqn.H" // --- PISO loop for (int corr=0; corrdivDevRhoReff(U) - == - rho.dimensionedInternalField()*g - + parcels.SU(U) - ); - - sources.constrain(UEqn); - - pZones.addResistance(UEqn); - - if (pimple.momentumPredictor()) - { - solve(UEqn == -fvc::grad(p) + sources(rho, U)); - } diff --git a/applications/solvers/lagrangian/LTSReactingParcelFoam/YEqn.H b/applications/solvers/lagrangian/LTSReactingParcelFoam/YEqn.H deleted file mode 100644 index d9ca1b09e0..0000000000 --- a/applications/solvers/lagrangian/LTSReactingParcelFoam/YEqn.H +++ /dev/null @@ -1,52 +0,0 @@ -tmp > mvConvection -( - fv::convectionScheme::New - ( - mesh, - fields, - phi, - mesh.divScheme("div(phi,Yi_h)") - ) -); - -combustion->correct(); -dQ = combustion->dQ(); - -if (solveSpecies) -{ - label inertIndex = -1; - volScalarField Yt(0.0*Y[0]); - - forAll(Y, i) - { - if (Y[i].name() != inertSpecie) - { - volScalarField& Yi = Y[i]; - - fvScalarMatrix YEqn - ( - fvm::ddt(rho, Yi) - + mvConvection->fvmDiv(phi, Yi) - - fvm::laplacian(turbulence->muEff(), Yi) - == - parcels.SYi(i, Yi) - + combustion->R(Yi) - + sources(rho, Yi) - ); - - sources.constrain(YEqn); - - YEqn.solve(mesh.solver("Yi")); - - Yi.max(0.0); - Yt += Yi; - } - else - { - inertIndex = i; - } - } - - Y[inertIndex] = scalar(1) - Yt; - Y[inertIndex].max(0.0); -} diff --git a/applications/solvers/lagrangian/LTSReactingParcelFoam/chemistry.H b/applications/solvers/lagrangian/LTSReactingParcelFoam/chemistry.H deleted file mode 100644 index ec92091de3..0000000000 --- a/applications/solvers/lagrangian/LTSReactingParcelFoam/chemistry.H +++ /dev/null @@ -1,28 +0,0 @@ -if (chemistry.chemistry()) -{ - Info<< "Solving chemistry" << endl; - - // update reaction rates - chemistry.calculate(); - - // turbulent time scale - if (turbulentReaction) - { - typedef DimensionedField dsfType; - - const dimensionedScalar e0("e0", sqr(dimLength)/pow3(dimTime), SMALL); - - const dsfType tk = - Cmix*sqrt(turbulence->muEff()/rho/(turbulence->epsilon() + e0)); - - const dsfType tc = chemistry.tc()().dimensionedInternalField(); - - kappa = tc/(tc + tk); - } - else - { - kappa = 1.0; - } - - chemistrySh = kappa*chemistry.Sh()(); -} diff --git a/applications/solvers/lagrangian/LTSReactingParcelFoam/createClouds.H b/applications/solvers/lagrangian/LTSReactingParcelFoam/createClouds.H deleted file mode 100644 index 954b74e069..0000000000 --- a/applications/solvers/lagrangian/LTSReactingParcelFoam/createClouds.H +++ /dev/null @@ -1,9 +0,0 @@ -Info<< "\nConstructing reacting cloud" << endl; -basicReactingMultiphaseCloud parcels -( - "reactingCloud1", - rho, - U, - g, - slgThermo -); diff --git a/applications/solvers/lagrangian/LTSReactingParcelFoam/createFields.H b/applications/solvers/lagrangian/LTSReactingParcelFoam/createFields.H deleted file mode 100644 index b78aeea132..0000000000 --- a/applications/solvers/lagrangian/LTSReactingParcelFoam/createFields.H +++ /dev/null @@ -1,121 +0,0 @@ - Info<< "Creating combustion model\n" << endl; - - autoPtr combustion - ( - combustionModels::rhoCombustionModel::New(mesh) - ); - - rhoReactionThermo& thermo = combustion->thermo(); - - SLGThermo slgThermo(mesh, thermo); - - basicMultiComponentMixture& composition = thermo.composition(); - PtrList& Y = composition.Y(); - - const word inertSpecie(thermo.lookup("inertSpecie")); - - if (!composition.contains(inertSpecie)) - { - FatalErrorIn(args.executable()) - << "Specified inert specie '" << inertSpecie << "' not found in " - << "species list. Available species:" << composition.species() - << exit(FatalError); - } - - volScalarField& p = thermo.p(); - volScalarField& hs = thermo.he(); - const volScalarField& T = thermo.T(); - const volScalarField& psi = thermo.psi(); - - volScalarField rho - ( - IOobject - ( - "rho", - runTime.timeName(), - mesh, - IOobject::NO_READ, - IOobject::AUTO_WRITE - ), - thermo.rho() - ); - - Info<< "\nReading field U\n" << endl; - volVectorField U - ( - IOobject - ( - "U", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh - ); - - #include "compressibleCreatePhi.H" - - dimensionedScalar rhoMax - ( - mesh.solutionDict().subDict("PIMPLE").lookup("rhoMax") - ); - - dimensionedScalar rhoMin - ( - mesh.solutionDict().subDict("PIMPLE").lookup("rhoMin") - ); - - Info<< "Creating turbulence model\n" << endl; - autoPtr turbulence - ( - compressible::turbulenceModel::New - ( - rho, - U, - phi, - thermo - ) - ); - - // Set the turbulence into the combustion model - combustion->setTurbulence(turbulence()); - - Info<< "Creating multi-variate interpolation scheme\n" << endl; - multivariateSurfaceInterpolationScheme::fieldTable fields; - - forAll(Y, i) - { - fields.add(Y[i]); - } - fields.add(hs); - - volScalarField dQ - ( - IOobject - ( - "dQ", - runTime.timeName(), - mesh, - IOobject::NO_READ, - IOobject::AUTO_WRITE - ), - mesh, - dimensionedScalar("dQ", dimEnergy/dimTime, 0.0) - ); - - - volScalarField rDeltaT - ( - IOobject - ( - "rDeltaT", - runTime.timeName(), - mesh, - IOobject::READ_IF_PRESENT, - IOobject::AUTO_WRITE - ), - mesh, - dimensionedScalar("one", dimless/dimTime, 1), - zeroGradientFvPatchScalarField::typeName - ); diff --git a/applications/solvers/lagrangian/LTSReactingParcelFoam/hsEqn.H b/applications/solvers/lagrangian/LTSReactingParcelFoam/hsEqn.H deleted file mode 100644 index 15698b6afa..0000000000 --- a/applications/solvers/lagrangian/LTSReactingParcelFoam/hsEqn.H +++ /dev/null @@ -1,25 +0,0 @@ -{ - fvScalarMatrix hsEqn - ( - fvm::ddt(rho, hs) - + mvConvection->fvmDiv(phi, hs) - - fvm::laplacian(turbulence->alphaEff(), hs) - == - - fvc::div(phi, 0.5*magSqr(U), "div(phi,K)") - + parcels.Sh(hs) - + radiation->Sh(thermo) - + combustion->Sh() - + sources(rho, hs) - ); - - sources.constrain(hsEqn); - - hsEqn.solve(); - - thermo.correct(); - - radiation->correct(); - - Info<< "T gas min/max = " << min(T).value() << ", " - << max(T).value() << endl; -} diff --git a/applications/solvers/lagrangian/LTSReactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/LTSReactingParcelFoam/pEqn.H deleted file mode 100644 index e6df2ca5cf..0000000000 --- a/applications/solvers/lagrangian/LTSReactingParcelFoam/pEqn.H +++ /dev/null @@ -1,66 +0,0 @@ -{ - rho = thermo.rho(); - - // Thermodynamic density needs to be updated by psi*d(p) after the - // pressure solution - done in 2 parts. Part 1: - thermo.rho() -= psi*p; - - volScalarField rAU(1.0/UEqn.A()); - volVectorField HbyA("HbyA", U); - HbyA = rAU*(UEqn == sources(rho, U))().H(); - - surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(HbyA) & mesh.Sf()); - if (pZones.size() == 0) - { - // ddtPhiCorr only used without porosity - phiHbyA += fvc::ddtPhiCorr(rAU, rho, U, phi); - } - - phiHbyA *= fvc::interpolate(rho); - - - fvScalarMatrix pDDtEqn - ( - fvc::ddt(rho) + psi*correction(fvm::ddt(p)) - + fvc::div(phiHbyA) - == - parcels.Srho() - + sources(psi, p, rho.name()) - ); - - while (pimple.correctNonOrthogonal()) - { - fvScalarMatrix pEqn - ( - pDDtEqn - - fvm::laplacian(rho*rAU, p) - ); - - sources.constrain(pDDtEqn, rho.name()); - - pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter()))); - - if (pimple.finalNonOrthogonalIter()) - { - phi = phiHbyA + pEqn.flux(); - } - } - - p.relax(); - - // Second part of thermodynamic density update - thermo.rho() += psi*p; - - #include "rhoEqn.H" // NOTE: flux and time scales now inconsistent - #include "compressibleContinuityErrs.H" - - U = HbyA - rAU*fvc::grad(p); - U.correctBoundaryConditions(); - sources.correct(U); - - rho = thermo.rho(); - rho = max(rho, rhoMin); - rho = min(rho, rhoMax); - - Info<< "p min/max = " << min(p).value() << ", " << max(p).value() << endl; -} diff --git a/applications/solvers/lagrangian/LTSReactingParcelFoam/readAdditionalSolutionControls.H b/applications/solvers/lagrangian/LTSReactingParcelFoam/readAdditionalSolutionControls.H deleted file mode 100644 index 8136af1ad3..0000000000 --- a/applications/solvers/lagrangian/LTSReactingParcelFoam/readAdditionalSolutionControls.H +++ /dev/null @@ -1,9 +0,0 @@ -dictionary additional = mesh.solutionDict().subDict("additional"); - -// pressure work term for enthalpy equation -bool pressureWork = additional.lookupOrDefault("pressureWork", true); -bool pressureWorkTimeDerivative = - additional.lookupOrDefault("pressureWorkTimeDerivative", true); - -// flag to activate solve transport for each specie (Y vector) -bool solveSpecies = additional.lookupOrDefault("solveSpecies", true); diff --git a/applications/solvers/lagrangian/coalChemistryFoam/EEqn.H b/applications/solvers/lagrangian/coalChemistryFoam/EEqn.H new file mode 100644 index 0000000000..006d2231f4 --- /dev/null +++ b/applications/solvers/lagrangian/coalChemistryFoam/EEqn.H @@ -0,0 +1,31 @@ +{ + volScalarField& he = thermo.he(); + + fvScalarMatrix EEqn + ( + fvm::ddt(rho, he) + mvConvection->fvmDiv(phi, he) + + fvc::ddt(rho, K) + fvc::div(phi, K) + + ( + he.name() == "e" + ? fvc::div(phi, volScalarField("Ep", p/rho)) + : -dpdt + ) + - fvm::laplacian(turbulence->alphaEff(), he) + == + combustion->Sh() + + coalParcels.Sh(he) + + limestoneParcels.Sh(he) + + radiation->Sh(thermo) + + sources(rho, he) + ); + + EEqn.relax(); + sources.constrain(EEqn); + EEqn.solve(); + + thermo.correct(); + radiation->correct(); + + Info<< "T gas min/max = " << min(T).value() << ", " + << max(T).value() << endl; +} diff --git a/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C b/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C index 9cdd2b7295..6832d01458 100644 --- a/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C +++ b/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C @@ -92,7 +92,7 @@ int main(int argc, char *argv[]) { #include "UEqn.H" #include "YEqn.H" - #include "hsEqn.H" + #include "EEqn.H" // --- Pressure corrector loop while (pimple.correct()) diff --git a/applications/solvers/lagrangian/coalChemistryFoam/createFields.H b/applications/solvers/lagrangian/coalChemistryFoam/createFields.H index ef241d2103..d76805c25a 100644 --- a/applications/solvers/lagrangian/coalChemistryFoam/createFields.H +++ b/applications/solvers/lagrangian/coalChemistryFoam/createFields.H @@ -6,6 +6,7 @@ ); psiReactionThermo& thermo = combustion->thermo(); + thermo.validate(args.executable(), "h", "e"); SLGThermo slgThermo(mesh, thermo); @@ -23,7 +24,6 @@ } volScalarField& p = thermo.p(); - volScalarField& hs = thermo.he(); const volScalarField& T = thermo.T(); const volScalarField& psi = thermo.psi(); @@ -33,7 +33,7 @@ { fields.add(Y[i]); } - fields.add(hs); + fields.add(thermo.he()); volScalarField rho ( diff --git a/applications/solvers/lagrangian/coalChemistryFoam/hsEqn.H b/applications/solvers/lagrangian/coalChemistryFoam/hsEqn.H deleted file mode 100644 index 37d29778c0..0000000000 --- a/applications/solvers/lagrangian/coalChemistryFoam/hsEqn.H +++ /dev/null @@ -1,29 +0,0 @@ -{ - fvScalarMatrix hsEqn - ( - fvm::ddt(rho, hs) - + mvConvection->fvmDiv(phi, hs) - - fvm::laplacian(turbulence->alphaEff(), hs) - == - dpdt - - (fvc::ddt(rho, K) + fvc::div(phi, K)) - + combustion->Sh() - + coalParcels.Sh(hs) - + limestoneParcels.Sh(hs) - + radiation->Sh(thermo) - + sources(rho, hs) - ); - - hsEqn.relax(); - - sources.constrain(hsEqn); - - hsEqn.solve(); - - thermo.correct(); - - radiation->correct(); - - Info<< "T gas min/max = " << min(T).value() << ", " - << max(T).value() << endl; -} diff --git a/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/Allwmake b/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/Allwmake new file mode 100755 index 0000000000..845cebd777 --- /dev/null +++ b/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/Allwmake @@ -0,0 +1,8 @@ +#!/bin/sh +cd ${0%/*} || exit 1 # run from this directory +set -x + +wmake +wmake icoUncoupledKinematicParcelDyMFoam + +# ----------------------------------------------------------------- end-of-file diff --git a/applications/solvers/lagrangian/icoUncoupledKinematicParcelDyMFoam/Make/files b/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelDyMFoam/Make/files similarity index 100% rename from applications/solvers/lagrangian/icoUncoupledKinematicParcelDyMFoam/Make/files rename to applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelDyMFoam/Make/files diff --git a/applications/solvers/lagrangian/icoUncoupledKinematicParcelDyMFoam/Make/options b/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelDyMFoam/Make/options similarity index 96% rename from applications/solvers/lagrangian/icoUncoupledKinematicParcelDyMFoam/Make/options rename to applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelDyMFoam/Make/options index 4535f635dc..942c217fa4 100644 --- a/applications/solvers/lagrangian/icoUncoupledKinematicParcelDyMFoam/Make/options +++ b/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelDyMFoam/Make/options @@ -1,5 +1,5 @@ EXE_INC = \ - -I../icoUncoupledKinematicParcelFoam \ + -I.. \ -I$(LIB_SRC)/lagrangian/basic/lnInclude \ -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \ diff --git a/applications/solvers/lagrangian/icoUncoupledKinematicParcelDyMFoam/icoUncoupledKinematicParcelDyMFoam.C b/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelDyMFoam/icoUncoupledKinematicParcelDyMFoam.C similarity index 100% rename from applications/solvers/lagrangian/icoUncoupledKinematicParcelDyMFoam/icoUncoupledKinematicParcelDyMFoam.C rename to applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelDyMFoam/icoUncoupledKinematicParcelDyMFoam.C diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/Make/files b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/Make/files deleted file mode 100644 index dfee7d7d19..0000000000 --- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/Make/files +++ /dev/null @@ -1,3 +0,0 @@ -porousExplicitSourceReactingParcelFoam.C - -EXE = $(FOAM_APPBIN)/porousExplicitSourceReactingParcelFoam diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/Make/options b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/Make/options deleted file mode 100644 index 1f776f0412..0000000000 --- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/Make/options +++ /dev/null @@ -1,54 +0,0 @@ -EXE_INC = \ - -I$(LIB_SRC)/finiteVolume/lnInclude \ - -I${LIB_SRC}/meshTools/lnInclude \ - -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \ - -I$(LIB_SRC)/lagrangian/basic/lnInclude \ - -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \ - -I$(LIB_SRC)/lagrangian/coalCombustion/lnInclude \ - -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \ - -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \ - -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ - -I$(LIB_SRC)/thermophysicalModels/properties/liquidProperties/lnInclude \ - -I$(LIB_SRC)/thermophysicalModels/properties/liquidMixtureProperties/lnInclude \ - -I$(LIB_SRC)/thermophysicalModels/properties/solidProperties/lnInclude \ - -I$(LIB_SRC)/thermophysicalModels/properties/solidMixtureProperties/lnInclude \ - -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \ - -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \ - -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \ - -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \ - -I$(LIB_SRC)/thermophysicalModels/radiationModels/lnInclude \ - -I$(LIB_SRC)/ODE/lnInclude \ - -I$(LIB_SRC)/regionModels/regionModel/lnInclude \ - -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \ - -I$(LIB_SRC)/combustionModels/lnInclude \ - -I$(LIB_SRC)/fieldSources/lnInclude \ - -I$(LIB_SRC)/sampling/lnInclude \ - -I$(FOAM_SOLVERS)/combustion/reactingFoam - - -EXE_LIBS = \ - -lfiniteVolume \ - -lmeshTools \ - -lcompressibleTurbulenceModel \ - -lcompressibleRASModels \ - -lcompressibleLESModels \ - -llagrangian \ - -llagrangianIntermediate \ - -lspecie \ - -lfluidThermophysicalModels \ - -lliquidProperties \ - -lliquidMixtureProperties \ - -lsolidProperties \ - -lsolidMixtureProperties \ - -lthermophysicalFunctions \ - -lreactionThermophysicalModels \ - -lSLGThermo \ - -lchemistryModel \ - -lradiationModels \ - -lODE \ - -lregionModels \ - -lsurfaceFilmModels \ - -lcombustionModels \ - -lfieldSources \ - -lsampling - diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/YEqn.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/YEqn.H deleted file mode 100644 index 60f6e1a55d..0000000000 --- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/YEqn.H +++ /dev/null @@ -1,55 +0,0 @@ - -tmp > mvConvection -( - fv::convectionScheme::New - ( - mesh, - fields, - phi, - mesh.divScheme("div(phi,Yi_h)") - ) -); - -combustion->correct(); -dQ = combustion->dQ(); - -if (solveSpecies) -{ - label inertIndex = -1; - volScalarField Yt(0.0*Y[0]); - - forAll(Y, i) - { - if (Y[i].name() != inertSpecie) - { - volScalarField& Yi = Y[i]; - - fvScalarMatrix YiEqn - ( - fvm::ddt(rho, Yi) - + mvConvection->fvmDiv(phi, Yi) - - fvm::laplacian(turbulence->muEff(), Yi) - == - parcels.SYi(i, Yi) - + combustion->R(Yi) - + sources(rho, Yi) - ); - - YiEqn.relax(); - - sources.constrain(YiEqn); - - YiEqn.solve(mesh.solver("Yi")); - - Yi.max(0.0); - Yt += Yi; - } - else - { - inertIndex = i; - } - } - - Y[inertIndex] = scalar(1) - Yt; - Y[inertIndex].max(0.0); -} diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createClouds.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createClouds.H deleted file mode 100644 index 954b74e069..0000000000 --- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createClouds.H +++ /dev/null @@ -1,9 +0,0 @@ -Info<< "\nConstructing reacting cloud" << endl; -basicReactingMultiphaseCloud parcels -( - "reactingCloud1", - rho, - U, - g, - slgThermo -); diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createExplicitSources.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createExplicitSources.H deleted file mode 100644 index b2e5775e58..0000000000 --- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createExplicitSources.H +++ /dev/null @@ -1,2 +0,0 @@ -Info<< "Creating sources\n" << endl; -IObasicSourceList sources(mesh); diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createPorousZones.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createPorousZones.H deleted file mode 100644 index 90506856d2..0000000000 --- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createPorousZones.H +++ /dev/null @@ -1,3 +0,0 @@ - Info<< "Creating porous zones" << nl << endl; - - porousZones pZones(mesh); diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/hsEqn.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/hsEqn.H deleted file mode 100644 index 15c598cd59..0000000000 --- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/hsEqn.H +++ /dev/null @@ -1,32 +0,0 @@ -{ - fvScalarMatrix hsEqn - ( - fvm::ddt(rho, hs) - + mvConvection->fvmDiv(phi, hs) - - fvm::laplacian(turbulence->alphaEff(), hs) - == - - (fvc::ddt(rho, K) + fvc::div(phi, K)) - + parcels.Sh(hs) - + radiation->Sh(thermo) - + combustion->Sh() - + sources(rho, hs) - ); - - if (pressureWorkTimeDerivative) - { - hsEqn -= dpdt; - } - - hsEqn.relax(); - - sources.constrain(hsEqn); - - hsEqn.solve(); - - thermo.correct(); - - radiation->correct(); - - Info<< "T gas min/max = " << min(T).value() << ", " - << max(T).value() << endl; -} diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H deleted file mode 100644 index 4401afc9fc..0000000000 --- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H +++ /dev/null @@ -1,66 +0,0 @@ -{ - rho = thermo.rho(); - - // Thermodynamic density needs to be updated by psi*d(p) after the - // pressure solution - done in 2 parts. Part 1: - thermo.rho() -= psi*p; - - volScalarField rAU(1.0/UEqn.A()); - volVectorField HbyA("HbyA", U); - HbyA = rAU*(UEqn == sources(rho, U))().H(); - - surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(HbyA) & mesh.Sf()); - if (pZones.size() == 0) - { - // ddtPhiCorr only used without porosity - phiHbyA += fvc::ddtPhiCorr(rAU, rho, U, phi); - } - - phiHbyA *= fvc::interpolate(rho); - - - fvScalarMatrix pDDtEqn - ( - fvc::ddt(rho) + psi*correction(fvm::ddt(p)) - + fvc::div(phiHbyA) - == - parcels.Srho() - + sources(psi, p, rho.name()) - ); - - - while (pimple.correctNonOrthogonal()) - { - fvScalarMatrix pEqn - ( - pDDtEqn - - fvm::laplacian(rho*rAU, p) - ); - - sources.constrain(pDDtEqn, rho.name()); - - pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter()))); - - if (pimple.finalNonOrthogonalIter()) - { - phi = phiHbyA + pEqn.flux(); - } - } - - // Second part of thermodynamic density update - thermo.rho() += psi*p; - - #include "rhoEqn.H" - #include "compressibleContinuityErrs.H" - - U = HbyA - rAU*fvc::grad(p); - U.correctBoundaryConditions(); - sources.correct(U); - - K = 0.5*magSqr(U); - - if (thermo.dpdt()) - { - dpdt = fvc::ddt(p); - } -} diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/porousExplicitSourceReactingParcelFoam.C b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/porousExplicitSourceReactingParcelFoam.C deleted file mode 100644 index 61533fcd6b..0000000000 --- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/porousExplicitSourceReactingParcelFoam.C +++ /dev/null @@ -1,126 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation - \\/ 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 3 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, see . - -Application - porousExplicitSourceReactingParcelFoam - -Description - Transient PIMPLE solver for compressible, laminar or turbulent flow with - reacting multiphase Lagrangian parcels for porous media, including explicit - sources for mass, momentum and energy - - The solver includes: - - reacting multiphase parcel cloud - - porous media - - mass, momentum and energy sources - - Note: ddtPhiCorr not used here when porous zones are active - - not well defined for porous calculations - -\*---------------------------------------------------------------------------*/ - -#include "fvCFD.H" -#include "turbulenceModel.H" -#include "basicReactingMultiphaseCloud.H" -#include "rhoCombustionModel.H" -#include "radiationModel.H" -#include "porousZones.H" -#include "IObasicSourceList.H" -#include "SLGThermo.H" -#include "pimpleControl.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -int main(int argc, char *argv[]) -{ - #include "setRootCase.H" - - #include "createTime.H" - #include "createMesh.H" - #include "readGravitationalAcceleration.H" - #include "createFields.H" - #include "createRadiationModel.H" - #include "createClouds.H" - #include "createExplicitSources.H" - #include "createPorousZones.H" - #include "initContinuityErrs.H" - #include "readTimeControls.H" - #include "compressibleCourantNo.H" - #include "setInitialDeltaT.H" - - pimpleControl pimple(mesh); - - // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - - Info<< "\nStarting time loop\n" << endl; - - while (runTime.run()) - { - #include "readTimeControls.H" - #include "readAdditionalSolutionControls.H" - #include "compressibleCourantNo.H" - #include "setDeltaT.H" - - runTime++; - - Info<< "Time = " << runTime.timeName() << nl << endl; - - parcels.evolve(); - - #include "rhoEqn.H" - - // --- Pressure-velocity PIMPLE corrector loop - while (pimple.loop()) - { - #include "UEqn.H" - #include "YEqn.H" - #include "hsEqn.H" - - // --- Pressure corrector loop - while (pimple.correct()) - { - #include "pEqn.H" - } - - if (pimple.turbCorr()) - { - turbulence->correct(); - } - } - - rho = thermo.rho(); - - runTime.write(); - - Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" - << " ClockTime = " << runTime.elapsedClockTime() << " s" - << nl << endl; - } - - Info<< "End\n" << endl; - - return(0); -} - - -// ************************************************************************* // diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/EEqn.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/EEqn.H new file mode 100644 index 0000000000..cdc475d42f --- /dev/null +++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/EEqn.H @@ -0,0 +1,29 @@ +{ + volScalarField& he = thermo.he(); + + fvScalarMatrix EEqn + ( + fvm::ddt(rho, he) + mvConvection->fvmDiv(phi, he) + + fvc::ddt(rho, K) + fvc::div(phi, K) + + ( + he.name() == "e" + ? fvc::div(phi, volScalarField("Ep", p/rho)) + : -dpdt + ) + - fvm::laplacian(turbulence->alphaEff(), he) + == + parcels.Sh(he) + + surfaceFilm.Sh() + + radiation->Sh(thermo) + + combustion->Sh() + ); + + EEqn.relax(); + EEqn.solve(); + + thermo.correct(); + radiation->correct(); + + Info<< "T gas min/max = " << min(T).value() << ", " + << max(T).value() << endl; +} diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/createFields.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/createFields.H index d5cd8e6aee..057978c8f8 100644 --- a/applications/solvers/lagrangian/reactingParcelFilmFoam/createFields.H +++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/createFields.H @@ -6,6 +6,7 @@ ); psiReactionThermo& thermo = combustion->thermo(); + thermo.validate(args.executable(), "h", "e"); SLGThermo slgThermo(mesh, thermo); @@ -29,7 +30,6 @@ ); volScalarField& p = thermo.p(); - volScalarField& hs = thermo.he(); const volScalarField& T = thermo.T(); const volScalarField& psi = thermo.psi(); @@ -107,7 +107,7 @@ { fields.add(Y[i]); } - fields.add(hs); + fields.add(thermo.he()); IOdictionary additionalControlsDict ( diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/hsEqn.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/hsEqn.H deleted file mode 100644 index dba268581e..0000000000 --- a/applications/solvers/lagrangian/reactingParcelFilmFoam/hsEqn.H +++ /dev/null @@ -1,25 +0,0 @@ -{ - fvScalarMatrix hsEqn - ( - fvm::ddt(rho, hs) - + mvConvection->fvmDiv(phi, hs) - - fvm::laplacian(turbulence->alphaEff(), hs) - == - dpdt - - (fvc::ddt(rho, K) + fvc::div(phi, K)) - + parcels.Sh(hs) - + surfaceFilm.Sh() - + radiation->Sh(thermo) - + combustion->Sh() - ); - - hsEqn.relax(); - - hsEqn.solve(); - - thermo.correct(); - - radiation->correct(); - - Info<< "min/max(T) = " << min(T).value() << ", " << max(T).value() << endl; -} diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C b/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C index 9badba12b7..06df5b1748 100644 --- a/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C +++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C @@ -87,7 +87,7 @@ int main(int argc, char *argv[]) { #include "UEqn.H" #include "YEqn.H" - #include "hsEqn.H" + #include "EEqn.H" // --- Pressure corrector loop while (pimple.correct()) diff --git a/applications/solvers/lagrangian/reactingParcelFoam/Allwmake b/applications/solvers/lagrangian/reactingParcelFoam/Allwmake new file mode 100755 index 0000000000..ecd3f26015 --- /dev/null +++ b/applications/solvers/lagrangian/reactingParcelFoam/Allwmake @@ -0,0 +1,8 @@ +#!/bin/sh +cd ${0%/*} || exit 1 # run from this directory +set -x + +wmake +wmake LTSReactingParcelFoam + +# ----------------------------------------------------------------- end-of-file diff --git a/applications/solvers/lagrangian/reactingParcelFoam/EEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/EEqn.H new file mode 100644 index 0000000000..14996cfdf5 --- /dev/null +++ b/applications/solvers/lagrangian/reactingParcelFoam/EEqn.H @@ -0,0 +1,30 @@ +{ + volScalarField& he = thermo.he(); + + fvScalarMatrix EEqn + ( + fvm::ddt(rho, he) + mvConvection->fvmDiv(phi, he) + + fvc::ddt(rho, K) + fvc::div(phi, K) + + ( + he.name() == "e" + ? fvc::div(phi, volScalarField("Ep", p/rho)) + : -dpdt + ) + - fvm::laplacian(turbulence->alphaEff(), he) + == + parcels.Sh(he) + + radiation->Sh(thermo) + + combustion->Sh() + + sources(rho, he) + ); + + EEqn.relax(); + sources.constrain(EEqn); + EEqn.solve(); + + thermo.correct(); + radiation->correct(); + + Info<< "T gas min/max = " << min(T).value() << ", " + << max(T).value() << endl; +} diff --git a/applications/solvers/lagrangian/LTSReactingParcelFoam/LTSReactingParcelFoam.C b/applications/solvers/lagrangian/reactingParcelFoam/LTSReactingParcelFoam/LTSReactingParcelFoam.C similarity index 98% rename from applications/solvers/lagrangian/LTSReactingParcelFoam/LTSReactingParcelFoam.C rename to applications/solvers/lagrangian/reactingParcelFoam/LTSReactingParcelFoam/LTSReactingParcelFoam.C index 6f41f328a3..0841d458b4 100644 --- a/applications/solvers/lagrangian/LTSReactingParcelFoam/LTSReactingParcelFoam.C +++ b/applications/solvers/lagrangian/reactingParcelFoam/LTSReactingParcelFoam/LTSReactingParcelFoam.C @@ -61,6 +61,7 @@ int main(int argc, char *argv[]) #include "readTimeControls.H" #include "readAdditionalSolutionControls.H" #include "createFields.H" + #include "createRDeltaT.H" #include "createRadiationModel.H" #include "createClouds.H" #include "createExplicitSources.H" @@ -93,7 +94,7 @@ int main(int argc, char *argv[]) #include "UEqn.H" #include "YEqn.H" - #include "hsEqn.H" + #include "EEqn.H" // --- Pressure corrector loop while (pimple.correct()) diff --git a/applications/solvers/lagrangian/LTSReactingParcelFoam/Make/files b/applications/solvers/lagrangian/reactingParcelFoam/LTSReactingParcelFoam/Make/files similarity index 100% rename from applications/solvers/lagrangian/LTSReactingParcelFoam/Make/files rename to applications/solvers/lagrangian/reactingParcelFoam/LTSReactingParcelFoam/Make/files diff --git a/applications/solvers/lagrangian/LTSReactingParcelFoam/Make/options b/applications/solvers/lagrangian/reactingParcelFoam/LTSReactingParcelFoam/Make/options similarity index 99% rename from applications/solvers/lagrangian/LTSReactingParcelFoam/Make/options rename to applications/solvers/lagrangian/reactingParcelFoam/LTSReactingParcelFoam/Make/options index f87b446266..fad2ee86fa 100644 --- a/applications/solvers/lagrangian/LTSReactingParcelFoam/Make/options +++ b/applications/solvers/lagrangian/reactingParcelFoam/LTSReactingParcelFoam/Make/options @@ -1,4 +1,5 @@ EXE_INC = \ + -I.. \ -I$(LIB_SRC)/finiteVolume/lnInclude \ -I${LIB_SRC}/meshTools/lnInclude \ -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \ diff --git a/applications/solvers/lagrangian/reactingParcelFoam/LTSReactingParcelFoam/createRDeltaT.H b/applications/solvers/lagrangian/reactingParcelFoam/LTSReactingParcelFoam/createRDeltaT.H new file mode 100644 index 0000000000..c445eaee71 --- /dev/null +++ b/applications/solvers/lagrangian/reactingParcelFoam/LTSReactingParcelFoam/createRDeltaT.H @@ -0,0 +1,14 @@ + volScalarField rDeltaT + ( + IOobject + ( + "rDeltaT", + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar("one", dimless/dimTime, 1), + zeroGradientFvPatchScalarField::typeName + ); diff --git a/applications/solvers/lagrangian/LTSReactingParcelFoam/readTimeControls.H b/applications/solvers/lagrangian/reactingParcelFoam/LTSReactingParcelFoam/readTimeControls.H similarity index 95% rename from applications/solvers/lagrangian/LTSReactingParcelFoam/readTimeControls.H rename to applications/solvers/lagrangian/reactingParcelFoam/LTSReactingParcelFoam/readTimeControls.H index 7b31fe9a84..73083f5fdd 100644 --- a/applications/solvers/lagrangian/LTSReactingParcelFoam/readTimeControls.H +++ b/applications/solvers/lagrangian/reactingParcelFoam/LTSReactingParcelFoam/readTimeControls.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License diff --git a/applications/solvers/lagrangian/LTSReactingParcelFoam/timeScales.H b/applications/solvers/lagrangian/reactingParcelFoam/LTSReactingParcelFoam/timeScales.H similarity index 97% rename from applications/solvers/lagrangian/LTSReactingParcelFoam/timeScales.H rename to applications/solvers/lagrangian/reactingParcelFoam/LTSReactingParcelFoam/timeScales.H index 24277a9f45..c1e9c20290 100644 --- a/applications/solvers/lagrangian/LTSReactingParcelFoam/timeScales.H +++ b/applications/solvers/lagrangian/reactingParcelFoam/LTSReactingParcelFoam/timeScales.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License diff --git a/applications/solvers/lagrangian/reactingParcelFoam/Make/options b/applications/solvers/lagrangian/reactingParcelFoam/Make/options index 6438a34c29..572cd1fb2c 100644 --- a/applications/solvers/lagrangian/reactingParcelFoam/Make/options +++ b/applications/solvers/lagrangian/reactingParcelFoam/Make/options @@ -1,10 +1,10 @@ EXE_INC = \ -I$(LIB_SRC)/finiteVolume/lnInclude \ -I${LIB_SRC}/meshTools/lnInclude \ - -I${LIB_SRC}/sampling/lnInclude \ -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \ -I$(LIB_SRC)/lagrangian/basic/lnInclude \ -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \ + -I$(LIB_SRC)/lagrangian/coalCombustion/lnInclude \ -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ @@ -22,13 +22,13 @@ EXE_INC = \ -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \ -I$(LIB_SRC)/combustionModels/lnInclude \ -I$(LIB_SRC)/fieldSources/lnInclude \ + -I$(LIB_SRC)/sampling/lnInclude \ -I$(FOAM_SOLVERS)/combustion/reactingFoam EXE_LIBS = \ -lfiniteVolume \ -lmeshTools \ - -lsampling \ -lcompressibleTurbulenceModel \ -lcompressibleRASModels \ -lcompressibleLESModels \ @@ -48,5 +48,6 @@ EXE_LIBS = \ -lODE \ -lregionModels \ -lsurfaceFilmModels \ + -lcombustionModels \ -lfieldSources \ - -lcombustionModels + -lsampling diff --git a/applications/solvers/lagrangian/reactingParcelFoam/UEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/UEqn.H index 2ddaa2c301..090922e572 100644 --- a/applications/solvers/lagrangian/reactingParcelFoam/UEqn.H +++ b/applications/solvers/lagrangian/reactingParcelFoam/UEqn.H @@ -1,5 +1,6 @@ fvVectorMatrix UEqn ( + //pZones.ddt(rho, U) fvm::ddt(rho, U) + fvm::div(phi, U) + turbulence->divDevRhoReff(U) @@ -12,6 +13,8 @@ sources.constrain(UEqn); + pZones.addResistance(UEqn); + if (pimple.momentumPredictor()) { solve(UEqn == -fvc::grad(p) + sources(rho, U)); diff --git a/applications/solvers/lagrangian/reactingParcelFoam/YEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/YEqn.H index 6cc645bde5..8ef6c610dd 100644 --- a/applications/solvers/lagrangian/reactingParcelFoam/YEqn.H +++ b/applications/solvers/lagrangian/reactingParcelFoam/YEqn.H @@ -5,14 +5,15 @@ tmp > mvConvection mesh, fields, phi, - mesh.divScheme("div(phi,Yi_hs)") + mesh.divScheme("div(phi,Yi_h)") ) ); +combustion->correct(); +dQ = combustion->dQ(); +if (solveSpecies) { - combustion->correct(); - dQ = combustion->dQ(); label inertIndex = -1; volScalarField Yt(0.0*Y[0]); @@ -22,22 +23,22 @@ tmp > mvConvection { volScalarField& Yi = Y[i]; - fvScalarMatrix YiEqn + fvScalarMatrix YEqn ( fvm::ddt(rho, Yi) + mvConvection->fvmDiv(phi, Yi) - fvm::laplacian(turbulence->muEff(), Yi) - == + == parcels.SYi(i, Yi) + combustion->R(Yi) + sources(rho, Yi) ); - YiEqn.relax(); + YEqn.relax(); - sources.constrain(YiEqn); + sources.constrain(YEqn); - YiEqn.solve(mesh.solver("Yi")); + YEqn.solve(mesh.solver("Yi")); Yi.max(0.0); Yt += Yi; diff --git a/applications/solvers/lagrangian/reactingParcelFoam/createClouds.H b/applications/solvers/lagrangian/reactingParcelFoam/createClouds.H index c568be12a1..954b74e069 100644 --- a/applications/solvers/lagrangian/reactingParcelFoam/createClouds.H +++ b/applications/solvers/lagrangian/reactingParcelFoam/createClouds.H @@ -1,5 +1,5 @@ Info<< "\nConstructing reacting cloud" << endl; -basicReactingCloud parcels +basicReactingMultiphaseCloud parcels ( "reactingCloud1", rho, diff --git a/applications/solvers/lagrangian/LTSReactingParcelFoam/createExplicitSources.H b/applications/solvers/lagrangian/reactingParcelFoam/createExplicitSources.H similarity index 100% rename from applications/solvers/lagrangian/LTSReactingParcelFoam/createExplicitSources.H rename to applications/solvers/lagrangian/reactingParcelFoam/createExplicitSources.H diff --git a/applications/solvers/lagrangian/reactingParcelFoam/createFields.H b/applications/solvers/lagrangian/reactingParcelFoam/createFields.H index ead3dac75f..d019704553 100644 --- a/applications/solvers/lagrangian/reactingParcelFoam/createFields.H +++ b/applications/solvers/lagrangian/reactingParcelFoam/createFields.H @@ -1,11 +1,12 @@ Info<< "Creating combustion model\n" << endl; - autoPtr combustion + autoPtr combustion ( - combustionModels::psiCombustionModel::New(mesh) + combustionModels::rhoCombustionModel::New(mesh) ); - psiReactionThermo& thermo = combustion->thermo(); + rhoReactionThermo& thermo = combustion->thermo(); + thermo.validate(args.executable(), "h", "e"); SLGThermo slgThermo(mesh, thermo); @@ -23,7 +24,6 @@ } volScalarField& p = thermo.p(); - volScalarField& hs = thermo.he(); const volScalarField& T = thermo.T(); const volScalarField& psi = thermo.psi(); @@ -56,6 +56,22 @@ #include "compressibleCreatePhi.H" + dimensionedScalar rhoMax + ( + "rhoMax", + dimDensity, + mesh.solutionDict().subDict("PIMPLE") + .lookupOrDefault("rhoMax", GREAT) + ); + + dimensionedScalar rhoMin + ( + "rhoMin", + dimDensity, + mesh.solutionDict().subDict("PIMPLE") + .lookupOrDefault("rhoMin", -GREAT) + ); + Info<< "Creating turbulence model\n" << endl; autoPtr turbulence ( @@ -86,13 +102,15 @@ Info<< "Creating field kinetic energy K\n" << endl; volScalarField K("K", 0.5*magSqr(U)); + + Info<< "Creating multi-variate interpolation scheme\n" << endl; multivariateSurfaceInterpolationScheme::fieldTable fields; forAll(Y, i) { fields.add(Y[i]); } - fields.add(hs); + fields.add(thermo.he()); volScalarField dQ ( @@ -107,6 +125,3 @@ mesh, dimensionedScalar("dQ", dimEnergy/dimTime, 0.0) ); - - Info<< "Creating sources\n" << endl; - IObasicSourceList sources(mesh); diff --git a/applications/solvers/lagrangian/LTSReactingParcelFoam/createPorousZones.H b/applications/solvers/lagrangian/reactingParcelFoam/createPorousZones.H similarity index 100% rename from applications/solvers/lagrangian/LTSReactingParcelFoam/createPorousZones.H rename to applications/solvers/lagrangian/reactingParcelFoam/createPorousZones.H diff --git a/applications/solvers/lagrangian/reactingParcelFoam/hsEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/hsEqn.H deleted file mode 100644 index c7c1ec0957..0000000000 --- a/applications/solvers/lagrangian/reactingParcelFoam/hsEqn.H +++ /dev/null @@ -1,28 +0,0 @@ -{ - fvScalarMatrix hsEqn - ( - fvm::ddt(rho, hs) - + mvConvection->fvmDiv(phi, hs) - - fvm::laplacian(turbulence->alphaEff(), hs) - == - dpdt - - (fvc::ddt(rho, K) + fvc::div(phi, K)) - + parcels.Sh(hs) - + radiation->Sh(thermo) - + combustion->Sh() - + sources(rho, hs) - ); - - hsEqn.relax(); - - sources.constrain(hsEqn); - - hsEqn.solve(); - - thermo.correct(); - - radiation->correct(); - - Info<< "T gas min/max = " << min(T).value() << ", " - << max(T).value() << endl; -} diff --git a/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H index d476a6bda8..5b90227323 100644 --- a/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H +++ b/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H @@ -1,65 +1,43 @@ -rho = thermo.rho(); - -volScalarField rAU(1.0/UEqn.A()); -volVectorField HbyA("HbyA", U); -HbyA = rAU*(UEqn == sources(rho, U))().H(); - -if (pimple.transonic()) { - surfaceScalarField phid - ( - "phid", - fvc::interpolate(psi) - *( - (fvc::interpolate(HbyA) & mesh.Sf()) - + fvc::ddtPhiCorr(rAU, rho, U, phi) - ) - ); + rho = thermo.rho(); - while (pimple.correctNonOrthogonal()) + // Thermodynamic density needs to be updated by psi*d(p) after the + // pressure solution - done in 2 parts. Part 1: + thermo.rho() -= psi*p; + + volScalarField rAU(1.0/UEqn.A()); + volVectorField HbyA("HbyA", U); + HbyA = rAU*(UEqn == sources(rho, U))().H(); + + surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(HbyA) & mesh.Sf()); + if (pZones.size() == 0) { - fvScalarMatrix pEqn - ( - fvm::ddt(psi, p) - + fvm::div(phid, p) - - fvm::laplacian(rho*rAU, p) - == - parcels.Srho() - + sources(psi, p, rho.name()) - ); - - pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter()))); - - if (pimple.finalNonOrthogonalIter()) - { - phi == pEqn.flux(); - } + // ddtPhiCorr only used without porosity + phiHbyA += fvc::ddtPhiCorr(rAU, rho, U, phi); } -} -else -{ - surfaceScalarField phiHbyA + + phiHbyA *= fvc::interpolate(rho); + + + fvScalarMatrix pDDtEqn ( - "phiHbyA", - fvc::interpolate(rho) - *( - (fvc::interpolate(HbyA) & mesh.Sf()) - + fvc::ddtPhiCorr(rAU, rho, U, phi) - ) + fvc::ddt(rho) + psi*correction(fvm::ddt(p)) + + fvc::div(phiHbyA) + == + parcels.Srho() + + sources(psi, p, rho.name()) ); while (pimple.correctNonOrthogonal()) { fvScalarMatrix pEqn ( - fvm::ddt(psi, p) - + fvc::div(phiHbyA) + pDDtEqn - fvm::laplacian(rho*rAU, p) - == - parcels.Srho() - + sources(psi, p, rho.name()) ); + sources.constrain(pDDtEqn, rho.name()); + pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter()))); if (pimple.finalNonOrthogonalIter()) @@ -67,17 +45,28 @@ else phi = phiHbyA + pEqn.flux(); } } -} - -#include "rhoEqn.H" -#include "compressibleContinuityErrs.H" - -U = HbyA - rAU*fvc::grad(p); -U.correctBoundaryConditions(); -sources.correct(U); -K = 0.5*magSqr(U); - -if (thermo.dpdt()) -{ - dpdt = fvc::ddt(p); + + p.relax(); + + // Second part of thermodynamic density update + thermo.rho() += psi*p; + + #include "rhoEqn.H" // NOTE: flux and time scales now inconsistent + #include "compressibleContinuityErrs.H" + + U = HbyA - rAU*fvc::grad(p); + U.correctBoundaryConditions(); + sources.correct(U); + K = 0.5*magSqr(U); + + if (thermo.dpdt()) + { + dpdt = fvc::ddt(p); + } + + rho = thermo.rho(); + rho = max(rho, rhoMin); + rho = min(rho, rhoMax); + + Info<< "p min/max = " << min(p).value() << ", " << max(p).value() << endl; } diff --git a/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C b/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C index 9b0ccc1963..d1c735f658 100644 --- a/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C +++ b/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C @@ -22,22 +22,32 @@ License along with OpenFOAM. If not, see . Application - reactingParcelFoam + porousExplicitSourceReactingParcelFoam Description Transient PIMPLE solver for compressible, laminar or turbulent flow with - reacting Lagrangian parcels. + reacting multiphase Lagrangian parcels for porous media, including explicit + sources for mass, momentum and energy + + The solver includes: + - reacting multiphase parcel cloud + - porous media + - mass, momentum and energy sources + + Note: ddtPhiCorr not used here when porous zones are active + - not well defined for porous calculations \*---------------------------------------------------------------------------*/ #include "fvCFD.H" #include "turbulenceModel.H" -#include "basicReactingCloud.H" -#include "psiCombustionModel.H" +#include "basicReactingMultiphaseCloud.H" +#include "rhoCombustionModel.H" #include "radiationModel.H" +#include "porousZones.H" +#include "IObasicSourceList.H" #include "SLGThermo.H" #include "pimpleControl.H" -#include "IObasicSourceList.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -48,16 +58,19 @@ int main(int argc, char *argv[]) #include "createTime.H" #include "createMesh.H" #include "readGravitationalAcceleration.H" + + pimpleControl pimple(mesh); + #include "createFields.H" - #include "createClouds.H" #include "createRadiationModel.H" + #include "createClouds.H" + #include "createExplicitSources.H" + #include "createPorousZones.H" #include "initContinuityErrs.H" #include "readTimeControls.H" #include "compressibleCourantNo.H" #include "setInitialDeltaT.H" - pimpleControl pimple(mesh); - // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Info<< "\nStarting time loop\n" << endl; @@ -65,6 +78,7 @@ int main(int argc, char *argv[]) while (runTime.run()) { #include "readTimeControls.H" + #include "readAdditionalSolutionControls.H" #include "compressibleCourantNo.H" #include "setDeltaT.H" @@ -81,7 +95,7 @@ int main(int argc, char *argv[]) { #include "UEqn.H" #include "YEqn.H" - #include "hsEqn.H" + #include "EEqn.H" // --- Pressure corrector loop while (pimple.correct()) diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/readAdditionalSolutionControls.H b/applications/solvers/lagrangian/reactingParcelFoam/readAdditionalSolutionControls.H similarity index 57% rename from applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/readAdditionalSolutionControls.H rename to applications/solvers/lagrangian/reactingParcelFoam/readAdditionalSolutionControls.H index 050c24fe7a..d04dfc406d 100644 --- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/readAdditionalSolutionControls.H +++ b/applications/solvers/lagrangian/reactingParcelFoam/readAdditionalSolutionControls.H @@ -1,8 +1,4 @@ dictionary additional = mesh.solutionDict().subDict("additional"); -// pressure work term for enthalpy equation -bool pressureWorkTimeDerivative = - additional.lookupOrDefault("pressureWorkTimeDerivative", false); - // flag to activate solve transport for each specie (Y vector) bool solveSpecies = additional.lookupOrDefault("solveSpecies", true); diff --git a/applications/solvers/lagrangian/reactingParcelFoam/rhoEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/rhoEqn.H index b6293f2c1f..f142f620be 100644 --- a/applications/solvers/lagrangian/reactingParcelFoam/rhoEqn.H +++ b/applications/solvers/lagrangian/reactingParcelFoam/rhoEqn.H @@ -42,6 +42,9 @@ Description sources.constrain(rhoEqn); rhoEqn.solve(); + + Info<< "rho min/max = " << min(rho).value() << ", " << max(rho).value() + << endl; } // ************************************************************************* // diff --git a/applications/solvers/lagrangian/sprayFoam/Make/options b/applications/solvers/lagrangian/sprayFoam/Make/options index 0d5e4b6dd4..bcbb2fb5cb 100644 --- a/applications/solvers/lagrangian/sprayFoam/Make/options +++ b/applications/solvers/lagrangian/sprayFoam/Make/options @@ -1,5 +1,5 @@ EXE_INC = \ - -I$(FOAM_SOLVERS)/lagrangian/reactingParcelFoam \ + -I../reactingParcelFoam \ -I$(LIB_SRC)/finiteVolume/lnInclude \ -I${LIB_SRC}/meshTools/lnInclude \ -I${LIB_SRC}/sampling/lnInclude \ diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/UEqn.H b/applications/solvers/lagrangian/sprayFoam/UEqn.H similarity index 85% rename from applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/UEqn.H rename to applications/solvers/lagrangian/sprayFoam/UEqn.H index ce25309532..2ddaa2c301 100644 --- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/UEqn.H +++ b/applications/solvers/lagrangian/sprayFoam/UEqn.H @@ -1,6 +1,5 @@ fvVectorMatrix UEqn ( - //pZones.ddt(rho, U) fvm::ddt(rho, U) + fvm::div(phi, U) + turbulence->divDevRhoReff(U) @@ -13,11 +12,8 @@ sources.constrain(UEqn); - pZones.addResistance(UEqn); - if (pimple.momentumPredictor()) { solve(UEqn == -fvc::grad(p) + sources(rho, U)); K = 0.5*magSqr(U); } - diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createFields.H b/applications/solvers/lagrangian/sprayFoam/createFields.H similarity index 85% rename from applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createFields.H rename to applications/solvers/lagrangian/sprayFoam/createFields.H index 41d5401f08..612ea4f3c9 100644 --- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createFields.H +++ b/applications/solvers/lagrangian/sprayFoam/createFields.H @@ -1,11 +1,12 @@ Info<< "Creating combustion model\n" << endl; - autoPtr combustion + autoPtr combustion ( - combustionModels::rhoCombustionModel::New(mesh) + combustionModels::psiCombustionModel::New(mesh) ); - rhoReactionThermo& thermo = combustion->thermo(); + psiReactionThermo& thermo = combustion->thermo(); + thermo.validate(args.executable(), "h", "e"); SLGThermo slgThermo(mesh, thermo); @@ -23,7 +24,6 @@ } volScalarField& p = thermo.p(); - volScalarField& hs = thermo.he(); const volScalarField& T = thermo.T(); const volScalarField& psi = thermo.psi(); @@ -68,32 +68,9 @@ ) ); - // Set the turbulence into the combustion model + // Set the turbulence into the combustion model combustion->setTurbulence(turbulence()); - Info<< "Creating multi-variate interpolation scheme\n" << endl; - multivariateSurfaceInterpolationScheme::fieldTable fields; - - forAll(Y, i) - { - fields.add(Y[i]); - } - fields.add(hs); - - volScalarField dQ - ( - IOobject - ( - "dQ", - runTime.timeName(), - mesh, - IOobject::NO_READ, - IOobject::AUTO_WRITE - ), - mesh, - dimensionedScalar("dQ", dimEnergy/dimTime, 0.0) - ); - Info<< "Creating field dpdt\n" << endl; volScalarField dpdt ( @@ -109,3 +86,27 @@ Info<< "Creating field kinetic energy K\n" << endl; volScalarField K("K", 0.5*magSqr(U)); + multivariateSurfaceInterpolationScheme::fieldTable fields; + + forAll(Y, i) + { + fields.add(Y[i]); + } + fields.add(thermo.he()); + + volScalarField dQ + ( + IOobject + ( + "dQ", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar("dQ", dimEnergy/dimTime, 0.0) + ); + + Info<< "Creating sources\n" << endl; + IObasicSourceList sources(mesh); diff --git a/applications/solvers/lagrangian/sprayFoam/pEqn.H b/applications/solvers/lagrangian/sprayFoam/pEqn.H new file mode 100644 index 0000000000..d476a6bda8 --- /dev/null +++ b/applications/solvers/lagrangian/sprayFoam/pEqn.H @@ -0,0 +1,83 @@ +rho = thermo.rho(); + +volScalarField rAU(1.0/UEqn.A()); +volVectorField HbyA("HbyA", U); +HbyA = rAU*(UEqn == sources(rho, U))().H(); + +if (pimple.transonic()) +{ + surfaceScalarField phid + ( + "phid", + fvc::interpolate(psi) + *( + (fvc::interpolate(HbyA) & mesh.Sf()) + + fvc::ddtPhiCorr(rAU, rho, U, phi) + ) + ); + + while (pimple.correctNonOrthogonal()) + { + fvScalarMatrix pEqn + ( + fvm::ddt(psi, p) + + fvm::div(phid, p) + - fvm::laplacian(rho*rAU, p) + == + parcels.Srho() + + sources(psi, p, rho.name()) + ); + + pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter()))); + + if (pimple.finalNonOrthogonalIter()) + { + phi == pEqn.flux(); + } + } +} +else +{ + surfaceScalarField phiHbyA + ( + "phiHbyA", + fvc::interpolate(rho) + *( + (fvc::interpolate(HbyA) & mesh.Sf()) + + fvc::ddtPhiCorr(rAU, rho, U, phi) + ) + ); + + while (pimple.correctNonOrthogonal()) + { + fvScalarMatrix pEqn + ( + fvm::ddt(psi, p) + + fvc::div(phiHbyA) + - fvm::laplacian(rho*rAU, p) + == + parcels.Srho() + + sources(psi, p, rho.name()) + ); + + pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter()))); + + if (pimple.finalNonOrthogonalIter()) + { + phi = phiHbyA + pEqn.flux(); + } + } +} + +#include "rhoEqn.H" +#include "compressibleContinuityErrs.H" + +U = HbyA - rAU*fvc::grad(p); +U.correctBoundaryConditions(); +sources.correct(U); +K = 0.5*magSqr(U); + +if (thermo.dpdt()) +{ + dpdt = fvc::ddt(p); +} diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/rhoEqn.H b/applications/solvers/lagrangian/sprayFoam/rhoEqn.H similarity index 95% rename from applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/rhoEqn.H rename to applications/solvers/lagrangian/sprayFoam/rhoEqn.H index 9eacd97597..b6293f2c1f 100644 --- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/rhoEqn.H +++ b/applications/solvers/lagrangian/sprayFoam/rhoEqn.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License diff --git a/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/Make/options b/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/Make/options index d9ae7414b2..d2d83e3920 100644 --- a/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/Make/options +++ b/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/Make/options @@ -1,6 +1,6 @@ EXE_INC = \ - -I$(FOAM_SOLVERS)/lagrangian/sprayFoam \ - -I$(FOAM_SOLVERS)/lagrangian/reactingParcelFoam \ + -I.. \ + -I../../reactingParcelFoam \ -I$(LIB_SRC)/finiteVolume/lnInclude \ -I${LIB_SRC}/meshTools/lnInclude \ -I${LIB_SRC}/sampling/lnInclude \ diff --git a/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/createClouds.H b/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/createClouds.H deleted file mode 100644 index ee0985ff70..0000000000 --- a/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/createClouds.H +++ /dev/null @@ -1,9 +0,0 @@ -Info<< "\nConstructing reacting cloud" << endl; -basicSprayCloud parcels -( - "sprayCloud", - rho, - U, - g, - slgThermo -); diff --git a/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/sprayEngineFoam.C b/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/sprayEngineFoam.C index baf1e40dc7..06ed6e1c2b 100644 --- a/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/sprayEngineFoam.C +++ b/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/sprayEngineFoam.C @@ -60,6 +60,7 @@ int main(int argc, char *argv[]) #include "startSummary.H" pimpleControl pimple(mesh); + bool solveSpecies = true; // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -86,7 +87,7 @@ int main(int argc, char *argv[]) { #include "UEqn.H" #include "YEqn.H" - #include "hsEqn.H" + #include "EEqn.H" // --- Pressure corrector loop while (pimple.correct()) diff --git a/applications/solvers/lagrangian/sprayFoam/sprayFoam.C b/applications/solvers/lagrangian/sprayFoam/sprayFoam.C index fce272d042..78713032f3 100644 --- a/applications/solvers/lagrangian/sprayFoam/sprayFoam.C +++ b/applications/solvers/lagrangian/sprayFoam/sprayFoam.C @@ -57,6 +57,7 @@ int main(int argc, char *argv[]) #include "setInitialDeltaT.H" pimpleControl pimple(mesh); + bool solveSpecies = true; // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -81,7 +82,7 @@ int main(int argc, char *argv[]) { #include "UEqn.H" #include "YEqn.H" - #include "hsEqn.H" + #include "EEqn.H" // --- Pressure corrector loop while (pimple.correct()) diff --git a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H index 80011d1ab5..3197a72bee 100644 --- a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H +++ b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H @@ -5,6 +5,27 @@ volScalarField rAU("rAU", 1.0/UEqn.A()); surfaceScalarField rAUf("Dp", fvc::interpolate(rAU)); + volVectorField HbyA("HbyA", U); + HbyA = rAU*UEqn.H(); + + surfaceScalarField phiHbyA + ( + "phiHbyA", + (fvc::interpolate(HbyA) & mesh.Sf()) + + fvc::ddtPhiCorr(rAU, rho, U, phi) + ); + phi = phiHbyA; + + surfaceScalarField phig + ( + ( + fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1) + - ghf*fvc::snGrad(rho) + )*rAUf*mesh.magSf() + ); + + phiHbyA += phig; + tmp p_rghEqnComp1; tmp p_rghEqnComp2; @@ -27,27 +48,6 @@ - fvc::Sp(fvc::div(phid2), p_rgh); } - volVectorField HbyA("HbyA", U); - HbyA = rAU*UEqn.H(); - - surfaceScalarField phiHbyA - ( - "phiHbyA", - (fvc::interpolate(HbyA) & mesh.Sf()) - + fvc::ddtPhiCorr(rAU, rho, U, phi) - ); - phi = phiHbyA; - - surfaceScalarField phig - ( - ( - fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1) - - ghf*fvc::snGrad(rho) - )*rAUf*mesh.magSf() - ); - - phiHbyA += phig; - // Thermodynamic density needs to be updated by psi*d(p) after the // pressure solution - done in 2 parts. Part 1: //thermo.rho() -= psi*p_rgh; diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/TEqns.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/TEqns.H index 861590b235..8f38ca9d91 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/TEqns.H +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/TEqns.H @@ -6,7 +6,6 @@ ( fvm::ddt(alpha1, T1) + fvm::div(alphaPhi1, T1) - - fvm::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), T1) - fvm::laplacian(kByCp1, T1) == heatTransferCoeff*T2/Cp1/rho1 @@ -18,7 +17,6 @@ ( fvm::ddt(alpha2, T2) + fvm::div(alphaPhi2, T2) - - fvm::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), T2) - fvm::laplacian(kByCp2, T2) == heatTransferCoeff*T1/Cp2/rho2 diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/UEqns.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/UEqns.H index 5d61f0eaf7..606b4573e4 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/UEqns.H +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/UEqns.H @@ -30,7 +30,6 @@ fvVectorMatrix U2Eqn(U2, U2.dimensions()*dimVol/dimTime); ( fvm::ddt(alpha1, U1) + fvm::div(alphaPhi1, U1) - - fvm::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), U1) + Cvm*rho2*alpha1*alpha2/rho1* ( @@ -61,7 +60,6 @@ fvVectorMatrix U2Eqn(U2, U2.dimensions()*dimVol/dimTime); ( fvm::ddt(alpha2, U2) + fvm::div(alphaPhi2, U2) - - fvm::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), U2) + Cvm*rho2*alpha1*alpha2/rho2* ( diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/UEqn.H b/applications/solvers/multiphase/interPhaseChangeFoam/UEqn.H index 062c5523c9..7cc250a66a 100644 --- a/applications/solvers/multiphase/interPhaseChangeFoam/UEqn.H +++ b/applications/solvers/multiphase/interPhaseChangeFoam/UEqn.H @@ -2,7 +2,6 @@ ( fvm::ddt(rho, U) + fvm::div(rhoPhi, U) - - fvm::Sp(fvc::ddt(rho) + fvc::div(rhoPhi), U) + turbulence->divDevRhoReff(rho, U) ); diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/TEqns.H b/applications/solvers/multiphase/multiphaseEulerFoam/TEqns.H index 26015141d2..427a5f0b4c 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/TEqns.H +++ b/applications/solvers/multiphase/multiphaseEulerFoam/TEqns.H @@ -6,7 +6,6 @@ ( fvm::ddt(alpha1, T1) + fvm::div(alphaPhi1, T1) - - fvm::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), T1) - fvm::laplacian(kByCp1, T1) == heatTransferCoeff*T2/Cp1/rho1 @@ -18,7 +17,6 @@ ( fvm::ddt(alpha2, T2) + fvm::div(alphaPhi2, T2) - - fvm::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), T2) - fvm::laplacian(kByCp2, T2) == heatTransferCoeff*T1/Cp2/rho2 diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/UEqns.H b/applications/solvers/multiphase/multiphaseEulerFoam/UEqns.H index d4b371915c..e6478e4d89 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/UEqns.H +++ b/applications/solvers/multiphase/multiphaseEulerFoam/UEqns.H @@ -17,7 +17,6 @@ forAllIter(PtrDictionary, fluid.phases(), iter) ( fvm::ddt(alpha, U) + fvm::div(phase.phiAlpha(), U) - - fvm::Sp(fvc::ddt(alpha) + fvc::div(phase.phiAlpha()), U) + (alpha/phase.rho())*fluid.Cvm(phase)* ( diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/newVirualMass/DDtU.H b/applications/solvers/multiphase/twoPhaseEulerFoam/newVirtualMass/DDtU.H similarity index 100% rename from applications/solvers/multiphase/twoPhaseEulerFoam/newVirualMass/DDtU.H rename to applications/solvers/multiphase/twoPhaseEulerFoam/newVirtualMass/DDtU.H diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/newVirualMass/UEqns.H b/applications/solvers/multiphase/twoPhaseEulerFoam/newVirtualMass/UEqns.H similarity index 100% rename from applications/solvers/multiphase/twoPhaseEulerFoam/newVirualMass/UEqns.H rename to applications/solvers/multiphase/twoPhaseEulerFoam/newVirtualMass/UEqns.H diff --git a/applications/test/dataEntry/Make/files b/applications/test/dataEntry/Make/files index c514002f99..a88ccad3bc 100644 --- a/applications/test/dataEntry/Make/files +++ b/applications/test/dataEntry/Make/files @@ -1,3 +1,6 @@ Test-DataEntry.C +interpolationWeights.C +splineInterpolationWeights.C +linearInterpolationWeights.C EXE = $(FOAM_USER_APPBIN)/Test-DataEntry diff --git a/applications/test/dataEntry/Make/options b/applications/test/dataEntry/Make/options index a071d9557c..72cae9f645 100644 --- a/applications/test/dataEntry/Make/options +++ b/applications/test/dataEntry/Make/options @@ -1,4 +1,5 @@ EXE_INC = \ + -DFULLDEBUG -g -O0 \ -I$(LIB_SRC)/finiteVolume/lnInclude \ -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/radiationModels/lnInclude \ diff --git a/applications/test/dataEntry/Test-DataEntry.C b/applications/test/dataEntry/Test-DataEntry.C index 92761b1d01..a47a29be24 100644 --- a/applications/test/dataEntry/Test-DataEntry.C +++ b/applications/test/dataEntry/Test-DataEntry.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -32,6 +32,8 @@ Description #include "fvCFD.H" #include "DataEntry.H" #include "IOdictionary.H" +#include "linearInterpolationWeights.H" +#include "splineInterpolationWeights.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -41,6 +43,62 @@ int main(int argc, char *argv[]) # include "createTime.H" # include "createMesh.H" +{ + scalarField samples(4); + samples[0] = 0; + samples[1] = 1; + samples[2] = 2; + samples[3] = 3; + scalarField values(4); + values = 1.0; + //values[0] = 0.0; + //values[1] = 1.0; + + //linearInterpolationWeights interpolator + splineInterpolationWeights interpolator + ( + samples, + interpolationWeights::WARN + ); + labelList indices; + scalarField weights; + + interpolator.integrationWeights(1.1, 1.2, indices, weights); + Pout<< "indices:" << indices << endl; + Pout<< "weights:" << weights << endl; + + scalar baseSum = interpolator.weightedSum + ( + weights, + UIndirectList(values, indices) + ); + Pout<< "baseSum=" << baseSum << nl << nl << endl; + + +// interpolator.integrationWeights(-0.01, 0, indices, weights); +// scalar partialSum = interpolator.weightedSum +// ( +// weights, +// UIndirectList(values, indices) +// ); +// Pout<< "partialSum=" << partialSum << nl << nl << endl; +// +// +// interpolator.integrationWeights(-0.01, 1, indices, weights); +// //Pout<< "samples:" << samples << endl; +// //Pout<< "indices:" << indices << endl; +// //Pout<< "weights:" << weights << endl; +// scalar sum = interpolator.weightedSum +// ( +// weights, +// UIndirectList(values, indices) +// ); +// Pout<< "integrand=" << sum << nl << nl << endl; + + + return 1; +} + IOdictionary dataEntryProperties ( IOobject diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict index 05dcb07286..8d7927c42b 100644 --- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict +++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict @@ -239,7 +239,8 @@ snapControls // Leave out altogether to disable. nFeatureSnapIter 10; - //- Detect (geometric) features by sampling the surface (default=false) + //- Detect (geometric only) features by sampling the surface + // (default=false). implicitFeatureSnap false; //- Use castellatedMeshControls::features (default = true) diff --git a/applications/utilities/mesh/manipulation/createPatch/createPatchDict b/applications/utilities/mesh/manipulation/createPatch/createPatchDict index f731f78df9..62e6cf619c 100644 --- a/applications/utilities/mesh/manipulation/createPatch/createPatchDict +++ b/applications/utilities/mesh/manipulation/createPatch/createPatchDict @@ -92,8 +92,10 @@ patches // Optional: explicitly set transformation tensor. // Used when matching and synchronising points. transform rotational; - rotationAxis ( 0 0 1 ); - rotationCentre ( 0.3 0 0 ); + rotationAxis (1 0 0); + rotationCentre (0 0 0); + // transform translational; + // separationVector (1 0 0); } // How to construct: either from 'patches' or 'set' diff --git a/applications/utilities/mesh/manipulation/moveDynamicMesh/moveDynamicMesh.C b/applications/utilities/mesh/manipulation/moveDynamicMesh/moveDynamicMesh.C index 00c97614de..052cd91210 100644 --- a/applications/utilities/mesh/manipulation/moveDynamicMesh/moveDynamicMesh.C +++ b/applications/utilities/mesh/manipulation/moveDynamicMesh/moveDynamicMesh.C @@ -107,6 +107,7 @@ void writeWeights(const polyMesh& mesh) int main(int argc, char *argv[]) { +# include "addRegionOption.H" argList::addBoolOption ( "checkAMI", @@ -115,7 +116,7 @@ int main(int argc, char *argv[]) # include "setRootCase.H" # include "createTime.H" -# include "createDynamicFvMesh.H" +# include "createNamedDynamicFvMesh.H" const bool checkAMI = args.optionFound("checkAMI"); diff --git a/applications/utilities/mesh/manipulation/topoSet/topoSetDict b/applications/utilities/mesh/manipulation/topoSet/topoSetDict index 7a2049c686..31fa348e97 100644 --- a/applications/utilities/mesh/manipulation/topoSet/topoSetDict +++ b/applications/utilities/mesh/manipulation/topoSet/topoSetDict @@ -359,6 +359,20 @@ FoamFile // cellSet c0; // name of cellSet of slave side // } // +// +// pointZoneSet +// ~~~~~~~~~~~~ +// (mirrors operations on a pointSet into a pointZone) +// +// // Select based on pointSet +// source setToPointZone; +// sourceInfo +// { +// set p0; // name of pointSet +// } +// +// +// actions ( diff --git a/applications/utilities/preProcessing/changeDictionary/changeDictionary.C b/applications/utilities/preProcessing/changeDictionary/changeDictionary.C index ad178b765c..3c664cfc3b 100644 --- a/applications/utilities/preProcessing/changeDictionary/changeDictionary.C +++ b/applications/utilities/preProcessing/changeDictionary/changeDictionary.C @@ -52,6 +52,7 @@ Description } } \endverbatim + Replacement entries starting with '~' will remove the entry. Usage @@ -172,6 +173,46 @@ bool addEntry } + +// List of indices into thisKeys +labelList findMatches +( + const HashTable& shortcuts, + const wordList& shortcutNames, + const wordList& thisKeys, + const keyType& key +) +{ + labelList matches; + + if (key.isPattern()) + { + // Wildcard match + matches = findStrings(key, thisKeys); + + } + else if (shortcuts.size()) + { + // See if patchGroups expand to valid thisKeys + labelList indices = findStrings(key, shortcutNames); + forAll(indices, i) + { + const word& name = shortcutNames[indices[i]]; + const wordList& keys = shortcuts[name]; + forAll(keys, j) + { + label index = findIndex(thisKeys, keys[j]); + if (index != -1) + { + matches.append(index); + } + } + } + } + return matches; +} + + // Dictionary merging/editing. // literalRE: // - true: behave like dictionary::merge, i.e. add regexps just like @@ -185,6 +226,8 @@ bool merge const HashTable& shortcuts ) { + const wordList shortcutNames(shortcuts.toc()); + bool changed = false; // Save current (non-wildcard) keys before adding items. @@ -203,7 +246,18 @@ bool merge { const keyType& key = mergeIter().keyword(); - if (literalRE || !(key.isPattern() || shortcuts.found(key))) + if (key[0] == '~') + { + word eraseKey = key(1, key.size()-1); + if (thisDict.remove(eraseKey)) + { + // Mark thisDict entry as having been match for wildcard + // handling later on. + thisKeysSet.erase(eraseKey); + } + changed = true; + } + else if (literalRE || !(key.isPattern() || shortcuts.found(key))) { entry* entryPtr = thisDict.lookupEntryPtr ( @@ -255,59 +309,69 @@ bool merge { const keyType& key = mergeIter().keyword(); - // List of indices into thisKeys - labelList matches; - - if (key.isPattern()) + if (key[0] == '~') { - // Wildcard match - matches = findStrings(key, thisKeys); + word eraseKey = key(1, key.size()-1); - } - else if (shortcuts.size()) - { - // See if patchGroups expand to valid thisKeys - const wordList shortcutNames = shortcuts.toc(); - labelList indices = findStrings(key, shortcutNames); - forAll(indices, i) - { - const word& name = shortcutNames[indices[i]]; - const wordList& keys = shortcuts[name]; - forAll(keys, j) - { - label index = findIndex(thisKeys, keys[j]); - if (index != -1) - { - matches.append(index); - } - } - } - } - - // Add all matches - forAll(matches, i) - { - const word& thisKey = thisKeys[matches[i]]; - - entry& thisEntry = const_cast + // List of indices into thisKeys + labelList matches ( - thisDict.lookupEntry(thisKey, false, false) + findMatches + ( + shortcuts, + shortcutNames, + thisKeys, + eraseKey + ) ); - if - ( - addEntry - ( - thisDict, - thisEntry, - mergeIter(), - literalRE, - HashTable(0) // no shortcuts - // at deeper levels - ) - ) + // Remove all matches + forAll(matches, i) { - changed = true; + const word& thisKey = thisKeys[matches[i]]; + thisKeysSet.erase(thisKey); + } + changed = true; + } + else + { + // List of indices into thisKeys + labelList matches + ( + findMatches + ( + shortcuts, + shortcutNames, + thisKeys, + key + ) + ); + + // Add all matches + forAll(matches, i) + { + const word& thisKey = thisKeys[matches[i]]; + + entry& thisEntry = const_cast + ( + thisDict.lookupEntry(thisKey, false, false) + ); + + if + ( + addEntry + ( + thisDict, + thisEntry, + mergeIter(), + literalRE, + HashTable(0) // no shortcuts + // at deeper levels + ) + ) + { + changed = true; + } } } } diff --git a/applications/utilities/thermophysical/adiabaticFlameT/adiabaticFlameT.C b/applications/utilities/thermophysical/adiabaticFlameT/adiabaticFlameT.C index 58de0b6f77..cfdb984b6d 100644 --- a/applications/utilities/thermophysical/adiabaticFlameT/adiabaticFlameT.C +++ b/applications/utilities/thermophysical/adiabaticFlameT/adiabaticFlameT.C @@ -36,15 +36,16 @@ Description #include "IFstream.H" #include "OSspecific.H" -#include "specieThermo.H" -#include "absoluteEnthalpy.H" -#include "janafThermo.H" +#include "specie.H" #include "perfectGas.H" +#include "thermo.H" +#include "janafThermo.H" +#include "absoluteEnthalpy.H" using namespace Foam; -typedef specieThermo, absoluteEnthalpy> thermo; - +typedef species::thermo >, absoluteEnthalpy> + thermo; // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/applications/utilities/thermophysical/equilibriumCO/equilibriumCO.C b/applications/utilities/thermophysical/equilibriumCO/equilibriumCO.C index dde7e8678a..8bf51d72c2 100644 --- a/applications/utilities/thermophysical/equilibriumCO/equilibriumCO.C +++ b/applications/utilities/thermophysical/equilibriumCO/equilibriumCO.C @@ -35,15 +35,18 @@ Description #include "OSspecific.H" #include "IOmanip.H" -#include "specieThermo.H" -#include "absoluteEnthalpy.H" -#include "janafThermo.H" +#include "specie.H" #include "perfectGas.H" +#include "thermo.H" +#include "janafThermo.H" +#include "absoluteEnthalpy.H" + #include "SLPtrList.H" using namespace Foam; -typedef specieThermo, absoluteEnthalpy> thermo; +typedef species::thermo >, absoluteEnthalpy> + thermo; // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // Main program: diff --git a/applications/utilities/thermophysical/equilibriumFlameT/equilibriumFlameT.C b/applications/utilities/thermophysical/equilibriumFlameT/equilibriumFlameT.C index 8cf2371898..40452749a4 100644 --- a/applications/utilities/thermophysical/equilibriumFlameT/equilibriumFlameT.C +++ b/applications/utilities/thermophysical/equilibriumFlameT/equilibriumFlameT.C @@ -38,14 +38,16 @@ Description #include "OSspecific.H" #include "IOmanip.H" -#include "specieThermo.H" -#include "absoluteEnthalpy.H" -#include "janafThermo.H" +#include "specie.H" #include "perfectGas.H" +#include "thermo.H" +#include "janafThermo.H" +#include "absoluteEnthalpy.H" using namespace Foam; -typedef specieThermo, absoluteEnthalpy> thermo; +typedef species::thermo >, absoluteEnthalpy> + thermo; // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/applications/utilities/thermophysical/mixtureAdiabaticFlameT/mixtureAdiabaticFlameT.C b/applications/utilities/thermophysical/mixtureAdiabaticFlameT/mixtureAdiabaticFlameT.C index 728897c88a..7f1def4476 100644 --- a/applications/utilities/thermophysical/mixtureAdiabaticFlameT/mixtureAdiabaticFlameT.C +++ b/applications/utilities/thermophysical/mixtureAdiabaticFlameT/mixtureAdiabaticFlameT.C @@ -35,16 +35,17 @@ Description #include "IFstream.H" #include "OSspecific.H" -#include "specieThermo.H" -#include "absoluteEnthalpy.H" -#include "janafThermo.H" +#include "specie.H" #include "perfectGas.H" +#include "thermo.H" +#include "janafThermo.H" +#include "absoluteEnthalpy.H" #include "mixture.H" using namespace Foam; -typedef specieThermo, absoluteEnthalpy> thermo; - +typedef species::thermo >, absoluteEnthalpy> + thermo; // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/etc/controlDict b/etc/controlDict index 816febf3ec..f9998bcdd7 100644 --- a/etc/controlDict +++ b/etc/controlDict @@ -510,30 +510,30 @@ DebugSwitches gradientUnburntEnthalpy 0; granularPressureModel 0; hCombustionThermo 0; - hMixtureThermo>>>> 0; - hMixtureThermo>>>> 0; - hMixtureThermo>>>> 0; - hMixtureThermo>>>> 0; - hMixtureThermo>>>> 0; + hMixtureThermo>>>> 0; + hMixtureThermo>>>> 0; + hMixtureThermo>>>> 0; + hMixtureThermo>>>> 0; + hMixtureThermo>>>> 0; hMixtureThermo 0; - hMixtureThermo>>>> 0; - hMixtureThermo>>>> 0; - hThermo>>>> 0; - hThermo>>>> 0; - hThermo>>>> 0; + hMixtureThermo>>>> 0; + hMixtureThermo>>>> 0; + hThermo>>>> 0; + hThermo>>>> 0; + hThermo>>>> 0; harmonic 0; heatTransferModel 0; hexCellLooper 0; hexRef8 0; hhuCombustionThermo 0; - hhuMixtureThermo>>>> 0; - hhuMixtureThermo>>>> 0; - hhuMixtureThermo>>>> 0; - hhuMixtureThermo>>>> 0; - hhuMixtureThermo>>>> 0; - hhuMixtureThermo>>>> 0; - hhuMixtureThermo>>>> 0; - hhuMixtureThermo>>>> 0; + hhuMixtureThermo>>>> 0; + hhuMixtureThermo>>>> 0; + hhuMixtureThermo>>>> 0; + hhuMixtureThermo>>>> 0; + hhuMixtureThermo>>>> 0; + hhuMixtureThermo>>>> 0; + hhuMixtureThermo>>>> 0; + hhuMixtureThermo>>>> 0; hierarchical 0; hollowConeInjector 0; iC3H8O 0; @@ -711,9 +711,9 @@ DebugSwitches processor 0; processorLduInterface 0; processorLduInterfaceField 0; - pureMixture>>> 0; - pureMixture>>> 0; - pureMixture>>> 0; + pureMixture>>> 0; + pureMixture>>> 0; + pureMixture>>> 0; quadratic 0; quadraticFit 0; quadraticLinearFit 0; @@ -965,4 +965,34 @@ DimensionedConstants } +DimensionSets +{ + unitSet SI; // USCS + + SICoeffs + { + // Basic units + kg kg [ 1 0 0 0 0 0 0 ] 1.0; + m m [ 0 1 0 0 0 0 0 ] 1.0; + s s [ 0 0 1 0 0 0 0 ] 1.0; + K K [ 0 0 0 1 0 0 0 ] 1.0; + mol mol [ 0 0 0 0 1 0 0 ] 1.0; + A A [ 0 0 0 0 0 1 0 ] 1.0; + Cd Cd [ 0 0 0 0 0 0 1 ] 1.0; + + // Derived units + Pa Pa [ kg^1 m^-2 ] 1.0; + + // Scaled units + mm mm [ kg^1 m^-2 ] 1e-3; + + + // Set of units used for printing. Can be any basic or derived + // but not scaled (only supported for dimensionedScalar, etc) + printUnits (kg m s K mol A Cd); + } +} + + + // ************************************************************************* // diff --git a/src/OSspecific/POSIX/fileMonitor.C b/src/OSspecific/POSIX/fileMonitor.C index 030e7639f4..f17e7b1dc9 100644 --- a/src/OSspecific/POSIX/fileMonitor.C +++ b/src/OSspecific/POSIX/fileMonitor.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -168,7 +168,7 @@ namespace Foam FatalErrorIn("fileMonitorWatcher(const bool, const label)") << "You selected inotify but this file was compiled" << " without FOAM_USE_INOTIFY" - << "Please select another fileModification test method" + << " Please select another fileModification test method" << exit(FatalError); #endif } diff --git a/src/OpenFOAM/primitives/functions/DataEntry/CompatibilityConstant/CompatibilityConstant.C b/src/OpenFOAM/primitives/functions/DataEntry/CompatibilityConstant/CompatibilityConstant.C index 78c6a4d7c5..953328f031 100644 --- a/src/OpenFOAM/primitives/functions/DataEntry/CompatibilityConstant/CompatibilityConstant.C +++ b/src/OpenFOAM/primitives/functions/DataEntry/CompatibilityConstant/CompatibilityConstant.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -30,7 +30,8 @@ License template Foam::CompatibilityConstant::CompatibilityConstant ( - const word& entryName, const dictionary& dict + const word& entryName, + const dictionary& dict ) : DataEntry(entryName), @@ -69,6 +70,7 @@ Foam::CompatibilityConstant::CompatibilityConstant dimensions_(cnst.dimensions_) {} + // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // template diff --git a/src/OpenFOAM/primitives/functions/DataEntry/Table/Table.H b/src/OpenFOAM/primitives/functions/DataEntry/Table/Table.H index 63d7b293e3..e89234cfa9 100644 --- a/src/OpenFOAM/primitives/functions/DataEntry/Table/Table.H +++ b/src/OpenFOAM/primitives/functions/DataEntry/Table/Table.H @@ -27,13 +27,13 @@ Class Description Templated table container data entry. Items are stored in a list of Tuple2's. First column is always stored as scalar entries. Data is read - in the form, e.g. for an entry \ that is (scalar, vector): + in Tuple2 form, e.g. for an entry \ that is (scalar, vector): \verbatim - table [0 1 0 0 0] //dimension set optional + table ( - 0.0 (1 2 3) - 1.0 (4 5 6) + (0.0 (1 2 3)) + (1.0 (4 5 6)) ); \endverbatim diff --git a/src/combustionModels/singleStepCombustion/singleStepCombustion.C b/src/combustionModels/singleStepCombustion/singleStepCombustion.C index e3754f0a58..8234e573c5 100644 --- a/src/combustionModels/singleStepCombustion/singleStepCombustion.C +++ b/src/combustionModels/singleStepCombustion/singleStepCombustion.C @@ -119,7 +119,7 @@ singleStepCombustion::R { const label fNorm = singleMixturePtr_->specieProd()[specieI]; const volScalarField fres(singleMixturePtr_->fres(specieI)); - wSpecie /= max(fNorm*(Y - fres), 1e-2); + wSpecie /= max(fNorm*(Y - fres), scalar(1e-2)); return -fNorm*wSpecie*fres + fNorm*fvm::Sp(wSpecie, Y); } diff --git a/src/dynamicFvMesh/include/createNamedDynamicFvMesh.H b/src/dynamicFvMesh/include/createNamedDynamicFvMesh.H new file mode 100644 index 0000000000..8ffc1490f5 --- /dev/null +++ b/src/dynamicFvMesh/include/createNamedDynamicFvMesh.H @@ -0,0 +1,32 @@ + Foam::word regionName; + + if (args.optionReadIfPresent("region", regionName)) + { + Foam::Info + << "Create mesh " << regionName << " for time = " + << runTime.timeName() << Foam::nl << Foam::endl; + } + else + { + regionName = Foam::fvMesh::defaultRegion; + Foam::Info + << "Create mesh for time = " + << runTime.timeName() << Foam::nl << Foam::endl; + } + + + autoPtr meshPtr + ( + dynamicFvMesh::New + ( + IOobject + ( + regionName, + runTime.timeName(), + runTime, + IOobject::MUST_READ + ) + ) + ); + + dynamicFvMesh& mesh = meshPtr(); diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files index 17abbc6363..6e19f299f0 100644 --- a/src/finiteVolume/Make/files +++ b/src/finiteVolume/Make/files @@ -308,6 +308,7 @@ $(ddtSchemes)/backwardDdtScheme/backwardDdtSchemes.C $(ddtSchemes)/boundedBackwardDdtScheme/boundedBackwardDdtScheme.C $(ddtSchemes)/boundedBackwardDdtScheme/boundedBackwardDdtSchemes.C $(ddtSchemes)/CrankNicholsonDdtScheme/CrankNicholsonDdtSchemes.C +$(ddtSchemes)/boundedDdtScheme/boundedDdtSchemes.C d2dt2Schemes = finiteVolume/d2dt2Schemes $(d2dt2Schemes)/d2dt2Scheme/d2dt2Schemes.C diff --git a/src/finiteVolume/fields/fvPatchFields/derived/flowRateInletVelocity/flowRateInletVelocityFvPatchVectorField.C b/src/finiteVolume/fields/fvPatchFields/derived/flowRateInletVelocity/flowRateInletVelocityFvPatchVectorField.C index bbd50fa4de..f5cd4413b1 100644 --- a/src/finiteVolume/fields/fvPatchFields/derived/flowRateInletVelocity/flowRateInletVelocityFvPatchVectorField.C +++ b/src/finiteVolume/fields/fvPatchFields/derived/flowRateInletVelocity/flowRateInletVelocityFvPatchVectorField.C @@ -40,8 +40,9 @@ flowRateInletVelocityFvPatchVectorField : fixedValueFvPatchField(p, iF), flowRate_(), - phiName_("phi"), - rhoName_("rho") + volumetric_(false), + rhoName_("rho"), + rhoInlet_(0.0) {} @@ -56,8 +57,9 @@ flowRateInletVelocityFvPatchVectorField : fixedValueFvPatchField(ptf, p, iF, mapper), flowRate_(ptf.flowRate_().clone().ptr()), - phiName_(ptf.phiName_), - rhoName_(ptf.rhoName_) + volumetric_(ptf.volumetric_), + rhoName_(ptf.rhoName_), + rhoInlet_(ptf.rhoInlet_) {} @@ -69,11 +71,47 @@ flowRateInletVelocityFvPatchVectorField const dictionary& dict ) : - fixedValueFvPatchField(p, iF, dict), - flowRate_(DataEntry::New("flowRate", dict)), - phiName_(dict.lookupOrDefault("phi", "phi")), - rhoName_(dict.lookupOrDefault("rho", "rho")) -{} + fixedValueFvPatchField(p, iF), + rhoInlet_(dict.lookupOrDefault("rhoInlet", -VGREAT)) +{ + if (dict.found("volumetricFlowRate")) + { + volumetric_ = true; + flowRate_ = DataEntry::New("volumetricFlowRate", dict); + rhoName_ = "rho"; + } + else if (dict.found("massFlowRate")) + { + volumetric_ = false; + flowRate_ = DataEntry::New("massFlowRate", dict); + rhoName_ = word(dict.lookupOrDefault("rho", "rho")); + } + else + { + FatalIOErrorIn + ( + "flowRateInletVelocityFvPatchVectorField::" + "flowRateInletVelocityFvPatchVectorField" + "(const fvPatch&, const DimensionedField&," + " const dictionary&)", + dict + ) << "Please supply either 'volumetricFlowRate' or" + << " 'massFlowRate' and 'rho'" << exit(FatalIOError); + } + + // Value field require if mass based + if (dict.found("value")) + { + fvPatchField::operator= + ( + vectorField("value", dict, p.size()) + ); + } + else + { + evaluate(Pstream::blocking); + } +} Foam::flowRateInletVelocityFvPatchVectorField:: @@ -84,8 +122,9 @@ flowRateInletVelocityFvPatchVectorField : fixedValueFvPatchField(ptf), flowRate_(ptf.flowRate_().clone().ptr()), - phiName_(ptf.phiName_), - rhoName_(ptf.rhoName_) + volumetric_(ptf.volumetric_), + rhoName_(ptf.rhoName_), + rhoInlet_(ptf.rhoInlet_) {} @@ -98,13 +137,44 @@ flowRateInletVelocityFvPatchVectorField : fixedValueFvPatchField(ptf, iF), flowRate_(ptf.flowRate_().clone().ptr()), - phiName_(ptf.phiName_), - rhoName_(ptf.rhoName_) + volumetric_(ptf.volumetric_), + rhoName_(ptf.rhoName_), + rhoInlet_(ptf.rhoInlet_) {} // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +void Foam::flowRateInletVelocityFvPatchVectorField::updateCoeffs +( + const scalar uniformRho +) +{ + if (updated()) + { + return; + } + + const scalar t = db().time().timeOutputValue(); + + // a simpler way of doing this would be nice + const scalar avgU = -flowRate_->value(t)/gSum(patch().magSf()); + + tmp n = patch().nf(); + + if (volumetric_ || rhoName_ == "none") + { + // volumetric flow-rate + operator==(n*avgU); + } + else + { + // mass flow-rate + operator==(n*avgU/uniformRho); + } +} + + void Foam::flowRateInletVelocityFvPatchVectorField::updateCoeffs() { if (updated()) @@ -119,40 +189,38 @@ void Foam::flowRateInletVelocityFvPatchVectorField::updateCoeffs() tmp n = patch().nf(); - const surfaceScalarField& phi = - db().lookupObject(phiName_); - - if (phi.dimensions() == dimVelocity*dimArea) + if (volumetric_ || rhoName_ == "none") { - // volumetric flow-rate + // volumetric flow-rate or density not given operator==(n*avgU); } - else if (phi.dimensions() == dimDensity*dimVelocity*dimArea) + else { - if (rhoName_ == "none") + // mass flow-rate + if + ( + patch().boundaryMesh().mesh().foundObject(rhoName_) + ) { - // volumetric flow-rate if density not given - operator==(n*avgU); - } - else - { - // mass flow-rate const fvPatchField& rhop = patch().lookupPatchField(rhoName_); operator==(n*avgU/rhop); } - } - else - { - FatalErrorIn - ( - "flowRateInletVelocityFvPatchVectorField::updateCoeffs()" - ) << "dimensions of " << phiName_ << " are incorrect" << nl - << " on patch " << this->patch().name() - << " of field " << this->dimensionedInternalField().name() - << " in file " << this->dimensionedInternalField().objectPath() - << nl << exit(FatalError); + else + { + // Use constant density + if (rhoInlet_ < 0) + { + FatalErrorIn + ( + "flowRateInletVelocityFvPatchVectorField::updateCoeffs()" + ) << "Did not find registered density field " << rhoName_ + << " and no constant density 'rhoInlet' specified" + << exit(FatalError); + } + operator==(n*avgU/rhoInlet_); + } } fixedValueFvPatchField::updateCoeffs(); @@ -163,8 +231,11 @@ void Foam::flowRateInletVelocityFvPatchVectorField::write(Ostream& os) const { fvPatchField::write(os); flowRate_->writeData(os); - writeEntryIfDifferent(os, "phi", "phi", phiName_); - writeEntryIfDifferent(os, "rho", "rho", rhoName_); + if (!volumetric_) + { + writeEntryIfDifferent(os, "rho", "rho", rhoName_); + writeEntryIfDifferent(os, "rhoInlet", -VGREAT, rhoInlet_); + } writeEntry("value", os); } diff --git a/src/finiteVolume/fields/fvPatchFields/derived/flowRateInletVelocity/flowRateInletVelocityFvPatchVectorField.H b/src/finiteVolume/fields/fvPatchFields/derived/flowRateInletVelocity/flowRateInletVelocityFvPatchVectorField.H index a4b48c7629..1a5d2f493d 100644 --- a/src/finiteVolume/fields/fvPatchFields/derived/flowRateInletVelocity/flowRateInletVelocityFvPatchVectorField.H +++ b/src/finiteVolume/fields/fvPatchFields/derived/flowRateInletVelocity/flowRateInletVelocityFvPatchVectorField.H @@ -28,21 +28,26 @@ Description Describes a volumetric/mass flow normal vector boundary condition by its magnitude as an integral over its area. - The basis of the patch (volumetric or mass) is determined by the - dimensions of the flux, phi. - - If the flux is mass-based - - the current density is used to correct the velocity - - volumetric flow rate can be applied by setting the 'rho' entry to 'none' + Either specify 'volumetricFlowRate' or 'massFlowRate' (requires additional + 'rho' or 'rhoInlet' entry). Example of the boundary condition specification: \verbatim inlet { - type flowRateInletVelocity; - flowRate 0.2; // Volumetric/mass flow rate [m3/s or kg/s] - rho rho; // none | rho [m3/s or kg/s] - value uniform (0 0 0); // placeholder + type flowRateInletVelocity; + volumetricFlowRate 0.2; // Volumetric [m3/s] + } + \endverbatim + + \verbatim + inlet + { + type flowRateInletVelocity; + massFlowRate 0.2; // mass flow rate [kg/s] + rho rho; // rho [m3/s or kg/s] + rhoInlet 1.0 // uniform rho if no rho field registered + // (e.g. at startup) } \endverbatim @@ -79,12 +84,15 @@ class flowRateInletVelocityFvPatchVectorField //- Inlet integral flow rate autoPtr > flowRate_; - //- Name of the flux transporting the field - word phiName_; + //- Is volumetric? + bool volumetric_; //- Name of the density field used to normalize the mass flux word rhoName_; + //- Rho initialisation value (for start; if value not supplied) + scalar rhoInlet_; + public: @@ -157,6 +165,10 @@ public: // Member functions + //- Update the coefficients associated with the patch field given + // uniform density field + void updateCoeffs(const scalar uniformRho); + //- Update the coefficients associated with the patch field virtual void updateCoeffs(); diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/boundedDdtScheme/boundedDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/boundedDdtScheme/boundedDdtScheme.C new file mode 100644 index 0000000000..a89629f6db --- /dev/null +++ b/src/finiteVolume/finiteVolume/ddtSchemes/boundedDdtScheme/boundedDdtScheme.C @@ -0,0 +1,170 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation + \\/ 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 3 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, see . + +\*---------------------------------------------------------------------------*/ + +#include "boundedDdtScheme.H" +#include "fvcDiv.H" +#include "fvcDdt.H" +#include "fvMatrices.H" +#include "fvmSup.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace fv +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +template +tmp > +boundedDdtScheme::fvcDdt +( + const dimensioned& dt +) +{ + return scheme_().fvcDdt(dt); +} + + +template +tmp > +boundedDdtScheme::fvcDdt +( + const GeometricField& vf +) +{ + return scheme_().fvcDdt(vf); +} + + +template +tmp > +boundedDdtScheme::fvcDdt +( + const dimensionedScalar& rho, + const GeometricField& vf +) +{ + return scheme_().fvcDdt(rho, vf); +} + + +template +tmp > +boundedDdtScheme::fvcDdt +( + const volScalarField& rho, + const GeometricField& vf +) +{ + return scheme_().fvcDdt(rho, vf) - fvc::ddt(rho)*vf; +} + + +template +tmp > +boundedDdtScheme::fvmDdt +( + const GeometricField& vf +) +{ + return scheme_().fvmDdt(vf); +} + + +template +tmp > +boundedDdtScheme::fvmDdt +( + const dimensionedScalar& rho, + const GeometricField& vf +) +{ + return scheme_().fvmDdt(rho, vf); +} + + +template +tmp > +boundedDdtScheme::fvmDdt +( + const volScalarField& rho, + const GeometricField& vf +) +{ + return scheme_().fvmDdt(rho, vf) - fvm::Sp(fvc::ddt(rho), vf); +} + + +template +tmp::fluxFieldType> +boundedDdtScheme::fvcDdtPhiCorr +( + const volScalarField& rA, + const GeometricField& U, + const fluxFieldType& phi +) +{ + return scheme_().fvcDdtPhiCorr(rA, U, phi); +} + + +template +tmp::fluxFieldType> +boundedDdtScheme::fvcDdtPhiCorr +( + const volScalarField& rA, + const volScalarField& rho, + const GeometricField& U, + const fluxFieldType& phi +) +{ + return scheme_().fvcDdtPhiCorr(rA, rho, U, phi); +} + + +template +tmp boundedDdtScheme::meshPhi +( + const GeometricField& vf +) +{ + return scheme_().meshPhi(vf); +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace fv + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/boundedDdtScheme/boundedDdtScheme.H b/src/finiteVolume/finiteVolume/ddtSchemes/boundedDdtScheme/boundedDdtScheme.H new file mode 100644 index 0000000000..0e5e0831f9 --- /dev/null +++ b/src/finiteVolume/finiteVolume/ddtSchemes/boundedDdtScheme/boundedDdtScheme.H @@ -0,0 +1,207 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation + \\/ 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 3 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, see . + +Class + Foam::fv::boundedDdtScheme + +Description + Bounded form of the selected ddt scheme. + + Boundedness is achieved by subtracting ddt(phi)*vf or Sp(ddt(rho), vf) + which is non-conservative if ddt(rho) != 0 but conservative otherwise. + + Can be used for the ddt of bounded scalar properties to improve stability + if insufficient convergence of the pressure equation causes temporary + divergence of the flux field. + +SourceFiles + boundedDdtScheme.C + +\*---------------------------------------------------------------------------*/ + +#ifndef boundedDdtScheme_H +#define boundedDdtScheme_H + +#include "ddtScheme.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace fv +{ + +/*---------------------------------------------------------------------------*\ + Class boundedDdtScheme Declaration +\*---------------------------------------------------------------------------*/ + +template +class boundedDdtScheme +: + public fv::ddtScheme +{ + // Private data + + tmp > scheme_; + + + // Private Member Functions + + //- Disallow default bitwise copy construct + boundedDdtScheme(const boundedDdtScheme&); + + //- Disallow default bitwise assignment + void operator=(const boundedDdtScheme&); + + +public: + + //- Runtime type information + TypeName("bounded"); + + + // Constructors + + //- Construct from mesh and Istream + boundedDdtScheme(const fvMesh& mesh, Istream& is) + : + ddtScheme(mesh, is), + scheme_ + ( + fv::ddtScheme::New(mesh, is) + ) + {} + + + // Member Functions + + //- Return mesh reference + const fvMesh& mesh() const + { + return fv::ddtScheme::mesh(); + } + + tmp > fvcDdt + ( + const dimensioned& + ); + + tmp > fvcDdt + ( + const GeometricField& + ); + + tmp > fvcDdt + ( + const dimensionedScalar&, + const GeometricField& + ); + + tmp > fvcDdt + ( + const volScalarField&, + const GeometricField& + ); + + tmp > fvmDdt + ( + const GeometricField& + ); + + tmp > fvmDdt + ( + const dimensionedScalar&, + const GeometricField& + ); + + tmp > fvmDdt + ( + const volScalarField&, + const GeometricField& + ); + + typedef typename ddtScheme::fluxFieldType fluxFieldType; + + tmp fvcDdtPhiCorr + ( + const volScalarField& rA, + const GeometricField& U, + const fluxFieldType& phi + ); + + tmp fvcDdtPhiCorr + ( + const volScalarField& rA, + const volScalarField& rho, + const GeometricField& U, + const fluxFieldType& phi + ); + + tmp meshPhi + ( + const GeometricField& + ); +}; + + +template<> +tmp boundedDdtScheme::fvcDdtPhiCorr +( + const volScalarField& rA, + const volScalarField& U, + const surfaceScalarField& phi +); + + +template<> +tmp boundedDdtScheme::fvcDdtPhiCorr +( + const volScalarField& rA, + const volScalarField& rho, + const volScalarField& U, + const surfaceScalarField& phi +); + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace fv + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "boundedDdtScheme.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/solvers/lagrangian/LTSReactingParcelFoam/rhoEqn.H b/src/finiteVolume/finiteVolume/ddtSchemes/boundedDdtScheme/boundedDdtSchemes.C similarity index 73% rename from applications/solvers/lagrangian/LTSReactingParcelFoam/rhoEqn.H rename to src/finiteVolume/finiteVolume/ddtSchemes/boundedDdtScheme/boundedDdtSchemes.C index 88b7a83221..8971a237e5 100644 --- a/applications/solvers/lagrangian/LTSReactingParcelFoam/rhoEqn.H +++ b/src/finiteVolume/finiteVolume/ddtSchemes/boundedDdtScheme/boundedDdtSchemes.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -21,30 +21,19 @@ License You should have received a copy of the GNU General Public License along with OpenFOAM. If not, see . -Global - rhoEqn - -Description - Solve the continuity for density. - \*---------------------------------------------------------------------------*/ +#include "boundedDdtScheme.H" +#include "fvMesh.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam { - fvScalarMatrix rhoEqn - ( - fvm::ddt(rho) - + fvc::div(phi) - == - parcels.Srho(rho) - + sources(rho) - ); - - sources.constrain(rhoEqn); - - rhoEqn.solve(); - - Info<< "rho min/max = " << min(rho).value() << ", " << max(rho).value() - << endl; +namespace fv +{ + makeFvDdtScheme(boundedDdtScheme) +} } // ************************************************************************* // diff --git a/src/fvMotionSolver/fvMotionSolvers/displacement/SBRStress/displacementSBRStressFvMotionSolver.H b/src/fvMotionSolver/fvMotionSolvers/displacement/SBRStress/displacementSBRStressFvMotionSolver.H index db98e1a7cc..6d31eec067 100644 --- a/src/fvMotionSolver/fvMotionSolvers/displacement/SBRStress/displacementSBRStressFvMotionSolver.H +++ b/src/fvMotionSolver/fvMotionSolvers/displacement/SBRStress/displacementSBRStressFvMotionSolver.H @@ -53,7 +53,6 @@ class motionDiffusivity; class displacementSBRStressFvMotionSolver : -// public displacementFvMotionSolver public displacementMotionSolver, public fvMotionSolverCore { diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C index 58725470e1..db246ba1f2 100644 --- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C +++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C @@ -99,76 +99,6 @@ void Foam::autoLayerDriver::sumWeights // Smooth field on moving patch -//void Foam::autoLayerDriver::smoothField -//( -// const motionSmoother& meshMover, -// const PackedBoolList& isMasterEdge, -// const labelList& meshEdges, -// const scalarField& fieldMin, -// const label nSmoothDisp, -// scalarField& field -//) const -//{ -// const indirectPrimitivePatch& pp = meshMover.patch(); -// const edgeList& edges = pp.edges(); -// const labelList& meshPoints = pp.meshPoints(); -// -// scalarField invSumWeight(pp.nPoints()); -// sumWeights -// ( -// isMasterEdge, -// meshEdges, -// meshPoints, -// edges, -// invSumWeight -// ); -// -// // Get smoothly varying patch field. -// Info<< "shrinkMeshDistance : Smoothing field ..." << endl; -// -// for (label iter = 0; iter < nSmoothDisp; iter++) -// { -// scalarField average(pp.nPoints()); -// averageNeighbours -// ( -// meshMover.mesh(), -// isMasterEdge, -// meshEdges, -// meshPoints, -// pp.edges(), -// invSumWeight, -// field, -// average -// ); -// -// // Transfer to field -// forAll(field, pointI) -// { -// //full smoothing neighbours + point value -// average[pointI] = 0.5*(field[pointI]+average[pointI]); -// -// // perform monotonic smoothing -// if -// ( -// average[pointI] < field[pointI] -// && average[pointI] >= fieldMin[pointI] -// ) -// { -// field[pointI] = average[pointI]; -// } -// } -// -// // Do residual calculation every so often. -// if ((iter % 10) == 0) -// { -// Info<< " Iteration " << iter << " residual " -// << gSum(mag(field-average)) -// /returnReduce(average.size(), sumOp