使用 Python 进行集总电容分析
当一个温度极高的物体突然掉入较冷的液体中时，如果假设固体的导电电阻与周围的对流电阻相比非常小，则传热分析称为集总电容分析(如下图所示)。在这里，我们将系统视为一个整体。在这种情况下，我们可以假设块体内部能量的变化率将等于与周围流体的热相互作用。
从数学上讲，这可以写成 −
$$\mathrm{pcV\frac{\partial T}{\partial t} \: = \: − hA(T \: − \: T_{\infty}) \: \dotso\dotso \: (1)}$$
$$\mathrm{\frac{\partial T}{\partial t} \: = \: − \frac{hA}{pcV}(T \: − \: T_{\infty}) \: \dotso\dotso \: (2)}$$
现在必须求解方程 2 以获得温度随时间的变化。
如果我们对方程进行积分。 2 物体的初始温度为 $\mathrm{T_{0}}$，在 0 到 t 的时间内变为任意温度 T，则最终答案为 −
$$\mathrm{\frac{T \: − \: T_{\infty}}{T_{0} \: − \: T_{\infty}} \: = \: exp(− \frac{hA}{pcV}t) \: \dotso\dotso \: (3)}$$
现在这个方程可以按以下方式使用 −
如果需要时间才能达到某个温度 $\mathrm{T_{f}}$ −
$$\mathrm{t \: = \: − \frac{pcV}{hA}In(\frac{T_{f} \: − \: T_{\infty}}{T_{0} \: − \: T_{\infty}}) \: \dotso\dotso \: (4)}$$
如果某个时间需要温度 $\mathrm{t_{f}}$ −
$$\mathrm{T \: = \: T_{\infty} \: + \: (T_{0} \: − \: T_{\infty}) \: \times \: exp (−\frac{hA}{pcV}t_{f}) \: \dotso\dotso \: (5)}$$
让我们借助以下示例来演示这一点 −
示例 1
将一个半径为 1 毫米的钢球置于温度为 $\mathrm{1200^{\circ}C}$ 的露天环境中，温度为 $\mathrm{25^{\circ}C}$。计算钢球冷却至 $\mathrm{100^{\circ}C}$ 所需的时间。 $\mathrm{(k_{b} \: = \: 50 \: W/mK, \: c_{b} \: = \: 500kJ/kgK, \: p_{b} \: = \: 8000kg/m^{3}, \: and \: h_{air} \: = \: 10000W/m^{2}K)}$.
解决方案
由于时间有限，我们将使用公式 4。
# 集总电容分析
from pylab import *
#Input Data
T_inf=25
r=1.E-3
T0=1200
Tf=100
k=50
c=500
ρ=8000
h=10000
#Evaluating Volume and Area
A=(4*pi*r**2)
V=(4/3)*pi*r**3
τ0=ρ*c*V/(h*A)
t=-τ0*log((Tf-T_inf)/(T0-T_inf))
print(f't = {round(t,3)} s')
输出
程序输出将是 −
t = 0.367 s
示例 2
在示例 1 中，0.1 秒后球的温度是多少。
解决方案
在这种情况下，必须求解方程 5。除了公式之外，其他一切都相同。
# 集总电容分析
from pylab import *
# 输入数据
T_inf=25
r=1.E-3
T0=1200
Tf=100
t=0.1
k=50
c=500
ρ=8000
h=10000
# 评估体积和面积
A=(4*pi*r**2)
V=(4/3)*pi*r**3
τ0=ρ*c*V/(h*A)
T=T_inf+(T0-T_inf)*exp(-t/τ0)
print(f'T = {round(T,3)} deg. C')
输出
输出结果为 −
T = 580.031 摄氏度
要理解该图，应绘制温度变化图。对于这种情况，公式 3 的使用方法如下 −
示例 3
from pylab import *
# 输入数据
T_inf=25
r=1.E-3
T0=1200
k=50
c=500
ρ=8000
h=10000
# 计算体积和面积
A=(4*pi*r**2)
V=(4/3)*pi*r**3
# 计算时间常数
τ0=ρ*c*V/(h*A)
t= linspace(0,1.0,50)
T=T_inf+(T0-T_inf)*exp(-t/τ0)
figure(1,dpi=300)
plot(t,T,'r-o')
xlabel('t (s)')
ylabel('T ($^\circ$C)')
savefig('plot1.jpg')
show()
输出
这是显示温度变化的图−
数值方法
如果想要遵循数值方法，那么您也可以这样做。让我们使用显式方法求解方程 2。
对于方程上方显示的网格。 2 可以离散化(使用前向差分)为 −
$$\mathrm{\frac{T_{i} \: − \: T_{i − 1}}{\Delta t} \: = \: −\tau_{0}(T_{i} \: − \: T_{\infty})}$$
$$\mathrm{T_{i} \: = \: T_{i − 1} \: − \: \Delta t \: \times \: (T_{i} \: − \: T_{\infty})/\tau_{0} \: \dotso\dotso \: (6)}$$
下标表示时间步骤。
我们将按时间前进，但在每个时间步骤中，我们将迭代直到解决方案收敛。
由于初始时间温度已知，即 1200，因此将从第二个时间步骤开始。首先，我们将猜测此时间步骤的温度值(例如 0)，然后我们将求解方程 6，然后比较猜测值和获得的值。如果绝对差小于收敛标准 $\mathrm{10^{− 5}}$，则我们将转到下一个时间步骤。否则，我们将将获得的温度值设置为猜测值并再次求解方程 6。该过程重复到最后一个时间步骤。此程序如下所示。我们还将结果与精确的分析值(公式 3)进行了比较。
示例
from pylab import *
# 输入数据
T_inf=25
r=1.E-3
T0=1200
k=50
c=500
ρ=8000
h=10000
# 评估体积和面积
A=(4*pi*r**2)
V=(4/3)*pi*r**3
τ0=ρ*c*V/(h*A)
n=35 # 时间步数
t = linspace(0,1.0,n)
Δt=1/(n-1)
# 所有时间步的初始猜测
T=zeros(n)
# 初始条件
T[0]=T0
# 猜测数组
Tg=T.copy()
# 时间推进
for i in range(1,n):
# 迭代循环
error=1
while error>1.E-5:
T[i]=T[i-1]-Δt*(T[i]-T_inf)/τ0
error=abs(T[i]-Tg[i])
Tg=T.copy()
# 精确解
T_exact=T_inf+(T0-T_inf)*exp(-t/τ0)
# 数据绘图
figure(2,dpi=300)
plot(t,T,'r-o',label='Numerical')
plot(t,T_exact,'b-o',label='Analytical')
xlabel('t (s)')
ylabel('T ($^\circ$C)')
legend()
savefig('plot2.jpg')
show()
输出
程序输出如下 −
生成上述图表所依据的数据如下所示 −
time T(analytical) T(numerical) 0.0 1200.0 1200.0 0.0294 987.6506 967.4051 0.0588 813.6776 780.853 0.0882 671.1455 631.2296 0.1176 554.3722 511.2245 0.1471 458.7025 414.9749 0.1765 380.3226 337.7781 0.2059 316.1076 275.8627 0.2353 263.4978 226.2036 0.2647 220.3958 186.3748 0.2941 185.0833 154.4301 0.3235 156.1526 128.809 0.3529 132.4503 108.2597 0.3824 113.0316 91.7782 0.4118 97.1223 78.5592 0.4412 84.0881 67.957 0.4706 73.4095 59.4535 0.5 64.6608 52.6334 0.5294 57.4932 47.1632 0.5588 51.6209 42.776 0.5882 46.8099 39.2572 0.6176 42.8684 36.4349 0.6471 39.6391 34.1713 0.6765 36.9935 32.3558 0.7059 34.826 30.8997 0.7353 33.0502 29.7319 0.7647 31.5954 28.7952 0.7941 30.4034 28.0439 0.8235 29.4269 27.4414 0.8529 28.6269 26.9581 0.8824 27.9714 26.5705 0.9118 27.4344 26.2596 0.9412 26.9945 26.0103 0.9706 26.634 25.8103 1.0 26.3387 25.6499
结论
本教程使用 Python 解释并建模了集总电容分析。讨论了分析方法和数值方法。
