Euler-Lagrange的一些解析，sprayFoam

[香柏树]

@星星星星晴 前辈，您好！我在使用sprayFoam喷射液滴的时候，在两个不同进气流量的case进行计算，case1为100%空气，case2为200%空气，所有的constant,system里的文件设置都一样，仅仅改变了进气流量。进行液滴后处理时发现，Case2的液滴数目（nParticles）比Case1少了近一半，请问您知道是什么原因么？液滴喷射模型设置具体如下：

model1
        {
            type            coneNozzleInjection;
            SOI             0.029;
            massTotal       0.0044166;
            parcelBasisType mass;
            injectionMethod disc;
            flowType        constantVelocity;
            UMag            33.0;
            outerDiameter   0.0007;
            innerDiameter   0;
            duration        1.0;
            position        (0 0 -0.03584797);
            direction       (0 0 1);
            parcelsPerSecond 2000000;
            massFlowRate    0.0044166;
            flowRateProfile table
            (
                (0              0.0044166)
                (1              0.0044166)
                (2              0.0044166)
                (3              0.0044166)
                (4              0.0044166)
                (5              0.0044166)
                (6              0.0044166)
                (7              0.0044166)
                (8              0.0044166)
                (9              0.0044166)
                (10              0.0044166)
                
            );

            Cd              constant 0.9;

            thetaInner      constant 25;
            thetaOuter      constant 30;

            sizeDistribution
            {
                type        RosinRammler;

                RosinRammlerDistribution
                {
                    minValue        1e-05;
                    maxValue        6e-05;
                    d               4e-05;
                    n               3;
                }
            }
        }

[星星星星晴]

@香柏树 请问100%？和200% 空气是什么意思？另外空气的速度有变化么？你有没有用碰撞模型？

[香柏树]

@星星星星晴 谢谢前辈！100%和200%代表计算域进口空气流量的相对值，比如100%是1kg/s的话，200%就是2kg/s。流场内空气速度增大了两倍左右，没有开启碰撞模型，sprayCloudProperties文件内容具体如下：

\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      binary;
    class       dictionary;
    location    "constant";
    object      SprayCloudProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solution
{
    active          true;
    coupled         true;
    transient       yes;
    cellValueSourceCorrection on;
    maxCo           0.3;

    sourceTerms
    {
        schemes
        {
            rho             explicit 1;
            U               explicit 1;
            Yi              explicit 1;
            h               explicit 1;
            radiation       explicit 1;
        }
    }

    interpolationSchemes
    {
        rho             cell;
        U               cellPoint;
        thermo:mu       cell;
        T               cell;
        Cp              cell;
        kappa           cell;
        p               cell;
    }

    integrationSchemes
    {
        U               Euler;
        T               analytical;
    }
}


constantProperties
{
    T0              300;

    // place holders for rho0 and Cp0
    // - reset from liquid properties using T0
    rho0            744;
    Cp0             2090;
    Tvap            298;
    Tbp             489;

    constantVolume  false;
}


subModels
{
    particleForces
    {
        sphereDrag;
    }

    injectionModels
    {
        model1
        {
            type            coneNozzleInjection;
            SOI             0.029;
            massTotal       0.0044166;
            parcelBasisType mass;
            injectionMethod disc;
            flowType        constantVelocity;
            UMag            33.0;
            outerDiameter   0.0007;
            innerDiameter   0;
            duration        1.0;
            position        (0 0 -0.03584797);
            direction       (0 0 1);
            parcelsPerSecond 2000000;
            massFlowRate    0.0044166;
            flowRateProfile table
            (
                (0              0.0044166)
                (1              0.0044166)
                (2              0.0044166)
                (3              0.0044166)
                (4              0.0044166)
                (5              0.0044166)
                (6              0.0044166)
                (7              0.0044166)
                (8              0.0044166)
                (9              0.0044166)
                (10              0.0044166)
                
            );

            Cd              constant 0.9;

            thetaInner      constant 25;
            thetaOuter      constant 30;

            sizeDistribution
            {
                type        RosinRammler;

                RosinRammlerDistribution
                {
                    minValue        1e-05;
                    maxValue        6e-05;
                    d               4e-05;
                    n               3;
                }
            }
        }

        

    dispersionModel none;

    patchInteractionModel standardWallInteraction;

    heatTransferModel RanzMarshall;

    compositionModel singlePhaseMixture;

    phaseChangeModel liquidEvaporationBoil;

    surfaceFilmModel none;

    atomizationModel none;

    breakupModel    ReitzDiwakar; // ReitzKHRT;

    stochasticCollisionModel none;

    radiation       off;

    standardWallInteractionCoeffs
    {
        type            rebound;
    }

    RanzMarshallCoeffs
    {
        BirdCorrection  true;
    }

    singlePhaseMixtureCoeffs
    {
        phases
        (
            liquid
            {
                C12H26               1;
            }
        );
    }

    liquidEvaporationBoilCoeffs
    {
        enthalpyTransfer enthalpyDifference;

        activeLiquids    ( C12H26 );
    }

    ReitzDiwakarCoeffs
    {
        solveOscillationEq yes;
        Cbag            6;
        Cb              0.785;
        Cstrip          0.5;
        Cs              10;
    }

/*
    ReitzKHRTCoeffs
    {
        solveOscillationEq yes;
        B0              0.61;
        B1              40;
        Ctau            1;
        CRT             0.1;
        msLimit         0.2;
        WeberLimit      6;
    }
*/
    TABCoeffs
    {
        y0              0;
        yDot0           0;
        Cmu             10;
        Comega          8;
        WeCrit          12;
    }
}


cloudFunctions
{
    voidFraction
    {
        type   voidFraction;
    }
}



}

[星星星星晴]

@香柏树

1. 首先你Euler那边速度发生了变化，自然Lagrange这边也会收到影响，我看你的parcel size 是10-60um, 然后parcel的初始速度就是33m/s，不知道你空气速度如何，自然drag的区别很大，然后parcel就可能很快的飞出去了。

2. 在此之上，你有破碎模型，如果气液相对速度过大，自然会破碎，你要具体看一下你的破碎模型是怎么处理质量的。

3. 其次，你的问题中“Case2的液滴数目（nParticles）比Case1少了近一半” 我没有get到，如果指的是particle per parcel的话，自然是受到其他模型的影响， 我看你这边有涉及到热，破碎等模型。这些模型也会影响nParticles，在使用这些模型的时候，你也是需要知道这些模型是怎么算的，算什么的，怎么处理mass/momentum balance的。

4. 另外flowrateprofile constant就行，我看你这个flowrate不是随时间变化的，而且就喷1秒。

[香柏树]

@星星星星晴 感谢前辈悉心指导！关于您说的问题，我今天查阅了一些公式，关于您提的四个问题，具体回复如下：

1. 根据图1中液滴阻力的公式，在100%和200%air中的阻力分别为3.7e-6N和2.3e-5N，也就是液滴受到的阻力变大了，是不是说明液滴更倾向于保持目前的球形状态呢？
2. 破碎模型用的是ReitzDiwakar model,具体公式如图中所示，计算得到100%和200%air中的rstable分别为18μm和1.2μm,也就是说后者的液滴粒径更小。这里您说的破碎模型如何处理质量是什么意思呢？是指公式中的密度么？
3. “Case2中的液滴数目（nParticles）”说错了，应该是流域内的parcels的数目，也就是图2中的8370.计算结果显示200%air case中的parcels数目比100%air case少了近一半。
4. Flowrateprofile constant已经改正。
   再次感谢前辈指点！
   
   

[星星星星晴]

@香柏树

1. spheredrag的前提就是假设所有的parcel 保持球形，所以根据阻力计算公式你可以自己推导一下，到底怎么回事。。严格来说drag 不一定是阻力吧，要看相对速度方向。

2. 你需要查一下关于破碎模型的原始论文，破碎同样也涉及到计算相对速度。在破碎模型中到底是如何假设破碎的，ReitzDiwakar 我不太了解，但是比如TAB，要计算一个y，KHRT也要计算波动的，而且是否生成了child parcel，生成的值是什么，你要看code的。单纯问，没办法解释，看code看原始论文，看of是如何implement的

3. 这个不难理解，如果你流场的速度增加了，那你喷雾的parcel也是收到流场影响的啊，自然可能会发生parcel被加速，如果被加速了的话，有些parcel离开的快也可以理解，你也没有具体描述你的case，所以我只能根据你当前的描述说了。

[香柏树]

@星星星星晴 好的 我查下原始文献。谢谢前辈！

[chengan.wang]

@星星星星晴 在 Euler-Lagrange的一些解析，sprayFoam 中说：

@chengan-wang 这也简单啊，同理，你不需要判断是否碰到面，而是做一个判断是否在你的区域内即可，比如你需要判断在两个面之间，z-dir

if(p.position(2)<1 & p.position(2)>0)
{
输出
}

但是这种方法会出现什么问题，比如一个parcel在这个区域待了3个$\mathrm{\Delta }{t}_L$, 那就输出3次，就是你的输出文件几何倍数的增长。。

您好，我按照您的提示，思路如下：
先是定义了三维空间的范围，以及坐标质量的数组。

        scalar xmin = -0.01;
        scalar xmax = 0.01;
        scalar ymin = 0;
        scalar ymax = 0.1;
        scalar zmin = -0.01;
        scalar zmax = 0.01;
        scalar dxyz = 5e-3;
        const int nx = floor((xmax - xmin)/dxyz);
        const int ny = floor((ymax - ymin)/dxyz);
        const int nz = floor((zmax - zmin)/dxyz);
        int nxc;
        int nyc;
        int nzc;
        scalar xc[nx];
        scalar yc[ny];
        scalar zc[nz];
        scalar mass3D[nx][ny][nz] = {0};

我只想最后时刻输出结果，所以先用if条件判断是否要输出数据

 
    if ((this->owner().db().time().value()) == 0.02)
    {
 
        for (int i=0; i<nx; i++)
        {
            xc[i] = i*dxyz+xmin;
        }
 
        for (int j=0; j<ny; j++)
        {
            yc[j] = j*dxyz+ymin;
        }
 
        for (int k=0; k<nz; k++)
        {
            zc[k] = k*dxyz+zmin;
        }
 
        nxc = floor((p.position().component(0) - xmin)/dxyz);
        nyc = floor((p.position().component(1) - ymin)/dxyz);
        nzc = floor((p.position().component(2) - zmin)/dxyz);
 
        mass3D[nxc][nyc][nzc] += p.nParticle()*p.mass();
 
 
        std::ofstream outfile;
        outfile.setf(ios_base::fixed, ios_base::floatfield);
        outfile.precision(8);
        outfile.open("postProcessingWCA/Mass3D.txt", ios_base::app);
 
        for (int i=0; i<nx; i++)
        {
            for (int j=0; j<ny; j++)
            {
                for (int k=0; k<nz; k++)
                {
                    outfile << xc[i] << tab << yc[j] << tab << zc[k] << mass3D[i][j][k] << tab << tab  << nl;
                }
            }
        }
    }

之后我的想法就是根据粒子的位置计算索引位置 nxc nyc nzc，带入到mass3D[nxc][nyc][nzc] 相应的位置。
但是存在两个问题：

1. 每一个粒子统计之后都输出到了Mass3D.txt，数据量非常大，而我是想把所有的粒子质量统计到相应的空间点上累积；
2. 另外就是定义方式scalar mass3D[nx][ny][nz] = {0}; 以及累积方法mass3D[nxc][nyc][nzc] += p.nParticle()*p.mass();感觉有问题，从结果上看感觉没有把相同位置的粒子质量加到一起。

想麻烦您帮忙看看，谢谢

[星星星星晴]

@chengan-wang

1. 数据量极大就是这个方法最大的问题，没办法。
2. 你的mass3D没有输出，在结束这一步的迭代以后，在下一步迭代会被清空的。
   你这一步已经很靠近方法1了，你不要用什么mass3D这个变量，你在solver中建立一个field， 比如叫summass，然后在cloudfunction的postmove 这个function中直接通过，就可以对一下field进行赋值了。

cellI = p.cell();
scalar& PPC = summass_->primitiveFieldRef() [cellI];
PPC += p.nParticle()*p.mass();

我的技能也是有限，所以没有办法给你解决所有问题，你可以自己试试，找找资料，如果发现更好的办法，也希望你分享出来。

对于方法2，我个人有个想法，你获得你所要区域的cell number， 然后仅仅对这些cell进行输出也可以的。我自己没做过，大家都是摸着石头过河，可能每个人的解决方法并不通用，很case- sensitive，这也是开源的一个弊端。。

祝你好运～

[chengan.wang]

@星星星星晴 在 Euler-Lagrange的一些解析，sprayFoam 中说：

@chengan-wang

1. 数据量极大就是这个方法最大的问题，没办法。
2. 你的mass3D没有输出，在结束这一步的迭代以后，在下一步迭代会被清空的。
   你这一步已经很靠近方法1了，你不要用什么mass3D这个变量，你在solver中建立一个field， 比如叫summass，然后在cloudfunction的postmove 这个function中直接通过，就可以对一下field进行赋值了。

cellI = p.cell();
scalar& PPC = summass_->primitiveFieldRef() [cellI];
PPC += p.nParticle()*p.mass();

我的技能也是有限，所以没有办法给你解决所有问题，你可以自己试试，找找资料，如果发现更好的办法，也希望你分享出来。

对于方法2，我个人有个想法，你获得你所要区域的cell number， 然后仅仅对这些cell进行输出也可以的。我自己没做过，大家都是摸着石头过河，可能每个人的解决方法并不通用，很case- sensitive，这也是开源的一个弊端。。

祝你好运～

谢谢。我在createFields.H建立了一个field:

    volScalarField summass
    (
        IOobject
        (
            "summass",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        mesh,
        dimensionedScalar ("zero", dimMass, 0.0)
    );

我没有理解“在cloudfunction的postmove 这个function中直接通过”是啥意思呢？麻烦了

[流体中二少年]

刚接触sprayfoam，请教三个个问题：

1. 初始场指定区域怎么布置液滴？setFields？
2. 液滴的喷射位置和喷射速度怎么设置为跟随inlet边界？
3. 我模拟的是预气化的液滴，比如C7H16既有液滴也有气体，我在算例里面该怎么区分他们？

谢谢前辈，希望解答或者提供一些解决问题的方向

[lyc]

@chengan-wang 前辈想问一下，GOFUN后来再没出particle simulation的第二讲了么，我看2019年本来有但是取消了

[youhaoyu]

@星星星星晴 你好，最近想向DPMFoam添加一些气泡长大方程，比如RP方程，如何去添加啊，有没有什么思路推荐啊，因为本人实在是太不会编程这些东西了。
最近看了篇博士论文，《水质对梢涡空化初生的影响研究》，上面有一些源文件不知道能不能用，意思是在运动包裹源文件中，添加一些颗粒体积更新代码，再植入自己方程什么的。这个思路能行吗？

[星星星星晴]

@youhaoyu
你好，
我不知道你说的RP方程是什么。

但是如果你是用欧拉拉格朗日法做气泡的话，其实很简单，就是调整parcel的直径就可以了p.d()。
当然你需要在合适的位置去调整。
另外如果欧拉拉格朗日的话，你要考虑cell中parcel的体积，这个是有限制的。我也不记得那个参数是什么了，论坛里之前我好像有聊过。

如果做OpenFOAM的话，不懂编程其实是不行的。
当然你也不需要懂得多深，多看，多想，有基础的方法就行。
祝收敛

[星星星星晴]

@流体中二少年

我刚看到，相信一年多过去你肯定已经解决了。

1. 在0/ 文件夹中，添加lagrangian的文件夹，提供parcel所需的相关信息即可。setFields 应该也行，但是我的第一反应是这个方案。
2. 随inlet边界？就是patchinjection, PatchFlowRateInjection 等喷射模型吧。
3. 在spraycloudproperties中设置的C7H16是parcel的，在外面的/constant文件夹中是外面的气体。在paraview中，parcel是parcel，流场是流场。你可以轻松的看出来。在时间步中的lagrangian文件夹中也能看到每个物质的量的

[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文件，在其中计算，还是通过函数拟合？