拉格朗日求解器代码求助

[youhaoyu]

各位大佬好，我想要在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_)
{}
 
 
 
 

我明明进行了定义，但是仍然报错，有没有相关操作过的大佬进行以下，这个方程植入我植入了两个月了还没成功

[李东岳]

把你这两个方程里面自定义参数的值都写一下我看看。

[youhaoyu]

@李东岳 我忽略了第一个方程最后那个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{array}{}\text{(1)}& r\frac{{\mathrm{\partial }}^{2}r}{\mathrm{\partial }{t}^2}+1.5\frac{\mathrm{\partial }r}{\mathrm{\partial }t}=\frac{1}{\rho }(p_{g0}{\left(r_0/r\right)}^{3\gamma }-p{\rho }_{\mathrm{c}}-\frac{2\sigma }{r}-\frac{4\mu }{r}\frac{\mathrm{\partial }r}{\mathrm{\partial }t})+\frac{({\mathbf{U}}_{\mathrm{c}}-\mathbf{U}{)}^2}{4}\end{array}$$

我把这个方程简化成这个行不行

$$\begin{array}{}\text{(2)}& r\frac{{\mathrm{\partial }}^{2}r}{\mathrm{\partial }{t}^2}+1.5\frac{\mathrm{\partial }r}{\mathrm{\partial }t}=A+\frac{({\mathbf{U}}_{\mathrm{c}}-\mathbf{U}{)}^2}{4}\end{array}$$

$A=0$

[youhaoyu]

@李东岳 这个把右侧括号部分考虑成0的话可能不大妥吧，因为我是个动态过程，气泡的入口和出口的压力差相差了两个数量级，还是挺大的。要看这个气泡慢慢变大过程。
如果你把A考虑成0，方程无法体现气泡内外压力差和表面张力，和我的模型可能也有一些问题。
方程的解法其实也不是很大问题，主要是这些变量如何去申明定义太恶心了。东岳老师您这种就简化的模型不大适合我目前的课题。 最开始我也想过，直接使用最简单的理想方程PV=NRT这种去修改，但是我觉得这种方程使用起来不是特别好

[李东岳]

我只是提供一个框架，而不是实际把r计算出来。后续你要在框架上自己进行修改，比如把A改成真实的。

[youhaoyu]

@李东岳 我的模型大概是1米高的管道，下方10000大气压，上方600Pa（液体是钢液，密度7000），底部通入气泡这种模型，对于气泡和液体之间压强差这个就是我们主要考虑的因素。这种来说A感觉不能直接考虑为0，所以我一直在看如何直接植入这个完整的方程

[杨英狄]

粒子性质的定义不仅在KinematicParcel.C和KinematicParcel.H，在其他文件中也涉及相关代码。你可以正则搜索一下，看看那些文件涉及粒子性质的声明和定义。

[youhaoyu]

@杨英狄 谢谢，这几天我我看了看，是我有些地方没有定义完成，我看您之前也在往DPMFoam添加RP方程，添加成功了吗？我最近添加编译成功了，库和求解器都成功了，但是计算出来结果颗粒直径没有变化。（您好久没登陆了，我还在尝试联系您，没成功）

[杨英狄]

@youhaoyu 如果计算过程中发现粒径没有变化，可能是添加的RP方程计算程序没有被调用。参考原生的lagrange库，粒子速度的更新函数calcVelocity（），是在函数calc()里面被调用的。你可以看看你的RP方程函数是否在calc()里面调用了。

[youhaoyu]

@杨英狄 是的，我就是这么编写的，我写了个calcDiameter，但是就是没发现哪儿有问题，现在我在一点点改看看哪儿有问题。
但是现在麻烦来了，修改了代码之后，编译能成功，但是使用求解器时候会报错，什么浮点数异常什么的，我初始条件都没改，哎麻烦死了。能不能给个联系方式，想咨询下您问题，因为刚好看到您也修改过

[mingyang]

@youhaoyu 你好，能学习一下您的植入方法吗，我最近也被困在这了，期待你的回复，感谢

[youhaoyu]

@mingyang 您好可以加我qq 2269609762一起探讨，我也没做成功

[李东岳]

我在几年前做了个一拉格朗日粒子直径变化的代码。我现在已经记不起来太多了。能想起来的是：

1. 粒子直径跟流场的湍流动能耗散率有关系
2. 植入的是一个跟PBM有关的粒径变化方程，考虑了coalescense和破碎
3. 代码可以编译并且测试过可以体现粒径变化

可以把下面的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);
}
 
 
// ************************************************************************* //