CFD中文网

@尚善若水 很显然你并没能成功构建matplotlib所定义的三角网格单元，可以看看我在13楼中针对cell变量的定义。直接用三角网格数据绘图的话，此处你直接设定cell = mesh.cells[0].data应该就行（我现在没电脑，无法验证）。

不正确设定cell和node的对应关系的话，tri.Triangulation会根据你输入的点集进行自动Delaunay三角化，造成你所描述的绘制错误。

@尚善若水 是不是有多边形网格？先尝试打开切片三角化吧

我以前对OF加过矢量化的编译指令，但是计算会出错。想用这套东西在设计代码时候就需要考虑，日常我写自己程序的时候用Eigen是支持矢量化的。

我不太清楚你的动网格求解器和相关案例能不能输出poindDisplacement或pointMotionU这类描述节点运动场量的数据。
如果需要凭空获得某个刚体表面运动速度的话，可以使用postProcess -func writeCellCentres获知每个时间步中动网格边界每个面心的坐标，然后根据时间步长去计算运动速度。

可以考虑用codedFunctionObject写一个针对表面的跟踪

@李东岳 哇，感谢李老师

感觉你的模型阻塞比有点大，我一般会压到2.5%以下。

ANSYS Fluent里，针对颗粒的运动方程写作：
$$
m_p \frac{d \vec{u}_p}{d t}=m_p \frac{\vec{u}-\vec{u}_p}{\tau_r}+m_p \frac{\vec{g}\left(\rho_p-\rho\right)}{\rho_p}+\vec{F}
$$
式中， 右侧第一项是曳力项，第二项是浮力+重力，第三项是其他外力。

理论手册里把这个公式简写成一个一阶线性非齐次常微分方程：
$$
\frac{d u_p}{d t}=\frac{1}{\tau_p}\left(u-u_p\right)+a
$$
式中，$\tau_p,a,u$被视为常数。

这个速度的方程如果用欧拉法或者龙格库塔做数值离散进行求解是比较好理解的，但是Fluent还给出了一个名为analytical integration（不定积分？）的求解方法，可以把上式离散成：
$$
u_p^{n+1}=u^n+e^{-\frac{\Delta t}{\tau_p}}\left(u_p^n-u^n\right)-a \tau_p\left(e^{-\frac{\Delta t}{\tau_p}}-1\right)
$$
这个公式看着有点像常微分的通解，但又不知道为什么又会出现$n,n+1$这种表示离散化的上标，有没有大佬知道如何推导？

@yuhxfoamer

1. 首先pathlib这个模块在最终的代码中并没有被使用到，所以你可以将这个模块和对应的代码移除，这并不会影响代码运行结果。
2. 然后谈谈你这个报错。readData.py这个文件是用pvpython执行的，这个python解释器是paraview自带的，对应的第三方库也是在安装paraview时配置好的，所以错误来源于你的paraview附带的pvpython解释器（可能是安装问题，可能是系统环境问题，也可能是软件版本问题，请自行排查），和系统python解释器无关。你使用pip和conda只是在安装你系统python解释器的第三方库。
3. 尽管最终版本这个模块被弃用，不过这里再补充一下：pathlib在本代码中的用途是创建空文件夹（pathlib.Path），在上面的讨论中也提到了，因此可以使用os.mkdir替代。

我不太清楚你表述的网格扭曲是什么意思，如果只是纵横比不正确的话，可以使用set_aspect进行调整。比如我上面代码段中的ax.set_aspect('equal')，表示x轴和y轴始终保持1:1的比例。

@尚善若水 画单纯的等值线可以参考 https://matplotlib.org/3.1.1/gallery/images_contours_and_fields/tricontour_smooth_delaunay.html#sphx-glr-gallery-images-contours-and-fields-tricontour-smooth-delaunay-py

@尚善若水 没能理解你的第二个问题

@尚善若水

其实是这么写的，就是指定读取的时间，后面改着改着就忘了

b = 10

os.system('pvpython readData.py -time_ %f' % b)

你想读取最后一步也很简单。注意第一段代码中的这一行：

times = reader.TimestepValues

这一行代码返回了一个包含所有可读时间的List，List的最后一个元素就是latestTime。因此接下来使用这个元素更新pipeline即可：

UpdatePipeline(time=times[-1])

@刀尔東 自己写代码生成文件应该是最直接的方法

@李东岳 感谢李老师！现在对这个求解器理解了一些，对于$\mathbf{U}_1$的控制方程, $K_d \mathbf{U}_2$一项作为显式项被直接纳入pEqn.H的界面通量的修正。

@李东岳 李老师你好，这边主要不是纠结符号正负，主要是疑惑为什么求解器中的曳力实现是$K_d \mathbf{U}_1$而不是$K_d (\mathbf{U}_1-\mathbf{U}_2)$？

最近研究了一下OpenFOAM-v2206的twoPhaseEulerFoam求解器，注意到在文件pU/Ueqns.H中，动量方程的代码为

U1Eqn =
        (
            fvm::ddt(alpha1, rho1, U1) + fvm::div(alphaRhoPhi1, U1)
          - fvm::Sp(contErr1, U1)
          + MRF.DDt(alpha1*rho1 + Vm, U1)
          + phase1.turbulence().divDevRhoReff(U1)
         ==
          - Vm
           *(
                fvm::ddt(U1)
              + fvm::div(phi1, U1)
              - fvm::Sp(fvc::div(phi1), U1)
              - DDtU2
            )
          + fvOptions(alpha1, rho1, U1)
        );
        U1Eqn.relax();
        U1Eqn += fvm::Sp(Kd, U1);
        fvOptions.constrain(U1Eqn);
        U1.correctBoundaryConditions();
        fvOptions.correct(U1);

我注意到曳力项被写作U1Eqn += fvm::Sp(Kd, U1);

然而按照我的经验，曳力似乎应该被写作
$$
F_D = K_d (\mathbf{U}_1 - \mathbf{U}_2)
$$
也就是后面的速度应该是两相速度的差值。但是好像OpenFOAM所有的多相流求解器都是按照开头代码那样实现的，非常困惑...恳请各位指点一下...

从报错信息看，你使用的编译器是Fluent 2022自带的Clang编译器，而不是MSVC。

离散格式用vanLeer01和linear01之类的格式可以把标量强制限制在0和1之间。