CFD中文网

地形网格cfmesh一键生成真的很方便，我前段时间做大尺度地形界面泥石流模拟，用这个真的省了不少事儿。

据我所知目前openfoam都是pointparticle，没有resolve方法，dev版本里有重叠网格法，适合单个或者几个颗粒的计算。resolve法的话cfdem可以实现，当然找找的话也有其他的开源求解器，不过不是基于openfoam的。

看这个散点图，每0.5s一个，是不是0.5s保存一组数据，到云图里面测量的？

openfoam自带的颗粒处理方法是点源(point model)，无法解析出颗粒的边界。

@杨英狄 目前是已经实现了，文章是前面那个题为"An optimized Eulerial-Lagrangian..."的文章。具体要请教什么呢？

@杨英狄 这个在肖恒老师的两篇文章详细讨论过吧，算是一种经典的处理粗颗粒的办法:
1：Diffusion-based coarse graining in hybrid continuum–discrete solvers: Theoretical formulation and a priori tests
2：Diffusion-based coarse graining in hybrid continuum–discrete solvers: Applications in CFD–DEM
然后也有很多这个方法的变种（基于OpenFOAM或者CFDEM）
1：An optimized Eulerian–Lagrangian method for two-phase flow with coarse particles: Implementation in open-source field operation and manipulation, verification, and validation
2: A semi-resolved CFD-DEM approach for particulate flows with kernel based approximation and Hilbert curve based searching strategy

@李东岳 感觉和相变类似，但是我理解来看相变是相间的物质交换，这里是边界处的输入，单纯的是泥石流相的物质增加。我在想能不能实现一个边界条件，先判断坡面上哪里有泥石流（边界面相分数判断），然后把边界面有相分数的的网格面提取出来，对这些网格面改变$\alpha$的值。

目前正在用OpenFOAM做坡面泥石流相关的数值模拟，主要是基于InterMixingFoam求解器（interFoam的衍生），将泥石流看作两相流体的混合物（颗粒与水），然后气相用来捕捉泥石流混合物运动过程的自由面。这是大致的示意图。

前期修改求解器实现了泥石流内部两相的相互作用（考虑速度滑移，interFoam本身未考虑速度滑移），预测了现场泥石流的运动过程，这是计算的几何域和计算的泥石流运动的速度云图。

图二 计算几何域

图三 泥石流运动速度云图

现在想要进一步延申。之前是直接将底部边界设置为noSlip边界，但是泥石流沿着坡面向下运动时会携带破面的物质向下走，导致泥石流相的体积会越来越大，文献中叫“坡面侵蚀”，想在数值上实现这个过程，我理解来看应该是下壁面上的一种边界条件，但是因为没有接触过类似的问题，现在不知道该用什么边界条件。

所以把问题发出来想咨询一下有没有研究类似问题或者类似过程的大佬可以给指个方向

@oitocfd 之前写过一个小函数，具体是在dpmfoam里用到的，是寻找一个颗粒所在的cell label以及这个cell周围的一圈cell label，不知道能不能帮到你。

template<class Type>
const Foam::labelList& Foam::interpolationCellAround<Type>::cellpointCells
(
    const label celli,
    DynamicList<label>& storage
) const
{
        const labelList& cPoints = this->psi_.mesh().cellPoints()[celli];

        storage.clear();

        forAll(cPoints, i)
        {
            const label pointi = cPoints[i];

            const labelList& pointcellsi = this->psi_.mesh().pointCells()[pointi];

            forAll(pointcellsi, j)
            {
                storage.append(pointcellsi[j]);
            }
        }

        // Filter duplicates
        if (storage.size() > 1)
        {
            sort(storage);

            label n = 1;
            for (label i = 1; i < storage.size(); i++)
            {
                if (storage[i-1] != storage[i])
                {
                    storage[n++] = storage[i];
                }
            }

            // truncate addressed list
            storage.setSize(n);
        }

        return storage;
}

template<class Type>
const Foam::labelList& Foam::interpolationCellAround<Type>::cellpointCells(const label celli) const
{
    return cellpointCells(celli, labels_);
}

interpolationCellAround是自定义的类，主要是函数的实现。

感谢两位大佬的帮助

计算滑块入水问题，初始想在计算域框选出一个三角形区域，不知道能否实现？
![Y}}FNTDQL8G3J4T)7OM9M2.png

不知道这篇论文能否给点帮助（同学的论文）：
"叶秉晟. 针对非定常空化流动的多相可压缩求解器开发及应用[D]. 中国科学院大学，2019."

@李东岳 谢谢李老师，我再详细看一下推导过程

@mechanicsdog 谢谢您，万分感谢

@李东岳 读了东岳老师的这篇文章，不知道公式（22）的量纲是否有点问题，按照公式(21)，弛豫时间的$\tau$量纲应该是$[M^{-1}L^{3}T^{1}]$,那么公式（22）右边第二项的两处量纲都是对不上的，大括号内$\mathbf{g}$那一串不是速度量纲，右边的$\tau$应该是少一个-1次方。另外，请问李老师，您文章onewaycoupledtestcase里面计算颗粒分离时，用的弛豫时间$10d^{2/3}$,请问一下这部分该怎么理解，2/3次方表示的是什么呢？另外，弛豫时间一般理解，量纲不是时间量纲么？谢谢李老师。
贴出您文章的公式(21)和(22)

@mechanicsdog 您好，我想再请教您一个问题，我按照您的方法实现颗粒上显示速度了，现在我在算一个带颗粒级配的流化床问题，颗粒直径是高斯分布，请问后处理的时候您知道如何显示颗粒的大小么，现在按照高斯点好像只能显示一样大小的颗粒，如果想这么显示您知道如何操作么？谢谢您

我也遇到过相同的问题，平均压力没问题，但是瞬时波动很大，暂时没有找出原因所在，难道是算法的问题么

想问作者一个问题，请问是如何在paraview里面把颗粒上色，显示不同的速度的大小的呢。
此外，是否压力本来就该取平均值，在controlDict里面添加平均库函数来获得平均场，在这里贴上我的用过的例子：

functions
{
    fieldAverage
    {
        type          fieldAverage;
        libs          ("libfieldFunctionObjects.so");
        writeControl  writeTime;
        timeStart     1;
        fields
        (
            UField
            {
                mean          on;
                prime2Mean    on;
                base          time;
            }
                     
            alpha
            {
                mean          on;
                prime2Mean    on;
                base          time;
            }
            
            alpha.air
            {
                mean          on;
                prime2Mean    on;
                base          time;
            }
            
            U.air
            {
                mean          on;
                prime2Mean    on;
                base          time;
            }
        );
    }
}

@BlookCFD 谢谢您的指点，我有看了一眼程序，我理解错了，OF里面是通过计算barycentric displacement来获得barycentric coordinates的，我一直被src/lagrangian/basic/particle/particle.C里面的一句代码给误导了，现贴出来：

particle.C 1062-1096 OpenFOAM6
void Foam::particle::correctAfterInteractionListReferral(const label celli)
 {
     // Get the position from the barycentric data
     const vector pos(coordinates_.b(), coordinates_.c(), coordinates_.d());
 
     // Create some arbitrary topology for the supplied cell
     celli_ = celli;
     tetFacei_ = mesh_.cells()[celli_][0];
     tetPti_ = 1;
     facei_ = -1;
 
     // Get the reverse transform and directly set the coordinates from the
     // position. This isn't likely to be correct; the particle is probably not
     // in this tet. It will, however, generate the correct vector when the
     // position method is called. A referred particle should never be tracked,
     // so this approximate topology is good enough. By using the nearby cell we
     // minimize the error associated with the incorrect topology.
     coordinates_ = barycentric(1, 0, 0, 0);
     if (mesh_.moving())
     {
         Pair<vector> centre;
         FixedList<scalar, 4> detA;
         FixedList<barycentricTensor, 3> T;
         movingTetReverseTransform(0, centre, detA, T);
         coordinates_ += (pos - centre[0]) & T[0]/detA[0];
     }
     else
     {
         vector centre;
         scalar detA;
         barycentricTensor T;
         stationaryTetReverseTransform(centre, detA, T);
         coordinates_ += (pos - centre) & T/detA;
     }
 }

其中的这句代码

coordinates_ += (pos - centre) & T/detA;

其实前面已经先定义了

coordinates_ = barycentric(1, 0, 0, 0);

后面计算单点的barycentric coordinates其实还是用的barycentric displacement的概念来转换的。

困扰了我两周的问题终于解决了，现在异常开心，在此对@东岳 和@BlookCFD 表示万分感谢。