提问人：shammas mohamed 提问时间：4/11/2023 最后编辑：beakershammas mohamed 更新时间：4/14/2023 访问量：87
我很难将流体动力学方程转换为MATLAB代码
I am having a hard time converting fluid dynamics equation into MATLAB code
问：
我之前将这两个方程转换为MATLAB代码，但是我丢失了文件，现在我正在尝试再次编写它们，但是有一个问题，我无法弄清楚。
这是主等式：
这是平衡半径方程：
这是我从丢失的MATLAB代码中得到的图形。
这是我现在正在编写的MATLAB代码：
clc
%Converting to SI units
Contact_angle = 142*pi/180;
Volume = 5* 1e-9;
ST_L = 0.072 ;
Density = 997 ;
Gravity = 9.807;
Shape_factor = 37.1;
T_zero = 0; Dynamic_viscosity = 8.9e-4 ;
R_e = ((4*Volume)/(pi*Contact_angle))^(1/3);
T = 0:0.1:5.2;
R = R_e * (1 - exp(((-(2*ST_L)/R_e^12) + ((Density*Gravity)/(R_e^10))) * ((24*Shape_factor*Volume^4 * (T+T_zero))/((pi^2)*Dynamic_viscosity)))).^(1/6);
%Convert R to mm R = R*1000;
plot(T,R);
grid on xlim([0 6]);
ylim([0 6]);
legend('Water','Location','northwest');
ylabel('Radius (mm)');
xlabel("Time (s)");
我试着玩弄形状因子（因为它是一个任意常量）并玩弄单位，以获得与我之前的代码类似的结果。然而，当我尝试使用不同的流体来测试我的数学模型时，我意识到改变单位或常数不是一个合适的解决方案。
答：
-1赞
duffymo
4/12/2023
#1
我尝试使用 Python 复制您的代码并使用 matplotlib 显示它，但没有成功。
y 数组中的值均为 。r 数组中的结果值从零跳到并保持在那里，而不会改变。我不确定我做错了什么。O(1.e-30)re = 1.3695297485934952
# see https://stackoverflow.com/questions/75985316/i-am-having-a-hard-time-converting-fluid-dynamics-equation-into-matlab-code
import math
import matplotlib.pyplot as plt
import numpy as np
if __name__ == '__main__':
theta = 142*math.pi/180
v = 5e-9
gamma_L = 0.072
rho = 997
g = 9.807
lamb = 37.1
t0 = 0
mu = 8.9e-4
t = np.arange(0.0, 2.0, 0.01)
re = (4*v/math.pi/theta)**(1/3)
x = 2*gamma_L/(re**12) + rho*g/9.0/(re**10)
y = [(24*lamb*(v**4)*(ti+t0)/(math.pi**2)/mu) for ti in t]
r = [(1000.0 * re * (1 - math.exp(-x*yi))**(1/6))for yi in y]
fig, ax = plt.subplots()
ax.plot(t, r, linewidth=2.0)
ax.set(xlim=(0, 5),
ylim=(0, max(r)))
plt.show()
评论
0赞
shammas mohamed
4/17/2023
谢谢大家的帮助。但我想通了
0赞
duffymo
4/17/2023
告诉我们你“想通了”什么。我想知道。
1赞
John BG
4/14/2023
#2
我得到与 duffymo 相同的稳定结果：
1.- 具有r_e近似值
theta_e = 142*pi/180;
V = 5* 1e-9;
gamma_L = 0.072 ;
rho = 997 ;
g = 9.807;
lambda = 37.1;
t_zero = 0;
etha = 8.9e-4 ;
% r_e=(6*V*(cos(theta_e/2))^3/(pi*(2+cos(theta_e))*sin(theta_e/2)))^(1/3)
r_e = (4*V/(pi*theta_e))^(1/3);
t = [0:0.1:5];
R =r_e*(1 - exp(-(2*gamma_L/r_e^12 + rho*g/(9*r_e^10)) *...
24*lambda*V^4*(t+t_zero)/(pi^2*etha))).^(1/6)
R =
Columns 1 through 4
0 0.0014 0.0014 0.0014
Columns 5 through 8
0.0014 0.0014 0.0014 0.0014
...
Columns 45 through 48
0.0014 0.0014 0.0014 0.0014
Columns 49 through 51
0.0014 0.0014 0.0014
2.- 具有精确的r_e表达式
theta_e = 142*pi/180;
V = 5* 1e-9;
gamma_L = 0.072 ;
rho = 997 ;
g = 9.807;
lambda = 37.1;
t_zero = 0;
etha = 8.9e-4 ;
r_e=(6*V*(cos(theta_e/2))^3/(pi*(2+cos(theta_e))*sin(theta_e/2)))^(1/3)
% r_e = (4*V/(pi*theta_e))^(1/3);
t = [0:0.1:5];
R =r_e*(1 - exp(-(2*gamma_L/r_e^12 + rho*g/(9*r_e^10)) *...
24*lambda*V^4*(t+t_zero)/(pi^2*etha))).^(1/6)
R =
1.0e-03 *
Columns 1 through 4
0 0.6600 0.6600 0.6600
Columns 5 through 8
0.6600 0.6600 0.6600 0.6600
...
Columns 45 through 48
0.6600 0.6600 0.6600 0.6600
Columns 49 through 51
0.6600 0.6600 0.6600
3.- 将 R 带入线性变化的变化
R =1e6* r_e*(1 - exp(-(2*gamma_L/r_e^12 + rho*g/(9*r_e^10)) *...
24*lambda*V^4*(1e-12)*(t+t_zero)/(pi^2*etha))) % .^(1/6)
%Convert R to mm R = R*1000;
plot(t,R);
grid on
ylim([0 6]);xlim([0 6]);
legend('Water','Location','northwest');
ylabel('Radius (mm)');
xlabel('Time (s)')
但这是线性的，而不是指数的，当加倍范围时，会不断增加，但似乎最终应该达到“带电电容器”状态。RtR
在我看来，输入参数中可能有一个或多个值可能需要检查，以防万一。
前面的 1e6 暗示值可能需要进一步评估。r_er_e
指数意味着问题中的整体表达式将其全部推离“电容器充电”斜率发生的地方很远。1e-12
希望对你有所帮助。
评论