Euler-Lagrange的一些解析，sprayFoam

chengan.wang

@星星星星晴我仔细想了一下，输出文件的方法可能数据量很大，感觉还是转换的方法更合适，
但还是有点没有搞清楚。

您之前告诉我可以放到cloudfunction中，我想知道您具体是放在哪个文件中了，能提供个思路么？我想看看安装包的原始文件，然后再比较一下您的代码，这样可能就清楚该怎么改和调用了。

另外，IVT8_StochasticCollisionModel_2022，PPC这是随便取得名字么？我在使用的时候是不是得改动呢？

星星星星晴

@chengan-wang
不能给你我的code，这个不是我自己的私产，我只能大致的说一下方法。
祝你好运。

chengan.wang

@星星星星晴不是要您的code，就是想知道应该改哪个文件而已，您误会了

星星星星晴

@chengan-wang
正如我前面所说的，你可以在任何Euler/Lagrange场的计算的地方建立这个field，solver部分也可以。
然后在任何Lagrange计算的地方对这个field进行操作。
只要是这个计算的方法是符合你的要求的，Euler 时间步还是Lagrange 时间步。
你可以参考具体下面的拉格朗日部分的流程图，在你需要的地方放置计算即可。

chengan.wang

@星星星星晴非常感谢，我再仔细研究一下。我在google以及源代码中没有查到IVT8_StochasticCollisionModel_2022，PPC，它们代表的意思您方便告诉我么？感觉这个PPC应该是粒子的相关信息，又感觉直接用这个名字不对

星星星星晴

@chengan-wang
ivt-xxxx 这个是我们自己的cloudfunction的名字而已，
ppc也就是这个cloudfunction中的一个function的名字，可以被调用。

你可以尝试在solver中建立一个field，叫任何名字都可以。

香柏树

@星星星星晴前辈，您好！我在使用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;
    }
}



}

星星星星晴

@香柏树

首先你Euler那边速度发生了变化，自然Lagrange这边也会收到影响，我看你的parcel size 是10-60um, 然后parcel的初始速度就是33m/s，不知道你空气速度如何，自然drag的区别很大，然后parcel就可能很快的飞出去了。
在此之上，你有破碎模型，如果气液相对速度过大，自然会破碎，你要具体看一下你的破碎模型是怎么处理质量的。
其次，你的问题中“Case2的液滴数目（nParticles）比Case1少了近一半” 我没有get到，如果指的是particle per parcel的话，自然是受到其他模型的影响，我看你这边有涉及到热，破碎等模型。这些模型也会影响nParticles，在使用这些模型的时候，你也是需要知道这些模型是怎么算的，算什么的，怎么处理mass/momentum balance的。
另外flowrateprofile constant就行，我看你这个flowrate不是随时间变化的，而且就喷1秒。

香柏树

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

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

星星星星晴

@香柏树

spheredrag的前提就是假设所有的parcel 保持球形，所以根据阻力计算公式你可以自己推导一下，到底怎么回事。。严格来说drag 不一定是阻力吧，要看相对速度方向。
你需要查一下关于破碎模型的原始论文，破碎同样也涉及到计算相对速度。在破碎模型中到底是如何假设破碎的，ReitzDiwakar 我不太了解，但是比如TAB，要计算一个y，KHRT也要计算波动的，而且是否生成了child parcel，生成的值是什么，你要看code的。单纯问，没办法解释，看code看原始论文，看of是如何implement的
这个不难理解，如果你流场的速度增加了，那你喷雾的parcel也是收到流场影响的啊，自然可能会发生parcel被加速，如果被加速了的话，有些parcel离开的快也可以理解，你也没有具体描述你的case，所以我只能根据你当前的描述说了。

香柏树

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

chengan.wang

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

@chengan-wang 这也简单啊，同理，你不需要判断是否碰到面，而是做一个判断是否在你的区域内即可，比如你需要判断在两个面之间，z-dir
if(p.position(2)<1 & p.position(2)>0)
{
输出
}
但是这种方法会出现什么问题，比如一个parcel在这个区域待了3个$\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] 相应的位置。
但是存在两个问题：

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

想麻烦您帮忙看看，谢谢

星星星星晴

@chengan-wang

数据量极大就是这个方法最大的问题，没办法。
你的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

数据量极大就是这个方法最大的问题，没办法。

你的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，请教三个个问题：

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

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

lyc

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

youhaoyu

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

星星星星晴

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

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

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

CFD中文网

Euler-Lagrange的一些解析，sprayFoam