问题引入
在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文件里即可。
三维散点图的绘制参考资料