Euler-Lagrange的一些解析，sprayFoam

星星星星晴

@chengan-wang
同样可以放到cloudfunction中，对应的改名什么的
第一步写这个function
第二步在constructor里面调用，生成这个field
第三步在你需要的地方调用

这个方法比较复杂，需要多个地方定义，我这仅仅是提供一个思路，可以在任何地方生成这个field，只要在计算的时候能够赋值即可，你也可以在你的solver里面生成这个field，只要在拉格朗日部分能调用到这个field就行。

为什么要插值？我没get到你的点。插值也有很多方法的，另外当你的数据量不够的时候自然会存在很大的误差，这是统计学的问题了。

方法一天然的会损失位置信息，因为网格中储存的数据都是在网格中心，如果做复杂的运算的话，这样会损失掉很多信息，anyway这是你的计算，你可以试试。

我个人觉得方法二足矣应对所有的统计问题，无非就是后处理计算复杂一点，无法paraview中直接连续性的可视化

chengan.wang

@星星星星晴谢谢回复，方法2我倒是能用了，只是我需要输出不仅仅是几个面上的数据，而是统计200x100x100这样的网格节点上的数据，不知道该怎么弄了

星星星星晴

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

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

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

chengan.wang

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

@chengan-wang 这也简单啊，同理，你不需要判断是否碰到面，而是做一个判断是否在你的区域内即可，比如你需要判断在两个面之间，z-dir
if(p.position(2)<1 & p.position(2)>0)
{
输出
}
但是这种方法会出现什么问题，比如一个parcel在这个区域待了3个$\Delta t_L$, 那就输出3次，就是你的输出文件几何倍数的增长。。

您好，如果加判断条件的话应该是在particleCollector代码中加，然后编译，对吧？这个圆柱面或者矩形面的定义就无所谓了。也就是说单独命名一个输出文件。
还有，我的网格数量挺大，200x100x100，我得定义好多个长方体区域，采用循环的办法，判断粒子与每一个区域的相对位置，如果符合条件就保存数据，跳出循环。但如果我的计算时间很长，每个时间点都保存数据，数据量肯定很大，能否只保存最后一个时间步的数据？也就是说我需要传递controlDict的参数，比如

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         5;

deltaT          0.01;

星星星星晴

@chengan-wang 自然可以这就需要你自己研究啦，也有许多其他的方法，比如你获得你要查的目标的cell number，然后只搞这部分就行。

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，这也是开源的一个弊端。。

祝你好运～

CFD中文网

Euler-Lagrange的一些解析，sprayFoam