Euler-Lagrange的一些解析，sprayFoam

[youhaoyu]

@星星星星晴 您好，大佬，我最近潜心修心了一个月，了解了拉格朗日源文件里面一些大概的东西，我目前打算根据一个方程更新我的拉格朗日颗粒的直径，所以我根据一篇名叫“水质对梢涡空化初生的影响研究”的博士论文在修改拉格朗日的源文件库，最近遇到了很多问题，不知道您能否给我指引一点方向来学习。
这个代码原理就是在KinematicParcel.C中按照他原本的calcVelocity，自己定义一个气泡长大方程取名为calcDimater，但是这篇博士论文的代码似乎有点不全，我一直在修改。目前改的文件有以下这6个基础文件

KinematicCloud.C
KinematicCloud.H
KinematicCloudI.H
KinematicParcel.C
KinematicParcel.H
KinematicParcelI.H

我能编译成功源文件库，但是在编译求解器时候会有报错，于是我又简单更改了这三个文件夹

CollidingCloud.C
CollidingCloud.H
CollidingCloudI.H

目前我的库能够编译成功，就是存在warning。但是当我编译我的新的newDPMFoam时候总是有很多奇奇怪怪报错，目前我的报错就很奇怪。总而言之就是说我有一些参数我定义错误

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

我这个原理上我感觉我的代码基本上都能说得通，但是报错经常就提醒我缺少一些东西，我现在就有一点迷茫，不知道怎么前进了，请大佬指点下迷津。我这个方法是否有可行性 附上一部分我的代码（//- myadd部分即为我添加部分）

KinematicParcel.C（这是方程的核心部分在里面，其他文件我只是按照原有的模板照猫画虎定义了一些参数，例如压力场的插值）

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_, Pc_);
    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 scalar Pc
) const
{
    // Saturated vapor pressure
    const scalar pv=0;
    //surface tension
    const scalar sigma=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*sigma)/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*sigma)/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*sigma)/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*sigma)/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;
}
 

KinematicParcel.H

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;
 
            // 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();
            //- myadd Return const access to the interpolator for continuous phase pressure field
            inline const interpolation<scalar>& PInterp() const;
    };
 
 
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 scalar Pc
        ) 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"
          + " tTurb"
          + " (UTurbx UTurby UTurbz)"
          //- myadd
          + " dRt"
          + " Ro"
          + " pgo"
        );

[星星星星晴]

@youhaoyu

hi，我已经好久没碰OF了。

但是我看起来 你可以先将calcDiameter 不要用function的方式植入，先用最简单的方式直接添加到 calc中。你把Ro()也转换成类似Re的那种简单的定义方式。

scalar Ro = xxxx

因为OF是一层一层的按照template编译的，所以你现在这个处理办法要一层一层往上走，改了底层的，你先按照简单植入办法植入，然后当你发现你不得不突破OF的限制了以后，再考虑其他的办法。

[dxl]

您好想请教一下，关于颗粒的碰撞OF中是怎么检测的呢，感觉不是通过两个循环对每两个颗粒做一次检测，您对OF中采用的方法了解吗，具体是在哪部分中体现的呢

[星星星星晴]

@dxl 首先要看你基于的哪个cloud做的，并不是所有的cloud都真实的包含了颗粒碰撞。要查一下code确认你的程序是不是真正的引用了碰撞模型。

其次，OF中的碰撞模型是O'Rourke的碰撞模型，是计算碰撞概率的，不是deteministic的（不知道中文怎么说）。即不是追踪所有颗粒，计算两个颗粒的路径，判断两个颗粒是否会发生碰撞的方法。

基于OF10的话，你可以在下面找到具体的源程序。

OpenFOAM-10/src/lagrangian/parcel/submodels/Spray/StochasticCollision
/ORourkeCollision/

[ann]

@星星星星晴 您好 想要请教一下flowRateProfile table里面的数值是怎么计算出来的，我的模型中将总喷射质量改为了3.6e-6，时长改为了1.54ms，之前想要简单地把时长均匀扩充到1.54ms，但是计算结果相差较大，想要问一下文档中提到的基于TimeFunction1是需要加一个boundaryField文件，在其中计算，还是通过函数拟合？

[星星星星晴]

@ann Hi 我好久没有碰OF了，没有办法给你精准的结果。
但是我建议你查一下源程序， 看看哪一项是控制质量流量，np是怎么计算的，尺寸分布是怎么算的等。
建议你可以先从最简单的情况，比如喷射多少个尺寸多大的particle，然后一点一点 你就知道什么是什么了。

[好学者]

大佬，你好！我有一些问题想问问，sprayFoam是基于欧拉-拉格朗日框架下进行喷雾模拟，我看了一些论文在用sprayFoam求解器计算喷雾时，其数值方法介绍时会介绍一些控制方程，其中会写 gas phase（论文里的小标题就写的气相控制方程）或Carrier phase的控制方程使用LES或者RANS的湍流模型，然后dispersed phase或Liquid phase的在拉格朗日框架下求解了喷雾包的质量、动量和能量控制方程。那么我想问一下，
1、那么只有在欧拉框架下的连续相才能使用LES或者RANS湍流模型对吗
2、这里的gas phase或Carrier phase的气相是指的什么？是喷雾室环境里的气体（N2）吗？是N2对液体喷雾工质造成的湍流的影响吗？
3、是不是在sprayFoam里的喷雾子模型里设置的喷射工质都是按照拉格朗日粒子进行计算的（一个parcel里有多个相同属性的液滴，然后喷射出大量parcel发展成雾化过程），即液体喷雾工质都是按照离散相计算吗？那么拉格朗日粒子不会产生连续的液柱吗？
4、如果上述问题3是正确的，那么喷射出来是不是就已经是雾化好的小液滴，这里的雾化小液滴和parcel是不是有所区别？
5、在工质喷射过程中，从喷嘴口区域到喷雾发展完全的底部区域如下图：

虚线的区域是不是都是parcel组成，那么这一块都是离散相吗？
6、如果我想在喷雾过程中增加闪沸（即有气泡成核、长大和破碎过程），那么怎么植入气泡呢？我想是在一个parcel里植入气泡，如果我上述说的正确，那么就是在一个parcel（因为里面是同种属性小液滴，是不是可以看做连续相）里增加了一个气相（这样是不是在一个parcel控制体内出现了两相流），那么这个气相是不是也是应该在拉格朗日框架下计算相应的控制方程。这样是不是相当于一个拉格朗日框架下嵌套了另一个拉格朗日框架进行计算？
大佬，可能问题有点密集，如果知道一二可以给我讲讲。感谢！感谢！感谢！路过的大佬有知道的也可以说说，指点一二，三人行，必有我师！再次感谢。

[youhaoyu]

@好学者 我使用DPMFoam，简单解答下你问题，可能有错误
问题1: DPM模型里面连续项也可以用湍流模型，具体在constant里面湍流模型文件可以设置。
问题2： gas phase应该是气体把，你可以看看设置文件他设置了什么，这个东西要么是离散相液滴，要么是连续项气体
问题3： 喷出来的就是parcel，of使用了parcel和particle两个概念，parcel包裹n个paticle，这个n你可以设置，当n为1时候particle就是parcel，你of设置的其实是parcel，不是particle，这是为了n不为1时候，简便计算。比如1w个颗粒，可以简便为5000个包裹，每个包裹2个颗粒。你用的拉格朗日求解器算出来的肯定是离散项啊，你可以在paraview看一下，是一个一个颗粒的，不会产生连续项，你想看连续项只能看看像分数，或者用Euler模型
问题4： 见问题3
问题5：你这一块肯定是有离散项有连续项，你可以通过像分数看看。openfoam的拉格朗日求解器中，离散项和连续项也有像分数的关系，你可以看看东岳老师的普适性方程看看。你这个虚线内部肯定是有连续项目有离散项，这和像分数有关，比如边长为1m的正方体网格，内部有20个1dm的颗粒，那么这个颗粒像分数就是20‰ 这个和双欧拉模型是异曲同工之妙。
问题6 移植气泡也是目前我的工作，这个太难了，我也在仔细研究。你在parcel里面添加气项是什么意思，如果你是想要添加一些离散项属性这个倒好添加。你是不是做类似于空化的啊。气泡成核我不知道，气泡长大我目前在做，包括气泡破碎这些。你这个要求太难了，需要修改openfoam 的源文件，这个玩意可以私下讨论下

[mingyang]

@youhaoyu 谢谢大佬你的解答，我再看看学习学习，对于您的工作气泡长大是否已经完成？很希望和你进行交流2055977489是我的qq号

[lyj_]

@chengan-wang 老师您好，请问您是如何只输出水雾的呀？具体文件是怎么设置的呀？谢谢老师！

[星星星星晴]

@好学者 在 Euler-Lagrange的一些解析，sprayFoam 中说：

大佬，你好！我有一些问题想问问，sprayFoam是基于欧拉-拉格朗日框架下进行喷雾模拟，我看了一些论文在用sprayFoam求解器计算喷雾时，其数值方法介绍时会介绍一些控制方程，其中会写 gas phase（论文里的小标题就写的气相控制方程）或Carrier phase的控制方程使用LES或者RANS的湍流模型，然后dispersed phase或Liquid phase的在拉格朗日框架下求解了喷雾包的质量、动量和能量控制方程。那么我想问一下，
1、那么只有在欧拉框架下的连续相才能使用LES或者RANS湍流模型对吗

• 可能你有点混淆流体力学中的欧拉观点和拉格朗日观点才有这个问题。至少在我个人的经历中，通常，流场都是在欧拉框架下求解的。求解的时候使用了 LES 和 RANS的湍流模型。所以才会有你所说的连续相下使用LES RANS湍流模型

2、这里的gas phase或Carrier phase的气相是指的什么？是喷雾室环境里的气体（N2）吗？是N2对液体喷雾工质造成的湍流的影响吗？

• 个人理解，Carrier Phase 可以算是总称。
• 什么是carrier phase，在多相流里面，最少涉及到2个相（气/液/固）。其中，主要的那一相就是Carrier phase， 另外一个离散的，被包裹的相也就是离散相（相不是项）。例子：bubble column中，水就是carrier phase，气泡就是被包裹的那一相，也就是离散相。比如旋风分离器中，气体就是carrier phase，颗粒就是离散相。所以这个就是相对的。在喷雾中，气相/气体就是carrierphase。当然也有喷雾工作在其他的条件下，所以不一定是你说的N_2

3、是不是在sprayFoam里的喷雾子模型里设置的喷射工质都是按照拉格朗日粒子进行计算的（一个parcel里有多个相同属性的液滴，然后喷射出大量parcel发展成雾化过程），即液体喷雾工质都是按照离散相计算吗？那么拉格朗日粒子不会产生连续的液柱吗？

• 如你所说，所有的离散相都是按照parcel进行计算的，一个parcel中包含np个particle。如同@youhaoyu 同学解释的那样。
• 如你所说，没有连续的液柱。所以在考虑雾化模型的时候就有所谓的TAB，KHRT等模型。如同@youhaoyu 所说。
• 另外，拉格朗日这边，也就是parcel这边，需要通过统计学的方式计算相关的参数。

4、如果上述问题3是正确的，那么喷射出来是不是就已经是雾化好的小液滴，这里的雾化小液滴和parcel是不是有所区别？

• 喷射出来的不一定是雾化好的小液滴。这个要看你的模拟到底位于喷雾的什么区，要看你要模拟什么物理现象。你喷出来的东西是雾化后的参数那就是雾化后的，如果你需要考虑atmozation model的话 那么你就不是雾化好的液滴。要根据你的具体设置具体分析
• 有区别，parcel代表的是np个有相同性质的雾化小液滴的集合。相同尺寸，质量相加

5、在工质喷射过程中，从喷嘴口区域到喷雾发展完全的底部区域如下图：
虚线的区域是不是都是parcel组成，那么这一块都是离散相吗？

• 首先你并没有解释说你这个是哪个field，流场的速度场？如果是流场的速度场的话，那么这个并不是离散相。

• 离散相是那些液滴，需要你额外在paraview中加载的。不知道你使用的那个版本的OF，在OF5以后，加载拉格朗日场会略微有点复杂。建议你在网上找一下解决办法，不同版本有不同或者近似的解决办法。

• 其次，在OF中并没有使用统计的方法把拉格朗日的速度场写出来。所以你需要使用统计的方法，比如在你的domain中添加一个面，统计在1s内通过这个面的颗粒的速度和大小，然后就可以其对应网格的颗粒的平均速度，平均尺寸等信息。不然你在of中查看到的只是瞬时值。

6、如果我想在喷雾过程中增加闪沸（即有气泡成核、长大和破碎过程），那么怎么植入气泡呢？我想是在一个parcel里植入气泡，如果我上述说的正确，那么就是在一个parcel（因为里面是同种属性小液滴，是不是可以看做连续相）里增加了一个气相（这样是不是在一个parcel控制体内出现了两相流），那么这个气相是不是也是应该在拉格朗日框架下计算相应的控制方程。这样是不是相当于一个拉格朗日框架下嵌套了另一个拉格朗日框架进行计算？

• parcel就是point-mass，就是一个点/一个变量/一个c语言中的sturct ，他包含了他的尺寸，密度等信息。但是他在计算的过程中并不会物理上的存在。例如，两个液滴半径都是1cm，在计算的时候，就是两个点p1, p2。 当两个液滴的质心距离效率2cm的时候，在真实世界中，两个液滴就发生了碰撞。但是在模拟中，如果你不考虑碰撞模型的话，两个液滴就是两个点。距离是小于2cm的两个点。没有任何物理现象自然的发生。

• 同理，如果你考虑在parcel中添加气相。你就要考虑到对其相关的物理性质进行修改。例如，原来一个液滴直径1cm，那么p1.d()=100mm, p1.rho() = 998.然后你想到，oh，这个液滴在什么条件下需要添加了气相，那么你就可以对这个parcel添加一个代表有气相的物理相。 还是p1，添加一个参数叫p1.y() = 0.2。当液滴中有了气相以后，相对的就要对其他对应的参数进行修改, p1.d()=105mm, p1.rho()=800。这个东西通常是使用spraycloudproperties文件中的compositionModel来调整和计算的的。可以调整你的parcel中都有什么相，每个相的百分比，相关的热力学参数。这部分我并没有太多的经验，mixture这个东西是很复杂的。

• 并没有所谓的嵌套了另外一个拉格朗日框架，都是在一个框架中的。

总体来说，所谓的欧拉拉格朗日的多相流模拟，简单的来说就是在普通的/常见模拟中添加另外的一个离散相，并将其视为质点。然后看这个离散相相是受到外部流体怎么样的影响。很基础，应用的就是牛二定律 F=ma。然后你就需要对这个parcel所受的力进行建模，常见的drag，lift，bouyancy等（BBO equation）。在不同的情况下应用不一样的模型，比如体积密度比较大， 那么drag就要使用不一样的模型。这个离散相是被流场所影响的就是所谓的1-way coupling，同时离散相也是会对流场产生影响的（反馈：外部流场速度因为parcel而减慢，由于parcel的运动产生的湍流等），如果考虑进去的话就是2-way coupling，即在NS方程中添加一个source term。如果考虑离散相相内部的关系的话（液滴碰撞）就是所谓的4-way coupling。有一点需要你注意，对于离散相的某些特征，比如在某个测量面的数据，请参照我在Q5所说的内容。需要大量的数据统计。

推荐一本书：Computational Fluid Dynamics of Dispersed Multi-Phase Flows - Martin Sommerfeld, ISBN - 978-91-633-3564-8.
你肯定能在什么地方找到。这里面基本上包含了所有你需要的基础知识。

@youhaoyu 感谢你的分享。

[好学者]

@星星星星晴 @youhaoyu 谢谢两位大佬前辈的讲解，我也有了新的认识，感谢感谢！

[youhaoyu]

@星星星星晴 您好大佬，我询问一些事情不知道您曾经是否遇到过。我在加（p.d）颗粒长大时候，假如我颗粒直径是个随着时间线性增长的函数时候，能够进行源文件编译，求解器运行，形如我代码红框那个部分（其他地方是我原始的方程代码，我没有进行使用罢了，storeddt就是我的d）

但是我修改一下代码，改成pv=nrt这种就无法使用，能够通过编译，就是无法进行求解器运行，一下子就报错

**我在想难道是他无法去识别我这个Ro_ pgo_，Pc_？？？**这个Pc是我在求解器里自定义的p*rho。但是我在其他文件下面仿造openfoam自带的一些颗粒属性都对这些进行了定义啊，这个问题困扰了我好几个月了

[星星星星晴]

@youhaoyu 应该是p.d()什么的吧。如果你在kinematicparcel.h等地方定义了的话，你应该通过p.xxx()赋值。另外如果你需要流场速度是可以通过两个方式获得，p.uc()插值的值或者通过mesh获得网格中心的数据吧

你的问题很具体，我也很久没接触OF了。所以建议你多去google。

[星星星星晴]

比如你的particle diameter 应该通过p.d()获得。你的rhoc 是打算使用网格中心值还是考虑是有parcel点的插值？都是有算法的。

整体来说，你这个算法不能这么直接从别的地方生搬硬套。要根据你用的cloud，parcel，根据怎么获取parcrl这个struct内的内容来处理。建议你一步一步来，一行一行写，编译，测试。不能一口吃个胖子

[youhaoyu]

@星星星星晴 比如of他自带Uc插值，是网格中心mesh的，但没有压力P插值，我人为在求解器中间写了个Pc（模仿createField.H的Uc写的），然后在源库里面使用跟着UInterp模仿写PInterp，我这个思路我感觉应该是正确的吧，包括加什么粒子新的属性，我都是模仿他们原有的p.age这些去编写，但是最后就是有问题。
google有什么更好的搜索方式吗，我之前尝试过，请问还有什么平台我可以去询问，cfdonline我打算也问问哈哈哈哈哈

[星星星星晴]

@youhaoyu 你这个想法是对的，我们也是那么植入p的interpolation的。
cfdonline一般没人回 就是多查 有空可以看看chalmers 的openfoam的课，里面有很多学生的project。很有帮助

[星星星星晴]

@youhaoyu 因为你是基于kinematic parcel做的，说不准你做的很多东西都可能在别的地方植入过。
有的时候你觉得你植入对了，但是很有可能你并没有，这边差一点，那边差一点。

而且有一些模型并没有在kinematiccloud中真的植入。有可能就是一个空的头文件。

最好你植入一点，查一点。我建议你捋顺整个of的library，在哪个地方call cloud，cloud在什么地方call 的什么方程，到底call了几次。慢慢查，整着整着就熟了。

我之前可能在某个帖子里面发过我当时对lagranian这部分的理解的流程图，你可以借鉴一下。
不要纠结于具体某个模型，先搞懂框架，然后思考你要植入的东西具体应该植入到哪边，然后通过看看相似的其他植入，具体问题具体分析吧。git有很多repo 有的论文作者会分享出来的。

[Thomas_Huster]

@星星星星晴 大佬你好，就是我按照那个gofun在OpenFOAM9中编译的时候其OpenFOAM4中basicKinematicCollidingCloud.H的文件改名为了collidingCloud.H，但是编译一直告诉我找不到collidingCloud.H

[星星星星晴]

@Thomas_Huster
你在这include了文件，但是在你的make文件里面并没有提供相关的目录和地址，所以程序自然找不到在哪。你需要调整你的make文件的路径