利用MATLAB计算颗粒在静止流体中的沉降速度

下面是MATLAB代码的内容：

clear;
clc;
close all;

% 流体粘度 kg/(m.s)
% Cate 等人：0.373
miu_f = 1.01e-03
% 液体密度
% Cate等人：960
rho_f = 1000
% 颗粒密度
% Cate等人：1120
rho_p = 2000
% 密度差 
delta_rho = rho_p - rho_f
% 重力
g = 9.81
% 参数 P
P = ( 3/4 * miu_f^2 / ( delta_rho * rho_f * g))^(1/3)
% 参数 Q
Q = ( 4/3 * delta_rho * miu_f * g / (rho_f * rho_f))^(1/3)
% 颗粒直径 d
% Cate等人：0.015
d = 0.030
% 无量纲颗粒粒径d_star
d_star = d / P
%  无量纲颗粒沉降速度 u_star
u_star = 20.52 / d_star * ((1+0.0921*d_star^1.5)^0.5-1)^2
% 颗粒沉降速度 u
u = u_star * Q

上面的单位均为SI单位。

将水的粘度1.01e-03，水的密度1000，颗粒的密度2000，颗粒的直径0.02，带入上面进行计算得到，d_star = 702.1711, u_star = 47.7176。结果接近Fig.6中的结果。

 参考文献：

[1] 刘德天,傅旭东,王光谦.CFD-DEM耦合计算中的体积分数算法[J].清华大学学报(自然科学版),2017,57(07):720-727.DOI:10.16511/j.cnki.qhdxxb.2017.25.028.

[2] Concha Arcil F. Settling Velocities of Particulate Systems [J]. KONA Powder and Particle Journal, 2009, 27: 18-37.