热物理过程的数值模拟Numerical Simulation of Thermophysics Process
讲稿
主讲：李隆键
第一章概论
1．1流动与传热过程的予测方法及特点
流动、传热、燃烧问题是热工类各专业和机械类动力机械专业所研究和解决的主要问题之一，燃烧问题实际上是有化学反应的流动与传热问题，推而广之，在所有热物理过程中，几乎都涉及到流动、传热问题。
预测的重要性：
①在规定设计参数的相应的结构下，热物理过程是否满足要求，达到预定的指
标？要预测；
②优化设计，不同方案的比较，要预测；
③减少设计、生产、再设计和再生产的费用；
④减少设计更改；
⑤减少试验和测量次数。
问题的核心：速度场、温度场（传热量）、浓度场等。
一、热物理问题的予测方法：理论分析法、实验测定、数值模拟
1、理论分析
以数学分析为基础，求解描述热物理过程的定解问题，获得函数形式的解，表示求解区域内物理量连续分布的场（速度场、温度场、浓度场……）。
控制方程+单值条件（数学模型）→理论解（分析解，解析解）
根据解的准确程度，又可再分为：
（1）精确分析解（严格解）
特点：函数形式的解；它在求解区域精确地满足定解问题。
具体解法：直接积分法、分离变量法、积分变换法、热源法、映射法。
（2）近似分析解法
特点：函数形式的解，在求解区域上近似地满足定解问题（但在总量上满足相应的守恒原理，动量守恒、动量守恒、能量守恒、质量守恒）。
具体解法：积分法（从积分方程出发）
变分近似解法
摄动法（从微分方程出发）
2、实验测定
（1）纯实验法
（2）相似理论实验法：同类相似，减少变量数目→减少工作量，得到规律性结
果，可直接应用。
（3）实验类比法：异类相似—物理现象不同，规律相同：微分方程形式相同，单值性条件类似
电热类比，水热类比……
3、数值模拟
以数值计算方法为基础，借助（利用）电子计算机求解物理过程的方法—热物理过程的数值模拟，对传热过程称为传热的数值模拟、数值传热、计算传热。如前述，传热过程函盖了流动、燃烧，所以计算传热学实质上就代表了热物理过理过程的数值模拟。
用电子计算机对热物理问题进行数值计算就象在实验室中对该现象进行实验测定一样，可称之为“数值实验”。随着高速、大容量电子计算机的发展，特别是微型计算机的普及和推广，这种数值实验的方法越来越被更多的科技人员掌握和应用，成为解决热物理过程的一种重要方法。
二、予测方法的比较与选择
1、分析解法
（1）精确预测了数学模型所控制的的热物理过程；
（2）函数形式的解使得可以确定区域中任意位置物理量的大小；
（3）以显函数的形式，展示各有关参量对该热物理过程的影响；
（4）由于是函数形式的解，便于进一步的运算、处理，例如求导、积分。
缺点：
（1）获得分析解的可能较小；
（2）即使能求得分析解，也常常是无究级数，特殊函数以及涉及特征值问题的超越函数，要得到具体的数值结果，也需要繁复的计算；
（3）数学模型的结果也需要有实验检验。
2、实验方法
（1）可以获得热物理过程可靠的数据资料；
（2）全比例设备实验可予测由它完全复制的同类设备在相同条件下将如何运行和变化；
（3）是研究一种新的基本现象的唯一方法；
（4）是检验其它预测方法准确程度的标准。
缺点：
（1）全比例实验代价大（投资，物力，人力，周期……）；
（2）缩小比例模型实验→结果的外推受准则数实验范围的限制，有些在全比例设备上才能出现的特征在缩小比例模型上并非总是能模拟（例如流动的涡），降低了模型试验的效果；
（3）测试困难及测量误差；
（4）有些过程无法预先进行试验（航天，气象预报……）。
3、数值模拟
（1）成本低：在大多数实际应用中，计算机运算的成本要比相应的实验研究成本低好几个数量级，对象愈庞大，过程愈复杂，此优点愈突出；同时，与大多数物品价格不断上涨的趋势相反，计算成本还会降低；
（2）速度快，周期短；不同方案的对比计算和优选，这对某些大型实验几乎是不可能的。
（3）信息完整：能提供计算区域内所有各个位置上有关变量的值（速度、压力、温度、浓度等），而实验则不可能测出整个区域各点处所有变量的值。
（4）具有模拟真实条件的能力：几何条件、边界条件、物性条件、初始条件……
很容易模拟真实条件，不需要采用缩小模型或冷态实验，无论大小、高位，低温、过程快、慢。
（5）具有模拟理想条件的能力：对于研究物理现象而不是工程问题时，注意力集中几个基本参数而要设法消除所有无关的因素。
几何条件（维数变化，尺寸→ ）、物性（常密度），BC（绝热表面……），ic（特定的初始温度分布……）。
缺点：
（1）数值模拟的对象是数学模型→简化处理，结果的准确性有特价检验；
（2）对一些十分复杂的问题（几何形状复杂，强烈非线性、物性变化大），数值解可能很难获得，或者即便可以获得，代价也是相当昂贵的，例如，对湍流问题，要想通过求解非稳态N-S方程来算出它们的全部与时间相关的结构，则仍然是计算所不能及的；
（3）对解的唯一判断力较弱
为了进一步讨论数值模拟的缺点，可以把所有的实际问题分成两大类：
A类：有完整数学模型的一类问题，如热传导、层流问题、简单的湍流边界层问题；
B类：迄今无完整数学模型的一类问题，如复杂湍流、某些非牛顿流体、某些两相流动等，问题的分类还有一标准问题，即描述到什么样的程度可以认为是“够了”，“合适”的。
A类缺点：对这类问题，用计算机求解的优越性远远大于实验研究。
有数值模拟缺点的（2）、（3）两条；在某些情况下也需要进行实验检
验。对于此类问题，研究计算方法的目的在于使这些计算方法理加可
靠、准确和有效，随着研究的进展，其缺点将被不断克服。
B类缺点：A类的缺点B类全有，此外，必须进行实验检验。
数学模型的研究不断地把B类问题转化为A类：试算与修正。
先提出一个模型→计算求解→与实验结果进行比较→修正模型，并不断完善、湍流模型的最新发展就是这种转换的一个典型例子，k-ε双方程模型最初建立在科尔莫戈洛夫（Kolmgorov）(1942)及普朗特（Prandtl）(1945)的工作基础上的，但并未，也不可能付诸实现，只有到了20世纪70年代，当计算机和计算方法变得更加强有力的是候，该模型才逐步趋于完整并付诸实际应用。
4、方法的选择
三种方法或三种手段相辅相成，互为补充。
（1）分析解可以为检验数值模拟结果的准确度提供比较依据；常常用有分析解的简单问题检验方法的准确度；
（2）简单的解极解可以为发展数值方法中的某些算法提供理论依据，λ→调和平均；
（3）物理规律、数学模型的正确建立必须通过对现象的充分观察和测定，-
=λ等；
?
n
tl
q?
（4）出现在数学模型中的物性参数只有通过实验测定才能获得。
数值模拟的对象是热物理过程的数学模型，所以其结果的准确度首先取决于数学模型反映实际热物理过程的准确度（包括所用的特性参数），然后才是所采用的数值方法，计算机并不能创造信息，发现规律，它只能把人们所送入的信息，按照计算者安排的程式对信息进行加工、处理，从而得到相应的结果。
但是，一旦确立了与实际物理过程相符合的物理模型、数学模型、数值模拟又可以发挥很大的作用，它可以减少实验工作量，拓宽实验研究的范围，实现对理想单值性条件的模拟；对那些耗资巨大，条件恶劣的实验，或者难于进行的实验来说，“数值实验”更是一种有吸引力的辅助或替代手段。
理论分析，实验测定和数值模拟有机而协调地结合，是研究热物理过程理想而有效的方法。
1．2本课程的内容及安排
一、内容
两个组成部分
1.理论部分：
①基础理论（数学物理、数值方法）
②热物理过程的数值模拟：通用性，并以热传导、对流换热、流场的计算为例，
推广到通用控制方程所描述的其它现象。
论述方式的特点：（1）强调物理的概念和方法，而不过分倚重纯数学的推导。（2）以一维为基础，推广（扩展）到二维、三维。
2、实践性环节：一个课程设计，程序设计，视情况而定。
二、授课方式
改变注入式，实行启发式，培养自学能力。
三、教材、参考书
1、S．V帕坦卡著、张政译（郭宽良译），“传热与流体流动的数值计算”。（Numerical Heat Transfer and Fouid Flow），科学出版社，1984。
2，陶文铨编著，数值传热学，西安交通学大出版社，1988。
第二章 物理现象的数学描述
数值计算的对象—过程的数学模型，核心是控制方程。 （数学模型：控制方程+单值性条件，？单值性条件）
2．1控制微分方程
1、控制微分方程的意义
控制微分方程是一定守恒原理的数学表达式，影响因变量的各因素之间必定存在某种联系。回忆导热方程：热力学第一定律，对任意控制容积V ，导入控制容积的热流量+控制容积中的内发热量=控制容积中物质能量的增量。
导入控制容积的净热流量?-=A
dA n q
.，
控制容积内的发热量=?v
v v d q ，
控制容积中物质内能的增量=???
v edv ρτ
。
??????+??-=??
+?-=??
v v v v v v v A dv
q dv q dv e dv
q dA n q dv e )()(ρτρτ
散度定理：?????==?v
v
v
dv q dv q diV dA n q
v v v v q q e d q q e +?-?=??
→=-??+??? )(0])([
ρτρτ 非稳态项 扩散项 源项
-??--=→dv t p T p dt c de t e v v ])/([:比容，对常密度（不可压缩）
dt
c de dv v == ,0v v q q t c +?-?=??
)(ρτ
ω
物质运动（流动）时，增加：物体宏观运动带入的能量=?
?-A
dA n w e ρ 压缩机械动 =???-v
dv w p .
粘性摩擦功 =?v
dv φ
2))(3
2(2w ij ij ??-'+=μμεμεφ
ij ε—流体的变形率张量
μ'—第二动力粘度，体积粘度
Φ++?-??-=??+??
v q q w p w e e .)()(ρρτ
Φ++??-??-=?+??+??+??v q p w p e w w e e e )()(ρρτ
ρ
τρ
S q q p w p w e e w e v +?-?=Φ++??-??-=??+??+??+?? )]([)(ρτ
ρ
τρ
φτ
ρ++??-?-=v q q w p D De
. 对简单可压缩系统
τ
ρ
ρτρττρρD d p D Dp D Dh D De p h e p e pv e h 211//,/,/+-=-=+=+= ∴ φτ
ρ
ρττρ
++??-??-=+-v q q w p D D p D Dp D Dh φρτ
ρ
ρττρ
++?-?+-=v q q w D D p D Dp D Dh .).( =??+??+??+??→++??-=)(z
h
w y h v x h u h q q D Dp D Dh v τρφττρ h w h ??+?? ρτρ
=)]([w h h w t h ρτ
ρ
ρρ
??+??+??+?? =)(w h h w h h ρρτ
ρ
τρ??+??+??+??
=)()
(h w h ρτ
ρ??+??
φτ
ρτρ++??-=+??∴
v q q D Dp
h w div h )()( 非稳 对流 扩散
对理想气体及恒密度（固液）dp T v v d c dh p ][β-+=， dt c dh p =，对固、液，忽略p 变化，0=dp
则)()()(h c div h c q p
p ?=???=???=??-λ
λτλ
各项则代表着各因素在单位容积时的作用效果.
τ
ρ??)
(h 单位容积焓的变化率
h w h w div ρρ)(单位时间、单位面积、传递的焓。对流流量密度J C 单位容积流出的净焓量s m J z
J y J x J divJ y x c ?++??=
32
/2222 d p
p
J h c h c div =?-
?λ
λ
),(
扩散流量密度)(d J div -
τφD DP q S v /++=… 源项s m J m W ?=33//
2、化学组分方程 l M M m l l /=的质量分量Mi M l ∑/
l d e l R J div m w div m +-=+??
)()()( ρρτ
扩散流密：law Fick m J l l d ?Γ-=
l l l l l R m div m w div m +?Γ=+??
)()()( ρρτ
3、能量方程
S h c div h w div h p
+?=+??)()()(λ
ρρτ 如果p c 为常数，OK 时h=0，则由T c h dt c dh p p =→=
S t div t w c div t c p p +?=+??
)()()(λρρτ
p
p c s
t c div t w div t +
?=+??)()()(λρρτ
4、动量方程：
某方向上的动量变化率u x F i -∑=，u ρ…单位体积质量x 方向的动量，uu ρ
x x v B x p u u div u div u ++??-?=+??
/)()()(ρωρτ
，其它粘性力项，u u j d ρρ
μ
μ?-
=?-=
5、湍流的时间平均方程
湍流：给定点处的物理量随时间而变化，随机性。
工程上：关心的是运动状态的时间平均特性…非稳定层流方程
的时间平均方程，滞流流动的时间平均方程，假设湍流中存在相对于时均值的脉动，平均化（时均化）运算→附加项，雷诺应力（湍流热流），湍流扩散流量密度…湍流模型，用平均性质来表示的附加项。
引入湍流粘度（湍流扩散系数）→滞流应力，流量必度→时均方程—层流方程，相应的层流交换系数用有效系数代替。
ε-R 模型
ρερρτ-+?Γ=+??
G k div k w div k k )()()(
6、通用方程 （1）向量形式
S div w div +?Γ=+??
)()()(φφρρφτ
对于不同的φ，有特定的Γ、S 与之间相对应。
大多数的传递过程，扩散流量密度Jd 因变量φ的梯度确定，φ?Γ-=d J 也有扩散流量密度不由φ梯度支配的情形，这时可将其（J d ）并入S 中。
连续性方程：
0,0,10)(==Γ=?=+??S w div φρτ
ρ
（2）直角张量形式
Φ
S x x u x j j j j +Γ=+??)22(22)(22)(φφρρφτ
0)(=??+??j j
u x ρτρ 和法则：一项某下标重复，则该下标依次取1、2、3然后做和运算。 如何把特定的控制微分方程改写成通用形式。
①把相关因变量的非稳态项，对流项及扩散项转化成标准形式； ②把扩散项内φ梯度前的系数取为Γ的表达式；
③把其余所有项之和置于方程右端，定义为源项。
通用方程可以是有量纲形式，也可以是无量纲形式，Γ、S 也相应无量纲化。 作用：通用方程→通用数值方法公式→通用计算机程度
2-2 坐标的性质及控制方程的类型
一、坐标的性质
坐标→自变量。坐标（系）或自变量的数目对问题的难易程度有很大影响，而对一定的问题而言，坐标（自变量）数目是可变的，既可以用这个坐标系来描述，又可以用别的坐标系来描述。
1、自变量的作用：),,,(),,,(τφφτz y x z y x t t ==或 ①在数值计算中，将选择用来计算φ值的自变量值，所需计算φ值的位置的多少与自变量的数目有关，自变量
数目↑所需计算φ的位置↑。
②因变量与自变量的相对性
),,,(),,,(t y x z z z y x t t ττ=?=，适用于温度场是坐标的单调函数的情形。
2、坐标的选择，原则：恰当，合适。自变量数目最少的坐标系→网格节点数少。 （1）选择最简单的坐标系
圆柱（简）中的轴对称导热：)(),(r t y x t ? 圆管内的轴对称流动：),(),,(z r t z y x ?φ
（2）运动坐标系→准稳态的概念，),,(),,,(z y x z y x '''?φτφ
)
（3）利用充分发展的概念：
存在这样一个坐标，当过程发展到一定深度后，因变量的无量纲分布与该坐标无关。
),(y x t t =，充分发展后
)(),(2
x t t y t t t t b b w
b ==--=
θθ 平面自由射流；
)(),(,/,/),,(x x u u y u u u y x u u c c c δδδη===≡=
但 )(ηu u =
（4）相似变换：减少自变量数目的变换统称相似变换。
半无限大物体（x ≥0）在1B.C.下的非稳态导热：τηη/),(),(cx t y x t =?
3、坐标的单、双向性
采用单向坐标和双向坐标的概念，可以形象地描绘不同类型控制方程的物理作用上的区别。单向坐标：在一个坐标轴上，如果扰动（或影响）只能向一个方向传递，则称此坐标为单向坐标。何谓向一个方
向传递？坐标上任意给定位置处因变量之值只受该位置一侧条件变化的影响，且该点因变量之值也只对其一侧位置上的固变量值发生影响。
（“抛物型”表示一种“单向作用”的概念）时间坐标是一个典型的单向坐标。 高温固体的冷却，其在某瞬时的温度只受该瞬时以前的条件的影响。
双向坐标：在一个坐标轴上扰动（或影响）可以向两侧传递，称为双向坐标。 何谓向两个方向传递？坐标上任意给定位置处因变量之值要受该位置两侧条件变化的影响，且该处因变量之值也会对其两侧位置上的因变量值发生影响。空间坐标是典型的双向坐标。
但在一定条件下，空间坐标也可以成为单向坐标；
在流体流动时，如果在某一个坐标上有很强的
单向流动，则重要的影响只能从上游传播到下游，
某处状态也主要受其上游条件的影响，受下游条件影响很小。
总结：对流是一种单向过程，而扩散是一种双向过程：在流体流动时，二者同时存在，仅当流作用很强（流量很大）时，扩散作用可忽略不计，空间坐标才近似成为单向坐标。
x
τ
τ=τ2
τ=τ1
τ=0
x
L 0 X
二、控制方程的类型
讨论控制方程的类型对数值解的影响
1、能量方程
导热：0)()(==+????+????v q S y
t
y x t x λλ S x t
x t c
+????=??)(λτρ 对流换热：
S h c div h w div h p
+?=+??)()()(λ
ρτρ （已假定：dt c dh p =） 若p c ≈常数，OK 时h=0，则T c h p = S t div t w c div t c p p +?=+??
)()()(λρρτ
p p p p c S
t c div w t t w t t c S t c div t w div t +?=??+??+??+??→+?=+??)()()()()(λρρτρτρλρρτ
p p p p c S z t c z y t c y x t c x wt z vt y ut x t +????+????+????=??+??+??+??)()()()()()()(λλλρρρρτp
p p c S
z t c y x t c x z t w y t v x t u t +
????+????=??+??+??+??)()(λλρρρτρ
组成：非稳态项，对流项，扩散项，源项 控制方程：微分容积中守恒原理的数学描述。
2、能量方程的守恒性质
数值计算中的控制容积无论多么小，总是有限容积，而不是微分容积。 问题：对任意有限容积，控制方程是否也一定满足守恒原理。
守恒型控制方程，对任意大小的有限容积能使守恒原理得到满足的控制方程。 非 不能
可以证明，在直角坐标中，当对流项为散度形式时，控制方程
是守恒的。
S t div t w c div t c t p p +?=+??)()()(λρρ 将上面的控制方程对任一有限容积V 作积分
ω
????+?=+??
v v v p v p Sdv dv t div dv t w c div dv t c t
)()()(λρρ 有限容积V 上的能量平衡原理的表达式： ①单位时间内有限容积中物质内能的增加量
②=-?A
p dA n w c
.ρ对流带入的净能量率
③=?-=????A
A
dA n q dA n t
)(λ导入的净能量率
④生成的能量率
若用非守恒型控制方程，有：
?????+?=??+??+??v v v p v p v p Sdv t div dv z t w c dv x t u c dv t c )(λρρτρ
??????+?=??+??+??++??v v v v p p v p v
p Sdv dv t div dv z t w c dv y t v c dv y t u c dv c t )()()()()(λρρρρ 简证：取平行六面体 ?
→??v
p dv x
t
u c )(ρ
3、控制方程的分类
（1）从数学上分类，限于二阶偏微分方程，对二元二阶线性偏微分方程，),(y x φ
),(y x g f e d c b a y x yy xy xx =+++++φφφφφφ
下标xy 表示φ对该变量的导数，,/,/,/22y x xx y xx x ??=??=?=φφφφφφ
y x y xy yy ???=??=/,/222φφφφ，
a 、
b 、
c 、
d 、
e 、
f 、是x 、y 的函数，对求解区域R ，方程的特性由系数ac b 42-之值决定
042>-ac b ，在区域内任意一点有两条实的特征线 双曲型 042=-ac b ，在区域内任意一点有一条实的特征线 抛物型 042<-ac b ，在区域内任意一点无实的特征线 椭圆型
三类方程在数学上的主要区别：影响区域和依赖区区域不相同。 影响区域、依赖区域；
R 中任一点P 的依赖区域，为了唯一定)(p φ之道，必须给函数φ值的点B 依赖区
'
x
域的集合；
P 的影响区域：)(P φ变化时，函数φ值发生变化的点的
集合。
椭圆型方程（即求解区域上每一点都是椭圆型的）：对于能量方程—无非稳态项，自变量与时间无关。
特点：无特征线，任一点P 的依赖区域是包围该点的区
域封闭边界曲线，而P 点的影响区域则是整个求解区域。
→边值问题，对应于物理学上的一类平衡问题，或稳态
问题。
→求解区域内各点处的因变量值是相互影响的，与双向坐标相联系。 结果：离散代数方程必须联立求解（直接解法或迭代解法），而不可能把区域中某一部分的值求得后再去确定其区域上的值。
抛物型方程：因变量与时间有关，或问题中存在类拟于时间的自变量，与单向坐标相联系，能量方程→有非稳态项。
对应于：物理学上一类非平衡问题，或非稳态问题，步进问题。 特点：过区域中任一点P 有一条实特征线，其方向
与单向坐标相垂直，如图，P 点的依赖区域和影响区域
以特征线为分界线。 对非稳态问题：某一瞬间物体中的温度（分布）取
决于该瞬时以前的情况及边界条件，而与该瞬时以后将要发生的情况无关，反之，某一时刻的温度只影响此后
的温度分布
边界层类型的流动与换热：忽略主流方向的扩散作
用，下游的物理量（u,v,t …）取决于上游，上游 只会影响下游的物理量。 结果：不必将单向坐标上所有位置处的离散方程联
立求解！只需从某一初始值出发，结合边界条件，沿单向坐标一步步向前推进→步进法→步进问题。大量节省计算机的内存及计算时间，二维→存储一维，三维→存储二维。
双曲型方程：依赖和影响区域与抛物型方程相同，即：依赖区域位于运动方向的上游，影响区域位于下游。
不同之处：某点的依赖区域或影响区域，不是其上游或下游区域的全部，而只是过该点的两条特征线之间的区域，如图。
例：无粘流体的非稳态流动，无粘流体的稳态超音速流动。 大多数工程导热、对流换热问题都属于椭圆型或抛物型方程，本课程只限于讨论这两类方程。 （2）从物理现象上划分
在传热相应于椭圆型与抛物型方程的物理过程，有专门的称谓。
抛物型→有一个空间坐标是单向坐标→边界层型（流动或换热）问题； 椭圆型→所有空间坐标都是双向的→回流型（ ）问题。
P ·
R
P 影响区域
x (y)
P τ
x 依赖区域 影响区域 (y)
x 0 x
P
P
影响区 特征线
依赖区
B
三、一维模型方程
控制方程由四部分组成：非稳态项、对流项、扩散项、源项。
在研究建立离散方程和方法时，为了避免复杂化，不必着眼于完全的方程，只需把同一类型的项取出一项作为代表来分析，得到一个维非稳态对流扩散方程，即在对流换热问题数值计算研究中广泛采用的所谓一维模型方程。
非守恒型：S x x x u +??Γ??=??+??)(φφρτφρ 守恒型：S x
x x u +??Γ??=??+??)()()(φ
φρτρφ φ—广义变量：u 、v 、w 、h 、c 等
Γ—广义扩散系数，对于传热c /λ=Γ
S —广义源项，代表一切不能归并到非稳、对流，扩散项中去的量，不定是物理上的真正源项。
对于导热问题：S x
h
c x h +????=??)()(λτρ C=常数，由cdT cdT dh →=
S x t x t c S x t x t c +????=????→?+????=??)()()(λτρλτρ常数 对流换热
S x
h
c x x uh t h +????=??+??)()()(λρρ C=常数，h=cT
S x
t c x x ut t S x t x x ut c t c '+????=??+??→+????=??+??)()()()()()(λρτρλρτρ
第三章 离散化方法
3-1传热问题数值求解的基本步骤
一、数值解法的基本思想
数值解法是一种离散近似的计算方法。这种方法得到的是求解区域中某些代表性位置上未知（待求）物理量（速度、温度、浓度等）的近似值，而不是象分析解或近似分析解那样的连续
函数，“数值解法”一词即由此而得名。 例：大平板（-L ≤x ≤L ）在第三类B.C.下的一维非稳态导热
f i i f t t t x t -=-=θτθ,),(
2
21x a ??=??θ
τθ B.C.：)0,0,(,=??==??±=x
x or x L x θαθλ
i.C.：i θθτ==,0
用分离变量法求解可得：
).cos()cos()sin()sin(])(exp[21
2L x
L L L L L Fo L m m m m m m m i ββββββθθ+-=∑∞
=
令)cos(cos sin sin )
exp(21
2L
x
Fo L m m m m m m m m m ξξξξξξβξ∑∞
=+-→==
式中，2L
a Fo τ
=
，m m ξβ,分别是下列（特征）方程的根：)/()(L Bi L tg m m ββ=， 21,/,/,/)(、m Bi ctg Bi tg Bi L L ctg m m m m m m ====ξξξξββ
λα/L Bi =毕渥数，物理意义，内部导热热阻/外部换热热阻
对不同的Bi 值，有不同的m ξ值系列。
2594
.162142.132003.102281.73058.44089.110
7713.156453.125293.94373.64256.38603.017143
.155743.124354.42991.61731.33111.00016
5
4321ξξξξξξBi
结论：
采用计算机进行数值计算不仅是求解
t
t i
t (X,τ)
L
-L
x
λ,a t t f
h h 0
c tg ξ
0 ξ/Bi
3π
ξ
3
ξ2
ξ1
2π 1π
偏微分方程的有力工具，而且对一些经验公式和用无穷级数表示的分析解，也常常需要用计算机来获得数值结果。
如何进行数值计算，基本思想？
把原来在时间，空间坐标上连续的物理量场（速度场，温度场，浓度场等），用求解区域中有限个离散点上的值的集合来代替，并按一定方式建立起关于这些值的代数方程，求解代数方程以获得物理量场的离散近拟解。
二、基本步骤 例：长方柱体中的导热，已知：稳态，四个侧表面各自维持均匀温度，用数值方法求柱体中的温度分布。 1．建立物理模型，对问题作必要的简化。 (1) 假定涉及的温度变化范围不大，常物性； (2) 柱体长度方向端部效应忽略不计，
→==??),(,0/y x t t z t 二维问题
稳态：二维、常物性、无内热源的导热问题，1st B.C.。
2. 建立相应的数学模型：控制方程+单值性条件 0//2222=??+??y t x t
2111, ;,:..t t L x t t L x C B =-===
4232, ;,t t L y t t L y =-===
有的工程实际问题的模型化工作较难，需要对实际问题进行仔细分析，还需要经验。
例：离散电阻片热源
t sat .a
q w
特点：温度分布在一定的深度的不均匀，畸变温度分布的不均匀性呈现周期性。所以只研究一个周期区间即可！
3. 区域离散化：在计算区域中配置需要计算温度的地点（位置）（称为节点）———这一步骤称为区域离散化。
4. 控制方程的离散化：按照一定的原则（通常是守恒原理），建立每个节点上未
x
知量（ ，温度）与其邻点上未知量之间的代数关系式（称为离散方程）——
—控制方程的离散化。例如，由能量方程可以得出每一节点温度与相邻节点温度之间的代数关系。
5.求解所得到的代数方程，获得节点上的未知量之值。
6.对所获得的数值结果进行分析、比较和讨论。
基本步骤框图（流程图）：
三、传热问题的主要数值方法
1、有限差分法（Finte Difference Method, FDM）
2、有限元法（Finite Element method, FEM）
3、边界元法（Boundary Element Method, BEM）
4、有限分析法（Finite Analysis Method, FAM）
应用最多的是有限差分法，与有限元法相比，差分法占有一定优势：
①方法发展的成熟程度；②实施的难易程度；③应用的广泛性。
3-2 区域离散化方法
定义：用一系列与坐标轴平行的曲线将计算区域划分成很多互不重叠的子区域，并选定每个子区域中的节点的过程。
每一个节点可视为相应微小容积（称为控制容积）的代表，控制容积的边界称为界面（常用虚线表示）。沿坐标轴方向联结相邻两节点而形成的曲线簇称为网格线（常用实线表示）。
一、分类根据节点在子区域中的位置进行
划分
1、外节点法
节点位于子区域的顶点
用直线簇划分出子区域。
节点：子区域的顶点即直线簇的突点
网线线：直线簇
界面：网格间距（相邻两节点间的网格线
P w
Δx B=0
段）的中垂线；
控制容积：不是子区域
特点：子区域与控制容积不重合，先确定节点，后确定界面，practiceA ，方法A
2、内节点法
节点不在控制容积中心，节点位于子区域的中心 界面：划分子区域的直线簇 控制容积：=子区域 节点：控制容积的中心 网格线：
特点：先定界面，后d 定节点,practiceB 方法B 。
二、网格系统的标记方法 以二维为例 1 i-j-n 系统
、j i-1,j
t
n i,j
n+1 i,j
2 P 、E 、W 、N 、S 系统
t 0 p
t
p
节点间距（网格间距）y x δδ,，界面间距y x ??,
对于均分网格：y y x x ?=?=δδ,，内、外节点法的节点分布在区域内部趋于一致。
例：一维网格系统
Δx
2
1+
i
i+1
E
e
P
i ω
W
i-1
2
1-i
(δx ）e +
(δx)w
三、两种区域离散化方法的比较
1、边界节点所代表的控制容积不相同
A ：0≠?
B x ，能更好地考虑边界节点之间的传热作用，更符合实际；
B ：0=?B x ，不能考虑边界节点之间的传热作用，（适合于今后处理2、3类B.C.的附加源项法）
2、网格不均匀时 ①界面位置
A ：界面位于相邻两节点正中间：
B 、否
对导热量计算的准确度不相同e
p E x t t x t )(δ-≈?? P E e ②节点位置
A ：不在控制容积中心
B ：在控制容积中心
３、对求解区域内材料物性突变的适应能力不同 A:
B: 易方便地将物性发生阶跃变化的交界面作为控制容积的界面，从而使同一控制容积内的物性保持均匀一致。
4、对B 、C 突变的适应能力不同
A ：不便于适应
B ：便于适应
结论：两种方法都得到应用，但当有材料物性突变时，建议采用内节点法。
1B.C 3B.C
四、非线笥问题迭代式解法的收敛性 每一层次上满足迭代法求解的收敛条件+相邻次间代数方程的系数变化不太大（亦即未知量 的变化不太大J多数情形下非线性问题迭代式解法是可以收敛的）。 使相邻两层次间未知量变化不太大的措施： 1欠松弛迭代常用逐次欠弛线迭法（SLUR）：一组临时系数下逐线迭代求解+对所得的解 施以欠松弛，再用欠松弛后的解去计算新的系数，常数，以进入下一层次的迭代。 实施：常把欠松弛处理纳入迭代过程，而不是在一个层次迭代完成后再行欠松弛。 .（n 1）川）.'a n bt n b t p =t p （t p ） a p （先）t p n1） = 7a n b t n b b （1一?）屯t p n） co o a'p t p n 9 、a n bt n b b' a'p -a^ ■, b' = b （^ ）（a p ）t p n），用交替方向线迭代法求解这一方程，就实现了SLUR 的迭代求解。为一般化起见，上式中t n b上没有标以迭代层次的符号（J, GS时不相同）。 2、采用拟非稳态法 前面已指出，稳态问题的迭代解法与非稳态问题的步进法十分相似。对于非线性稳态问题， 从代数方程的一组临时系数进入到另一组临时系数亦好象非稳态问题前进了一个时间层，非稳态问题的物理特性：系数热惯性越大（a； = PM v/也I ），温度变化越慢，仿此，对稳态非线性 问题，可在离散方程中加入拟非稳态项，以减小未知量托两个层次间的变化，即 由 （=a n b -S p：V）t p n。= Ua n bt n b b=（3a n b - S p：V a；）t p n。=二a n bt n b b a p tf Za n bt n b - b - a；t p n） （n 1） t p o Ea n b -S p心V +a p 一直进行到t p,t n b收敛，虚拟时间步的大小通过计算实践确定。 3、采用Jacobi点迭代法 中止迭代的判据（该层次迭代）除前述变化率判据外，还可以规定迭代的轮数，例如规定进 行4-6次ADI线迭代就结束该层次上的计算。此时，用收敛速度低的丁迭代也就起到了欠松弛的作用。 五、迭代法的收敛速度 1收敛速度 对给定的代数方程组（包括是临时系数的情形），采用不同的迭代方法求解时，使一定的初始误差缩小成：?倍所需要的迭代轮数K是不相的
传热学 二维稳态导热问题的数值解法 杨达文2011151419 赵树明2011151427 杨文晓2011151421 吴鸿毅2011151416
第一题： a=linspace(0,0.6,121); t1=[60+20*sin(pi*a/0.6)]; t2=repmat(60,[80 121]); s=[t1;t2]; %构造矩阵 for k=1:10000000 %理论最大迭代次数，想多大就设置多大S=s; for j=2:120 for i=2:80 S(i,j)=0.25*(S(i-1,j)+S(i+1,j)+S(i,j-1)+S(i,j+1)); end end if norm(S-s)<0.0001 break; %如果符合精度要求，提前结束迭代else s=S; end end S %输出数值解 数值解数据量太大，这里就不打印出来，只画出温度分布。 画出温度分布： figure(1) xx=linspace(0,0.6,121); yy=linspace(0.4,0,81); [x,y]=meshgrid(xx,yy); surf(x,y,S) axis([0 0.6 0 0.4 60 80]) grid on xlabel('L1') ylabel('L2') zlabel('t（温度）')
.60.66666777778L 1 L 2t （温度）
A0=[S(:,61)]; for k=1:81 B1(k)=A0(81-k+1); end B1 %x=L1/2时y方向的温度 A1=[S(41,:)] %y=L2/2时x方向的温度 x=0:0.005:0.6; y=0:0.005:0.4; A2=60+20*sin(pi*x/0.6)*((exp(pi*0.2/0.6)-exp(-pi*0.2/0.6))/2)/((exp(pi*0.4/0.6)-exp(-pi*0.4/0.6) )/2) %计算y=L2/2时x方向的解析温度 B2=60+20*sin(pi*0.3/0.6)*((exp(pi*y/0.6)-exp(-pi*y/0.6))/2)/((exp(pi*0.4/0.6)-exp(-pi*0.4/0.6))/ 2) %计算x=L1/2时y方向的解析温度 figure(2) subplot(2,2,1); plot(x,A1,'g-.',x,A2,'k:x'); %画出x=L1/2时y方向的温度场、画出x=L1/2时y方向的解析温度场曲线 xlabel('L1');ylabel('t温度'); title('y=L2/2'); legend('数值解','解析解'); subplot(2,2,2); plot(x,A1-A2); %画出具体温度场与解析温度场的差值曲线 xlabel('L1');ylabel('差值'); title('y=L2/2时，比较=数值解-解析解'); subplot(2,2,3); plot(y,B1,'g-.',y,B2,'k:x'); %画出y=L2/2时x方向的温度场、画出y=L2/2时x方向的解析温度场曲线 xlabel('L2');ylabel('t温度'); title('x=L1/2'); legend('数值解','解析解'); subplot(2,2,4); plot(y,B1-B2); %画出具体温度场与解析温度场的差值曲线 xlabel('L2');ylabel('差值'); title('x=L1/2时，比较=数值解-解析解'); y=L2/2时x方向的温度： 60 60.1635347276130 60.3269574318083 60.4901561107239 60.6530189159961 60.8154342294146 60.9772907394204 61.1384775173935 61.2988840936779 61.4584005332920 61.6169175112734 61.7743263876045 61.9305192816696 62.0853891461909 62.2388298405943 62.3907362037523 62.5410041260577 62.6895306207746 62.8362138946214 62.9809534175351 63.1236499915702 63.2642058188844 63.4025245687647 63.5385114436490 63.6720732440951 63.8031184326565 63.9315571966177 64.0573015095482 64.1802651916318 64.3003639687311 64.4175155301449 64.5316395850212 64.6426579173846 64.7504944397430 64.8550752452343 64.9563286582797 65.0541852837075
1、已知：一块厚度为0.1mm 的无限大平板，具有均匀内热源，q ＝50×103W/m 3，，导热系数K ＝10W/m.℃，一侧边界给定温度为75℃，另一侧对流换热，T f ＝25℃，，h=50W/m 2.℃,求解稳态分布。（边界条件用差分代替微分和能量平衡法），画图。(内,外节点) 2、试以下述一维非稳态导热问题为模型，编写求解一维非稳态扩散型问题的通用程序： 00 00000()()()() L L f x x x x L fL L x x x x T T k s c x x T k h T T W x T k h T T W x T T x τρτ =====???+=????=-+??-=-+?= 其中，x 是空间坐标变量，τ是时间坐标变量，T 是温度（分布），k 是材料的导热系数，s 是内热源强度，ρ是材料的密度，c 是材料的比热，h 0和h L 分别是x 0和x L 处流体与固体壁面间的换热系数，而T f0和T fL 分别是固体壁两侧流体的温度，W 0和W L 是x 0和x L 处(非对流换热)热流密度，T 0(x )是固体壁内初始温度分布。注意k 、ρ、c 、s 、h 0 、h L 、W 0和W L 均可以是温度T 和/或空间坐标x 的函数。 具体要求： 1) 将数学模型无量纲化； 2) 考虑各种可能的边界条件和初始条件组合 3) 提供完整的程序设计说明，包括数学推导过程和程序使用说明 3、对于有源项的一维稳态方程， s dx d T dx d u dx d +=)()(φφρ 已知 x=0,φ=0,x=1, φ=1.源项S=0.5-X 利用迎风格式、混合格式、乘方格式求解φ的分布.
数值计算大作业 一、用数值方法求解尺度为100mm×100mm 的二维矩形物体的稳态导热问题。物体的导热系数λ为1.0w/m·K。边界条件分别为： 1、上壁恒热流q=1000w/m2; 2、下壁温度t1=100℃； 3、右侧壁温度t2=0℃; 4、左侧壁与流体对流换热，流体温度tf=0℃，表面传热系数 h 分别为1w/m2·K、10 w/m2·K、100w/m2·K 和1000 w/m2·K; 要求： 1、写出问题的数学描述； 2、写出内部节点和边界节点的差分方程； 3、给出求解方法； 4、编写计算程序(自选程序语言)； 5、画出4个工况下的温度分布图及左、右、下三个边界的热流密度分布图； 6、就一个工况下（自选）对不同网格数下的计算结果进行讨论； 7、就一个工况下（自选）分别采用高斯迭代、高斯——赛德尔迭代及松弛法（亚松弛和超松弛）求解的收敛性（cpu 时间，迭代次数）进行讨论； 8、对4个不同表面传热系数的计算结果进行分析和讨论。 9、自选一种商业软件（fluent 、ansys 等）对问题进行分析，并与自己编程计算结果进行比较验证（一个工况）。（自选项） 1、写出问题的数学描述 设H=0.1m 微分方程 22220t t x y ??+=?? x=0，0 y=H ，0 物理模拟设备的发展综述 摘要：物理模拟技术，作为材料成形工艺的简单实验，可以对复杂成形技术提供可靠的支持，在材料的加工领域里面有不可取代的作用。早期使用橡皮泥，铅块，石蜡等塑性较好的材料来进行复杂成形过程的模拟，以提供合理的设计参数，这种方法浪费大，时间长，效率较低，随着计算机技术的发展，目前更多的模拟同在在电脑上进行，先在热物理模拟机上进行的简单的模拟，得到材料的性能参数，然后在电脑上利用专门的商业软件进行模拟，这样不仅花费小，开发周期短，而且可以使材料的数据得到最大的用途。因此，热物理模拟设备的发展对物理模拟的进步有着举足轻重的作用。 关键词：物理模拟，热物理模拟机，Gleeble 前言 “物理模拟”是一个内涵十分丰富的广义概念，也是一种重要的科学方法和工程手段。通常，“物理模拟”是指缩小或放大比例，或简化条件，或待用材料，用实验的模型来代替原型的研究。对材料和热加工工艺来说，物理模拟通常指利用小试样，借助某种实验装置在线材料制备或热加工过程中受热火受力的物理过程，充分而准确的揭示材料或工件在制备和热加工过程中的组织和性能变化规律，用这些来评定或预测材料制备或加工过程中可能出现的问题，为制定合理的加工工艺和参数，以及研制新材料提供理论指导和技术支持。物理实验可以分为以下两种，一种是在模拟过程中进行的实验，另一种是模拟完成后进行的实验。 以往我们在进行科学研究或者工件的生产过程，为评价工艺方案对材料性能或产品质量的影响，多采用实验的方法，这种简单直接的实验不仅仅要消耗大量的时间，材料和金钱，而且得到结果仅仅能够表示在该工艺下的结果，并不能对其他工艺有太多的指导意义，因此我们必须在实验工艺和方法上进行有一定的创新和改造。 近些年来，随着计算机技术和工程检测技术的迅速发展，物理模拟，数值模拟以及与模拟相关的专业软件都有了长足的进步，相关软件在材料科学和工程领域的运用都取得了非常好的效果，材料学科的研究开始从“经验”走向“科学”。新模拟技术的应用使得人们不仅可以对变形过程有了更加直观的认识，对模具的设计参数好坏有了更加直观的评价，为工艺的制定和工艺参数的设计提供了更加可靠的依据，从而大大减少了新产品和新材料的开发周期和开发费用，降低了企业的成本，提高企业的竞争力。 取步长δx=0.02。已知x=0,Φ=0;x=1,Φ=1.令k=ρｕ/Γ计算结果图表： 程序及数据结果： 追赶法： #include a[i]=2+0.02*k; b[i]=4; c[i]=2-0.02*k; f[i]=0; } tdma(a,b,c,f,x); for(i=0;i 第四章 复习题 1、 试简要说明对导热问题进行有限差分数值计算的基本思想与步骤。 2、 试说明用热平衡法建立节点温度离散方程的基本思想。 3、 推导导热微分方程的步骤和过程与用热平衡法建立节点温度离散方程的过程十分相似， 为什么前者得到的是精确描述，而后者解出的确实近似解。 4、 第三类边界条件边界节点的离散那方程，也可用将第三类边界条件表达式中的一阶导数 用差分公式表示来建立。试比较这样建立起来的离散方程与用热平衡建立起来的离散方程的异同与优劣。 5．对绝热边界条件的数值处理本章采用了哪些方法？试分析比较之． 6．什么是非稳态导热问题的显示格式？什么是显示格式计算中的稳定性问题？ 7．用高斯－塞德尔迭代法求解代数方程时是否一定可以得到收敛德解？不能得出收敛的解时是否因为初场的假设不合适而造成？ 8．有人对一阶导数()()()2 21,253x t t t x t i n i n i n i n ?-+-≈ ??++ 你能否判断这一表达式是否正确，为什么？ 一般性数值计算 4-1、采用计算机进行数值计算不仅是求解偏微分方程的有力工具，而且对一些复杂的经验公式及用无穷级数表示的分析解，也常用计算机来获得数值结果。试用数值方法对Bi=0.1,1,10的三种情况计算下列特征方程的根:)6,2,1( =n n μ 3,2,1,tan == n Bi n n μμ 并用计算机查明，当2 .02≥=δτ a Fo 时用式（3-19）表示的级数的第一项代替整个级数（计 算中用前六项之和来替代）可能引起的误差。 解：Bi n n =μμtan ，不同Bi 下前六个根如下表所示： Bi μ 1 μ2 μ3 μ 4 μ 5 μ 6 0.1 0.3111 3.1731 6.2991 9.4354 12.5743 15.7143 1.0 0.8603 3.4256 6.4373 9.5293 12.6453 15.7713 10 1.4289 4.3058 7.2281 10.2003 13.2142 16.2594 Fo=0.2及0.24时计算结果的对比列于下表： Fo=0.2 δ=x Bi=0.1 Bi=1 Bi=10 第一项的值 0.94879 0.62945 0.11866 前六和的值 0.95142 0.64339 0.12248 比值 0.99724 0.97833 0.96881 Fo=0.2 0=x Bi=0.1 Bi=1 Bi=10 第一项的值 0.99662 0.96514 0.83889 前六项和的值 0.994 0.95064 0.82925 比值 1.002 1.01525 1.01163 Fo=0.24 δ=x 《热物理过程数值模拟》考试题 一、名词解释 1、差分方程的截差 2、差分格式的稳定性 3、时间项差分的C-N 格式 二、问答题 1、试简要分析为什么要对源项进行线性化处理？ 2、试问区域离散化方法有哪几种？它们的异同如何？ 3、何为双向坐标？在什么情况下可转换为单向坐标？ 三、理论推导 1、如图1所示的一维导热问题，给定 边界 热条件，请采用控制容积积分法确定边界节点B 的温度的离散方程。 2、如图2所示的二维导热问题，采用离散化方法B 求第二类边界条件下的边界节点w 和中心节点p 的温度方程。 四、试对如下源项线性化处理方式进行评价。 （1）S=3+7t 7 ,3==p c S S 0 ,73* =+=p p c S t S ① ② ③ 2 ,93* -=+=p p c S t S f t ,α x 图1 图2 A B w q B ·N P E S ω Δy ωωωδδδ)()(,0x x x ==+ - 1、 研究传热问题基本都可用二元二阶线性偏微分方程表达出来，如： ，这里， ，a 、b 、c 、d 、e 、f 、是x 、y 的函数。试描述该控制方程求解域内任一点控制方程可能存在的类型，并简要分析各控制方程类型之差异性。(15分) 2、 试证明控制微分方程 扩散项采用中心差分格式进行离散能保持其迁移特 性之本质（可考虑成常物性、均匀网格）。 （15分） 3、 对于一个双层平板（P 平板与E 平板）的一维导热问题， 如图1所示。P 平板为λP 材料组成，E 平板为λE 材料组成，试简要分析确定合理的P 、E 板界面导热系数的求取方法及表达式。（15分） ① （2）S=5-4t 4 ,5-==p c S S 0,45*=-=p p c S t S 11 ,75*-=+=p p c S t S ② ③ 354t S -=0 ,543 *=-=p p c S t S 2 *5,4p p c t S S -==（3） ① ② ③ 2 *3*25,204p p p c t S t S -=+=),(y x g f e d c b a y x yy xy xx =+++++φφφφφφ,/,/,/22y x x y xx x ??=??=??=φφφφφφy x y xy yy ???=??=/,/222φφφφx e x )( 图1 P/E 双层平板示意 传热学大作业报告二维稳态计算 院系：能源与环境学院 专业：核工程与核技术 姓名：杨予琪 学号：03311507 一、原始题目及要求 计算要求： 1. 写出各未知温度节点的代数方程 2. 分别给出G-S 迭代和Jacobi 迭代程序 3. 程序中给出两种自动判定收敛的方法 4. 考察三种不同初值时的收敛快慢 5. 上下边界的热流量（λ=1W/(m ℃）） 6. 绘出最终结果的等值线 报告要求： 1. 原始题目及要求 2. 各节点的离散化的代数方程 3. 源程序 4. 不同初值时的收敛快慢 5. 上下边界的热流量（λ=1W/(m ℃）） 6. 计算结果的等温线图 7. 计算小结 二、各节点的离散化的代数方程 左上角节点 )(21 1,22,11,1t t t += 右上角节点 )(2 15,24,15,1t t t += 左下角节点 C t ?=1001,5 右下角节点 )2(211,24,55,5λ λ x h t t x h t ?++?+= 左边界节点 C t i ?=1001,，42≤≤i 上边界节点 C t j ?=200,1，42≤≤j 右边界节点 )2(415,15,14,5,+-++= i i i i t t t t ，42≤≤i 下边界节点 )42()2(211,51,5,4,5∞+-?+++?+=t x h t t t x h t j j j j λλ ，42≤≤j 内部节点 )(2 1,1,11,1,,j i j i j i j i j i t t t t t +-+-+++= ，4,2≤≤j i 三、源程序 1、G-S 迭代法 t=zeros(5,5); t0=zeros(5,5); dteps=0.0001; for i=2:5 %左边界节点 t(i,1)=100; end for j=2:4 %上边界节点 t(1,j)=200; end t(1,1)=(t(1,2)+t(2,1))/2; t for k=1:100 for i=2:4 %内部节点 for j=2:4 t(i,j)=(t(i-1,j)+t(i+1,j)+t(i,j-1)+t(i,j+1))/4; end end t(1,5)=(t(1,4)+t(2,5))/2;%右上角节点 for i=2:4;%右边界节点 t(i,5)=(2*t(i,4)+t(i-1,5)+t(i+1,5))/4; end for j=2:4; %下边界节点 1、Jacobi 迭代 在Jacobi 迭代法中任一点上未知值的更新是用上一轮迭代中所获得的各邻 点之值来计算的，即 kk k k l l n l k n k a b T a T /)(1)1()(+=∑≠=- k=1,2,...,L 1×M 1 这里带括号的上角标表示迭代轮数。所谓一轮是指把求解区域中每一节点之值都更新一次的运算环节。显然，采用Jacobi 迭代式，迭代前进的方向（又称扫描方向）并不影响迭代收敛速度。这种迭代法收敛速度很慢，一般较少采用。但对强烈的非线性问题，如果两个层次的迭代之间未知量的变化过大，容易引起非线性问题迭代的发散。在规定每一层次计算的迭代轮次数的情况下，有利于Jacobi 迭代有利于非线性问题迭代的收敛。 2、Gauss-Seidel 迭代 在这种迭代法中，每一种计算总是取邻点的最新值来进行。如果每一轮迭代按T 的下角标由小到大的方式进行，则可表示为： kk k M L k l n l kl k l l n l kl n k a b T a T a T /)(1 11 ) 1(1 1) ()(++ =∑∑?+=--≠= 此时迭代计算进行的方向（即扫描方向）会影响到收敛速度，这是与边界条件的影响传入到区域内部的快慢有关的。 3、例题： 一矩形薄板几何尺寸如图所示，薄板左侧的边界温度T L =100K ，右侧温度T R =300K ，上侧温度T T =200K ，下侧温度T B =200K ，其余各面绝热，求板上个节点的温度。要求节点数目可以变化，写出程序。 解析： ⑴列出描述问题的微分方程和定解条件。 22 220t t x y ??+=??；对于离散化的问题，其微分方程根据热平衡原理得到： 粘土的热物理参数和冻结过程中的温度场演变研究冻土是一种由固、液、气、冰四相物质组成的复合材料。冻土中的冰相对其物理、力学和工程性质具有决定影响。 因此,研究冰相的存在形式和形成过程,是揭示冻土复杂性质的基础和条件。而认识冻土的热参数及其导热特性,则是揭示冻土形成规律必不可少的条件。 本论文从冻土的材料组成入手,研究了冻结过程中的热量构成,建立了冻土 热物理参数的计算新模型,研发了室内模型试验装置,开展室内模型试验,通过模型试验结合数值模拟分析,验证了本文提出的冻土热物理参数计算新模型的有效性。为寒冷地区的工程建设和冻结法施工,提供了设计理论基础。 获得的成果主要包括以下几个方面。(1)揭示了土在冻结过程中的热量构成。 将未冻土当作颗粒、孔隙水和孔隙气组成的复合体,则冻结过程实质上是孔隙水的降温过程、孔隙水的结冰过程、孔隙冰的降温过程、土颗粒的降温过程共四个物理过程的组合。因此,土在冻结过程中的热量消耗主要用于颗粒的温度变化、水(冰)的温度变化和冰水相变。 (2)揭示了土在冻结过程中,比热和导热系数的演变特点。由于冰和水的含量在冻结阶段是不断变化的,而冰的比热和导热系数完全不同于水的比热和导热系数,因此土体的比热和导热系数是一个动态变化过程。 因此,研究土冻结阶段的导热系数和比热,必须与孔隙水的相变过程和相应 潜热相联系。(3)根据是否存在冰水相变或潜热,提出了一种冻结过程的划分方法。 随温度降低,若土中水一直以液态形式存在则称之为未冻阶段;若存在固态 冰且液态水随温度降低逐渐减少则称之为冻结阶段;若存在固态冰但液态水不再随温度降低而减少则称之为冻实阶段。(4)建立了完整的冻结过程三阶段比热计 二维导热物体温度场的数值模拟 一、物理问题 有一个用砖砌成的长方形截面的冷空气通道， 于纸面方向上用冷空气及砖墙的温度变化很小， 可以近似地予以忽略。 在下列两种情况下试计算： 砖墙横截面上的温度分布；垂直于纸面方向的每 米长度上通过砖墙的导热量。 第一种情况：内外壁分别均匀维持在 0℃及 30℃； 第二种情况：内外壁均为第三类边界条 件， 且已知： t 1 30 C,h 1 10.35W / m 2 K 2 t 2 10 C, h 2 3.93W / m 2 K 砖墙导热系数 0.35/ m K 二、数学描写 由对称的界面必是绝热面， 态、无内热源的导热问题。 控制方程： 22 tt 22 xy 边界条件： 第一种情况： 由对称性知边界 1 绝热： 边界 2 为等温边界，满足第一类边界条件： t w 0 C ； 边界 3 为等温边界，满足第一类边界条件： t w 30 C 。 第一种情况： 由对称性知边界 1 绝热： q w 0； 边界 2 为对流边界，满足第三类边界条件： q w ( t )w h 2(t w 可取左上方的四分之一墙角为研究对象， 该问题为二维、 稳 图1- t f )； n t 边界3 为对流边界，满足第三类边界条件：q w ( ) w h 2 (t w t f )。 w n w 2 w f 0,m 6,n 1~ 7;m 7 ~ 16,n 7 30,m 1,n 1~12;m 2 ~ 16,n 12 三、方程离散 用一系列与坐标轴平行的间隔 0.1m 的二维网格线 将温度区域划分为若干子区域，如图 1-3 所示。 采用热平衡法， 利用傅里叶导热定律和能量守恒定 律，按照以导入元体（ m,n ）方向的热流量为正，列写 每个节点代表的元体的代数方程， 第一种情况： 边界点： 1 边界 绝热边界） ： 边界 图1-3 t m ,1 t 16,n 等温内边界） ： 14 (2t m,2 1 4 (2t 15,n t m 1,1 t m 1,1)，m 2 ~ 5 t 16,n 1 t 16,n 1)， n 8 ~ 11 边界 等温外边界） ： 内节 点： 1 (t t t t ) 4 m 1,n m 1,n m ,n 1 m,n 1 m 2 ~ 5,n 2 ~11;m 6 ~ 15,n 8 ~ 11 t m,n 第二种情况 边界点： 边界 1（绝热边界） ： t m ,1 1 4 (2t m,2 t m 1,1 t m 1,1)，m 2 ~ 5 t 16,n 1 4 (2t 15,n t 16,n 1 t 16,n 1)， n 8 ~11 4 边界 2（内对流边界） ： t6,n 2t 5,n t 6,n 1 t 6,n 1 2Bi 1t 1 ，n 1~ 6 6,n 2(Bi 2) t m,n t m,n 1. 传热学的发展概述 18世纪30年代首先从英国开始的工业革命促进了生产力的空前发展。生产力的发展为自然科学的发展成长开辟了广阔的道路。传热学这一门学科就是在这种大背景下发展成长起来的。导热和对流两种基本热量传递方式早为人们所认识，第三种热量传递方式则是在1803年发现了红外线才确认的，它就是热辐射方式。在批判“热素说”确认热是一种运动的过程中，科学史上的两个著名实验起着关键作用。其一是1798年伦福特(B ．T ．Rumford)钻炮筒大量发热的实验，其二是 1799年戴维(H ．Davy)两块冰块摩擦生热化为水的实验。确认热来源于物体本身内部的运动开辟了探求导热规律的途径。1804年毕渥根据实验提出了一个公式，认为每单位时间通过每单位面积的导热热量正比例于两侧表面温差，反比例于壁厚，比例系数是材料的物理性质。傅里叶于1822年发表了他的著名论著“热的解析理论”，成功地完成了创建导热理论的任务。他提出的导热定律正确概括了导热实验的结果，现称为傅里叶定律，奠定了导热理论的基础。他从傅里叶定律和能量守恒定律推出的导热微分方程是导热问题正确的数学描写，成为求解大多数工程导热问题的出发点。他所提出的采用无穷级数表示理论解的方法开辟了数学求解的新途径。傅里叶被公认为导热理论的奠基人。在傅里叶之后，导热理论求解的领域不断扩大。同样，自1823年M. Navier 提出流动方程以来，通过1845 年 G.G. Stokes 的改进，完成了流体流动基本方程的创建任务。流体流动理论是更加复杂的对流换热理论的必要前提，1909和1915年W. Nusselt 开辟了在无量纲数原则关系正确指导下，通过实验研究对流换热问题的一种基本方法。1904 年，L. Prandtl 提出的对流边界层理论使流动微分方程得到了简化，1921年 E. Pohlhausen 基于流动边界层理论引进了热边界层的概念，为对流传热微分方程的理论求解建立了基础。在辐射传热研究方面，19世纪J. Stefan 根据实验确定了黑体辐射力正比于它的绝对温度的四次方的规律，1900年M.Planck 提出的量子假说奠定了热辐射传热理论基础。上述传热理论为传热分析解析、数值以及实验研究奠定了理论基础。还要特别提到的是，由于计算机的迅速发展，用数值方法对传热问题的分析研究取得了重大进展，在20世纪70年代已经形成一个新兴分支—数值传热学。近年来，数值传热学得到了蓬勃的发展[2-4]。 2. 传热分析计算理论 热量传递主要有三种传递形式，分别是热传导、热对流和热辐射。热传导是指两个相互接触良好的物体之间的能量交换或一个物体由于其自身温度梯度而 引起的内部能量的传递。其遵循傅里叶定律[5]：dT q dx λ=-，其中λ是热导率， dT dx 是温度梯度，q 是热流密度。热对流是指在物体与其周围介质之间发生的热量交换。热对流分为自然对流和强制对流，用牛顿冷却方程描述为()w f q h t t =-，其中h 为表面传热系数，w t 为物体表面的温度，f t 为物体周围流体的温度。一个 物体或两个物体之间通过电磁波形式进行的能量传递交换称为热辐射，通常由斯 热物理过程的数值模拟Numerical Simulation of Thermophysics Process 讲稿 主讲：李隆键 第一章概论 1．1流动与传热过程的予测方法及特点 流动、传热、燃烧问题是热工类各专业和机械类动力机械专业所研究和解决的主要问题之一，燃烧问题实际上是有化学反应的流动与传热问题，推而广之，在所有热物理过程中，几乎都涉及到流动、传热问题。 预测的重要性： ①在规定设计参数的相应的结构下，热物理过程是否满足要求，达到预定的指 标？要预测； ②优化设计，不同方案的比较，要预测； ③减少设计、生产、再设计和再生产的费用； ④减少设计更改； ⑤减少试验和测量次数。 问题的核心：速度场、温度场（传热量）、浓度场等。 一、热物理问题的予测方法：理论分析法、实验测定、数值模拟 1、理论分析 以数学分析为基础，求解描述热物理过程的定解问题，获得函数形式的解，表示求解区域内物理量连续分布的场（速度场、温度场、浓度场……）。 控制方程+单值条件（数学模型）→理论解（分析解，解析解） 根据解的准确程度，又可再分为： （1）精确分析解（严格解） 特点：函数形式的解；它在求解区域精确地满足定解问题。 具体解法：直接积分法、分离变量法、积分变换法、热源法、映射法。 （2）近似分析解法 特点：函数形式的解，在求解区域上近似地满足定解问题（但在总量上满足相应的守恒原理，动量守恒、动量守恒、能量守恒、质量守恒）。 具体解法：积分法（从积分方程出发） 变分近似解法 摄动法（从微分方程出发） 2、实验测定 （1）纯实验法 （2）相似理论实验法：同类相似，减少变量数目→减少工作量，得到规律性结 、室内一根水平放置的无限长的蒸汽管道，其保温层外径d=583 mm，外表面 实测平均温度及空气温度分别为，此时空气与管道外 表面间的自然对流换热的表面传热系数h=3.42 W /(m2 K),墙壁的温度近似取为 室内空气的温度，保温层外表面的发射率 问：(1)此管道外壁的换热必须考虑哪些热量传递方式； (2)计算每米长度管道外壁的总散热量。(12分) 解： (1)此管道外壁的换热有辐射换热和自然对流换热两种方式。 (2)把管道每米长度上的散热量记为qi 当仅考虑自然对流时，单位长度上的自然对流散热 q i,c =二d h t =二dh (j - t f ) = 3.14 0.583 3.42 (48 - 23 ) 二156 .5(W / m) 近似地取墙壁的温度为室内空气温度，于是每米长度管道外表面与室内物体及墙壁 之间的辐射为： q i厂d (T； -T；) = 3.14 0.583 5.67 10》0.9 [(48 273)4-(23 273)4] = 274.7(W /m) 总的散热量为q i = q i,c +q i,r = 156.5 +274.7 = 431.2(W/m) 2、如图所示的墙壁，其导热系数为50W/(m- K),厚度为50mm在稳态情况下的 墙壁内的一维温度分布为：t=200-2000x 2，式中t的单位为°C, x单位为m 试 求： t (1) 墙壁两侧表面的热流密度； (2) 墙壁内单位体积的内热源生成的热量 2 t =200 —2000x 解：(1)由傅立叶定律: ① dt W q ' (―4000x) = 4000二x A dx 所以墙壁两侧的热流密度： q x _. =4000 50 0.05 =10000 (1)由导热微分方程 茫?生=0得: dx 扎 3、一根直径为1mm 勺铜导线，每米的电阻为2.22 10 。导线外包有厚度为 0.5mm 导热系数为0.15W/(m ? K)的绝缘层。限定绝缘层的最高温度为 65°C,绝 缘层的外表面温度受环境影响，假设为40°C 。试确定该导线的最大允许电流为多 少？ 解：(1)以长度为L 的导线为例，导线通电后生成的热量为I 2RL ,其中的一部分 热量用于导线的升温，其热量为心务中：一部分热量通过绝热层的 导热传到大气中，其热量为：门二 1 , d In 2 L d 1 根据能量守恒定律知：l 2RL -门 述二厶E = I 2RL -门 即 E = — L dT m = I 2RL - t w1 _tw2 4 di 1 , d 2 In 2 L d 1 q v 、d 2t ——' 2 dx =-(7000)= 4000 50 二 200000 W/m 3 t w1 - t w2 。 2 q x 卫=4000.： 0 = 0 2．迁移性 传递过程的两种机制：扩散传递、对流传递 两种机制在物理特上的差异：对信息或扰动的传递性质上有很大的区别 扩散传递：物质分子不规则热运动所致，这种分子的不规则热运动对空间不同方向的几率是一样，所以扩扩散作用可以把发生在某一位置处的扰动影响向各个方向传递。 对流传递：是流体微团的宏观定向运动，带有强烈的方向性。对流作用只能将发生在某一位置处的扰动向其下游方向传递，而不会逆向传播。 图示 ε扩散 对流 ε 扩散与对流作用在物理本质上的这种差异，应在其各自的差分格式中反映出来。 （1）扩散项的中心差分把扰动向四周均匀传递 一堆非稳态扩散方程： )()(x x ??Γ??=??φ τρφ 对于常物性222x ??Γ=?φτφρ 差分格式：时间导数向前差分，空间导数中心差分（显式），均匀网格x x δ=? 2 1 11) (2x n i n i n i n i n i ?+-Γ=?--++φφφτφφρ 为简化起见，假定初始时刻物理量场已均匀化，且0=φ，在某一时刻（例如第n 时层），节点i 处突然有一个扰动ε，而其余各节点的扰动均匀为零，如图所示，随着时间的推移，这一扰 动传递的情形可由上述差分方程来确定，（n+1）时层： 2 1 11) (2x n i n i n i n i n i ?+-Γ=?--++φφφτφφρ 其中011==-+n i n i φφ ∴ )21())(21(2 21x x n i n i ?Γ?-=Γ??- =+ρτ ερτφφ 在这里，网格傅里叶数)/(2 x F ?Γ?=?ρτ，按稳定性要求， 1210,2/12 2≤?Γ?-≤∴≤?Γ?x x ρτ ρτ， 《传热学》上机大作业 二维导热物体温度场的数值模拟 学校：西安交通大学 姓名：张晓璐 学号:10031133 班级：能动A06 一．问题（4-23） 有一个用砖砌成的长方形截面的冷空气通道，形状和截面尺寸如下图所示，假设在垂直纸面方向冷空气和砖墙的温度变化很小，差别可以近似的予以忽略。在下列两种情况下计算：砖墙横截面上的温度分布；垂直于纸面方向上的每米长度上通过墙砖上的导热量。 第一种情况：内外壁分别维持在10C ?和30C ? 第二种情况：内外壁与流体发生对流传热，且有C t f ?=101， )/(2021k m W h ?=，C t f ?=302，)/(422k m W h ?=，K m W ?=/53.0λ 二．问题分析 1.控制方程 02222=??+??y t x t 2.边界条件 所研究物体关于横轴和纵轴对称，所以只研究四分之一即可，如下图： 对上图所示各边界： 边界1：由对称性可知：此边界绝热，0=w q 。 边界2：情况一：第一类边界条件 C t w ?=10 情况二：第三类边界条件 )()( 11f w w w t t h n t q -=??-=λ 边界3：情况一：第一类边界条件 C t w ?=30 情况二：第三类边界条件 )()( 22f w w w t t h n t q -=??-=λ 三：区域离散化及公式推导 如下图所示，用一系列和坐标抽平行的相互间隔cm 10的网格线将所示区域离散化，每个交点可以看做节点，该节点的温度近似看做节点所在区域的平均温度。利用热平衡法列出各个节点温度的代数方程。 第一种情况： 内部角点： 计算机在材料科学中的应用 课程设计 材料0707班 组长：赵宇 组员：杨林波 阴晓宁 王昆 陈晓辉 周琦 2011/5/27 目录 一、课程设计及团队介绍 (2) 二、课程设计内容 (3) 模型1：蒸汽管道保温层的稳态温度场 (3) 1、建立模型 (3) 2、边界条件 (3) 3、差分方程 (4) 模型2：台式电脑CPU散热片的二维稳态温度场 (4) 1、建立模型 (4) 2、热传导方程 (5) 3、边界条件 (5) 4、初始条件 (6) 5、差分方程 (6) 模型3：无限大钢板淬火冷却过程一维非稳定温度场 (7) 1．建立模型 (7) 2. 热传导方程 (8) 3．边界条件 (8) 4．初始条件 (8) 5．差分方程 (9) 模型4：平板焊接过程的二维温度场 (11) 1、建立模型 (11) 2、边界条件 (12) 3、初始条件 (12) 4、差分方程 (13) 模型5：立方体钢锭三维非稳定温度场 (14) 1、建立模型 (14) 2、边界条件 (14) 3、初始条件 (15) 4、差分方程 (15) 模型6：T型钢淬火 (18) 1、模型建立 (18) 2、边界条件 (19) 3、初始条件 (19) 4、差分方程 (20) 5、C语言源程序。（略） (21) 6、数据处理 (21) 三、课程总结与感想 (25) 1 / 26 一、课程设计及团队介绍 我们小组共有6人，包括组长赵宇，组员杨林波、阴晓宁、王昆、陈晓辉、周琦。经过组员的共同努力，我们顺利的完成了本次作业。我们的作业内容包括，6个实际问题的模型建立与差分方程推导，以及T型钢淬火模型的编程计算与数据处理等。 以下是我们小组各组员的任务完成情况： 以下是我们小组各成员的排序情况：热物理模拟设备的发展
计算传热学中国石油大学(华东)第四章大作业
传热学第四版课后题答案第四章
热物理过程数值模拟考试题(2009冉)
传热学大作业报告 二维稳态导热
计算传热学数值模拟
粘土的热物理参数和冻结过程中的温度场演变研究
西安交通大学传热学大作业二维温度场热电比拟实验1
计算传热学-传热基本原理及其有限元应用
热物理过程的数值模拟-计算传热学1
传热学计算例题
热物理过程的数值模拟-计算传热学2
西安交通大学传热学大作业
传热学数值模拟
相关主题
文本预览