坐标变换法计算不规则图形导热
以lua脚本为例
一、不规则图形导热问题
《数值传热学》课本上大多讨论正交坐标系下的导热问题,网格划分和控制方程的求解相对简单。而对于不规则的几何体内部的导热问题,由于无法直接建立正交坐标系并进行网格划分,即使划分了网格也无法适应几何结构,从而降低了求解精度。
以如下简单几何体的导热问题为例,介绍坐标变换法。下图为一个梯形几何结构,顶边长1米,底边长2米,高1米,定边、底边及左侧边均为第一类边界条件,恒温10℃,右侧斜边给定一个第三类边界条件,与温度40℃的环境换热,换热系数为30W/m2。
针对以上几何结构,常规算法也可以实现导热的数值计算,网格划分如下:
正交网格下,右侧斜边部分仅占半个网格。当然也好处理,按照半个网格来求解热平衡就行,只是编程处理变复杂了,也无法适应右侧边界。
二、坐标变换法
坐标变换法就是将不规则的图形通过一定的坐标变换,变换为规则的几何图形,而在新的坐标系下网格划分得到简化。对于以上图形,网格划分如下:
假定原坐标系为(x,y),变换后为(
η
\eta
η,
ξ
\xi
ξ),变换关系如下:
ξ
=
x
(
2
−
y
)
η
=
y
\xi=x(2-y) \\ \eta=y
ξ=x(2−y)η=y
新坐标系下的网格如下:
正交的坐标系下,网格划分及离散相对简单,但是,右侧边界条件及控制方程需要调整。
1.一阶导数
原边界条件如下:
λ
∂
T
∂
n
=
h
a
(
T
a
−
T
)
=
λ
2
(
∂
T
∂
x
+
∂
T
∂
y
)
\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} )
λ∂n∂T=ha(Ta−T)=2λ(∂x∂T+∂y∂T)
新坐标系下:
d
T
=
∂
T
∂
x
d
x
+
∂
T
∂
y
d
y
(
∗
1
)
d
T
=
∂
T
∂
ξ
d
ξ
+
∂
T
∂
η
d
η
(
∗
2
)
d
ξ
=
(
2
−
y
)
d
x
−
x
d
y
(
∗
3
)
d
η
=
d
y
(
∗
4
)
dT=\frac{\partial T}{\partial x} dx+\frac{\partial T}{\partial y} dy (*1)\\ dT=\frac{\partial T}{\partial \xi} d\xi+\frac{\partial T}{\partial \eta} d\eta (*2)\\ d\xi =(2-y)dx-xdy (*3)\\ d\eta=dy (*4)
dT=∂x∂Tdx+∂y∂Tdy(∗1)dT=∂ξ∂Tdξ+∂η∂Tdη(∗2)dξ=(2−y)dx−xdy(∗3)dη=dy(∗4)
将上式(*3)(*4)反代回(*2),对比(*1)可得:
∂
T
∂
x
=
∂
T
∂
ξ
(
2
−
y
)
=
T
ξ
(
2
−
η
)
∂
T
∂
y
=
∂
T
∂
η
−
∂
T
∂
ξ
x
=
T
η
−
ξ
2
−
η
T
ξ
\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
∂x∂T=∂ξ∂T(2−y)=Tξ(2−η)∂y∂T=∂η∂T−∂ξ∂Tx=Tη−2−ηξTξ
矩阵
J
=
[
2
−
y
0
−
x
1
]
=
[
2
−
η
0
−
ξ
2
−
η
1
]
J= \begin{bmatrix} 2-y & 0 \\ -x & 1 \end{bmatrix}= \begin{bmatrix} 2-\eta & 0 \\ -\frac{\xi}{2-\eta} & 1 \end{bmatrix}
J=[2−y−x01]=[2−η−2−ηξ01]
为该变换的雅可比矩阵。总结如下:
V
o
r
g
′
=
J
V
n
e
w
′
(
∗
4
)
V_{org}'=JV_{new}' (*4)
Vorg′=JVnew′(∗4)
其中V为该问题中任意变量。
基于此,新坐标系下的边界条件变为:
h
a
(
T
a
−
T
)
=
λ
2
(
∂
T
∂
x
+
∂
T
∂
y
)
=
λ
2
[
∂
T
∂
ξ
(
2
−
x
−
y
)
+
∂
T
∂
η
]
=
λ
2
[
∂
T
∂
ξ
(
2
−
ξ
2
−
η
−
η
)
+
∂
T
∂
η
]
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} ]
ha(Ta−T)=2λ(∂x∂T+∂y∂T)=2λ[∂ξ∂T(2−x−y)+∂η∂T]=2λ[∂ξ∂T(2−2−ηξ−η)+∂η∂T]
2.二阶导数
原坐标下,该问题的控制方程如下:
∂
2
T
∂
x
2
+
∂
2
T
∂
y
2
=
1
α
∂
T
∂
τ
\frac{\partial^2 T}{\partial x^2}+\frac{\partial^2 T}{\partial y^2}=\frac{1}{\alpha} \frac{\partial T}{\partial \tau}
∂x2∂2T+∂y2∂2T=α1∂τ∂T
需要在新的坐标系下获得改问题的控制方程,为此需要获得温度关于新坐标系的二阶导数。
由(*4)可知:
(
T
x
)
o
r
g
′
=
J
(
T
x
)
n
e
w
′
(
T
y
)
o
r
g
′
=
J
(
T
y
)
n
e
w
′
(T_x)_{org}'=J(T_x)_{new}'\\ (T_y)_{org}'=J(T_y)_{new}'
(Tx)org′=J(Tx)new′(Ty)org′=J(Ty)new′
其中:
T
x
=
∂
T
∂
x
,
T
y
=
∂
T
∂
y
T_x=\frac{\partial T}{\partial x}, T_y=\frac{\partial T}{\partial y}
Tx=∂x∂T,Ty=∂y∂T
所以,通过化简可以得到二阶导数:
∂
2
T
∂
x
2
=
(
T
x
)
x
=
(
T
x
)
ξ
(
2
−
η
)
=
(
T
ξ
(
2
−
η
)
)
ξ
(
2
−
η
)
=
(
2
−
η
)
2
T
ξ
ξ
∂
2
T
∂
y
2
=
(
T
y
)
y
=
(
T
y
)
η
−
ξ
2
−
η
(
T
y
)
ξ
=
(
T
η
−
ξ
2
−
η
T
ξ
)
η
−
ξ
2
−
η
(
T
η
−
ξ
2
−
η
T
ξ
)
ξ
=
T
η
η
−
2
ξ
2
−
η
T
ξ
η
+
ξ
2
(
2
−
η
)
2
T
ξ
ξ
\frac{\partial^2 T}{\partial x^2}=(T_x)_x=(T_x)_\xi(2-\eta)=(T_\xi(2-\eta))_\xi(2-\eta)=(2-\eta)^2T_{\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}
∂x2∂2T=(Tx)x=(Tx)ξ(2−η)=(Tξ(2−η))ξ(2−η)=(2−η)2Tξξ∂y2∂2T=(Ty)y=(Ty)η−2−ηξ(Ty)ξ=(Tη−2−ηξTξ)η−2−ηξ(Tη−2−ηξTξ)ξ=Tηη−2−η2ξTξη+(2−η)2ξ2Tξξ
从而,控制方程简化为:
(
(
2
−
η
)
2
+
ξ
2
(
2
−
η
)
2
)
T
ξ
ξ
−
2
ξ
2
−
η
T
ξ
η
+
T
η
η
=
1
α
T
τ
((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
((2−η)2+(2−η)2ξ2)Tξξ−2−η2ξTξη+Tηη=α1Tτ
三、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循环,预测向量化后程序计算速度将得到进一步提升,目前这部分工作笔者还未进行。