多孔介质

import numpy as np
import matplotlib.pyplot as plt

########## Begin User Inputs ##########
fluid_viscosity    = 1.781e-5  # 流体的动态粘滞度 [N·s/m2]
fluid_density      = 1.225     # 流体密度 [kg/m^3]
porous_zone_length = 0.1       # 多孔区域长度 [m]

# 压力与速度数据
velocity = np.array([0.0, 0.1, 0.2, 0.3])          # 流体速度 u [m/s]
pressure_loss = np.array([0.0, 0.01, 0.02, 0.04])  # 多孔区域压力差 ΔP [N/m2]

########## End User Inputs ##########

# 使用 ΔP = A * u + B * u^2 进行抛物线拟合
coefficients = np.polyfit(velocity, pressure_loss, 2)
B = coefficients[0]
A = coefficients[1]

delta_p_fit = A * velocity + B * velocity**2
plt.scatter(velocity, pressure_loss, color='red', label='data')
plt.plot(velocity, delta_p_fit, label=f'fitting: ΔP = {A:.4f} * u + {B:.4f} * u^2', color='blue')

plt.xlabel('Velocity u [m/s]')
plt.ylabel('Pressure Loss ΔP [Pa]')
plt.title('Pressure Loss vs. Velocity')
plt.legend()
plt.grid(True)
plt.show()

# 计算 d 时将 A 限定为正值
d = max(0, A) / (fluid_viscosity * porous_zone_length)
f = 2 * B / (fluid_density * porous_zone_length)

# 打印系数
print(f"Curve-fit Coefficients A = {A:.4f}, B = {B:.4f}")
print(f"Viscous  Resistance = {d:.6f} 1/m^{2}")
print(f"Inertial Resistance = {f:.6f} 1/m")

然后，用自己的数据替换 fluid_viscosity、fluid_density、porous_zone_length、velocity = np.array 系列，和pressure_loss = np.array 系列。

或者

输入以下代码，运行穿孔板计算：

import numpy as np

########## Begin User Inputs ##########

# 穿孔板的物理参数
porosity              = 0.8  # porosity, dimensionless, 0 < eps <= 1. Do not use 0.
hole_diameter         = 0.003  # 圆孔直径 [m]
porous_zone_thickness = 0.025  # 多孔区域厚度 [m]

########## End User Inputs ##########

K = porosity * hole_diameter * hole_diameter * porous_zone_thickness / (32 * porous_zone_thickness + 15 * hole_diameter)

laminar_flow = True
if (laminar_flow): 
  # 如果流动是层流
  alpha = 3 * (1 - porosity) / (4 * porosity * porosity * porous_zone_thickness)
else: 
  # 如果流动是湍流
  delta_over_D = porous_zone_thickness / hole_diameter
  alpha = 9 * (6 * delta_over_D - 5 * delta_over_D * delta_over_D) / (40 * porosity * porosity * porous_zone_thickness)

viscous_resistance = 1 / K
inertial_resistance = 2 * max(0, alpha)

# 打印系数
print(f"Coefficients K = {K:.8f}, alpha = {alpha:.8f}")
print(f"Viscous Resistance = {viscous_resistance:.6f} 1/m^{2}")
print(f"Inertial Resistance = {inertial_resistance:.6f} 1/m")

然后，用自己的数据替换 porosity, hole_diameter，和 porous_zone_thickness。