plt python 自己制定cmap_Python数值优化：使用Euler法求解二维热传导方程2

上一篇实践了简单的二维热传导抛物线方程求解，对于简单的热源和等温边界条件进行了模拟，下面来探究一些更有趣的现象。

pizh12thu：Python数值优化：使用Euler法求解二维热传导方程1​zhuanlan.zhihu.com

移动热源

移动热源在焊接分析、3D增材制造等瞬态分析中十分常见。先前计算式中给了恒定常数的热源

,由于Numpy的广播机制，

与矩阵进行相加操作时会被扩展成和U形状相同的矩阵进行运算，为了控制热源的移动，需要对矩阵形式的热源进行操作。

U = U + rx*np.dot(U,A) + ry*np.dot(B,U) + heat*ft

在Python代码中实现，指定一个具有空间位置的热源

import numpy as np
import time
import matplotlib.pyplot as plt

if __name__ == "__main__":
    a = 1
    D = np.array([0,15,0,15])
    Mx = 150
    My = 150
    T,Tn = 20,5
    Nt = 20000

    dx = (D[1]-D[0])/Mx
    dy = (D[3]-D[2])/My
    dt = T/Nt
    Frame = 100
    fq = 10
    u0 = 20
    u_env = 20

    ### 初始化U、A、B矩阵，各点初始温度u0
    U = u0*np.ones((My+1,Mx+1))
    # x方向二阶导系数矩阵A
    A = (-2)*np.eye(Mx+1,k=0) + (1)*np.eye(Mx+1,k=-1) + (1)*np.eye(Mx+1,k=1)
    # y方向二阶导系数矩阵B 
    B = (-2)*np.eye(My+1,k=0) + (1)*np.eye(My+1,k=-1) + (1)*np.eye(My+1,k=1)
    
    #用网格化映射矩阵坐标信息
    x=np.linspace(D[0],D[1],U.shape[1])
    y=np.linspace(D[2],D[3],U.shape[0])
    X,Y = np.meshgrid(x,y)
    x0,y0 = 0,7 #热源中心位置

    ### 按时间增量逐次计算
    for k in range(Nt+1):
        tt = k*dt
        #热源沿着X方向移动
        x0 += 2*dt
        Qv = 80*fq*np.exp(-8*((X-x0)**2+(Y-y0)**2))
        # solve inside nodes
        U = U + rx*np.dot(U,A) + ry*np.dot(B,U) + Qv*dt
        # solve boundary nodes 
        U[:,0] = U[:,-1] = u_env 
        U[0,:] = U[-1,:] = u_env 

用等温云图显示一下计算结果

        if k%Frame == 0:
            end = time.time()
            print('T = {:.3f} s  process time = {:.1f}'.format(tt,end-start))
            showcontourf(U,D,vmax=120)


def showcontourf(mat,D,cmap=plt.cm.get_cmap('jet'),fsize=(12,12),vmin=0, vmax=100):
    plt.clf()
    levels = np.arange(vmin,vmax,1)
    x=np.linspace(D[0],D[1],mat.shape[1])
    y=np.linspace(D[2],D[3],mat.shape[0])
    X,Y=np.meshgrid(x,y)
    z_max = np.max(mat)
    i_max,j_max = np.where(mat==z_max)[0][0], np.where(mat==z_max)[1][0]
    show_max = "U_max: {:.1f}".format(z_max)
    plt.plot(x[j_max],y[i_max],'go')  
    plt.contourf(X, Y, mat, 100, cmap = cmap, origin = 'lower', levels = levels)
    plt.annotate(show_max,xy=(x[j_max],y[i_max]),xytext=(x[j_max],y[i_max]),fontsize=14)
    plt.colorbar()
    plt.xlabel('x', fontsize=20)
    plt.ylabel('y', fontsize=20)
    plt.axis('equal')
    plt.draw()
    plt.pause(0.001)
    plt.clf()

移动内热源的实现

---

Neumann边界条件

可以看到热源向右匀速移动的效果，但是左右边界上似乎有点奇怪，边界没有受到热源的影响。这是因为我们给定了边界上的节点温度，即Dirichlet边界条件，边界温度与环境温度始终相同，那么如何处理热流边界条件呢

        # solve boundary nodes 
        U[:,0] = U[:,-1] = u_env 
        U[0,:] = U[-1,:] = u_env 

热流边界条件是第二类边界条件即Neumann边界条件，给定边界节点处的一阶导（斜率）

同样采用差分格式离散化，则有

通常边界给定对流换热系数，这样就可以建立热流平衡方程，使边界点内外两边的温度梯度达到平衡，即

这样边界温度值

就可以根据环境温度

和内部节点温度

的值加权平均求得，

代表导热系数，

表示边界处的对流换热系数，系数

就反映了热传导和对流的相对大小，也决定了边界节点温度更接近内部还是外部温度。假定环境温度保持不变，边界均为对流散热边界，则有

将 solve boundary nodes 部分改为如下形式

        # solve boundary nodes
        ch = 0.1 
        U[:,0] = (u_env + ch*U[:,1]/dx)/(1+ch/dx)
        U[:,-1] = (u_env + ch*U[:,-2]/dx)/(1+ch/dx) 
        U[0,:] = (u_env + ch*U[1,:]/dy)/(1+ch/dy)
        U[-1,:] = (u_env + ch*U[-2,:]/dy)/(1+ch/dy) 

对流换热边界条件，边界热量流出速度下降

如果

，

，那么边界节点温度将等于环境温度，与等温边界结果相同。如果

，

，那么将没有热量可以从边界流出，此时就形成了绝热边界条件，边界温度等于内部相邻节点温度，结果如下图

绝热边界条件，等温线垂直于左右边界

---

边界形状

如果边界并非矩形，该如何用矩阵处理呢？关键在于确定矩阵中哪些节点是边界节点，可以使用Python中的Boolean逻辑矩阵。这里构造一个半径为7的圆形边界，用plt.imshow()来显示Boolean矩阵

    x=np.linspace(D[0],D[1],U.shape[1])
    y=np.linspace(D[2],D[3],U.shape[0])
    X,Y = np.meshgrid(x,y)
    
    Mat_e = (X-7.5)**2 + (Y-7.5)**2 >= 7**2
    Mat_b = abs(((X-7.5)**2 + (Y-7.5)**2)**0.5 - 7) <= 0.5*min(dx,dy)
    plt.imshow(Mat_b, cmap=plt.cm.get_cmap('jet'), origin = 'lower',vmin=0, vmax=1)
    plt.show()

Mat_e代表外部环境，Mat_b代表边界。离散的边界点不可能完全在解析曲线上，只要节点到边界曲线（圆周）的距离小于半个单元长度，就认为是边界节点。元素值为Ture的地方显示红色，False显示蓝色，

边界节点

外部节点

求解边界时令外部节点温度等于环境温度，边界节点采用等效的热流边界，即以内部节点温度的最小值作为梯度参考，省略了计算法向的步骤。

        # solve boundary nodes
        ch = 0.1e6 
        U[Mat_e] = u_env
        U[Mat_b] = (u_env + ch*np.min(U[~(Mat_e+Mat_b)])/dx)/(1+ch/dx)
    

可以看到变成了圆形边界，但温差太小，分界不是很明显。为了使结果表达更清晰，可以让外部节点不显示呢，Python中可以通过将对应位置的点赋值为np.nan实现

#在函数showcontourf增加下列语句，且多传入一个参数M_e
def showcontourf(mat,D,M_e,cmap=plt.cm.get_cmap('jet'),fsize=(12,12),vmin=0, vmax=100):
    ...
    ...
    Z = mat.copy()
    Z[M_e] = np.nan 
    plt.contourf(X, Y, Z, 100, cmap = cmap, origin = 'lower', levels = levels)

#在main函数中,传入外节点Boolean矩阵
    showcontourf(U,D,Mat_e,vmax=100)

注意这里需要将mat复制一份（mat.copy()）再进行赋值操作，否则将会改变原数组mat，后续计算无法进行。显示结果如下

不显示外部环境节点

---

周期边界条件

假设矩形计算区域代表的不是一个平面，而是一个首位相接的圆环外壁，那么这时候会出现一个十分有趣的周期边界。最左边的节点和最右边的节点实际上是相同的节点，这个周期边界该如何实现呢？代码上十分简单，只需要将二阶偏导系数矩阵的左下角和右上角的值置为1就可以了（以x方向为例），这样左右边界点的二阶导就可以由完整三项差分求得。

    # x方向二阶导系数矩阵A
    A = (-2)*np.eye(Mx+1,k=0) + (1)*np.eye(Mx+1,k=-1) + (1)*np.eye(Mx+1,k=1)
    A[0,-1] = A[-1,0] = 1

    # solve boundary nodes 去掉X方向的热流边界约束
    ch = 0.1e2
    U[0,:] = (u_env + ch*U[1,:]/dy)/(1+ch/dy)
    U[-1,:] = (u_env + ch*U[-2,:]/dy)/(1+ch/dy) 

周期边界，X方向连续

可以看到x方向的左右边界是连续的、周期变化的，类似于圆环上的圆周切线方向。