MATLAB程序初始化OpenFOAM颗粒位置

问题引入

在OpenFOAM的颗粒两相流求解器中，我们可以采用manualInjection的方式进行自定义颗粒的初始位置，这个命令十分方便，在CFDEM中也有类似的命令，不过CFDEM中的命令更加强大，我们不仅可以定义颗粒的初始位置，而且还可以定义颗粒的初始速度,唯一的缺点就是CFDEM学起来资料比较少。。。
anyway，这里我们介绍如何在OpenFOAM中实现这个操作。
其核心步骤是在cosntant文件夹下创建kinematicCloudPositions文件，其大致内容如下：

于是自定义粒子的位置的思路就非常简单，我们利用MATLAB创建上面的文件内容即可。

MATLAB编程

这里不对程序做解说，直接放在这里。

%% write by rqli 2023/8/10
%% MATLAB程序初始化OpenFOAM颗粒位置（颗粒位置为均匀分布）
clc;clear;format long%% 各个方向颗粒的数目
NumXofP = 5;NumYofP = 5;NumZofP = 1;
%% 计算域参数
Lx = 0.278;Ly = 0.278;%这里是二维计算域，因此只有高度和宽度
Lz = 0.278;
%% 生成初始颗粒位置，每个点处速度的标准形式为(0.0，0.0，0.0)
position = zeros(NumXofP,NumYofP);
x = linspace(Lx/(2*NumXofP),Lx-Lx/(2*NumXofP),NumXofP);
y = linspace(Ly/(2*NumYofP),Ly-Ly/(2*NumYofP),NumYofP);
[X,Y] = meshgrid(x,y);
if NumZofP == 1postion_x = reshape(X,NumXofP*NumYofP*NumZofP,1);postion_y = reshape(Y,NumXofP*NumYofP*NumZofP,1);postion_z = zeros(NumXofP*NumYofP*NumZofP,1);
elseif NumZofP >1z = linspace(Lz/(2*NumZofP),Lz-Lz/(2*NumZofP),NumZofP); X3d = [];Y3d = [];Z3d = [];for i = 1:length(z)X3d = [X3d;X];Y3d = [Y3d;Y];Z3d = [Z3d;z(i).*ones(NumXofP,NumYofP)];endpostion_x = reshape(X3d,NumXofP*NumYofP*NumZofP,1);postion_y = reshape(Y3d,NumXofP*NumYofP*NumZofP,1);postion_z = reshape(Z3d,NumXofP*NumYofP*NumZofP,1);disp('done')
end%% 将速度场转为txt文件
position_xyz = zeros(NumXofP*NumYofP*NumZofP,1);
position_xyz = num2str(position_xyz);
position_xyz = string(position_xyz);
for i = 1:NumXofP*NumYofP*NumZofPposition_xyz(i) = strcat("(",num2str(postion_x(i))," ",...num2str(postion_y(i))," ",num2str(postion_z(i)),")");
end 
position_file = 'position';
writematrix(position_xyz,position_file);

二维粒子初始位置如下(采用scatter函数)：

三维粒子位置排布如下(采用scatter3函数)：

最后生成的文件如下：

我们把上面的坐标粘贴到kinematicCloudPositions文件里即可。

三维散点图的绘制参考资料