diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/UEqn.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/UEqn.H index 2e8f979be4..61d805e01d 100644 --- a/applications/solvers/lagrangian/reactingParcelFilmFoam/UEqn.H +++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/UEqn.H @@ -4,13 +4,23 @@ + fvm::div(phi, U) + turbulence->divDevRhoReff(U) == - rho.dimensionedInternalField()*g - + parcels.SU(U) + parcels.SU(U) ); UEqn.relax(); if (momentumPredictor) { - solve(UEqn == -fvc::grad(p)); + solve + ( + UEqn + == + fvc::reconstruct + ( + ( + - ghf*fvc::snGrad(rho) + - fvc::snGrad(p_rgh) + )*mesh.magSf() + ) + ); } diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/YEqn.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/YEqn.H index 8526b590a2..db79376614 100644 --- a/applications/solvers/lagrangian/reactingParcelFilmFoam/YEqn.H +++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/YEqn.H @@ -5,7 +5,7 @@ tmp > mvConvection mesh, fields, phi, - mesh.divScheme("div(phi,Yi_h)") + mesh.divScheme("div(phi,Yi_hs)") ) ); @@ -23,7 +23,7 @@ tmp > mvConvection ( fvm::ddt(rho, Yi) + mvConvection->fvmDiv(phi, Yi) - - fvm::laplacian(turbulence->muEff(), Yi) + - fvm::laplacian(turbulence->alphaEff(), Yi) == parcels.SYi(i, Yi) + surfaceFilm.Srho(i) diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/chemistry.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/chemistry.H index 5dd7ce9cb0..ce36fa0d1b 100644 --- a/applications/solvers/lagrangian/reactingParcelFilmFoam/chemistry.H +++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/chemistry.H @@ -1,3 +1,4 @@ +if (chemistry.chemistry()) { Info << "Solving chemistry" << endl; diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/createFields.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/createFields.H index 2e3208b0ee..a5a87c73ef 100644 --- a/applications/solvers/lagrangian/reactingParcelFilmFoam/createFields.H +++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/createFields.H @@ -15,11 +15,6 @@ const word inertSpecie(thermo.lookup("inertSpecie")); - volScalarField& p = thermo.p(); - volScalarField& hs = thermo.hs(); - const volScalarField& T = thermo.T(); - const volScalarField& psi = thermo.psi(); - Info<< "Creating field rho\n" << endl; volScalarField rho ( @@ -34,6 +29,11 @@ thermo.rho() ); + volScalarField& p = thermo.p(); + volScalarField& hs = thermo.hs(); + const volScalarField& T = thermo.T(); + const volScalarField& psi = thermo.psi(); + Info<< "\nReading field U\n" << endl; volVectorField U ( @@ -84,6 +84,28 @@ fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p) ); + + Info<< "Calculating field g.h\n" << endl; + volScalarField gh("gh", g & mesh.C()); + + surfaceScalarField ghf("gh", g & mesh.Cf()); + + volScalarField p_rgh + ( + IOobject + ( + "p_rgh", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + // Force p_rgh to be consistent with p + p_rgh = p - rho*gh; + multivariateSurfaceInterpolationScheme::fieldTable fields; forAll(Y, i) diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/hsEqn.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/hsEqn.H index 0cb2318252..feb112f652 100644 --- a/applications/solvers/lagrangian/reactingParcelFilmFoam/hsEqn.H +++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/hsEqn.H @@ -19,4 +19,6 @@ thermo.correct(); radiation->correct(); + + Info<< "min/max(T) = " << min(T).value() << ", " << max(T).value() << endl; } diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H index 9d028d4502..ff7c716968 100644 --- a/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H +++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H @@ -1,74 +1,58 @@ rho = thermo.rho(); volScalarField rAU(1.0/UEqn.A()); +surfaceScalarField rhorAUf(rAU.name(), fvc::interpolate(rho*rAU)); U = rAU*UEqn.H(); -if (transonic) +surfaceScalarField phiU +( + fvc::interpolate(rho) + *( + (fvc::interpolate(U) & mesh.Sf()) + + fvc::ddtPhiCorr(rAU, rho, U, phi) + ) +); + +phi = phiU - rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf(); + +for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) { - surfaceScalarField phid + fvScalarMatrix p_rghEqn ( - "phid", - fvc::interpolate(psi) - *( - (fvc::interpolate(U) & mesh.Sf()) - + fvc::ddtPhiCorr(rAU, rho, U, phi) + fvc::ddt(psi, rho)*gh + + fvc::div(phi) + + fvm::ddt(psi, p_rgh) + - fvm::laplacian(rhorAUf, p_rgh) + == + parcels.Srho() + + surfaceFilm.Srho() + ); + + p_rghEqn.solve + ( + mesh.solver + ( + p_rgh.select + ( + pimpleCorr.finalIter() + && corr == nCorr-1 + && nonOrth == nNonOrthCorr + ) ) ); - for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) + if (nonOrth == nNonOrthCorr) { - fvScalarMatrix pEqn - ( - fvm::ddt(psi, p) - + fvm::div(phid, p) - - fvm::laplacian(rho*rAU, p) - == - parcels.Srho() - + surfaceFilm.Srho() - ); - - pEqn.solve(); - - if (nonOrth == nNonOrthCorr) - { - phi == pEqn.flux(); - } + phi += p_rghEqn.flux(); } } -else -{ - phi = - fvc::interpolate(rho) - *( - (fvc::interpolate(U) & mesh.Sf()) - + fvc::ddtPhiCorr(rAU, rho, U, phi) - ); - for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) - { - fvScalarMatrix pEqn - ( - fvm::ddt(psi, p) - + fvc::div(phi) - - fvm::laplacian(rho*rAU, p) - == - parcels.Srho() - + surfaceFilm.Srho() - ); - - pEqn.solve(); - - if (nonOrth == nNonOrthCorr) - { - phi += pEqn.flux(); - } - } -} +p = p_rgh + rho*gh; #include "rhoEqn.H" #include "compressibleContinuityErrs.H" -U -= rAU*fvc::grad(p); +U += rAU*fvc::reconstruct((phi - phiU)/rhorAUf); U.correctBoundaryConditions(); DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p); diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C b/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C index f788da57e4..9ec100f3fd 100644 --- a/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C +++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2009-2010 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2009-2011 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -39,6 +39,7 @@ Description #include "chemistrySolver.H" #include "radiationModel.H" #include "SLGThermo.H" +#include "pimpleLoop.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -66,8 +67,9 @@ int main(int argc, char *argv[]) while (runTime.run()) { #include "readTimeControls.H" - #include "readPISOControls.H" + #include "readPIMPLEControls.H" #include "compressibleCourantNo.H" + #include "setMultiRegionDeltaT.H" #include "setDeltaT.H" runTime++; @@ -84,20 +86,22 @@ int main(int argc, char *argv[]) #include "rhoEqn.H" // --- PIMPLE loop - for (int ocorr=1; ocorr<=nOuterCorr; ocorr++) + for + ( + pimpleLoop pimpleCorr(mesh, nOuterCorr); + pimpleCorr.loop(); + pimpleCorr++ + ) { #include "UEqn.H" #include "YEqn.H" + #include "hsEqn.H" // --- PISO loop for (int corr=1; corr<=nCorr; corr++) { - #include "hsEqn.H" #include "pEqn.H" } - - Info<< "T gas min/max = " << min(T).value() << ", " - << max(T).value() << endl; } turbulence->correct(); diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/setMultiRegionDeltaT.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/setMultiRegionDeltaT.H new file mode 100644 index 0000000000..037a09f93d --- /dev/null +++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/setMultiRegionDeltaT.H @@ -0,0 +1,57 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2010-2011 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 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 . + +Global + setMultiRegionDeltaT + +Description + Reset the timestep to maintain a constant maximum Courant numbers. + Reduction of time-step is immediate, but increase is damped to avoid + unstable oscillations. + +\*---------------------------------------------------------------------------*/ + +if (adjustTimeStep) +{ + if (CoNum == -GREAT) + { + CoNum = SMALL; + } + + const scalar TFactorFluid = maxCo/(CoNum + SMALL); + const scalar TFactorFilm = maxCo/(surfaceFilm.CourantNumber() + SMALL); + + const scalar dt0 = runTime.deltaTValue(); + + runTime.setDeltaT + ( + min + ( + dt0*min(min(TFactorFluid, TFactorFilm), 1.2), + maxDeltaT + ) + ); +} + + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/bubbleFoam/wallFunctions.H b/applications/solvers/multiphase/bubbleFoam/wallFunctions.H index a998d4c3dd..bc9e93c63d 100644 --- a/applications/solvers/multiphase/bubbleFoam/wallFunctions.H +++ b/applications/solvers/multiphase/bubbleFoam/wallFunctions.H @@ -34,7 +34,7 @@ { const scalarField& nuw = nutb.boundaryField()[patchi]; - scalarField magFaceGradU(mag(U.boundaryField()[patchi].snGrad())); + scalarField magFaceGradU(mag(Ub.boundaryField()[patchi].snGrad())); forAll(currPatch, facei) { diff --git a/applications/solvers/multiphase/interMixingFoam/incompressibleThreePhaseMixture/threePhaseMixture.C b/applications/solvers/multiphase/interMixingFoam/incompressibleThreePhaseMixture/threePhaseMixture.C index 3cb03fc6d8..5edf003719 100644 --- a/applications/solvers/multiphase/interMixingFoam/incompressibleThreePhaseMixture/threePhaseMixture.C +++ b/applications/solvers/multiphase/interMixingFoam/incompressibleThreePhaseMixture/threePhaseMixture.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -36,6 +36,10 @@ Class //- Calculate and return the laminar viscosity void Foam::threePhaseMixture::calcNu() { + nuModel1_->correct(); + nuModel2_->correct(); + nuModel3_->correct(); + // Average kinematic viscosity calculated from dynamic viscosity nu_ = mu()/(alpha1_*rho1_ + alpha2_*rho2_ + alpha3_*rho3_); } diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C b/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C index c8561a6d67..90e6a379e0 100644 --- a/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C +++ b/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C @@ -179,7 +179,7 @@ Foam::label Foam::checkTopology if (nTwoCells > 0) { Info<< " < tokenCompound##Name##_; \ + typedef token::Compound tokenCompound##Name##_; \ defineTemplateTypeNameAndDebugWithName(tokenCompound##Name##_, #Type, 0); #define addCompoundToRunTimeSelectionTable(Type, Name) \ - token::compound::addIstreamConstructorToTable > \ + token::compound::addIstreamConstructorToTable > \ add##Name##IstreamConstructorToTable_; diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C index 1c59cae4a6..80ec5ae74c 100644 --- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C +++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C @@ -100,13 +100,7 @@ template template void Foam::KinematicCloud::solve(TrackData& td) { - if (solution_.transient()) - { - td.cloud().preEvolve(); - - evolveCloud(td); - } - else + if (solution_.steadyState()) { td.cloud().storeState(); @@ -114,7 +108,21 @@ void Foam::KinematicCloud::solve(TrackData& td) evolveCloud(td); - td.cloud().relaxSources(td.cloud().cloudCopy()); + if (solution_.coupled()) + { + td.cloud().relaxSources(td.cloud().cloudCopy()); + } + } + else + { + td.cloud().preEvolve(); + + evolveCloud(td); + + if (solution_.coupled()) + { + td.cloud().scaleSources(); + } } td.cloud().info(); @@ -258,6 +266,7 @@ void Foam::KinematicCloud::cloudReset(KinematicCloud& c) injectionModel_.reset(c.injectionModel_.ptr()); patchInteractionModel_.reset(c.patchInteractionModel_.ptr()); postProcessingModel_.reset(c.postProcessingModel_.ptr()); + surfaceFilmModel_.reset(c.surfaceFilmModel_.ptr()); UIntegrator_.reset(c.UIntegrator_.ptr()); } @@ -556,11 +565,23 @@ void Foam::KinematicCloud::relax ) const { const scalar coeff = solution_.relaxCoeff(name); - field = field0 + coeff*(field - field0); } +template +template +void Foam::KinematicCloud::scale +( + DimensionedField& field, + const word& name +) const +{ + const scalar coeff = solution_.relaxCoeff(name); + field *= coeff; +} + + template void Foam::KinematicCloud::relaxSources ( @@ -568,6 +589,15 @@ void Foam::KinematicCloud::relaxSources ) { this->relax(UTrans_(), cloudOldTime.UTrans(), "U"); + this->relax(UCoeff_(), cloudOldTime.UCoeff(), "U"); +} + + +template +void Foam::KinematicCloud::scaleSources() +{ + this->scale(UTrans_(), "U"); + this->scale(UCoeff_(), "U"); } diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H index 765abfedd7..88eff6449f 100644 --- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H +++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H @@ -498,9 +498,20 @@ public: const word& name ) const; + //- Scale field + template + void scale + ( + DimensionedField& field, + const word& name + ) const; + //- Apply relaxation to (steady state) cloud sources void relaxSources(const KinematicCloud& cloudOldTime); + //- Apply scaling to (transient) cloud sources + void scaleSources(); + //- Evolve the cloud void evolve(); diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H index 5987192bfd..568d089030 100644 --- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H +++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H @@ -370,42 +370,22 @@ Foam::KinematicCloud::theta() const false ), mesh_, - dimensionedScalar("zero", dimless, 0.0) + dimensionedScalar("zero", dimless, 0.0), + zeroGradientFvPatchScalarField::typeName ) ); volScalarField& theta = ttheta(); - theta.boundaryField() == 0; - forAllConstIter(typename KinematicCloud, *this, iter) { const parcelType& p = iter(); const label cellI = p.cell(); - if ((p.face() != -1)) - { - const label patchI = p.patch(p.face()); - if (patchI != -1) - { - scalarField& thetap = theta.boundaryField()[patchI]; - const label faceI = p.patchFace(patchI, p.face()); - thetap[faceI] += p.nParticle()*p.areaP(); - } - } - theta[cellI] += p.nParticle()*p.volume(); } theta.internalField() /= mesh_.V(); - - forAll(theta.boundaryField(), patchI) - { - scalarField& thetap = theta.boundaryField()[patchI]; - if (thetap.size() > 0) - { - thetap /= mesh_.magSf().boundaryField()[patchI]; - } - } + theta.correctBoundaryConditions(); return ttheta; } diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.C b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.C index 165a7be0b6..75b987c24b 100644 --- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.C +++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.C @@ -116,40 +116,47 @@ void Foam::cloudSolution::read() dict_.lookup("calcFrequency") >> calcFrequency_; dict_.lookup("maxCo") >> maxCo_; dict_.lookup("maxTrackTime") >> maxTrackTime_; - dict_.subDict("sourceTerms").lookup("resetOnStartup") - >> resetSourcesOnStartup_; + + if (coupled_) + { + dict_.subDict("sourceTerms").lookup("resetOnStartup") + >> resetSourcesOnStartup_; + } } - const dictionary& - schemesDict(dict_.subDict("sourceTerms").subDict("schemes")); - - wordList vars(schemesDict.toc()); - schemes_.setSize(vars.size()); - forAll(vars, i) + if (coupled_) { - // read solution variable name - schemes_[i].first() = vars[i]; + const dictionary& + schemesDict(dict_.subDict("sourceTerms").subDict("schemes")); - // set semi-implicit (1) explicit (0) flag - Istream& is = schemesDict.lookup(vars[i]); - const word scheme(is); - if (scheme == "semiImplicit") + wordList vars(schemesDict.toc()); + schemes_.setSize(vars.size()); + forAll(vars, i) { - schemes_[i].second().first() = true; - } - else if (scheme == "explicit") - { - schemes_[i].second().first() = false; - } - else - { - FatalErrorIn("void cloudSolution::read()") - << "Invalid scheme " << scheme << ". Valid schemes are " - << "explicit and semiImplicit" << exit(FatalError); - } + // read solution variable name + schemes_[i].first() = vars[i]; - // read under-relaxation factor - is >> schemes_[i].second().second(); + // set semi-implicit (1) explicit (0) flag + Istream& is = schemesDict.lookup(vars[i]); + const word scheme(is); + if (scheme == "semiImplicit") + { + schemes_[i].second().first() = true; + } + else if (scheme == "explicit") + { + schemes_[i].second().first() = false; + } + else + { + FatalErrorIn("void cloudSolution::read()") + << "Invalid scheme " << scheme << ". Valid schemes are " + << "explicit and semiImplicit" << exit(FatalError); + } + + // read under-relaxation factor + is >> schemes_[i].second().second(); + } } } diff --git a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C index eb1f66e85a..c536cd84d6 100644 --- a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C +++ b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C @@ -295,10 +295,10 @@ void Foam::ReactingCloud::relaxSources const ReactingCloud& cloudOldTime ) { - typedef DimensionedField dsfType; - CloudType::relaxSources(cloudOldTime); + typedef DimensionedField dsfType; + forAll(rhoTrans_, fieldI) { dsfType& rhoT = rhoTrans_[fieldI]; @@ -308,6 +308,21 @@ void Foam::ReactingCloud::relaxSources } +template +void Foam::ReactingCloud::scaleSources() +{ + CloudType::scaleSources(); + + typedef DimensionedField dsfType; + + forAll(rhoTrans_, fieldI) + { + dsfType& rhoT = rhoTrans_[fieldI]; + this->scale(rhoT, "rho"); + } +} + + template void Foam::ReactingCloud::evolve() { diff --git a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.H b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.H index 6d46d08bb7..324636cca5 100644 --- a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.H +++ b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.H @@ -287,6 +287,9 @@ public: //- Apply relaxation to (steady state) cloud sources void relaxSources(const ReactingCloud& cloudOldTime); + //- Apply scaling to (transient) cloud sources + void scaleSources(); + //- Evolve the cloud void evolve(); diff --git a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C index e0365b54e2..91cb53aff0 100644 --- a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C +++ b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C @@ -290,6 +290,17 @@ void Foam::ThermoCloud::relaxSources CloudType::relaxSources(cloudOldTime); this->relax(hsTrans_(), cloudOldTime.hsTrans(), "hs"); + this->relax(hsCoeff_(), cloudOldTime.hsCoeff(), "hs"); +} + + +template +void Foam::ThermoCloud::scaleSources() +{ + CloudType::scaleSources(); + + this->scale(hsTrans_(), "hs"); + this->scale(hsCoeff_(), "hs"); } diff --git a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.H b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.H index 44656d07fc..f725dca73e 100644 --- a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.H +++ b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.H @@ -301,6 +301,9 @@ public: //- Apply relaxation to (steady state) cloud sources void relaxSources(const ThermoCloud& cloudOldTime); + //- Apply scaling to (transient) cloud sources + void scaleSources(); + //- Evolve the cloud void evolve(); diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C index 6ee0d77474..0d605da15d 100644 --- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C +++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C @@ -327,7 +327,6 @@ void Foam::ReactingMultiphaseParcel::calc - // Heat transfer // ~~~~~~~~~~~~~ @@ -369,7 +368,7 @@ void Foam::ReactingMultiphaseParcel::calc d0, U0, rho0, - mass0, + 0.5*(mass0 + mass1), Su, dUTrans, Spu @@ -386,15 +385,11 @@ void Foam::ReactingMultiphaseParcel::calc { label gid = composition.localToGlobalCarrierId(GAS, i); td.cloud().rhoTrans(gid)[cellI] += np0*dMassGas[i]; -// td.cloud().hsTrans()[cellI] += -// np0*dMassGas[i]*composition.carrier().Hs(gid, T0); } forAll(YLiquid_, i) { label gid = composition.localToGlobalCarrierId(LIQ, i); td.cloud().rhoTrans(gid)[cellI] += np0*dMassLiquid[i]; -// td.cloud().hsTrans()[cellI] += -// np0*dMassLiquid[i]*composition.carrier().Hs(gid, T0); } /* // No mapping between solid components and carrier phase @@ -402,15 +397,11 @@ void Foam::ReactingMultiphaseParcel::calc { label gid = composition.localToGlobalCarrierId(SLD, i); td.cloud().rhoTrans(gid)[cellI] += np0*dMassSolid[i]; -// td.cloud().hsTrans()[cellI] += -// np0*dMassSolid[i]*composition.carrier().Hs(gid, T0); } */ forAll(dMassSRCarrier, i) { td.cloud().rhoTrans(i)[cellI] += np0*dMassSRCarrier[i]; -// td.cloud().hsTrans()[cellI] += -// np0*dMassSRCarrier[i]*composition.carrier().Hs(i, T0); } // Update momentum transfer @@ -458,8 +449,8 @@ void Foam::ReactingMultiphaseParcel::calc } */ td.cloud().UTrans()[cellI] += np0*mass1*U1; - td.cloud().hsTrans()[cellI] += - np0*mass1*HEff(td, pc, T1, idG, idL, idS); // using total h + + // enthalpy transfer accounted for via change in mass fractions } } diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C index 56f24c97f1..b5a2706732 100644 --- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C +++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C @@ -379,13 +379,12 @@ void Foam::ReactingParcel::calc d0, U0, rho0, - mass0, + 0.5*(mass0 + mass1), Su, dUTrans, Spu ); - dUTrans += 0.5*(mass0 - mass1)*(U0 + U1); // Accumulate carrier phase source terms // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -396,8 +395,6 @@ void Foam::ReactingParcel::calc { label gid = composition.localToGlobalCarrierId(0, i); td.cloud().rhoTrans(gid)[cellI] += np0*dMassPC[i]; -// td.cloud().hsTrans()[cellI] += -// np0*dMassPC[i]*composition.carrier().Hs(gid, T0); } // Update momentum transfer @@ -429,8 +426,8 @@ void Foam::ReactingParcel::calc td.cloud().rhoTrans(gid)[cellI] += np0*mass1*Y_[i]; } td.cloud().UTrans()[cellI] += np0*mass1*U1; - td.cloud().hsTrans()[cellI] += - np0*mass1*composition.H(0, Y_, pc_, T1); + + // enthalpy transfer accounted for via change in mass fractions } } diff --git a/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.C b/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.C index 07ec7cc21e..cb0bc500a3 100644 --- a/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.C +++ b/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.C @@ -148,7 +148,8 @@ void Foam::SurfaceFilmModel::inject(TrackData& td) const labelList& filmPatches = filmModel.intCoupledPatchIDs(); const labelList& primaryPatches = filmModel.primaryPatchIDs(); - const polyBoundaryMesh& pbm = this->owner().mesh().boundaryMesh(); + const fvMesh& mesh = this->owner().mesh(); + const polyBoundaryMesh& pbm = mesh.boundaryMesh(); forAll(filmPatches, i) { @@ -163,6 +164,10 @@ void Foam::SurfaceFilmModel::inject(TrackData& td) cacheFilmFields(filmPatchI, primaryPatchI, distMap, filmModel); + const vectorField& Cf = mesh.C().boundaryField()[primaryPatchI]; + const vectorField& Sf = mesh.Sf().boundaryField()[primaryPatchI]; + const scalarField& magSf = mesh.magSf().boundaryField()[primaryPatchI]; + forAll(injectorCellsPatch, j) { if (diameterParcelPatch_[j] > 0) @@ -179,7 +184,15 @@ void Foam::SurfaceFilmModel::inject(TrackData& td) const label tetFaceI = this->owner().mesh().cells()[cellI][0]; const label tetPtI = 1; - const point& pos = this->owner().mesh().C()[cellI]; +// const point& pos = this->owner().mesh().C()[cellI]; + + const scalar offset = + max + ( + diameterParcelPatch_[j], + deltaFilmPatch_[primaryPatchI][j] + ); + const point pos = Cf[j] - 1.1*offset*Sf[j]/magSf[j]; // Create a new parcel parcelType* pPtr = @@ -217,11 +230,11 @@ void Foam::SurfaceFilmModel::cacheFilmFields const regionModels::surfaceFilmModels::surfaceFilmModel& filmModel ) { - massParcelPatch_ = filmModel.massForPrimary().boundaryField()[filmPatchI]; + massParcelPatch_ = filmModel.cloudMassTrans().boundaryField()[filmPatchI]; distMap.distribute(massParcelPatch_); diameterParcelPatch_ = - filmModel.diametersForPrimary().boundaryField()[filmPatchI]; + filmModel.cloudDiameterTrans().boundaryField()[filmPatchI]; distMap.distribute(diameterParcelPatch_); UFilmPatch_ = filmModel.Us().boundaryField()[filmPatchI]; diff --git a/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C b/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C index 7249846b13..541e983155 100644 --- a/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C +++ b/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C @@ -391,8 +391,8 @@ void Foam::ThermoSurfaceFilm::splashInteraction const scalar dBarSplash = 1/cbrt(6.0)*cbrt(mRatio/Ns)*d + ROOTVSMALL; // cumulative diameter splash distribution - const scalar dMax = cbrt(mRatio)*d; - const scalar dMin = 0.001*dMax; + const scalar dMax = 0.9*cbrt(mRatio)*d; + const scalar dMin = 0.1*dMax; const scalar K = exp(-dMin/dBarSplash) - exp(-dMax/dBarSplash); // surface energy of secondary parcels [J] @@ -437,7 +437,7 @@ void Foam::ThermoSurfaceFilm::splashInteraction // magnitude of the normal velocity of the first splashed parcel const scalar magUns0 = - sqrt(2.0*parcelsPerSplash_*EKs/mSplash/(1 + coeff1/sqr(coeff2))); + sqrt(2.0*parcelsPerSplash_*EKs/mSplash/(1.0 + coeff1/sqr(coeff2))); // Set splashed parcel properties forAll(dNew, i) @@ -467,7 +467,7 @@ void Foam::ThermoSurfaceFilm::splashInteraction // Apply correction to velocity for 2-D cases meshTools::constrainDirection(mesh, mesh.solutionD(), pPtr->U()); - +Info<< "NEW PARTICLE: " << *pPtr << endl; // Add the new parcel this->owner().addParticle(pPtr); diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C index 3787156a53..6c1925f064 100644 --- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C +++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C @@ -722,10 +722,15 @@ void Foam::autoSnapDriver::preSmoothPatch if (debug) { const_cast(mesh.time())++; - Pout<< "Writing patch smoothed mesh to time " << meshRefiner_.timeName() - << endl; - - mesh.write(); + Pout<< "Writing patch smoothed mesh to time " + << meshRefiner_.timeName() << '.' << endl; + meshRefiner_.write + ( + debug, + mesh.time().path()/meshRefiner_.timeName() + ); + Pout<< "Dumped mesh in = " + << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl; } Info<< "Patch points smoothed in = " @@ -956,25 +961,6 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface vector(GREAT, GREAT, GREAT) // null value (note: cannot use VGREAT) ); - - // Check for displacement being outwards. - outwardsDisplacement(pp, patchDisp); - - // Set initial distribution of displacement field (on patches) from - // patchDisp and make displacement consistent with b.c. on displacement - // pointVectorField. - meshMover.setDisplacement(patchDisp); - - if (debug) - { - dumpMove - ( - mesh.time().path()/"patchDisplacement.obj", - pp.localPoints(), - pp.localPoints() + patchDisp - ); - } - return patchDisp; } @@ -1128,8 +1114,8 @@ Foam::autoPtr Foam::autoSnapDriver::repatchToSurface indirectPrimitivePatch& pp = ppPtr(); // Divide surfaces into zoned and unzoned - labelList zonedSurfaces = meshRefiner_.surfaces().getNamedSurfaces(); - labelList unzonedSurfaces = meshRefiner_.surfaces().getUnnamedSurfaces(); + labelList zonedSurfaces = surfaces.getNamedSurfaces(); + labelList unzonedSurfaces = surfaces.getUnnamedSurfaces(); // Faces that do not move @@ -1344,19 +1330,6 @@ void Foam::autoSnapDriver::doSnap // Pre-smooth patch vertices (so before determining nearest) preSmoothPatch(snapParams, nInitErrors, baffles, meshMover); - if (debug) - { - Pout<< "Writing patch smoothed mesh to time " - << meshRefiner_.timeName() << '.' << endl; - meshRefiner_.write - ( - debug, - mesh.time().path()/meshRefiner_.timeName() - ); - Pout<< "Dumped mesh in = " - << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl; - } - for (label iter = 0; iter < nFeatIter; iter++) { @@ -1382,6 +1355,15 @@ void Foam::autoSnapDriver::doSnap ); } + // Check for displacement being outwards. + outwardsDisplacement(ppPtr(), disp); + + // Set initial distribution of displacement field (on patches) + // from patchDisp and make displacement consistent with b.c. + // on displacement pointVectorField. + meshMover.setDisplacement(disp); + + if (debug&meshRefinement::OBJINTERSECTIONS) { dumpMove @@ -1467,7 +1449,7 @@ void Foam::autoSnapDriver::doSnap meshRefiner_.write ( debug, - mesh.time().path()/meshRefiner_.timeName() + meshRefiner_.timeName() ); } } diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H index 4b0e1a30fd..59cd772a98 100644 --- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H +++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H @@ -117,12 +117,12 @@ class autoSnapDriver const List& constraints, vectorField& disp ) const; - void calcNearest - ( - const pointField& points, - vectorField& disp, - vectorField& surfaceNormal - ) const; + //void calcNearest + //( + // const pointField& points, + // vectorField& disp, + // vectorField& surfaceNormal + //) const; void calcNearestFace ( const label iter, diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C index cc73efbe51..0227234822 100644 --- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C +++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C @@ -135,63 +135,6 @@ void Foam::autoSnapDriver::smoothAndConstrain } -void Foam::autoSnapDriver::calcNearest -( - const pointField& points, - vectorField& disp, - vectorField& surfaceNormal -) const -{ - const refinementSurfaces& surfaces = meshRefiner_.surfaces(); - - // Displacement and orientation per pp face. - disp.setSize(points.size()); - disp = vector::zero; - surfaceNormal.setSize(points.size()); - surfaceNormal = vector::zero; - - { - // Divide surfaces into zoned and unzoned - labelList zonedSurfaces = - meshRefiner_.surfaces().getNamedSurfaces(); - labelList unzonedSurfaces = - meshRefiner_.surfaces().getUnnamedSurfaces(); - - { - List hitInfo; - labelList hitSurface; - labelList hitRegion; - surfaces.findNearestRegion - ( - unzonedSurfaces, - points, - sqr(scalarField(points.size(), GREAT)),// sqr of attract dist - hitSurface, - hitInfo, - hitRegion, - surfaceNormal - ); - - forAll(hitInfo, i) - { - if (hitInfo[i].hit()) - { - disp[i] = - hitInfo[i].hitPoint() - - points[i]; - } - else - { - WarningIn("calcNearest(..)") - << "Did not hit anything from face:" << i - << " at:" << points[i] << endl; - } - } - } - } -} - - void Foam::autoSnapDriver::calcNearestFace ( const label iter, @@ -202,6 +145,7 @@ void Foam::autoSnapDriver::calcNearestFace ) const { const fvMesh& mesh = meshRefiner_.mesh(); + const refinementSurfaces& surfaces = meshRefiner_.surfaces(); // Displacement and orientation per pp face. faceDisp.setSize(pp.size()); @@ -209,7 +153,148 @@ void Foam::autoSnapDriver::calcNearestFace faceSurfaceNormal.setSize(pp.size()); faceSurfaceNormal = vector::zero; - calcNearest(pp.faceCentres(), faceDisp, faceSurfaceNormal); + // Divide surfaces into zoned and unzoned + labelList zonedSurfaces = surfaces.getNamedSurfaces(); + labelList unzonedSurfaces = surfaces.getUnnamedSurfaces(); + + // Per pp face the current surface snapped to + labelList snapSurf(pp.size(), -1); + + + // Do zoned surfaces + // ~~~~~~~~~~~~~~~~~ + // Zoned faces only attract to corresponding surface + + // Extract faces per zone + const wordList& faceZoneNames = surfaces.faceZoneNames(); + + forAll(zonedSurfaces, i) + { + label zoneSurfI = zonedSurfaces[i]; + + // Get indices of faces on pp that are also in zone + label zoneI = mesh.faceZones().findZoneID(faceZoneNames[zoneSurfI]); + if (zoneI == -1) + { + FatalErrorIn + ( + "autoSnapDriver::calcNearestFace(..)" + ) << "Problem. Cannot find zone " << faceZoneNames[zoneSurfI] + << exit(FatalError); + } + const faceZone& fZone = mesh.faceZones()[zoneI]; + PackedBoolList isZonedFace(mesh.nFaces()); + forAll(fZone, i) + { + isZonedFace[fZone[i]] = 1; + } + + DynamicList