@尚善若水 很显然你并没能成功构建matplotlib所定义的三角网格单元,可以看看我在13楼中针对cell
变量的定义。直接用三角网格数据绘图的话,此处你直接设定cell = mesh.cells[0].data
应该就行(我现在没电脑,无法验证)。
不正确设定cell和node的对应关系的话,tri.Triangulation
会根据你输入的点集进行自动Delaunay三角化,造成你所描述的绘制错误。
@尚善若水 很显然你并没能成功构建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$这种表示离散化的上标,有没有大佬知道如何推导?
readData.py
这个文件是用pvpython
执行的,这个python解释器是paraview自带的,对应的第三方库也是在安装paraview时配置好的,所以错误来源于你的paraview附带的pvpython解释器(可能是安装问题,可能是系统环境问题,也可能是软件版本问题,请自行排查),和系统python解释器无关。你使用pip
和conda
只是在安装你系统python解释器的第三方库。pathlib.Path
),在上面的讨论中也提到了,因此可以使用os.mkdir
替代。我不太清楚你表述的网格扭曲是什么意思,如果只是纵横比不正确的话,可以使用set_aspect
进行调整。比如我上面代码段中的ax.set_aspect('equal')
,表示x轴和y轴始终保持1:1的比例。
@尚善若水 没能理解你的第二个问题
其实是这么写的,就是指定读取的时间,后面改着改着就忘了
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之间。
@李东岳 macOS官方的C++编译器是Clang,和GCC有一些差异。然后现在苹果产品全线基本都转到基于ARM64架构的Apple Silicon处理器了,导致编译参数和Intel/AMD这些x86/x64的CPU不太一样。
比如有个大坑就是基于苹果M1/M2处理器的Clang无法用传统方法自动捕获浮点异常(https://developer.apple.com/forums/thread/689159),编译OpenFOAM的时候要去掉对应参数,不然会冒出来一大堆警告。
@vbcwl 并行计算下处理这种数据,需要对MPI节点的数据进行归约,比如求和可以写作:
Info << "V = " << V << endl; // 输出第?个MPI节点的V值
Pout << "V = " << V << endl; // 输出每个MPI节点的V值
reduce(V, sumOp<scalar>()); // 对标量求和的归约操作
Info << "V = " << V << endl; // 输出所有MPI节点的V的总和
@一颗鸭蛋 默认情况下,你编译求解器生成的可执行文件会保存到$FOAM_APPBIN
,比如我的默认在<上级文件夹>/OpenFOAM-v2206/platforms/darwinARM64ClangDPInt32Opt/bin
,但是这个App的加载磁盘是只读的,所以会报错。
解决方案有两个:
Make/files
中你的求解器保存路径EXE = $(FOAM_APPBIN)/xxxFoam -> EXE = <新路径>/xxxFoam
然后在环境变量中为$FOAM_APPBIN附加这个路径,以便让OpenFOAM在执行时能搜索到这个文件
export $FOAM_APPBIN=<新路径>:$FOAM_APPBIN
pointMotionU
和pointDisplacement
这些控制动网格节点的边界条件文件的数据类型是pointVectorField
,对应的patch类是Foam::pointPatch
而不是Foam::fvPatch
,这应该就是报错的根本原因。
我的思路如下:
// 获取pointPatch的ID,它和相同边界上的fvPatch的ID是一致的
label patchIndex = patch().index();
// 在objectRegistry随意获取一个存在的物理场
const volScalarField & p
(
this->db().objectRegistry::lookupObject<volScalarField>("p");
);
// 从上述物理场中访问fvMesh中的体心/面心场
const vectorField& Cfp = p.mesh().Cf().boundaryField()[patchIndex];
一般放外面就行,如果要放pimple.loop()里的话,标量方程的时间步长需要除以nOuterCorrector。
直接将壁面的速度修改成运动的速度是否是正确的?目前正在尝试这个方法,但是收敛性似乎很糟糕。
找到的案例似乎都是旋转运动...
MRF1
{
cellZone rotor;
active yes;
origin (0 0 0);
axis (0 0 1);
omega constant 2.2831853;
}
特征尺度是L的话,物体到出口一般取10L,其余包括物体到入口、顶部、左右侧取5L。
@落花风 代码只是一个早期雏形,很多调试功能最后并没有用到,我也懒得删干净了。不必理会。
caseType
参数的作用在上一行的注释中已经给出。caseType=0
指读写decompose
算例,即包含processor*
文件夹的并行算例。caseType=1
指读写reconstruct
算例,即串行或使用reconstructPar
重建后的算例。
driftDensity
是我自己写的某个求解器里的标量场的名称。在标准求解器里,你想读压力就是p
,想读速度就是U
,想读温度就是T
...
我的这个算例用的也是六面体网格,通过一次slice
操作转换成了四边形切片,然后对每个四边形单元进行三角化。请注意理解我的这段代码:
cell_0 = mesh.cells[0].data[:, :3]
cell_1 = np.concatenate(
(mesh.cells[0].data[:, 2:], mesh.cells[0].data[:, [0]]), axis=1)
cell = np.concatenate((cell_0, cell_1), axis=0)
triang = tri.Triangulation(mesh.points[:, 0], mesh.points[:, 2], cell)
就是取每个四边形的0 1 2
和0 2 3
两组顶点,组成两个互不重叠的三角形,然后生成matplotlib自己的三角网格。
Triangulatetheslice=True
。这个功能和paraview图形操作里是一样的,这样出来的切片就直接是三角网格了。就不需要进行上述的额外的三角化操作了。@李东岳 就是下面两段python代码绘制的,用的是matplotlib库。
代码分两个文件,一个读数据,一个绘制。
readData.py
from paraview.simple import *
from pathlib import Path
import sys
from paraview.vtk.numpy_interface import dataset_adapter as dsa
# CaseType: 0 = decomposed case, 1 = reconstructed case.
reader = OpenFOAMReader(FileName="PEE.foam", CaseType=0)
if reader:
print("Successful read.")
else:
print("Failed read.")
time = reader.TimestepValues
# print(time)
argList = sys.argv
# print(argList)
try:
timeIndex = argList.index('-time_') + 1
timeI = float(argList[timeIndex])
except ValueError:
print('ValueError: undefined time parameter.')
exit()
UpdatePipeline(time=timeI)
data = paraview.servermanager.Fetch(reader)
data = dsa.WrapDataObject(data)
slicer = Slice(Input=reader, SliceType="Plane", Triangulatetheslice=False)
slicer.SliceType.Origin = [0, 5.5, 0]
slicer.SliceType.Normal = [0, 1, 0]
Path("slice").mkdir(parents=True, exist_ok=True)
SaveData("/Volumes/RAM Disk/slicePlane_%.2f.e" % timeI, slicer)
print('Successful write.')
plotData.py
import meshio
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.tri as tri
import os
#from matplotlib import ticker, cm
plt.rcParams.update({
# "text.usetex": True,
"font.family": "serif",
"font.serif": ['SimSun'],
# "axes.unicode_minus": False,
"mathtext.fontset": "stix",
"font.size": 16,
"figure.figsize": (12, 5)})
ticklabel_style = {"fontname": "Times New Roman", "fontsize": 14}
b = 10
os.system('pvpython readData.py -time_ 10')
mesh = meshio.read('/Volumes/RAM Disk/slicePlane_%.2f.e' % b)
os.system('rm -r "/Volumes/RAM Disk/slicePlane_%.2f.e"' % b)
logData = np.where(mesh.point_data['driftDensity'] <
1e-12, -12., np.log10(mesh.point_data['driftDensity']))
fig, ax = plt.subplots()
ax.set_aspect('equal')
cell_0 = mesh.cells[0].data[:, :3]
print(mesh.cells[0].data[:, 2:])
print(mesh.cells[0].data[:, 0])
cell_1 = np.concatenate(
(mesh.cells[0].data[:, 2:], mesh.cells[0].data[:, [0]]), axis=1)
cell = np.concatenate((cell_0, cell_1), axis=0)
triang = tri.Triangulation(mesh.points[:, 0], mesh.points[:, 2], cell)
#triang = tri.Triangulation(mesh.points[:,0], mesh.points[:,2])
tpc = ax.tricontourf(triang, logData, cmap='jet', levels=16)
#ax.triplot(triang, 'k-', linewidth=0.5)
for cellI in cell:
ax.plot(mesh.points[:, 0][cellI], mesh.points[:, 2][cellI],'k-', linewidth=0.5)
cbar = fig.colorbar(tpc, format='$10^{%.1f}$', fraction=0.0172, pad=0.06)
cbar.ax.tick_params(labelsize=14)
plt.savefig('test.jpg', dpi=400)
plt.show()
# print(mesh.cell_data)
注意我这里的存储路径/Volumes/RAM Disk/
是macOS风格的,Linux下面可以把临时文件存在/dev/shm/
,这样可以回避外部存储,直接在内存中实现IO。
比较方便的方法是联合paraview的pvpython和你自己的python进行处理:
@疏影横斜水清浅 OpenFOAM不可压求解器的密度是被约掉的,我这里的空气密度取的是1。当时的写法不太严谨,没检查量纲。
算是解决了吧,其实大可不必费力去猜,先写个Info看看这玩意执行的时间点合不合适,然后直接访问fvMesh对象,拿个引用出来就完事了
codeCorrect
#{
volVectorField &U = mesh().lookupObjectRef<volVectorField>("U");
#}
有一项工作需要在每个时间步对物理场进行修正,之前采用的方法是参考fvOptions/corrections
的源码编译自己的库(OpenFOAM-v2206),但是这种方法比较麻烦。现在我注意到coded source的参考文档有这么一句话:
The coded option is available to all primitive field types, i.e.
scalarCodedSource: scalar equations
...
It provides hooks to:include sources: codeAddSup
constrain values before the equation is solved: codeSetValue
apply corrections after the equation has been solved: codeCorrect
看上去似乎能在求解控制方程的前后使用codeSetValue
或codeCorrect
调整场量,从源码看似乎传入了一个场的引用:
template<class Type>
void Foam::fv::CodedSource<Type>::correct
(
GeometricField<Type, fvPatchField, volMesh>& field
)
{
DebugInfo
<< "fv::CodedSource<" << pTraits<Type>::typeName
<< ">::correct for source " << name_ << endl;
updateLibrary(name_);
redirectOption().correct(field);
}
但是最后并没有找到用法...在codeCorrect
代码段里试了试变量名称和field
都显示变量未声明,不知哪位大佬有正确的使用方法?还是我理解有误,这个函数并不能实现修改场量的功能?
@李东岳 可以的,李老师!不过做这部分内容的时候我对OF的认知还比较浅,很多细节可能还比较糙..
我写过类似的东西,网格可以根据当前的场值进行变动,用于研究风场作用下的积雪沉积/侵蚀。代码和算例可以参考[https://github.com/fightingxiaoxiao/driftScalarDyFoam],如果有用的话,也可以引下我的论文[https://doi.org/10.3389/feart.2022.822140]。
大致的思路是这样的:
dynamicMeshDict
中定义速度相关的动网格:dynamicFvMesh dynamicMotionSolverFvMesh;
motionSolverLibs (fvMotionSolvers);
motionSolver velocityLaplacian; //displacementLaplacian;
diffusivity quadratic inverseDistance 3(bottom.snow roof top);
pointMotionU
边界条件为codedFixedValue
bottom.snow
{
type codedFixedValue;
value uniform (0 0 0);
name erosionDeposition;
}
codeDict
嵌入代码,这里比较大的一个坑就是动网格使用节点值进行位移,但是实际计算拿到的都是面心值和体心值,需要做一个插值。另外,deltaH
是我自己求解器里定义的场,需要替换成你自己的内容。erosionDeposition
{
code
#{
label patchIndex = patch().index();
const volScalarField& deltaH = this->db().objectRegistry::lookupObject<volScalarField>("deltaH");
primitivePatchInterpolation facePointInterp(deltaH.mesh().boundaryMesh()[patchIndex]); //初始化插值类
const scalarField& deltaHp = deltaH.boundaryField()[patchIndex];
auto deltaHpp = facePointInterp.faceToPointInterpolate(deltaHp); //面心向节点插值
Info << "Successfully interpolate "<< deltaHpp().size() << " points." << endl;
Info << "Local points number = " << patch().localPoints().size() << endl;
pointField pVec(deltaHpp().size());
const scalar lowerLimit = 0.;
forAll(pVec, i)
{
pVec[i][0] = 0;
pVec[i][1] = 0;
pVec[i][2] = deltaHpp()[i];
if(patch().localPoints()[i][2] + pVec[i][2] * this->db().time().deltaTValue() < lowerLimit)
{
pVec[i][2] = (lowerLimit - patch().localPoints()[i][2])/this->db().time().deltaTValue();
}
}
(*this) == pVec;
#};
codeInclude
#{
#include "fvCFD.H"
#include "primitivePatchInterpolation.H"
#};
codeOptions
#{
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude
#};
@李东岳 是的,现在消费级平台已经开始普及DDR5内存了。苹果这个芯片里的CPU和GPU是共享内存的,所以理论内存带宽很大,但是似乎8个CPU核心无法饱和那么大的带宽。
入手了一台MacBook,在macOS上原生编译的OpenFOAM-v2206居然和Ubuntu ARM虚拟机上的OpenFOAM-10半斤八两。
CPU型号:Apple M1 Max 8P+2E (MacBook Pro 16 2021)
系统:macOS Monterey 12.3 真系统
内存:64GB LPDDR5 6400MHz
OpenFOAM版本:OpenFOAM-v2206
8 166.85
4 277.05
2 467.17
1 798.8
CPU型号:Apple M1 Max 8P+2E (MacBook Pro 16 2021)
系统:Ubuntu ARM 20.04 虚拟机
内存:32GB LPDDR5 6400MHz
OpenFOAM版本:OpenFOAM-10
8 154
4 289.836
2 453.656
1 803.381
如果颗粒刚度比较大的话,可以用重叠网格(overset mesh),也就是背景网格+颗粒周边网格,两套网格重合的地方会进行插值处理。这样可以回避网格重画的难点,不过这东西算二维还行,三维想算快算好比较折磨。
Cyclic似乎只适合几何拓扑一致的网格边界吧,你这种情况用cyclicAMI边界可能会好一些。
@星星星星晴 非常感谢!主要是考虑到便携和静音,拿来开会或者以后授课似乎不错,现在手头的普通笔记本略吵...不过新Air拉满配置和旧Pro14丐版差不多了,目前还在纠结哈哈
@星星星星晴 大佬你好,请问这个Apple M1是Macbook Air 13上的么?最近看到有人已经成功在M1上编译了OpenFOAM,想买个M2的macbook air跑跑小算例顺便远程工作站使用。
Liggghts是一个源自知名分子动力学软件lammps的DEM软件,CFDEM是一个第三方实现的接口,用于实现OpenFOAM和Liggghts的耦合。这应该算是最经典的CFD-DEM方法,和商软里的Fluent-EDEM耦合差不多。但是似乎现在已经停止开发了,对OpenFOAM的支持停留在OpenFOAM 5。
DPMFoam大概算是对DEM的简化,给出了一个Parcel的概念,每个Parcel可以包含N个particle,在计算过程中,计算每对相近Parcel之间的接触并修正速度。具体的碰撞模型好像挺多,我不是特别熟悉。但是由于OpenFOAM的MPI并行框架很难对颗粒做负载平衡,而且本身自带的颗粒搜索算法效率似乎也挺低下的,所以对稠密颗粒流的计算效率不是很好。
MPPICFoam是对DPMFoam的简化(事实上你从源码可以看到,MPPICFoam的代码是直接包含DPMFoam的,只是两者的Parcel不同),在MP-PIC方法的求解过程中,不会对每对Parcel的接触进行判断,而是把拉格朗日颗粒映射进所在的Cell,然后根据每个Cell里所有颗粒的速度/浓度来修正颗粒碰撞后的速度。MP-PIC模型里,粒子的接触模型大概细分为三类,Packing(计算颗粒稠密区的堆积压力梯度),Damping(计算颗粒之间产生的碰撞), Return-to-isotropy(计算颗粒碰撞产生的随机散射)。
因此从计算精度上讲 CFDEM >= DPMFoam > MPPCFoam
从计算效率上讲 MPPICFoam > DPMFoam > CFDEM (没用过CFDEM,属于猜测)
具体使用的时候还是要看需求。
我能想到的一个方法是先用blockMesh画一个比较粗的网格,然后执行decomposePar做分区,最后并行执行refineMesh/refineHexMesh在每个计算节点上细分网格。
不过这个操作我没尝试过,毕竟孩子长这么大最多也就画过1000万网格
只使用物理核心的话差异似乎不大。个人觉得个人PC超线程的好处是在物理核心满负载跑计算任务的时候还能运行点别的东西,比如看看网页写写文档…
我试过用双路EPYC 7742在Windows Server系统下面跑ANSYS Fluent,出现过同样的问题。本质上就是现在双CPU一般都是两个NUMA节点(开了超线程的话Windows任务管理器会显示成4个),每个NUMA节点对应你的一半内存,当你NUMA0节点的CPU想访问NUMA1节点对应的内存,就需要先通信NUMA1节点的CPU。而当你的MPI不支持跨NUMA通信时,你的数据在哪一侧的内存里,哪一侧的CPU才会运行,就会出现双路只能跑一路的情况。当然我不太熟悉计算机的底层知识,这些言论可能不太严谨。
解决这个问题的方法很简单,用支持NUMA通信的MPI协议就行。默认的ibmmpi和msmpi都不行,换成intelmpi就可以。