5 热传导问题的数值方法
5.1一维稳态导热
一维稳态导热在直角坐标系下的控制方程可表示为：
0)(=+s dx
dT k dx d (5-1) 式中k 为导热系数，T 是温度，s 是单位容积的热产生率。
首先选定控制体和网格，如图5.1所示，并对方程（5-1）在所选定的控制体进行积分，即得：
0)()(=+-?dx s dx
dT
k dx dT k e w w e (5-2)
图5.1 控制体和网格
然后进行离散化。如果用分线段性分布来计算方程（5-2）中的微商dx
dT
，那么最终的方程为:
0)()
()()(=?+---x s x T T k x T T k w
W P w e P E e δδ (5-3)
假设源项s 在任一控制体中之值可以表示为温度的线性函数，即P P c T s s s +=，则导出的离散化方程为：
b T a T a T a W W E E P P ++= (5-4)
式中
x s b x
s a a a x k a x k a c P W E P w w
W e
e E ?=?-+=δ=
δ=
)()( (5-5) 式(5-4)就是一维稳态导热方程的离散形式，系数a E 和a W 分别代表了节点P 与E 间及W 与P 间导热阻力的倒数，它们的大小反映了节点W 和E 处的温度对P 点的影响程度。式中的k e 和k w 是控制容积中的e 和w 界面上的当量导热系数。进行计算时，物理参数值存储在节点的位置上。为了确定k e 和k w ，还需规定由节点上的物理量来计算相应界面上的量的方法。常用的方法由两种，即算术平均法与调和平均法。
1、算术平均法
假定k 与x 呈线性关系，由P 与E 点的导数系数确定e k 的公式为：
e
e
E e e P e x x k x x k k )()()()(δδ+δδ=-+ （5-6）
2、调和平均法
利用传热学的基本公式可以导出确定界面上当量导热系数的调和平均公式。控制容积中P 和E 的导热系数不相等，但界面上热流密度应该连续，则由Fourier 定律可得：
()
()
()
()
E
e
P
e
P
E E
e
e
E P
e
P
e e k x k x T T k x T T k x T T q +-+-δ+
δ-=
δ-=
δ-=
（5-7）
而
()P
e P
E e k x T T q δ-=
则
()()()E
e P
e e
e
k x k x k x +
-
+
=
δδδ （5-8）
这就是确定界面上当量导热系数的调和平均公式，它反映了串联过程热阻的迭加原则。
3、两种方法比较
设E P k k >>，网格均分时，22P
E P e k k k k ≈
+=
（5-9） 即P 和E 两点间的导热阻力为
()()P
e
P e
k x k x δ=
δ22
/，表明此时P 和E 间的热阻主要是由导热系数大的物体所决定，这显然不符合传热学的基本原理。
实际上，此时控制体E 构成了热阻的主要部分。P 和E 间的热阻可表示为：
()()()E
e E
e P
e k x k x k x +
+
-≈
+δδδ （5-10）
从中可以看出与调和平均一致。
令 ()()e
e e
x x f δδ-
=， 则 1
1-????
??+-=E e P
e e
k f k f k （5-11） 因此，总体上看，调和平均要比算术平均更好一些。
5.2 边界条件与源项的处理
式（5-4）导出的离散方程只适用于内部节点。为了对某个特定的导热问题进行计算，还应加上边界条件。传热的边界条件有三类，即
（1）给定边界温度； （2）给出边界热流量；
（3）通过换热系数以及周围流体温度给定边界热流量。
如果是第一类边界条件，则就比较简单。但如果是第二或第三类边界，则需要对边界条件进行处理。在有限差分范围内，有两种处理方法，即对边界点补充代数方程和附加源项法。
首先介绍边界点补充代数方程的方法。先讨论区域离散外节点法的情形。在边界给定热流量，如图5-2所示，即给定第二类边界条件，可表示为：
B L
x q dx
dT
k
== （5-12）
上式可化为： B M
M q x
T T k
=--δ1 （5-13）
即 k
q x T T B
M M ?+=-δ1 (5-14)
图5-2 边界节点控制体
上式的截断误差为一阶，而内节点上如采用中心差分，则截断误差为二阶。一般希望内节点与边界点离散方程截断误差等级保持一致，否则会影响计算结果的准确度。为得出具有二阶截断误差的公式，可采用虚拟点法。即在右边界外虚设M+1点，这样节点M 就可视为内节点，其一阶导数就可采用中心差分，即 B M M q x
T T k
=--+δ21
1 （5-15）
为消去1+M T ，由一维稳态含内热源的控制方程（5-1）可得到在 M 点的离散形式，即 ()
022
1
1=++--+s x T T T k
M M M δ （5-16）
从以上两式消去1+M T 可得 ()k
x q k
s x x T T B M M δδ+???+=-1 (5-17)
其中2
x
x δ=
?，是节点M 所代表的控制体的厚度。
对于第三类边界条件，以()
M f B T T q -=α代入式(5-14)和(5-17)可得相应于一阶与二阶截断误差的节点离散方程：
???
?
?+?
??
?????? ??+=-k x T k x T T f M M αδαδ11 （5-18）
()??
? ??αδ+?
?????
αδ+??δ+=-k x T k x k s x x T T f M M 11 （5-19） 用控制容积平衡法来推导。对图5-2所示边界节点的控制体作能量平衡，即可得
02
1=?+-+-x
s x T T k
q M M B δδ （5-20）
解出M T 即得(5-17)。
为了使第二和第三类边界条件的离散方程具有统一的形式，把边界上的热流密度表示成以下形式W B bT a q -=，其中W T 为边界温度。a 等于给定的热流密度(进入为+)，b =0相当于第二类边界条件；a =f T α，b =α相当于第三类边界条件。
??
?????+??????
?+=-k b x k a x T T M M δδ11，0()x δ （5-21） ()??? ???δ+?????
?
?δ+??δ+=-k b x k a x k s x x T T M M 11，0()
2x δ （5-22） 当区域离散化采用第二种方法（即内节点法）时，边界节点可以看成是第一种区域离散法中当边界条件所代表的控制容积厚度x ?趋近于0时的极限，式（5-22）变为
???????+?????
??+=-k b x k a x T T M M δδ11，0()
2
x δ （5-23） 其中x δ为边界节点与第一个内节点之间的距离。上式虽然在形式上与区域离散外节点方法具有一阶截断误差的公式一样，但它却是区域离散内节点方法中具有二阶截断误差的公式。
离散化过程中源项的处理问题：
如果源项是常数，则在离散方程的建立过程中不带来任何困难。当源项是新求解未知量的函数时，源项处理显得十分重要，有时甚至是数值求解成败的关键所在。目前应用较广泛的一种处理方法是把源项局部线性化，即把源项写成P P C T S S S +=，其中C S 为常数，P S 为S 随T 而变化的曲线在P 点的斜率。源项的线性化处理需要注意以下几点：
1、当源项为未知量函数时，线性化的处理比假定源项为常数更为合理。因为如果S=f(T)，那么把各个控制体中的S 作为常数处理就是将上一次迭代计算得到的T*计算S ，也就是说源
项的计算相对于当前T 的计算有一个滞后，而如果按线性化处理，就不存在这个问题。
2、线性化处理是建立线性代数方程所必须的。因为如果采用更高阶的多项式，则所得到的离散化方程就不可能成为线性代数方程。
3、为了保证代表方程迭代求解的收敛，要求P S ≤0。因为b T a
T a nb
nb P P +=
∑，
V ?-=∑P nb P S a a （下标nb 表示邻点，V ?为控制体的体积），线性代数迭代求解收敛
的一个充分条件是对角占优，这就要求P S ≤0。
4、由代数方程迭代求解的公式为
∑∑?-+=
V
S a
b
T a T P nb
nb
nb P
(5-24)
P S 的大小影响到迭代过程中温度的变化速度，P S 的绝对值越大，所对应系统的惯性越大，
相邻两次迭代之间的变化越小，因而收敛速度下降，但有利于克服迭代过程的发散。
例题：设有一导热型方程，02
2=-T dx T d ，边界条件为x =0，T =0；x =1，1=dx dT 。试
将该区域三等分，分别用外节点法和内节点法求解该问题。
解：（1）采用区域离散方法外节点法时，网格划分如图5-3(a)所示。内节点采用中心差分，即(
)
2
1n 1i 2x n
i n i ?Φ+Φ-Φ++。右端点采用一阶截差时，离散方程为：
图5-3 两种离散方法的网格 对2T 节点：
0)3
1(222123=-+-T T T T => 091
2123=-+-T T T
对3T ：
0)3
1(232234=-+-T T T T => 091
2234=-+-T T T
边界条件：01=T (x =0)
13
134=-T T (x =1) => 31
34=-T T 右端点采用二阶截差时,k
x
q k xs x T T B M M δδ+
?+
=-)(1，即 3116131434?+??-=T T T => 32
291234=-T T
此问题精确解为： )(1
2x x e e e e
T --+=
边界条件采用不同离散表达式，对计算结果会产生影响，由表5-1可知，右端点采用二阶截差时的结果远比采用一阶截差时的结果准确。
表5-1 采用外节点法时的数值解
（2）采用离散方法内节点法求解此问题。网格选取如图5-3(b)，则对于节点2，3，4和5的离散化方程分别为：
0289
32=-
T T 0199
199432=+-T T T 028
18
289543=+-T T T 61
54-=-T T
数值解与精确解的比较见表5-2。
表5-2 采用内节点法时的数值解
5.3 一维非稳态导热
一维非稳态导热控制方程可写为: S x T k x t T c
+??
?
??????=??ρ (5-25)
为了建立其离散化方程，在[t , t +Δt ]时间间隔内对控制体P 作积分，为简便起见，设c ρ与时间无关，则可得到:
()()()()()()()xdt T S S dt x T T k x T T k T T x c t t t P P c t
t t
w W P w e P E e t P t t P P ?++??
?
???δ--δ-=-?ρ???+?+?+ (5-26)
为了完成上式的积分，需对上式右端项中T 如何随时间而变的方式作出选择。常用的有显式、隐式和Crank-Nicolson 三种，可用下式来表示：
[]
[]
t
T f fT t
T f fT Tdt t t t t
t t
?-+=?-+=?+?+?
0)1()1( （5-27）
式中f 是在0与1之间的加权因子。上述积分式最后可化为:
()()
()()()()()()()()
()()()()[]
x
T S S f T S S f x T T k x T T k f x T T k x T T k f T T t x c P P c P P c w W P e e P E e w W P w e P E e P P ?+-+++??
??
??δ--δ--+??????δ--δ-=-??ρ00
0000
11 （5-28） 整理上式可得:
[][
][
]x
S x S f a f a f a T T f fT a T f fT a T a c
P
W E P P W W W E E E P P ?+?-+----+
-++-+=)1()1()1()1()1(0
00
0 （5-29）
其中()e e E x k a δ=
，()w w W x k a δ=，()t
x c a P P ??ρ=0
式(5-29)就是一维非稳态导热的离散化方程，取f =0，1，以及1/2可依次得到显式、隐式以及C-N 格式。在直角坐标系中，当网格均匀时，无内热源、常物性导热问题有三种格式，它们分别为： 1、显式：
()
()
00002
002P W E P W W E E P P W P o E P P
T a a a T a T a T a x
T T T a t T T --++=??+-=?- (5-30) 2、隐式：
()
2
2x
T T T a t T T W P E P P
?+-=?- （5-31） 3、C-N 格式：
()
???
?
???+-+?+-=?-20020
222x
T T T x T T T a t T T
W P o E
W P E P
P （5-32） 采用von Neumann 分析方法可以证明，对于源项不随时间而变的问题，对于式(5-29)，
当0.5≤f≤1时，绝对稳定；而当0≤f <1时，稳定的条件则为
()f x
t a 21212
-≤???，其中c k
a ρ=。 当采用显式时，由式（5-30）可知，T P 并不与T E 和T W 等其它未知数有关。任何f≠0的格式都将是隐式，T P 必定和未知数T E 和T W 有关，并且必须求解一组联立方程。相比之下，显式就比较方便，但这一点将被它局限性所抵消。因为式(5-30)中0
P T 的系数可能为负。事实上，为了使这个系数为正，必须使时间步长t ?小到足以使0
P a 大于W E a a +。对于均匀导热，
均匀网格，可以表达为：t ?≤a
x 22
?。从中可以看出，当为了提高计算精度而减小x ?时，
只能采用更小的t ?。
至此，可以总结出保持离散化方程计算收敛和不失真的四个基本规则:
1、控制容积界面上的连续性原则。在两个控制容积的离散化方程中，通过界面的流量必须用相同的表达式来表示。如对于扩散项，由于取了分段线性分布假设，因而保证了在控制容积界面上热量和热通量的连续性。如果如图5-4所示的二次曲线来计算界面上的热通量（-k d T /d x ），则在公共界面上会造成梯度的不连续性。
图5-4 由二次曲线分布引发的热通量不连续示意图
1-右边斜率；2-左边斜率
2、系数为正原则。某个网格节点处的因变量值只是通过对流和扩散的过程才受到相邻节点上值的影响。如对于一维导热问题，如果T E 的增加必然导致T P 的增加，则要求T E 的系数a E 和T P 的系数a P 具有相同的符号。换言之，中心系数P a 与相邻系数nb a 的符号必须相同。
3、源项的负斜率线性化原则。由式（5-5）可知，当S P 为正值时，中心节点的系数P
a 仍有可能为负值。只有当S P 不为正值时，才能完全避免这种情况的出现。因此，当源项线性化Φ+=P c S S S 处理时，0
过程中的不稳定以及导致产生物理意义上的不真实解，这是至关重要的。
4、相邻节点系数加和原则。基本微分方程往往只包含有因变量的导数项，于是，如果T 代表因变量，则函数T 与T+C （C 为任一常数）两者均满足微分方程，微分方程所具有的这个特性必定要反映到与之相对应的离散化方程中。因此，当T P 以及所有的T nb 都增加同一常数时，离散化方程应该仍然适用，这就要求P a =
∑nb
a
。但当源项与T 有关时，特别进
行线性化处理时，相邻节点系数不服从此原则，这只是一个特例，如S P 为零，仍满足此原则。
需要指出的是在数学上认为是稳定的初值问题的差分格式，未必能保证在所有的时间步长下均获得具有物理意义的解。C-N 格式就属于这种情况。当5.0=f 时，方程中0
P T 的系数变成()2/0
w E P a a a +-。对于均匀导热系数和均匀网格问题，可以得出该系数为
x
k
t x c ?-
??ρ。从中可以看出只要时间步长不是足够小，该系数就有可能变为负，从而有可能出现物理上不真实的解。如果要求方程中0
P T 的系数永远不为负的话，只有f 为1才能确保这个要求。因此，全隐身格式(f =1)能满足格式简单且在物理上不失真的要求。正是因为此，我们通常采用全隐式格式。但必须承认，对小的时间步长，全隐式格式不如C-N 格式精确。
全隐式离散化方程可表示为:
b T a T a T a w w E E P P ++= （5-33）
式中()e e E x k a δ=
, ()w w w x k a δ=, t
x c a P ??ρ=0
, x S a a a a c P w E P ?-++=0, 全隐式格式的主要原理是新值P T 在整个时间步长占控制地位。
5.4二维非稳态导热问题的全隐式格式
在直角坐标系下二维非稳态导热方程为: S y T k y x T k x t T c
+???
?
??????+??? ??????=??ρ （5-34） 在时间间隔[t, t+Δt]内，对图5-4所示的控制容积P 做积分，假设在控制容积界面上热流密度均匀，采用全隐式格式，则有：
1、时间项：
()()
y x T T c dxdydt t
T
c
n s
e w t
t t
P P P ??-ρ=??ρ???
?+0 （5-35）
2、扩散项:
()()()()t
x x T T k x T T k t y x T T k x T T k dxdydt y T k y dxdydt x T k t s S P s n P E n t t t
n s
e w t t t e w n s w W P w e P E e ???????
?δ--δ-+????
?
???δ--δ-=???? ??????+???
????????
?????+?+ （5-36）
图5-5 二维直角坐标系下的网格体系 3、源项:
()t y x T S S Sdxdydt n s
e w t
t t
P P c ???+=???
?+ （5-37）
将（5-35）、（5-36）和（5-37）代入（5-34）并整理可得:
b T a T a T a T a T a S S N N W W E E P P ++++= （5-38）
其中: ()e e E x y k a δ?=
, ()w w W x y k a δ?=, ()n n N y x k a δ?=, ()s s S y x k a δ?=, t
y x c a P P ???ρ=)(0
, 00P P C T a y x S b +??=, y x S a a a a a a P P S N w E P ??-++++=0
y x ??为控制体的容积，界面上当量导热系数按调解和平均计算。
轴对称的圆柱坐标中，非稳态导热问题的控制方程为:
S r T rk r r x T k x t T c
+??
? ??????+??? ??????=??ρ1 （5-39） 取一弧度为中心角所包含的范围作为研究对象，网格体系见图5-6所示，采用类似的推导方法可得形式上与上式完全相同的离散方程式，即
b T a T a T a T a T a S S N N W W E E P P ++++= （5-40）
()e e P E k x r r a δ?=
, ()w w P w k x r
r a δ?=, ()n n n N k r x r a δ?=, ()s
s
s S k r x r a δ?=,
V S a a a a a a P P S N W E P ?-++++=0, t
V c a P P ??ρ=
)(0
, 0
0P P C T a V S b +?=, ()x r r r V s n ??+=?5.0
图5-6 二维柱坐标系下的网格体系
在极坐标θ,r 中, 非稳态导热问题的控制方程为：
S T r k r r T kr r r t T c
+??
?
??θ??θ??+??? ??????=??ρ11 （5-41）
b T a T a T a T a T a S S N N W W E E P P ++++= （5-42）
()e e e E r r k a δθ?=
, ()w w w w r r k a δθ?=, ()n n n N r r k a δθ?=, ()s s s S r r k a δθ?=, t
V c a P ??ρ=0
,
00P P C T a V S b +?=, V S a a a a a a P P
S N W E P ?-++++=0
,()r r r V s n ?θ?+=?5.0
图5-7 极坐标系下的网格体系
5.5三维问题的离散化方程
在z 方向上再加两个相邻近节点T 和B(顶和底),即可完成三维情况的离散化，离散化方程的形式与二维的类似，即
b T a T a T a T a T a T a T a B B T T S S N N W W E E P P ++++++= （5-43）
式中()e e E x z y k a δ??=
, ()w w W x z y k a δ??=, ()n n N y x z k a δ??=, ()s s S y x
z k a δ??=,
()t t T z y x k a δ??=
, ()b b B z y x k a δ??=, t
z y x c a P ????ρ=0
, 0
0P P C T a z y x S b +???=,
z y x S a a a a a a a a P P B T S N W E P ???-++++++=0
5.6多维导热问题边界条件的处理
对于第一类边界条件的问题，由于边界温度为已知值，这样就不必补充边界节点方程，求解的区域仅限于内部节点。对于二维、三维问题的第二、三类边界条件问题，如果也能做到这样，则不仅能节省计算时间，而且有助于编程。目前比较成功的处理方法就是附加源项法。所谓附加源项法就是把由第二类或第三类边界条件所规定的进入或导出计算区域的热量作为与边界相邻的控制容积的当量源项。以直角坐标中的情形为例，如图5-8所示。与边界相邻的控制容积中的节点为P 。此控制容积的离散化方程可表达为：
b T a T a T a T a T a S S N N W W E E P P ++++= （5-44）
图5-8 二维边界控制体
为了在P T 的代数方程中不出现未知的边界温度，就需要利用已知的边界条件把W T 消去，为此将上式作如下变换：
b T T a T a T a T a T a a P W W S S N N E E P W P +-+++=-)()(
由于y q x T T y k T T a B w
P W B P W W ?=δ-?=
-)()
()(
其中B k 为边界节点的导热系数，B q 为进入该控制容积的热流密度，以进入为正。这样关于P 点的方程可化为：
b y q T a T a T a T a B S S N N E E P P +?+++=' （5-45）
对于第二类边界条件，B q 为已知，故可把它与b 组成一个新项：
()y x S S y x y x y q S b y q ad C C B C B ??+=?????
? ?????+=+?,
同时，y x S a a a a a a P S N E W P P ??-+++=-=0'
这就是说，对第二类边界条件，如果把)/(y x y q B ???作为与边界相邻的控制容积的附加常数源项，记为ad C S ,，同时令0=W a ，则所得的离散方程既满足能量守恒关系，又可把未知的边界温度排除在外。
当边界条件为第三类时，即()
W f B T T q -α=。由Fourier 定律可得：
()
()w
P W B B x T T k q δ-=
则可得：
()()B
w P f B w P
W W f B k x T T k x T T T T q //1//1δ+α-=δ-=
α
-=
（5-46）
将上式代入式（5-45）经整理后可得：
()()()y
x k x y x yT S T a T a T a T k
x y a B w f C S S N N E E P B w P ??????
?
?δ+α???++++=???
? ?
?δ+α?+//1//1' （5-47）
上式表明，对第三类边界条件，如果在边界控制容积中加入以下附加源项： ()B
w f ad
C k x T y x y
S //1,δ+α????= ()B
w ad P k x y x y S //11
,δ+α????-
= 同时令0=W a 就可实现使未知的边界温度不进入离散方程的目的。为了编程方便，在采用
附加源项法时，边界节点控制容积的系数W a 是以节点上的导热系数B k 来计算，这样按已知的边界导热系数确定ad ad C S S ,,,后，只要令0 B k ，就可实现附加源项法。在或内部节点的解后，由式（5-46）可确定边界节点温度。
需要注意的是附加源项法仅适用于区域离散化方法的内节点法。
三、编程题 4.16 设计工程，已知圆的半径r，求圆面积S。 【解答】设圆半径为r，圆面积为S。根据数学知识，已知圆半径r，求圆面积S的公式为：2r Sπ =。 设计步骤如下。 （1）建立应用程序用户界面，如图4-1所示。 （2）设置对象属性： Label1的Caption属性为“已知圆半径r=”； Text1的Text属性为空； Command1的Caption属性为“圆面积为：”； Label2的Caption属性为空； Label2的BorderStyle属性为1-Fixed Single。 各控件的属性设置如图4-2所示。 图4-1 建立用户界面图4-2 设置各控件的属性（3）编写程序代码。 写出“圆面积为：”命令按钮Command1的Click事件代码为： Private Sub Command1_Click( ) Const pi = 3.14 Dim r As Single, S As Single r = V al(Text1.Text) S = pi * r ^ 2 Label2.Caption = S End Sub 运行程序时，在文本框输入圆半径的值，单击“圆面积为：”按钮后，输出结果如图4-3所示。 也可以不用文本框接收输入值，改用InputBox函数接收圆的半径r，求圆面积S，代码如下。 图4-3 程序运行结果 Private Sub Form_Load( ) Show Const pi = 3.1415926
Dim r As Single, S As Single r = V al(InputBox("输入半径：", "计算圆面积", "10")) FontSize = 18 S = pi * r ^ 2 Print "圆面积："; S End Sub 程序运行时，首先显示如图4-4所示的对话框，在该对话框的文本框中输入数字，按Enter 键或单击“确定”按钮后，才能显示窗体。 图4-4 输入对话框 用InputBox 函数输入文本虽然很方便，但是由于输入框弹出后将暂停程序的运行，直到用户响应，因此输入框不符合VB 自由环境的精神。输入框适合于像要求用户输入口令等这样不常见的输入方式。还可以用更好的用户输入方式，如文本框、选项按钮等。 4.17 已知平面坐标系中两点的坐标，求两点间的距离。 【解答】 由数学知识可知，已知两点坐标（x A , y A ）、（x B , y B ），求两点间距离的计算公式为 2 A B 2 A B )()(y y x x s -+-= 建立用户界面如图4-5所示。在该界面中用TextBox 控件输入数据，用Label 控件输出数据。为了形象地表示两点之间的距离，可用Picture 控件插入一幅图，该图用画图软件绘制。 命令按钮Command1的Click 事件代码为： Private Sub Command1_Click( ) Dim xa As Single, xb As Single Dim ya As Single, yb As Single Dim s As Single xa = Val(Text1.Text) ya = V al(Text2.Text) xb = V al(Text3.Text) yb = V al(Text4.Text) s = Sqr((xb - xa) ^ 2 + (yb - ya) ^ 2) Label6.Caption = s End Sub 程序运行结果如图4-6所示。
第四章 数值积分与数值微分 1.确定下列求积公式中的特定参数，使其代数精度尽量高，并指明所构造出的求积公式所具有的代数精度： 101210121 12120 (1)()()(0)(); (2)()()(0)(); (3)()[(1)2()3()]/3; (4)()[(0)()]/2[(0)()]; h h h h h f x dx A f h A f A f h f x dx A f h A f A f h f x dx f f x f x f x dx h f f h ah f f h -----≈-++≈-++≈-++''≈++-?? ?? 解： 求解求积公式的代数精度时，应根据代数精度的定义，即求积公式对于次数不超过m 的多项式均能准确地成立，但对于m+1次多项式就不准确成立，进行验证性求解。 （1）若101(1) ()()(0)()h h f x dx A f h A f A f h --≈-++? 令()1f x =，则 1012h A A A -=++ 令()f x x =，则 110A h Ah -=-+ 令2 ()f x x =，则 3 221123 h h A h A -=+ 从而解得 01 1431313A h A h A h -?=?? ?=?? ?=?? 令3 ()f x x =，则 3()0h h h h f x dx x dx --==? ? 101()(0)()0A f h A f A f h --++= 故 101()()(0)()h h f x dx A f h A f A f h --=-++? 成立。 令4 ()f x x =，则
取步长δx=0.02。已知x=0,Φ=0;x=1,Φ=1.令k=ρｕ/Γ计算结果图表： 程序及数据结果： 追赶法： #include a[i]=2+0.02*k; b[i]=4; c[i]=2-0.02*k; f[i]=0; } tdma(a,b,c,f,x); for(i=0;i 4-1 解：采用区域离散方法A 时；网格划分如右图。内点采用中心差分 23278.87769.9 T T T === 22d T T=0dx - 有 i+1i 12 2+T 0i i T T T x ---=? 将2点,3点带入 32122 2+T 0T T T x --=? 即321 209T T -+= 432322+T 0T T T x --=?4321322+T 0T T T x --=? 即4 321 209 T T T -+-= 边界点4 （1）一阶截差 由x=1 1dT dx =，得 431 3 T T -= （2）二阶截差 11B M M q x x x T T S δδλλ -=++ 所以 434111. 1. 36311 T T T =++ 即 431 22293 T T -= 采用区域离散方法B 22d T T=0dx - 由控制容积法 0w e dT dT T x dT dT ????--?= ? ????? 所以代入2点4点有 322121011336 T T T T T ----= 即 239 028T T -= 544431011363 T T T T T ----= 即 34599 02828T T T -+= 对3点采用中心差分有 432 32 2+T 013T T T --=?? ??? 即 23499 01919 T T T -+= 对于点5 由x=1 1dT dx =，得 541 6 T T -= （1）精确解求左端点的热流密度 由 ()2 1 x x e T e e e -= -+ 所以有 ()22 20.64806911x x x x dT e e q e e dx e e λ -====- +=-=++ （2）由A 的一阶截差公式 21 0.247730.743113 x T T dT q dx λ =-=-= =?= （3）由B 的一阶截差公式 0 0.21640 0.649213 x dT q dx λ =-=-= = （4）由区域离散方法B 中的一阶截差公式： 210.108460.6504()B B T T dT dx x δ-?? ==?= ? ?? 通过对上述计算结果进行比较可得：区域离散B 有控制容积平衡法建立的离散方程与区域离散方程A 中具有二阶精度的格式精确度相当！ 4-3 解：将平板沿厚度方向3等分，如图 第四章 数值积分与数值微分 1.确定下列求积公式中的特定参数，使其代数精度尽量高，并指明所构造出的求积公式所具有的代数精度： 101210121 12120 (1)()()(0)(); (2)()()(0)(); (3)()[(1)2()3()]/3; (4)()[(0)()]/2[(0)()]; h h h h h f x dx A f h A f A f h f x dx A f h A f A f h f x dx f f x f x f x dx h f f h ah f f h -----≈-++≈-++≈-++''≈++-?? ?? 解： 求解求积公式的代数精度时，应根据代数精度的定义，即求积公式对于次数不超过m 的多项式均能准确地成立，但对于m+1次多项式就不准确成立，进行验证性求解。 （1）若101(1) ()()(0)()h h f x dx A f h A f A f h --≈-++? 令()1f x =，则 1012h A A A -=++ 令()f x x =，则 110A h Ah -=-+ 令2 ()f x x =，则 3 221123 h h A h A -=+ 从而解得 011431313A h A h A h -?=?? ? =?? ?=?? 令3 ()f x x =，则 3()0h h h h f x dx x dx --==? ? 101()(0)()0A f h A f A f h --++= 令4()f x x =，则 455 1012()5 2 ()(0)()3 h h h h f x dx x dx h A f h A f A f h h ---== -++=? ? 故此时， 101()()(0)()h h f x dx A f h A f A f h --≠-++? 故 101()()(0)()h h f x dx A f h A f A f h --≈-++? 具有3次代数精度。 （2）若 21012()()(0)()h h f x dx A f h A f A f h --≈-++? 令()1f x =，则 1014h A A A -=++ 令()f x x =，则 110A h Ah -=-+ 令2 ()f x x =，则 3 2211163 h h A h A -=+ 从而解得 1143 8383A h A h A h -?=-?? ? =?? ?=?? 令3 ()f x x =，则 22322()0h h h h f x dx x dx --==? ? 101()(0)()0A f h A f A f h --++= 导热问题的数值解法 1 、重点内容：① 掌握导热问题数值解法的基本思路； ② 利用热平衡法和泰勒级数展开法建立节点的离散方程。 2 、掌握内容：数值解法的实质。 3 、了解内容：了解非稳态导热问题的两种差分格式及其稳定性。 由前述3 可知，求解导热问题实际上就是对导热微分方程在定解条件下的积分求解，从而获得分析解。但是，对于工程中几何形状及定解条件比较复杂的导热问题，从数学上目前无法得出其分析解。随着计算机技术的迅速发展，对物理问题进行离散求解的数值方法发展得十分迅速，并得到广泛应用，并形成为传热学的一个分支——计算传热学（数值传热学），这些数值解法主要有以下几种：（１）有限差分法（ 2 ）有限元方法（ 3 ）边界元方法 数值解法能解决的问题原则上是一切导热问题，特别是分析解方法无法解决的问题。如：几何形状、边界条件复杂、物性不均、多维导热问题。 分析解法与数值解法的异同点： 相同点：根本目的是相同的，即确定① t=f(x ，y ，z) ；②。不同点：数值解法求解的是区域或时间空间坐标系中离散点的温度分布代替连续的温度场；分析解法求解的是连续的温度场的分布特征，而不是分散点的数值。§4-1 导热问题数值求解的基本思想及内节点离散方程的建立 实质 对物理问题进行数值解法的基本思路可以概括为：把原来在时间、空间坐标系中连续的物理量的场，如导热物体的温度场等，用有限个离散点上的值的集合来代替，通过求解按一定方法建立起来的关于这些值的代数方程，来获得离散点上被求物理量的值。该方法称为数值解法。 这些离散点上被求物理量值的集合称为该物理量的数值解。 2 、基本思路：数值解法的求解过程可用框图4-1 表示。 由此可见： 1 ）物理模型简化成数学模型是基础； 2 ）建立节点离散方程是关键； 3 ）一般情况微分方程中，某一变量在某一坐标方向所需边界条件的个数等于该变量在该坐标方向最高阶导数的阶数。 一数值求解的步骤 如图4-2 （a ），二维矩形域内无内热源、稳态、常物性的导热问题采用数值解法的步骤如下： 1 建立控制方程及定解条件 控制方程：是指描写物理问题的微分方程 针对图示的导热问题，它的控制方程（即导热微分方程）为：（a ）边界条件：x=0 时， x=H 时， 当y=0 时， 当y=W 时， 区域离散化（确立节点） 第四章导热问题的数值解法 1 、重点内容：①掌握导热问题数值解法的基本思路； ②利用热平衡法和泰勒级数展开法建立节点的离散方程。 2 、掌握内容：数值解法的实质。 3 、了解内容：了解非稳态导热问题的两种差分格式及其稳定性。 §4—1导热问题数值求解的基本思想及内节点方程的建立由前述 3 可知，求解导热问题实际上就是对导热微分方程在定解条件下的积分求解，从而获得分析解。但是，对于工程中几何形状及定解条件比较复杂的导热问题，从数学上目前无法得出其分析解。随着计算机技术的迅速发展，对物理问题进行离散求解的数值方法发展得十分迅速，并得到广泛应用，并形成为传热学的一个分支——计算传热学（数值传热学），这些数值解法主要有以下几种： （１）有限差分法（ 2 ）有限元方法（ 3 ）边界元方法 数值解法能解决的问题原则上是一切导热问题，特别是分析解方法无法解决的问题。如：几何形状、边界条件复杂、物性不均、多维导热问题。 一．分析解法与数值解法的异同点： ?相同点：根本目的是相同的，即确定① t=f(x ， y ， z) ；② 。 ?不同点：数值解法求解的是区域或时间空间坐标系中离散点的温度分布代替连续的温度场；分析解法求解的是连续的温度场的分布特征，而不是分散点的数值。 数值求解的基本思路及稳态导热内节点离散方程的建立 二．解法的基本概念 ?实质 对物理问题进行数值解法的基本思路可以概括为：把原来在时间、空间坐标系中连续的物理量的场，如导热物体的温度场等，用有限个离散点上的值的集合来代替，通过求解按一定方法建立起来的关于这些值的代数方程，来获得离散点上被求物理量的值。该方法称为数值解法。 这些离散点上被求物理量值的集合称为该物理量的数值解。 2 、基本思路：数值解法的求解过程可用框图 4-1 表示。 由此可见： 1 ）物理模型简化成数学模型是基础； 2 ）建立节点离散方程是关键； 3 ）一般情况微分方程中，某一变量在某一坐标方向所需边界条件的个数等于该变量在该坐标方向最高阶导数的阶数。 ?数值求解的步骤 如图 4-2 （ a ），二维矩形域内无内热源、稳态、常物性的导热问题采用数值解法的步骤如下：（1）建立控制方程及定解条件 控制方程：是指描写物理问题的微分方程 针对图示的导热问题，它的控制方程（即导热微分方程）为：（ a ） 边界条件： x=0 时， x=H 时， 当 y=0 时， 导热问题数值解法初次研究 对物理物体的数值求解的基本思想可以概括为：把原有的时间、空间坐标系中连续的物理量的场，如导热问题的温度场，用有限个离散点上的值的集合来代替，通过求解按一定方法建立起来的关于这些值的代数方程，来获得离散点上的值。这些离散点上的被求解物理量的值的集合称为该物理量的数值解。 物理模型 在四个输气的管道中间有一个各边长为10厘米的薄铁片，求导热其达到稳态后，这块铁片的温度分布。四个输气管道里的气体温度是恒值分别为100℃、200℃、500℃、1000℃。因此，可以看成是二维矩形域内的稳态、无内热源、常物性的导热问题。 建立数学模型 描写物理问题的微分方程称为控制方程，导热微分方程为： 2 2 22 0t t x y ??+ =?? （1） 其四个边界分别为第一类边界条件，1234t =1005002001000===℃、t ℃、t ℃、t ℃。 区域离散化 用一系列与坐标轴平行的网格线把求解区域划分成许多子区域，以网格线的交点作为需要确定温度值的空间位置，称为节点。相邻两节点的距离称为步长，记为x ?、y ?。本模型x 、y 方向是各自均分的，各自为100个子区域。节点的位置以该节点在两个方向上的标号m 、n 来表示。 每一个节点都可以看成是以它为中心的一个小区域的代表，由相邻的两节点连接的中垂线构成。为叙述方便，我们把节点所代表的小区域称为元体。 数学模型离散化 它的建立是数值求解过程中的重要环节，主要有泰勒级数展开法及热平衡法两种，取节点（m ，n ）及其临点为例。 泰勒级数展开法 以节点（m ，n ）处的二阶偏导数为例用这种方法来导出其差分表达式。对节点(1,)m n +及(1,)m n -分别写出函数t 对（m ，n ）的泰勒级数展开式： 2 2 33 4 4 1,,,,,,2 3 4 2624m n m n m n m n m n m n t x t x t x t t t x x x x x +???????=+?+ + + +???? （2） 2 2 3 34 4 1,,,,,,2 3 4 2 624m n m n m n m n m n m n t x t x t x t t t x x x x x -???????=-?+ - + +???? （3） 将式（2）、（3）相加得 第四章 习题 1. 采用数值计算方法，画出dt t t x y x ?= 0sin )(在]10 ,0[区间曲线，并计算)5.4(y 。 〖答案〗 1.6541 2. 求函数 x e x f 3sin )(=的数值积分?=π 0 )(dx x f s ，并请采用符号计算尝试复算。 〖答案〗 s = 5.1354 Warning: Explicit integral could not be found. > In sym.int at 58 s = int(exp(sin(x)^3),x = 0 .. pi) 3. 用quad 求取dx x e x sin 7.15? --ππ的数值积分，并保证积分的绝对精度为910-。 〖答案〗 1.08784943754779 4. 求函数 5.08.12cos 5.1)5(sin )(20 6.02++-=t t t e t t f t 在区间]5,5[-中的最小值点。 〖答案〗 最小值点是 -1.28498111480531 相应目标值是 -0.18604801006545 5. 设 0)0(,1)0(,1)(2)(3)(22===+-dt dy y t y dt t dy dt t y d ，用数值法和符号法求5.0)(=t t y 。 〖答案〗 数值解 y_05 = 0.78958020790127 符号解 ys = 1/2-1/2*exp(2*t)+exp(t) ys_05 = .78958035647060552916850705213780 6. 求矩阵b Ax =的解，A 为3阶魔方阵，b 是)13(?的全1列向量。 〖答案〗 x = 0.0667 0.0667 0.0667 7. 求矩阵b Ax =的解，A 为4阶魔方阵，b 是)14(?的全1列向量。 〖答案〗 解不唯一 x = -0.0074 -0.0809 0.1397 0.0662 0.0588 0.1176 -0.0588 4-5迭代法求解节点温度。 说明：此处给出的是C++程序代码，使用牛顿迭代法，迭代收敛精度1.0e-6；程序运行结果附后。 /*NHT 4-5 newton *created on 2012-10-19 by Sanye */ #include 数值传热学陶文铨第四章作业
数值分析第四章数值积分与数值微分习题答案
传热学导热问题的数值解法
第四章导热题的数值解法
罗大雷-导热问题的数值解法
数值分析第四章习题
数值传热学第四章编程题