数值传热学部分习题答案
传热学课后答案(完结版)
2
tw2
3
tw1 tw 2 q2 1 2 3 1 2 3
再由：
tw1
λ
λ 3
tw2
q1
q2 0.2q1 ，有
tw1 tw 2 t t 0.2 w1 w 2 1 2 1 2 3 1 2 1 2 3
得：
3 43 (
'2 3 2 5 6 2 R 0.265m k / W 2 3 0.65 0.024
"
由计算可知，双 Low-e 膜双真空玻璃的导热热阻高于中空玻璃，也就是说双 Low-e 膜双真 空玻璃的保温性能要优于中空玻璃。 3． 4．略 5 ．
m2
(m 2 K )
、 h2 85W
(m 2 K )
、 t1 45 ℃
t2 500 ℃、 k ' h2 85W
求： k 、 、
(m 2 K )
、 1mm 、 398 W
(m K )
解：由于管壁相对直径而言较小，故可将此圆管壁近似为平壁 即： k
tw1 t w 2 x
(设 tw1 tw 2 ) ， 否则 t 与平壁 coust （即常物性假设）
其与平壁的材料无关的根本原因在 的材料有关 （2）由 4.略
q
dt dx
知，q 与平壁的材料即物性有关
5.解：
d 2 dt (r )0 dr dr r r1 , t tw1 (设tw1 t w 2 ) r r2 , t tw 2
绪论
思考题与习题（ P89 ）答案： 1． 冰雹落体后溶化所需热量主要是由以下途径得到：
数值传热课后学习题及答案(权威版)
习题 4-14
充分发展区的温度控制方程如下：
c
p
u
T x
1 r
r
(r
T r
)
对于三种无量纲定义 T Tw 、 T T 、 T Tw 进行分析如下
Tb Tw
Tw T
T Tw
1）由 T Tw 得： Tb Tw
T (Tb Tw ) Tw
由 T 可得：
T [(Tb Tw ) Tw ] Tb (1 ) Tw
x
x
x
x
T r
[(Tb
Tw ) Tw ] r
(Tb
Tw
)
r
(1 ) Tw r
由 Tb 与 r
无关、
与 x 无关以及
T x
、
T r
的表达式可知，除了 Tw 均匀的情况外，该无量
纲温度定义在一般情况下是不能用分离变量法的；
2）由 T T 得： Tw T
T (Tw T ) T
由 T 可得：
节点 3：
T2 4T3 75
求解结果：
T2 85 ， T3 40
对整个控制容积作能量平衡，有：
qB Sx h f (T f T3 ) Sx 15 (20 40) 150 2 0
即：计算区域总体守恒要求满足
习题 4－5
在 4－2 习题中，如果 h 10 (T3 T f )0.25 ，则各节点离散方程如下：
coematrix(n,n)=A(1,n); if n>=2
coematrix(n,n-1)=-1*B(1,n); end if n<mdim
coematrix(n,n+1)=-1*C(1,n); end end %计算 D 矢量 D=(coematrix*T')'; %由已知的 A、B、C、D 用 TDMA 方法求解 T %消元 P(1,1)=C(1,1)/A(1,1); Q(1,1)=D(1,1)/A(1,1); for n=2:mdim P(1,n)=C(1,n)/(A(1,n)-B(1,n)*P(1,n-1)); Q(1,n)=(D(1,n)+B(1,n)*Q(1,n-1))/(A(1,n)-B(1,n)*P(1,n-1)); end %回迭 Tcal(1,mdim)=Q(1,mdim); for n=(mdim-1):-1:1 Tcal(1,n)=P(1,n)*Tcal(1,n+1)+Q(1,n);
数值传热学(陶文铨)第二章习题答案
p1 p2
（1）
p 2 p x2 x 2 x i 2 x i 2 2!
p 2 p x2 p3 p2 x 2 x i 2 x i 2 2!
（2） 式（2）-（1）得
p p1 p 3 x i 2 2 x
2-11.解：对于均分网格用泰勒级数展开法，用 2 点分别表示 1,3,4 处热流量得
（i+1,n）= （i,n）+
(1)
x
( x)
i ,n
2 x 2
i ,n
(( x) )2 2!
（i-1,n）= （i,n）+
(2) 式（1）-（2）得：
x
( x)
i ,n
2 x 2
i ,n
(( x) )2 2!
（i+1,n）- （i-1,n）=
即
krp k krp k 2k rpTp ( )Te ( )T r s r r r 2 r 2 w p
化简后得
2krp
r
Tp
kre kr S Te w Tw rp r r r 2
计算结果与控制容积积分法一致。 2-9.解：对于均分网格用泰勒级数展开法分别表示 1 和 3 点处的压力值
（1） 对 T 随 r 由 Tw 变到 Te 的过程进行积分
（2） 可化为
rk
dT dr
e
w
S 2 r 2 w
e
（3） 取 T 随 r 呈分段线性的变化，则（3）式中
Tw
Tp Tw 2 Tp Te 2
（4）
Te
（5）
dT Te Tp r dr e
（6）
数值传热学答案范文
数值传热学是热力学的重要分支之一，研究物质中热量的传递和分布规律。
与传统的实验方法相比，数值传热学采用计算机模拟技术，通过数学模型和计算实验方法，能够更加深入、系统地研究热传递现象的规律和特性，为工程设计和实际生产提供重要的技术支持。
数值传热学的本质是热传递方程的数值求解。
热传递方程是描述物质中热量传递和分布的方程，它包含了热传导、热对流和热辐射三种传热方式。
热传导是指热量沿着物质内部的温度梯度传递，主要发生在固体和液体中；热对流是指热量随物质的流动而传递，主要发生在液体和气体中；热辐射是指热量通过辐射传递，主要发生在光学和辐射热转换材料中。
通过数值方法求解热传递方程，可以得到物体的温度分布、热传递速率和热流密度等参数，为材料和工程设计提供准确的数据支持。
数值传热学的核心是数值方法，主要包括有限差分、有限元和边界元等方法。
有限差分法是一种利用离散化方法求解微分方程的数值方法，它将微分方程中的连续变量离散化，将求解微分方程转化为求解线性方程组。
有限元法是一种利用有限元逼近方法解决偏微分方程的数值方法，采用对物体进行简单的几何划分，将问题离散化，通过数学建模来表示物体的温度分布和热流密度分布。
边界元法是一种较新的有限元法补充，它能够快速解决边界值问题，并且可以减少问题的维数。
数值传热学的应用范围广泛，包括热工和物理问题的研究、能源系统分析和设计、建筑工程中的热传递和能源效率研究等。
例如，在太阳能发电系统设计中，数值传热学可以帮助设计人员确定集热器表面温度和吸收率等参数，提高太阳能效率并减少系统成本。
在建筑工程中，数值传热学可以帮助设计师分析建筑物的保温性能，合理评估保温材料的性能和使用效果，确保建筑节能和环保。
在机械加工领域中，数值传热学可以帮助工程师分析材料切削过程中的热量和温度分布，挑选适合材料和刀具的加工工艺，提高机械切削效率。
数值传热学是现代科学技术的重要分支之一，是研究物质中热传递和分布规律的重要工具。
传热学部分习题答案(第五版)
6.一厚度为50mm的无限大平壁，其稳态温度分布为 ℃。式中a=200℃，b=－2000℃/m2。若平壁材料导热系数为45W/m.℃，试求：(1)平壁两侧表面处的热流通量；(2)平壁中是否有内热源?为什么？若有的话，它的强度应是多大？
解：（1）利用傅立叶定律确定热流密度：
（2）
或：根据 ，得：
补充题：
1.为什么潮湿的多孔材料的导热系数不但比干多孔材料的导热系数大而且还比水的导热系数大?
答：因为水分的渗入替代了相当一部分空气，而且更重要的是，当潮湿材料有温度梯度时，水和湿汽顺热流方向迁移(即从高温区向低温区迁移)，同时携带很多热量，温度梯度越大，这种“迁移”量就越大，致使潮湿多孔材料的导热系数不仅比干多孔材料的导热系数值大，而且比水的导热系数值还大。
9.一半径为R的实心球，初始温度均匀并等于 ，突然将其放入一温度恒定并等于 的液体槽内冷却。已知球的热物性参数 、 和 ，球壁表面的对流换热系数为h，试写出描写球体冷却过程的完整数学描述。
解：该球经历了一维、非稳态、无内热源、导热系数为常数的导热现象，其完整的数学描述为：
补充题：一个质量为M、具有常物性(导热系数λ，比热c，密度ρ)、无内热源的长方体，该物体开始处于均匀温度t0，然后将一电加热器突然通电，在表面x＝0处提供均匀热流密度q0，而在x＝L处和其他边界都完全绝热。试写出该导热现象完整的数学描述。
冬季挂上窗帘减少了通过窗户的热辐射散热，因此人感觉暖和。
9.一般保温瓶胆为真空玻璃夹层，夹层内两侧镀银，为什么它能较长时间地保持热水的温度?并分析热水的热量是如何通过胆壁传到外界的?什么情况下保温性能会变得很差？
答：保温瓶胆为真空玻璃夹层，其目的是保证夹层散热方式仅是热辐射而没有对流换热方式，同时夹层内两测镀银是为了提高表面反射率，以降低热辐射散热，因此保温瓶可以较长时间地保持热水温度。
《传热学》课后习题答案-第三章
第三章思考题1. 试说明集总参数法的物理概念及数学处理的特点答：当内外热阻之比趋于零时，影响换热的主要环节是在边界上的换热能力。
而内部由于热阻很小而温度趋于均匀，以至于不需要关心温度在空间的分布，温度只是时间的函数， 数学描述上由偏微分方程转化为常微分方程、大大降低了求解难度。
2. 在用热电偶测定气流的非稳态温度场时,怎么才能改善热电偶的温度响应特性?答：要改善热电偶的温度响应特性，即最大限度降低热电偶的时间常数,形状上要降低体面比，要选择热容小的材料，要强化热电偶表面的对流换热。
3. 试说明”无限大平板”物理概念,并举出一二个可以按无限大平板处理的非稳态导热问题 答；所谓“无限大”平板，是指其长宽尺度远大于其厚度，从边缘交换的热量可以忽略 不计，当平板两侧换热均匀时，热量只垂直于板面方向流动。
如薄板两侧均匀加热或冷却、 炉墙或冷库的保温层导热等情况可以按无限大平板处理。
4. 什么叫非稳态导热的正规状态或充分发展阶段?这一阶段在物理过程及数学处理上都有些什么特点?答：非稳态导热过程进行到一定程度,初始温度分布的影响就会消失，虽然各点温度仍 随时间变化,但过余温度的比值已与时间无关，只是几何位置()和边界条件(Bi 数) 的函数，亦即无量纲温度分布不变，这一阶段称为正规状况阶段或充分发展阶段。
这一阶段的数学处理十分便利，温度分布计算只需取无穷级数的首项进行计算。
5. 有人认为,当非稳态导热过程经历时间很长时,采用图3-7记算所得的结果是错误的.理由是： 这个图表明,物体中各点的过余温度的比值与几何位置及Bi 有关,而与时间无关.但当时间趋于无限大时,物体中各点的温度应趋近流体温度,所以两者是有矛盾的。
你是否同意这种看法,说明你的理由。
答：我不同意这种看法，因为随着时间的推移，虽然物体中各点过余温度的比值不变 但各点温度的绝对值在无限接近。
这与物体中各点温度趋近流体温度的事实并不矛盾。
6. 试说明Bi 数的物理意义。
数值传热学 第六章答案 (2)
数值传热学第六章答案简介本文档将为读者提供《数值传热学》第六章的答案。
第六章主要涉及热对流传热的数值计算方法，包括网格划分、边界条件、离散方法等内容。
通过本文档，读者将了解如何使用数值方法解决热对流传热问题，并学会应用这些方法进行实际计算。
问题回答1. 简述热对流传热的数值计算方法。
热对流传热的数值计算方法主要包括三个步骤：网格划分、边界条件设置和离散方法。
网格划分是指将传热区域划分为若干个离散的小单元，每个单元内部温度变化均匀。
常见的网格划分方法有结构化网格和非结构化网格。
结构化网格适用于简单几何形状，易于处理；非结构化网格则适用于复杂几何形状。
边界条件设置是指给定物体表面的边界条件，如温度或热流密度。
边界条件的设置需要根据实际问题来确定，可以通过实验或经验公式来获取。
离散方法是指将传热控制方程进行离散化，通常使用有限差分法或有限元法。
有限差分法将控制方程离散化为代数方程组，而有限元法则通过近似方法将方程离散化。
2. 什么是结构化网格和非结构化网格？它们在热对流传热计算中有何不同？结构化网格是指由规则排列的矩形或立方体单元组成的网格。
在结构化网格中，每个单元与其相邻单元之间的联系都是固定的，因此易于处理。
结构化网格适用于简单几何形状，如长方体或圆柱体。
非结构化网格是指由不规则形状的三角形、四边形或多边形组成的网格。
在非结构化网格中，每个单元与其相邻单元之间的联系可能是不确定的，需要使用邻接表来表示网格拓扑关系。
非结构化网格适用于复杂几何形状，如复杂流体流动中的腔体或障碍物。
在热对流传热计算中，结构化网格和非结构化网格的主要区别在于网格的配置方式和计算复杂度。
结构化网格由正交单元组成，计算稳定性较高，但对于复杂几何形状的处理能力较差。
非结构化网格可以灵活地适应复杂几何形状，但计算复杂度较高。
3. 如何设置边界条件？边界条件的设置是热对流传热计算中非常重要的一步，它决定了计算结果的准确性和可靠性。
数值传热学习题答案(汇总版)
2-4-9
= rP rS
式（2-4-9）也可以写成 a PTP = a E TE + aW TW + b 的形式。而且两种结果是一致的。
2—6：
n n TE −TW dT P , n = 解：将 ， dx 2x n n TE −2TPn + TW d 2T P , n = , dx2 x 2
dk = f (x ) 代入原方程，得： dx
令
2-4-4
rk rk a E = ， aW = ， a P = a E + aW ， b x w x e
= SrP r ，
式（2-4-4）可以写成 a PTP = a E TE + aW TW + b 的形式。 2. 再用 Taylor 展开法导出 k
2 2 uE + uP u = ， 2 2 e
2 2 uW + uP u = 2 2 w
t u ut N − uP y = (y ) ， n n
t
t ut u p − uS y = (y ) 。 s s
t
(y ) n = (y ) s = y
n n n n TE −TW TE −2TPn + TW k + f (x ) +S=0 整理得： 2x x 2
4kT P= 2k + xf ( x)T E+2k − xf ( x)T W +2x 2 S
− 2k 时， a E 会成为负值， x 2k 当 f(x）> 时， aW 会成为负值。 x
rk dr = rk r r dr dr dr
w
e
1 d
传热学课后答案【第五版】
绪论思考题与习题（89P -）答案：1． 冰雹落体后溶化所需热量主要是由以下途径得到：Q λ—— 与地面的导热量 f Q ——与空气的对流换热热量注：若直接暴露于阳光下可考虑辐射换热，否则可忽略不计。
6． 夏季：在维持20℃的室内，人体通过与空气的对流换热失去热量，但同时又与外界和内墙面通过辐射换热得到热量，最终的总失热量减少。
（T T 〉外内）冬季：在与夏季相似的条件下，一方面人体通过对流换热失去部分热量，另一方面又与外界和内墙通过辐射换热失去部分热量，最终的总失热量增加。
（T T 〈外内）挂上窗帘布阻断了与外界的辐射换热，减少了人体的失热量。
7．热对流不等于对流换热，对流换热 = 热对流 + 热传导 热对流为基本传热方式，对流换热为非基本传热方式 8．门窗、墙壁、楼板等等。
以热传导和热对流的方式。
9．因内、外两间为真空，故其间无导热和对流传热，热量仅能通过胆壁传到外界，但夹层 两侧均镀锌，其间的系统辐射系数降低，故能较长时间地保持热水的温度。
当真空被破坏掉后，1、2两侧将存在对流换热，使其保温性能变得很差。
10．t R R A λλ= ⇒ 1t R R Aλλ==2218.331012m --=⨯ 11．q t λσ=∆ c o n s t λ=→直线 c o n s t λ≠ 而为λλ=（t ）时→曲线12. i R α 1R λ 3R λ 0R α 1f t −−→ q首先通过对流换热使炉子内壁温度升高，炉子内壁通过热传导，使内壁温度生高，内壁与空气夹层通过对流换热继续传递热量，空气夹层与外壁间再通过热传导，这样使热量通过空气夹层。
（空气夹层的厚度对壁炉的保温性能有影响，影响a α的大小。
） 13．已知：360mm σ=、0.61()Wm K λ=∙ 118f t =℃ 2187()Wh m K =∙210f t =-℃ 22124()Wh m K =∙ 墙高2.8m ，宽3m求：q 、1w t 、2w t 、φ 解：1211t q h h σλ∆=++＝18(10)45.9210.361870.61124--=++2W m111()f w q h t t =-⇒ 11137.541817.5787w f q t t h =-=-=℃222()w f q h t t =-⇒ 22237.54109.7124w f q t t h =+=-+=-℃ 45.92 2.83385.73q A W φ=⨯=⨯⨯=14.已知：3H m =、0.2m σ=、2L m =、45λ=()W m K ∙ 1150w t =℃、2285w t =℃求：t R λ、R λ、q 、φ解：40.27.407104532t K R W A HL λσσλλ-====⨯⨯⨯30.24.4441045t R λσλ-===⨯2m K W ∙ 3232851501030.44.44410t KW q m R λ--∆-==⨯=⨯ 3428515010182.37.40710t t KW R λφ--∆-==⨯=⨯ 15．已知：50i d mm =、 2.5l m =、85f t =℃、273()Wh m K =∙、25110Wq m =求：i w t 、φ()i w f q h t h t t =∆=-⇒iw f qt t h =+51108515573=+=℃0.05 2.551102006.7i Aq d lq Wφππ===⨯⨯=16.已知:150w t =℃、220w t =℃、241.2 3.96()W c m K =∙、1'200w t =℃求： 1.2q 、'1.2q 、 1.2q ∆ 解：12441.2 1.2()()100100w w t t q c ⎡⎤=-⎢⎥⎣⎦44227350273203.96()()139.2100100W m ++⎡⎤=⨯-=⎢⎥⎣⎦12''441.21.2()()100100w w t t qc ⎡⎤=-⎢⎥⎢⎥⎣⎦442273200273203.96()()1690.3100100W m ++⎡⎤=⨯-=⎢⎥⎣⎦'21.2 1.2 1.21690.3139.21551.1Wq q q m ∆=-=-=17．已知：224A m =、215000()Wh m K =∙、2285()Wh m K =∙、145t =℃2500t =℃、'2285()Wk h m K ==∙、1mm σ=、398λ=()Wm K ∙求：k 、φ、∆解：由于管壁相对直径而言较小，故可将此圆管壁近似为平壁 即：12111k h h σλ=++＝3183.5611101500039085-=⨯++2()W m k ∙ 383.5624(50045)10912.5kA t KW φ-=∆=⨯⨯-⨯=若k ≈2h'100k kk-∆=⨯％8583.56 1.7283.56-==％ 因为：1211h h，21h σλ 即：水侧对流换热热阻及管壁导热热阻远小于燃气侧对流换热热阻，此时前两个热阻均可以忽略不记。
数值传热学部分习题答案
习题4－2一维稳态导热问题的控制方程：022=+∂∂S xTλ 依据本题给定条件，对节点2节点3采用第三类边界条件具有二阶精度的差分格式，最后得到各节点的离散方程： 节点1： 1001=T节点2： 1505105321-=+-T T T 节点3：75432=+-T T求解结果：852=T ，403=T对整个控制容积作能量平衡，有：02150)4020(15)(3=⨯--⨯=∆+-=∆+x S T T h x S q f f B即：计算区域总体守恒要求满足习题4－5在4－2习题中，如果25.03)(10f T T h -⨯=，则各节点离散方程如下：节点1： 1001=T节点2： 1505105321-=+-T T T节点3：25.03325.032)20(4015])20(21[-⨯+=-⨯++-T T T T对于节点3中的相关项作局部线性化处理，然后迭代计算； 求解结果：818.822=T ，635.353=T （迭代精度为10-4）迭代计算的Matlab 程序如下： x=30; x1=20;while abs(x1-x)>0.0001a=[1 0 0;5 -10 5;0 -1 1+2*(x-20)^(0.25)]; b=[100;-150; 15+40*(x-20)^(0.25)]; t=a^(-1)*b; x1=x; x=t(3,1);endtcal=t习题4-12的Matlab程序%代数方程形式A i T i=C i T i+1+B i T i-1+D imdim=10;%计算的节点数x=linspace(1,3,mdim);%生成A、C、B、T数据的基数；A=cos(x);%TDMA的主对角元素B=sin(x);%TDMA的下对角线元素C=cos(x)+exp(x); %TDMA的上对角线元素T=exp(x).*cos(x); %温度数据%由A、B、C构成TDMAcoematrix=eye(mdim,mdim);for n=1:mdimcoematrix(n,n)=A(1,n);if n>=2coematrix(n,n-1)=-1*B(1,n);endif n<mdimcoematrix(n,n+1)=-1*C(1,n);endend%计算D矢量D=(coematrix*T')';%由已知的A、B、C、D用TDMA方法求解T%消元P(1,1)=C(1,1)/A(1,1);Q(1,1)=D(1,1)/A(1,1);for n=2:mdimP(1,n)=C(1,n)/(A(1,n)-B(1,n)*P(1,n-1));Q(1,n)=(D(1,n)+B(1,n)*Q(1,n-1))/(A(1,n)-B(1,n)*P(1,n-1));end%回迭Tcal(1,mdim)=Q(1,mdim);for n=(mdim-1):-1:1Tcal(1,n)=P(1,n)*Tcal(1,n+1)+Q(1,n);endTcom=[T;Tcal];%绘图比较给定T值和计算T值plot(Tcal,'r*')hold onplot(T)结果比较如下，由比较可知两者值非常切合（在小数点后8位之后才有区别）：习题4-14充分发展区的温度控制方程如下：)(1rTr r r x T uc p ∂∂∂∂=∂∂λρ 对于三种无量纲定义w b w T T T T --=Θ、∞∞--=ΘT T T T w 、w w T T T T --=Θ∞进行分析如下1）由wb wT T T T --=Θ得：w w b T T T T +Θ-=)(由T 可得：x T x T x T T T x T w b w w b ∂∂Θ-+∂∂Θ=∂+Θ-∂=∂∂)1(])[(rT r T T r T T T r T w w b w w b ∂∂Θ-+∂Θ∂-=∂+Θ-∂=∂∂)1()(])[( 由b T 与r 无关、Θ与x 无关以及x T ∂∂、rT∂∂的表达式可知，除了w T 均匀的情况外，该无量纲温度定义在一般情况下是不能用分离变量法的； 2）由∞∞--=ΘT T T T w 得： ∞∞+Θ-=T T T T w )(由T 可得：xT x T T T x T w w ∂∂Θ=∂+Θ-∂=∂∂∞∞])[(rT r T T r T T T r T w w w ∂∂Θ+∂Θ∂-=∂+Θ-∂=∂∂∞∞∞)(])[( 由b T 与r 无关、Θ与x 无关以及x T ∂∂、rT∂∂的表达式可知，在常见的四种边界条件中除了轴向及周向均匀热流const q w =的情况外，有0=∂∂rT w，则该无量纲温度定义是可以用分离变量法的； 3）由wwT T T T --=Θ∞得： w w T T T T +Θ-=∞)(由T 可得：xT x T T T x T w w w ∂∂Θ-=∂+Θ-∂=∂∂∞)1(])[(rT r T T r T T T r T w w w w ∂∂Θ-+∂Θ∂-=∂+Θ-∂=∂∂∞∞)1()(])[( 同2）分析可知，除了轴向及周向均匀热流const q w =的情况外，有0=∂∂rT w，该无量纲温度定义是可以用分离变量法的；习题4－181）采用柱坐标分析，写出统一的稳态柱坐标形式动量方程：r r x x w r v r r r u x ∂∂+∂∂∂∂=∂∂+∂∂+∂∂1)()(1)(1)(φλφρθφρφρx 、r 和θ分别是圆柱坐标的3个坐标轴，u 、v 和w 管内的流动方向；对于管内的层流充分发展有：0=v 、0=w ，0=∂∂xu； 并且x 方向的源项：x pS ∂∂-=r 方向的源项：r pS ∂∂-= θ方向的源项：θ∂∂-=pr S 1 由以上分析可得到圆柱坐标下的动量方程： x 方向： 0)(1)(1=∂∂-∂∂∂∂+∂∂∂∂x pu r r r u r r r θλθλ r 方向：0=∂∂r pθ方向：0=∂∂θp边界条件： R r =，0=u0=r ，0=∂∂r u ；对称线上，0=∂∂θu 不考虑液体的轴向导热，并简化分析可以得到充分发展的能量方程为：)(1)(1θλθλρ∂∂∂∂+∂∂∂∂=∂∂Tr r r T r r r x T uc p 边界条件： R r =，w q r T =∂∂λ；0=r ，0=∂∂rTπθ/0=，0=∂∂-θλT2）定义无量纲流速：dxdp R uU 2-=λ并定义无量纲半径：R r /=η；将无量纲流速和无量纲半径代入x 方向的动量方程得：0))1((1))1((122=∂∂-∂-∂∂∂+∂-∂∂∂xp U dx dp R R R R U dx dp R RR R θληλθηηλληηη 上式化简得：01)1(1)(1=+∂∂∂∂+∂∂∂∂θηθηηηηηU U 边界条件：1=η，0=U0=η，0=∂∂ηU ；对称线上，0=∂∂θU定义无量纲温度：λ/0R q T T b-=Θ其中，0q 是折算到管壁表面上的平均热流密度，即：Rq q wπ=0； 由无量纲温度定义可得：b T Rq T +Θ=λ将T 表达式和无量纲半径η代入能量方程得：)(1)(100θληλθηηλληηηρ∂Θ∂∂∂+∂Θ∂∂∂=∂∂R q R R R R q R R R x T uc b p 化简得：)1(1)(10θηθηηηηηρ∂Θ∂∂∂+∂Θ∂∂∂=∂∂x T u c q R b p （1）由热平衡条件关系可以得：mm m b m p b p p RU U q R u u R q A u u dx dT A u c x T u c x T uc 020221221)(===∂∂=∂∂ππρρρ 将上式代入式（1）可得：)1(1)(12θηθηηηηη∂Θ∂∂∂+∂Θ∂∂∂=m U U 边界条件：0=η，0=∂Θ∂η；1=η，R q q w πη10==∂Θ∂0=θ，0=∂Θ∂θ；πθ=，0=∂Θ∂θ单值条件： 由定义可知：0/0=-=ΘλR q T T b b b 且： ⎰⎰Θ=ΘAAb UdAUdA即得单值性条件：0=Θ⎰⎰AA UdAUdA 3）由阻力系数f 及Re 定义有：228)(21/Re ⎪⎭⎫ ⎝⎛=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-=D D U D u u dx dp D f e m e m me νρ 且：m W b m W b m W R q T T D T T q Nu ,0,,0~2)/(2Θ=-=-=λλ5－21．一维稳态无源项的对流－扩散方程如下所示：xx u 22∂∂Γ=∂∂φφρ （取常物性）边界条件如下：L L x x φφφφ====,;,00上述方程的精确解如下：11)/(00--=--⋅PeL x Pe L e e φφφφ Γ=/uL Pe ρ 2．将L 分成20等份，所以有：∆=P Pe 201 2 3 4 5 6 ………… …………… 17 18 19 20 21 对于中心差分、一阶迎风、混合格式和QUICK 格式分别分析如下： 1） 中心差分中间节点： 2)5.01()5.01(11-∆+∆++-=i i i P P φφφ 20,2 =i2） 一阶迎风中间节点： ∆-∆++++=P P i i i 2)1(11φφφ 20,2 =i3） 混合格式当1=∆P 时，中间节点：2)5.01()5.01(11-∆+∆++-=i i i P P φφφ20,2 =i当10,5=∆P 时，中间节点： 1-=i i φφ 20,2 =i 4） QUICK 格式*12111)35(8122121⎥⎦⎤⎢⎣⎡---++++++=+--∆∆-∆∆+∆i i i i i i i P P P P P φφφφφφφ 2≠i *1111)336(8122121⎥⎦⎤⎢⎣⎡--++++++=+-∆∆-∆∆+∆i i i i i i P P P P P φφφφφφ 2=i数值计算结果与精确解的计算程序如下：%except for HS, any other scheme doesnt take Pe<0 into consideration %expression of exact solutiony=dsolve('a*b*Dy=c*D2y','y(0)=y0,y(L)=yL','x')y=subs(y,'L*a*b/c','t')y=simple(subs(y,'a*b/c*x','t*X'));ysim=simple(sym(strcat('(',char(y),'-y0)','/(yL-y0)')))y=sym(strcat('(',char(ysim),')*(yL-y0)','+y0'))% in the case of Pe=0y1=dsolve('D2y=0','y(0)=y0,y(L)=yL','x')y1=subs(y1,'-(y0-yL)/L*x','(-y0+yL)*X')%grid Pe numbertt=[1 5 10];%dimensionless lengthm=20;%mdim is the number of inner nodemdim=m-1;X=linspace(0,1,m+1);%initial value of variable during calculationy0=1;yL=2;%cal exact solutionfor n=1:size(tt,2)t=m*tt(1,n);if t==0yval1(n,:)=eval(y1);elseyval1(n,:)=eval(y);endend%extra treatment because max number in MATLAB is 10^308if max(isnan(yval1(:)))yval1=yval1';yval1=yval1(:);indexf=find(isnan(yval1));for n=1:size(indexf,1)if rem(indexf(n,1),size(X,2))==0yval1(indexf(n),1)=yL;elseyval1(indexf(n),1)=y0;endendyval1=reshape(yval1,size(X,2),size(yval1,1)/size(X,2));yval1=yval1';end%CD solutiond=zeros(size(tt,2),mdim);a=repmat([1],size(tt,2),mdim);for n=1:size(tt,2)t=tt(1,n);b(n,:)=repmat([0.5*(1-0.5*t)],1,mdim);c(n,:)=repmat([0.5*(1+0.5*t)],1,mdim);d(n,1)=0.5*(1+0.5*tt(1,n))*y0;d(n,mdim)=0.5*(1-0.5*tt(1,n))*yL;endc(:,1)=0;b(:,mdim)=0;%numerical cal by using TDMA subfuctionyval2=TDMA(a,b,c,d,mdim);yval2=[repmat([1],size(tt,2),1),yval2,repmat([2],size(tt,2),1)]; Fig(1,X,yval1,yval2,tt);title('CD Vs. Exact Solution')% FUS solutiond=zeros(size(tt,2),mdim);a=repmat([1],size(tt,2),mdim);for n=1:size(tt,2)t=tt(1,n);b(n,:)=repmat([1/(2+t)],1,mdim);c(n,:)=repmat([(1+t)/(2+t)],1,mdim);d(n,1)=(1+tt(1,n))/(2+tt(1,n))*y0;d(n,mdim)=1/(2+tt(1,n))*yL;endc(:,1)=0;b(:,mdim)=0;%numerical cal by using TDMA subfuctionyval3=TDMA(a,b,c,d,mdim);yval3=[repmat([1],size(tt,2),1),yval3,repmat([2],size(tt,2),1)]; Fig(2,X,yval1,yval3,tt);title('FUS Vs. Exact Solution')% HS solutiond=zeros(size(tt,2),mdim);a=repmat([1],size(tt,2),mdim);for n=1:size(tt,2)t=tt(1,n);if t>2b(n,:)=repmat([0],1,mdim);c(n,:)=repmat([1],1,mdim);d(n,1)=y0;elseif t<-2b(n,:)=repmat([1],1,mdim);c(n,:)=repmat([0],1,mdim);d(n,mdim)=yL;elseb(n,:)=repmat([0.5*(1-0.5*t)],1,mdim);c(n,:)=repmat([0.5*(1+0.5*t)],1,mdim);d(n,1)=0.5*(1+0.5*t)*y0;d(n,mdim)=0.5*(1-0.5*t)*yL;endendc(:,1)=0;b(:,mdim)=0;% numerical cal by using TDMA subfuctionyval4=TDMA(a,b,c,d,mdim);yval4=[repmat([1],size(tt,2),1),yval4,repmat([2],size(tt,2),1)]; Fig(3,X,yval1,yval4,tt);title('HS Vs. Exact Solution')%QUICK Solutiond=zeros(size(tt,2),mdim);a=repmat([1],size(tt,2),mdim);for n=1:size(tt,2)t=tt(1,n);b(n,:)=repmat([1/(2+t)],1,mdim);c(n,:)=repmat([(1+t)/(2+t)],1,mdim);d(n,1)=(1+tt(1,n))/(2+tt(1,n))*y0;d(n,mdim)=1/(2+tt(1,n))*yL;endc(:,1)=0;b(:,mdim)=0;%numerical cal by using TDMA subfuctionyval5=zeros(size(tt,2),mdim);yval5com=yval5+1;counter=1;%iterativewhile max(max(abs(yval5-yval5com)))>10^-10if counter==1yval5com=TDMA(a,b,c,d,mdim);endfor nn=1:size(tt,2)for nnn=1:mdimif nnn==1d(nn,nnn)=((6*yval5com(nn,nnn)-3*y0-3*yval5com(nn,nnn+1))*tt(1,nn))/(8*(2+tt(1, nn)))+((1+tt(1,nn))/(2+tt(1,nn))*y0);elseif nnn==2d(nn,nnn)=((5*yval5com(nn,nnn)-3*yval5com(nn,nnn+1)-yval5com(nn,nnn-1)-y0)*tt (1,nn))/(8*(2+tt(1,nn)));elseif nnn==mdimd(nn,nnn)=((5*yval5com(nn,nnn)-3*yL-yval5com(nn,nnn-1)-yval5com(nn,nnn-2))*tt (1,nn))/(8*(2+tt(1,nn)))+(1/(2+tt(1,nn))*yL);elsed(nn,nnn)=((5*yval5com(nn,nnn)-3*yval5com(nn,nnn+1)-yval5com(nn,nnn-1)-yval5 com(nn,nnn-2))*tt(1,nn))/(8*(2+tt(1,nn)));endendendyval5=TDMA(a,b,c,d,mdim);temp=yval5;yval5=yval5com;yval5com=temp;counter=counter+1;endyval5=yval5com;yval5=[repmat([1],size(tt,2),1),yval5,repmat([2],size(tt,2),1)];Fig(4,X,yval1,yval5,tt);title('QUICK Vs. Exact Solution')%-------------TDMA SubFunction------------------function y=TDMA(a,b,c,d,mdim)%form a b c d resolve yval2 by using TDMA%eliminationp(:,1)=b(:,1)./a(:,1);q(:,1)=d(:,1)./a(:,1);for n=2:mdimp(:,n)=b(:,n)./(a(:,n)-c(:,n).*p(:,n-1));q(:,n)=(d(:,n)+c(:,n).*q(:,n-1))./(a(:,n)-c(:,n).*p(:,n-1));end%iterativey(:,mdim)=q(:,mdim);for n=(mdim-1):-1:1y(:,n)=p(:,n).*y(:,n+1)+q(:,n);end%-------------ResultCom SubFunction------------------ function y=ResultCom (a,b,c)for n=1:max(size(c,2))y(2*n-1,:)=a(n,:);y(2*n,:)=b(n,:);end%-------------Fig SubFunction------------------ function y=Fig(n,a,b,c,d)figure(n);plot(a,b);hold onplot(a,c,'*');str='''legend(';for n=1:size(d,2)if n==size(d,2)str=strcat(str,'''''Pe=',num2str(d(1,n)),''''')''');elsestr=strcat(str,'''''Pe=',num2str(d(1,n)),''''',');endendeval(eval(str));精确解与数值解的对比图，其中边界条件给定10=φ，2=L φ。
(完整版)数值传热学陶文铨主编第二版习题答案
数值传热学4-9章习题答案习题4－2一维稳态导热问题的控制方程：022=+∂∂S xTλ依据本题给定条件，对节点2节点3采用第三类边界条件具有二阶精度的差分格式，最后得到各节点的离散方程：节点1：1001=T 节点2：1505105321-=+-T T T 节点3：75432=+-T T 求解结果：，852=T 403=T 对整个控制容积作能量平衡，有：2150)4020(15)(3=⨯+-⨯=∆+-=∆+x S T T h x S q f f B 即：计算区域总体守恒要求满足习题4－5在4－2习题中，如果，则各节点离散方程如下：25.03)(10f T T h -⨯=节点1：1001=T 节点2：1505105321-=+-T T T 节点3：25.03325.032)20(4015])20(21[-⨯+=-⨯++-T T T T 对于节点3中的相关项作局部线性化处理，然后迭代计算；求解结果：，（迭代精度为10-4）818.822=T 635.353=T 迭代计算的Matlab 程序如下：x=30;x1=20;while abs(x1-x)>0.0001a=[1 0 0;5 -10 5;0 -1 1+2*(x-20)^(0.25)]; b=[100;-150; 15+40*(x-20)^(0.25)]; t=a^(-1)*b;x1=x;x=t(3,1);endtcal=t习题4-12的Matlab程序%代数方程形式A i T i=C i T i+1+B i T i-1+D imdim=10;%计算的节点数x=linspace(1,3,mdim);%生成A、C、B、T数据的基数；A=cos(x);%TDMA的主对角元素B=sin(x);%TDMA的下对角线元素C=cos(x)+exp(x); %TDMA的上对角线元素T=exp(x).*cos(x); %温度数据%由A、B、C构成TDMAcoematrix=eye(mdim,mdim);for n=1:mdimcoematrix(n,n)=A(1,n);if n>=2coematrix(n,n-1)=-1*B(1,n);endif n<mdimcoematrix(n,n+1)=-1*C(1,n);endend%计算D矢量D=(coematrix*T')';%由已知的A、B、C、D用TDMA方法求解T%消元P(1,1)=C(1,1)/A(1,1);Q(1,1)=D(1,1)/A(1,1);for n=2:mdimP(1,n)=C(1,n)/(A(1,n)-B(1,n)*P(1,n-1));Q(1,n)=(D(1,n)+B(1,n)*Q(1,n-1))/(A(1,n)-B(1,n)*P(1,n-1)); end%回迭Tcal(1,mdim)=Q(1,mdim);for n=(mdim-1):-1:1Tcal(1,n)=P(1,n)*Tcal(1,n+1)+Q(1,n);endTcom=[T;Tcal];%绘图比较给定T值和计算T值plot(Tcal,'r*')hold onplot(T)n gin th a r e 结果比较如下，由比较可知两者值非常切合（在小数点后8位之后才有区别）：习题4-14充分发展区的温度控制方程如下：)(1rTr r r x T uc p ∂∂∂∂=∂∂λρ对于三种无量纲定义、、进行分析如下w b w T T T T --=Θ∞∞--=ΘT T T T w ww T T T T --=Θ∞1）由得：wb wT T T T --=Θww b T T T T +Θ-=)(由可得：T x T x T x T T T x T w b w w b ∂∂Θ-+∂∂Θ=∂+Θ-∂=∂∂)1(])[(rT r T T r T T T r T w w b w w b ∂∂Θ-+∂Θ∂-=∂+Θ-∂=∂∂)1()(])[(由与无关、与无关以及、的表达式可知，除了均匀的情况外，该无量b T r Θx x T ∂∂rT∂∂w T 纲温度定义在一般情况下是不能用分离变量法的；2）由得：∞∞--=ΘT T T T w ∞∞+Θ-=T T T T w )(由可得：T xT x T T T x T w w ∂∂Θ=∂+Θ-∂=∂∂∞∞])[(rT r T T r T T T r T w w w ∂∂Θ+∂Θ∂-=∂+Θ-∂=∂∂∞∞∞)(])[(由与无关、与无关以及、的表达式可知，在常见的四种边界条件中除了b T r Θx x T ∂∂rT ∂∂轴向及周向均匀热流的情况外，有，则该无量纲温度定义是可以用分const q w =0=∂∂rT w离变量法的；3）由得：wwT T T T --=Θ∞ww T T T T +Θ-=∞)(由可得：T xT x T T T x T w w w ∂∂Θ-=∂+Θ-∂=∂∂∞)1(])[(r T T r T T T r T w w w -+∂Θ∂-=∂+Θ-∂=∂∂∞∞1()(])[(同2）分析可知，除了轴向及周向均匀热流const q w =温度定义是可以用分离变量法的；习题4－181）采用柱坐标分析，写出统一的稳态柱坐标形式动量方程：S r r r r r r x x w r v r r r u x +∂∂∂∂+∂∂∂∂+∂∂∂∂=∂∂+∂∂+∂∂(1)(1)()(1)(1)(θφλθφλφλφρθφρφρ、和分别是圆柱坐标的3个坐标轴，、和分别是其对应的速度分量，其中x r θu v w 是管内的流动方向；x 对于管内的层流充分发展有：、，；0=v 0=w 0=∂∂xu并且方向的源项：x x pS ∂∂-=方向的源项：r r pS ∂∂-=方向的源项：θθ∂∂-=pr S 1由以上分析可得到圆柱坐标下的动量方程：方向：x 0)(1)(1=∂∂-∂∂∂∂+∂∂∂∂x pu r r r u r r r θλθλ方向：r 0=∂∂r p 方向：θ0=∂∂θp 边界条件：，R r =0=u ，；对称线上，0=r 0=∂∂r u 0=∂∂θu 不考虑液体的轴向导热，并简化分析可以得到充分发展的能量方程为：)(1(1θλθλρ∂∂∂∂+∂∂∂∂=∂∂Tr r r T r r r x T uc p 边界条件：，；，R r =w q r T =∂∂λ0=r 0=∂∂rT，πθ/0=0=∂∂-θλT2）定义无量纲流速：dxdp R uU 2-=λ并定义无量纲半径：；将无量纲流速和无量纲半径代入方向的动量方程得：R r /=ηx 0))1((1)1((122=∂∂-∂-∂∂∂+∂-∂∂∂xp U dx dp R R R R U dx dp R RR R θληλθηηλληηη上式化简得：011(1(1=+∂∂∂∂+∂∂∂∂θηθηηηηηU U 边界条件：，1=η0=U ，；对称线上，0=η0=∂∂ηU 0=∂∂θU定义无量纲温度：λ/0R q T T b-=Θ其中，是折算到管壁表面上的平均热流密度，即：；0q Rq q wπ=0由无量纲温度定义可得：bT Rq T +Θ=λ0将表达式和无量纲半径代入能量方程得：T η(1)(100θληλθηηλληηηρ∂Θ∂∂∂+∂Θ∂∂∂=∂∂R q R R R R q R R R x T uc b p 化简得：（1）)1(1)(10θηθηηηηηρ∂Θ∂∂∂+∂Θ∂∂∂=∂∂x T u c q R b p 由热平衡条件关系可以得：mm m b m p b p p RU U q R u u R q A u u dx dT A u c x T u c x T uc 020221221)(===∂∂=∂∂ππρρρ将上式代入式（1）可得：)1(1)(12θηθηηηηη∂Θ∂∂∂+∂Θ∂∂∂=m U U 边界条件：，；，0=η0=∂Θ∂η1=ηR q q w πη10==∂Θ∂，；，0=θ0=∂Θ∂θπθ=0=∂Θ∂θ单值条件：由定义可知： 且： 0/0=-=ΘλR q T T b b b ⎰⎰Θ=ΘAAb UdAUdA 即得单值性条件：=Θ⎰⎰AA UdAUdA 3）由阻力系数及定义有：f Re 228)(21/Re ⎪⎭⎫ ⎝⎛=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-=D D U D u u dx dp D f e m e m me νρ且：m W b m W b m W R q T T D T T q Nu ,0,,0~2)/(2Θ=-=-=λλ5－21．一维稳态无源项的对流－扩散方程如下所示： （取常物性）xx u 22∂∂Γ=∂∂φφρ边界条件如下：LL x x φφφφ====,;,00上述方程的精确解如下： 11)/(00--=--⋅Pe L x Pe L e e φφφφΓ=/uL Pe ρ2．将分成20等份，所以有：L ∆=P Pe 20 1 2 3 4 5 6……………………… 17 18 19 20 21对于中心差分、一阶迎风、混合格式和QUICK 格式分别分析如下：1）中心差分中间节点： 2)5.01()5.01(11-∆+∆++-=i i i P P φφφ20,2 =i 2）一阶迎风中间节点： ∆-∆++++=P P i i i 2)1(11φφφ20,2 =i 3）混合格式当时，中间节点： 1=∆P 2)5.01()5.01(11-∆+∆++-=i i i P P φφφ 20,2 =i 当时，中间节点： 10,5=∆P 1-=i i φφ20,2 =i 4）QUICK 格式*12111)35(8122121⎥⎦⎤⎢⎣⎡---++++++=+--∆∆-∆∆+∆i i i i i i i P P P P P φφφφφφφ2≠i*1111)336(8122121⎥⎦⎤⎢⎣⎡--++++++=+-∆∆-∆∆+∆i i i i i i P P P P P φφφφφφ2=i 数值计算结果与精确解的计算程序如下：%except for HS, any other scheme doesnt take Pe<0 into consideration %expression of exact solutiony=dsolve('a*b*Dy=c*D2y','y(0)=y0,y(L)=yL','x')y=subs(y,'L*a*b/c','t')y=simple(subs(y,'a*b/c*x','t*X'));ysim=simple(sym(strcat('(',char(y),'-y0)','/(yL-y0)')))y=sym(strcat('(',char(ysim),')*(yL-y0)','+y0'))% in the case of Pe=0y1=dsolve('D2y=0','y(0)=y0,y(L)=yL','x')y1=subs(y1,'-(y0-yL)/L*x','(-y0+yL)*X')%grid Pe number tt=[1 5 10];%dimensionless length m=20;%mdim is the number of inner node mdim=m-1;X=linspace(0,1,m+1);%initial value of variable during calculation y0=1;yL=2;%cal exact solution for n=1:size(tt,2) t=m*tt(1,n); if t==0 yval1(n,:)=eval(y1); else yval1(n,:)=eval(y); end end%extra treatment because max number in MATLAB is 10^308if max(isnan(yval1(:))) yval1=yval1'; yval1=yval1(:);indexf=find(isnan(yval1)); for n=1:size(indexf,1) if rem(indexf(n,1),size(X,2))==0 yval1(indexf(n),1)=yL; else yval1(indexf(n),1)=y0; endendyval1=reshape(yval1,size(X,2),size(yval1,1)/size(X,2));yval1=yval1';end%CD solutiond=zeros(size(tt,2),mdim);a=repmat([1],size(tt,2),mdim);for n=1:size(tt,2)t=tt(1,n);b(n,:)=repmat([0.5*(1-0.5*t)],1,mdim);c(n,:)=repmat([0.5*(1+0.5*t)],1,mdim);d(n,1)=0.5*(1+0.5*tt(1,n))*y0;d(n,mdim)=0.5*(1-0.5*tt(1,n))*yL;endc(:,1)=0;b(:,mdim)=0;%numerical cal by using TDMA subfuctionyval2=TDMA(a,b,c,d,mdim);yval2=[repmat([1],size(tt,2),1),yval2,repmat([2],size(tt,2),1)]; Fig(1,X,yval1,yval2,tt);title('CD Vs. Exact Solution')% FUS solutiond=zeros(size(tt,2),mdim);a=repmat([1],size(tt,2),mdim);for n=1:size(tt,2)t=tt(1,n);b(n,:)=repmat([1/(2+t)],1,mdim);c(n,:)=repmat([(1+t)/(2+t)],1,mdim);d(n,1)=(1+tt(1,n))/(2+tt(1,n))*y0;d(n,mdim)=1/(2+tt(1,n))*yL;endc(:,1)=0;b(:,mdim)=0;%numerical cal by using TDMA subfuctionyval3=TDMA(a,b,c,d,mdim);yval3=[repmat([1],size(tt,2),1),yval3,repmat([2],size(tt,2),1)]; Fig(2,X,yval1,yval3,tt);title('FUS Vs. Exact Solution')% HS solutiond=zeros(size(tt,2),mdim);a=repmat([1],size(tt,2),mdim);for n=1:size(tt,2)t=tt(1,n);if t>2b(n,:)=repmat([0],1,mdim);c(n,:)=repmat([1],1,mdim);d(n,1)=y0;elseif t<-2b(n,:)=repmat([1],1,mdim);c(n,:)=repmat([0],1,mdim);d(n,mdim)=yL;elseb(n,:)=repmat([0.5*(1-0.5*t)],1,mdim);c(n,:)=repmat([0.5*(1+0.5*t)],1,mdim);d(n,1)=0.5*(1+0.5*t)*y0;d(n,mdim)=0.5*(1-0.5*t)*yL;endendc(:,1)=0;b(:,mdim)=0;% numerical cal by using TDMA subfuctionyval4=TDMA(a,b,c,d,mdim);yval4=[repmat([1],size(tt,2),1),yval4,repmat([2],size(tt,2),1)]; Fig(3,X,yval1,yval4,tt);title('HS Vs. Exact Solution')%QUICK Solutiond=zeros(size(tt,2),mdim);a=repmat([1],size(tt,2),mdim);for n=1:size(tt,2)t=tt(1,n);b(n,:)=repmat([1/(2+t)],1,mdim);c(n,:)=repmat([(1+t)/(2+t)],1,mdim);d(n,1)=(1+tt(1,n))/(2+tt(1,n))*y0;d(n,mdim)=1/(2+tt(1,n))*yL;endc(:,1)=0;b(:,mdim)=0;%numerical cal by using TDMA subfuctionyval5=zeros(size(tt,2),mdim);yval5com=yval5+1;counter=1;%iterativewhile max(max(abs(yval5-yval5com)))>10^-10if counter==1yval5com=TDMA(a,b,c,d,mdim);endfor nn=1:size(tt,2)for nnn=1:mdimif nnn==1d(nn,nnn)=((6*yval5com(nn,nnn)-3*y0-3*yval5com(nn,nnn+1))*tt(1,nn))/(8*(2+tt(1,nn)))+((1+tt(1,nn))/(2+tt(1,nn))*y0);elseif nnn==2d(nn,nnn)=((5*yval5com(nn,nnn)-3*yval5com(nn,nnn+1)-yval5com(nn,nnn-1)-y0)*tt(1,nn))/(8*(2+tt(1,nn)));elseif nnn==mdimd(nn,nnn)=((5*yval5com(nn,nnn)-3*yL-yval5com(nn,nnn-1)-yval5com(nn,nnn-2))*tt(1,nn))/(8*(2+tt(1,nn)))+(1/(2+tt(1,nn))*yL);elsed(nn,nnn)=((5*yval5com(nn,nnn)-3*yval5com(nn,nnn+1)-yval5com(nn,nnn-1)-yval5com(nn,nnn-2))*tt(1,nn))/(8*(2+tt(1,nn)));endendendyval5=TDMA(a,b,c,d,mdim);temp=yval5;yval5=yval5com;yval5com=temp;counter=counter+1;endyval5=yval5com;yval5=[repmat([1],size(tt,2),1),yval5,repmat([2],size(tt,2),1)];Fig(4,X,yval1,yval5,tt);title('QUICK Vs. Exact Solution')%-------------TDMA SubFunction------------------function y=TDMA(a,b,c,d,mdim)%form a b c d resolve yval2 by using TDMA%eliminationp(:,1)=b(:,1)./a(:,1);q(:,1)=d(:,1)./a(:,1);for n=2:mdimp(:,n)=b(:,n)./(a(:,n)-c(:,n).*p(:,n-1));q(:,n)=(d(:,n)+c(:,n).*q(:,n-1))./(a(:,n)-c(:,n).*p(:,n-1));end%iterativey(:,mdim)=q(:,mdim);for n=(mdim-1):-1:1y(:,n)=p(:,n).*y(:,n+1)+q(:,n);end%-------------ResultCom SubFunction------------------function y=ResultCom (a,b,c)for n=1:max(size(c,2))y(2*n-1,:)=a(n,:);y(2*n,:)=b(n,:);end%-------------Fig SubFunction------------------function y=Fig(n,a,b,c,d)figure(n);plot(a,b);hold onplot(a,c,'*');str='''legend(';for n=1:size(d,2)if n==size(d,2)str=strcat(str,'''''Pe=',num2str(d(1,n)),''''')''');elsestr=strcat(str,'''''Pe=',num2str(d(1,n)),''''',');endendeval(eval(str));a n d A l l t h i n g s i n t h ei r b e i n g a r e g 13精确解与数值解的对比图，其中边界条件给定，。
数值传热学第二章部分习题参考答案
习题2－4 [解]1．先用控制容积积分法得出离散方程： 以r 乘式01=+⎪⎭⎫⎝⎛S dr dT rk dr d r ，并对图2-2所示的控制容积P 作积分： wewe dr dT rk dr dT rk dr dr dT rk dr d r r⎪⎭⎫⎝⎛-⎪⎭⎫ ⎝⎛=⎪⎭⎫ ⎝⎛⎰1 2-4-1 ()E P e eT T dT dr r δ-⎛⎫=⎪⎝⎭ 2-4-1-1 ()P W w wT T dT dr r δ-⎛⎫= ⎪⎝⎭2-4-1-2将式（2-4-1-1）、式（2-4-1-2）代入式（2-4-1）可以得到：()()W P wP E e ewT T x rk T T x rk dr dr dT rk dr d r r-⎪⎭⎫⎝⎛--⎪⎭⎫ ⎝⎛=⎪⎭⎫ ⎝⎛⎰δδ1 2-4-2 222ee wP w r r rSdr S Sr r -==∆⎰2-4-3根据式（2-4-2）、式（2-4-3）可以得到：P E W P e w e w rk rk rk rk T T T Sr r x x x x δδδδ⎡⎤⎛⎫⎛⎫⎛⎫⎛⎫+=++∆ ⎪ ⎪ ⎪ ⎪⎢⎥⎝⎭⎝⎭⎝⎭⎝⎭⎣⎦2-4-4令e E x rk a ⎪⎭⎫ ⎝⎛=δ，wW x rk a ⎪⎭⎫⎝⎛=δ，W E P a a a +=，P b Sr r =∆，式（2-4-4）可以写成b T a T a T a W W E E P P ++=的形式。
2. 再用Taylor 展开法导出022=++S drdTr k dr T d k的离散方程。
将点E T 对点P T 作Taylor 展开，有：()() +++=!2222e Pe P P E x dr Td x drdTT T δδ 2-4-5再将点W T 对点P T 作Taylor 展开，有：()() ++-=!2222wPw P P W x dr Td x drdTT T δδ 2-4-6根据式（2-4-5）、式（2-4-6）可以计算出dr dT ，22drTd ()[]()()[]()[]()()[]()()[]222222w e e w We P w e E w x x x x T x T x x T x dr dT δδδδδδδδ+---= 2-4-7 ()()()[]()()[]()()[]()222222we e w We P w e E w x x x x T x T x x T x dr T d δδδδδδδδ+++-= 2-4-8 将式（2-4-7）、式（2-4-8）代入上面的非守恒型方程，整理成（并考虑到常物性、均分网格）：222P P P P E W P kr kr kr k k T T T r rS r r r ⎡⎤⎡⎤=++-+∆⎢⎥⎢⎥∆∆∆⎣⎦⎣⎦2-4-9 令12e P E kr ra k r r ⎡⎤=+=⎢⎥∆∆⎣⎦，12w P W kr r a k r r⎡⎤=-=⎢⎥∆∆⎣⎦，W E P a a a +=，P b r rS =∆式（2-4-9）也可以写成b T a T a T a W W E E P P ++=的形式。
数值传热学部分习题问题详解2
习题4－2一维稳态导热问题的控制方程：022=+∂∂S xTλ 依据本题给定条件，对节点2采用二阶精度的中心差分格式，节点3采用第三类边界条件具有二阶精度的差分格式，最后得到各节点的离散方程： 节点1： 1001=T节点2： 1505105321-=+-T T T 节点3：75432=+-T T求解结果：852=T ，403=T对整个控制容积作能量平衡，有：02150)4020(15)(3=⨯--⨯=∆+-=∆+x S T T h x S q f f B即：计算区域总体守恒要求满足习题4－5在4－2习题中，如果25.03)(10f T T h -⨯=，则各节点离散方程如下：节点1： 1001=T节点2： 1505105321-=+-T T T节点3：25.03325.032)20(4015])20(21[-⨯+=-⨯++-T T T T对于节点3中的相关项作局部线性化处理，然后迭代计算； 求解结果：818.822=T ，635.353=T （迭代精度为10-4）迭代计算的Matlab 程序如下：x=30; x1=20;while abs(x1-x)>0.0001a=[1 0 0;5 -10 5;0 -1 1+2*(x-20)^(0.25)]; b=[100;-150; 15+40*(x-20)^(0.25)]; t=a^(-1)*b; x1=x; x=t(3,1); end tcal=t习题4-14充分发展区的温度控制方程如下：)(1rTr r r x T uc p ∂∂∂∂=∂∂λρ 对于三种无量纲定义w b w T T T T --=Θ、∞∞--=ΘT T T T w 、w w T T T T --=Θ∞进行分析如下1）由wb wT T T T --=Θ得：w w b T T T T +Θ-=)(由T 可得：x T x T x T T T x T w b w w b ∂∂Θ-+∂∂Θ=∂+Θ-∂=∂∂)1(])[(rT r T T r T T T r T w w b w w b ∂∂Θ-+∂Θ∂-=∂+Θ-∂=∂∂)1()(])[( 由b T 与r 无关、Θ与x 无关以及x T ∂∂、rT∂∂的表达式可知，除了w T 均匀的情况外，该无量纲温度定义在一般情况下是不能用分离变量法的； 2）由∞∞--=ΘT T T T w 得：∞∞+Θ-=T T T T w )(由T 可得：xT x T T T x T w w ∂∂Θ=∂+Θ-∂=∂∂∞∞])[(rT r T T r T T T r T w w w ∂∂Θ+∂Θ∂-=∂+Θ-∂=∂∂∞∞∞)(])[( 由b T 与r 无关、Θ与x 无关以及x T ∂∂、rT∂∂的表达式可知，在常见的四种边界条件中除了轴向及周向均匀热流const q w =的情况外，有0=∂∂rT w离变量法的；3）由wwT T T T --=Θ∞得：w w T T T T +Θ-=∞)(由T 可得： xT x T T T x T w w w ∂∂Θ-=∂+Θ-∂=∂∂∞)1(])[(rT r T T r T T T r T w w w w ∂∂Θ-+∂Θ∂-=∂+Θ-∂=∂∂∞∞)1()(])[( 同2）分析可知，除了轴向及周向均匀热流const q w =的情况外，有0=∂∂rT w，该无量纲温度定义是可以用分离变量法的；习题4－181）采用柱坐标分析，写出统一的稳态柱坐标形式动量方程：S r r r r r r x x w r v r r r u x +∂∂∂∂+∂∂∂∂+∂∂∂∂=∂∂+∂∂+∂∂)(1)(1)()(1)(1)(θφλθφλφλφρθφρφρ x 、r 和θ分别是圆柱坐标的3个坐标轴，u 、v 和w 分别是其对应的速度分量，其中x 是管内的流动方向；对于管内的层流充分发展有：0=v 、0=w ，0=∂∂xu；并且x 方向的源项：x p S ∂∂-= r 方向的源项：r pS ∂∂-=θ方向的源项：θ∂∂-=pr S 1 由以上分析可得到圆柱坐标下的动量方程： x 方向： 0)(1)(1=∂∂-∂∂∂∂+∂∂∂∂x pu r r r u r r r θλθλ r 方向：0=∂∂r pθ方向：0=∂∂θp边界条件： R r =，0=u0=r ，0=∂∂r u ；对称线上，0=∂∂θu 不考虑液体的轴向导热，并简化分析可以得到充分发展的能量方程为：)(1)(1θλθλρ∂∂∂∂+∂∂∂∂=∂∂Tr r r T r r r x T uc p 边界条件： R r =，w q r T =∂∂λ；0=r ，0=∂∂rTπθ/0=，0=∂∂-θλT2）定义无量纲流速：dxdp R uU 2-=λ并定义无量纲半径：R r /=η；将无量纲流速和无量纲半径代入x 方向的动量方程得：0))1((1))1((122=∂∂-∂-∂∂∂+∂-∂∂∂xp U dx dp R R R R U dx dp R RR R θληλθηηλληηη 上式化简得：01)1(1)(1=+∂∂∂∂+∂∂∂∂θηθηηηηηU U 边界条件：1=η，0=U0=η，0=∂∂ηU ；对称线上，0=∂∂θU定义无量纲温度：λ/0R q T T b-=Θ其中，0q 是折算到管壁表面上的平均热流密度，即：Rq q wπ=0； 由无量纲温度定义可得： b T Rq T +Θ=λ将T 表达式和无量纲半径η代入能量方程得：)(1)(100θληλθηηλληηηρ∂Θ∂∂∂+∂Θ∂∂∂=∂∂R q R R R R q R R R x T uc b p 化简得：)1(1)(10θηθηηηηηρ∂Θ∂∂∂+∂Θ∂∂∂=∂∂x T u c q R b p （1）由热平衡条件关系可以得：mm m b m p b p p RU U q R u u R q A u u dx dT A u c x T u c x T uc 020221221)(===∂∂=∂∂ππρρρ 将上式代入式（1）可得：)1(1)(12θηθηηηηη∂Θ∂∂∂+∂Θ∂∂∂=m U U 边界条件： 0=η，0=∂Θ∂η；1=η，R q q w πη10==∂Θ∂0=θ，0=∂Θ∂θ；πθ=，0=∂Θ∂θ单值条件：由定义可知：0/0=-=ΘλR q T T b b b 且： ⎰⎰Θ=ΘAAb UdAUdA即得单值性条件：0=Θ⎰⎰AA UdAUdA 3）由阻力系数f 及Re 定义有：228)(21/Re ⎪⎭⎫ ⎝⎛=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-=D D U D u u dx dp D f e m e m me νρ 且：m W b m W b m W R q T T D T T q Nu ,0,,0~2)/(2Θ=-=-=λλ5－21．一维稳态无源项的对流－扩散方程如下所示：xx u 22∂∂Γ=∂∂φφρ （取常物性）边界条件如下：L L x x φφφφ====,;,00上述方程的精确解如下：11)/(00--=--⋅PeL x Pe L e e φφφφ Γ=/uL Pe ρ 2．将L 分成20等份，所以有：∆=Pe 20图示如下：1 2 3 4 5 6 ………… …………… 17 18 19 20 21 对于中心差分、一阶迎风、混合格式和QUICK 格式分别分析如下： 1） 中心差分中间节点： 2)5.01()5.01(11-∆+∆++-=i i i P P φφφ 20,2Λ=i2） 一阶迎风中间节点： ∆-∆++++=P P i i i 2)1(11φφφ 20,2Λ=i3） 混合格式当1=∆P 时，中间节点：2)5.01()5.01(11-∆+∆++-=i i i P P φφφ20,2Λ=i当10,5=∆P 时，中间节点： 1-=i i φφ 20,2Λ=i 4） QUICK 格式*12111)35(8122121⎥⎦⎤⎢⎣⎡---++++++=+--∆∆-∆∆+∆i i i i i i i P P P P P φφφφφφφ 2≠i *1111)336(8122121⎥⎦⎤⎢⎣⎡--++++++=+-∆∆-∆∆+∆i i i i i i P P P P P φφφφφφ 2=i 5－3乘方格式：⎪⎪⎩⎪⎪⎨⎧<-≤≤--+≤≤->=∆∆∆∆∆∆∆∆10,010,)1.01(100,)1.01(10,055P P P P P P P P D a e E当1.0=∆P 时有：951.0)1.01.01()1.01(55=⨯-=-=∆P D a eE因为：301.0/3)()()()()()(===Γ=Γ=∆eee e e e e e e P u x u u x D ρδρρδ 所以：5297.2830951.0951.0=⨯==e E D a由系数关系式∆=-P D a D a eEw W 可得： 53.3130)951.01.0()(=⨯+=⨯+=∆w eEW D D a P a 且： 205.01.010=⨯=∆∆=txa P p ρ 当采用隐式时1=f ，因此可得：0597.62253.315297.280=++=++=P W E P a fa fa a同理可得当10=∆P 时有：0=E a ，3=W a ，5=P a5－5二维稳态无源项的对流－扩散问题的控制方程：)()()()(yy x x y v x u ∂∂Γ∂∂+∂∂Γ∂∂=∂∂+∂∂φφφρφρφφ 对于一阶迎风、混合、乘方格式的通用离散方程：S S N N W W E E P P a a a a a φφφφφ+++=其中：[]0,)(e e e E F P A D a -+=∆ []0,)(w w w W F P A D a +=∆ []0,)(n n n N F P A D a -+=∆ []0,)(s s s S F P A D a +=∆5－71）QUICK 格式的界面值定义如下：⎪⎪⎩⎪⎪⎨⎧-+=-+=)36(81)36(81WW P W w W E P e φφφφφφφφ0>u 对（5－1）式dxdx d d dx u d )()(φφρΓ=积分可得： w e w e dxd dx d u u )()()()(φφφρφρΓ-Γ=-对流项采用QUICK 格式的界面插值，扩散项采用线性界面插值，对于0>u 及均分网格有：)]()([]))(36())(36[(81x x u u W P w P E e w WW P W e W E P ∆-Γ-∆-Γ=-+--+φφφφρφφφρφφφ 整理得：WW w W w e w E e e P w e w e u u u x u x x x u u φρφρρφρφρρ)(81])(43)(81[])(83[)]()(83)(43[-++∆Γ+-∆Γ=∆Γ+∆Γ+-上式即为QUICK 格式离散得到的离散方程；2）要分析QUICK 格式的稳定性，则应考虑非稳平流方程：xut ∂∂-=∂∂φφ 在t ∆时间间隔内对控制容积作积分：⎰⎰⎰⎰∆+∆+∂∂-=∂∂t t t e w e w tt tdxdt x u dtdx xφφ得：dt u dx tt tw e e wttt ⎰⎰∆+∆+--=-)()(φφφφφ随时间变化采用阶梯显式，随空间变化采用QUICK 格式得：t u x WW P W W E P tP t t P ∆+---+-=∆-∆+)]3636(81[)(φφφφφφφφ整理得：xu t ni n i n i n i ni n i ∆+-+-∆---++87332111φφφφφφ对于初始均匀零场，假设在),(n i 点有一个扰动n i ε； 对1+i 点写出QUICK 格式的离散方程：xu tni n i n i n i n i n i ∆+-+-∆--+++++8733121111φφφφφφ可得：ni n i xt u εφ∆∆=++8711 对1-i 点分析可得：ni n i xt u εφ∆∆-=+-8311 由于扩散对扰动的传递恒为正，其值为ni x t ερ2∆Γ∆，所以根据符号不变原则有： 0)/)83(2≥∆Γ∆+∆∆-ni n i n i xt x t u εερε 整理得到QUICK 格式的稳定性条件为：38≤∆P 5－91）三阶迎风格式采用上游两个节点和下游一个节点的值来构造函数界面插值形式，所以定义如下：⎩⎨⎧<++=>++=00u c b a u c b a EEE P e W P E e φφφφφφφφ根据上述定义，在0>u 时对控制容积内的对流项作积分平均可得：])()([1)(11WW W P E e w w e c b c a b a xxdx x x φφφφφφφ--+-+∆=-∆=∂∂∆⎰由表2－1式可知三阶迎风格式的差分格式：xxni n i n i n i ni ∆+-+=∂∂--+1221264211,φφφφφ 由控制容积积分法得到的对流项离散格式应与Taylor 离散展开得到的离散格式具有相同的形式和精度，所以比较可得：61,65,31-===c b a所以三阶迎风格式的函数插值定义为：⎪⎪⎩⎪⎪⎨⎧<-+=>-+=06165310616531u u EE E P e W P E e φφφφφφφφ2）由上述分析可知，得到的三阶迎风格式的插值定义与给出节点上导数表达式的定义在形式上显然是一致的；6－1二维直角坐标中不可压缩流体的连续方程及动量方程如下：⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎨⎧+∂∂∂∂+∂∂∂∂+∂∂-=∂∂+∂∂+∂∂+∂∂∂∂+∂∂∂∂+∂∂-=∂∂+∂∂+∂∂=∂∂+∂∂)3()()()()()()2()()()()()()1(0vu S y y v x x v y p y vv x vu t v S y y u x x u x p y uv x uu tu y v x u ηηρρρηηρρρ假设常粘性，则0==v u S S ；对公式（2）及（3）分别对y x ,求偏导得：⎪⎪⎩⎪⎪⎨⎧∂∂+⎪⎪⎭⎫ ⎝⎛∂∂∂∂+⎪⎪⎭⎫ ⎝⎛∂∂∂∂-=⎪⎪⎭⎫ ⎝⎛∂∂∂∂+⎪⎭⎫ ⎝⎛∂∂∂∂+⎪⎭⎫ ⎝⎛∂∂∂∂⎪⎪⎭⎫⎝⎛∂∂∂∂+∂∂+⎪⎭⎫ ⎝⎛∂∂∂∂-=⎪⎪⎭⎫ ⎝⎛∂∂∂∂+⎪⎭⎫ ⎝⎛∂∂∂∂+⎪⎭⎫ ⎝⎛∂∂∂∂33222233)()()()()()(y v x v y y p y y vv y x vu y t v y y u x x u x p x y uv x x uu x t u x ηηρρρηηρρρ 两式相加得并变换积分顺序有：⎥⎦⎤⎢⎣⎡⎪⎪⎭⎫ ⎝⎛∂∂+∂∂∂∂+⎥⎦⎤⎢⎣⎡⎪⎪⎭⎫ ⎝⎛∂∂+∂∂∂∂+⎪⎪⎭⎫⎝⎛∂∂+∂∂-=⎥⎦⎤⎢⎣⎡⎪⎪⎭⎫ ⎝⎛∂∂+∂∂+∂∂∂∂+⎪⎪⎭⎫ ⎝⎛∂∂+∂∂+∂∂∂∂+⎪⎪⎭⎫ ⎝⎛∂∂+∂∂∂∂y v x u yy v x u xy p x p x v u x u v y v v y y v u y u v x u u x y v x u t 2222222222ηηρρ利用连续方程有：⎪⎪⎭⎫ ⎝⎛∂∂+∂∂-=⎥⎦⎤⎢⎣⎡⎪⎪⎭⎫ ⎝⎛∂∂+∂∂∂∂+⎪⎪⎭⎫ ⎝⎛∂∂+∂∂∂∂2222y p x p x v u y v v y y u v x u u x ρ ⎪⎪⎭⎫ ⎝⎛∂∂+∂∂-=⎥⎦⎤⎢⎣⎡∂∂∂∂-∂∂∂∂+∂∂+∂∂+∂∂∂∂22222222222y p xp y v x u y v x u y v x u x v y u ρ 最后即得：⎥⎦⎤⎢⎣⎡∂∂∂∂-∂∂∂∂=⎪⎪⎭⎫ ⎝⎛∂∂+∂∂x v y u y v x u y p x p ρ222226－4假设5*=P p ，则有：5105*-=-=e u 5.3)05(7.0*=-⨯=n v由连续性条件有：s w n e v u v u +=+按SIMPLE 算法有：'''*5)(P E P e e e p p p d u u +-=-+= '''*7.05.3)(P n P n n n p p p d v v +=-+=将上两式代入连续性方程中有：20507.05.35''+=+++-P P p p计算得：06.42'=P p所以：06.4706.425'*=+=+=P P P p p p06.371006.47=-=-=E P e p p u 94.32)006.47(7.0)(7.0=-⨯=-=N P n p p v6－5假设250*3=p ，150*6=p ，所以各点的流量为：⎪⎪⎪⎩⎪⎪⎪⎨⎧-=-⨯==-⨯=-=-⨯=-=-⨯==-⨯=11)15040(1.020)150250(2.024)25010(1.04)270250(2.010)250275(4.0*****E DC B A Q Q Q Q Q 上述流量满足动量方程，但并不满足连续性方程，所以对流量修正：⎪⎪⎪⎩⎪⎪⎪⎨⎧-⨯+-=-⨯+=-⨯+-=-⨯+-=-⨯+=)(1.011)(2.020)(1.024)(2.04)(4.010'6'5'6'3'3'4'2'3'3'1p p Q p p Q p p Q p p Q p p Q ED C B A 对节点3作质量守恒有：B DC A Q Q Q Q +=+即得：)(2.04)(2.020)(1.024)(4.010'2'3'6'3'3'4'3'1p p p p p p p p -⨯+--⨯+=-⨯+--⨯+对节点3作质量守恒有：F E D Q Q Q =+即得：20)(1.011)(2.020'6'5'6'3=-⨯+--⨯+p p p p联立求解上两式有：70.48'3-=p ，13.69'6-=p修正后的压力为：3.20170.48250'3*33=-=+=p p p 87.8013.69150'6*66=-=+=p p p修正后的流量为：⎪⎪⎪⎩⎪⎪⎪⎨⎧-=-⨯==-⨯=-=-⨯=-=-⨯==-⨯=09.4)87.8040(1.009.24)87.803.201(2.013.19)3.20110(1.074.13)2703.201(2.048.29)3.201275(4.0ED C B A Q Q Q Q Q由)(76p p C Q F F -=。
数值传热学答案
例题3-3：算例的matlab程序clear%************一维非稳态热传导方程的数值解法与解析解对比程序**************** %****************************前处理（定义网格与节点）********************x=linspace(0,1,11);t=linspace(0,0.2,101);%****************************前处理（定义边界）************************** % 初值条件for i=1:11if i*0.1<=0.5T(i+1,1)=2*i*0.1;elseT(i+1,1)=2*(1-i*0.1);endend% 边界条件T(1,:)=0;T(11,:)=0;%****************************求解器（差分算法）*************************F=input('please input 网格fourier数：F=')%采用时间一阶向前差分，空间的的二阶中心差分的显式离散格式for j=2:101for i=2:10T(i,j)=F*(T(i+1,j-1)+T(i-1,j-1))+(1-2*F)*T(i,j-1);endend%****************************后处理（结果输出）*************************%输出x和T的关系的二维图形axis([0 1 0 0.6]);plot(x,T([1:11],41),':*r',x,T([1:11],21),':*b',x,T([1:11],11),':*k')hold on%**************************与解析解进行比较**************************%方程的解析解：T(x,t)=(8/pi^2)*exp(-((pi^2)*t))*sin(pi*x)%*******************************************************************axis([0 1 0 0.6]);y1=(8/pi^2)*exp(-(pi^2)*0.192)*sin(pi*x);% t=0.192y2=(8/pi^2)*exp((-pi^2)*0.096)*sin(pi*x);% t=0.096y3=(8/pi^2)*exp((-pi^2)*0.048)*sin(pi*x);% t=0.048xlabel('x');ylabel('T');title('一维非稳态热传导方程的数值解法与解析解');plot(x,y1,'r',x,y2,'b',x,y3,'k')3-9.试证明扩散项的中心差分格式具有守恒性。
传热计算习题附详细答案
传热计算题1.在一内径为0.25cm的管轴心位置上，穿一直径为 0.005cm的细导线，用以测定气体的导热系数。
当导线以0.5A 的电流时，产生的电压降为0.12V/cm，测得导线温度为167℃，空心管内壁温度为150℃。
试求充入管内的气体的导热系数试分析仪器精度以外造成结果误差的客观原因。
2．有两个铜质薄球壳，内球壳外径为0。
015m，外球壳内径为 0.1m，两球壳间装入一种其导热系数待测的粉粒料。
内球用电加热，输入功率为 50w，热量稳定地传向外球，然后散发到周围大气中。
两球壁上都装有热电偶，侧得内球壳的平均温度为120℃，外求壳的平均温度为50℃，周围大气环境温度为20℃；设粉粒料与球壁贴合，试求：（1）待测材料的导热系数（2）外球壁对周围大气的传热系数3．有一面积为10cm2带有保护套的热电偶插入一输送空气的长管内，用来测量空气的温度。
已知热电偶的温度读数为300℃，输气管的壁温为 200℃，空气对保护套的对流传热系数为60w/m2.k，该保护套的黑度为 0.8，试估算由于辐射造成的气体温度测量误差。
并叙述减小测量误差的途径。
已知 Stefan-Bohzman常数σ=5.67×10-9w/m2k 。
4．用两个结构尺寸相同的列管换热器按并联方式加热某中料液。
换热器的管束由32根长 3m 的Ф25×3mm 的钢管组成。
壳程为120℃的饱和蒸汽。
料液总流量为20m3/h，按相等流量分配到两个换热器中作湍流流动，由 25℃加热到 80℃。
蒸汽冷凝对流传热系数为8Kw/m2.℃，管壁及污垢热阻可不记，热损失为零，料液比热为 4.1KJ/kg.℃，密度为 1000kg/m3。
试求：（1）管壁对料液的对流传热系数（2）料液总流量不变，将两个换热器串联，料液加热程度有何变化？（3）此时蒸汽用量有无变化？若有变化为原来的多少倍？（两者情况下蒸汽侧对流传热系数和料液物性不变）5．某厂现有两台单壳程单管程的列管式空气加热器，每台传热面积为A0=20m2（管外面积），均由128根Ф25×2.5mm的钢管组成。
数值传热学习题解答(汇总版)
习题1－7解：由于对称性，取半个通道作为求解区域。
常物性不可压缩流体，二维层流、稳态对流换热的控制方程组为： 质量守恒方程0=∂∂+∂∂yv x u 动量守恒方程 ()()⎪⎪⎭⎫ ⎝⎛∂∂+∂∂+∂∂−=∂∂+∂∂22221y u x u x py vu x uu νρ ()()⎪⎪⎭⎫ ⎝⎛∂∂+∂∂+∂∂−=∂∂+∂∂22221y v x v y p y vv x uv νρ 能量守恒方程 ()()⎪⎪⎭⎫ ⎝⎛∂∂+∂∂=∂∂+∂∂2222y T xT a y vT x uT 边界条件：进口截面 ()0,,===v c T y u u in ； 平板通道上（下）壁面 0,0=∂∂==yTv u ； 中心线上对称条件： 0,0u T v y y∂∂===∂∂； 出口截面0,0,0=∂∂=∂∂=∂∂xT x v x u ； 或者写：采用数值传热学的处理方法。
图1-10 习题1-7的图示本题如果采用整个通道作为计算区域，应该扣除0.5 分2-3.解：由u x u ∂∂=()xuu ∂∂21=η22y u ∂∂得： 其守恒形式为：()xuu ∂∂=2η22y u ∂∂ 对方程两端在t ∆时间间隔内对其控制容积积分得：()dxdydt x uu t t t nsew ⎰⎰⎰∆+∂∂=⎰⎰⎰∆+∂∂t t t e w n s dydxdt y u 222η()()[]dxdt y u y u dydt uu uu s n t t t ewtt t w e n s ][2⎪⎪⎭⎫ ⎝⎛∂∂−⎪⎪⎭⎫ ⎝⎛∂∂=−⎰⎰⎰⎰∆+∆+η 将()()2)(PE e uu uu uu +=, ()()()2P W w uu uu uu +=,()n PN n y u u y u δ−=⎪⎪⎭⎫ ⎝⎛∂∂，()sSP s y u u y u δ−=⎪⎪⎭⎫ ⎝⎛∂∂。
y y y s n ∆==)()(δδ 带入，得：xdt y u u u ydt uu uu t t t S P N tt tW E ∆∆+−=∆⎥⎦⎤⎢⎣⎡−⎰⎰∆+∆+]2[22)()(η t x yu u u t y uu uu tSt P t N t W t E ∆∆∆+−=∆∆−222)()(η整理得离散方程为：()()0242=∆−+−∆−yu u u xuu uu t P t S t N tWt E η2—3：解：由2221()u 2u u ux x y η∂∂∂===∂∂∂得：原方程的守恒形式为： 222()2u ux yη∂∂=∂∂ 对方程两端在t ∆时间间隔内对其控制容积积分，把可积的部分积出后得：22()t tsne wtu u dtdy +∆−⎰⎰= 2t te wtn s u u dtdx y y η+∆⎡⎤⎛⎫⎛⎫∂∂−⎢⎥ ⎪ ⎪∂∂⎝⎭⎝⎭⎣⎦⎰⎰选定2u 随y 而变化的型线，这里取为阶梯式，即在控制容积内沿y 方向不变，则2222()=y ()t tt ts ne we w ttu u dtdy u u dt +∆+∆−∆−⎰⎰⎰选定2u 随t 而变化的规律，这里采用阶梯式显式，则22()t tewty u u dt +∆∆−⎰= ()()22t t e w u u t y ⎡⎤−∆∆⎢⎥⎣⎦选定uy∂∂随x 而变化的型线，这里取为阶梯式，即在控制容积内沿x 方向不变，则22t tt t e wtt n s n s u u u u dtdx x dt y y y y ηη+∆+∆⎡⎤⎡⎤⎛⎫⎛⎫⎛⎫⎛⎫∂∂∂∂−=∆−⎢⎥⎢⎥ ⎪ ⎪ ⎪ ⎪∂∂∂∂⎝⎭⎝⎭⎝⎭⎝⎭⎣⎦⎣⎦⎰⎰⎰ 选定uy∂∂随t 而变化的规律，这里采用阶梯显式，则 2t ttn s u u x dt y y η+∆⎡⎤⎛⎫⎛⎫∂∂∆−⎢⎥ ⎪ ⎪∂∂⎝⎭⎝⎭⎣⎦⎰= 2t t n s u u t x y y η⎡⎤⎛⎫⎛⎫∂∂−∆∆⎢⎥ ⎪ ⎪∂∂⎢⎥⎝⎭⎝⎭⎣⎦进一步选取u 随x,y 分段线性变化，则2222E Pe u u u += ， 222w 2W P u u u +=()nt PtN ty uu y u δ−=⎪⎪⎭⎫ ⎝⎛∂∂n ， ()stSt p ts y u u y u δ−=⎪⎪⎭⎫ ⎝⎛∂∂。
数值传热学习题答案
数值传热学习题答案【篇一：习题解答】s=txt>1.1 什么是汇编语言？汇编语言的特点是什么？答：为了克服机器语言难以记忆、表达和阅读的缺点，人们采用具有一定含义的符号作为助忆符，用指令助忆符、符号地址等组成的符号指令称为汇编格式指令(或汇编指令)。
汇编语言是汇编指令集、伪指令集和使用它们规则的统称。
① 55h，70h，70h，65h，72h② 53h，6ch，6fh，77h③ 43h，6fh，6dh，70h，75h，74h，65h，72h ④ 57h，68h，61h，74h 1.7 求下列带符号十进制数的8位基2码补码。
① +127② ?2③ ?128④ +2 答：汇编语言的特点是：（1）执行速度快。
（2）程序短小。
（3）可以直接控制硬件。
（4）可以方便地编译。
（5）辅助计算机工作者掌握计算机体系结构。
（6）程序编制耗时，可读性差。
（7）程序可移植性差。
1.2 把下列十进制数转换成二进制数、八进制数、十六进制数。
①127② 1021 ③ 0.875④ 6.25 答：① 1111111b；177q；7fh ② 1111111101；1775q；3fdh ③0.111 b；0.7q；0.eh ④ 110.01b；6.2q；6.4h1.3 把下列二进制数转换成十进制数。
① 1001.11 ② 101011.10011 ③ 111.011④1011.1答：① 9.75d ② 43.59375d③ 7.375d④ 11.5d1.4 把下列八进制数转换成十进制数。
① 573.06 ② 75.23 ③ 431.7④ 123.45 答：① 379.09375d ② 61.296875d③ 281.875④ 83.5781251.5 把下列十六进制数转换成十进制数。
① 0d5.f4 ② 8ba.7c ③ 0b2e.3a④ 6ec.2d 答：① 213.953125d ② 2234.484375③ 2862.2265625 ④1772.175781251.6 把下列英文单词转换成ascii编码的字符串。
传热学课后答案(完整版)
绪论思考题与习题（89P -）答案：1．冰雹落体后溶化所需热量主要是由以下途径得到： Q λ—— 与地面的导热量 f Q ——与空气的对流换热热量注：若直接暴露于阳光下可考虑辐射换热，否则可忽略不计。
2．略 3．略 4．略 5．略6．夏季：在维持20℃的室内，人体通过与空气的对流换热失去热量，但同时又与外界和内墙面通过辐射换热得到热量，最终的总失热量减少。
（T T 〉外内）冬季：在与夏季相似的条件下，一方面人体通过对流换热失去部分热量，另一方面又与外界和内墙通过辐射换热失去部分热量，最终的总失热量增加。
（T T 〈外内）挂上窗帘布阻断了与外界的辐射换热，减少了人体的失热量。
7．热对流不等于对流换热，对流换热 = 热对流 + 热传导 热对流为基本传热方式，对流换热为非基本传热方式 8．门窗、墙壁、楼板等等。
以热传导和热对流的方式。
9．因内、外两间为真空，故其间无导热和对流传热，热量仅能通过胆壁传到外界，但夹层两侧均镀锌，其间的系统辐射系数降低，故能较长时间地保持热水的温度。
当真空被破坏掉后，1、2两侧将存在对流换热，使其保温性能变得很差。
10．t R R A λλ= ⇒ 1t R R A λλ== 2218.331012m --=⨯11．q t λσ=∆ const λ=→直线 const λ≠ 而为λλ=（t ）时→曲线12、略13．解：1211t q h h σλ∆=++＝18(10)45.9210.361870.61124--=++2W m111()f w q h t t =-⇒ 11137.541817.5787w f q t t h =-=-=℃222()w f q h t t =-⇒ 22237.54109.7124w f q t t h =+=-+=-℃ 45.92 2.83385.73q A W φ=⨯=⨯⨯= 14. 解：40.27.407104532t K R W A HL λσσλλ-====⨯⨯⨯30.2 4.4441045t R λσλ-===⨯2m K W • 3232851501030.44.44410t KW q m R λ--∆-==⨯=⨯ 3428515010182.37.40710t t KW R λφ--∆-==⨯=⨯ 15．()i w f q h t h t t =∆=-⇒i w f qt t h=+51108515573=+=℃0.05 2.551102006.7i Aq d lq W φππ===⨯⨯=16.解：12441.2 1.2()()100100w w t t q c ⎡⎤=-⎢⎥⎣⎦ 44227350273203.96()()139.2100100W m ++⎡⎤=⨯-=⎢⎥⎣⎦12''441.21.2()()100100w w t t qc ⎡⎤=-⎢⎥⎢⎥⎣⎦442273200273203.96()()1690.3100100W m ++⎡⎤=⨯-=⎢⎥⎣⎦'21.2 1.2 1.21690.3139.21551.1Wq q q m ∆=-=-=17．已知：224A m =、215000()Wh m K =•、2285()Wh m K =•、145t =℃2500t =℃、'2285()Wk h m K ==•、1mm σ=、398λ=()W m K •求：k 、φ、∆解：由于管壁相对直径而言较小，故可将此圆管壁近似为平壁即：12111k h h σλ=++＝3183.5611101500039085-=⨯++2()W m k • 383.5624(50045)10912.5kA t KW φ-=∆=⨯⨯-⨯= 若k ≈2h'100k k k -∆=⨯％8583.561.7283.56-==％ 因为：1211h h ，21h σλ 即：水侧对流换热热阻及管壁导热热阻远小于燃气侧对流换热热阻，此时前两个热阻均可以忽略不记。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
习题4－2一维稳态导热问题的控制方程：022=+∂∂S xTλ 依据本题给定条件，对节点2节点3采用第三类边界条件具有二阶精度的差分格式，最后得到各节点的离散方程： 节点1： 1001=T节点2： 1505105321-=+-T T T 节点3：75432=+-T T求解结果：852=T ，403=T对整个控制容积作能量平衡，有：02150)4020(15)(3=⨯--⨯=∆+-=∆+x S T T h x S q f f B即：计算区域总体守恒要求满足习题4－5在4－2习题中，如果25.03)(10f T T h -⨯=，则各节点离散方程如下：节点1： 1001=T节点2： 1505105321-=+-T T T节点3：25.03325.032)20(4015])20(21[-⨯+=-⨯++-T T T T对于节点3中的相关项作局部线性化处理，然后迭代计算； 求解结果：818.822=T ，635.353=T （迭代精度为10-4）迭代计算的Matlab 程序如下： x=30; x1=20;while abs(x1-x)>0.0001a=[1 0 0;5 -10 5;0 -1 1+2*(x-20)^(0.25)]; b=[100;-150; 15+40*(x-20)^(0.25)]; t=a^(-1)*b; x1=x; x=t(3,1);endtcal=t习题4-12的Matlab程序%代数方程形式A i T i=C i T i+1+B i T i-1+D imdim=10;%计算的节点数x=linspace(1,3,mdim);%生成A、C、B、T数据的基数；A=cos(x);%TDMA的主对角元素B=sin(x);%TDMA的下对角线元素C=cos(x)+exp(x); %TDMA的上对角线元素T=exp(x).*cos(x); %温度数据%由A、B、C构成TDMAcoematrix=eye(mdim,mdim);for n=1:mdimcoematrix(n,n)=A(1,n);if n>=2coematrix(n,n-1)=-1*B(1,n);endif n<mdimcoematrix(n,n+1)=-1*C(1,n);endend%计算D矢量D=(coematrix*T')';%由已知的A、B、C、D用TDMA方法求解T%消元P(1,1)=C(1,1)/A(1,1);Q(1,1)=D(1,1)/A(1,1);for n=2:mdimP(1,n)=C(1,n)/(A(1,n)-B(1,n)*P(1,n-1));Q(1,n)=(D(1,n)+B(1,n)*Q(1,n-1))/(A(1,n)-B(1,n)*P(1,n-1));end%回迭Tcal(1,mdim)=Q(1,mdim);for n=(mdim-1):-1:1Tcal(1,n)=P(1,n)*Tcal(1,n+1)+Q(1,n);endTcom=[T;Tcal];%绘图比较给定T值和计算T值plot(Tcal,'r*')hold onplot(T)结果比较如下，由比较可知两者值非常切合（在小数点后8位之后才有区别）：习题4-14充分发展区的温度控制方程如下：)(1rTr r r x T uc p ∂∂∂∂=∂∂λρ 对于三种无量纲定义w b w T T T T --=Θ、∞∞--=ΘT T T T w 、w w T T T T --=Θ∞进行分析如下1）由wb wT T T T --=Θ得：w w b T T T T +Θ-=)(由T 可得：x T x T x T T T x T w b w w b ∂∂Θ-+∂∂Θ=∂+Θ-∂=∂∂)1(])[(rT r T T r T T T r T w w b w w b ∂∂Θ-+∂Θ∂-=∂+Θ-∂=∂∂)1()(])[( 由b T 与r 无关、Θ与x 无关以及x T ∂∂、rT∂∂的表达式可知，除了w T 均匀的情况外，该无量纲温度定义在一般情况下是不能用分离变量法的； 2）由∞∞--=ΘT T T T w 得： ∞∞+Θ-=T T T T w )(由T 可得：xT x T T T x T w w ∂∂Θ=∂+Θ-∂=∂∂∞∞])[(rT r T T r T T T r T w w w ∂∂Θ+∂Θ∂-=∂+Θ-∂=∂∂∞∞∞)(])[( 由b T 与r 无关、Θ与x 无关以及x T ∂∂、rT∂∂的表达式可知，在常见的四种边界条件中除了轴向及周向均匀热流const q w =的情况外，有0=∂∂rT w，则该无量纲温度定义是可以用分离变量法的； 3）由wwT T T T --=Θ∞得： w w T T T T +Θ-=∞)(由T 可得：xT x T T T x T w w w ∂∂Θ-=∂+Θ-∂=∂∂∞)1(])[(rT r T T r T T T r T w w w w ∂∂Θ-+∂Θ∂-=∂+Θ-∂=∂∂∞∞)1()(])[( 同2）分析可知，除了轴向及周向均匀热流const q w =的情况外，有0=∂∂rT w，该无量纲温度定义是可以用分离变量法的；习题4－181）采用柱坐标分析，写出统一的稳态柱坐标形式动量方程：r r x x w r v r r r u x ∂∂+∂∂∂∂=∂∂+∂∂+∂∂1)()(1)(1)(φλφρθφρφρx 、r 和θ分别是圆柱坐标的3个坐标轴，u 、v 和w 管内的流动方向；对于管内的层流充分发展有：0=v 、0=w ，0=∂∂xu； 并且x 方向的源项：x pS ∂∂-=r 方向的源项：r pS ∂∂-= θ方向的源项：θ∂∂-=pr S 1 由以上分析可得到圆柱坐标下的动量方程： x 方向： 0)(1)(1=∂∂-∂∂∂∂+∂∂∂∂x pu r r r u r r r θλθλ r 方向：0=∂∂r pθ方向：0=∂∂θp边界条件： R r =，0=u0=r ，0=∂∂r u ；对称线上，0=∂∂θu 不考虑液体的轴向导热，并简化分析可以得到充分发展的能量方程为：)(1)(1θλθλρ∂∂∂∂+∂∂∂∂=∂∂Tr r r T r r r x T uc p 边界条件： R r =，w q r T =∂∂λ；0=r ，0=∂∂rTπθ/0=，0=∂∂-θλT2）定义无量纲流速：dxdp R uU 2-=λ并定义无量纲半径：R r /=η；将无量纲流速和无量纲半径代入x 方向的动量方程得：0))1((1))1((122=∂∂-∂-∂∂∂+∂-∂∂∂xp U dx dp R R R R U dx dp R RR R θληλθηηλληηη 上式化简得：01)1(1)(1=+∂∂∂∂+∂∂∂∂θηθηηηηηU U 边界条件：1=η，0=U0=η，0=∂∂ηU ；对称线上，0=∂∂θU定义无量纲温度：λ/0R q T T b-=Θ其中，0q 是折算到管壁表面上的平均热流密度，即：Rq q wπ=0； 由无量纲温度定义可得：b T Rq T +Θ=λ将T 表达式和无量纲半径η代入能量方程得：)(1)(100θληλθηηλληηηρ∂Θ∂∂∂+∂Θ∂∂∂=∂∂R q R R R R q R R R x T uc b p 化简得：)1(1)(10θηθηηηηηρ∂Θ∂∂∂+∂Θ∂∂∂=∂∂x T u c q R b p （1）由热平衡条件关系可以得：mm m b m p b p p RU U q R u u R q A u u dx dT A u c x T u c x T uc 020221221)(===∂∂=∂∂ππρρρ 将上式代入式（1）可得：)1(1)(12θηθηηηηη∂Θ∂∂∂+∂Θ∂∂∂=m U U 边界条件：0=η，0=∂Θ∂η；1=η，R q q w πη10==∂Θ∂0=θ，0=∂Θ∂θ；πθ=，0=∂Θ∂θ单值条件： 由定义可知：0/0=-=ΘλR q T T b b b 且： ⎰⎰Θ=ΘAAb UdAUdA即得单值性条件：0=Θ⎰⎰AA UdAUdA 3）由阻力系数f 及Re 定义有：228)(21/Re ⎪⎭⎫ ⎝⎛=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-=D D U D u u dx dp D f e m e m me νρ 且：m W b m W b m W R q T T D T T q Nu ,0,,0~2)/(2Θ=-=-=λλ5－21．一维稳态无源项的对流－扩散方程如下所示：xx u 22∂∂Γ=∂∂φφρ （取常物性）边界条件如下：L L x x φφφφ====,;,00上述方程的精确解如下：11)/(00--=--⋅PeL x Pe L e e φφφφ Γ=/uL Pe ρ 2．将L 分成20等份，所以有：∆=P Pe 201 2 3 4 5 6 ………… …………… 17 18 19 20 21 对于中心差分、一阶迎风、混合格式和QUICK 格式分别分析如下： 1） 中心差分中间节点： 2)5.01()5.01(11-∆+∆++-=i i i P P φφφ 20,2 =i2） 一阶迎风中间节点： ∆-∆++++=P P i i i 2)1(11φφφ 20,2 =i3） 混合格式当1=∆P 时，中间节点：2)5.01()5.01(11-∆+∆++-=i i i P P φφφ20,2 =i当10,5=∆P 时，中间节点： 1-=i i φφ 20,2 =i 4） QUICK 格式*12111)35(8122121⎥⎦⎤⎢⎣⎡---++++++=+--∆∆-∆∆+∆i i i i i i i P P P P P φφφφφφφ 2≠i *1111)336(8122121⎥⎦⎤⎢⎣⎡--++++++=+-∆∆-∆∆+∆i i i i i i P P P P P φφφφφφ 2=i数值计算结果与精确解的计算程序如下：%except for HS, any other scheme doesnt take Pe<0 into consideration %expression of exact solutiony=dsolve('a*b*Dy=c*D2y','y(0)=y0,y(L)=yL','x')y=subs(y,'L*a*b/c','t')y=simple(subs(y,'a*b/c*x','t*X'));ysim=simple(sym(strcat('(',char(y),'-y0)','/(yL-y0)')))y=sym(strcat('(',char(ysim),')*(yL-y0)','+y0'))% in the case of Pe=0y1=dsolve('D2y=0','y(0)=y0,y(L)=yL','x')y1=subs(y1,'-(y0-yL)/L*x','(-y0+yL)*X')%grid Pe numbertt=[1 5 10];%dimensionless lengthm=20;%mdim is the number of inner nodemdim=m-1;X=linspace(0,1,m+1);%initial value of variable during calculationy0=1;yL=2;%cal exact solutionfor n=1:size(tt,2)t=m*tt(1,n);if t==0yval1(n,:)=eval(y1);elseyval1(n,:)=eval(y);endend%extra treatment because max number in MATLAB is 10^308if max(isnan(yval1(:)))yval1=yval1';yval1=yval1(:);indexf=find(isnan(yval1));for n=1:size(indexf,1)if rem(indexf(n,1),size(X,2))==0yval1(indexf(n),1)=yL;elseyval1(indexf(n),1)=y0;endendyval1=reshape(yval1,size(X,2),size(yval1,1)/size(X,2));yval1=yval1';end%CD solutiond=zeros(size(tt,2),mdim);a=repmat([1],size(tt,2),mdim);for n=1:size(tt,2)t=tt(1,n);b(n,:)=repmat([0.5*(1-0.5*t)],1,mdim);c(n,:)=repmat([0.5*(1+0.5*t)],1,mdim);d(n,1)=0.5*(1+0.5*tt(1,n))*y0;d(n,mdim)=0.5*(1-0.5*tt(1,n))*yL;endc(:,1)=0;b(:,mdim)=0;%numerical cal by using TDMA subfuctionyval2=TDMA(a,b,c,d,mdim);yval2=[repmat([1],size(tt,2),1),yval2,repmat([2],size(tt,2),1)]; Fig(1,X,yval1,yval2,tt);title('CD Vs. Exact Solution')% FUS solutiond=zeros(size(tt,2),mdim);a=repmat([1],size(tt,2),mdim);for n=1:size(tt,2)t=tt(1,n);b(n,:)=repmat([1/(2+t)],1,mdim);c(n,:)=repmat([(1+t)/(2+t)],1,mdim);d(n,1)=(1+tt(1,n))/(2+tt(1,n))*y0;d(n,mdim)=1/(2+tt(1,n))*yL;endc(:,1)=0;b(:,mdim)=0;%numerical cal by using TDMA subfuctionyval3=TDMA(a,b,c,d,mdim);yval3=[repmat([1],size(tt,2),1),yval3,repmat([2],size(tt,2),1)]; Fig(2,X,yval1,yval3,tt);title('FUS Vs. Exact Solution')% HS solutiond=zeros(size(tt,2),mdim);a=repmat([1],size(tt,2),mdim);for n=1:size(tt,2)t=tt(1,n);if t>2b(n,:)=repmat([0],1,mdim);c(n,:)=repmat([1],1,mdim);d(n,1)=y0;elseif t<-2b(n,:)=repmat([1],1,mdim);c(n,:)=repmat([0],1,mdim);d(n,mdim)=yL;elseb(n,:)=repmat([0.5*(1-0.5*t)],1,mdim);c(n,:)=repmat([0.5*(1+0.5*t)],1,mdim);d(n,1)=0.5*(1+0.5*t)*y0;d(n,mdim)=0.5*(1-0.5*t)*yL;endendc(:,1)=0;b(:,mdim)=0;% numerical cal by using TDMA subfuctionyval4=TDMA(a,b,c,d,mdim);yval4=[repmat([1],size(tt,2),1),yval4,repmat([2],size(tt,2),1)]; Fig(3,X,yval1,yval4,tt);title('HS Vs. Exact Solution')%QUICK Solutiond=zeros(size(tt,2),mdim);a=repmat([1],size(tt,2),mdim);for n=1:size(tt,2)t=tt(1,n);b(n,:)=repmat([1/(2+t)],1,mdim);c(n,:)=repmat([(1+t)/(2+t)],1,mdim);d(n,1)=(1+tt(1,n))/(2+tt(1,n))*y0;d(n,mdim)=1/(2+tt(1,n))*yL;endc(:,1)=0;b(:,mdim)=0;%numerical cal by using TDMA subfuctionyval5=zeros(size(tt,2),mdim);yval5com=yval5+1;counter=1;%iterativewhile max(max(abs(yval5-yval5com)))>10^-10if counter==1yval5com=TDMA(a,b,c,d,mdim);endfor nn=1:size(tt,2)for nnn=1:mdimif nnn==1d(nn,nnn)=((6*yval5com(nn,nnn)-3*y0-3*yval5com(nn,nnn+1))*tt(1,nn))/(8*(2+tt(1, nn)))+((1+tt(1,nn))/(2+tt(1,nn))*y0);elseif nnn==2d(nn,nnn)=((5*yval5com(nn,nnn)-3*yval5com(nn,nnn+1)-yval5com(nn,nnn-1)-y0)*tt (1,nn))/(8*(2+tt(1,nn)));elseif nnn==mdimd(nn,nnn)=((5*yval5com(nn,nnn)-3*yL-yval5com(nn,nnn-1)-yval5com(nn,nnn-2))*tt (1,nn))/(8*(2+tt(1,nn)))+(1/(2+tt(1,nn))*yL);elsed(nn,nnn)=((5*yval5com(nn,nnn)-3*yval5com(nn,nnn+1)-yval5com(nn,nnn-1)-yval5 com(nn,nnn-2))*tt(1,nn))/(8*(2+tt(1,nn)));endendendyval5=TDMA(a,b,c,d,mdim);temp=yval5;yval5=yval5com;yval5com=temp;counter=counter+1;endyval5=yval5com;yval5=[repmat([1],size(tt,2),1),yval5,repmat([2],size(tt,2),1)];Fig(4,X,yval1,yval5,tt);title('QUICK Vs. Exact Solution')%-------------TDMA SubFunction------------------function y=TDMA(a,b,c,d,mdim)%form a b c d resolve yval2 by using TDMA%eliminationp(:,1)=b(:,1)./a(:,1);q(:,1)=d(:,1)./a(:,1);for n=2:mdimp(:,n)=b(:,n)./(a(:,n)-c(:,n).*p(:,n-1));q(:,n)=(d(:,n)+c(:,n).*q(:,n-1))./(a(:,n)-c(:,n).*p(:,n-1));end%iterativey(:,mdim)=q(:,mdim);for n=(mdim-1):-1:1y(:,n)=p(:,n).*y(:,n+1)+q(:,n);end%-------------ResultCom SubFunction------------------ function y=ResultCom (a,b,c)for n=1:max(size(c,2))y(2*n-1,:)=a(n,:);y(2*n,:)=b(n,:);end%-------------Fig SubFunction------------------ function y=Fig(n,a,b,c,d)figure(n);plot(a,b);hold onplot(a,c,'*');str='''legend(';for n=1:size(d,2)if n==size(d,2)str=strcat(str,'''''Pe=',num2str(d(1,n)),''''')''');elsestr=strcat(str,'''''Pe=',num2str(d(1,n)),''''',');endendeval(eval(str));精确解与数值解的对比图，其中边界条件给定10=φ，2=L φ。