数值传热学陶文铨第四章作业
计算传热学_高等教育-实验设计
坐标系不同，控制方程的形式不尽相同
必要的简化与化简
2.1 控制方程
传热的三种模式（Modes of heat transfer)
热传导（Thermal conduction) 热对流（Thermal advection)
对流换热（Convection heat transfer) 辐射换热（Radiation heat transfer)
数值方法
分析解法与实验研究
分析解法
成本最低 结果最理想 影响因素表达清楚 缺点：局限与非常简单的问题 成本较低：数值实验 适用范围宽 缺点：可靠性差，表达困难 可靠 成本高
数值方法
实验研究
将三种方法有 机结合，互为 补充，必然会 取得相得益彰 的效果
第2讲
传热问题的数学描述
1) 2)
将上面的数学模型无量纲化，并给出其分析解； 取β=1， 就 PeL=(ρuL)/Γ=1、10、100 三种情况分别用三点中心差分格式、迎风格式、幂律格式和 QUICK 格式进行计算，并与分析解比较（计算时节点数目可取为 10 ~ 20） ； 3) 改变参数β，譬如取β=10，重复 2）中的计算； 分析 2）和 3）中得到的结果，对各种格式进行比较。
Tf h A B
Tf h
δ
δ
计算传热学习题之三
考虑下述一维稳态对流-扩散问题,
d d dU ( ρuU ) = (Γ )+s dx dx dx U x=0 = U 0 U
x=L
= UL
其中 u 是流速，Γ和ρ均为常数，而 s 是 x 的单值函数，
s = 0.5 β U 0 −U L L
数值传热学(陶文铨)第二章习题答案
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．什么是非稳态导热问题的显示格式？什么是显示格式计算中的稳定性问题？
7．用高斯－塞德尔迭代法求解代数方程时是否一定可以得到收敛德解？不能得出收敛的解时是否因为初场的假设不合适而造成？
8．有人对一阶导数
你能否判断这一表达式是否正确，为什么？
Fo=0.24
Bi=0.1
Bi=1
Bi=10
第一项的值
0.99277
0.93698
0.77311
前六项和的值
0.99101
0.92791
0.76851
比值
1.00177
1.00978
1.00598
4-2、试用数值计算证实，对方程组
用高斯-赛德尔迭代法求解，其结果是发散的，并分析其原因。
解：将上式写成下列迭代形式
其中第五次与第六次相对偏差已小于 迭代终止。
4-4、试对附图所示的等截面直肋的稳态导热问题用数值方法求解节点2，3的温度。图中 .肋高H=4cm,纵剖面面积 导热系数 。
解：对于2点可以列出：
节点2：
节点3： 。
由此得：
， ，
，于是有： ，
，代入得：
， ，
， ，
，
。
离散方程的建立
4-5、试将直角坐标中的常物性无内热源的二维稳态导热微分方程化为显式差分格式，并指出其稳定性条件（ 。
3△
200
116.98
42.63
42.23
4△
200
125.51
52.57
51.94
4-16、一厚为2.54cm的钢板，初始温度为650℃，后置于水中淬火，其表面温度突然下降为93.5℃并保持不变。试用数值方法计算中心温度下降到450℃所需的时间。已知 。建议将平板8等分，取9个节点，并把数值计算的结果与按海斯勒计算的结果作比较。
传热学第4章热传导问题的数值解法重点习题
t1 t5 y t9 t5 x t 6 t5 1 y xy yh t5 t f 0 y 2 x 2 节点 5： y 2 ； t 2 t6 t7 t6 t10 t5 t5 t 6 x y x y xy 0 y x y x 节点、一等截面直肋，高 H,厚 ，肋根温度为 t 0 ，流体温度为 t f ，表 面传热系数为 h,肋片导热系数为 。将它均分成 4 个节点（见附图） ， 并对肋端为绝热及为对流边界条件（h 同侧面）的两种情况列出节点 2 ， 3 ， 4 的 离 散 方 程 式 。 设
节点 2： 节点 3：
t3 t 2
x
2hx t2 t f 0 2hx t3 t f 0
t 2 t3
x
t 4 t3
x x
； ；
t3 t 4
节点 4：肋端绝热 肋端对流
0 0 0 由此解得：肋端绝热 t2 92.2 C ， t3 87.7 C ， t4 86.2 C ；
肋端对流 t2 91.5 C ， t3 86.2 C ， t4 83.8 C 。 肋端对流换热的条件使肋端温度更接近于流体温度。
0 0 0
传热学第4章热传导问题的数值解法重点习题数值传热学传热学课后习题答案数值传热学答案数值传热学第二版答案数值传热学陶文铨数值传热学第二版pdf传热学习题解答数值传热学pdf传热学课后习题
第 4 章热传导问题的数值解法
一般性数值计算
4-4、试对附图所示的等截面直肋的稳态导热问题用数值方法求解节 点 2， 3 的温度。 图中
2 H=45cm, 10mm, h 50W /(m .K ) , =50W/(m.K), t 0 100 ℃， t f 20 ℃， 计算节点 2，3，4 的温度（对于肋端的两种边界条件） 。
第一章数值传热学
(uT ) (vT ) T T a( 2 2 ) x y x y
2 2
19/80
MOE KLTFSE
3. 边界条件
定u,v,T随 y 的分布；
（1）进口边界条件：给
u T （3）中心线： 0; v 0 y y
y x
界：数学上要 求给定u,v,T或 其导数随 y 的 分布；实际上 做不到；数值 上近似处理。
cp
c p
( ) c p
Pr
12/80
MOE KLTFSE
4. 通用控制方程
( ) * * div( U ) div( grad ) S t
瞬态项 对流项 扩散项 广义源项 不同求解变量之间的区别： （1）边界条件与初始条件不同； （2）广义源项表达式不同； （3）广义扩散系数不同。 文献中常以表格形式给出所求解变量的源项与 广义扩散系数的表达式。
常物性不可压缩流体动量方程源项中显含速度部分 为零。
11/80
MOE KLTFSE
3. 能量守恒方程
[微元体内热力学能的增加率]=[进入微元体内的净热 流量]+[体积力与表面力对微元体所做的功] 引入导热Fourier定律，忽略力所作的功， 设hc
pT ;
c p 为常数
( T ) div( T U ) div( gradT ) ST cp t
4/80
MOE KLTFSE
绪论教学目录
1.1 传热与流动问题的数学描写 1.2 传热与流动问题数值计算的基本思想及应 用举例 1.3 传热与流动问题的数学描写的分类及其对 数值解的影响 1.4 传热与流动问题的数值计算的近代发展
5/80
传热学第四版课后思考题答案(杨世铭-陶文铨)]之欧阳物创编
![传热学第四版课后思考题答案(杨世铭-陶文铨)]之欧阳物创编](https://img.taocdn.com/s3/m/cb874f98168884868662d685.png)
第一章思考题1． 试用简练的语言说明导热、对流换热及辐射换热三种热传递方式之间的联系和区别。
答：导热和对流的区别在于：物体内部依靠微观粒子的热运动而产生的热量传递现象，称为导热；对流则是流体各部分之间发生宏观相对位移及冷热流体的相互掺混。
联系是：在发生对流换热的同时必然伴生有导热。
导热、对流这两种热量传递方式，只有在物质存在的条件下才能实现，而辐射可以在真空中传播，辐射换热时不仅有能量的转移还伴有能量形式的转换。
2． 以热流密度表示的傅立叶定律、牛顿冷却公式及斯忒藩－玻耳兹曼定律是应当熟记的传热学公式。
试写出这三个公式并说明其中每一个符号及其意义。
答：①傅立叶定律：dx dt q λ-=，其中，q －热流密度；λ－导热系数；dx dt －沿x 方向的温度变化率，“－”表示热量传递的方向是沿着温度降低的方向。
②牛顿冷却公式：)(f w t t h q -=，其中，q －热流密度；h －表面传热系数；w t －固体表面温度；f t －流体的温度。
③斯忒藩－玻耳兹曼定律：4T=，其中，q－热流密qσ度；σ－斯忒藩－玻耳兹曼常数；T－辐射物体的热力学温度。
3．导热系数、表面传热系数及传热系数的单位各是什么？哪些是物性参数，哪些与过程有关？答：①导热系数的单位是：W/(m.K)；②表面传热系数的单位是：W/(m2.K)；③传热系数的单位是：W/(m2.K)。
这三个参数中，只有导热系数是物性参数，其它均与过程有关。
4．当热量从壁面一侧的流体穿过壁面传给另一侧的流体时，冷、热流体之间的换热量可以通过其中任何一个环节来计算（过程是稳态的），但本章中又引入了传热方程式，并说它是“换热器热工计算的基本公式”。
试分析引入传热方程式的工程实用意义。
答：因为在许多工业换热设备中，进行热量交换的冷、热流体也常处于固体壁面的两侧，是工程技术中经常遇到的一种典型热量传递过程。
5．用铝制的水壶烧开水时，尽管炉火很旺，但水壶仍然安然无恙。
Chapter_4(2)_
主讲陶文铨西安交通大学能源与动力工程学院热流中心CFD-NHT-EHT CENTER 2010年9月27日, 西安数值传热学第四章扩散方程的数值解及其应用(2)4.4 求解代数方程的TDMA 及ADI 方法4.4.1求解一维导热问题代数方程的三对角阵算法4.4.2求解多维非稳态导热全隐格式的ADI 方法1.求解方法概述1.一维导热问题代数方程通用形式2.Peaceman-Rachford 的ADI 迭代2.Thomas 算法3.第一类边界条件的处理4.4 方法4.4.1 求解一维导热问题代数方程的三对角阵算法Thomas算法的一般形式将上式改写为：所谓消元就是要找出系数间的关系。
(b)乘以CiP=11;i−−1,i=B D3.第一类边界条件下Thomas 算法的实施第一类边界条件下，求解区域为i=2,….M1-1=M2。
将消元公式用于i=1, 注意T 1是给定的：1121T PT Q =+10;P =11Q T =因T M1已知，消元从M 2开始：2212M M M T P T Q =+注意：采用附加源项法来处理第二类，第三类边界条件时，均将第二类，第三类边界条件问题视为第一类边界条件问题，数学上的处理与此相同。
4.4.2 求解多维导热问题代数方程的方法求解二维非稳态导热全隐格式代数方程的方法(1) 五对角阵算法(Penta-diagonal ,PDMA)(2) 交替方向隐式方法(Alternative-directionImplicit, ADI)2. 3-D Peaceman-Rachford方法将tΔ三等分：tΔX方向为隐式，第一个/3y,z方向为显式式；第二，三个/3tΔ分别在y,z方向实施隐式;2－D交替方向隐式u i,j,k , v 表示方向二阶导数的中心差分；nT用von Neumann分析方法可以证明稳定性条件为：4.5 管道内充分发展对流换热概说4.5.1管道内充分发展对流换热的定义1. 简单的充分发展对流换热2. 复杂的充分发展对流换热4.5.2能实现充分发展对流换热的边界条件4.5.3部分算例汇总4.5 管道内充分发展对流换热概说4.5.1管道内充分发展对流换热的定义平直通道中的充分发展对流换热属于这一类。
数值传热学第4章作业
习题4－2图4-22 习题4-2插图[解]一维稳态导热问题的控制方程为：0=+⎪⎭⎫ ⎝⎛S dx dT dx d λ 4-2-1 该问题的边界条件为：()⎪⎩⎪⎨⎧=-=-==2,0,1001x T T h dx dT x T f λ 4-2-2分别对节点2，3进行离散，将已知数据代入离散格式中，得到方程组：130232=-T T 4-2-375432=+-T T 4-2-4 联立式（4-2-3）、式（4-2-4），可以解出2T ，3T ： 852=T ，403=T 。
下面验证总体守恒性：4-2-5右端3放出的热量为：()()30020401533=-⨯=-=f T T h Q 4-2-6在总体容积内部产生的热量为：2.0150 2.0300S Q S x =⨯∆=⨯=还需要证明左端是绝热条件： 节点2的热平衡为：21851000.550.5150757501T T xS x λ--+∆=+⨯=-+=∆ 左端绝热，所以计算结果符合总体能量守恒。
习题 4-5[解] 根据习题4-2的分析，可以得到节点2的离散方程：130232+=T T 4-5-1对于节点3，应用边界条件：()()1324330.510f f T T S T T T T xλδ--+⨯=-- 4-5-2式（4-5-2）可以整理成：()5432355751020T T T =+-- 4-5-3采用局部线性化方法，可以得到：()()()()515***444333331020102012.520T T T T T -=-+-- 4-5-4节点3的离散方程表示成：()()()51**44323335575 2.52012.52020T T T T T =++---- 4-5-5迭代求解得出：2382.82;35.64T T == 检验热平衡：内热源生成热1300φ=； 右端散热5/4210(35.6420)311.0h T φ=∆=-=左端散热382.821000.5150510.91φ-=⨯+⨯=-所以123()30031110.90φφφ-+=-+≅不作热平衡扣0.5 分。
计算流体力学与传热学大作业
########学院计算流体力学与传热学学号：专业：学生姓名：任课教师：教授2013年12月目录第一章验证显式格式的稳定性 (4)1.1 概述 (4)1.2 数学推导 (4)1.3 问题描述 (4)1.4 数值模拟 (4)1.5 结果及分析 (5)第二章判断肋片可以按一维问题处理的主要依据 (6)2.1 概述 (6)2.2 问题描述及算法 (6)2.3 数值模拟 (7)2.4 结果及分析 (8)第三章三层墙导热 (9)3.1 概述 (9)3.2 问题描述 (9)3.3 TDMA算法 (9)3.4 结果 (10)第四章一维无源稳态对流扩散问题 (11)4.1 公式及初值 (11)4.2 情况一 (11)4.3 情况二 (12)4.4 情况三 (13)第五章用ADI算法计算长方肋内的温度分布 (14)5.1 问题描述 (14)5.2 初始参数 (14)5.3 情况一，一列列扫 (14)5.4 情况二，一行行扫 (14)5.5 情况三，采用ADI算法 (15)5.6 结果分析 (15)参考文献 (16)第一章 验证显式格式的稳定性1.1 概述将一维非稳态热传导方程用显式格式差分化为代数方程，在求解的迭代过程中必须满足一定的条件，才能使方程收敛且结果正确。
此处即验证β≤½。
1.2 数学推导方程： 22T t T x α∂∂=∂∂（1）显式离散格式： 此处时间向前差分，空间中心差分11122n n n n ni i i i i T T T T T t x α+-+--+=∆∆1112(2)n n n n ni i i i i t T T T T T xα+-+∆-=-+∆ 令β=2tx α∆∆则： 111(2)n n n n ni i i i i T T T T T β+-+-=-+ （2）误差也应该满足上式，故：()()1()()()2()()iiiiiIkx Ikx Ik x x Ikx Ik x x n n n n nT e T e T e T e T e ψψβψψψ----∆--+∆+⎡⎤-=-+⎣⎦()()()1()12()()()iiiiIkx Ikx Ik x x Ik x x n n n nT e T e T e T e ψβψβψψ----∆-+∆+⎡⎤=-++⎣⎦()()1()12()()iiiIkx Ikx Ikxn n Ik x Ik x n T e T e e e T e ψβψβψ---+-∆∆=-++()()1()121()n Ik x Ik x nT e e T ψββψ+-∆∆=-++≤ 因此 β≤½。
传热学思考题参考答案(陶文铨第四版)
答：放大镜的中间厚，边缘薄，光线在透过放大镜时会产生折射，因此会把物图像放大。要点： 值越大则温度变化率越小，在图上标示出来就是斜率越小（具体可参考换热器原理一书）。当相等时，顺流为对称的两曲线，而逆流时则为平行线。
答：在圆管外敷设保温层和设置肋片都使表面换热热阻降低而导热热阻增加，而一般情况下保温使导热热阻增加较多，使换热热阻降低较少，使总热阻增加，起到削弱传热的效果；设置肋片使导热热阻增加较少，而换热热阻降低较多，使总热阻下降，起到强化传热的作用。但当外径小于临界直径时，增加保温层厚度反而会强化传热。理论上只有当肋化系数与肋面总效率的乘积小于1时，肋化才会削弱传热。
答：条件：（1）材料的导热系数，表面传热系数以及沿肋高方向的横截面积均各自为常数（2）肋片温度在垂直纸面方向（即长度方向）不发生变化，因此可取一个截面（即单位长度）来分析（3）表面上的换热热阻远远大于肋片中的导热热阻，因而在任一截面上肋片温度可认为是均匀的（4）肋片顶端可视为绝热。并不是扩展表面细长就可以按一维问题处理，必须满足上述四个假设才可视为一维问题。
第八章：
1、选择太阳能集热器的表面涂层时，该涂层表面吸收率随波长的变化最佳曲线是什么？有人认为取暖用的辐射采暖片也需要涂上这种材料，你认为合适吗？
分析：太阳辐射的主要能量集中在0.2～2μm，该涂层表面吸收率随波长的变化最佳曲线是当波长小于2μm时，吸收率大，当波长大于2μm时，吸收率要小。
不合适。因为如果暖片在高温（波长小）时有很大的吸收比，那么暖片将有很大的辐射换热量，减小了对流换热量，因此不适合。
答：虽然黑体表面与重辐射面均具有J=Eb的特点，但二者具有不同的性质。黑体表面的温度不依赖于其他参与辐射的表面，相当于源热势。而重辐射面的温度则是浮动的，取决于参与辐射的其他表面。
数值传热学陶文铨主编第二版习题答案
数值传热学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对整个控制容积作能量平衡，有：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(])[(T r T T r T T T r T w 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)(θφλθφλφλφρθφρφρ x 、r 和θ分别是圆柱坐标的3个坐标轴，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=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精确解与数值解的对比图，其中边界条件给定，。
传热学第四版课后思考题答案(杨世铭-陶文铨)]之欧阳地创编
![传热学第四版课后思考题答案(杨世铭-陶文铨)]之欧阳地创编](https://img.taocdn.com/s3/m/a2de67350975f46526d3e121.png)
第一章思考题1． 试用简练的语言说明导热、对流换热及辐射换热三种热传递方式之间的联系和区别。
答：导热和对流的区别在于：物体内部依靠微观粒子的热运动而产生的热量传递现象，称为导热；对流则是流体各部分之间发生宏观相对位移及冷热流体的相互掺混。
联系是：在发生对流换热的同时必然伴生有导热。
导热、对流这两种热量传递方式，只有在物质存在的条件下才能实现，而辐射可以在真空中传播，辐射换热时不仅有能量的转移还伴有能量形式的转换。
2． 以热流密度表示的傅立叶定律、牛顿冷却公式及斯忒藩－玻耳兹曼定律是应当熟记的传热学公式。
试写出这三个公式并说明其中每一个符号及其意义。
答：①傅立叶定律：dx dt q λ-=，其中，q －热流密度；λ－导热系数；dx dt－沿x 方向的温度变化率，“－”表示热量传递的方向是沿着温度降低的方向。
②牛顿冷却公式：)(f w t t h q -=，其中，q －热流密度；h －表面传热系数；w t －固体表面温度；f t －流体的温度。
③斯忒藩－玻耳兹曼定律：4T q σ=，其中，q －热流密度；σ－斯忒藩－玻耳兹曼常数；T －辐射物体的热力学温度。
3． 导热系数、表面传热系数及传热系数的单位各是什么？哪些是物性参数，哪些与过程有关？答：①导热系数的单位是：W/(m.K)；②表面传热系数的单位是：W/(m2.K)；③传热系数的单位是：W/(m2.K)。
这三个参数中，只有导热系数是物性参数，其它均与过程有关。
4．当热量从壁面一侧的流体穿过壁面传给另一侧的流体时，冷、热流体之间的换热量可以通过其中任何一个环节来计算（过程是稳态的），但本章中又引入了传热方程式，并说它是“换热器热工计算的基本公式”。
试分析引入传热方程式的工程实用意义。
答：因为在许多工业换热设备中，进行热量交换的冷、热流体也常处于固体壁面的两侧，是工程技术中经常遇到的一种典型热量传递过程。
5．用铝制的水壶烧开水时，尽管炉火很旺，但水壶仍然安然无恙。
数值传热学陶文铨第四章作业
4-1解：采用区域离散方法A 时；网格划分如右图。
内点采用中心差分123278.87769.9T T T ===22d T T=0dx - 有 i+1i 122+T 0i i T T T x ---=∆将2点,3点带入321222+T 0T T T x --=∆ 即321209T T -+= 432322+T 0T T T x --=∆4321322+T 0T T T x --=∆ 即4321209T T T -+-= 边界点4（1）一阶截差 由x=11dT dx =，得 4313T T -= （2）二阶截差 11B M M q x x xT T S δδλλ-=++V所以 434111. 1.36311T T T =++即 43122293T T -= 采用区域离散方法B22d T T=0dx - 由控制容积法 0w edT dT T x dT dT ⎛⎫⎛⎫--∆= ⎪ ⎪⎝⎭⎝⎭所以代入2点4点有322121011336T T T T T ----= 即 239028T T -=544431011363T T T T T ----= 即3459902828T T T -+=对3点采用中心差分有432322+T 013T T T --=⎛⎫⎪⎝⎭即2349901919T T T -+= 对于点5 由x=11dT dx =，得 5416T T -= （1）精确解求左端点的热流密度 由 ()21x x eT e e e -=-+ 所以有 ()2220.64806911x xx x dT e e q e e dxe e λ-====-+=-=++ （2）由A 的一阶截差公式210.247730.743113x T T dT q dxλ=-=-==⨯= （3）由B 的一阶截差公式0.216400.649213x dTq dxλ=-=-== （4）由区域离散方法B 中的一阶截差公式： 210.108460.6504()B BT T dT dx x δ-⎛⎫==⨯=⎪⎝⎭ 通过对上述计算结果进行比较可得：区域离散B 有控制容积平衡法建立的离散方程与区域离散方程A 中具有二阶精度的格式精确度相当！ 4-3解：将平板沿厚度方向3等分，如图由题可知该导热过程可看作无限大平板的一维稳态有源导热问题，则控制方程为22d T+S=0dxλx=0， T 0=75℃ x=0.1 dT =h(T-T )dxf λ- 1点 ，2点采用中心差分有21022+T 0T T S x λ-+=∆ （1） 32122+T 0T T S x λ-+=∆ （2）右端点采用一阶截差的离散231f hx T T T x h λλ⎡⎤+⎢⎥⎣⎦=⎛⎫+ ⎪⎝⎭V （3）右端点采用二阶截差的离散232.1f x S hx T x T T x h λλλ⎡⎤⎢⎥++⎢⎥⎢⎥⎣⎦=⎛⎫+ ⎪⎝⎭V V V 代入（1）（2）（3）得1223132280.62 5.67625T T T T T T T -=--=-= 解得123278.87769.9T T T ===代入（4）得12380.6380.6675.1T T T === 3221T 18125T -=解得 12380.6380.6675.1T T T ===精确解 22d T+S=0dxλ （4）x=0， T 0=75℃ （5） x=0.1 dT =h(T-T )dxf λ- （6） 代入数据积分的2250025075T x x =-++ 将 x 1=10.13⨯，x 2=20.13⨯, x 3=0.1 T 1=80.56 T 2=80.56 T 3=75.1通过比较可得右端点采用二阶截差的离散更接近真实值。
传热学杨世铭 陶文铨 第四章热传导问题数值解法
(m+1,n)
y
(m,n-1)
y
x
o
第四章 导热问题的数值解法
x
x
15
以二维、稳态、有内热源的导热问题为例 此时：
Φ上 Φ下 Φ左＋Φ右 Φv 0
dt dt 左 A y dx dx
可见：当温度场还没有求出来之前，我们并不知道 dt dx 所以，必须假设相邻节点间的温度分布形式，这里我们
无内热源时变为：
x 2
Φ
4tm,n tm1,n tm1,n tm,n1 tm,n1
重要说明：所求节点的温度前的系数一定等于其他 所有相邻节点温度前的系数之和。这一结论也适用 于边界节点。但这里不包括热流(或热流密度)前的
系数。
第四章 导热问题的数值解法
19
4-2 边界节点离散方程的建立及代数 方程的求解
第四章 导热问题的数值解法
4
§4-1 导热问题数值求解的基本思想 及内部节点离散方程的建立
1 物 理 问 题 的 数 值 求 解 过 程
建立控制方程及定解条件 确定节点（区域离散化）
设立温度场的迭代初值
建立节点物理量的代数方程
求解代数方程
改进初场
是否收敛 是 解的分析
否
第四章 导热问题的数值解法
5
针对二维矩形区域内的稳态无内热源导热问题外部与温度为t的流体对流换热换热系数为h请建立外部角点的温度离散方程并化简到最后的形式inclassproblems第四章导热问题的数值解法26quickreviewquickreview1导热数值解法的重要意义2导热数值解法的基本思想3网格划分区域离散的过程及涉及的基本概念4代数方程离散方程的建立方法和过程第四章导热问题的数值解法2742代数方程组的求解代数方程组的求解一维无限大平板稳态常物性无内热源左侧第一类边条右侧第三类如右图所示将其均匀分成三个控制体试建立离散方程内部节点内部节点边界节点第四章导热问题的数值解法28形成如下代数方程组
数值传热学第四章课件陶文铨
主讲陶文铨西安交通大学能源与动力工程学院热流中心CFD-NHT-EHT CENTER 2010年9月27日, 西安数值传热学第四章扩散方程的数值解及其应用(1)4.1 一维导热问题4.1.1一维稳态导热的通用控制方程4.1.3界面导热系数的确定方法4.1.4 一维非稳态导热控制方程的离散化4.1.2通用控制方程控制容积积分法的离散4.1.5 数学上的稳定未必导致物理上有意义的解一维稳态导热问题不同坐标系通用控制方程0 P P()0P x x Δ=i调和平均已经广泛为国内外学术界所接受。
≤1数学上的稳定未必导致物理上有意义的解无内热源一维非稳态导热，初场均匀，两表面0]T +代入下式：P（全隐格式）才能满足。
结论：数学上的稳定未必导致物理上有意义的解；推=xΔa TP P极坐标均可以表示成为：2.解决通用化的一种方案为写出适合于三种坐标系中系数的通用表达式，特引进两个辅助变量：（1）x –方向标尺因子，scaling factor ，x-方向的距离表示成为sx x δi 。
对直角、圆柱坐标规定1;sx ≡（2）y-方向引入一个名义半径，R 。
对直角坐标R ＝1，据此，东西导热距离为：sx xδi 东西导热面积为：R /y sxΔ对极坐标取;sx r =对圆柱与极坐标R ＝r三种二维正交坐标系中离散方程的统一表达式按这种方式编制程序时，只要设置一个变量MODE，4.3 源项与边界条件的处理4.3.1非常数源项的线性化处理1. 线性化方法4.3.2第二、三类边界条件使方程组封闭的处理2. 线性化方法讨论3. 线性化方法应用实例1. 补充以边界节点代数方程的方法2. 附加源项法S= P2. 线性化方法讨论（1）对与被求解变量有关的非常数源项，线性化比假定为常数更合理：用*()PS f T =来表示P 的源项比落后一个迭代步；P C P T S S S =+（2）任何复杂的函数总可以用线性函数来近似逼近；线性又是建立线性代数方程所必须的；（3）是为保证代数方程迭代求解收敛所必须；0P S ≤P P nb nb a a b φφ=+∑P nb a a ≥∑P nb P a a S V =−Δ∑代数方程迭代求解收敛的充分条件是，因为可以确保代数方程迭代求解收敛。
传热学第四版课后作业答案(杨世铭-陶文铨)]
![传热学第四版课后作业答案(杨世铭-陶文铨)]](https://img.taocdn.com/s3/m/5b92aceb9b89680203d8251d.png)
1-9 一砖墙的表面积为122m ，厚为260mm ，平均导热系数为1.5W/（m.K ）。
设面向室内的表面温度为25℃，而外表面温度为-5℃，试确定次砖墙向外界散失的热量。
解：根据傅立叶定律有：WtA9.207626.05)(25125.1=--⨯⨯=∆=Φδλ1-12 在一次测定空气横向流过单根圆管的对流换热实验中，得到下列数据：管壁平均温度t w =69℃，空气温度t f =20℃，管子外径 d=14mm ，加热段长 80mm ，输入加热段的功率8.5w ，如果全部热量通过对流换热传给空气，试问此时的对流换热表面传热系数多大？ 解：根据牛顿冷却公式()f w t t rlh q -=π2所以()f w t t d qh -=π＝49.33W/(m 2.k)1-20 半径为0.5 m 的球状航天器在太空中飞行，其表面发射率为0.8。
航天器内电子元件的散热总共为175W 。
假设航天器没有从宇宙空间接受任何辐射能量，试估算其表面的平均温度。
解：电子原件的发热量＝航天器的辐射散热量即：4T Q εσ=4A QT εσ=∴ =187K 热阻分析 ;;2-4 一烘箱的炉门由两种保温材料A 及B 组成，且B A δδ2=(见附图)。
已知)./(1.0K m W A =λ,)./(06.0K m W B =λ,烘箱内空气温度4001=f t ℃，内壁面的总表面传热系数)./(501K m W h =。
为安全起见，希望烘箱炉门的 外表面温度不得高于50℃。
设可把炉门导热作为一维问题处理，试决定所需保温材料的厚度。
环境温度=2f t 25℃，外表面总传热系数)./(5.922K m W h =。
解：热损失为()()22111f f BBA A fwf t t h t t h t t q -+-=+-=λδλδ又50=fw t ℃；B A δδ=联立得m m B A 039.0;078.0==δδ2-9 双层玻璃窗系由两层厚为6mm 的玻璃及其间的空气隙所组成，空气隙厚度为8mm 。
传热学第四版课后思考题目解析（杨世铭-陶文铨）]
![传热学第四版课后思考题目解析（杨世铭-陶文铨）]](https://img.taocdn.com/s3/m/19ca3113f02d2af90242a8956bec0975f465a403.png)
传热学第四版课后思考题目解析（杨世铭-陶文铨）]第一章思考题1．试用简练的语言说明导热、对流换热及辐射换热三种热传递方式之间的联系和区别。
答：导热和对流的区别在于：物体内部依靠微观粒子的热运动而产生的热量传递现象，称为导热；对流则是流体各部分之间发生宏观相对位移及冷热流体的相互掺混。
联系是：在发生对流换热的同时必然伴生有导热。
导热、对流这两种热量传递方式，只有在物质存在的条件下才能实现，而辐射可以在真空中传播，辐射换热时不仅有能量的转移还伴有能量形式的转换。
2．以热流密度表示的傅立叶定律、牛顿冷却公式及斯忒藩－玻耳兹曼定律是应当熟记的传热学公式。
试写出这三个公式并说明其中每一个符号及其意义。
答：① 傅立叶定律：dx dtq λ-=，其中，q －热流密度；λ－导热系数；dx dt －沿x 方向的温度变化率，“－”表示热量传递的方向是沿着温度降低的方向。
② 牛顿冷却公式：)(f w t t h q -=，其中，q －热流密度；h －表面传热系数；w t －固体表面温度；f t －流体的温度。
③ 斯忒藩－玻耳兹曼定律：4T q σ=，其中，q －热流密度；σ－斯忒藩－玻耳兹曼常数；T－辐射物体的热力学温度。
3．导热系数、表面传热系数及传热系数的单位各是什么？哪些是物性参数，哪些与过程有关？答：① 导热系数的单位是：W/(m.K)；② 表面传热系数的单位是：W/(m 2.K)；③ 传热系数的单位是：W/(m 2.K)。
这三个参数中，只有导热系数是物性参数，其它均与过程有关。
4．当热量从壁面一侧的流体穿过壁面传给另一侧的流体时，冷、热流体之间的换热量可以通过其中任何一个环节来计算（过程是稳态的），但本章中又引入了传热方程式，并说它是“换热器热工计算的基本公式”。
试分析引入传热方程式的工程实用意义。
答：因为在许多工业换热设备中，进行热量交换的冷、热流体也常处于固体壁面的两侧，是工程技术中经常遇到的一种典型热量传递过程。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
4-1
解：采用区域离散方法A 时；网格划分如右图。
内点采用中心差分123278.8
77
69.9
T T T === 22d T T=0dx - 有 i+1i 122+T 0i i T T T x ---=∆ 将2点,3点带入
321222+T 0T T T x --=∆ 即321209
T T -+= 432322+T 0T T T x --=∆4321322+T 0T T T x --=∆ 即4321209T T T -+-= 边界点4
（1）一阶截差 由x=1 1dT dx =，得 4313
T T -= （2）二阶截差 11B M M q x x x T T S δδλλ
-=++V 所以 434111. 1.36311
T T T =++ 即 43122293
T T -=
采用区域离散方法B 22d T T=0dx - 由控制容积法 0w e dT dT T x dT dT ⎛⎫⎛⎫--∆= ⎪ ⎪⎝⎭⎝⎭ 所以代入2点4点有
322121011336
T T T T T ----= 即 239028
T T -= 544431011363T T T T T ----= 即 34599 02828
T T T -+=
对3点采用中心差分有 432
322+T 013T T T --=⎛⎫ ⎪⎝⎭ 即 2349901919
T T T -+= 对于点5 由x=1 1dT dx =，得 5416
T T -= （1）精确解求左端点的热流密度
由 ()21
x x e T e e e -=-+ 所以有 ()2200
20.64806911x x x x dT
e e q e e dx e e λ-====-
+=-=++ （2）由A 的一阶截差公式
（3）由B 的一阶截差公式
（4）由区域离散方法B 中的一阶截差公式：
通过对上述计算结果进行比较可得：区域离散B 有控制容积平衡法建立的离散方程与区域离散方程A 中具有二阶精度的格式精确度相当！ 4-3
解：将平板沿厚度方向3等分，如图
由题可知该导热过程可看作无限大平板的一维稳态有源导热问题，则控制方程为
x=0， T 0=75℃
x=0.1 dT =h(T-T )dx f λ-
1点 ，2点采用中心差分有
21022+T 0T T S x
λ
-+=∆ （1） 3
2122+T 0T T S x λ-+=∆ （2） 右端点采用一阶截差的离散
231f hx T T T x h λλ⎡⎤+⎢⎥⎣⎦=⎛⎫+ ⎪⎝
⎭V （3） 右端点采用二阶截差的离散
代入（1）（2）（3）得
1223132280.6
2 5.67625T T T T T T T -=--=-= 解得123278.877
69.9T T T ===
代入（4）得
解得 12380.63
80.6675.1
T T T ===
精确解 22d T +S=0dx
λ （4） x=0， T 0=75℃ （5） x=0.1 dT =h(T-T )dx
f λ- （6）
代入数据积分的
将 x 1=10.13⨯，x 2=20.13⨯, x 3=0.1
T 1=80.56 T 2=80.56 T 3=75.1
通过比较可得右端点采用二阶截差的离散更接近真实值。
4-4
解:采用区域离散方法B 进行离散，如图
控制方程为
x=0， T 0=75℃
x=0.1 dT =h(T-T )dx f λ-
对1点进行离散得1对点进行离散得32
43482.935/2T T T T T x x --=
=∆∆10
21
02
T T T T S x x x λλ---+∆=∆∆
对2点进行离散得
对右端点采用附加源法的
本题中()p w e
a a x λ
δ== C S S =
代入数据，
T 1= 82.4℃ T 2= 84.87 ℃ T 3=81.7℃
由Fourier 导热定理
得 482.935T =
4-12
function x=zhuiganfa
A=[1 2 3 4 5 6 7 8 9 10];
B=[0 1 2 3 4 5 6 7 8 9];
C=[1 2 3 4 5 6 7 8 10 0];
D=[3;11;25;45;71;103;141;185;235;190];
n=length(A);
u0=0;y0=0;B(1)=0;
%追得过程
L(1)=A(1)-B(1)*u0;
y(1)=(D(1)-y0*B(1))/L(1);u(1)=C(1)/L(1); for i=2:(n-1)
L(i)=A(i)-B(i)*u(i-1);
y(i)=(D(i)-y(i-1)*B(i))/L(i);
u(i)=C(i)/L(i);
end
L(n)=A(n)-B(n)*u(n-1);
y(n)=(D(n)-y(n-1)*B(n))/L(n);
%赶的过程
x(n)=y(n);
for i=(n-1):-1:1 x(i)=y(i)-u(i)*x(i+1); end。