使用 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}}$ 到某个任意温度 T 在时间 0 到 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。

# lumped capacitance analysis
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。除了公式外,其他一切都是相同的。

# lumped capacitance analysis
from pylab import *

# Input Data
T_inf=25
r=1.E-3
T0=1200
Tf=100
t=0.1

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=T_inf+(T0-T_inf)*exp(-t/τ0)

print(f'T = {round(T,3)} deg. C')

输出

输出将是:

T = 580.031 deg. C

要理解图像,应该绘制温度的变化曲线。在这种情况下,将使用方程 3,如下所示:

示例 3

from pylab import *

# Input Data
T_inf=25
r=1.E-3
T0=1200

k=50
c=500
ρ=8000
h=10000

# Evaluating Volume and Area
A=(4*pi*r**2)
V=(4/3)*pi*r**3

# Evaluating time constant

τ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 *

# Input Data
T_inf=25
r=1.E-3
T0=1200

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)

n=35 # Number of time steps

t = linspace(0,1.0,n)

Δt=1/(n-1)

# initial guess for all time steps
T=zeros(n)

# initial condition
T[0]=T0

# Array of guess
Tg=T.copy()

# Time marching
for i in range(1,n):

   # Iteration loop
   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()

# Exact solution
T_exact=T_inf+(T0-T_inf)*exp(-t/τ0)
# Data Plotting
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 解释并建模了集总电容分析。讨论了分析和数值方法。

更新于: 2023 年 10 月 3 日

199 次查看

开启您的 职业生涯

通过完成课程获得认证

开始学习
广告