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/compressible/rhoPimpleFoam/EEqn.H b/applications/solvers/compressible/rhoPimpleFoam/EEqn.H new file mode 100644 index 0000000000..3ad3c0f03e --- /dev/null +++ b/applications/solvers/compressible/rhoPimpleFoam/EEqn.H @@ -0,0 +1,20 @@ +{ + 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(turbulence->alphaEff(), he) + ); + + EEqn.relax(); + EEqn.solve(); + + thermo.correct(); +} 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/EEqn.H b/applications/solvers/compressible/rhoSimpleFoam/EEqn.H new file mode 100644 index 0000000000..e4c79b40cc --- /dev/null +++ b/applications/solvers/compressible/rhoSimpleFoam/EEqn.H @@ -0,0 +1,19 @@ +{ + 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(turbulence->alphaEff(), he) + ); + + EEqn.relax(); + EEqn.solve(); + + thermo.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/compressible/rhoSimpleFoam/createFields.H b/applications/solvers/compressible/rhoSimpleFoam/createFields.H index d8fdff55ae..dea35b7657 100644 --- a/applications/solvers/compressible/rhoSimpleFoam/createFields.H +++ b/applications/solvers/compressible/rhoSimpleFoam/createFields.H @@ -5,6 +5,7 @@ psiThermo::New(mesh) ); psiThermo& thermo = pThermo(); + thermo.validate(args.executable(), "h", "e"); volScalarField rho ( @@ -20,7 +21,6 @@ ); volScalarField& p = thermo.p(); - volScalarField& e = thermo.he(); const volScalarField& psi = thermo.psi(); Info<< "Reading field U\n" << endl; diff --git a/applications/solvers/compressible/rhoSimpleFoam/eEqn.H b/applications/solvers/compressible/rhoSimpleFoam/eEqn.H deleted file mode 100644 index a1ea771573..0000000000 --- a/applications/solvers/compressible/rhoSimpleFoam/eEqn.H +++ /dev/null @@ -1,18 +0,0 @@ -{ - // Kinetic + pressure energy - volScalarField Ekp("Ekp", 0.5*magSqr(U) + p/rho); - - fvScalarMatrix eEqn - ( - fvm::div(phi, e) - - fvm::Sp(fvc::div(phi), e) - - fvm::laplacian(turbulence->alphaEff(), e) - == - fvc::div(phi)*Ekp - fvc::div(phi, Ekp) - ); - - eEqn.relax(); - eEqn.solve(); - - thermo.correct(); -} diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/EEqn.H b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/EEqn.H new file mode 100644 index 0000000000..ff467c0382 --- /dev/null +++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/EEqn.H @@ -0,0 +1,21 @@ +{ + 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(turbulence->alphaEff(), he) + ); + + pZones.addEnergySource(thermo, rho, EEqn); + + EEqn.relax(); + EEqn.solve(); + + thermo.correct(); +} diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/createFields.H b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/createFields.H index 96bdb2a122..4fff74d224 100644 --- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/createFields.H +++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/createFields.H @@ -5,6 +5,7 @@ rhoThermo::New(mesh) ); rhoThermo& thermo = pThermo(); + thermo.validate(args.executable(), "h", "e"); volScalarField rho ( @@ -20,7 +21,6 @@ ); volScalarField& p = thermo.p(); - volScalarField& e = thermo.he(); Info<< "Reading field U\n" << endl; volVectorField U diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/eEqn.H b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/eEqn.H deleted file mode 100644 index 5d0f174623..0000000000 --- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/eEqn.H +++ /dev/null @@ -1,20 +0,0 @@ -{ - // Kinetic + pressure energy - volScalarField Ekp("Ekp", 0.5*magSqr(U) + p/rho); - - fvScalarMatrix eEqn - ( - fvm::div(phi, e) - - fvm::Sp(fvc::div(phi), e) - - fvm::laplacian(turbulence->alphaEff(), e) - == - fvc::div(phi)*Ekp - fvc::div(phi, Ekp) - ); - - pZones.addEnergySource(thermo, rho, eEqn); - - eEqn.relax(); - eEqn.solve(); - - thermo.correct(); -} diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/rhoPorousMRFSimpleFoam.C b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/rhoPorousMRFSimpleFoam.C index 22a58da1cc..3898295b91 100644 --- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/rhoPorousMRFSimpleFoam.C +++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/rhoPorousMRFSimpleFoam.C @@ -63,7 +63,7 @@ int main(int argc, char *argv[]) // Pressure-velocity SIMPLE corrector { #include "UEqn.H" - #include "eEqn.H" + #include "EEqn.H" #include "pEqn.H" } diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C b/applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C index 0eee9129b5..312196583e 100644 --- a/applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C +++ b/applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C @@ -59,7 +59,7 @@ int main(int argc, char *argv[]) // Pressure-velocity SIMPLE corrector { #include "UEqn.H" - #include "eEqn.H" + #include "EEqn.H" #include "pEqn.H" } diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/rhoSimplecFoam.C b/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/rhoSimplecFoam.C index 561c5e2dc9..f707570f60 100644 --- a/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/rhoSimplecFoam.C +++ b/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/rhoSimplecFoam.C @@ -61,7 +61,7 @@ int main(int argc, char *argv[]) // Velocity-pressure-enthalpy SIMPLEC corrector { #include "UEqn.H" - #include "eEqn.H" + #include "EEqn.H" #include "pEqn.H" } diff --git a/applications/solvers/compressible/sonicFoam/EEqn.H b/applications/solvers/compressible/sonicFoam/EEqn.H new file mode 100644 index 0000000000..15979381a9 --- /dev/null +++ b/applications/solvers/compressible/sonicFoam/EEqn.H @@ -0,0 +1,10 @@ +{ + solve + ( + fvm::ddt(rho, e) + fvm::div(phi, e) + + fvc::ddt(rho, K) + fvc::div(phi, volScalarField("Ekp", K + p/rho)) + - fvm::laplacian(turbulence->alphaEff(), e) + ); + + thermo.correct(); +} diff --git a/applications/solvers/compressible/sonicFoam/createFields.H b/applications/solvers/compressible/sonicFoam/createFields.H index 3c39d3e7aa..c3421369e7 100644 --- a/applications/solvers/compressible/sonicFoam/createFields.H +++ b/applications/solvers/compressible/sonicFoam/createFields.H @@ -5,6 +5,7 @@ psiThermo::New(mesh) ); psiThermo& thermo = pThermo(); + thermo.validate(args.executable(), "e"); volScalarField& p = thermo.p(); volScalarField& e = thermo.he(); diff --git a/applications/solvers/compressible/sonicFoam/eEqn.H b/applications/solvers/compressible/sonicFoam/eEqn.H deleted file mode 100644 index cfd908ded8..0000000000 --- a/applications/solvers/compressible/sonicFoam/eEqn.H +++ /dev/null @@ -1,12 +0,0 @@ -{ - solve - ( - fvm::ddt(rho, e) - + fvm::div(phi, e) - - fvm::laplacian(turbulence->alphaEff(), e) - == - - (fvc::ddt(rho, K) + fvc::div(phi, volScalarField("Ekp", K + p/rho))) - ); - - thermo.correct(); -} diff --git a/applications/solvers/compressible/sonicFoam/pEqn.H b/applications/solvers/compressible/sonicFoam/pEqn.H index 99ca26fce3..4d700b3a41 100644 --- a/applications/solvers/compressible/sonicFoam/pEqn.H +++ b/applications/solvers/compressible/sonicFoam/pEqn.H @@ -17,7 +17,8 @@ surfaceScalarField phid volScalarField Dp("Dp", rho*rAU); -for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) +// Non-orthogonal pressure corrector loop +while (pimple.correctNonOrthogonal()) { fvScalarMatrix pEqn ( @@ -28,7 +29,7 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) pEqn.solve(); - if (nonOrth == nNonOrthCorr) + if (pimple.finalNonOrthogonalIter()) { phi = pEqn.flux(); } diff --git a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/EEqn.H b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/EEqn.H new file mode 100644 index 0000000000..994c2862b8 --- /dev/null +++ b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/EEqn.H @@ -0,0 +1,11 @@ +{ + solve + ( + fvm::ddt(rho, e) + fvm::div(phi, e) + + fvc::ddt(rho, K) + fvc::div(phi, K) + + fvc::div(phi/fvc::interpolate(rho) + mesh.phi(), p, "div(phiv,p)") + - fvm::laplacian(turbulence->alphaEff(), e) + ); + + thermo.correct(); +} diff --git a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/eEqn.H b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/eEqn.H deleted file mode 100644 index 0c17353ec4..0000000000 --- a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/eEqn.H +++ /dev/null @@ -1,12 +0,0 @@ -{ - solve - ( - fvm::ddt(rho, e) - + fvm::div(phi, e) - - fvm::laplacian(turbulence->alphaEff(), e) - == - - p*fvc::div(phi/fvc::interpolate(rho) + mesh.phi()) - ); - - thermo.correct(); -} diff --git a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/sonicDyMFoam.C b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/sonicDyMFoam.C index 163baa5b69..2f7ac9962a 100644 --- a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/sonicDyMFoam.C +++ b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/sonicDyMFoam.C @@ -34,6 +34,7 @@ Description #include "psiThermo.H" #include "turbulenceModel.H" #include "motionSolver.H" +#include "pimpleControl.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -45,6 +46,8 @@ int main(int argc, char *argv[]) #include "createFields.H" #include "initContinuityErrs.H" + pimpleControl pimple(mesh); + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Info<< "\nStarting time loop\n" << endl; @@ -62,19 +65,23 @@ int main(int argc, char *argv[]) #include "rhoEqn.H" - #include "UEqn.H" - - #include "eEqn.H" - - - // --- PISO loop - - for (int corr=0; corrcorrect(); + // --- Pressure corrector loop + while (pimple.correct()) + { + #include "pEqn.H" + } + + if (pimple.turbCorr()) + { + turbulence->correct(); + } + } rho = thermo.rho(); diff --git a/applications/solvers/compressible/sonicFoam/sonicFoam.C b/applications/solvers/compressible/sonicFoam/sonicFoam.C index 954404b6a4..71d032d422 100644 --- a/applications/solvers/compressible/sonicFoam/sonicFoam.C +++ b/applications/solvers/compressible/sonicFoam/sonicFoam.C @@ -33,6 +33,7 @@ Description #include "fvCFD.H" #include "psiThermo.H" #include "turbulenceModel.H" +#include "pimpleControl.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -44,6 +45,8 @@ int main(int argc, char *argv[]) #include "createFields.H" #include "initContinuityErrs.H" + pimpleControl pimple(mesh); + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Info<< "\nStarting time loop\n" << endl; @@ -52,20 +55,28 @@ int main(int argc, char *argv[]) { Info<< "Time = " << runTime.timeName() << nl << endl; - #include "readPISOControls.H" + #include "readTimeControls.H" #include "compressibleCourantNo.H" #include "rhoEqn.H" - #include "UEqn.H" - // --- PISO loop - for (int corr=0; corrcorrect(); + // --- Pressure corrector loop + while (pimple.correct()) + { + #include "pEqn.H" + } + + if (pimple.turbCorr()) + { + turbulence->correct(); + } + } rho = thermo.rho(); diff --git a/applications/solvers/compressible/sonicFoam/sonicLiquidFoam/sonicLiquidFoam.C b/applications/solvers/compressible/sonicFoam/sonicLiquidFoam/sonicLiquidFoam.C index 19a530a6aa..1ac1ae4ae3 100644 --- a/applications/solvers/compressible/sonicFoam/sonicLiquidFoam/sonicLiquidFoam.C +++ b/applications/solvers/compressible/sonicFoam/sonicLiquidFoam/sonicLiquidFoam.C @@ -31,6 +31,7 @@ Description \*---------------------------------------------------------------------------*/ #include "fvCFD.H" +#include "pimpleControl.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -44,6 +45,8 @@ int main(int argc, char *argv[]) #include "createFields.H" #include "initContinuityErrs.H" + pimpleControl pimple(mesh); + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Info<< "\nStarting time loop\n" << endl; @@ -52,57 +55,59 @@ int main(int argc, char *argv[]) { Info<< "Time = " << runTime.timeName() << nl << endl; - #include "readPISOControls.H" + #include "readTimeControls.H" #include "compressibleCourantNo.H" #include "rhoEqn.H" - fvVectorMatrix UEqn - ( - fvm::ddt(rho, U) - + fvm::div(phi, U) - - fvm::laplacian(mu, U) - ); - - solve(UEqn == -fvc::grad(p)); - - - // --- PISO loop - - for (int corr=0; corralphaEff(), h) - == - dpdt - - (fvc::ddt(rho, K) + fvc::div(phi, K)) - ); - - hEqn.relax(); - hEqn.solve(); - - thermo.correct(); -} diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/Allwmake b/applications/solvers/heatTransfer/buoyantSimpleFoam/Allwmake new file mode 100755 index 0000000000..0fe8e8f4ad --- /dev/null +++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/Allwmake @@ -0,0 +1,8 @@ +#!/bin/sh +cd ${0%/*} || exit 1 # run from this directory +set -x + +wmake +wmake buoyantSimpleRadiationFoam + +# ----------------------------------------------------------------- end-of-file diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/EEqn.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/EEqn.H new file mode 100644 index 0000000000..e4c79b40cc --- /dev/null +++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/EEqn.H @@ -0,0 +1,19 @@ +{ + 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(turbulence->alphaEff(), he) + ); + + EEqn.relax(); + EEqn.solve(); + + thermo.correct(); +} diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C b/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C index 2a355a21de..4a9387c188 100644 --- a/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C +++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C @@ -59,7 +59,7 @@ int main(int argc, char *argv[]) // Pressure-velocity SIMPLE corrector { #include "UEqn.H" - #include "hEqn.H" + #include "EEqn.H" #include "pEqn.H" } diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleRadiationFoam/EEqn.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleRadiationFoam/EEqn.H new file mode 100644 index 0000000000..0d1f41d0c3 --- /dev/null +++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleRadiationFoam/EEqn.H @@ -0,0 +1,22 @@ +{ + 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(turbulence->alphaEff(), he) + == + radiation->Sh(thermo) + ); + + EEqn.relax(); + EEqn.solve(); + + thermo.correct(); + radiation->correct(); +} diff --git a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/Make/files b/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleRadiationFoam/Make/files similarity index 98% rename from applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/Make/files rename to applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleRadiationFoam/Make/files index 79ebcd90aa..a8347c5525 100644 --- a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/Make/files +++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleRadiationFoam/Make/files @@ -1,4 +1,3 @@ buoyantSimpleRadiationFoam.C EXE = $(FOAM_APPBIN)/buoyantSimpleRadiationFoam - diff --git a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/Make/options b/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleRadiationFoam/Make/options similarity index 94% rename from applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/Make/options rename to applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleRadiationFoam/Make/options index f2ebc095ff..f26046adb2 100644 --- a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/Make/options +++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleRadiationFoam/Make/options @@ -1,5 +1,5 @@ EXE_INC = \ - -I../buoyantSimpleFoam \ + -I.. \ -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/radiationModels/lnInclude \ -I$(LIB_SRC)/turbulenceModels \ diff --git a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C b/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C similarity index 98% rename from applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C rename to applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C index 07718aa41d..d02853591e 100644 --- a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C +++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C @@ -62,7 +62,7 @@ int main(int argc, char *argv[]) // Pressure-velocity SIMPLE corrector { #include "UEqn.H" - #include "hEqn.H" + #include "EEqn.H" #include "pEqn.H" } diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H index d2fff48837..bb7a65cb1d 100644 --- a/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H +++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H @@ -5,6 +5,7 @@ psiThermo::New(mesh) ); psiThermo& thermo = pThermo(); + thermo.validate(args.executable(), "h", "e"); volScalarField rho ( @@ -20,7 +21,6 @@ ); volScalarField& p = thermo.p(); - volScalarField& h = thermo.he(); const volScalarField& psi = thermo.psi(); Info<< "Reading field U\n" << endl; diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/hEqn.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/hEqn.H deleted file mode 100644 index 24ed135c08..0000000000 --- a/applications/solvers/heatTransfer/buoyantSimpleFoam/hEqn.H +++ /dev/null @@ -1,15 +0,0 @@ -{ - fvScalarMatrix hEqn - ( - fvm::div(phi, h) - - fvm::Sp(fvc::div(phi), h) - - fvm::laplacian(turbulence->alphaEff(), h) - == - - fvc::div(phi, 0.5*magSqr(U), "div(phi,K)") - ); - - hEqn.relax(); - hEqn.solve(); - - thermo.correct(); -} diff --git a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/hEqn.H b/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/hEqn.H deleted file mode 100644 index 2c422a5cab..0000000000 --- a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/hEqn.H +++ /dev/null @@ -1,19 +0,0 @@ -{ - fvScalarMatrix hEqn - ( - fvm::div(phi, h) - - fvm::Sp(fvc::div(phi), h) - - fvm::laplacian(turbulence->alphaEff(), h) - == - - fvc::div(phi, 0.5*magSqr(U), "div(phi,K)") - + radiation->Sh(thermo) - ); - - hEqn.relax(); - - hEqn.solve(); - - thermo.correct(); - - radiation->correct(); -} 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/porousExplicitSourceReactingParcelFoam/readAdditionalSolutionControls.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/readAdditionalSolutionControls.H deleted file mode 100644 index 050c24fe7a..0000000000 --- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/readAdditionalSolutionControls.H +++ /dev/null @@ -1,8 +0,0 @@ -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/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 96% rename from applications/solvers/lagrangian/LTSReactingParcelFoam/LTSReactingParcelFoam.C rename to applications/solvers/lagrangian/reactingParcelFoam/LTSReactingParcelFoam/LTSReactingParcelFoam.C index 6f41f328a3..fabb0bba60 100644 --- a/applications/solvers/lagrangian/LTSReactingParcelFoam/LTSReactingParcelFoam.C +++ b/applications/solvers/lagrangian/reactingParcelFoam/LTSReactingParcelFoam/LTSReactingParcelFoam.C @@ -59,8 +59,8 @@ int main(int argc, char *argv[]) pimpleControl pimple(mesh); #include "readTimeControls.H" - #include "readAdditionalSolutionControls.H" #include "createFields.H" + #include "createRDeltaT.H" #include "createRadiationModel.H" #include "createClouds.H" #include "createExplicitSources.H" @@ -73,7 +73,6 @@ int main(int argc, char *argv[]) while (runTime.run()) { - #include "readAdditionalSolutionControls.H" #include "readTimeControls.H" runTime++; @@ -93,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/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..a22267d9be 100644 --- a/applications/solvers/lagrangian/reactingParcelFoam/YEqn.H +++ b/applications/solvers/lagrangian/reactingParcelFoam/YEqn.H @@ -5,11 +5,10 @@ tmp > mvConvection mesh, fields, phi, - mesh.divScheme("div(phi,Yi_hs)") + mesh.divScheme("div(phi,Yi_h)") ) ); - { combustion->correct(); dQ = combustion->dQ(); @@ -22,22 +21,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..8ef6838c59 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; @@ -81,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/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..1b4f46289f 100644 --- a/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/sprayEngineFoam.C +++ b/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/sprayEngineFoam.C @@ -86,7 +86,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..02739bda71 100644 --- a/applications/solvers/lagrangian/sprayFoam/sprayFoam.C +++ b/applications/solvers/lagrangian/sprayFoam/sprayFoam.C @@ -81,7 +81,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/config/settings.csh b/etc/config/settings.csh index 0896420bcc..069624f513 100644 --- a/etc/config/settings.csh +++ b/etc/config/settings.csh @@ -151,15 +151,17 @@ setenv FOAM_LIBBIN $WM_PROJECT_DIR/platforms/$WM_OPTIONS/lib # external (ThirdParty) libraries setenv FOAM_EXT_LIBBIN $WM_THIRD_PARTY_DIR/platforms/$WM_OPTIONS/lib +# site-specific directory +if ( $?WM_PROJECT_SITE ) then + set siteDir=$WM_PROJECT_SITE +else + set siteDir=$WM_PROJECT_INST_DIR/site +endif + # shared site executables/libraries # similar naming convention as ~OpenFOAM expansion -if ( $?WM_PROJECT_SITE ) then - setenv FOAM_SITE_APPBIN $WM_PROJECT_SITE/$WM_PROJECT_VERSION/platforms/$WM_OPTIONS/bin - setenv FOAM_SITE_LIBBIN $WM_PROJECT_SITE/$WM_PROJECT_VERSION/platforms/$WM_OPTIONS/lib -else - setenv FOAM_SITE_APPBIN $WM_PROJECT_INST_DIR/site/$WM_PROJECT_VERSION/platforms/$WM_OPTIONS/bin - setenv FOAM_SITE_LIBBIN $WM_PROJECT_INST_DIR/site/$WM_PROJECT_VERSION/platforms/$WM_OPTIONS/lib -endif +setenv FOAM_SITE_APPBIN $siteDir/$WM_PROJECT_VERSION/platforms/$WM_OPTIONS/bin +setenv FOAM_SITE_LIBBIN $siteDir/$WM_PROJECT_VERSION/platforms/$WM_OPTIONS/lib # user executables/libraries setenv FOAM_USER_APPBIN $WM_PROJECT_USER_DIR/platforms/$WM_OPTIONS/bin @@ -182,6 +184,15 @@ if ( -d "${WM_DIR}" ) setenv PATH ${WM_DIR}:${PATH} # add OpenFOAM scripts to the path setenv PATH ${WM_PROJECT_DIR}/bin:${PATH} +# add site-specific scripts to path - only if they exist +if ( -d "$siteDir/bin" ) then # generic + _foamAddPath "$siteDir/bin" +endif +if ( -d "$siteDir/$WM_PROJECT_VERSION/bin" ) then # version-specific + _foamAddPath "$siteDir/$WM_PROJECT_VERSION/bin" +endif +unset siteDir + _foamAddPath ${FOAM_USER_APPBIN}:${FOAM_SITE_APPBIN}:${FOAM_APPBIN} # Make sure to pick up dummy versions of external libraries last _foamAddLib ${FOAM_USER_LIBBIN}:${FOAM_SITE_LIBBIN}:${FOAM_LIBBIN}:${FOAM_EXT_LIBBIN}:${FOAM_LIBBIN}/dummy diff --git a/etc/config/settings.sh b/etc/config/settings.sh index 1a7e8237f9..1fc2136530 100644 --- a/etc/config/settings.sh +++ b/etc/config/settings.sh @@ -179,16 +179,13 @@ export FOAM_LIBBIN=$WM_PROJECT_DIR/platforms/$WM_OPTIONS/lib # external (ThirdParty) libraries export FOAM_EXT_LIBBIN=$WM_THIRD_PARTY_DIR/platforms/$WM_OPTIONS/lib +# site-specific directory +siteDir="${WM_PROJECT_SITE:-$WM_PROJECT_INST_DIR/site}" + # shared site executables/libraries # similar naming convention as ~OpenFOAM expansion -if [ -n "$WM_PROJECT_SITE" ] -then - export FOAM_SITE_APPBIN=$WM_PROJECT_SITE/$WM_PROJECT_VERSION/platforms/$WM_OPTIONS/bin - export FOAM_SITE_LIBBIN=$WM_PROJECT_SITE/$WM_PROJECT_VERSION/platforms/$WM_OPTIONS/lib -else - export FOAM_SITE_APPBIN=$WM_PROJECT_INST_DIR/site/$WM_PROJECT_VERSION/platforms/$WM_OPTIONS/bin - export FOAM_SITE_LIBBIN=$WM_PROJECT_INST_DIR/site/$WM_PROJECT_VERSION/platforms/$WM_OPTIONS/lib -fi +export FOAM_SITE_APPBIN=$siteDir/$WM_PROJECT_VERSION/platforms/$WM_OPTIONS/bin +export FOAM_SITE_LIBBIN=$siteDir/$WM_PROJECT_VERSION/platforms/$WM_OPTIONS/lib # user executables/libraries export FOAM_USER_APPBIN=$WM_PROJECT_USER_DIR/platforms/$WM_OPTIONS/bin @@ -211,6 +208,17 @@ export FOAM_RUN=$WM_PROJECT_USER_DIR/run # add OpenFOAM scripts to the path export PATH=$WM_PROJECT_DIR/bin:$PATH +# add site-specific scripts to path - only if they exist +if [ -d "$siteDir/bin" ] # generic +then + _foamAddPath "$siteDir/bin" +fi +if [ -d "$siteDir/$WM_PROJECT_VERSION/bin" ] # version-specific +then + _foamAddPath "$siteDir/$WM_PROJECT_VERSION/bin" +fi +unset siteDir + _foamAddPath $FOAM_USER_APPBIN:$FOAM_SITE_APPBIN:$FOAM_APPBIN # Make sure to pick up dummy versions of external libraries last _foamAddLib $FOAM_USER_LIBBIN:$FOAM_SITE_LIBBIN:$FOAM_LIBBIN:$FOAM_EXT_LIBBIN:$FOAM_LIBBIN/dummy @@ -232,6 +240,12 @@ fi case "${foamCompiler}" in OpenFOAM | ThirdParty) case "$WM_COMPILER" in + Gcc463) + gcc_version=gcc-4.6.3 + gmp_version=gmp-5.0.2 + mpfr_version=mpfr-3.0.1 + mpc_version=mpc-0.9 + ;; Gcc | Gcc++0x | Gcc46 | Gcc46++0x) gcc_version=gcc-4.6.1 gmp_version=gmp-5.0.4 diff --git a/etc/controlDict b/etc/controlDict index 04ec43aca0..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; @@ -884,6 +884,7 @@ DebugSwitches wallHeatTransfer 0; wallLayerCells 0; wallModel 0; + warnUnboundedGauss 1; waveTransmissive 0; wedge 0; weighted 0; @@ -964,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/Allwmake b/src/Allwmake index a61c7a2e3e..94afb0d985 100755 --- a/src/Allwmake +++ b/src/Allwmake @@ -67,7 +67,6 @@ turbulenceModels/Allwmake $* wmake $makeType combustionModels regionModels/Allwmake $* lagrangian/Allwmake $* -postProcessing/Allwmake $* mesh/Allwmake $* fvAgglomerationMethods/Allwmake $* @@ -77,4 +76,6 @@ wmake $makeType engine wmake $makeType fieldSources +postProcessing/Allwmake $* + # ----------------------------------------------------------------- end-of-file 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/db/dictionary/dictionary.C b/src/OpenFOAM/db/dictionary/dictionary.C index 5e67648d17..40abb17d47 100644 --- a/src/OpenFOAM/db/dictionary/dictionary.C +++ b/src/OpenFOAM/db/dictionary/dictionary.C @@ -465,7 +465,8 @@ const Foam::entry* Foam::dictionary::lookupScopedEntryPtr *this ) << "keyword " << keyword << " is undefined in dictionary " - << name() + << name() << endl + << "Valid keywords are " << keys() << exit(FatalIOError); } if (!entPtr->isDict()) 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 1acc4bcf01..7a354e0c64 100644 --- a/src/combustionModels/singleStepCombustion/singleStepCombustion.C +++ b/src/combustionModels/singleStepCombustion/singleStepCombustion.C @@ -118,7 +118,7 @@ tmp 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/fieldSources/basicSource/basicSource/basicSourceList.C b/src/fieldSources/basicSource/basicSource/basicSourceList.C index 819d7b8f6b..2d484868a4 100644 --- a/src/fieldSources/basicSource/basicSource/basicSourceList.C +++ b/src/fieldSources/basicSource/basicSource/basicSourceList.C @@ -62,6 +62,22 @@ Foam::basicSourceList::basicSourceList PtrList(), mesh_(mesh), checkTimeIndex_(mesh_.time().startTimeIndex() + 2) +{ + reset(dict); +} + + +Foam::basicSourceList::basicSourceList(const fvMesh& mesh) +: + PtrList(), + mesh_(mesh), + checkTimeIndex_(mesh_.time().startTimeIndex() + 2) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void Foam::basicSourceList::reset(const dictionary& dict) { label count = 0; forAllConstIter(dictionary, dict, iter) @@ -85,15 +101,13 @@ Foam::basicSourceList::basicSourceList this->set ( i++, - basicSource::New(name, sourceDict, mesh) + basicSource::New(name, sourceDict, mesh_) ); } } } -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - bool Foam::basicSourceList::read(const dictionary& dict) { checkTimeIndex_ = mesh_.time().timeIndex() + 2; diff --git a/src/fieldSources/basicSource/basicSource/basicSourceList.H b/src/fieldSources/basicSource/basicSource/basicSourceList.H index 7d77de55d1..6a7a0126ed 100644 --- a/src/fieldSources/basicSource/basicSource/basicSourceList.H +++ b/src/fieldSources/basicSource/basicSource/basicSourceList.H @@ -84,7 +84,10 @@ public: // Constructors - //- Construct from components with list of field names + //- Construct null + basicSourceList(const fvMesh& mesh); + + //- Construct from mesh and dictionary basicSourceList(const fvMesh& mesh, const dictionary& dict); @@ -95,6 +98,9 @@ public: // Member Functions + //- Reset the source list + void reset(const dictionary& dict); + //- Correct template void correct(GeometricField& fld); diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files index 22988f3216..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 @@ -345,6 +346,7 @@ convectionSchemes = finiteVolume/convectionSchemes $(convectionSchemes)/convectionScheme/convectionSchemes.C $(convectionSchemes)/gaussConvectionScheme/gaussConvectionSchemes.C $(convectionSchemes)/multivariateGaussConvectionScheme/multivariateGaussConvectionSchemes.C +$(convectionSchemes)/boundedConvectionScheme/boundedConvectionSchemes.C laplacianSchemes = finiteVolume/laplacianSchemes $(laplacianSchemes)/laplacianScheme/laplacianSchemes.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 0b78dddc22..8051261154 100644 --- a/src/finiteVolume/fields/fvPatchFields/derived/flowRateInletVelocity/flowRateInletVelocityFvPatchVectorField.H +++ b/src/finiteVolume/fields/fvPatchFields/derived/flowRateInletVelocity/flowRateInletVelocityFvPatchVectorField.H @@ -32,13 +32,12 @@ Description from the flux (volumetric or mass-based), whose direction is assumed to be normal to the patch. - The basis of the patch (volumetric or mass) is determined by the - dimensions of the flux, \c phi, and the value of the density field name, - \c rhoName entry. - For a mass-based flux: - - if \c rhoName is a valid density field name, the flow rate is in kg/s + - the flow rate should be provided in kg/s - if \c rhoName is "none" the flow rate is in m3/s + - otherwise \c rhoName should correspond to the name of the density field + - if the density field cannot be found in the database, the user must + specify the inlet density using the \c rhoInlet entry For a volumetric-based flux: - the flow rate is in m3/s @@ -47,24 +46,38 @@ Description \table Property | Description | Required | Default value - flowRate | volumetric [m3/s] OR mass flow rate [kg/s] | yes | + massFlowRate | mass flow rate [kg/s] | no | + volumetricFlowRate | volumetric flow rate [m3/s]| no | + rhoInlet | inlet density | no | \endtable - Example of the boundary condition specification: + Example of the boundary condition specification for a volumetric flow rate: \verbatim myPatch { type flowRateInletVelocity; - flowRate 0.2; - rho rho; + volumetricFlowRate 0.2; value uniform (0 0 0); // placeholder } \endverbatim + Example of the boundary condition specification for a mass flow rate: + \verbatim + myPatch + { + type flowRateInletVelocity; + massFlowRate 0.2; + rho rho; + rhoInlet 1.0; + } + \endverbatim + The \c flowRate entry is a \c DataEntry type, meaning that it can be specified as constant, a polynomial fuction of time, and ... Note + - \c rhoInlet is required for the case of a mass flow rate, where the + density field is not available at start-up - the value is positive into the domain (as an inlet) - may not work correctly for transonic inlets - strange behaviour with potentialFoam since the U equation is not solved @@ -101,12 +114,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: @@ -179,6 +195,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/convectionSchemes/boundedConvectionScheme/boundedConvectionScheme.C b/src/finiteVolume/finiteVolume/convectionSchemes/boundedConvectionScheme/boundedConvectionScheme.C new file mode 100644 index 0000000000..1e5c911763 --- /dev/null +++ b/src/finiteVolume/finiteVolume/convectionSchemes/boundedConvectionScheme/boundedConvectionScheme.C @@ -0,0 +1,103 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 "boundedConvectionScheme.H" +#include "fvcSurfaceIntegrate.H" +#include "fvMatrices.H" +#include "fvmSup.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace fv +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +template +tmp > +boundedConvectionScheme::interpolate +( + const surfaceScalarField& phi, + const GeometricField& vf +) const +{ + return scheme_().interpolate(phi, vf); +} + + +template +tmp > +boundedConvectionScheme::flux +( + const surfaceScalarField& faceFlux, + const GeometricField& vf +) const +{ + return scheme_().flux(faceFlux, vf); +} + + +template +tmp > +boundedConvectionScheme::fvmDiv +( + const surfaceScalarField& faceFlux, + const GeometricField& vf +) const +{ + return + scheme_().fvmDiv(faceFlux, vf) + - fvm::Sp(fvc::surfaceIntegrate(faceFlux), vf); +} + + +template +tmp > +boundedConvectionScheme::fvcDiv +( + const surfaceScalarField& faceFlux, + const GeometricField& vf +) const +{ + return + scheme_().fvcDiv(faceFlux, vf) + - fvc::surfaceIntegrate(faceFlux)*vf; +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace fv + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/convectionSchemes/boundedConvectionScheme/boundedConvectionScheme.H b/src/finiteVolume/finiteVolume/convectionSchemes/boundedConvectionScheme/boundedConvectionScheme.H new file mode 100644 index 0000000000..9426104d33 --- /dev/null +++ b/src/finiteVolume/finiteVolume/convectionSchemes/boundedConvectionScheme/boundedConvectionScheme.H @@ -0,0 +1,150 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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::boundedConvectionScheme + +Description + Bounded form of the selected convection scheme. + + Boundedness is achieved by subtracting div(phi)*vf or Sp(div(phi), vf) + which is non-conservative if div(phi) != 0 but conservative otherwise. + + Can be used for convection of bounded scalar properties in steady-state + solvers to improve stability if insufficient convergence of the pressure + equation causes temporary divergence of the flux field. + +SourceFiles + boundedConvectionScheme.C + +\*---------------------------------------------------------------------------*/ + +#ifndef boundedConvectionScheme_H +#define boundedConvectionScheme_H + +#include "convectionScheme.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace fv +{ + +/*---------------------------------------------------------------------------*\ + Class boundedConvectionScheme Declaration +\*---------------------------------------------------------------------------*/ + +template +class boundedConvectionScheme +: + public fv::convectionScheme +{ + // Private data + + tmp > scheme_; + + + // Private Member Functions + + //- Disallow default bitwise copy construct + boundedConvectionScheme(const boundedConvectionScheme&); + + //- Disallow default bitwise assignment + void operator=(const boundedConvectionScheme&); + + +public: + + //- Runtime type information + TypeName("bounded"); + + + // Constructors + + //- Construct from flux and Istream + boundedConvectionScheme + ( + const fvMesh& mesh, + const surfaceScalarField& faceFlux, + Istream& is + ) + : + convectionScheme(mesh, faceFlux), + scheme_ + ( + fv::convectionScheme::New(mesh, faceFlux, is) + ) + {} + + + // Member Functions + + tmp > interpolate + ( + const surfaceScalarField&, + const GeometricField& + ) const; + + tmp > flux + ( + const surfaceScalarField&, + const GeometricField& + ) const; + + tmp > fvmDiv + ( + const surfaceScalarField&, + const GeometricField& + ) const; + + tmp > fvcDiv + ( + const surfaceScalarField&, + const GeometricField& + ) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace fv + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "boundedConvectionScheme.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/convectionSchemes/boundedConvectionScheme/boundedConvectionSchemes.C b/src/finiteVolume/finiteVolume/convectionSchemes/boundedConvectionScheme/boundedConvectionSchemes.C new file mode 100644 index 0000000000..33385d4ecf --- /dev/null +++ b/src/finiteVolume/finiteVolume/convectionSchemes/boundedConvectionScheme/boundedConvectionSchemes.C @@ -0,0 +1,39 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 "boundedConvectionScheme.H" +#include "fvMesh.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace fv +{ + makeFvConvectionScheme(boundedConvectionScheme) +} +} + +// ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/convectionSchemes/gaussConvectionScheme/gaussConvectionScheme.H b/src/finiteVolume/finiteVolume/convectionSchemes/gaussConvectionScheme/gaussConvectionScheme.H index 1b27e361a9..61d6c22d37 100644 --- a/src/finiteVolume/finiteVolume/convectionSchemes/gaussConvectionScheme/gaussConvectionScheme.H +++ b/src/finiteVolume/finiteVolume/convectionSchemes/gaussConvectionScheme/gaussConvectionScheme.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 @@ -47,6 +47,10 @@ namespace Foam namespace fv { +//- Temporary debug switch to provide warning about backward-compatibility +// issue with setting div schemes for steady-state +extern int warnUnboundedGauss; + /*---------------------------------------------------------------------------*\ Class gaussConvectionScheme Declaration \*---------------------------------------------------------------------------*/ @@ -103,7 +107,29 @@ public: ( surfaceInterpolationScheme::New(mesh, faceFlux, is) ) - {} + { + is.rewind(); + word bounded(is); + + if + ( + warnUnboundedGauss + && word(mesh.ddtScheme("default")) == "steadyState" + && bounded != "bounded" + ) + { + fileNameList controlDictFiles(findEtcFiles("controlDict")); + + IOWarningIn("gaussConvectionScheme", is) + << "Unbounded 'Gauss' div scheme used in " + "steady-state solver, use 'bounded Gauss' " + "to ensure boundedness.\n" + << " To remove this warning switch off " + << "'boundedGauss' in " + << controlDictFiles[controlDictFiles.size()-1] + << endl; + } + } // Member Functions diff --git a/src/finiteVolume/finiteVolume/convectionSchemes/gaussConvectionScheme/gaussConvectionSchemes.C b/src/finiteVolume/finiteVolume/convectionSchemes/gaussConvectionScheme/gaussConvectionSchemes.C index 55a5e6049d..c5ea8dbac7 100644 --- a/src/finiteVolume/finiteVolume/convectionSchemes/gaussConvectionScheme/gaussConvectionSchemes.C +++ b/src/finiteVolume/finiteVolume/convectionSchemes/gaussConvectionScheme/gaussConvectionSchemes.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,11 @@ namespace Foam { namespace fv { + int warnUnboundedGauss + ( + Foam::debug::debugSwitch("warnUnboundedGauss", true) + ); + makeFvConvectionScheme(gaussConvectionScheme) } } 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/src/finiteVolume/finiteVolume/ddtSchemes/boundedDdtScheme/boundedDdtSchemes.C b/src/finiteVolume/finiteVolume/ddtSchemes/boundedDdtScheme/boundedDdtSchemes.C new file mode 100644 index 0000000000..8971a237e5 --- /dev/null +++ b/src/finiteVolume/finiteVolume/ddtSchemes/boundedDdtScheme/boundedDdtSchemes.C @@ -0,0 +1,39 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 . + +\*---------------------------------------------------------------------------*/ + +#include "boundedDdtScheme.H" +#include "fvMesh.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace fv +{ + makeFvDdtScheme(boundedDdtScheme) +} +} + +// ************************************************************************* // diff --git a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C index 0cd882fe54..d83681cba1 100644 --- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C +++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C @@ -65,6 +65,15 @@ Foam::solverPerformance Foam::fvMatrix::solve << endl; } + label maxIter = -1; + if (solverControls.readIfPresent("maxIter", maxIter)) + { + if (maxIter == 0) + { + return solverPerformance(); + } + } + word type(solverControls.lookupOrDefault("type", "segregated")); if (type == "segregated") 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/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H index 321fba33ad..29e696dc81 100644 --- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H +++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H @@ -470,8 +470,8 @@ public: //- Total rotational kinetic energy in the system inline scalar rotationalKineticEnergyOfSystem() const; - //- Penetration for percentage of the current total mass - inline scalar penetration(const scalar& prc) const; + //- Penetration for fraction [0-1] of the current total mass + inline scalar penetration(const scalar& fraction) const; //- Mean diameter Dij inline scalar Dij(const label i, const label j) const; diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H index a6bdf47cc4..2e793ac71e 100644 --- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H +++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H @@ -24,6 +24,7 @@ License \*---------------------------------------------------------------------------*/ #include "fvmSup.H" +#include "SortableList.H" // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // @@ -336,106 +337,130 @@ inline Foam::scalar Foam::KinematicCloud::Dmax() const reduce(d, maxOp()); - return d; + return max(0.0, d); } template inline Foam::scalar Foam::KinematicCloud::penetration ( - const scalar& prc + const scalar& fraction ) const { - scalar distance = 0.0; - scalar mTot = 0.0; - - label np = this->size(); - - // arrays containing the parcels mass and - // distance from injector in ascending order - scalarField mass(np); - scalarField dist(np); - - if (np > 0) + if ((fraction < 0) || (fraction > 1)) { - label n = 0; - - // first arrange the parcels in ascending order - // the first parcel is closest to its injection position - // and the last one is most far away. - forAllConstIter(typename KinematicCloud, *this, iter) - { - const parcelType& p = iter(); - scalar mi = p.nParticle()*p.mass(); - scalar di = mag(p.position() - p.position0()); - mTot += mi; - - // insert at the last place - mass[n] = mi; - dist[n] = di; - - label i = 0; - bool found = false; - - // insert the parcel in the correct place - // and move the others - while ((i < n) && (!found)) - { - if (di < dist[i]) - { - found = true; - for (label j=n; j>i; j--) - { - mass[j] = mass[j-1]; - dist[j] = dist[j-1]; - } - mass[i] = mi; - dist[i] = di; - } - i++; - } - n++; - } + FatalErrorIn + ( + "inline Foam::scalar Foam::KinematicCloud::penetration" + "(" + "const scalar&" + ") const" + ) + << "fraction should be in the range 0 < fraction < 1" + << exit(FatalError); } - reduce(mTot, sumOp()); + scalar distance = 0.0; - if (np > 0) + const label nParcel = this->size(); + globalIndex globalParcels(nParcel); + const label nParcelSum = globalParcels.size(); + + if (nParcelSum == 0) { - scalar mLimit = prc*mTot; - scalar mOff = (1.0 - prc)*mTot; + return distance; + } - if (np > 1) + // lists of parcels mass and distance from initial injection point + List mass(nParcel, 0.0); + List dist(nParcel, 0.0); + + label i = 0; + scalar mSum = 0.0; + forAllConstIter(typename KinematicCloud, *this, iter) + { + const parcelType& p = iter(); + scalar m = p.nParticle()*p.mass(); + scalar d = mag(p.position() - p.position0()); + mSum += m; + + mass[i] = m; + dist[i] = d; + + i++; + } + + // calculate total mass across all processors + reduce(mSum, sumOp()); + + // flatten the mass list + List allMass(nParcelSum, 0.0); + SubList + ( + allMass, + globalParcels.localSize(Pstream::myProcNo()), + globalParcels.offset(Pstream::myProcNo()) + ).assign(mass); + + // flatten the distance list + SortableList allDist(nParcelSum, 0.0); + SubList + ( + allDist, + globalParcels.localSize(Pstream::myProcNo()), + globalParcels.offset(Pstream::myProcNo()) + ).assign(dist); + + // sort allDist distances into ascending order + // note: allMass masses are left unsorted + allDist.sort(); + + if (nParcelSum > 1) + { + const scalar mLimit = fraction*mSum; + const labelList& indices = allDist.indices(); + + if (mLimit > (mSum - allMass[indices.last()])) { - // 'prc' is large enough that the parcel most far - // away will be used, no need to loop... - if (mLimit > mTot - mass[np-1]) - { - distance = dist[np-1]; - } - else - { - scalar mOffSum = 0.0; - label i = np; - - while ((mOffSum < mOff) && (i>0)) - { - i--; - mOffSum += mass[i]; - } - distance = - dist[i+1] - + (dist[i] - dist[i+1])*(mOffSum - mOff) - /mass[i+1] ; - } + distance = allDist.last(); } else { - distance = dist[0]; + // assuming that 'fraction' is generally closer to 1 than 0, loop + // through in reverse distance order + const scalar mThreshold = (1.0 - fraction)*mSum; + scalar mCurrent = 0.0; + label i0 = 0; + + forAllReverse(indices, i) + { + label indI = indices[i]; + + mCurrent += allMass[indI]; + + if (mCurrent > mThreshold) + { + i0 = i; + break; + } + } + + if (i0 == indices.size() - 1) + { + distance = allDist.last(); + } + else + { + // linearly interpolate to determine distance + scalar alpha = (mCurrent - mThreshold)/allMass[indices[i0]]; + distance = allDist[i0] + alpha*(allDist[i0+1] - allDist[i0]); + } } } - - reduce(distance, maxOp()); + else + { + distance = allDist.first(); + } return distance; } diff --git a/src/lagrangian/intermediate/clouds/baseClasses/kinematicCloud/kinematicCloud.H b/src/lagrangian/intermediate/clouds/baseClasses/kinematicCloud/kinematicCloud.H index 01c80c16e2..d7a0fcf0d0 100644 --- a/src/lagrangian/intermediate/clouds/baseClasses/kinematicCloud/kinematicCloud.H +++ b/src/lagrangian/intermediate/clouds/baseClasses/kinematicCloud/kinematicCloud.H @@ -89,7 +89,7 @@ public: virtual scalar rotationalKineticEnergyOfSystem() const = 0; //- Penetration for percentage of the current total mass -// virtual scalar penetration(const scalar& prc) const = 0; +// virtual scalar penetration(const scalar& fraction) const = 0; //- Mean diameter Dij virtual scalar Dij(const label i, const label j) const = 0; diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C index ff2f3a6e84..32473b8acf 100644 --- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C +++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C @@ -585,16 +585,23 @@ void Foam::InjectionModel::inject(TrackData& td) pPtr->rho() ); - if ((pPtr->nParticle() >= 1.0) && (pPtr->move(td, dt))) + if (!pPtr->move(td, dt)) { - td.cloud().addParticle(pPtr); - massAdded += pPtr->nParticle()*pPtr->mass(); - parcelsAdded++; + delete pPtr; } else { - delayedVolume += pPtr->nParticle()*pPtr->volume(); - delete pPtr; + if (pPtr->nParticle() >= 1.0) + { + td.cloud().addParticle(pPtr); + massAdded += pPtr->nParticle()*pPtr->mass(); + parcelsAdded++; + } + else + { + delayedVolume += pPtr->nParticle()*pPtr->volume(); + delete pPtr; + } } } } diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/PressureGradient/PressureGradientForce.C b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/PressureGradient/PressureGradientForce.C index a1d121ee29..c9d61f6d34 100644 --- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/PressureGradient/PressureGradientForce.C +++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/PressureGradient/PressureGradientForce.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) 2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -24,6 +24,7 @@ License \*---------------------------------------------------------------------------*/ #include "PressureGradientForce.H" +#include "fvcDdt.H" #include "fvcGrad.H" // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // @@ -33,12 +34,13 @@ Foam::PressureGradientForce::PressureGradientForce ( CloudType& owner, const fvMesh& mesh, - const dictionary& dict + const dictionary& dict, + const word& forceType ) : - ParticleForce(owner, mesh, dict, typeName, true), - UName_(this->coeffs().lookup("U")), - gradUPtr_(NULL) + ParticleForce(owner, mesh, dict, forceType, true), + UName_(this->coeffs().template lookupOrDefault("U", "U")), + DUcDtInterpPtr_(NULL) {} @@ -50,7 +52,7 @@ Foam::PressureGradientForce::PressureGradientForce : ParticleForce(pgf), UName_(pgf.UName_), - gradUPtr_(NULL) + DUcDtInterpPtr_(NULL) {} @@ -66,18 +68,48 @@ Foam::PressureGradientForce::~PressureGradientForce() template void Foam::PressureGradientForce::cacheFields(const bool store) { + static word fName("DUcDt"); + + bool fieldExists = this->mesh().template foundObject(fName); + if (store) { - const volVectorField& U = this->mesh().template - lookupObject(UName_); - gradUPtr_ = fvc::grad(U).ptr(); + if (!fieldExists) + { + const volVectorField& Uc = this->mesh().template + lookupObject(UName_); + + volVectorField* DUcDtPtr = new volVectorField + ( + fName, + fvc::ddt(Uc) + (Uc & fvc::grad(Uc)) + ); + + DUcDtPtr->store(); + } + + const volVectorField& DUcDt = this->mesh().template + lookupObject(fName); + + DUcDtInterpPtr_.reset + ( + interpolation::New + ( + this->owner().solution().interpolationSchemes(), + DUcDt + ).ptr() + ); } else { - if (gradUPtr_) + DUcDtInterpPtr_.clear(); + + if (fieldExists) { - delete gradUPtr_; - gradUPtr_ = NULL; + const volVectorField& DUcDt = this->mesh().template + lookupObject(fName); + + const_cast(DUcDt).checkOut(); } } } @@ -95,11 +127,24 @@ Foam::forceSuSp Foam::PressureGradientForce::calcCoupled { forceSuSp value(vector::zero, 0.0); - const volTensorField& gradU = *gradUPtr_; - value.Su() = mass*p.rhoc()/p.rho()*(p.U() & gradU[p.cell()]); + vector DUcDt = + DUcDtInterp().interpolate(p.position(), p.currentTetIndices()); + + value.Su() = mass*p.rhoc()/p.rho()*DUcDt; return value; } +template +Foam::scalar Foam::PressureGradientForce::massAdd +( + const typename CloudType::parcelType&, + const scalar +) const +{ + return 0.0; +} + + // ************************************************************************* // diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/PressureGradient/PressureGradientForce.H b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/PressureGradient/PressureGradientForce.H index 663723a820..b1986348fb 100644 --- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/PressureGradient/PressureGradientForce.H +++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/PressureGradient/PressureGradientForce.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) 2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -38,6 +38,7 @@ SourceFiles #include "ParticleForce.H" #include "volFields.H" +#include "interpolation.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -53,13 +54,15 @@ class PressureGradientForce : public ParticleForce { - // Private data +protected: + + // Protected data //- Name of velocity field const word UName_; - //- Velocity gradient field - const volTensorField* gradUPtr_; + //- Rate of change of carrier phase velocity interpolator + autoPtr > DUcDtInterpPtr_; public: @@ -75,7 +78,8 @@ public: ( CloudType& owner, const fvMesh& mesh, - const dictionary& dict + const dictionary& dict, + const word& forceType = typeName ); //- Construct copy @@ -99,8 +103,8 @@ public: // Access - //- Return const access to the velocity gradient field - inline const volTensorField& gradU() const; + //- Return the rate of change of carrier phase velocity interpolator + inline const interpolation& DUcDtInterp() const; // Evaluation @@ -117,6 +121,13 @@ public: const scalar Re, const scalar muc ) const; + + //- Return the added mass + virtual scalar massAdd + ( + const typename CloudType::parcelType& p, + const scalar mass + ) const; }; diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/PressureGradient/PressureGradientForceI.H b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/PressureGradient/PressureGradientForceI.H index c9bd24c62c..6c085241de 100644 --- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/PressureGradient/PressureGradientForceI.H +++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/PressureGradient/PressureGradientForceI.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) 2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -26,23 +26,20 @@ License // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // template -const Foam::volTensorField& Foam::PressureGradientForce::gradU() -const +inline const Foam::interpolation& +Foam::PressureGradientForce::DUcDtInterp() const { - if (gradUPtr_) - { - return *gradUPtr_; - } - else + if (!DUcDtInterpPtr_.valid()) { FatalErrorIn ( - "const volTensorField& PressureGradientForce::gradU()" - "const" - ) << "gradU field not allocated" << abort(FatalError); - - return *reinterpret_cast(0); + "inline const Foam::interpolation&" + "Foam::PressureGradientForce::DUcDtInterp() const" + ) << "Carrier phase DUcDt interpolation object not set" + << abort(FatalError); } + + return DUcDtInterpPtr_(); } diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/VirtualMass/VirtualMassForce.C b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/VirtualMass/VirtualMassForce.C index 3dca602dd4..c701e07a00 100644 --- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/VirtualMass/VirtualMassForce.C +++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/VirtualMass/VirtualMassForce.C @@ -24,8 +24,6 @@ License \*---------------------------------------------------------------------------*/ #include "VirtualMassForce.H" -#include "fvcDdt.H" -#include "fvcGrad.H" // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // @@ -34,14 +32,12 @@ Foam::VirtualMassForce::VirtualMassForce ( CloudType& owner, const fvMesh& mesh, - const dictionary& dict + const dictionary& dict, + const word& forceType ) : - ParticleForce(owner, mesh, dict, typeName, true), - UName_(this->coeffs().template lookupOrDefault("U", "U")), - Cvm_(readScalar(this->coeffs().lookup("Cvm"))), - DUcDtPtr_(NULL), - DUcDtInterpPtr_(NULL) + PressureGradientForce(owner, mesh, dict, forceType), + Cvm_(readScalar(this->coeffs().lookup("Cvm"))) {} @@ -51,11 +47,8 @@ Foam::VirtualMassForce::VirtualMassForce const VirtualMassForce& vmf ) : - ParticleForce(vmf), - UName_(vmf.UName_), - Cvm_(vmf.Cvm_), - DUcDtPtr_(NULL), - DUcDtInterpPtr_(NULL) + PressureGradientForce(vmf), + Cvm_(vmf.Cvm_) {} @@ -71,36 +64,7 @@ Foam::VirtualMassForce::~VirtualMassForce() template void Foam::VirtualMassForce::cacheFields(const bool store) { - if (store && !DUcDtPtr_) - { - const volVectorField& Uc = this->mesh().template - lookupObject(UName_); - - DUcDtPtr_ = new volVectorField - ( - "DUcDt", - fvc::ddt(Uc) + (Uc & fvc::grad(Uc)) - ); - - DUcDtInterpPtr_.reset - ( - interpolation::New - ( - this->owner().solution().interpolationSchemes(), - *DUcDtPtr_ - ).ptr() - ); - } - else - { - DUcDtInterpPtr_.clear(); - - if (DUcDtPtr_) - { - delete DUcDtPtr_; - DUcDtPtr_ = NULL; - } - } + PressureGradientForce::cacheFields(store); } @@ -114,12 +78,10 @@ Foam::forceSuSp Foam::VirtualMassForce::calcCoupled const scalar muc ) const { - forceSuSp value(vector::zero, 0.0); + forceSuSp value = + PressureGradientForce::calcCoupled(p, dt, mass, Re, muc); - vector DUcDt = - DUcDtInterp().interpolate(p.position(), p.currentTetIndices()); - - value.Su() = mass*p.rhoc()/p.rho()*Cvm_*DUcDt; + value.Su() *= Cvm_; return value; } diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/VirtualMass/VirtualMassForce.H b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/VirtualMass/VirtualMassForce.H index c5277ec696..ea0d580932 100644 --- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/VirtualMass/VirtualMassForce.H +++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/VirtualMass/VirtualMassForce.H @@ -28,7 +28,6 @@ Description Calculates particle virtual mass force SourceFiles - VirtualMassForceI.H VirtualMassForce.C \*---------------------------------------------------------------------------*/ @@ -36,9 +35,7 @@ SourceFiles #ifndef VirtualMassForce_H #define VirtualMassForce_H -#include "ParticleForce.H" -#include "volFields.H" -#include "interpolation.H" +#include "PressureGradientForce.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -52,22 +49,13 @@ namespace Foam template class VirtualMassForce : - public ParticleForce + public PressureGradientForce { // Private data - //- Name of velocity field - const word UName_; - //- Virtual mass coefficient - typically 0.5 scalar Cvm_; - //- Rate of change of carrier phase velocity - volVectorField* DUcDtPtr_; - - //- Rate of change of carrier phase velocity interpolator - autoPtr > DUcDtInterpPtr_; - public: @@ -82,7 +70,8 @@ public: ( CloudType& owner, const fvMesh& mesh, - const dictionary& dict + const dictionary& dict, + const word& forceType = typeName ); //- Construct copy @@ -104,15 +93,6 @@ public: // Member Functions - // Access - - //- Return the rate of change of carrier phase velocity - inline const volVectorField& DUcDt() const; - - //- Return the rate of change of carrier phase velocity interpolator - inline const interpolation& DUcDtInterp() const; - - // Evaluation //- Cache fields @@ -143,8 +123,6 @@ public: // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -#include "VirtualMassForceI.H" - #ifdef NoRepository #include "VirtualMassForce.C" #endif diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/VirtualMass/VirtualMassForceI.H b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/VirtualMass/VirtualMassForceI.H deleted file mode 100644 index e156c58d14..0000000000 --- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/VirtualMass/VirtualMassForceI.H +++ /dev/null @@ -1,66 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / 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 . - -\*---------------------------------------------------------------------------*/ - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -template -const Foam::volVectorField& Foam::VirtualMassForce::DUcDt() const -{ - if (DUcDtPtr_) - { - return *DUcDtPtr_; - } - else - { - FatalErrorIn - ( - "const volVectorField& VirtualMassForce::DUcDt()" - "const" - ) << "DUcDt field not allocated" << abort(FatalError); - - return *reinterpret_cast(0); - } -} - - -template -inline const Foam::interpolation& -Foam::VirtualMassForce::DUcDtInterp() const -{ - if (!DUcDtInterpPtr_.valid()) - { - FatalErrorIn - ( - "inline const Foam::interpolation&" - "Foam::VirtualMassForce::DUcDtInterp() const" - ) << "Carrier pahase DUcDt interpolation object not set" - << abort(FatalError); - } - - return DUcDtInterpPtr_(); -} - - -// ************************************************************************* // 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