使用Python建模二维热传导问题
在本教程中,我们将学习如何使用Python模拟二维热传导方程。带有热生成的二维稳态热传导方程可以用笛卡尔坐标表示如下:
▽2T+qgk=∂2T∂x2+∂2T∂y2+qgk=0……(1)
这需要离散化才能得到有限差分方程。让我们考虑如下所示的矩形网格。

索引𝑖垂直运行,即行方向,而索引𝑗水平运行,即列方向。任何内部节点(i,j)都被四个节点包围,因此该节点的信息只能从这四个节点传输。而对于边界节点,信息将来自边界条件。
为了离散化方程1,将使用二阶精度的中心差分格式。此外,水平方向和垂直方向的网格点间距保持相同,即Δx=Δy=Δ。单变量函数二阶导数的中心差分公式为:
d2f(x)dx2=fi−1−2fi+fi+1Δx2……(2)
同样,方程1的离散化可以写成:
Ti−1,j−2Ti,j+Ti+1,jΔx2+Ti,j−1−2Ti,j+Ti,j+1Δ62+qgk=0……(3)
方程3的简化形式为:
Ti,j=14(Ti−1,j+Ti+1,j+Ti,j−1+Ti,j+1+Δ2qgk)……(4)
方程4是必须针对域中每个内部节点求解的方程。
对于这类问题,信息从边界传播到内部,因此边界条件的知识至关重要。基本上,在任何数值和分析分析中经常使用三种类型的边界条件:
狄利克雷边界条件:在此条件下,边界上因变量的值已知。例如
x=0→T=500K
诺依曼边界条件:在此条件下,因变量的导数给出为零或自变量的函数。例如
x=0→∂T∂x=0
罗宾边界条件:在此条件下,因变量的导数是因变量本身的某个函数。例如
−k(∂T∂x)x=0=h(T−T∞)
算法
选择区域
选择x和y方向的网格点数
定义Δx和Δy。
定义网格。
提供温度的初始猜测值。
提供边界条件
对方程4中的内部节点求解
评估误差。
如果误差大于收敛标准,则将当前T值设置为新的猜测值,并遵循步骤4-9。否则,答案已达到,并绘制结果。
在本教程中,我们将解决两种情况,即有热量产生和无热量产生。此外,我们将假设狄利克雷边界条件。
待解决的问题如下:
问题1

问题2

对于这两个问题,除了热量产生参数外,其他所有参数都保持不变。这就是为什么本文只显示一次代码的原因。
Importing modules from pylab import * from numpy import* # Defining thermal properties # case 1 # qg=0 # case 2 qg= 1000000 k=100 # Define domain ℓx=1.0 ℓy=1.0 # Number of grid points nx=20 ny=20 # Grid spacing Δx=ℓx/(nx-1) Δy=ℓy/(ny-1) Δ=Δx # Grid generation x=linspace(0,ℓx,nx) y=linspace(0,ℓy,ny) X,Y=meshgrid(x,flipud(y)) # Initial Guess T=zeros(shape=(nx,ny)) # Boundary conditions #left T[:,0]=300 #right T[:,-1]=300 #top T[0,:]=800 #bottom T[-1,:]=300 # Guess array for comparison Tg=T.copy() # Initial error or entry in loop error=1 # Iteration counter count=1 # Comparison loop while error>1.E-5: # Sweeping in the domain for i in range(1,nx-1): for j in range(1,ny-1): T[i,j]=(T[i,j-1]+T[i,j+1]+T[i-1,j]+T[i+1,j]+Δ**2*qg/k)/4 # Evaluating and printing error error=sqrt(sum(abs(T-Tg)**2)) figure(1,dpi=300) if count%100==0: print(error) semilogy(count,error,'ko') xlabel("Iteration Counter") ylabel("Error") savefig("Error.jpg") # Updating guess for next iteration Tg=T.copy() # Incrementing counter count+= # Result Plotting (Temperature contour plot) figure(3,dpi=300) cp1=contourf(X,Y,T,10,cmap='jet') colorbar() cp1=contour(X,Y,T,10,colors='k') clabel(cp1,inline=True, fontsize=10) savefig("temp.jpg")
第一种情况的误差和温度图


第二种情况的误差和温度图


现在您可以观察到,热量的产生导致核心温度升高,以至于边界条件的影响无法与之匹配。在问题一中,您观察到影响从顶部开始,向底部和侧面移动。然而,在第二种情况下,边界条件试图施加影响,但是由于热量产生过高,热流的方向现在从中心向外移动。然而,由于问题必须匹配边界条件,因此边界温度保持不变。
结论
在本教程中,二维热传导方程已在Python中建模。有限差分法已被用于解决两个问题:一个是有热量产生,另一个是无热量产生。已经观察到,代码中采用的热量产生大小对最终稳态解的影响不大。