数值传热学陶文铨第三章

热流问题的数值计算Numerical Simulations of Thermal & Fluid Problems第一章 绪论主讲 陶文铨西安交通大学能源与动力工程学院 热流中心 CFD-NHT-EHT CENTER 2007年10月16日, 西安1/88物理问题数值解的基本思想 把原来在空间与时间坐标中连续的物理量的场 (如速度场,温度场,浓度场等),用一系列有限 个离散点(称为节点,node)上的值的集合来代替; 通过一定的原则建立起这些离散点上变量值之间关 系的代数方程 (称为离散方程,discretizationequation);求解所建立起来的代数方程以获得所求解变量的近似解.2/88大规模科学计算的重要性 传热与流动问题数值计算是应用计算机求解热量传 递过程中的速度场,温度场等的分支学科,是大规模 科学计算的重要组成部分,其重要性不言而喻. 2005年美国总统顾问委员会向美国总统提出要大 力发展计算科学以确保美国在世界上的竞争能力. 波音公司实现了对航空发动机的网格数达10亿量 级的直接数值模拟,以研究所设计发动机的性能.3/88现代科学研究的三大基本方法及其关系理论分析Analytical实验研究Experimental数值模拟Numerical4/88课程简介1. 学时- 30学时理论教学;6学时计算机作业 2. 考核- 平时作业/计算机大作业/考试: 20/30/50 3. 方法- 理解,参与,应用 努力将与数学处理相对应的物理背景联系起来理解. 4. 助手- 于乐 5. 参考教材-《计算流体力学与传热学》,中国建筑 工业出版社,19915/88学习方法建议1. 善于从物理过程基本特性来掌握理解数值方法; 2. 对数值方法-明其全而析其微:明其全-了解基本原理;析其微-掌握实施细节;3. 努力上机实践; 4. 学会分析计算结果: 合理性,规律性; 5. 应用商业软件与自编程序相结合.6/88《热流问题的数值计算》 主要教学内容第一章 绪论(物理与数学基础) 第二章 一维导热问题的数值解 第三章 多维导热问题的数值解 第四章 势流及管道内充分发展流动与换热的数值解 第五章 有回流的动与换热问题的数值解 第六章 二维涡量-流函数法通用程序介绍 第七章 原始变量法与湍流数值模拟简介7/88绪论1.1 流动与传热问题控制方程的基本类型 1.2 流动与传热问题数值计算的基本步骤 1.3 建立离散方程的方法 1.4 离散方程数学与物理特性分析简介8/881.1 流动与传热问题控制方程的基本类型1.1.1 流动与传热问题完整的数学描写 1.1.2 控制方程 1. 质量守恒方程 3. 能量守恒方程 1.1.3 单值性条件 1.1.4 建立数学描写举例 1.1.5 控制方程式的分类9/882. 动量守恒方程1.1 流动与传热问题控制方程的基本类型1.1.1 流动与传热问题完整的数学描写 1. 有关的守恒定律的偏微分方程(控制方程)一切宏观的流动与传热问题都由三个守恒定律所 支配:质量,动量与能量守恒(conservation law).2. 与表述守恒定律的偏微分方程相关的单值性条件.不同问题的区别主要在于单值性条件 (conditions for unique solution) 的不同:初始条件以,边界条件 以及物性数据.10/881.1.2 控制方程(Governing equations) Mass conservation1. 质量守恒方程r ( r u ) ( r v) ( r w) + + + =0 t x y z单位时间 内质量的 增加 单位时间内流 进微元体的净 质量物理意义:单位时间内空 间某一微元容积质量的增 加等于流入该微元容积的 净质量.11/88对不可压缩流体: r = const 对二维不可压缩流体:u v + =0 x yu v w + + =0 x y z对二维问题,速度矢量:ur u v 数学上称: + = div(U ) x yur r ur U =ui+v j为速度矢量的散度,因此对二维不可压流体有:ur div(U ) = 0下面只讨论不可压缩流体(incompressible flow).12/882. 动量守恒方程(Momentum conservation)对上图所示的微元体分别在三个坐标方向上应用 Newton第2定律(F=ma)在流体中的表现形式: [微元体内动量的增加率]=[作用在微元体上各种力之和] 可得出三个坐标方向的动量方程:u uu uv uw 1 p 2u 2u 2u + + + =+ n ( 2 + 2 + 2 ) + Fx t x y z r x x y z 1 p v vu vv vw 2v 2v 2v + + + =+ n ( 2 + 2 + 2 ) + Fy t x y z r y x y z 1 p w wu wv ww 2 w 2 w 2 w + + + =+ n ( 2 + 2 + 2 ) + Fz t x y z r z x y z微元体内动 量的增加率压力粘性力体积力13/883. 能量守恒方程(Energy conservation)[微元体内热力学能的增加率]=[通过流动与导热进入 微元体内的净热流量]+[体积力与表面力对微元体所做 的功率] 引入导热Fourier定律,假定热物性为常数,可得T (uT ) (vT ) ( wT ) 2T 2T 2T rcp[ + + + ] = l( 2 + 2 + 2 ) + S t x y z x y z微元体 内能增 加率 由于流动被带出 微元体的净功率 由于导热而进入 源项 微元体的净功率 生成 热14/88l =a rcp流体的热扩散率(thermal diffusivity)4. 对于二维稳态对流换热问题控制方程汇总u v + =0 x yuu uv 2u 2u 1 p + =+ n ( 2 + 2 ) + Fx y z r x x yvu vv 2v 2v 1 p + =+ n ( 2 + 2 ) + Fy y z r y x y(uT ) (vT ) 2T 2T + = a( 2 + 2 ) + ST x y x y对流项扩散项源项数值计算中常用的术语.15/88不同的二维,稳态求解问题之间的区别在于: (1)边界条件不同; (2)源项与扩散系数不同.5. 二点说明1. 所导出的三维非稳态Navier-Stokes方程,无论对 层流或是湍流都是适用的. 2. 辐射换热需要用积分方程来描述,课程中将不涉及 这类问题.16/881.1.3 单值性条件 1. 初始条件 2. 边界条件 (1) 第一类 (Dirichlet):t = 0, T = f ( x, y, z )TB = Tgiven(2) 第二类 (Neumann): qB = -l (T ) B = qgiven n(3) 第三类 (Rubin):规定了边界上被求函数的一阶导数与函数之间的关系: -l ( T ) B = h(TB - T f )n数值计算中计算区域的出口边界条件常常最难 确定,要做近似处理.17/881.1.4 建立数学描写举例 1. 问题与假设条件突扩区域中的对流传热:二维,稳态,不可压缩, 常物性,不计重力与黏性耗散.18/882. 控制方程u v + =0 x y1 p u u u u u +v =+n ( 2 + 2 ) r x x y x y 2 2 v v 1 p v v u +v =+n ( 2 + 2 ) x y r y x y2 2T T T T u +v = a( 2 + 2 ) x y x y2 219/883. 边界条件 (1)进口边界条件:给定u,v,T随y 的分布; (3)中心线: u = T = 0; v = 0 y y(4)出口边y x界:数学上要 求给定u,v,T 或其导数随y 的分布;实际 上做不到;数 值上近似处理20/88(2)固体边界条件:速度无滑移,温度无跳跃1.1.5 传热与流动问题的数学描写的分类 1. 从数学角度分类-椭圆型与抛物型椭圆型 (Elliptic)椭圆型方程数学上的特点是:所求解的因变量对每个 空间自变量均存在二阶导数项: 导热方程-所求解的因变量为温度T ,空间自变量x,y; 动量方程-所求解的因变量为速度u ,空间自变量x,y.21/88抛物型(Parabolic)抛物型方程数学上的特点是:所求解的因变量对某个 个自变量只存在一阶导数项: 非稳态导热方程-因变量T 对时间t仅有一阶导数; 边界层动量方程-u对空间自变量x仅有一阶导数. 仅存在一阶导数的自变量在物理过程上的重要特 点:过程只能沿该坐标的单个方向进行而不能逆向进 行.22/88抛物型与椭圆型流动的例子椭圆型方程的求解必须全场联立进行,而抛物性 方程的求解可以沿坐标正向逐步推进, 大大节省时间.23/88(1)椭圆型问题: 流动有回流,必须 全场同时求解; (2)抛物型问题:流动无回流,可以沿主流方向步 步逼进,不必全场同时求解,大大节省时间.Marching method24/882. 从物理角度分类-守恒型与非守恒型守恒型( Conservative)-对任意大小容积守恒特性 都能得到满足的方程; 凡对流项表示成散度形式的方程具有守恒性 . 非守恒型方程+u v v u u v u ++ u = 0= 0 u ( + ) = 0 x x y y x y (uu ) (uv) 1 p 2u 2 v =+n ( 2 + 2 ) + r x x x y x守恒型方程凡是从守恒型控制方程推导得到的用于数值求解 的代数方程也具有守恒特性.25/881.2 流动与传热问题数值求解的基本步骤1.2.1 流动与传热问题数值求解步骤 1. 建立数理模型 3. 方程的离散化 5.代数方程求解 1.2.2 区域离散化方法 2.区域的离散化 4. 边界条件离散 6. 求解结果分析1.区域离散化的任务 2. 区域离散方法1.2.3 网格系统标记方法26/881) 外节点法2. 内节点法1.2.1 流动与传热问题数值求解步骤把原来在空间与时间坐标中连续的物理量的场 (如速度场,温度场,浓度场等),用一系列有限个 离散点(称为节点,node)上的值的集合来代替;通过 一定的原则建立起这些离散点上变量值之间关系的代 数方程(称为离散方程,discretization equation);求 解所建立起来的代数方程以获得所求解变量的近似解.27/88(1) 区域离散 (2) (3) (4) (5) 代数求解 (6)28/88方程离散结果分析1.2.2 区域离散化1.区域离散化的任务将所计算的区域分割成许多不重叠的子区域,确 定每个子区域中节点的位置以及所代表的控制容积. 离散结果得出四种几何要素: (1) 节点(node):所求解未知量的位置; (2) 控制容积(control volume):实施守恒定律的最 小几何单位; (3) 界面(interface):控制容积的分界位置; (4) 网格线(grid lines):沿坐标方向相邻节点连接 成的曲线簇.29/882. 区域离散方法 (a) 外节点法:节点位于子区域的角顶;控制容积界 面位于两节点之间;生成过程:先节点后界面;又 称 Practice A.子区域控制容积30/88YPractice A-外节点法 x31/88(b) 内节点法:节点位于子区域的中心;子区域即为 控制容积;生成过程:先界面,后节点,又称 Practice B.子区域即为控制容积32/88YPractice B-内节点法 x33/88 1.2.3 内接点与外节点法的比较 (a)边界节点所代表的控制容积不同 方法A 边界节点代表半个CV方法B 边界节点代表零个CV(b)网格非均分时,节点作为控制容积的代表方法B 更合理 方法A 方法B34/881.2.3 网格系统表示方法 网格线-节点间连线,用实线表示;界面为虚线; 节点间距离-dx;界面间距离-Dx .35/881.2.4 网格独立解 当网格足够细密以至于再进一步加密网格已对 数值计算结果基本上没有影响时所得到的数值解称 为网格独立解(grid-independent solution).Int. Journal Numerical Methods in Fluids, 1998, 28: 1371-1387.36/881.3 建立离散方程的方法 1.3.1 一维模型方程( 1-D model equation ) 1.3.2 由Taylor 展开法导出导数的差分表示式 1.3.3 控制容积积分法导出导数的差分表示式 1.3.4 讨论37/881.3 建立离散方程的方法 1.3.1 一维模型方程( 1-D model equation ) 一维模型方程是一维非稳态有源项的对流-扩 散方程,具有四个特征项,便于离散方法的研讨. 非守恒型 守恒型 ( rf ) f f + ru = (G ) + Sf t t x xFDM采用 ( rf ) ( r uf ) f + = (G ) + Sf FVM采用 t t x x 瞬态 对流 扩散 源项38/88"麻雀虽小,五脏俱全!"1.3.2 由Taylor 展开法导出导数的差分表示式 1. 一阶导数的差分表达式的导出 将函数f ( x, t ) 在(i+1,n)的值对(i,n)点做Taylor展开:f 2f Dx 2 2 f (i + 1, n) = f (i, n) + )i ,n Dx + 2 )i ,n Dx + ..... x x 2!f f (i + 1, n) - f (i, n) Dx 2f ) i ,n = - ( 2 )i ,n + ... x Dx 2 x39/88O ( Dx ) 称为截断误差, truncation error,表示:随 Dx 的趋于零,用 f (i + 1, n) - f (i, n) 代替 f )i ,n 的误差 x Dxf f (i + 1, n) - f (i, n) )i ,n = + O(Dx) x Dx KD x, K 与 Dx 无关.D x 的方次称为截差的阶数(order of TE).用数值计算的近似解 fin 代替精确解 f (i, n)fin 1 - fin f )i ,n @ + , O(Dx) 得向前差分: x Dx40/88f -f f )i ,n @ 向后差分: x Dxn in i -1, O (Dx )fin 1 - fin 1 f )i , n @ + , O(Dx 2 ) 中心差分: x 2Dx2. 一,二阶导数的各种差分表达式. 表达差分结构的格式图案o构筑差分表达式的位置; 构筑差分表达式所用到的节点.41/88一阶导数的 常用差分表达式42/88二阶导数的常用差分表达式定性判别导数的差分表达式正确与否的方法: (1)量纲是否正确-与导数本身一致; (2)均匀场的各阶导数应为零.43/883. 一维模型方程的有限差分显式离散表示式 微分方程形式: 假设 ( rf ) f f + ru = (G ) t t x xr , u, G均为常数,显式差分表达式:fin +1 - fin fin 1 - fin 1 r + ru + = Dt 2Dx fin 1 - 2fin + fin 1 G + , O (Dt , Dx 2 ) Dx 2差分方程 截断误差44/88显式(Explicit)-空间导数均以初 始时刻之值计算.1.3.3 控制容积积分法导出导数的差分表示式 1. 控制容积积分法实施步骤 1. 将守恒型的方程对控制容积做积分; 2. 选定被求函数及其一阶导数对时间,空间的变化 曲线-型线; 3. 完成积分,整理成相邻节点间未知量的代数方程. 2. 两种常用型线 型线-被求函数随自变量的局部变化方式,本是 所求内容,近似求解需先假定.45/88随空间自变量的变化型线 型线 型线分段线性阶梯逼近46/88piece-wise linear step-wise approximation随时间自变量的变化型线分段线性 piece-wise linear阶梯逼近 step-wise approximation47/883. 一维模型方程的控制容积积分法离散 将守恒型控制方程对控制容积P 在[t, t+ Dt ]内 做积分, ( rf ) ( r uf ) ft立即可得e+xt +Dt t=xe(Gx)r ò (ft +Dt -ft )dx +rwò [(uf)òt- (uf)w ]dt =t +Dt=Gf f [( )e - ( ) w ]dt x xf 以及 x48/88继续积分,需要知道:f对空间与时间的变化型线.1. 非稳态项假设 f 对空间呈阶梯型变化:t t r ò (f t +Dt - f t )dx = r (f P+Dt - f P )Dx w e2. 对流项假设 f 对时间呈显示阶梯型变化:rt +Dtòt[(uf )e - (uf ) w ]dt = r[(uf )te - (uf )tw ]Dt49/88假设 f 对空间呈分段线性变化:fE + fP fP + fW fE - fW r[(uf ) - (uf ) ]Dt = r uDt ( ) = r uDt 2 2 2t e t w均分网格3. 扩散项f 假设 对时间呈显式阶梯型变化: xt +DtGòtf f f t f t [( )e - ( ) w ]dt = G[( )e - ( ) w ]Dt x x x x50/88假设 f 对空间呈分段线性变化:。

2020年硕士研究生考试
同等学力加试传热学科目考试大纲
一、考查目标
按全国硕士研究生入学考试要求为沈阳建筑大学招收建筑设备与环境、供暖通风与空调专业硕士研究生而设置的专业课程考试科目。

其中，传热学是属招生学校自行命题的性质。

它的考查目标是高等学校优秀本科毕业生能达到的及格或及格以上水平，以保证被录取者具有基本的传热理论知识并有利于招生学校在专业上择优选拔。

传热学考试的目标在于考查考生对传热学基本概念、基本理论的掌握和分析求解基本问题的能力。

考生应能：
1. 准确地把握定义的物理量以及它们的量纲；
2. 正确理解基本概念和基本规律；
3. 正确应用基本理论知识分析和处理实际传热问题；
4. 掌握基本计算方法，准确完成传热问题的定量计算。

二、考试形式与试卷结构
（一）试卷满分及考试时间
传热学满分为100分，考试时间为2小时。

（二）答题方式
答题方式为闭卷、笔试。

（三）试卷内容结构。

数值传热学陶文铨课后答案第三章A.陶文铨,B.周伟俊D.黄仁龙，黄国生，周伟俊，周国生，周大白。

在计算流体力学中，流场是一个基本问题。

计算流体力学中常用到的主要方法是()。

一、对于固体，流场和速度的分布问题的计算是通过计算流体的物理量来进行的。

在计算流体力学中，流场和速度的分布问题的计算是通过计算流体的物理量来进行的。

A.质量流场 F:计算流体力学中流场和速度分布问题方法，其目的是用物理量来表示（质量）。

D.速度分布方程 F:计算流体力学中非稳态速度分布方程。

二、对于固体流动，常用的方法有流线理论法、流体压力分析法、边界层理论法、粒子质量守恒定律法、自由流动技术法、流动力学理论以及流体力学仿真技术。

解析：本题考查的是流线理论中的流场理论。

固体流动方向上分为轴流和向后流场。

在轴流流动中，固体的扩散与温度成正比；向后流一般不流动。

三、流体粘性为固液相物之间的粘滞系数，且在固态中较为稳定，因此研究固液相物之间结合力和粘滞系数问题的方法通常称为粘性流场方法。

【参考答案】 C解析：流体粘性指流体相物之间形成的粘性，根据粘性对固液相物有粘滞力或阻力，流体粘性小于固体流体之间粘滞力或阻力时固体粘滞力>固体之间粘性；流体粘性大于固体之间粘性时固体粘滞力>固体之间粘性。

固液相混是指固液相混在一起。

流体特性参数是用于定量计算和评价流体性能的重要参数。

流体特性参数如流速、温度、压力、密度、粘度等。

A选项; D选项四、流体流动分为粘性流动和非粘性流动。

A.粘性流动是指流体在流动过程中具有一定的粘性而非粘性，其流动过程主要分为四个阶段：B.非粘性流动是指流体在流动过程中不具有任何粘性，其流场主要分为湍流、不稳定、流动速度等多种形式。

C.湍流是指流体在湍流状态下不能充分发挥其流动能力，其流场主要分为两种形式：湍流和不稳定。

D.不稳定是指流体在流动过程中没有任何不稳定状态，其流场主要分为两种形式：一种是无规则流动（在流动过程中无规则流动),另一种是不规则流动时无规则流动）。

主讲陶文铨西安交通大学能源与动力工程学院热流中心CFD-NHT-EHT CENTER2007年11月20日,西安第三章多维导热问题热流问题的数值计算Numerical Simulations of Thermal & Fluid Problems=xΔP Pa T+极坐标均可以表示成为：2.解决的一种方案为写出适合于三种坐标系中系数的通用表达式，特引进两个辅助变量：（1）x –方向标尺因子，scaling factor ，x-方向的距离表示成为sx x δi 。

对直角、圆柱坐标规定1;sx ≡（2）y-方向引入一个名义半径，R 。

对直角坐标R ＝1，据此，东西导热距离为：sx xδi 东西导热面积为：R /y sxΔ对圆柱与极坐标R ＝r 。

3.1.3三种二维正交坐标系中离散方程的统一表达式按这种方式编制程序时，只要设置一个变量MODE，它按以下方式取值，程序即可自动处理三种坐标系。

3.2附加源项法3.2.1第二、三类边界条件的处理方法1. 补充以边界节点代数方程的方法2. 附加源项法3.2.2附加源项法的实施细则1. 处理第二类边界条件的附加源项法2. 处理第三类边界条件的附加源项法3. 附加源项法的实施步骤3.2.3 附加源项法与补充节点法的对比对第二，第三类边界条件问题，边界温度未知，x缺点：对多维问题显著增加了求解的节点数目。

对20X20的二维区域，增加了内点的23％。

寻找不增加结点数目而又能使内节点的代数方程封闭的方法对多维问题具有重要意义。

2. 附加源项法（additional source term method,ASTM)区域的热量折算成与边界相邻的第一个控制容积的源将由第二类、第三类边界条件所规定的进入计算项；切断内点与边界点的联系，从而将未知的边界点温度从内点离散方程中排除。

W =B q y VΔΔfTC S V Δ(2)令该边界上的导热系数为零；(3)按常规方法建立内接点的离散方程，并在内接点区域求解方程组；(4)在内点区域求解方程组; 获得收敛解后按Newton 冷却公式或Fourier 定律确定边界温度。

第三章思考题1. 试说明集总参数法的物理概念及数学处理的特点答：当内外热阻之比趋于零时，影响换热的主要环节是在边界上的换热能力。

而内部由于热阻很小而温度趋于均匀，以至于不需要关心温度在空间的分布，温度只是时间的函数， 数学描述上由偏微分方程转化为常微分方程、大大降低了求解难度。

2. 在用热电偶测定气流的非稳态温度场时,怎么才能改善热电偶的温度响应特性?答：要改善热电偶的温度响应特性，即最大限度降低热电偶的时间常数hA cvc ρτ=,形状上要降低体面比，要选择热容小的材料，要强化热电偶表面的对流换热。

3. 试说明”无限大平板”物理概念,并举出一二个可以按无限大平板处理的非稳态导热问题 答；所谓“无限大”平板，是指其长宽尺度远大于其厚度，从边缘交换的热量可以忽略 不计，当平板两侧换热均匀时，热量只垂直于板面方向流动。

如薄板两侧均匀加热或冷却、 炉墙或冷库的保温层导热等情况可以按无限大平板处理。

4. 什么叫非稳态导热的正规状态或充分发展阶段?这一阶段在物理过程及数学处理上都有些什么特点?答：非稳态导热过程进行到一定程度,初始温度分布的影响就会消失，虽然各点温度仍 随时间变化,但过余温度的比值已与时间无关，只是几何位置(δ/x )和边界条件(Bi 数) 的函数，亦即无量纲温度分布不变，这一阶段称为正规状况阶段或充分发展阶段。

这一阶段的数学处理十分便利，温度分布计算只需取无穷级数的首项进行计算。

5. 有人认为,当非稳态导热过程经历时间很长时,采用图3-7记算所得的结果是错误的.理由是： 这个图表明,物体中各点的过余温度的比值与几何位置及Bi 有关,而与时间无关.但当时间趋于无限大时,物体中各点的温度应趋近流体温度,所以两者是有矛盾的。

你是否同意这种看法,说明你的理由。

答：我不同意这种看法，因为随着时间的推移，虽然物体中各点过余温度的比值不变 但各点温度的绝对值在无限接近。

这与物体中各点温度趋近流体温度的事实并不矛盾。

数值传热学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精确解与数值解的对比图，其中边界条件给定，。

课程编号：S201E045数值传热学学时：32 学分： 2 开课时间：春季授课单位：机械工程任课教师：李汛一、课程内容简介本教程的目的是使学生掌握一种能够预测传热与传质、流体流动过程的数值方法，用以解决工程实际中大量存在的，而用解析方法难以解决的传热与流体流动问题。

本课程的特点是：它强烈地以物理上的依据为基础，而不只是以数学推演为基础。

学生在学习数值方法的同时，可以加深对基本物理过程的认识和理解。

物理的手段将使学生掌握通用的评定准则。

他们应用这些准则就可以对现有的以及未来的数值方法做出评判。

课程由三部分组成，每部分分成三章，共九章。

前三章是预备性的知识，其中包括对数学与数值方法的基本讨论。

此外，这一部分还概述了本课程所特有的方法。

第四到第六章包含着数值方法的主要推演。

最后的三章则致力于解释和应用。

第一章是绪论。

第二章中概述并讨论了有关的物理现象和微分方程。

在这一章中尤为重要的是，从物理意义的观点来分析这些方程的抛物型或椭圆型特性。

在第三章中提出了数值解的概念，其中描述了构成数值方法的一般步骤。

以这一部分内容为基础，本书进一步系统地提出了形式为四项基本法则的一般准则。

这些基本法则构成本书其余部分中数值方法推演过程的准绳。

数值方法的构成开始于第四章，它通过三个步骤进行。

第四章处理热传导。

第五章集中讨论对流与热传导的相互作用；这时认为：流场是已知的。

最后，在第六章中处理速度场本身的计算。

最后三章用于对前面的方法的解释和应用。

This course is primarily aimed at to help the students to mastery a numerical method which is able to predict the process of heat and mass transfer and fluid flow and with which to solve the complex engineering problems, to which the analytical method can find no way out.An important characteristic of the numerical method to be developed here is that they strongly based on the physical consideration, not just on mathematical manipulation. A significant of this strategy is that the student, while learning about the numerical method, develops a deeper understanding of, and insight into, the underlying physical process. Further, the physical approach will equip the students with general criteria with which to judge other existing and future numerical methods.The course is consisted of three different parts with three chapters each. The first three chapters constitute the preparatory phase. Here,a preliminary discussion about the mathematical and numerical aspect is included, and the particular philosophy of the course is outlined. Chapters 4-6 contain the main development of the numerical method. The last three chapters are devoted to elucidations and applications.The first chapter is introduction. In Chapter 2, the related physical phenomena and differential equations are outlined and discussed. Of special importance in that chapter is the examination of the parabolic or elliptic nature of these equations from a physically meaningful viewpoint.The concept of numerical solution is developed in Chapter 3, where the common procedures of constructing numerical methods are described. This introductory material is used to formulate general criteria in the form of four basic rules. These rules form the guideposts for the development of the numerical method in the rest of this book.The construction of the numerical method begins in Chapter 4. It is carried out in three stages. Heat conduction (i.e., the general problem without the convection term) is treated in Chapter 4. Chapter 5 concentrates on the interaction of convection and conduction with the flow field regarded as given. Finally, the calculation of the velocity field itself is dealt with in Chapter 6.As to the last three chapters, as mentioned above, are devoted to elucidations and applications of the previous methods.二、先修课程流体力学、传热学三、教材Numerical Heat Transfer and Fluid Flow, by S.V. Patankar, McGraw-Hill Company, 1980四、主要参考书目及文献1、《数值传热学》（第二版），陶文铨著，西安交通大学出版社，20012、《计算传热学的近代进展》，陶文铨著，科学出版社，2000在燃烧问题中，高温气流和与之相邻的液体或固体物质之间存在着一个相分界面。