代码供参考
g = 9.8;
rho_S = 2000;
rho_L = 1030;
d_s = 30; % 10~50
C_D = 0.42;
C_S = 0.2; % 0.05~0.3
mu_L = 0.0015;
v_S_tn = sqrt(4*g*d_s*(rho_S-rho_L)/(3*C_D*rho_L));
v_S_th = v_S_tn;
while true
Re_p = rho_L*v_S_th*d_s/mu_L;
alpha_h = 4.7*(1+0.15*Re_p^0.687)/(1+0.253*Re_p^0.687);
v_S_th_old = v_S_th; % 把原来的v_S_th记作v_S_th_old
v_S_th = v_S_tn*(1-C_S)^alpha_h;
if(abs(v_S_th-v_S_th_old)<1e-10)%如果前后两步的v_S_th相差不大就退出循环
break;
end
end
fprintf('v_S_th = %f \n',v_S_th)
结果(本算例):
v_S_th = 15.919848