坐标变换法计算不规则图形导热

坐标变换法计算不规则图形导热

以lua脚本为例

---

一、不规则图形导热问题

《数值传热学》课本上大多讨论正交坐标系下的导热问题，网格划分和控制方程的求解相对简单。而对于不规则的几何体内部的导热问题，由于无法直接建立正交坐标系并进行网格划分，即使划分了网格也无法适应几何结构，从而降低了求解精度。

以如下简单几何体的导热问题为例，介绍坐标变换法。下图为一个梯形几何结构，顶边长1米，底边长2米，高1米，定边、底边及左侧边均为第一类边界条件，恒温10℃，右侧斜边给定一个第三类边界条件，与温度40℃的环境换热，换热系数为30W/m²。

针对以上几何结构，常规算法也可以实现导热的数值计算，网格划分如下：

正交网格下，右侧斜边部分仅占半个网格。当然也好处理，按照半个网格来求解热平衡就行，只是编程处理变复杂了，也无法适应右侧边界。

---

二、坐标变换法

坐标变换法就是将不规则的图形通过一定的坐标变换，变换为规则的几何图形，而在新的坐标系下网格划分得到简化。对于以上图形，网格划分如下：

假定原坐标系为（x,y），变换后为（$\eta$,$\xi$)，变换关系如下：
$$\begin{gathered} \xi =x(2-y) \\ \eta =y \end{gathered}$$
新坐标系下的网格如下：

正交的坐标系下，网格划分及离散相对简单，但是，右侧边界条件及控制方程需要调整。

1.一阶导数

原边界条件如下：
$$\lambda \frac{\partial T}{\partial n}=h_a(T_a-T)=\frac{\lambda }{\sqrt{2}}(\frac{\partial T}{\partial x}+\frac{\partial T}{\partial y})$$
新坐标系下：
$$\begin{gathered} dT=\frac{\partial T}{\partial x}dx+\frac{\partial T}{\partial y}dy(\ast 1) \\ dT=\frac{\partial T}{\partial \xi }d\xi +\frac{\partial T}{\partial \eta }d\eta (\ast 2) \\ d\xi =(2-y)dx-xdy(\ast 3) \\ d\eta =dy(\ast 4) \end{gathered}$$
将上式（*3）（*4）反代回（*2），对比（*1）可得：
$$\begin{gathered} \frac{\partial T}{\partial x}=\frac{\partial T}{\partial \xi }(2-y)=T_{\xi }(2-\eta ) \\ \frac{\partial T}{\partial y}=\frac{\partial T}{\partial \eta }-\frac{\partial T}{\partial \xi }x=T_{\eta }-\frac{\xi }{2-\eta }{T}_{\xi } \end{gathered}$$
矩阵
$$J=\left[\begin{array}{cc}2-y& 0\\ -x& 1\end{array}\right]=\left[\begin{array}{cc}2-\eta & 0\\ -\frac{\xi }{2-\eta }& 1\end{array}\right]$$
为该变换的雅可比矩阵。总结如下：
$$V_{org}^{\prime }=J{V}_{new}^{\prime }（\ast 4）$$
其中V为该问题中任意变量。
基于此，新坐标系下的边界条件变为：
$$\begin{gathered} h_a(T_a-T)=\frac{\lambda }{\sqrt{2}}(\frac{\partial T}{\partial x}+\frac{\partial T}{\partial y})=\frac{\lambda }{\sqrt{2}}[\frac{\partial T}{\partial \xi }(2-x-y)+\frac{\partial T}{\partial \eta }] \\ =\frac{\lambda }{\sqrt{2}}[\frac{\partial T}{\partial \xi }(2-\frac{\xi }{2-\eta }-\eta )+\frac{\partial T}{\partial \eta }] \end{gathered}$$

2.二阶导数

原坐标下，该问题的控制方程如下：
$$\frac{{\partial }^{2}T}{\partial x^2}+\frac{{\partial }^{2}T}{\partial y^2}=\frac{1}{\alpha }\frac{\partial T}{\partial \tau }$$
需要在新的坐标系下获得改问题的控制方程，为此需要获得温度关于新坐标系的二阶导数。

由（*4）可知：
$$\begin{gathered} (T_x{)}_{org}^{\prime }=J(T_x{)}_{new}^{\prime } \\ (T_y{)}_{org}^{\prime }=J(T_y{)}_{new}^{\prime } \end{gathered}$$
其中：
$$T_x=\frac{\partial T}{\partial x},T_y=\frac{\partial T}{\partial y}$$
所以，通过化简可以得到二阶导数：
$$\begin{gathered} \frac{{\partial }^{2}T}{\partial x^2}=(T_x{)}_x=(T_x{)}_{\xi }(2-\eta )=(T_{\xi }(2-\eta ){)}_{\xi }(2-\eta )=(2-\eta {)}^2{T}_{\xi \xi } \\ \frac{{\partial }^{2}T}{\partial y^2}=(T_y{)}_y=(T_y{)}_{\eta }-\frac{\xi }{2-\eta }(T_y{)}_{\xi }= \\ (T_{\eta }-\frac{\xi }{2-\eta }{T}_{\xi }{)}_{\eta }-\frac{\xi }{2-\eta }(T_{\eta }-\frac{\xi }{2-\eta }{T}_{\xi }{)}_{\xi }= \\ {T}_{\eta \eta }-\frac{2\xi }{2-\eta }{T}_{\xi \eta }+\frac{{\xi }^2}{(2-\eta {)}^2}{T}_{\xi \xi } \end{gathered}$$
从而，控制方程简化为：
$$((2-\eta {)}^2+\frac{{\xi }^2}{(2-\eta {)}^2})T_{\xi \xi }-\frac{2\xi }{2-\eta }{T}_{\xi \eta }+T_{\eta \eta }=\frac{1}{\alpha }{T}_{\tau }$$

---

三、lua脚本实现

function evenval(imax,jmax,num)
    local array = {}
    local i
    local j
    for i=1,imax do
        array[i] = {}
        for j=1,jmax do
            array[i][j] = num
        end
    end
    return array
end
-- Matrix
function val(imax,jmax,num)
    local array = {}
    local i
    local j
    for i=1,imax do
        array[i] = {}
        for j=1,jmax do
            array[i][j] = num[i][j]
        end
    end
    return array
end
-- Print a matrix
function printMatrix(imax,jmax,num,filename)
    local i
    local j
    local f = io.open(filename, "w")
    for i=1,imax do
        for j=1,jmax do
            f:write(num[i][j],"\t")
        end
        f:write("\n")
    end
    f:close(f)
end
-- Test Program
function Main(imax,jmax)
    local T=evenval(imax,jmax,10.)    
    --
    local x=evenval(imax,jmax,0.)
    local y=evenval(imax,jmax,0.)
    local a1=evenval(imax,jmax,0.)
    local a2=evenval(imax,jmax,0.)
    local a3
    local a4=evenval(imax,jmax,0.)
    local m=evenval(imax,jmax,0.)
    local n=evenval(imax,jmax,0.)
    --
    local dx=1./imax
    local dy=1./jmax
    --
    local i
    local j
    --
    for i=1,imax do
        for j=1,jmax do
            x[i][j]=(j-0.5)*dx
            y[i][j]=(i-0.5)*dy

            a1[i][j]=(x[i][j]^2+1)/(2-y[i][j])^2
            a2[i][j]=2*x[i][j]/(2-y[i][j])
            a3=1
            a4[i][j]=(2*x[i][j]^2+1)/(2-y[i][j])^2/x[i][j]
        end
    end
    --
    local dt=4
    local alph=2e-7
    local Ta=40
    local ha=30
    local lamd=1.5
    local dx0,dx1,dx2,dx3,dx4,dx5
    local d1,d2,d3,d4
    --
    local k
    ---
    local Tn=val(imax,jmax,T)
    ----
    --
    for k=1,100000 do --You can change the time here.
        for i=2,imax-1 do
            for j=2,jmax-1 do
                d1=(T[i][j-1]+T[i][j+1]-2*T[i][j])/dx^2
                d2=(T[i+1][j+1]-T[i-1][j+1]-T[i+1][j-1]+T[i-1][j-1])/4/dy/dx
                d3=(T[i-1][j]+T[i+1][j]-2*T[i][j])/dy^2
                d4=(T[i][j+1]-T[i][j-1])/2/dx
                Tn[i][j]=T[i][j]+dt*alph*(a1[i][j]*d1+a2[i][j]*d2+a3*d3+a4[i][j]*d4)
            end
        end
        j=jmax

        for i=2,imax-1 do
            dx0=-ha*(2-y[i][j])*(T[i][j]-Ta)/math.sqrt(2)/lamd
            dx1=(T[i][j]-T[i][j-1])/dx
            dx2=-ha*(2-y[i-1][j])*(T[i-1][j]-Ta)/math.sqrt(2)/lamd
            dx3=(T[i-1][j]-T[i-1][j-1])/dx
            dx4=-ha*(2-y[i+1][j])*(T[i+1][j]-Ta)/math.sqrt(2)/lamd
            dx5=(T[i+1][j]-T[i+1][j-1])/dx

            d1=(dx0-dx1)/dx
            d2=(dx4+dx5-dx2-dx3)/4/dy
            d3=(T[i-1][j]+T[i+1][j]-2*T[i][j])/dy^2
            d4=(dx0+dx1)/2
            Tn[i][j]=T[i][j]+dt*alph*(a1[i][j]*d1+a2[i][j]*d2+a3*d3+a4[i][j]*d4)
        end
        T=val(imax,jmax,Tn)
    end

    for i=1,imax do
        for j=1,jmax do
            m[i][j]=x[i][j]*(2-y[i][j])
            n[i][j]=y[i][j]
        end
    end
    printMatrix(imax,jmax,m,"./m.csv")
    printMatrix(imax,jmax,n,"./n.csv")
    printMatrix(imax,jmax,T,"./T.csv")
end
--[[
Only approx. 20+ s used.
]]
Main(100,100)

---

几秒内出结果，生成m.csv、n.csv、T.csv三个文件，通过python读取画图，可以得到温度场的分布情况，代码如下：

import numpy as np
import matplotlib.pyplot as plt

m=np.loadtxt('m.csv')
n=np.loadtxt('n.csv')
T=np.loadtxt('T.csv')

fig, ax = plt.subplots()

ctf=ax.contourf(m,n,T,10,cmap='jet')

ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_aspect('equal')

ax.set_xlim(0,2)
ax.set_ylim(0,1)

fig.colorbar(ctf)

plt.show()

可以得到温度分布如下：

---

扩展

笔者采用同样的方法分别编写了matlab、python、octave、scilab的程序，测试了不同编程语言在计算速度上的差距，其中matlab也可以快速得到温度分布，与luajit相差不大，python、octave、scilab则无法在有限时间内给出结果。笔者又测试了python开JIT的情况，也可以快速得到结果，但是numba配置似乎比较奇怪，总是报错。
产生以上差距的主要原因是本程序并未向量化，大量使用for循环，预测向量化后程序计算速度将得到进一步提升，目前这部分工作笔者还未进行。