传热学课件--导热问题的数值求解
第六章 导热问题的数值求解随着计算机的普及应用和性能的不断改善，以及相关的数值计算方法的发展和应用程序的开发，传热学数值计算方法作为数值求解传热问题的有效工具也得到了相应的发展，利用计算机求解传热学问题愈来愈受到人们的普遍重视，而且在计算复杂传热问题中显示出它的优越性，因而成为传热学的一个重要的分支。
数值传热的相关内容也很自然地成为工程类学生学习传热学课程的不可缺少的部分。
为了使学生能简要地掌握传热学数值计算的基本方法，在这里我们以导热问题为例对传热学数值计算方法做一个简单的介绍。
6-1导热问题数值解概述在第二章和第三章中我们对较为简单的导热问题，如一维、二维简单几何形状和边界条件的稳态导热和非稳态导热、以及一些特殊导热问题，象通过肋片的导热和忽略内阻的集总导热系统，进行了分析求解。
然而对于一些更为复杂的导热问题，如复杂的几何形状和边界条件以及物性变化较大的情况，分析求解往往很复杂或者根本不可能。
此时求解问题的唯一途径是利用数值分析的办法获得数值解。
数值求解通常是对微分方程直接进行数值积分或者把微分方程转化为一组代数方程组再求解。
这里要介绍的是后一种方法。
如何实现从微分方程到代数方程的转化又可以采用不同的数学方法，如有限差分法、有限元法和边界元法等。
作为一本入门的教材，这里仅向读者简要地介绍用有限差分析方法从微分方程确立代数方程的处理过程。
有限差分法的基本思想是把原来在时间和空间坐标中连续变化的物理量（如温度、压力、速度和热流等），用有限个离散点上的数值集合来近似表示。
有限差分的数学基础是用差商代替微商（导数），而几何意义是用函数在某区域内的平均变化率代替函数的真实变化率。
在图6-1中可以看出有限差分表示的温度场与真实温度场的区别。
图中用T 0、T 1、T 2…表示连续的温度场T ；Δx 为步长，它将区域的x 方向划分为有限个数的区域，Δx 0、Δx 1、Δx 2…，它们可以相等，也可以不相等。
当Δx 相等时，T 1处的真实变化率a 可以用平均变化率b 、c 或d 来表示，其中b 、c 和d 分别表示三种不同差分格式下的温度随时间的变化率，即：T 图4-1温度场的有限差分表示b 为向后差分格式xx x T x T dx dT ∆∆--≈)()(111； c 为向前差分格式xx T x T dx dT ∆-∆≈)()(111＋；d 为中心差分格式xx x T x x T dx dT ∆∆--∆+≈)()(111。
这种差分格式也可以推广到高阶微商的情形。
对于二阶微商的差分格式可以在一阶差分格式的基础上得出：2111212)()()(2)(x x x T x T x x T dx T d ∆∆-+-∆+≈。
采用这样的处理之后，反映温度场随时间、空间连续变化的微分方程就可以用反映离散点间温度线性变化规律的代数方程来表示。
当利用相应的数学办法求解这些代数方程组之后，我们就能获得离散点上的温度值。
这些温度值就可以近似表示温度场的连续的温度分布。
从上面的分析不难看出，当我们要对导热问题进行数值求解时一定要采取三个大的步骤，即研究区域的离散化；离散点（节点）差分方程的建立；节点方程（代数方程）的求解。
下面我们将导热问题的数值求解进行较为详细的讨论。
6-2研究区域温度场的离散化导热问题的温度场是假设为时间和空间的连续函数，当进行数值求解时首先要做的事情是在所研究的时间和空间区域内把时间和空间分割成为有限大小的小区域，尤如地球被人为地划分为不同的地域且冠以不同的名称，时间被年、月、日和时、分、秒分割。
如果在所分割的每一个时间间隔和空间区域内均用同一个温度值来表示，那么原来连续变化的温度场就被一个离散的阶跃变化的温度分布所代替。
这就是连续变化的温度场离散化处理的基本思路。
这里我们以一个矩形长柱体的非稳态导热过程为例来讨论区域离散化问题。
如果不考虑矩形长柱体长度方向上的温度变化，那么它是一个二维非稳态导热问题，图6-2表示了长柱体矩形截面上区域离散化的情况。
图中可见，对于给定的空间区域，在x 方向上的步长为Δx ，在y 方向上的步长为Δy ，用它们作为空间尺度可以将矩形区域划分成纵横交错的网格，交点称为节点。
然后以节点为中心，在两个节点的中心处划分界限，定出节点的控制面积，对于三维情况则为控制体积或控制容积，因而常在一般意义上称之为节点的控制体。
控制体的形状是随着坐标系的不同而改变的，这里的控制体是一个个的矩形面积。
网格的步长在每一个方向上可以均匀划分，也可以不均匀的划分。
因此，x图4-2矩形长柱体截面区域离散化选用不同的步长和不同的划分方法，可以将同一区域划分出不同大小、不同数目的控制区域，以及不同数目的节点数。
显然，随着步长的不断减小，节点数目的不断增加，由节点温度表示的离散的温度场就会更加接近连续的温度场，但计算工作量也会随之增加。
在时间方向上离散化的步长常用Δτ来表示，Δτ的选取也是可大可小的，也可以随时间的进程而变化。
显然，无限小的时间步长Δτ亦会使得离散温度变化接近连续的温度改变，但随之而来的是相应的计算工作量的增加。
6-3温度场节点方程的建立为了得出所研究区域的节点温度，必须建立相应的节点方程。
建立节点方程可以采用不同的方法，为了更好地理解节点方程的物理意义和掌握节点方程的建立方法，我们采用控制体热平衡法来建立节点方程。
下面我们以实例的形式介绍不同节点的节点方程的建立过程。
1 控制体的内节点方程控制体热平衡法建立节点方程的过程是将能量守恒方程应用于控制体，建立该节点与周围节点之间的能量平衡关系式，再利用傅立叶的导热定律，最后获得控制体节点温度与周围节点温度之间的关系式。
考察图6-2中的节点P 及其控制体，由能量平衡关系应有E Q Q Q Q Q V N S E W ∆=＋＋＋＋，6-1式中，Q W 、Q E 、Q S 和Q N 分别为邻近节点W 、E 、S 和N 通过传导方式传给节点P 的热流量；Q V 为单位时间控制体内热源的发热量；ΔΕ为控制体单位时间内热能的增加量。
由导热傅立叶定律，在线性温度分布的假设下，时刻K 周围节点传给节点P 的热流量分别为：，1.)(;1.)(;1.)(;1.)(x T T yQ x T T yQ y T T xQ y T T xQ K P K N N K P K S S K P K E E K P KW W ∆-∆=∆-∆=∆-∆=∆-∆=λλλλ以及控制体的发热流量1⋅∆∆=y x q Q V V ,（q V 为内热源强度，即单位时间单位体积的内热源发热量。
）控制体单位时间的内能增加量为1111⋅∆∆∆-=∆⋅∆∆∆-=∆y x T T c E y x T T c E K P K P K P K P τρτρ－＋或，前者为时间上的向前差分，而后者为时间上向后差分。
以上关系式中温度T 的上标为所在时刻，下标为所在空间位置。
将以上关系式一并代入方程6-1中，且假设Δx=Δy ，经整理可以得出二维非稳态导热问题的内节点的两种差分格式的差分方程，即显示差分格式c q T xa T T T T x a T V K P KN K S K E K W K P ρτττ∆+∆∆-++++∆∆=)41()(221＋ 6-2和隐示差分格式])([11122c q T T T T T xa xa T VK P K N K S K E K W K P ρτττ∆+++++∆∆∆∆=-＋。
6-3比较上面两种差分格式可以看出，显示差分格式最突出的优点是节点温度表达式的右边只涉及K 时刻的节点温度值，那么只要知道这一时刻周围节点的温度值就可以求出该节点的下一时刻的温度值；而隐示差分格式却相反，温度表达式的两端都是同一K 时刻的节点温度值，这就意味着必须同时计算同一时刻所有节点的温度值，即必须联立求解K 时刻所有节点的差分方程组，增大计算工作量是显而易见的。
虽然显示差分格式计算比较方便，但它却存在着一个缺点，即计算式中a Δτ/Δx 2值必须满足一定的条件才不至于引起数值计算出现不收敛的问题，这在数值计算中称为差分格式的不稳定性。
这里差分方程稳定性的条件是方程6-2中的变量T 前面的系数必须大于或等于零，分析一下差分方程中的各项系数，在T P K前的系数应为0412≥∆∆x a τ－，改写成为 4142≤∆∆x a τ。
6-4此式称为显示差分格式的稳定性判据，从中看出时间步长和空间步长是相互制约的。
为了获得较为精确的节点温度值，空间步长Δx 的选择不能太小，按照稳定性判据的要求势必会使时间步长Δτ也要相应地不能太大，因而必须在增加节点数目的同时增多时间间隔，从而使计算工作量加大。
与显示差分格式相反，由于隐示差分格式的节点方程中没有会使方程系数成为负值的系数项，因而不存在方程求解的不稳定性的问题。
也就是说，对于隐示差分格式，无论a Δτ/Δx 2中的Δx 和Δτ取什么样的数值，均不会出现数值计算结果的不收敛问题，因而是无条件稳定的。
这样就使得我们能在满足一定精度的情况下尽可能地加大时间步长或空间步长，亦可以在计算过程中随意改变步长，而不必担心会造成计算结果的不收敛。
2 控制体的边界节点方程在数值计算中所研究的区域的边界条件是通过边界节点的节点方程来反映的，因而边界节点的差分方程的建立十分重要。
这里同样采用边界节点的控制体热平衡来确立边界节点的差分方程。
下面将具体讨论一些典型的边界节点的节点方程。
对流换热边界条件下的节点方程。
如图6-3所示的边界节点P ，其控制体的热平衡关系式为E Q Q Q Q Q C V N S W ∆=++++，6-5式中各项可以写成：（隐示格式）。
（显示格式）；1212;1)(;12;12;12;111⋅∆∆∆-=∆⋅∆∆∆-=∆⋅∆-=⋅∆∆=⋅∆∆-=⋅∆∆-=⋅∆∆-=-+∞y xT T cE y xT T cE y T T Q y xq Q xy T T Q xy T T Q y xT T Q K P K P KP K P K P K c V V K P K N N K P K S S KP K W W τρτραλλλ将上述关系式一并代入6-5式中，经整理可得出在二维非稳态导热过程中处于平直边界上的边界节点在对流换热边界条件下的节点差分方程（当Δx=Δy ），即：对于显示差分格式有c q T x a x xa T x a x T T T x a T V K P K KNK S K W K P ρττλαττλατ∆+∆∆∆-∆∆-+∆∆∆+++∆∆=∞)241(2)2(22221＋； 6-6对于隐示差分格式有]2)2([21112222c q T T xa x T T T x a x a x xa T V K P K KNK S K W K P ρττλαττλατ∆++∆∆∆+++∆∆∆∆∆+∆∆=-∞＋。
内科大传热学课件04导热问题的数值解法
x t x
m,n
x2 2
2t x2
m,n
x3 6
3t x3
m,n
x4 24
4t x4
用节点(i,j)的温度ti,j来表示节点(i-1,j)的
温度ti-1,j
t
x2 2t
x3 3t
x4 4t
tm1,n tm,n x x m,n
2
x2 m,n
6
x3 24 x4 m,n
将上两式相加可得
tm1,n
对元体 (m,n). 根据能量守恒定律可知：
e w n s 0
其中，规定：导入元体（ m,n ）的热流量 为正；导出元体（ m,n ）的热流量为负。
说明： ① 上述分析与推导是在笛卡儿坐标系中进行 的； ② 热平衡法概念清晰，过程简捷； ③ 热平衡法与建立微分方程的思路与过程一 致，但不同的是前者是有限大小的元体，后 者是微元体。
的最低阶数为2
根据导热问题的控制方程 ( 导热微分方程 )
2t x2
2t y 2
0
得
tm1,n
2tm,n x2
tm1,n
tm,n1
2tm,n y 2
tm,n1
0
若 △x=△y 则有
1 tm,n 4 (tm1,n tm1,n tm,n1 tm,n1)
(2) 控制容积平衡法(热平衡法)
h1t f
x
（a）
3 基本概念：控制容积、网格线、节点、
界面线、步长
(m,n) N
二维矩
形域内
稳态无
n
内热源，
常物性
y
的导热
y
问题
x x
M m
ቤተ መጻሕፍቲ ባይዱ
传热学60-第四章 导热问题的数值解法
B (i,j 1)
第四章 导热问题的数值解法 9
根据傅里叶定律， L,R,T ,B各节点向P节点的导热量：
T (i,j 1)
LP RP BP LP
t i 1 , j t i , j ti 1 , j ti , j t i , j 1 t i , j ti , j 1 ti , j
tik 1 tik
的一阶导数采用向前差分，则
tik 1 2tik tik 1 a x 2
第四章 导热问题的数值解法 28
上式移项整理 k k 1 k k ti a ( t i 1 t i 1 ) ( 1 2 a )ti 2 2 x x
x
y 1 y 1 x 1 x 1
(i 1,j )
(i 1,j )
L
R
x
Z方向取单位长度 (i,j 1) B
y
y
第四章 导热问题的数值解法
10
若有内热源， v ,i , j 为P节点所在网格单元的内热源强度
则内热源发热量
在稳态导热下：
v , p v ,i , j x y 1
ti 1, j ti , j y x y 3 1 h(t f ti , j ) xy v ,i , j 0 x 2 2 4
第四章 导热问题的数值解法
19
不规则区域的处理
用阶梯形的折线来模拟真实边界
t2
t1
t2
t2
t2
t2
第四章 导热问题的数值解法 20
第四章 导热问题的数值解法
17
3.
内部角点
内部角点
P(i ,j ) ti ,j
传热学-学习课件-4-1 导热问题数值求解基本思想
有限差分法（FDM），主要介绍方法 有限元法（FEM） 边界元法（BEM）
传热学 Heat Transfer
三、数值解法的基本步骤
传热学 Heat Transfer
以一个二维稳态导热问题为例，介绍数值求解导热问题的具体 过程，重点是节点离散方程的建立和代数方程组的迭代求解。
（1）物理问题 二维矩形域稳态无内热源，常物性的导热问题 y
传热学 Heat Transfer
4-1 导热问题数值求解的基本思想
一、数值解法的本质
数值解法是用物理问题所 涉及的空间和时间区域内有 限个离散点（称为节点）的 物理量近似值来代替物体内 实际连续的物理量分布，将 连续物理量分布函数的求解 问题转化为各节点物理量值 的求解问题。
传热学 Heat Transfer
传热学 Heat Transfer
传热学 Heat Transfer
主讲老师：王舫 适用专业：能源与动力工程专业
传热学 Heat Transfer
第四章 热传导问题的数值解法
§4-0 引言
1 求解导热问题的三种基本方法： (1) 理论分析 (2)实验 (3)数值计算
2 三种方法的特点 3 三种方法的基本求解过程
传热学heattransfer?本章教学内容41导热问题数值求解基本思想42内节点离散方程的建立43边界节点离散方程的建立及代数方程的求解44非稳态导热问题的数值解法传热学heattransfer41导热问题数值求解的基本思想一数值解法的本质数值解法是用物理问题所涉及的空间和时间区域内有限个离散点称为节点的物理量近似值来代替物体内实际连续的物理量分布将连续物理量分布函数的求解问题转化为各节点物理量值的求解问题
传热学 Heat Transfer
四章导热问题的数值解法-PPT文档资料
fi 2 fi1 h2
fi 2
二 阶 导 数 向 后 差 分 ： f i
fi2 2 fi1 h2
fi
二 阶 导 数 中 心 差 分 ： f i
fi1 2 fi h2
fi1
式中：
f i
df dx
i
fi
d2 f dx 2
i
图4-2 有限差分表达 式的几何意义
• 解此代数方程组，得到节点上温度的近似值
2、函数 f（x）在点 x 的导数的有限差分表达式：
函数f (x)在点 x 0 的泰勒级数展开形式为：
f( x ) f( x 0 ) ( x x 0 )f( x 0 ) ( x 2 ! x 0 ) 2f( x 0 ) ( x 3 ! x 0 ) 3 f( x 0 )
导数向前、向后及中心差分公式为：
、
一阶导数向前差分：
f(x)f(xh)f(x)
h
一阶导数向后差分： f(x)f(x)f(xh) h
一阶导数中心差分：
f(x)f(xh)f(xh) 2h
二阶导数向前差分： f(x)f(x)f(xh 22 h)2f(xh)
二阶导数向后差分： f(x)f(x)f(xh 22 h)2f(xh) 二阶导数中心差分： f(x)f(xh)f(hx2h)2f(x)
由式(b)和式(d)消去f (x) 得：
f(x)f(x )f(x 2 h ) 2f(x h ) O (h 2) h 2
(i)
由式(a)和式(b)消去f (x) 得： f(x )f(x h )f( h x 2 h ) 2f(x ) O (h 3 ) (j)
由(e)式～(j)式分别略去 h 、h 2 及 h 3 以上各项得一阶、二阶
传热学第五章导热问题数值解法
2 h∆y
=
h∆y 2 2 + λ t m ,n+1 + t m ,n−1 + 2(t m−1,n + Bi ⋅ t f ) 2(2 + Bi )
9
λ
tf
3）方程组的求解 ） 对应每个未知量（一个节点温度） 一条方程 一条方程（ 对应每个未知量（一个节点温度）→一条方程（一 个节点方程） 方程组有唯一解 个节点方程）→方程组有唯一解 解法：（ ：（1） 解法：（ ）矩阵法 （2）迭代法：高斯 赛德尔迭代 ）迭代法：高斯-赛德尔迭代
λ∆y
t m−1,n − t m ,n ∆x ∆x t m ,n−1 − t m ,n +λ + ∆y ⋅ h(t f − t m ,n ) = 0 2 ∆y ∆y ∆x t m ,n+1 − t m ,n +λ ∆y 2
如取正方形网络
∆x = ∆y
上式简化为： 上式简化为：
t m ,n =
2t m−1,n + t m ,n+1 + t m ,n−1 +
对于非稳态导热问题，除了空间上进行网格划分外， 对于非稳态导热问题，除了空间上进行网格划分外， 还要把时间分割成许多间隔。 还要把时间分割成许多间隔。
4
2）有限元法 ）
把整个求解域离散成为有限个子域， 把整个求解域离散成为有限个子域，每一子域内运 用变分法， 用变分法，即利用与原问题中微分方程相等价的变 分原理来进行推导，从而使原问题的微分方程组退 分原理来进行推导， 化到代数联立方程组，得到数值解。 化到代数联立方程组，得到数值解。 有限元法和差分法都是常用的数值计算方法， 有限元法和差分法都是常用的数值计算方法，差分 法计算模型对于不规则的几何形状难以应用。 法计算模型对于不规则的几何形状难以应用。有限 元法能够很好地适应复杂的几何形状、 元法能够很好地适应复杂的几何形状、复杂的材料 特性和复杂的边界条件。 特性和复杂的边界条件。
传热学课件第四章 导热问题数值解法基础
i , j
t x
t i 1 , j t i , j x
0 x
2.一阶导级的向后差分表达式：舍去<2>式△x2后各项，则有：
i , j
t x
t i , j t i 1 , j x
0 x
第一节 建立离散方程的方法
二、泰勒级数展开法（有限差分法）
k 2 k 1
对 流 h t f t1 A
k k
显式
△x
C.内能增量△u：
u c
x 2
A t1
k
k 1
t1 /
k
△x/2
k hx
据热平衡A+B=C并整理得：
k f
t 2 t1
k
t
t1
k
1 2
c
x
2
t1
k 1
LP
△y
t i 1 , j t i , j x
t i , j 1 t i , j y
y 2
x 2
1
BP
1
x 2
y 2
EP h t f t i , j
△x
1
FP h t f t i , j
t x
t
2
2
x i , j 2!
2
t x
3
x i , j 3!
3
3.一阶导级的中心差分表达式：<1>-<2>式且忽略后项，则有：
i , j
t x
热传导问题的数值解法
1. 空间离散化
01
将求解区域划分为一系列小的网格或节点，用离散的差分代替
微分。
2. 时间离散化
02
将时间轴划分为一系列小的时段，用离散的差分代替微分。
3. 初始条件和边界条件的离散化
03
将初始条件和边界条件转化为离散形式。
差分方程的求解
01
1. 迭代法
2. 直接法
02
03
3. 松弛法
通过迭代逐步逼近解，常用的有 Jacobi迭代法和Gauss-Seidel迭 代法。
02
根据问题的几何特性和求解精度要求，选择合适的单元类型和
划分方式。
单元划分应尽量保证求解精度和计算效率。
03
建立系统方程
01
根据热传导的物理定律和边界条件，建立每个单元的热平衡 方程。
02
将各个单元的方程联立起来，形成整个求解域上的系统方程 。
03
系统方程通常为线性方程组，可以使用不同的求解方法进行 求解。
步骤
首先将求解区域划分为一系列离散点，然后根据泰勒级数展开，将偏微分方程 中的导数项用离散点上的函数值之差代替，从而得到离散化的差分方程。
特点
有限差分法简单直观，适用于规则区域，但对不规则区域 处理较为复杂。
有限元法
定义
有限元法是一种将连续的求解区域离散化为有限个小的子域（即有限元），然后在每个子 域上应用数学方法进行求解的方法。
热传导定律也称为傅里叶定律，指出热流密度与温度梯度成正比，方向由高温指向低温。数学表达式为：q = -k * grad(T)，其中 q为热流密度，k为导热系数，T为温度，grad表示梯度。
热传导定律是热传导过程的基本规律，描述了热量传递的方向和大小，是数值解法的基础。
传热学第四章导热问题的数值解法
三.灰体表面间的辐射换热
由于灰体表面存在着多次的反射和吸收，计算起来比 黑体复杂得多，为了计算方便我们引进了两个概念： 1.投入辐射G和有效辐射J G—单位时间内投射到表面单位面积上的总辐射能。 J—单位时间内离开该表面单位面积的总辐射能。
G(投入辐射) (1-α)G(反射辐射) αG (吸收辐射) J(有效辐射) E= ε Eb (本身辐射)
定义：表面1发出的辐射能直接落到表面2上的百分数， 称表面1对表面2的角系数，记为X1,2，同理有X2,1。 角系数是一个无量纲能量百分比。引入角系数是为了说 明两个表面之间的辐射换热量与它们之间的相对位置有 很大关系。 角系数为几何因子，其值取决于物体的几何特性（形状、 尺寸及物体的相对位置）而与物体的种类和温度无关。
i =1 n
表面为凸面或平面时，有何性质？
3.可加性
从表面1上发出落到表面2上的总能量，等于落到表 面2各部分的能量之和。于是有：
X 1， 2 = X 1， 2 a + + X 1， 2 n = ∑ X 1， 2i
i =a n
设表面由a、b两部分组成，写出其可加性表 达式。
三.角系数的计算
或写成：
Eb − J q = ， 1− ε
ε
1− ε 称为表面辐射热阻 εA
2.灰体表面间的辐射换热
考虑两个等温的漫灰表面组成的二维封闭系统，表 面1、2间的Φ为： φ1, 2 = A1 J 1 X 1, 2 − A2 J 2 X 2,1
J 1 A1 = A1 Eb1 − (1 / ε 1 − 1)φ1, 2 J 2 A2 = A2 Eb 2 − (1 / ε 2 − 1)φ 2,1
本节讨论的是被透热介质隔开的两固体表面间的辐射换 热。透热介质是指不参与热辐射的介质，最常见的是空 气。
传热学-第4章-热传导问题的数值解珐
若步长∆x=∆y，有： ， 若步长
t m ,n = 1 ( 2 t m −1 , n + t m , n + 1 + t m , n −1 + 4 ∆2 x Φ m , n
λ
+
2 ∆ xq w
λ
)
2. 外部角点 控制容积的热平衡为： 控制容积的热平衡为：
∆y tm−1,n − tm,n ∆x tm,n−1 − tm,n ∆x∆y ∆x + ∆y λ +λ + Φ m, n + qw = 0 ∆x 2 2 ∆y 4 2
4. 边界热流密度的三种情况
q （1）绝热边界： w = 0 ）绝热边界：
（2） qw 值不为零：代入给定的 qw 值。 ） 值不为零： （3）对流边界：qw = h(t f ）对流边界： 平直边界节点： 平直边界节点：
2( h∆x
− t m n = 2 t m − 1 , n + t m , n + 1 + t m , n −1 +
第一类边界条件 — 边界温度已知 m-1,n 第二类边界条件 需建立边界节点温度 ∆y 第三类边界条件 的差分方程 n 1. 位于平直边界上的节点
λ∆y
tm−1,n − tm,n ∆x +λ
m m,n+1
qw
m,n m,n-1
∆x
∆x tm,n+1 − tm,n ∆x tm,n−1 − tm,n ∆x∆y +λ + Φm,n + ∆yqw = 0 2 ∆y 2 ∆y 2
若步长∆x=∆y，有： ， 若步长
t m ,n = 1 ( t m −1 , n + t m , n −1 + 2