拉格朗日求解器代码求助
-
各位大佬好，我想要在DPMFoam植入我的颗粒张大的RP方程，简单如下
这个里面需要对颗粒进行一些定义，所以我修改了 OpenFOAM/OpenFOAM-3.0.0/src/lagrangian/intermediate这个源文件库，再去进行我的求解器编译，但是编译求解器中我遇到了如下报错g++ -m64 -Dlinux64 -DWM_ARCH_OPTION=64 -DWM_DP -DWM_LABEL_SIZE=32 -Wall -Wextra -Wold-style-cast -Wnon-virtual-dtor -Wno-unused-parameter -Wno-invalid-offsetof -O3 -DNoRepository -ftemplate-depth-100 -I./DPMTurbulenceModels/lnInclude -I/home/zly/OpenFOAM/OpenFOAM-3.0.0/src/lagrangian/basic/lnInclude -I/home/zly/OpenFOAM/OpenFOAM-3.0.0/src/lagrangian/intermediate/lnInclude -I/home/zly/OpenFOAM/OpenFOAM-3.0.0/src/thermophysicalModels/specie/lnInclude -I/home/zly/OpenFOAM/OpenFOAM-3.0.0/src/transportModels/compressible/lnInclude -I/home/zly/OpenFOAM/OpenFOAM-3.0.0/src/thermophysicalModels/basic/lnInclude -I/home/zly/OpenFOAM/OpenFOAM-3.0.0/src/thermophysicalModels/reactionThermo/lnInclude -I/home/zly/OpenFOAM/OpenFOAM-3.0.0/src/thermophysicalModels/radiation/lnInclude -I/home/zly/OpenFOAM/OpenFOAM-3.0.0/src/transportModels -I/home/zly/OpenFOAM/OpenFOAM-3.0.0/src/transportModels/incompressible/singlePhaseTransportModel -I/home/zly/OpenFOAM/OpenFOAM-3.0.0/src/TurbulenceModels/turbulenceModels/lnInclude -I/home/zly/OpenFOAM/OpenFOAM-3.0.0/src/TurbulenceModels/incompressible/lnInclude -I/home/zly/OpenFOAM/OpenFOAM-3.0.0/src/TurbulenceModels/phaseIncompressible/lnInclude -I/home/zly/OpenFOAM/OpenFOAM-3.0.0/src/regionModels/regionModel/lnInclude -I/home/zly/OpenFOAM/OpenFOAM-3.0.0/src/regionModels/surfaceFilmModels/lnInclude -I/home/zly/OpenFOAM/OpenFOAM-3.0.0/src/finiteVolume/lnInclude -I/home/zly/OpenFOAM/OpenFOAM-3.0.0/src/meshTools/lnInclude -IlnInclude -I. -I/home/zly/OpenFOAM/OpenFOAM-3.0.0/src/OpenFOAM/lnInclude -I/home/zly/OpenFOAM/OpenFOAM-3.0.0/src/OSspecific/POSIX/lnInclude -fPIC -Xlinker --add-needed -Xlinker --no-as-needed Make/linux64GccDPInt32Opt/ppDPMFoam.o -L/home/zly/OpenFOAM/OpenFOAM-3.0.0/platforms/linux64GccDPInt32Opt/lib \ -llagrangian -llagrangianIntermediate -llagrangianTurbulence -lthermophysicalFunctions -lspecie -lradiationModels -lincompressibleTransportModels -lturbulenceModels -lincompressibleTurbulenceModels -lDPMTurbulenceModels -lregionModels -lsurfaceFilmModels -lsampling -lfiniteVolume -lmeshTools -lOpenFOAM -ldl \ -lm -o /home/zly/OpenFOAM/zly-3.0.0/platforms/linux64GccDPInt32Opt/bin/ppDPMFoam Make/linux64GccDPInt32Opt/ppDPMFoam.o: In function `Foam::KinematicCloud<Foam::Cloud<Foam::CollidingParcel<Foam::KinematicParcel<Foam::particle> > > >::setParcelThermoProperties(Foam::CollidingParcel<Foam::KinematicParcel<Foam::particle> >&, double) [clone .isra.656]': ppDPMFoam.C:(.text+0x1377): undefined reference to `Foam::KinematicParcel<Foam::particle>::Ro()' ppDPMFoam.C:(.text+0x13c9): undefined reference to `Foam::KinematicParcel<Foam::particle>::pgo()' Make/linux64GccDPInt32Opt/ppDPMFoam.o: In function `void Foam::KinematicParcel<Foam::particle>::setCellValues<Foam::KinematicParcel<Foam::particle>::TrackingData<Foam::CollidingCloud<Foam::KinematicCloud<Foam::Cloud<Foam::CollidingParcel<Foam::KinematicParcel<Foam::particle> > > > > > >(Foam::KinematicParcel<Foam::particle>::TrackingData<Foam::CollidingCloud<Foam::KinematicCloud<Foam::Cloud<Foam::CollidingParcel<Foam::KinematicParcel<Foam::particle> > > > > >&, double, int)': ppDPMFoam.C:(.text._ZN4Foam15KinematicParcelINS_8particleEE13setCellValuesINS2_12TrackingDataINS_14CollidingCloudINS_14KinematicCloudINS_5CloudINS_15CollidingParcelIS2_EEEEEEEEEEEEvRT_di[_ZN4Foam15KinematicParcelINS_8particleEE13setCellValuesINS2_12TrackingDataINS_14CollidingCloudINS_14KinematicCloudINS_5CloudINS_15CollidingParcelIS2_EEEEEEEEEEEEvRT_di]+0x23c): undefined reference to `Foam::KinematicParcel<Foam::particle>::TrackingData<Foam::CollidingCloud<Foam::KinematicCloud<Foam::Cloud<Foam::CollidingParcel<Foam::KinematicParcel<Foam::particle> > > > > >::PInterp() const' collect2: error: ld returned 1 exit status make: *** [/home/zly/OpenFOAM/zly-3.0.0/platforms/linux64GccDPInt32Opt/bin/ppDPMFoam] Error 1似乎在说我的Ro，pgo和PInterp没有定义，其中第一个是初始的气泡半径，第二个是初始气泡内压，第三个是颗粒受到的压强插值。
但是我明明在KinematicParcel.H和.C文件进行了定义，如下其中//-myadd即为我的自己添加)
KinematicParcel.H\*---------------------------------------------------------------------------*/ #ifndef KinematicParcel_H #define KinematicParcel_H #include "particle.H" #include "IOstream.H" #include "autoPtr.H" #include "interpolation.H" #include "demandDrivenEntry.H" // #include "ParticleForceList.H" // TODO // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // namespace Foam { template<class ParcelType> class KinematicParcel; // Forward declaration of friend functions template<class ParcelType> Ostream& operator<< ( Ostream&, const KinematicParcel<ParcelType>& ); /*---------------------------------------------------------------------------*\ Class KinematicParcel Declaration \*---------------------------------------------------------------------------*/ template<class ParcelType> class KinematicParcel : public ParcelType { // Private data //- Size in bytes of the fields static const std::size_t sizeofFields_; //- Number of particle tracking attempts before we assume that it stalls static label maxTrackAttempts; public: //- Class to hold kinematic particle constant properties class constantProperties { protected: // Protected data //- Constant properties dictionary const dictionary dict_; private: // Private data //- Parcel type id - used for post-processing to flag the type // of parcels issued by this cloud demandDrivenEntry<label> parcelTypeId_; //- Minimum density [kg/m3] demandDrivenEntry<scalar> rhoMin_; //- Particle density [kg/m3] (constant) demandDrivenEntry<scalar> rho0_; //- Minimum parcel mass [kg] demandDrivenEntry<scalar> minParcelMass_; //- myadd Particle initial radius [m] demandDrivenEntry<scalar> Ro0_; //- myadd Particle initial pressure inside [kg/m/s2] demandDrivenEntry<scalar> pgo0_; public: // Constructors //- Null constructor constantProperties(); //- Copy constructor constantProperties(const constantProperties& cp); //- Construct from dictionary constantProperties(const dictionary& parentDict); // Member functions //- Return const access to the constant properties dictionary inline const dictionary& dict() const; //- Return const access to the parcel type id inline label parcelTypeId() const; //- Return const access to the minimum density inline scalar rhoMin() const; //- Return const access to the particle density inline scalar rho0() const; //- Return const access to the minimum parcel mass inline scalar minParcelMass() const; //- myadd Return const access to the particle intitial radius inline scalar Ro0() const; //- myadd Return const access to the particle initial pressure inline scalar pgo0() const; }; template<class CloudType> class TrackingData : public ParcelType::template TrackingData<CloudType> { public: enum trackPart { tpVelocityHalfStep, tpLinearTrack, tpRotationalTrack }; private: // Private data // Interpolators for continuous phase fields //- Density interpolator autoPtr<interpolation<scalar> > rhoInterp_; //- Velocity interpolator autoPtr<interpolation<vector> > UInterp_; //- Dynamic viscosity interpolator autoPtr<interpolation<scalar> > muInterp_; //-myadd Pressure interpolator autoPtr<interpolation<scalar> > PInterp_; //- Local gravitational or other body-force acceleration const vector& g_; // label specifying which part of the integration // algorithm is taking place trackPart part_; public: // Constructors //- Construct from components inline TrackingData ( CloudType& cloud, trackPart part = tpLinearTrack ); // Member functions //- Return conat access to the interpolator for continuous // phase density field inline const interpolation<scalar>& rhoInterp() const; //- Return conat access to the interpolator for continuous // phase velocity field inline const interpolation<vector>& UInterp() const; //- Return conat access to the interpolator for continuous // phase dynamic viscosity field inline const interpolation<scalar>& muInterp() const; //- myadd Return const access to the interpolator for continuous phase pressure field inline const interpolation<scalar>& PInterp() const; // Return const access to the gravitational acceleration vector inline const vector& g() const; //- Return the part of the tracking operation taking place inline trackPart part() const; //- Return access to the part of the tracking operation taking place inline trackPart& part(); }; protected: // Protected data // Parcel properties //- Active flag - tracking inactive when active = false bool active_; //- Parcel type id label typeId_; //- Number of particles in Parcel scalar nParticle_; //- Diameter [m] scalar d_; //- Target diameter [m] scalar dTarget_; //- Velocity of Parcel [m/s] vector U_; //- Density [kg/m3] scalar rho_; //- Age [s] scalar age_; //- Time spent in turbulent eddy [s] scalar tTurb_; //- Turbulent velocity fluctuation [m/s] vector UTurb_; //- myadd change rate of bubble radius [m/s] scalar dRt_; //- myadd initial bubble radius [m] scalar Ro_; //- myadd initial bubble pressuer [Pa] scalar pgo_; // Cell-based quantities //- Density [kg/m3] scalar rhoc_; //- Velocity [m/s] vector Uc_; //- Viscosity [Pa.s] scalar muc_; //- myadd Pressure [Pa] scalar Pc_; // Protected Member Functions //- Calculate new particle velocity template<class TrackData> const vector calcVelocity ( TrackData& td, const scalar dt, // timestep const label cellI, // owner cell const scalar Re, // Reynolds number const scalar mu, // local carrier viscosity const scalar mass, // mass const vector& Su, // explicit particle momentum source vector& dUTrans, // momentum transfer to carrier scalar& Spu // linearised drag coefficient ) const; //- myadd Calculate new particle Diameter template<class TrackData> const vector calcDiameter ( TrackData& td, const scalar dt, const label cellI, const scalar mu, const scalar rhoc ) const; public: // Static data members //- Runtime type information TypeName("KinematicParcel"); //- String representation of properties AddToPropertyList ( ParcelType, " active" + " typeId" + " nParticle" + " d" + " dTarget " + " (Ux Uy Uz)" + " rho" + " age" //- myadd + " dRt" + " Ro" + " pgo" //- end + " tTurb" + " (UTurbx UTurby UTurbz)" ); // Constructors //- Construct from owner, position, and cloud owner // Other properties initialised as null inline KinematicParcel ( const polyMesh& mesh, const vector& position, const label cellI, const label tetFaceI, const label tetPtI ); //- Construct from components inline KinematicParcel ( const polyMesh& mesh, const vector& position, const label cellI, const label tetFaceI, const label tetPtI, const label typeId, const scalar nParticle0, const scalar d0, const scalar dTarget0, const vector& torque0, //- myadd const scalar dRt0, const constantProperties& constProps ); //- Construct from Istream KinematicParcel ( const polyMesh& mesh, Istream& is, bool readFields = true ); //- Construct as a copy KinematicParcel(const KinematicParcel& p); //- Construct as a copy KinematicParcel(const KinematicParcel& p, const polyMesh& mesh); //- Construct and return a (basic particle) clone virtual autoPtr<particle> clone() const { return autoPtr<particle>(new KinematicParcel(*this)); } //- Construct and return a (basic particle) clone virtual autoPtr<particle> clone(const polyMesh& mesh) const { return autoPtr<particle>(new KinematicParcel(*this, mesh)); } //- Factory class to read-construct particles used for // parallel transfer class iNew { const polyMesh& mesh_; public: iNew(const polyMesh& mesh) : mesh_(mesh) {} autoPtr<KinematicParcel<ParcelType> > operator()(Istream& is) const { return autoPtr<KinematicParcel<ParcelType> > ( new KinematicParcel<ParcelType>(mesh_, is, true) ); } }; // Member Functions // Access //- Return const access to active flag inline bool active() const; //- Return const access to type id inline label typeId() const; //- Return const access to number of particles inline scalar nParticle() const; //- Return const access to diameter inline scalar d() const; //- Return const access to target diameter inline scalar dTarget() const; //- Return const access to velocity inline const vector& U() const; //- Return const access to density inline scalar rho() const; //- Return const access to the age inline scalar age() const; //- myadd Return const access to the change rate of bubble radius inline scalar dRt() const; //- myadd Return const access to intial bubble radius inline scalar Ro() const; //- myadd Return const access to initial pressure inside inline scalar pgo() const; //- Return const access to time spent in turbulent eddy inline scalar tTurb() const; //- Return const access to turbulent velocity fluctuation inline const vector& UTurb() const; //- Return const access to carrier density [kg/m3] inline scalar rhoc() const; //- Return const access to carrier velocity [m/s] inline const vector& Uc() const; //- Return const access to carrier viscosity [Pa.s] inline scalar muc() const; //- myadd Return const carrier pressure inline scalar Pc() const; // Edit //- Return const access to active flag inline bool& active(); //- Return access to type id inline label& typeId(); //- Return access to number of particles inline scalar& nParticle(); //- Return access to diameter inline scalar& d(); //- Return access to target diameter inline scalar& dTarget(); //- Return access to velocity inline vector& U(); //- Return access to density inline scalar& rho(); //- Return access to the age inline scalar& age(); //- myadd Return access to the change rate of bubble radius inline scalar& dRt(); //- myadd Return access to intial bubble radius inline scalar& Ro(); //- myadd Return access to intial pressure inside inline scalar& pgo(); //- Return access to time spent in turbulent eddy inline scalar& tTurb(); //- Return access to turbulent velocity fluctuation inline vector& UTurb(); // Helper functions //- Return the index of the face used in the interpolation routine inline label faceInterpolation() const; //- Cell owner mass inline scalar massCell(const label cellI) const; //- Particle mass inline scalar mass() const; //- Particle moment of inertia around diameter axis inline scalar momentOfInertia() const; //- Particle volume inline scalar volume() const; //- Particle volume for a given diameter inline static scalar volume(const scalar d); //- Particle projected area inline scalar areaP() const; //- Projected area for given diameter inline static scalar areaP(const scalar d); //- Particle surface area inline scalar areaS() const; //- Surface area for given diameter inline static scalar areaS(const scalar d); //- Reynolds number inline scalar Re ( const vector& U, // particle velocity const scalar d, // particle diameter const scalar rhoc, // carrier density const scalar muc // carrier dynamic viscosity ) const; //- Weber number inline scalar We ( const vector& U, // particle velocity const scalar d, // particle diameter const scalar rhoc, // carrier density const scalar sigma // particle surface tension ) const; //- Eotvos number inline scalar Eo ( const vector& a, // acceleration const scalar d, // particle diameter const scalar sigma // particle surface tension ) const; // Main calculation loop //- Set cell values template<class TrackData> void setCellValues ( TrackData& td, const scalar dt, const label cellI ); //- Correct cell values using latest transfer information template<class TrackData> void cellValueSourceCorrection ( TrackData& td, const scalar dt, const label cellI ); //- Update parcel properties over the time interval template<class TrackData> void calc ( TrackData& td, const scalar dt, const label cellI ); // Tracking //- Move the parcel template<class TrackData> bool move(TrackData& td, const scalar trackTime); // Patch interactions //- Overridable function to handle the particle hitting a face // without trackData void hitFace(int& td); //- Overridable function to handle the particle hitting a face template<class TrackData> void hitFace(TrackData& td); //- Overridable function to handle the particle hitting a patch // Executed before other patch-hitting functions template<class TrackData> bool hitPatch ( const polyPatch& p, TrackData& td, const label patchI, const scalar trackFraction, const tetIndices& tetIs ); //- Overridable function to handle the particle hitting a // processorPatch template<class TrackData> void hitProcessorPatch ( const processorPolyPatch&, TrackData& td ); //- Overridable function to handle the particle hitting a wallPatch template<class TrackData> void hitWallPatch ( const wallPolyPatch&, TrackData& td, const tetIndices& ); //- Overridable function to handle the particle hitting a polyPatch template<class TrackData> void hitPatch ( const polyPatch&, TrackData& td ); //- Transform the physical properties of the particle // according to the given transformation tensor virtual void transformProperties(const tensor& T); //- Transform the physical properties of the particle // according to the given separation vector virtual void transformProperties(const vector& separation); //- The nearest distance to a wall that the particle can be // in the n direction virtual scalar wallImpactDistance(const vector& n) const; // I-O //- Read template<class CloudType> static void readFields(CloudType& c); //- Write template<class CloudType> static void writeFields(const CloudType& c); // Ostream Operator friend Ostream& operator<< <ParcelType> ( Ostream&, const KinematicParcel<ParcelType>& ); }; // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace Foam // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // #include "KinematicParcelI.H" #include "KinematicParcelTrackingDataI.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // #ifdef NoRepository #include "KinematicParcel.C" #endif // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // #endif // ************************************************************************* //以及KinematicParcel.C
/*---------------------------------------------------------------------------*\ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright (C) 2011-2014 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/>. \*---------------------------------------------------------------------------*/ #include "KinematicParcel.H" #include "forceSuSp.H" #include "IntegrationScheme.H" #include "meshTools.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // template<class ParcelType> Foam::label Foam::KinematicParcel<ParcelType>::maxTrackAttempts = 1; // * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * * // template<class ParcelType> template<class TrackData> void Foam::KinematicParcel<ParcelType>::setCellValues ( TrackData& td, const scalar dt, const label cellI ) { tetIndices tetIs = this->currentTetIndices(); rhoc_ = td.rhoInterp().interpolate(this->position(), tetIs); if (rhoc_ < td.cloud().constProps().rhoMin()) { if (debug) { WarningIn ( "void Foam::KinematicParcel<ParcelType>::setCellValues" "(" "TrackData&, " "const scalar, " "const label" ")" ) << "Limiting observed density in cell " << cellI << " to " << td.cloud().constProps().rhoMin() << nl << endl; } rhoc_ = td.cloud().constProps().rhoMin(); } Uc_ = td.UInterp().interpolate(this->position(), tetIs); muc_ = td.muInterp().interpolate(this->position(), tetIs); // Apply dispersion components to carrier phase velocity Uc_ = td.cloud().dispersion().update ( dt, cellI, U_, Uc_, UTurb_, tTurb_ ); Pc_ = td.PInterp().interpolate(this->position(), tetIs); //- myadd } template<class ParcelType> template<class TrackData> void Foam::KinematicParcel<ParcelType>::cellValueSourceCorrection ( TrackData& td, const scalar dt, const label cellI ) { Uc_ += td.cloud().UTrans()[cellI]/massCell(cellI); } template<class ParcelType> template<class TrackData> void Foam::KinematicParcel<ParcelType>::calc ( TrackData& td, const scalar dt, const label cellI ) { // Define local properties at beginning of time step // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ const scalar np0 = nParticle_; const scalar mass0 = mass(); // Reynolds number const scalar Re = this->Re(U_, d_, rhoc_, muc_); // Sources //~~~~~~~~ // Explicit momentum source for particle vector Su = vector::zero; // Linearised momentum source coefficient scalar Spu = 0.0; // Momentum transfer from the particle to the carrier phase vector dUTrans = vector::zero; //- myadd if(Ro_>0) { vector storeddRt = calcDiameter(td, dt, cellI, muc_, rhoc_); this->dRt_=storeddRt[0]; this->d_=storeddRt[1]; } // Motion // ~~~~~~ // Calculate new particle velocity this->U_ = calcVelocity(td, dt, cellI, Re, muc_, mass0, Su, dUTrans, Spu); // Accumulate carrier phase source terms // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if (td.cloud().solution().coupled()) { // Update momentum transfer td.cloud().UTrans()[cellI] += np0*dUTrans; // Update momentum transfer coefficient td.cloud().UCoeff()[cellI] += np0*Spu; } } template<class ParcelType> template<class TrackData> const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity ( TrackData& td, const scalar dt, const label cellI, const scalar Re, const scalar mu, const scalar mass, const vector& Su, vector& dUTrans, scalar& Spu ) const { typedef typename TrackData::cloudType cloudType; typedef typename cloudType::parcelType parcelType; typedef typename cloudType::forceType forceType; const forceType& forces = td.cloud().forces(); // Momentum source due to particle forces const parcelType& p = static_cast<const parcelType&>(*this); const forceSuSp Fcp = forces.calcCoupled(p, dt, mass, Re, mu); const forceSuSp Fncp = forces.calcNonCoupled(p, dt, mass, Re, mu); const forceSuSp Feff = Fcp + Fncp; const scalar massEff = forces.massEff(p, mass); // New particle velocity //~~~~~~~~~~~~~~~~~~~~~~ // Update velocity - treat as 3-D const vector abp = (Feff.Sp()*Uc_ + (Feff.Su() + Su))/massEff; const scalar bp = Feff.Sp()/massEff; Spu = dt*Feff.Sp(); IntegrationScheme<vector>::integrationResult Ures = td.cloud().UIntegrator().integrate(U_, dt, abp, bp); vector Unew = Ures.value(); // note: Feff.Sp() and Fc.Sp() must be the same dUTrans += dt*(Feff.Sp()*(Ures.average() - Uc_) - Fcp.Su()); // Apply correction to velocity and dUTrans for reduced-D cases const polyMesh& mesh = td.cloud().pMesh(); meshTools::constrainDirection(mesh, mesh.solutionD(), Unew); meshTools::constrainDirection(mesh, mesh.solutionD(), dUTrans); return Unew; } //- myadd template<class ParcelType> template<class TrackData> const Foam::vector Foam::KinematicParcel<ParcelType>::calcDiameter ( TrackData& td, const scalar dt, const label cellI, const scalar mu, const scalar rhoc ) const { // Saturated vapor pressure const scalar pv=0; //surface tension const scalar siggma=0.0712; scalar R=d_/2; //The fourth-order Rungekuta method solves second-order differential equations const scalar k1=dRt_; const scalar h1=-3*pow(k1,2)/2/R-(4*mu*k1+2*siggma)/rhoc/pow(R,2)+(pv-Pc_+pgo_*pow(Ro_,4.2)*pow(R,-4.2))/rhoc/R; const scalar k2=k1+0.5*dt*h1; const scalar h2=-3*pow(k2,2)/2/(R+0.5*dt*k1)-(4*mu*k2+2*siggma)/rhoc/pow(R+0.5*dt*k1,2)+(pv-Pc_+pgo_*pow(Ro_,4.2)*pow(R+0.5*dt*k1,-4.2))/rhoc/(R+0.5*dt*k1); const scalar k3=k1+0.5*dt*h2; const scalar h3=-3*pow(k3,2)/2/(R+0.5*dt*k2)-(4*mu*k3+2*siggma)/rhoc/pow(R+0.5*dt*k2,2)+(pv-Pc_+pgo_*pow(Ro_,4.2)*pow(R+0.5*dt*k2,-4.2))/rhoc/(R+0.5*dt*k2); const scalar k4=k1+dt*h3; const scalar h4=-3*pow(k4,2)/2/(R+0.5*dt*k3)-(4*mu*k4+2*siggma)/rhoc/pow(R+0.5*dt*k3,2)+(pv-Pc_+pgo_*pow(Ro_,4.2)*pow(R+0.5*dt*k3,-4.2))/rhoc/(R+0.5*dt*k3); vector Rtvector(0,0,0); Rtvector[0]=k1+dt/6*(h1+2*h2+2*h3+h4); Rtvector[1]=2*(R+dt/6*(k1+2*k2+2*k3+k4)); return Rtvector; } // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // template<class ParcelType> Foam::KinematicParcel<ParcelType>::KinematicParcel ( const KinematicParcel<ParcelType>& p ) : ParcelType(p), active_(p.active_), typeId_(p.typeId_), nParticle_(p.nParticle_), d_(p.d_), dTarget_(p.dTarget_), U_(p.U_), rho_(p.rho_), age_(p.age_), tTurb_(p.tTurb_), UTurb_(p.UTurb_), rhoc_(p.rhoc_), Uc_(p.Uc_), muc_(p.muc_), //- myadd dRt_(p.dRt_), Ro_(p.Ro_), pgo_(p.pgo_), Pc_(p.Pc_) {} template<class ParcelType> Foam::KinematicParcel<ParcelType>::KinematicParcel ( const KinematicParcel<ParcelType>& p, const polyMesh& mesh ) : ParcelType(p, mesh), active_(p.active_), typeId_(p.typeId_), nParticle_(p.nParticle_), d_(p.d_), dTarget_(p.dTarget_), U_(p.U_), rho_(p.rho_), age_(p.age_), tTurb_(p.tTurb_), UTurb_(p.UTurb_), rhoc_(p.rhoc_), Uc_(p.Uc_), muc_(p.muc_), //- myadd dRt_(p.dRt_), Ro_(p.Ro_), pgo_(p.pgo_), Pc_(p.Pc_) {}我明明进行了定义，但是仍然报错，有没有相关操作过的大佬进行以下，这个方程植入我植入了两个月了还没成功
-
@李东岳 我忽略了第一个方程最后那个u和ub
R — 气泡半径
··R — 气泡半径对时间二阶导
·R — 气泡半径对时间一阶导
pv — 气泡内饱和蒸汽压（我这里设置的0，因为这个方程考虑了空化现象，我直接取的0）
pg — 气泡内部气压（根据第二个方程在计算）
p∞ — 气泡受到的液体压强（我这里使用PInterp插值，在DPMFoam求解器也自定义了个P=p*rhoc）
σ — 表面张力
pgo — 初始气泡内部压强
Ro — 初始气泡直径值得一提的是，修改完求解器之后，这个pgo和Ro我看我现在阅读的论文代码上是直接在KinematicCloudProperties的颗粒参数上进行了赋值。
在我上面那个KinematicParcel.C是主要算法，使用龙格库朗四阶对于这个二阶的方程进行的求解（他直接把方程2-49代入了2-48计算）。 -
\begin{equation}
r\frac{\p^2 r}{\p t^2}+1.5\frac{\p r}{\p t}=
\frac{1}{\rho}\left(p_{g0} \left(r_0/r\right)^{3\gamma} -p\rho_\rc-\frac{2\sigma}{r} - \frac{4\mu}{r}\frac{\p r}{\p t}\right)+\frac{(\bfU_\rc-\bfU)^2}{4}
\end{equation}我把这个方程简化成这个行不行
\begin{equation}
r\frac{\p^2 r}{\p t^2}+1.5\frac{\p r}{\p t}=A+\frac{(\bfU_\rc-\bfU)^2}{4}
\end{equation}$A=0$
-
-
-
我在几年前做了个一拉格朗日粒子直径变化的代码。我现在已经记不起来太多了。能想起来的是：
- 粒子直径跟流场的湍流动能耗散率有关系
- 植入的是一个跟PBM有关的粒径变化方程，考虑了coalescense和破碎
- 代码可以编译并且测试过可以体现粒径变化
可以把下面的
MomentumCloud.C：替换到OpenFOAM原来的代码然后编译。搞明白之后看看思路然后适配到你们的算法上。MomentumCloud.C：
/*---------------------------------------------------------------------------*\ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org \\ / A nd | Copyright (C) 2011-2021 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/>. \*---------------------------------------------------------------------------*/ #include "MomentumCloud.H" #include "integrationScheme.H" #include "interpolation.H" #include "subCycleTime.H" #include "InjectionModelList.H" #include "DispersionModel.H" #include "PatchInteractionModel.H" #include "StochasticCollisionModel.H" #include "SurfaceFilmModel.H" // * * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * // template<class CloudType> void Foam::MomentumCloud<CloudType>::setModels() { dispersionModel_.reset ( DispersionModel<MomentumCloud<CloudType>>::New ( subModelProperties_, *this ).ptr() ); patchInteractionModel_.reset ( PatchInteractionModel<MomentumCloud<CloudType>>::New ( subModelProperties_, *this ).ptr() ); stochasticCollisionModel_.reset ( StochasticCollisionModel<MomentumCloud<CloudType>>::New ( subModelProperties_, *this ).ptr() ); surfaceFilmModel_.reset ( SurfaceFilmModel<MomentumCloud<CloudType>>::New ( subModelProperties_, *this ).ptr() ); UIntegrator_.reset ( integrationScheme::New ( "U", solution_.integrationSchemes() ).ptr() ); } template<class CloudType> template<class TrackCloudType> void Foam::MomentumCloud<CloudType>::solve ( TrackCloudType& cloud, typename parcelType::trackingData& td ) { if (solution_.steadyState()) { cloud.storeState(); cloud.preEvolve(); evolveCloud(cloud, td); if (solution_.coupled()) { cloud.relaxSources(cloud.cloudCopy()); } } else { ///////////////////////////////// Info<< "Dyfluid: Evolve function in MomentumCloud.C" << nl; List<DynamicList<label>> PartLabelPre(U_.size()); List<DynamicList<label>> PartLabelPost(U_.size()); DynamicList<label> cellWithPartPre; DynamicList<label> cellWithPartPost; pNumber_.primitiveFieldRef() = 0.0; forAllIter(typename MomentumCloud<CloudType>, *this, iter) { parcelType& p = iter(); pNumber_[p.cell()] += p.nParticle(); PartLabelPre[p.cell()].append(p.origId()); } forAll(U_, cell) { if (PartLabelPre[cell].size() != 0) { cellWithPartPre.append(cell); } } //Info<< cellWithPartPre << nl; ///////////////////////////////// cloud.preEvolve(); evolveCloud(cloud, td); if (solution_.coupled()) { cloud.scaleSources(); } ///////////////////////////////// forAllIter(typename MomentumCloud<CloudType>, *this, iter) { parcelType& p = iter(); PartLabelPost[p.cell()].append(p.origId()); } const volScalarField& epsi = U_.mesh().lookupObject<volScalarField>("epsilon.water"); scalar C1 = 0.00481; scalar C2 = 0.08; scalar sigma = 0.072; scalar rhoc = 998.0; scalar D1 = 0.88; scalar D2 = 9e6; scalar muc = 1e-3; const scalar dt = this->mesh().time().deltaTValue(); forAll(U_, cell) { // if there are particles in some cell if (PartLabelPost[cell].size() != 0) { cellWithPartPost.append(cell); scalar allVolume = 0.0; scalar allnParticle = 0.0; scalar allnParticleNew = 0.0; scalar diamAve = 0.0; //Info<< "Cell [" << cell << "] has " << PartLabelPost[cell].size() // << " particles" << nl; // loop all over particles forAllIter(typename MomentumCloud<CloudType>, *this, iter) { parcelType& p = iter(); // loop particles label in this cell forAll(PartLabelPost[cell], i) { // if this particle belongs to this cell if (p.origId() == PartLabelPost[cell](i)) { allVolume += p.nParticle()*M_PI/6.0*pow3(p.d()); allnParticle += p.nParticle(); diamAve = pow(allVolume/allnParticle*6.0/M_PI, 1.0/3.0); //Info<< " Labels are " << PartLabelPost[cell](i) // << ", its nParticle is " << p.nParticle() // << ", volume is " << allVolume << nl; } } } const scalar d = diamAve; // cell manipulation, calculate average diameter if (d > SMALL) { const scalar epsilon = epsi[cell]; scalar gN = allnParticle*C1*pow(epsilon, 1.0/3.0)/pow(d, 2.0/3.0) *exp(-C2*sigma/(rhoc*pow(epsilon, 2.0/3.0)*pow(d, 5.0/3.0))); scalar aSqrN = sqr(allnParticle)*D1*pow(epsilon, 1.0/3.0)*4.0*sqrt(2.0)*pow(d, 7.0/6.0) *exp(-D2*muc*rhoc*epsilon/sqr(sigma)*pow4(d)/16.0); allnParticleNew = allnParticle + (gN - aSqrN)*dt; diamAve = diamAve*pow(allnParticle/allnParticleNew, 1.0/3.0); //Info<< "gN = " << gN << nl; //Info<< "aSqrN = " << aSqrN << nl; //Info<< "Particle increased by " << allnParticleNew - allnParticle << nl; //Info<< "diamAve is " << diamAve << nl; } scalar newVolume = 0.0; forAllIter(typename MomentumCloud<CloudType>, *this, iter) { parcelType& p = iter(); forAll(PartLabelPost[cell], i) { if (p.origId() == PartLabelPost[cell](i)) { p.nParticle() = allnParticleNew/PartLabelPost[cell].size(); p.d() = diamAve; newVolume += p.nParticle()*M_PI/6.0*pow3(p.d()); //Info<< " after brecoa, volume is " << newVolume // << "p.nParticle is " << p.nParticle() << nl; } } } } } ///////////////////////////////// } cloud.info(); cloud.postEvolve(); if (solution_.steadyState()) { cloud.restoreState(); } } template<class CloudType> void Foam::MomentumCloud<CloudType>::buildCellOccupancy() { if (cellOccupancyPtr_.empty()) { cellOccupancyPtr_.reset ( new List<DynamicList<parcelType*>>(this->mesh().nCells()) ); } else if (cellOccupancyPtr_().size() != this->mesh().nCells()) { // If the size of the mesh has changed, reset the // cellOccupancy size cellOccupancyPtr_().setSize(this->mesh().nCells()); } List<DynamicList<parcelType*>>& cellOccupancy = cellOccupancyPtr_(); forAll(cellOccupancy, cO) { cellOccupancy[cO].clear(); } forAllIter(typename MomentumCloud<CloudType>, *this, iter) { cellOccupancy[iter().cell()].append(&iter()); } } template<class CloudType> void Foam::MomentumCloud<CloudType>::updateCellOccupancy() { // Only build the cellOccupancy if the pointer is set, i.e. it has // been requested before. if (cellOccupancyPtr_.valid()) { buildCellOccupancy(); } } template<class CloudType> template<class TrackCloudType> void Foam::MomentumCloud<CloudType>::evolveCloud ( TrackCloudType& cloud, typename parcelType::trackingData& td ) { if (solution_.coupled()) { cloud.resetSourceTerms(); } if (solution_.transient()) { label preInjectionSize = this->size(); this->surfaceFilm().inject(cloud); // Update the cellOccupancy if the size of the cloud has changed // during the injection. if (preInjectionSize != this->size()) { updateCellOccupancy(); preInjectionSize = this->size(); } injectors_.inject(cloud, td); // Assume that motion will update the cellOccupancy as necessary // before it is required. cloud.motion(cloud, td); stochasticCollision().update(td, solution_.trackTime()); } else { // this->surfaceFilm().injectSteadyState(cloud); injectors_.injectSteadyState(cloud, td, solution_.trackTime()); CloudType::move(cloud, td, solution_.trackTime()); } } template<class CloudType> void Foam::MomentumCloud<CloudType>::postEvolve() { Info<< endl; if (debug) { this->writePositions(); } this->dispersion().cacheFields(false); forces_.cacheFields(false); functions_.postEvolve(); solution_.nextIter(); if (this->db().time().writeTime()) { outputProperties_.writeObject ( IOstream::ASCII, IOstream::currentVersion, this->db().time().writeCompression(), true ); } } template<class CloudType> void Foam::MomentumCloud<CloudType>::cloudReset(MomentumCloud<CloudType>& c) { CloudType::cloudReset(c); rndGen_ = c.rndGen_; forces_.transfer(c.forces_); functions_.transfer(c.functions_); injectors_.transfer(c.injectors_); dispersionModel_.reset(c.dispersionModel_.ptr()); patchInteractionModel_.reset(c.patchInteractionModel_.ptr()); stochasticCollisionModel_.reset(c.stochasticCollisionModel_.ptr()); surfaceFilmModel_.reset(c.surfaceFilmModel_.ptr()); UIntegrator_.reset(c.UIntegrator_.ptr()); } // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // template<class CloudType> Foam::MomentumCloud<CloudType>::MomentumCloud ( const word& cloudName, const volScalarField& rho, const volVectorField& U, const volScalarField& mu, const dimensionedVector& g, const bool readFields ) : CloudType(cloudName, rho, U, mu, g, false), cloudCopyPtr_(nullptr), particleProperties_ ( IOobject ( cloudName + "Properties", this->mesh().time().constant(), this->mesh(), IOobject::MUST_READ_IF_MODIFIED, IOobject::NO_WRITE ) ), outputProperties_ ( IOobject ( cloudName + "OutputProperties", this->mesh().time().timeName(), "uniform"/cloud::prefix/cloudName, this->mesh(), IOobject::READ_IF_PRESENT, IOobject::NO_WRITE ) ), solution_(this->mesh(), particleProperties_.subDict("solution")), constProps_(particleProperties_), subModelProperties_ ( particleProperties_.subOrEmptyDict("subModels", true) ), rndGen_(0), cellOccupancyPtr_(), cellLengthScale_(mag(cbrt(this->mesh().V()))), pNumber_ ( IOobject ( "pNumbers", this->mesh().time().timeName(), this->mesh(), IOobject::NO_READ, IOobject::AUTO_WRITE ), this->mesh(), dimensionedScalar(dimless, 0.0) ), rho_(rho), U_(U), mu_(mu), g_(g), pAmbient_(0.0), forces_ ( *this, this->mesh(), subModelProperties_.subOrEmptyDict ( "particleForces", true ), true ), functions_ ( *this, particleProperties_.subOrEmptyDict("cloudFunctions"), true ), injectors_ ( subModelProperties_.subOrEmptyDict("injectionModels"), *this ), dispersionModel_(nullptr), patchInteractionModel_(nullptr), stochasticCollisionModel_(nullptr), surfaceFilmModel_(nullptr), UIntegrator_(nullptr), UTrans_ ( new volVectorField::Internal ( IOobject ( this->name() + ":UTrans", this->db().time().timeName(), this->db(), IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE ), this->mesh(), dimensionedVector(dimMass*dimVelocity, Zero) ) ), UCoeff_ ( new volScalarField::Internal ( IOobject ( this->name() + ":UCoeff", this->db().time().timeName(), this->db(), IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE ), this->mesh(), dimensionedScalar( dimMass, 0) ) ) { setModels(); if (readFields) { parcelType::readFields(*this); this->deleteLostParticles(); } if (solution_.resetSourcesOnStartup()) { resetSourceTerms(); } } template<class CloudType> Foam::MomentumCloud<CloudType>::MomentumCloud ( const word& cloudName, const volScalarField& rho, const volVectorField& U, const dimensionedVector& g, const fluidThermo& carrierThermo, const bool readFields ) : MomentumCloud(cloudName, rho, U, carrierThermo.mu(), g, readFields) {} template<class CloudType> Foam::MomentumCloud<CloudType>::MomentumCloud ( MomentumCloud<CloudType>& c, const word& name ) : CloudType(c, name), cloudCopyPtr_(nullptr), particleProperties_(c.particleProperties_), outputProperties_(c.outputProperties_), solution_(c.solution_), constProps_(c.constProps_), subModelProperties_(c.subModelProperties_), rndGen_(c.rndGen_), cellOccupancyPtr_(nullptr), cellLengthScale_(c.cellLengthScale_), pNumber_(c.pNumber_), rho_(c.rho_), U_(c.U_), mu_(c.mu_), g_(c.g_), pAmbient_(c.pAmbient_), forces_(c.forces_), functions_(c.functions_), injectors_(c.injectors_), dispersionModel_(c.dispersionModel_->clone()), patchInteractionModel_(c.patchInteractionModel_->clone()), stochasticCollisionModel_(c.stochasticCollisionModel_->clone()), surfaceFilmModel_(c.surfaceFilmModel_->clone()), UIntegrator_(c.UIntegrator_->clone()), UTrans_ ( new volVectorField::Internal ( IOobject ( this->name() + ":UTrans", this->db().time().timeName(), this->db(), IOobject::NO_READ, IOobject::NO_WRITE, false ), c.UTrans_() ) ), UCoeff_ ( new volScalarField::Internal ( IOobject ( name + ":UCoeff", this->db().time().timeName(), this->db(), IOobject::NO_READ, IOobject::NO_WRITE, false ), c.UCoeff_() ) ) {} template<class CloudType> Foam::MomentumCloud<CloudType>::MomentumCloud ( const fvMesh& mesh, const word& name, const MomentumCloud<CloudType>& c ) : CloudType(mesh, name, c), cloudCopyPtr_(nullptr), particleProperties_ ( IOobject ( name + "Properties", mesh.time().constant(), mesh, IOobject::NO_READ, IOobject::NO_WRITE, false ) ), outputProperties_ ( IOobject ( name + "OutputProperties", this->mesh().time().timeName(), "uniform"/cloud::prefix/name, this->mesh(), IOobject::NO_READ, IOobject::NO_WRITE, false ) ), solution_(mesh), constProps_(), subModelProperties_(dictionary::null), rndGen_(0), cellOccupancyPtr_(nullptr), cellLengthScale_(c.cellLengthScale_), pNumber_(c.pNumber_), rho_(c.rho_), U_(c.U_), mu_(c.mu_), g_(c.g_), pAmbient_(c.pAmbient_), forces_(*this, mesh), functions_(*this), injectors_(*this), dispersionModel_(nullptr), patchInteractionModel_(nullptr), stochasticCollisionModel_(nullptr), surfaceFilmModel_(nullptr), UIntegrator_(nullptr), UTrans_(nullptr), UCoeff_(nullptr) {} // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // template<class CloudType> Foam::MomentumCloud<CloudType>::~MomentumCloud() {} // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // template<class CloudType> void Foam::MomentumCloud<CloudType>::setParcelThermoProperties ( parcelType& parcel, const scalar lagrangianDt ) { parcel.rho() = constProps_.rho0(); } template<class CloudType> void Foam::MomentumCloud<CloudType>::checkParcelProperties ( parcelType& parcel, const scalar lagrangianDt, const bool fullyDescribed ) { const scalar carrierDt = this->mesh().time().deltaTValue(); parcel.stepFraction() = (carrierDt - lagrangianDt)/carrierDt; if (parcel.typeId() == -1) { parcel.typeId() = constProps_.parcelTypeId(); } } template<class CloudType> void Foam::MomentumCloud<CloudType>::storeState() { cloudCopyPtr_.reset ( static_cast<MomentumCloud<CloudType>*> ( clone(this->name() + "Copy").ptr() ) ); } template<class CloudType> void Foam::MomentumCloud<CloudType>::restoreState() { cloudReset(cloudCopyPtr_()); cloudCopyPtr_.clear(); } template<class CloudType> void Foam::MomentumCloud<CloudType>::resetSourceTerms() { UTransRef().field() = Zero; UCoeffRef().field() = 0.0; } template<class CloudType> template<class Type> void Foam::MomentumCloud<CloudType>::relax ( DimensionedField<Type, volMesh>& field, const DimensionedField<Type, volMesh>& field0, const word& name ) const { const scalar coeff = solution_.relaxCoeff(name); field = field0 + coeff*(field - field0); } template<class CloudType> template<class Type> void Foam::MomentumCloud<CloudType>::scale ( DimensionedField<Type, volMesh>& field, const word& name ) const { const scalar coeff = solution_.relaxCoeff(name); field *= coeff; } template<class CloudType> void Foam::MomentumCloud<CloudType>::relaxSources ( const MomentumCloud<CloudType>& cloudOldTime ) { this->relax(UTrans_(), cloudOldTime.UTrans_(), "U"); this->relax(UCoeff_(), cloudOldTime.UCoeff_(), "U"); } template<class CloudType> void Foam::MomentumCloud<CloudType>::scaleSources() { this->scale(UTrans_(), "U"); this->scale(UCoeff_(), "U"); } template<class CloudType> void Foam::MomentumCloud<CloudType>::preEvolve() { // force calculation of mesh dimensions - needed for parallel runs // with topology change due to lazy evaluation of valid mesh dimensions label nGeometricD = this->mesh().nGeometricD(); Info<< "\nSolving " << nGeometricD << "-D cloud " << this->name() << endl; this->dispersion().cacheFields(true); forces_.cacheFields(true); updateCellOccupancy(); pAmbient_ = constProps_.dict().template lookupOrDefault<scalar>("pAmbient", pAmbient_); functions_.preEvolve(); } template<class CloudType> void Foam::MomentumCloud<CloudType>::evolve() { if (solution_.canEvolve()) { typename parcelType::trackingData td(*this); solve(*this, td); } } template<class CloudType> template<class TrackCloudType> void Foam::MomentumCloud<CloudType>::motion ( TrackCloudType& cloud, typename parcelType::trackingData& td ) { CloudType::move(cloud, td, solution_.trackTime()); updateCellOccupancy(); } template<class CloudType> void Foam::MomentumCloud<CloudType>::patchData ( const parcelType& p, const polyPatch& pp, vector& nw, vector& Up ) const { p.patchData(nw, Up); Up /= p.mesh().time().deltaTValue(); // If this is a wall patch, then there may be a non-zero tangential velocity // component; the lid velocity in a lid-driven cavity case, for example. We // want the particle to interact with this velocity, so we look it up in the // velocity field and use it to set the wall-tangential component. if (!this->mesh().moving() && isA<wallPolyPatch>(pp)) { const label patchi = pp.index(); const label patchFacei = pp.whichFace(p.face()); // We only want to use the boundary condition value only if it is set // by the boundary condition. If the boundary values are extrapolated // (e.g., slip conditions) then they represent the motion of the fluid // just inside the domain rather than that of the wall itself. if (U_.boundaryField()[patchi].fixesValue()) { const vector Uw1 = U_.boundaryField()[patchi][patchFacei]; const vector& Uw0 = U_.oldTime().boundaryField()[patchi][patchFacei]; const scalar f = p.currentTimeFraction(); const vector Uw = Uw0 + f*(Uw1 - Uw0); const tensor nnw = nw*nw; Up = (nnw & Up) + Uw - (nnw & Uw); } } } template<class CloudType> void Foam::MomentumCloud<CloudType>::updateMesh() { updateCellOccupancy(); injectors_.updateMesh(); cellLengthScale_ = mag(cbrt(this->mesh().V())); } template<class CloudType> void Foam::MomentumCloud<CloudType>::autoMap(const mapPolyMesh& mapper) { Cloud<parcelType>::autoMap(mapper); updateMesh(); } template<class CloudType> void Foam::MomentumCloud<CloudType>::info() { vector linearMomentum = linearMomentumOfSystem(); reduce(linearMomentum, sumOp<vector>()); scalar linearKineticEnergy = linearKineticEnergyOfSystem(); reduce(linearKineticEnergy, sumOp<scalar>()); Info<< "Cloud: " << this->name() << nl << " Current number of parcels = " << returnReduce(this->size(), sumOp<label>()) << nl << " Current mass in system = " << returnReduce(massInSystem(), sumOp<scalar>()) << nl << " Linear momentum = " << linearMomentum << nl << " |Linear momentum| = " << mag(linearMomentum) << nl << " Linear kinetic energy = " << linearKineticEnergy << nl; injectors_.info(Info); this->surfaceFilm().info(Info); this->patchInteraction().info(Info); } // ************************************************************************* //
