使用Python建模二维热传导问题

在本教程中，我们将学习如何使用Python模拟二维热传导方程。带有热生成的二维稳态热传导方程可以用笛卡尔坐标表示如下：

$$\mathrm{\triangledown^{2} T \: + \: \frac{q_{g}}{k} \: = \: \frac{\partial^{2}T}{\partial x^{2}} \: + \: \frac{\partial^{2}T}{\partial y^{2}} \: + \: \frac{q_{g}}{k} \: = \: 0 \:\:\dotso\dotso (1)}$$

这需要离散化才能得到有限差分方程。让我们考虑如下所示的矩形网格。

索引𝑖垂直运行，即行方向，而索引𝑗水平运行，即列方向。任何内部节点(i,j)都被四个节点包围，因此该节点的信息只能从这四个节点传输。而对于边界节点，信息将来自边界条件。

为了离散化方程1，将使用二阶精度的中心差分格式。此外，水平方向和垂直方向的网格点间距保持相同，即$\mathrm{\Delta x \: = \: \Delta y \: = \: \Delta}$。单变量函数二阶导数的中心差分公式为：

$$\mathrm{\frac{d^{2}f(x)}{dx^{2}} \: = \: \frac{f_{i-1} \: − \: 2f_{i} \: + \: f_{i+1}}{\Delta x^{2}}\:\:\dotso\dotso (2)}$$

同样，方程1的离散化可以写成：

$$\mathrm{\frac{T_{i-1, \: j} \: − \: 2T_{i,\: j} \: + \: T_{i+1, \: j}}{\Delta x^{2}} \: + \: \frac{T_{i, \: j-1} \: − \: 2T_{i,\: j} \: + \: T_{i, \: j+1}}{\Delta 6^{2}} \: + \: \frac{q_{g}}{k} \: = \: 0 \:\:\dotso\dotso (3)}$$

方程3的简化形式为：

$$\mathrm{T_{i,\: j} \: = \: \frac{1}{4}(T_{i-1, \: j}\: + \: T_{i+1, \: j} \: + \: T_{i,\: j-1} \: + \: T_{i, \: j+1} \: + \: \Delta^{2} \frac{q_{g}}{k}) \:\:\dotso\dotso (4)}$$

方程4是必须针对域中每个内部节点求解的方程。

对于这类问题，信息从边界传播到内部，因此边界条件的知识至关重要。基本上，在任何数值和分析分析中经常使用三种类型的边界条件：

狄利克雷边界条件：在此条件下，边界上因变量的值已知。例如
$$\mathrm{x \: = \: 0 \: \rightarrow \: T \: = \: 500 \: K }$$
诺依曼边界条件：在此条件下，因变量的导数给出为零或自变量的函数。例如
$$\mathrm{x \: = \: 0 \: \rightarrow \: \frac{\partial T}{\partial x} \: = \: 0 }$$
罗宾边界条件：在此条件下，因变量的导数是因变量本身的某个函数。例如
$$\mathrm{− \: k \: (\frac{\partial T}{\partial x})_{x=0} \: = \: h(T \: − \: T_{\infty}) }$$

算法

选择区域
选择x和y方向的网格点数
定义$\mathrm{\Delta x}$和$\mathrm{\Delta 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")