diff --git a/src/OpenFOAM/db/dictionary/dictionary.C b/src/OpenFOAM/db/dictionary/dictionary.C index 5e67648d17aca1e61ded720d947772e9bdbef972..40abb17d470ed5358919c4f46c9139a4577caed9 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/finiteVolume/fields/fvPatchFields/derived/uniformTotalPressure/uniformTotalPressureFvPatchScalarField.C b/src/finiteVolume/fields/fvPatchFields/derived/uniformTotalPressure/uniformTotalPressureFvPatchScalarField.C index e9bd3270727c8050c5fad99ec892cbe0ea3e236b..c3b35384cdaa0cd3ab1a8e06d972eea42131bb4a 100644 --- a/src/finiteVolume/fields/fvPatchFields/derived/uniformTotalPressure/uniformTotalPressureFvPatchScalarField.C +++ b/src/finiteVolume/fields/fvPatchFields/derived/uniformTotalPressure/uniformTotalPressureFvPatchScalarField.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 @@ -44,7 +44,6 @@ uniformTotalPressureFvPatchScalarField rhoName_("none"), psiName_("none"), gamma_(0.0), - p0_(0.0), pressure_() {} @@ -63,7 +62,6 @@ uniformTotalPressureFvPatchScalarField rhoName_(dict.lookupOrDefault<word>("rho", "none")), psiName_(dict.lookupOrDefault<word>("psi", "none")), gamma_(readScalar(dict.lookup("gamma"))), - p0_(readScalar(dict.lookup("p0"))), pressure_(DataEntry<scalar>::New("pressure", dict)) { if (dict.found("value")) @@ -75,7 +73,8 @@ uniformTotalPressureFvPatchScalarField } else { - fvPatchField<scalar>::operator=(p0_); + scalar p0 = pressure_->value(this->db().time().timeOutputValue()); + fvPatchField<scalar>::operator=(p0); } } @@ -95,7 +94,6 @@ uniformTotalPressureFvPatchScalarField rhoName_(ptf.rhoName_), psiName_(ptf.psiName_), gamma_(ptf.gamma_), - p0_(ptf.p0_), pressure_(ptf.pressure_().clone().ptr()) {} @@ -112,7 +110,6 @@ uniformTotalPressureFvPatchScalarField rhoName_(tppsf.rhoName_), psiName_(tppsf.psiName_), gamma_(tppsf.gamma_), - p0_(tppsf.p0_), pressure_(tppsf.pressure_().clone().ptr()) {} @@ -130,7 +127,6 @@ uniformTotalPressureFvPatchScalarField rhoName_(tppsf.rhoName_), psiName_(tppsf.psiName_), gamma_(tppsf.gamma_), - p0_(tppsf.p0_), pressure_(tppsf.pressure_().clone().ptr()) {} @@ -147,14 +143,14 @@ void Foam::uniformTotalPressureFvPatchScalarField::updateCoeffs return; } - p0_ = pressure_->value(this->db().time().timeOutputValue()); + scalar p0 = pressure_->value(this->db().time().timeOutputValue()); const fvsPatchField<scalar>& phip = patch().lookupPatchField<surfaceScalarField, scalar>(phiName_); if (psiName_ == "none" && rhoName_ == "none") { - operator==(p0_ - 0.5*(1.0 - pos(phip))*magSqr(Up)); + operator==(p0 - 0.5*(1.0 - pos(phip))*magSqr(Up)); } else if (rhoName_ == "none") { @@ -167,7 +163,7 @@ void Foam::uniformTotalPressureFvPatchScalarField::updateCoeffs operator== ( - p0_ + p0 /pow ( (1.0 + 0.5*psip*gM1ByG*(1.0 - pos(phip))*magSqr(Up)), @@ -177,7 +173,7 @@ void Foam::uniformTotalPressureFvPatchScalarField::updateCoeffs } else { - operator==(p0_/(1.0 + 0.5*psip*(1.0 - pos(phip))*magSqr(Up))); + operator==(p0/(1.0 + 0.5*psip*(1.0 - pos(phip))*magSqr(Up))); } } else if (psiName_ == "none") @@ -185,7 +181,7 @@ void Foam::uniformTotalPressureFvPatchScalarField::updateCoeffs const fvPatchField<scalar>& rho = patch().lookupPatchField<volScalarField, scalar>(rhoName_); - operator==(p0_ - 0.5*rho*(1.0 - pos(phip))*magSqr(Up)); + operator==(p0 - 0.5*rho*(1.0 - pos(phip))*magSqr(Up)); } else { @@ -219,7 +215,6 @@ void Foam::uniformTotalPressureFvPatchScalarField::write(Ostream& os) const os.writeKeyword("rho") << rhoName_ << token::END_STATEMENT << nl; os.writeKeyword("psi") << psiName_ << token::END_STATEMENT << nl; os.writeKeyword("gamma") << gamma_ << token::END_STATEMENT << nl; - os.writeKeyword("p0") << p0_ << token::END_STATEMENT << nl; pressure_->writeData(os); writeEntry("value", os); } diff --git a/src/finiteVolume/fields/fvPatchFields/derived/uniformTotalPressure/uniformTotalPressureFvPatchScalarField.H b/src/finiteVolume/fields/fvPatchFields/derived/uniformTotalPressure/uniformTotalPressureFvPatchScalarField.H index aa77ffe2ca4952ba3c4fb64a44d0f76fcd631e2f..325ed6e47191b8b74c9407c23ac6e663c47249fa 100644 --- a/src/finiteVolume/fields/fvPatchFields/derived/uniformTotalPressure/uniformTotalPressureFvPatchScalarField.H +++ b/src/finiteVolume/fields/fvPatchFields/derived/uniformTotalPressure/uniformTotalPressureFvPatchScalarField.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 @@ -75,9 +75,6 @@ class uniformTotalPressureFvPatchScalarField //- Heat capacity ratio scalar gamma_; - //- Total pressure - scalar p0_; - //- Table of time vs total pressure, including the bounding treatment autoPtr<DataEntry<scalar> > pressure_; @@ -178,18 +175,6 @@ public: return gamma_; } - //- Return the total pressure - scalar p0() const - { - return p0_; - } - - //- Return reference to the total pressure to allow adjustment - scalar p0() - { - return p0_; - } - // Evaluation functions diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H index 321fba33ad63dd0c4dc09975a6ea8f097f7e9bdd..29e696dc81a734dd6b14150203d9dc6ff6425bc1 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 a6bdf47cc46e750f79091d14201fad830c9e4a2a..2e793ac71e3ff9211c6761bc1176c17c1f2fbcaa 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<CloudType>::Dmax() const reduce(d, maxOp<scalar>()); - return d; + return max(0.0, d); } template<class CloudType> inline Foam::scalar Foam::KinematicCloud<CloudType>::penetration ( - const scalar& prc + const scalar& fraction ) const { - scalar distance = 0.0; - scalar mTot = 0.0; + if ((fraction < 0) || (fraction > 1)) + { + FatalErrorIn + ( + "inline Foam::scalar Foam::KinematicCloud<CloudType>::penetration" + "(" + "const scalar&" + ") const" + ) + << "fraction should be in the range 0 < fraction < 1" + << exit(FatalError); + } - label np = this->size(); + scalar distance = 0.0; - // arrays containing the parcels mass and - // distance from injector in ascending order - scalarField mass(np); - scalarField dist(np); + const label nParcel = this->size(); + globalIndex globalParcels(nParcel); + const label nParcelSum = globalParcels.size(); - if (np > 0) + if (nParcelSum == 0) { - label n = 0; + return distance; + } - // 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<CloudType>, *this, iter) - { - const parcelType& p = iter(); - scalar mi = p.nParticle()*p.mass(); - scalar di = mag(p.position() - p.position0()); - mTot += mi; + // lists of parcels mass and distance from initial injection point + List<scalar> mass(nParcel, 0.0); + List<scalar> dist(nParcel, 0.0); - // insert at the last place - mass[n] = mi; - dist[n] = di; + label i = 0; + scalar mSum = 0.0; + forAllConstIter(typename KinematicCloud<CloudType>, *this, iter) + { + const parcelType& p = iter(); + scalar m = p.nParticle()*p.mass(); + scalar d = mag(p.position() - p.position0()); + mSum += m; - label i = 0; - bool found = false; + mass[i] = m; + dist[i] = d; - // 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++; - } + i++; } - reduce(mTot, sumOp<scalar>()); + // calculate total mass across all processors + reduce(mSum, sumOp<scalar>()); + + // flatten the mass list + List<scalar> allMass(nParcelSum, 0.0); + SubList<scalar> + ( + allMass, + globalParcels.localSize(Pstream::myProcNo()), + globalParcels.offset(Pstream::myProcNo()) + ).assign(mass); + + // flatten the distance list + SortableList<scalar> allDist(nParcelSum, 0.0); + SubList<scalar> + ( + 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 (np > 0) + if (nParcelSum > 1) { - scalar mLimit = prc*mTot; - scalar mOff = (1.0 - prc)*mTot; + const scalar mLimit = fraction*mSum; + const labelList& indices = allDist.indices(); - if (np > 1) + 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 + distance = allDist.last(); + } + else + { + // 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) { - scalar mOffSum = 0.0; - label i = np; + label indI = indices[i]; + + mCurrent += allMass[indI]; - while ((mOffSum < mOff) && (i>0)) + if (mCurrent > mThreshold) { - i--; - mOffSum += mass[i]; + i0 = i; + break; } - distance = - dist[i+1] - + (dist[i] - dist[i+1])*(mOffSum - mOff) - /mass[i+1] ; } - } - else - { - distance = dist[0]; + + 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<scalar>()); + 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 01c80c16e21ee5a1bc8c8a518b8a7addeb37e6a0..d7a0fcf0d03056df7411dac483549f62240cbca0 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/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C index 84942b034932601a40b18b1899c71203422c5fc6..461ec5b285c6e4eac516534ca1eb7180cdcf7f15 100644 --- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C +++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C @@ -636,7 +636,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions *(sum(dMassSRGas) + sum(dMassSRLiquid) + sum(dMassSRSolid)) ); - const scalar xsi = min(T/5000.0, 1.0); + const scalar xsi = min(T/td.cloud().constProps().TMax(), 1.0); const scalar coeff = (1.0 - xsi*xsi)*td.cloud().constProps().hRetentionCoeff(); diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/PressureGradient/PressureGradientForce.C b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/PressureGradient/PressureGradientForce.C index a1d121ee2990e0bbd078f8f882f22de9e1d99dc8..c9d61f6d34d606eb2f72618b435f17cea7bcdaa2 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<CloudType>::PressureGradientForce ( CloudType& owner, const fvMesh& mesh, - const dictionary& dict + const dictionary& dict, + const word& forceType ) : - ParticleForce<CloudType>(owner, mesh, dict, typeName, true), - UName_(this->coeffs().lookup("U")), - gradUPtr_(NULL) + ParticleForce<CloudType>(owner, mesh, dict, forceType, true), + UName_(this->coeffs().template lookupOrDefault<word>("U", "U")), + DUcDtInterpPtr_(NULL) {} @@ -50,7 +52,7 @@ Foam::PressureGradientForce<CloudType>::PressureGradientForce : ParticleForce<CloudType>(pgf), UName_(pgf.UName_), - gradUPtr_(NULL) + DUcDtInterpPtr_(NULL) {} @@ -66,18 +68,48 @@ Foam::PressureGradientForce<CloudType>::~PressureGradientForce() template<class CloudType> void Foam::PressureGradientForce<CloudType>::cacheFields(const bool store) { + static word fName("DUcDt"); + + bool fieldExists = this->mesh().template foundObject<volVectorField>(fName); + if (store) { - const volVectorField& U = this->mesh().template - lookupObject<volVectorField>(UName_); - gradUPtr_ = fvc::grad(U).ptr(); + if (!fieldExists) + { + const volVectorField& Uc = this->mesh().template + lookupObject<volVectorField>(UName_); + + volVectorField* DUcDtPtr = new volVectorField + ( + fName, + fvc::ddt(Uc) + (Uc & fvc::grad(Uc)) + ); + + DUcDtPtr->store(); + } + + const volVectorField& DUcDt = this->mesh().template + lookupObject<volVectorField>(fName); + + DUcDtInterpPtr_.reset + ( + interpolation<vector>::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<volVectorField>(fName); + + const_cast<volVectorField&>(DUcDt).checkOut(); } } } @@ -95,11 +127,24 @@ Foam::forceSuSp Foam::PressureGradientForce<CloudType>::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<class CloudType> +Foam::scalar Foam::PressureGradientForce<CloudType>::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 663723a82093a8bf9526af9e84a227d9a8c0d2ff..b1986348fba94f22f64ada538c6008935d6032b0 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<CloudType> { - // 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<interpolation<vector> > 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<vector>& 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 c9bd24c62c66a0df23ebf75ad245fc3e2c571b52..6c085241dedf4289c9cc7b183141a57a0e401e8e 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<class CloudType> -const Foam::volTensorField& Foam::PressureGradientForce<CloudType>::gradU() -const +inline const Foam::interpolation<Foam::vector>& +Foam::PressureGradientForce<CloudType>::DUcDtInterp() const { - if (gradUPtr_) - { - return *gradUPtr_; - } - else + if (!DUcDtInterpPtr_.valid()) { FatalErrorIn ( - "const volTensorField& PressureGradientForce<CloudType>::gradU()" - "const" - ) << "gradU field not allocated" << abort(FatalError); - - return *reinterpret_cast<const volTensorField*>(0); + "inline const Foam::interpolation<Foam::vector>&" + "Foam::PressureGradientForce<CloudType>::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 3dca602dd443f59548e4a671fd8686b827ae17dc..c701e07a00d17d64de5fe69c9b0d7f8bbfb19ba6 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<CloudType>::VirtualMassForce ( CloudType& owner, const fvMesh& mesh, - const dictionary& dict + const dictionary& dict, + const word& forceType ) : - ParticleForce<CloudType>(owner, mesh, dict, typeName, true), - UName_(this->coeffs().template lookupOrDefault<word>("U", "U")), - Cvm_(readScalar(this->coeffs().lookup("Cvm"))), - DUcDtPtr_(NULL), - DUcDtInterpPtr_(NULL) + PressureGradientForce<CloudType>(owner, mesh, dict, forceType), + Cvm_(readScalar(this->coeffs().lookup("Cvm"))) {} @@ -51,11 +47,8 @@ Foam::VirtualMassForce<CloudType>::VirtualMassForce const VirtualMassForce& vmf ) : - ParticleForce<CloudType>(vmf), - UName_(vmf.UName_), - Cvm_(vmf.Cvm_), - DUcDtPtr_(NULL), - DUcDtInterpPtr_(NULL) + PressureGradientForce<CloudType>(vmf), + Cvm_(vmf.Cvm_) {} @@ -71,36 +64,7 @@ Foam::VirtualMassForce<CloudType>::~VirtualMassForce() template<class CloudType> void Foam::VirtualMassForce<CloudType>::cacheFields(const bool store) { - if (store && !DUcDtPtr_) - { - const volVectorField& Uc = this->mesh().template - lookupObject<volVectorField>(UName_); - - DUcDtPtr_ = new volVectorField - ( - "DUcDt", - fvc::ddt(Uc) + (Uc & fvc::grad(Uc)) - ); - - DUcDtInterpPtr_.reset - ( - interpolation<vector>::New - ( - this->owner().solution().interpolationSchemes(), - *DUcDtPtr_ - ).ptr() - ); - } - else - { - DUcDtInterpPtr_.clear(); - - if (DUcDtPtr_) - { - delete DUcDtPtr_; - DUcDtPtr_ = NULL; - } - } + PressureGradientForce<CloudType>::cacheFields(store); } @@ -114,12 +78,10 @@ Foam::forceSuSp Foam::VirtualMassForce<CloudType>::calcCoupled const scalar muc ) const { - forceSuSp value(vector::zero, 0.0); + forceSuSp value = + PressureGradientForce<CloudType>::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 c5277ec6963a7b09c0b20a91ef575d4f1f22ad17..ea0d5809324da2d5bd2024ea6a40e34aafa5c527 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 CloudType> class VirtualMassForce : - public ParticleForce<CloudType> + public PressureGradientForce<CloudType> { // 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<interpolation<vector> > 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<vector>& 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 e156c58d141b23e834db6602b3ac96564dd07493..0000000000000000000000000000000000000000 --- 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 <http://www.gnu.org/licenses/>. - -\*---------------------------------------------------------------------------*/ - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -template<class CloudType> -const Foam::volVectorField& Foam::VirtualMassForce<CloudType>::DUcDt() const -{ - if (DUcDtPtr_) - { - return *DUcDtPtr_; - } - else - { - FatalErrorIn - ( - "const volVectorField& VirtualMassForce<CloudType>::DUcDt()" - "const" - ) << "DUcDt field not allocated" << abort(FatalError); - - return *reinterpret_cast<const volVectorField*>(0); - } -} - - -template<class CloudType> -inline const Foam::interpolation<Foam::vector>& -Foam::VirtualMassForce<CloudType>::DUcDtInterp() const -{ - if (!DUcDtInterpPtr_.valid()) - { - FatalErrorIn - ( - "inline const Foam::interpolation<Foam::vector>&" - "Foam::VirtualMassForce<CloudType>::DUcDtInterp() const" - ) << "Carrier pahase DUcDt interpolation object not set" - << abort(FatalError); - } - - return DUcDtInterpPtr_(); -} - - -// ************************************************************************* // diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C index abe3bc487c336d67a38b8b3c213efffb449e37d4..26b4fb0187a26eae754ec845156f3f32c863cb28 100644 --- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C +++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C @@ -48,6 +48,10 @@ Description #include "globalIndex.H" #include "DynamicField.H" #include "PatchTools.H" +#include "slipPointPatchFields.H" +#include "fixedValuePointPatchFields.H" +#include "calculatedPointPatchFields.H" +#include "cyclicSlipPointPatchFields.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -71,7 +75,7 @@ void Foam::autoLayerDriver::dumpDisplacement ) { OFstream dispStr(prefix + "_disp.obj"); - Info<< "Writing all displacements to " << dispStr.name() << nl << endl; + Info<< "Writing all displacements to " << dispStr.name() << endl; label vertI = 0; @@ -87,7 +91,7 @@ void Foam::autoLayerDriver::dumpDisplacement OFstream illStr(prefix + "_illegal.obj"); - Info<< "Writing invalid displacements to " << illStr.name() << nl << endl; + Info<< "Writing invalid displacements to " << illStr.name() << endl; vertI = 0; @@ -801,6 +805,81 @@ void Foam::autoLayerDriver::setNumLayers } +// Construct pointVectorField with correct boundary conditions for adding +// layers +Foam::tmp<Foam::pointVectorField> +Foam::autoLayerDriver::makeLayerDisplacementField +( + const pointMesh& pMesh, + const labelList& numLayers +) +{ + // Construct displacement field. + const pointBoundaryMesh& pointPatches = pMesh.boundary(); + + wordList patchFieldTypes + ( + pointPatches.size(), + slipPointPatchVectorField::typeName + ); + + forAll(numLayers, patchI) + { + // 0 layers: do not allow lslip so fixedValue 0 + // >0 layers: fixedValue which gets adapted + if (numLayers[patchI] >= 0) + { + patchFieldTypes[patchI] = fixedValuePointPatchVectorField::typeName; + } + } + + forAll(pointPatches, patchI) + { + if (isA<processorPointPatch>(pointPatches[patchI])) + { + patchFieldTypes[patchI] = calculatedPointPatchVectorField::typeName; + } + else if (isA<cyclicPointPatch>(pointPatches[patchI])) + { + patchFieldTypes[patchI] = cyclicSlipPointPatchVectorField::typeName; + } + } + + +//Pout<< "*** makeLayerDisplacementField : boundary conditions:" << endl; +//forAll(patchFieldTypes, patchI) +//{ +// Pout<< "\t" << patchI << " name:" << pointPatches[patchI].name() +// << " type:" << patchFieldTypes[patchI] +// << " nLayers:" << numLayers[patchI] +// << endl; +//} + + const polyMesh& mesh = pMesh(); + + // Note: time().timeName() instead of meshRefinement::timeName() since + // postprocessable field. + tmp<pointVectorField> tfld + ( + new pointVectorField + ( + IOobject + ( + "pointDisplacement", + mesh.time().timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + pMesh, + dimensionedVector("displacement", dimLength, vector::zero), + patchFieldTypes + ) + ); + return tfld; +} + + void Foam::autoLayerDriver::growNoExtrusion ( const indirectPrimitivePatch& pp, @@ -2392,10 +2471,10 @@ void Foam::autoLayerDriver::addLayers mesh, pp(), patchIDs, - meshRefinement::makeDisplacementField + makeLayerDisplacementField ( pointMesh::New(mesh), - patchIDs + layerParams.numLayers() ), motionDict ) @@ -3186,7 +3265,7 @@ void Foam::autoLayerDriver::doLayers // Merge coplanar boundary faces mergePatchFacesUndo(layerParams, motionDict); - // Per patch the number of layers (0 if no layer) + // Per patch the number of layers (-1 or 0 if no layer) const labelList& numLayers = layerParams.numLayers(); // Patches that need to get a layer diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H index cb720dc64fe2de82cd000036181d622dd4a86116..d36d15a4cb15fe0138f8e087ce0b31b91b28c1a5 100644 --- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H +++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H @@ -195,6 +195,19 @@ class autoLayerDriver List<extrudeMode>& extrudeStatus ) const; + //- Helper function to make a pointVectorField with correct + // bcs for layer addition: + // - numLayers > 0 : fixedValue + // - numLayers == 0 : fixedValue (always zero) + // - processor : calculated (so free to move) + // - cyclic/wedge/symmetry : slip + // - other : slip + static tmp<pointVectorField> makeLayerDisplacementField + ( + const pointMesh& pMesh, + const labelList& numLayers + ); + //- Grow no-extrusion layer. void growNoExtrusion ( @@ -444,7 +457,6 @@ class autoLayerDriver const PackedBoolList& isMasterEdge, const labelList& meshEdges, const scalar minCosLayerTermination, - scalarField& field, List<extrudeMode>& extrudeStatus, pointField& patchDisp, labelList& patchNLayers diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C index 35133f062d10f9ea37ef43f81c69058c334fa8be..58725470e1363a641c00da979c5b0c9f2c39887f 100644 --- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C +++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C @@ -99,6 +99,76 @@ 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<label>()) +// << endl; +// } +// } +//} +//XXXXXXXXX void Foam::autoLayerDriver::smoothField ( const motionSmoother& meshMover, @@ -126,9 +196,15 @@ void Foam::autoLayerDriver::smoothField // Get smoothly varying patch field. Info<< "shrinkMeshDistance : Smoothing field ..." << endl; - for (label iter = 0; iter < nSmoothDisp; iter++) + + const scalar lambda = 0.33; + const scalar mu = -0.34; + + for (label iter = 0; iter < 90; iter++) { scalarField average(pp.nPoints()); + + // Calculate average of field averageNeighbours ( meshMover.mesh(), @@ -141,23 +217,37 @@ void Foam::autoLayerDriver::smoothField average ); - // Transfer to field - forAll(field, pointI) + forAll(field, i) { - //full smoothing neighbours + point value - average[pointI] = 0.5*(field[pointI]+average[pointI]); + if (field[i] >= fieldMin[i]) + { + field[i] = (1-lambda)*field[i]+lambda*average[i]; + } + } - // perform monotonic smoothing - if - ( - average[pointI] < field[pointI] - && average[pointI] >= fieldMin[pointI] - ) + + // Calculate average of field + averageNeighbours + ( + meshMover.mesh(), + isMasterEdge, + meshEdges, + meshPoints, + pp.edges(), + invSumWeight, + field, + average + ); + + forAll(field, i) + { + if (field[i] >= fieldMin[i]) { - field[pointI] = average[pointI]; + field[i] = (1-mu)*field[i]+mu*average[i]; } } + // Do residual calculation every so often. if ((iter % 10) == 0) { @@ -168,7 +258,7 @@ void Foam::autoLayerDriver::smoothField } } } - +//XXXXXXXXX // Smooth normals on moving patch. void Foam::autoLayerDriver::smoothPatchNormals @@ -480,7 +570,6 @@ void Foam::autoLayerDriver::findIsolatedRegions const PackedBoolList& isMasterEdge, const labelList& meshEdges, const scalar minCosLayerTermination, - scalarField& field, List<extrudeMode>& extrudeStatus, pointField& patchDisp, labelList& patchNLayers @@ -661,7 +750,6 @@ void Foam::autoLayerDriver::findIsolatedRegions ) ) { - field[f[fp]] = 0.0; nPointCounter++; } } @@ -670,7 +758,8 @@ void Foam::autoLayerDriver::findIsolatedRegions } reduce(nPointCounter, sumOp<label>()); - Info<< "Number isolated points extrusion stopped : "<< nPointCounter<< endl; + Info<< "Number isolated points extrusion stopped : "<< nPointCounter + << endl; } @@ -875,11 +964,14 @@ void Foam::autoLayerDriver::medialAxisSmoothingInfo forAll(patches, patchI) { const polyPatch& pp = patches[patchI]; + //const pointPatchVectorField& pvf = + // meshMover.displacement().boundaryField()[patchI]; if ( !pp.coupled() && !isA<emptyPolyPatch>(pp) + //&& pvf.constraintType() != word::null //exclude fixedValue && !adaptPatches.found(patchI) ) { @@ -1168,12 +1260,22 @@ void Foam::autoLayerDriver::shrinkMeshMedialDistance isMasterEdge, meshEdges, minCosLayerTermination, - thickness, + extrudeStatus, patchDisp, patchNLayers ); + // Update thickess for changed extrusion + forAll(thickness, patchPointI) + { + if (extrudeStatus[patchPointI] == NOEXTRUDE) + { + thickness[patchPointI] = 0.0; + } + } + + // smooth layer thickness on moving patch smoothField ( @@ -1182,6 +1284,7 @@ void Foam::autoLayerDriver::shrinkMeshMedialDistance meshEdges, minThickness, nSmoothThickness, + thickness ); diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C index 7f612246f42c33c395a5d5b4e678b1dd0de90ec9..edbd6a04742889fd0930ebf80adaca5299355661 100644 --- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C +++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C @@ -506,7 +506,7 @@ void Foam::autoSnapDriver::dumpMove ) { // Dump direction of growth into file - Pout<< nl << "Dumping move direction to " << fName << endl; + Info<< "Dumping move direction to " << fName << endl; OFstream nearestStream(fName); @@ -721,14 +721,14 @@ void Foam::autoSnapDriver::preSmoothPatch if (debug) { const_cast<Time&>(mesh.time())++; - Pout<< "Writing patch smoothed mesh to time " + Info<< "Writing patch smoothed mesh to time " << meshRefiner_.timeName() << '.' << endl; meshRefiner_.write ( debug, mesh.time().path()/meshRefiner_.timeName() ); - Pout<< "Dumped mesh in = " + Info<< "Dumped mesh in = " << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl; } @@ -998,7 +998,7 @@ void Foam::autoSnapDriver::smoothDisplacement if (debug) { const_cast<Time&>(mesh.time())++; - Pout<< "Writing smoothed mesh to time " << meshRefiner_.timeName() + Info<< "Writing smoothed mesh to time " << meshRefiner_.timeName() << endl; // Moving mesh creates meshPhi. Can be cleared out by a mesh.clearOut @@ -1010,13 +1010,12 @@ void Foam::autoSnapDriver::smoothDisplacement debug, mesh.time().path()/meshRefiner_.timeName() ); - - Pout<< "Writing displacement field ..." << endl; + Info<< "Writing displacement field ..." << endl; disp.write(); tmp<pointScalarField> magDisp(mag(disp)); magDisp().write(); - Pout<< "Writing actual patch displacement ..." << endl; + Info<< "Writing actual patch displacement ..." << endl; vectorField actualPatchDisp(disp, pp.meshPoints()); dumpMove ( @@ -1068,11 +1067,11 @@ bool Foam::autoSnapDriver::scaleMesh if (debug) { const_cast<Time&>(mesh.time())++; - Pout<< "Writing scaled mesh to time " << meshRefiner_.timeName() + Info<< "Writing scaled mesh to time " << meshRefiner_.timeName() << endl; mesh.write(); - Pout<< "Writing displacement field ..." << endl; + Info<< "Writing displacement field ..." << endl; meshMover.displacement().write(); tmp<pointScalarField> magDisp(mag(meshMover.displacement())); magDisp().write(); @@ -1521,14 +1520,14 @@ void Foam::autoSnapDriver::doSnap if (debug) { const_cast<Time&>(mesh.time())++; - Pout<< "Writing scaled mesh to time " + Info<< "Writing scaled mesh to time " << meshRefiner_.timeName() << endl; meshRefiner_.write ( debug, mesh.time().path()/meshRefiner_.timeName() ); - Pout<< "Writing displacement field ..." << endl; + Info<< "Writing displacement field ..." << endl; meshMover.displacement().write(); tmp<pointScalarField> magDisp(mag(meshMover.displacement())); magDisp().write(); @@ -1564,7 +1563,7 @@ void Foam::autoSnapDriver::doSnap if (nChanged > 0 && debug) { const_cast<Time&>(mesh.time())++; - Pout<< "Writing patchFace merged mesh to time " + Info<< "Writing patchFace merged mesh to time " << meshRefiner_.timeName() << endl; meshRefiner_.write ( diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C index e526f8faf96546aaaa5e133e7d3a813f22d31dcc..88ec55b1c8b452ce9f2cb3cb25ac272f72101cf3 100644 --- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C +++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C @@ -2565,6 +2565,76 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature } } + + // Collect additionally 'normal' boundary faces for boundaryPoints of pp + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + // points on the boundary of pp should pick up non-pp normals + // as well for the feature-reconstruction to behave correctly. + // (the movement is already constrained outside correctly so it + // is only that the unconstrained attraction vector is calculated + // correctly) + { + const polyBoundaryMesh& pbm = mesh.boundaryMesh(); + labelList patchID(pbm.patchID()); + + // Unmark all non-coupled boundary faces + forAll(pbm, patchI) + { + const polyPatch& pp = pbm[patchI]; + + if (pp.coupled() || isA<emptyPolyPatch>(pp)) + { + forAll(pp, i) + { + label meshFaceI = pp.start()+i; + patchID[meshFaceI-mesh.nInternalFaces()] = -1; + } + } + } + + // Remove any meshed faces + PackedBoolList ppFaces(mesh.nFaces()); + forAll(pp.addressing(), i) + { + label meshFaceI = pp.addressing()[i]; + patchID[meshFaceI-mesh.nInternalFaces()] = -1; + } + + // See if pp point uses any non-meshed boundary faces + + const labelList& boundaryPoints = pp.boundaryPoints(); + forAll(boundaryPoints, i) + { + label pointI = boundaryPoints[i]; + label meshPointI = pp.meshPoints()[pointI]; + const point& pt = mesh.points()[meshPointI]; + const labelList& pFaces = mesh.pointFaces()[meshPointI]; + + List<point>& pNormals = pointFaceSurfNormals[pointI]; + List<point>& pDisp = pointFaceDisp[pointI]; + List<point>& pFc = pointFaceCentres[pointI]; + labelList& pFid = pointFacePatchID[pointI]; + + forAll(pFaces, i) + { + label meshFaceI = pFaces[i]; + if (!mesh.isInternalFace(meshFaceI)) + { + label patchI = patchID[meshFaceI-mesh.nInternalFaces()]; + + if (patchI != -1) + { + vector fn = mesh.faceAreas()[meshFaceI]; + pNormals.append(fn/mag(fn)); + pDisp.append(mesh.faceCentres()[meshFaceI]-pt); + pFc.append(mesh.faceCentres()[meshFaceI]); + pFid.append(patchI); + } + } + } + } + } + syncTools::syncPointList ( mesh, diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C index 95d9bf0dd467ba9269b4759605e5fb869cd087a0..cc466df4be1ac0385c8c6905b03d7829b9a6f409 100644 --- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C +++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C @@ -43,7 +43,7 @@ Foam::layerParameters::layerParameters const polyBoundaryMesh& boundaryMesh ) : - numLayers_(boundaryMesh.size(), 0), + numLayers_(boundaryMesh.size(), -1), expansionRatio_ ( boundaryMesh.size(), diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.H index 1a7913a8efd8fd6788a6b034e78f7d86c1a0f232..4a48784ee9821e9157cf9b3f197628e27f0f7a9a 100644 --- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.H +++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.H @@ -131,6 +131,10 @@ public: // Per patch information //- How many layers to add. + // -1 : no specification. Assume 0 layers but allow sliding + // to make layers + // 0 : specified to have 0 layers. No sliding allowed. + // >0 : number of layers const labelList& numLayers() const { return numLayers_; diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C index f2ca0857dcd992e9268b85fe60e20a613d8c0156..5561b21829417d528a2d5762b4bbdc32e2506ad7 100644 --- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C +++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C @@ -571,7 +571,7 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells hitInfo ); - // Start of from current points + // Start off from current points newPoints = mesh_.points(); forAll(hitInfo, i) @@ -581,6 +581,17 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells newPoints[meshPoints[i]] = hitInfo[i].hitPoint(); } } + + if (debug) + { + const_cast<Time&>(mesh_.time())++; + pointField oldPoints(mesh_.points()); + mesh_.movePoints(newPoints); + Pout<< "Writing newPoints mesh to time " << timeName() + << endl; + write(debug, mesh_.time().path()/"newPoints"); + mesh_.movePoints(oldPoints); + } } diff --git a/src/turbulenceModels/compressible/LES/GenEddyVisc/GenEddyVisc.H b/src/turbulenceModels/compressible/LES/GenEddyVisc/GenEddyVisc.H index 5f19ee925b224157eb39d8ffe807f2533f5fa9ab..7e32d57396f99f0ea751223df915dcb2e3a63453 100644 --- a/src/turbulenceModels/compressible/LES/GenEddyVisc/GenEddyVisc.H +++ b/src/turbulenceModels/compressible/LES/GenEddyVisc/GenEddyVisc.H @@ -51,7 +51,7 @@ namespace LESModels { /*---------------------------------------------------------------------------*\ - Class GenEddyVisc Declaration + Class GenEddyVisc Declaration \*---------------------------------------------------------------------------*/ class GenEddyVisc @@ -108,7 +108,21 @@ public: //- Return sub-grid disipation rate virtual tmp<volScalarField> epsilon() const { - return ce_*k()*sqrt(k())/delta(); + return tmp<volScalarField> + ( + new volScalarField + ( + IOobject + ( + "epsilon", + runTime_.timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + ce_*k()*sqrt(k())/delta() + ) + ); } //- Return viscosity diff --git a/src/turbulenceModels/compressible/LES/GenSGSStress/GenSGSStress.H b/src/turbulenceModels/compressible/LES/GenSGSStress/GenSGSStress.H index e69d2253c9cea305747f103d4d5ffe850852667c..202a85371ead0719fad8e24d40b1ddfca8baec53 100644 --- a/src/turbulenceModels/compressible/LES/GenSGSStress/GenSGSStress.H +++ b/src/turbulenceModels/compressible/LES/GenSGSStress/GenSGSStress.H @@ -52,7 +52,7 @@ namespace LESModels { /*---------------------------------------------------------------------------*\ - Class GenSGSStress Declaration + Class GenSGSStress Declaration \*---------------------------------------------------------------------------*/ class GenSGSStress @@ -109,14 +109,43 @@ public: //- Return the SGS turbulent kinetic energy virtual tmp<volScalarField> k() const { - return 0.5*tr(B_); + return tmp<volScalarField> + ( + new volScalarField + ( + IOobject + ( + "k", + runTime_.timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + 0.5*tr(B_) + ) + ); } //- Return the SGS turbulent dissipation virtual tmp<volScalarField> epsilon() const { const volScalarField K(k()); - return ce_*K*sqrt(K)/delta(); + + return tmp<volScalarField> + ( + new volScalarField + ( + IOobject + ( + "epsilon", + runTime_.timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + ce_*K*sqrt(K)/delta() + ) + ); } //- Return the SGS viscosity diff --git a/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/turbulentHeatFluxTemperature/turbulentHeatFluxTemperatureFvPatchScalarField.C b/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/turbulentHeatFluxTemperature/turbulentHeatFluxTemperatureFvPatchScalarField.C index 29757f1ea6eb381255799639c79c6c4352cc54a0..fbae992c25fd4f773a95087dfa40ca4fe3bef3e7 100644 --- a/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/turbulentHeatFluxTemperature/turbulentHeatFluxTemperatureFvPatchScalarField.C +++ b/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/turbulentHeatFluxTemperature/turbulentHeatFluxTemperatureFvPatchScalarField.C @@ -79,7 +79,8 @@ turbulentHeatFluxTemperatureFvPatchScalarField fixedGradientFvPatchScalarField(p, iF), temperatureCoupledBase(patch(), "undefined", "undefined-K"), heatSource_(hsPower), - q_(p.size(), 0.0) + q_(p.size(), 0.0), + QrName_("undefinedQr") {} @@ -95,7 +96,8 @@ turbulentHeatFluxTemperatureFvPatchScalarField fixedGradientFvPatchScalarField(ptf, p, iF, mapper), temperatureCoupledBase(patch(), ptf.KMethod(), ptf.kappaName()), heatSource_(ptf.heatSource_), - q_(ptf.q_, mapper) + q_(ptf.q_, mapper), + QrName_(ptf.QrName_) {} @@ -110,7 +112,8 @@ turbulentHeatFluxTemperatureFvPatchScalarField fixedGradientFvPatchScalarField(p, iF), temperatureCoupledBase(patch(), dict), heatSource_(heatSourceTypeNames_.read(dict.lookup("heatSource"))), - q_("q", dict, p.size()) + q_("q", dict, p.size()), + QrName_(dict.lookupOrDefault<word>("Qr", "none")) { fvPatchField<scalar>::operator=(patchInternalField()); gradient() = 0.0; @@ -126,7 +129,8 @@ turbulentHeatFluxTemperatureFvPatchScalarField fixedGradientFvPatchScalarField(thftpsf), temperatureCoupledBase(patch(), thftpsf.KMethod(), thftpsf.kappaName()), heatSource_(thftpsf.heatSource_), - q_(thftpsf.q_) + q_(thftpsf.q_), + QrName_(thftpsf.QrName_) {} @@ -140,7 +144,8 @@ turbulentHeatFluxTemperatureFvPatchScalarField fixedGradientFvPatchScalarField(thftpsf, iF), temperatureCoupledBase(patch(), thftpsf.KMethod(), thftpsf.kappaName()), heatSource_(thftpsf.heatSource_), - q_(thftpsf.q_) + q_(thftpsf.q_), + QrName_(thftpsf.QrName_) {} @@ -183,17 +188,25 @@ void turbulentHeatFluxTemperatureFvPatchScalarField::updateCoeffs() const scalarField& Tp = *this; + scalarField qr(this->size(), 0.0); + + //- qr is negative going into the domain + if (QrName_ != "none") + { + qr = patch().lookupPatchField<volScalarField, scalar>(QrName_); + } + switch (heatSource_) { case hsPower: { const scalar Ap = gSum(patch().magSf()); - gradient() = q_/(Ap*kappa(Tp)); + gradient() = (q_/Ap + qr)/kappa(Tp); break; } case hsFlux: { - gradient() = q_/kappa(Tp); + gradient() = (q_ + qr)/kappa(Tp); break; } default: @@ -220,6 +233,7 @@ void turbulentHeatFluxTemperatureFvPatchScalarField::write << token::END_STATEMENT << nl; temperatureCoupledBase::write(os); q_.writeEntry("q", os); + os.writeKeyword("Qr")<< QrName_ << token::END_STATEMENT << nl; writeEntry("value", os); } diff --git a/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/turbulentHeatFluxTemperature/turbulentHeatFluxTemperatureFvPatchScalarField.H b/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/turbulentHeatFluxTemperature/turbulentHeatFluxTemperatureFvPatchScalarField.H index 254ed52b756f1ed0521f1d50b33514a6af9b6c89..740183e35902541cc29b212efcd879cdf71266ce 100644 --- a/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/turbulentHeatFluxTemperature/turbulentHeatFluxTemperatureFvPatchScalarField.H +++ b/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/turbulentHeatFluxTemperature/turbulentHeatFluxTemperatureFvPatchScalarField.H @@ -37,6 +37,7 @@ Description heatSource flux; // power [W]; flux [W/m2] q uniform 10; // heat power or flux kappa fluidThermo; // calculate kappa=alphaEff*thermo.Cp + Qr none; // name of the radiative flux value uniform 300; // initial temperature value } @@ -93,6 +94,9 @@ private: //- Heat power [W] or flux [W/m2] scalarField q_; + //- Name of radiative in flux field + word QrName_; + public: diff --git a/src/turbulenceModels/incompressible/LES/GenEddyVisc/GenEddyVisc.H b/src/turbulenceModels/incompressible/LES/GenEddyVisc/GenEddyVisc.H index 62c216ee3ec06fd13b84ef69c3baff5d2991d977..c40adc383be052813551a0ad60bf89b373ba27bb 100644 --- a/src/turbulenceModels/incompressible/LES/GenEddyVisc/GenEddyVisc.H +++ b/src/turbulenceModels/incompressible/LES/GenEddyVisc/GenEddyVisc.H @@ -51,7 +51,7 @@ namespace LESModels { /*---------------------------------------------------------------------------*\ - Class GenEddyVisc Declaration + Class GenEddyVisc Declaration \*---------------------------------------------------------------------------*/ class GenEddyVisc @@ -104,7 +104,21 @@ public: //- Return sub-grid disipation rate virtual tmp<volScalarField> epsilon() const { - return ce_*k()*sqrt(k())/delta(); + return tmp<volScalarField> + ( + new volScalarField + ( + IOobject + ( + "epsilon", + runTime_.timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + ce_*k()*sqrt(k())/delta() + ) + ); } //- Return the SGS viscosity. diff --git a/src/turbulenceModels/incompressible/LES/GenSGSStress/GenSGSStress.H b/src/turbulenceModels/incompressible/LES/GenSGSStress/GenSGSStress.H index 0456d69a7178193ffa6b0f05cea10f521796c247..f7b6040f377701637cffcb50a58d53a4a723960a 100644 --- a/src/turbulenceModels/incompressible/LES/GenSGSStress/GenSGSStress.H +++ b/src/turbulenceModels/incompressible/LES/GenSGSStress/GenSGSStress.H @@ -52,7 +52,7 @@ namespace LESModels { /*---------------------------------------------------------------------------*\ - Class GenSGSStress Declaration + Class GenSGSStress Declaration \*---------------------------------------------------------------------------*/ class GenSGSStress @@ -104,14 +104,43 @@ public: //- Return the SGS turbulent kinetic energy. virtual tmp<volScalarField> k() const { - return 0.5*tr(B_); + return tmp<volScalarField> + ( + new volScalarField + ( + IOobject + ( + "k", + runTime_.timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + 0.5*tr(B_) + ) + ); } //- Return the SGS turbulent dissipation. virtual tmp<volScalarField> epsilon() const { const volScalarField K(k()); - return ce_*K*sqrt(K)/delta(); + + return tmp<volScalarField> + ( + new volScalarField + ( + IOobject + ( + "epsilon", + runTime_.timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + ce_*K*sqrt(K)/delta() + ) + ); } //- Return the SGS viscosity. diff --git a/tutorials/mesh/cvMesh/flange/system/cvMeshDict b/tutorials/mesh/cvMesh/flange/system/cvMeshDict index cbd86c379d4f08a8008b6e76ad3b062b93305dd8..e0dace7c28388a91b23c8f0ad12e712d419d8594 100644 --- a/tutorials/mesh/cvMesh/flange/system/cvMeshDict +++ b/tutorials/mesh/cvMesh/flange/system/cvMeshDict @@ -217,7 +217,9 @@ polyMeshFiltering } -#include "meshQualityControls" - +meshQualityControls +{ + #include "meshQualityDict" +} // ************************************************************************* // diff --git a/tutorials/mesh/cvMesh/flange/system/meshQualityControls b/tutorials/mesh/cvMesh/flange/system/meshQualityControls deleted file mode 100644 index 09ebbb2458e42e23f5b07ed411ce74493f1d7d2b..0000000000000000000000000000000000000000 --- a/tutorials/mesh/cvMesh/flange/system/meshQualityControls +++ /dev/null @@ -1,76 +0,0 @@ -/*--------------------------------*- C++ -*----------------------------------*\ -| ========= | | -| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | -| \\ / O peration | Version: dev | -| \\ / A nd | Web: www.OpenFOAM.org | -| \\/ M anipulation | | -\*---------------------------------------------------------------------------*/ - -FoamFile -{ - version 2.0; - format ascii; - - root ""; - case ""; - instance ""; - local ""; - - class dictionary; - object meshQualityControls; -} - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -meshQualityControls -{ - //- Maximum non-orthogonality allowed. Set to 180 to disable. - maxNonOrtho 65; - - //- Max skewness allowed. Set to <0 to disable. - maxBoundarySkewness 50; - maxInternalSkewness 10; - - //- Max concaveness allowed. Is angle (in degrees) below which concavity - // is allowed. 0 is straight face, <0 would be convex face. - // Set to 180 to disable. - maxConcave 80; - - //- Minimum quality of the tet formed by the face-centre - // and variable base point minimum decomposition triangles and - // the cell centre. This has to be a positive number for tracking - // to work. Set to very negative number (e.g. -1E30) to - // disable. - // <0 = inside out tet, - // 0 = flat tet - // 1 = regular tet - minTetQuality 1e-30; - - //- Minimum pyramid volume. Is absolute volume of cell pyramid. - // Set to a sensible fraction of the smallest cell volume expected. - // Set to very negative number (e.g. -1E30) to disable. - minVol 0; - - //- Minimum face area. Set to <0 to disable. - minArea -1; - - //- Minimum face twist. Set to <-1 to disable. dot product of face normal - //- and face centre triangles normal - minTwist 0.001; - - //- minimum normalised cell determinant - //- 1 = hex, <= 0 = folded or flattened illegal cell - minDeterminant 0.001; - - //- minFaceWeight (0 -> 0.5) - minFaceWeight 0.02; - - //- minVolRatio (0 -> 1) - minVolRatio 0.01; - - //must be >0 for Fluent compatibility - minTriangleTwist -1; -} - - -// ************************************************************************* // diff --git a/tutorials/mesh/cvMesh/flange/system/meshQualityDict b/tutorials/mesh/cvMesh/flange/system/meshQualityDict new file mode 100644 index 0000000000000000000000000000000000000000..1b83e2bd87b647e97cf8b39d2dccbb30af267c8e --- /dev/null +++ b/tutorials/mesh/cvMesh/flange/system/meshQualityDict @@ -0,0 +1,73 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ + +FoamFile +{ + version 2.0; + format ascii; + + root ""; + case ""; + instance ""; + local ""; + + class dictionary; + object meshQualityDict; +} + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +//- Maximum non-orthogonality allowed. Set to 180 to disable. +maxNonOrtho 65; + +//- Max skewness allowed. Set to <0 to disable. +maxBoundarySkewness 50; +maxInternalSkewness 10; + +//- Max concaveness allowed. Is angle (in degrees) below which concavity +// is allowed. 0 is straight face, <0 would be convex face. +// Set to 180 to disable. +maxConcave 80; + +//- Minimum quality of the tet formed by the face-centre +// and variable base point minimum decomposition triangles and +// the cell centre. This has to be a positive number for tracking +// to work. Set to very negative number (e.g. -1E30) to +// disable. +// <0 = inside out tet, +// 0 = flat tet +// 1 = regular tet +minTetQuality 1e-30; + +//- Minimum pyramid volume. Is absolute volume of cell pyramid. +// Set to a sensible fraction of the smallest cell volume expected. +// Set to very negative number (e.g. -1E30) to disable. +minVol 0; + +//- Minimum face area. Set to <0 to disable. +minArea -1; + +//- Minimum face twist. Set to <-1 to disable. dot product of face normal +//- and face centre triangles normal +minTwist 0.001; + +//- minimum normalised cell determinant +//- 1 = hex, <= 0 = folded or flattened illegal cell +minDeterminant 0.001; + +//- minFaceWeight (0 -> 0.5) +minFaceWeight 0.02; + +//- minVolRatio (0 -> 1) +minVolRatio 0.01; + +//must be >0 for Fluent compatibility +minTriangleTwist -1; + + +// ************************************************************************* // diff --git a/tutorials/mesh/cvMesh/flange/system/snappyHexMeshDict b/tutorials/mesh/cvMesh/flange/system/snappyHexMeshDict index 034e3c0e749d95c36508975773bcddabbd49a150..9da7fabcafea5f7ef0e7951a190148200f5a701c 100644 --- a/tutorials/mesh/cvMesh/flange/system/snappyHexMeshDict +++ b/tutorials/mesh/cvMesh/flange/system/snappyHexMeshDict @@ -278,12 +278,12 @@ addLayersControls -// Generic mesh quality settings. At any undoable phase these determine -// where to undo. -#include "meshQualityControls" - meshQualityControls { + // Generic mesh quality settings. At any undoable phase these determine + // where to undo. + #include "meshQualityDict" + //- Number of error distribution iterations nSmoothScale 4; //- amount to scale back displacement at error points