diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H b/applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H index c3da8d972a..a382b769b5 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H @@ -78,10 +78,6 @@ - fvc::div(phi2)*U2 ); - - Info<< "Calculating field g.h\n" << endl; - volScalarField gh("gh", g & mesh.C()); - volScalarField rAU1 ( IOobject diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H index b097abf745..ea3f6bdb33 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H @@ -34,6 +34,7 @@ fvc::interpolate((1.0/rho1)*rAU1*phase1.turbulence().pPrime()) *fvc::snGrad(alpha1)*mesh.magSf() ); + phiP1.boundaryField() == 0; // Phase-2 pressure flux (e.g. due to particle-particle pressure) surfaceScalarField phiP2 @@ -42,6 +43,7 @@ fvc::interpolate((1.0/rho2)*rAU2*phase2.turbulence().pPrime()) *fvc::snGrad(alpha2)*mesh.magSf() ); + phiP2.boundaryField() == 0; surfaceScalarField phiHbyA1 ( @@ -126,9 +128,11 @@ ); pEqnComp1 = - fvc::ddt(rho1) - + fvc::div(phi1, rho1) - fvc::Sp(fvc::div(phi1), rho1) - + correction + ( + fvc::ddt(alpha1, rho1) + fvc::div(alphaPhi1, rho1) + - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1) + )/rho1 + + (alpha1/rho1)*correction ( psi1*fvm::ddt(p) + fvm::div(phid1, p) - fvm::Sp(fvc::div(phid1), p) @@ -137,9 +141,11 @@ pEqnComp1().relax(); pEqnComp2 = - fvc::ddt(rho2) - + fvc::div(phi2, rho2) - fvc::Sp(fvc::div(phi2), rho2) - + correction + ( + fvc::ddt(alpha2, rho2) + fvc::div(alphaPhi2, rho2) + - fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2) + )/rho2 + + (alpha2/rho2)*correction ( psi2*fvm::ddt(p) + fvm::div(phid2, p) - fvm::Sp(fvc::div(phid2), p) @@ -150,12 +156,18 @@ else { pEqnComp1 = - fvc::ddt(rho1) + psi1*correction(fvm::ddt(p)) - + fvc::div(phi1, rho1) - fvc::Sp(fvc::div(phi1), rho1); + ( + fvc::ddt(alpha1, rho1) + fvc::div(alphaPhi1, rho1) + - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1) + )/rho1 + + (alpha1*psi1/rho1)*correction(fvm::ddt(p)); pEqnComp2 = - fvc::ddt(rho2) + psi2*correction(fvm::ddt(p)) - + fvc::div(phi2, rho2) - fvc::Sp(fvc::div(phi2), rho2); + ( + fvc::ddt(alpha2, rho2) + fvc::div(alphaPhi2, rho2) + - fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2) + )/rho2 + + (alpha2*psi2/rho2)*correction(fvm::ddt(p)); } // Cache p prior to solve for density update @@ -171,11 +183,7 @@ solve ( - ( - (alpha1/rho1)*pEqnComp1() - + (alpha2/rho2)*pEqnComp2() - ) - + pEqnIncomp, + pEqnComp1() + pEqnComp2() + pEqnIncomp, mesh.solver(p.select(pimple.finalInnerIter())) ); @@ -201,8 +209,8 @@ fluid.dgdt() = ( - pos(alpha2)*(pEqnComp2 & p)/rho2 - - pos(alpha1)*(pEqnComp1 & p)/rho1 + pos(alpha2)*(pEqnComp2 & p)/max(alpha2, scalar(1e-3)) + - pos(alpha1)*(pEqnComp1 & p)/max(alpha1, scalar(1e-3)) ); p.relax(); diff --git a/applications/utilities/postProcessing/dataConversion/foamToEnsight/cellSets.H b/applications/utilities/postProcessing/dataConversion/foamToEnsight/cellSets.H index 1f66fa276e..54887f662a 100644 --- a/applications/utilities/postProcessing/dataConversion/foamToEnsight/cellSets.H +++ b/applications/utilities/postProcessing/dataConversion/foamToEnsight/cellSets.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-2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -78,6 +78,26 @@ public: hexes(nCells), polys(nCells) {} + + + // Member Functions + + void setSize(const label nCells) + { + nTets = 0; + nPyrs = 0; + nPrisms = 0; + nHexesWedges = 0; + nPolys = 0; + + tets.setSize(nCells); + pyrs.setSize(nCells); + prisms.setSize(nCells); + wedges.setSize(nCells); + hexes.setSize(nCells); + polys.setSize(nCells); + } + }; diff --git a/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightMesh.C b/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightMesh.C index 0332572eeb..267b8dc844 100644 --- a/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightMesh.C +++ b/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightMesh.C @@ -47,7 +47,8 @@ License void Foam::ensightMesh::correct() { patchPartOffset_ = 2; - meshCellSets_ = mesh_.nCells(); + meshCellSets_.setSize(mesh_.nCells()); + boundaryFaceSets_.setSize(mesh_.boundary().size()); allPatchNames_.clear(); patchNames_.clear(); @@ -194,44 +195,41 @@ void Foam::ensightMesh::correct() { forAll(mesh_.boundary(), patchi) { - if (mesh_.boundary()[patchi].size()) + const polyPatch& p = mesh_.boundaryMesh()[patchi]; + + labelList& tris = boundaryFaceSets_[patchi].tris; + labelList& quads = boundaryFaceSets_[patchi].quads; + labelList& polys = boundaryFaceSets_[patchi].polys; + + tris.setSize(p.size()); + quads.setSize(p.size()); + polys.setSize(p.size()); + + label nTris = 0; + label nQuads = 0; + label nPolys = 0; + + forAll(p, faceI) { - const polyPatch& p = mesh_.boundaryMesh()[patchi]; + const face& f = p[faceI]; - labelList& tris = boundaryFaceSets_[patchi].tris; - labelList& quads = boundaryFaceSets_[patchi].quads; - labelList& polys = boundaryFaceSets_[patchi].polys; - - tris.setSize(p.size()); - quads.setSize(p.size()); - polys.setSize(p.size()); - - label nTris = 0; - label nQuads = 0; - label nPolys = 0; - - forAll(p, faceI) + if (f.size() == 3) { - const face& f = p[faceI]; - - if (f.size() == 3) - { - tris[nTris++] = faceI; - } - else if (f.size() == 4) - { - quads[nQuads++] = faceI; - } - else - { - polys[nPolys++] = faceI; - } + tris[nTris++] = faceI; + } + else if (f.size() == 4) + { + quads[nQuads++] = faceI; + } + else + { + polys[nPolys++] = faceI; } - - tris.setSize(nTris); - quads.setSize(nQuads); - polys.setSize(nPolys); } + + tris.setSize(nTris); + quads.setSize(nQuads); + polys.setSize(nPolys); } } @@ -242,12 +240,9 @@ void Foam::ensightMesh::correct() if (patchNames_.empty() || patchNames_.found(patchName)) { - if (mesh_.boundary()[patchi].size()) - { - nfp.nTris = boundaryFaceSets_[patchi].tris.size(); - nfp.nQuads = boundaryFaceSets_[patchi].quads.size(); - nfp.nPolys = boundaryFaceSets_[patchi].polys.size(); - } + nfp.nTris = boundaryFaceSets_[patchi].tris.size(); + nfp.nQuads = boundaryFaceSets_[patchi].quads.size(); + nfp.nPolys = boundaryFaceSets_[patchi].polys.size(); } reduce(nfp.nTris, sumOp