热物理过程的数值模拟-计算传热学3汇总
四、非线笥问题迭代式解法的收敛性
每一层次上满足迭代法求解的收敛条件+相邻次间代数方程的系数变化不太大（亦即未知量的变化不太大←多数情形下非线性问题迭代式解法是可以收敛的）。
使相邻两层次间未知量变化不太大的措施： 1、欠松弛迭代 常用逐次欠弛线迭法（SLUR ）：一组临时系数下逐线迭代求解+对所得的解施以欠松弛，再用欠松弛后的解去计算新的系数，常数，以进入下一层次的迭代。
实施：常把欠松弛处理纳入迭代过程，而不是在一个层次迭代完成后再行欠松弛。
)(
)
()()1(n p p
n n n p n p t a b bt a t t -∑+=+ω )()1()
1()(
n p p
n n n p p
t a b b t b a t a ω
ωω
-+++∑=+
∑+=+')1('b b bt a t a n n n p p
)('))(1(',n p p p p t a b b a a ωωω-+==，用交替方向线迭代法求解这一方程，就实现了SLUR
的迭代求解。为一般化起见，上式中b t n 上没有标以迭代层次的符号（J ，GS 时不相同）。 2、采用拟非稳态法
前面已指出，稳态问题的迭代解法与非稳态问题的步进法十分相似。对于非线性稳态问题，从代数方程的一组临时系数进入到另一组临时系数亦好象非稳态问题前进了一个时间层，非稳态问题的物理特性：系数热惯性越大（↑??=τρ/v c a o
p ），温度变化越慢，仿此，对稳态非线性问题，可在离散方程中加入拟非稳态项，以减小未知量托两个层次间的变化，即
由
)()1()1()()(n p o p n n n p o p p n n n n p p n t a b b bt a t a V S b a b b bt a t V S b a ++∑=+?-∑?+∑=?-∑++
o p
p n n p
o p n n n p a V S b a t a b b bt a t +?-∑++∑=
+)
()1(
一直进行到b t t n p ,收敛，虚拟时间步τ?的大小通过计算实践确定。
3、采用Jacobi 点迭代法
中止迭代的判据（该层次迭代）除前述变化率判据外，还可以规定迭代的轮数，例如规定进行4-6次ADI 线迭代就结束该层次上的计算。此时，用收敛速度低的丁迭代也就起到了欠松弛的作用。
五、迭代法的收敛速度 1、收敛速度
对给定的代数方程组（包括是临时系数的情形），采用不同的迭代方法求解时，使一定的初始误差缩小成α倍所需要的迭代轮数K 是不相的。1<α
)0()0()(e G e G e k k k ≤=
因为G 是对称的，所以)()()(22
G G G G G T ρρρ===
所以 2
)0(2
)0(2
2
)
())(e G e G e k k k ρ=≤
)(,))((/2
)
0(2
)
(G kl l G e e n n k k ραρα≤≤=
即 ))(/()(/G l l G l l k n n n n ραρα--=≥ )(,G ρα均<1 对于给定的α，所需迭代轮数k 与)(G l n ρ-成反比，规定用
)(G l R n ρ-=
表示迭代法的收敛速度，则
k l R or R l k n n //αα-≥-≥
即所需迭代轮数与收敛速度成反比，收敛速度又与谱半径成反比，收敛速度愈快，迭代轮数愈少。
注意：不同的迭代方法每进行一轮迭代所需的运算次数不同，最终所需的计算时间的多少取决于迭代轮数及每一轮迭代所需的时间。
2、收敛性的定性分析
为什么不同的迭代方法的收敛速度不同，亦即为达到满足一定精度要求所需的迭代轮数不同？以二维常物性、无内热源、稳态导热问题来进行讨论。
0'2
22=??+??y t
x t 1B 、C S S N N w w E E p p t a t a t a t a t a +++=
迭代法需要假定一个初场，例如假定一个均场，从微分方程为看，均匀是其一个解，但却不是所研究问题的解，为什么？因为它虽然满足内部节点上的离散方程，却不满足与边界有关的节点的离散方程（图中红点），即不满足边界条件。所以迭代法的实质是要通过迭代，尽快建立起与边界条件相适应的φ变量场，
关键：必须使B 、C 的影响迅速传入计算区域内部，以改进节点φ变量值，尽快与B 、C 相应，B 、C 的影响传入愈快，逼近真解就愈快，收敛就越快！
B 、
C 的影响传入计算区域内部的快慢与哪些因素有关？ （1）与迭代方法有关
J 迭代：节点温度的更新均用上轮迭代所得的“旧”值来计算，所以完成一轮迭代后，B 、C
t
TB
t RB
t BB
t wB
i
j y
x
的影响只能传入与边界相邻的一批节点上，即仅可传入一个网格，且扫描方向与收敛快慢无关。要在以后各轮迭代中，B 、C 的影响才由这些节点逐步向内渗透，所以收敛慢。
GS 迭代：假设从左向右扫描，则每做完一轮迭代，左边界和下边界的影响传遍全区域，而右边界的影响只能传入一个网格，且收敛速度受迭代扫描方向的影响。
线迭代：GS 线迭代。自左→右扫描，完成一轮迭代不仅左边界的影响逐步传入，而且在每一列的直接求解中，上、下边界的影响全部传入到该列的各节点上，即一轮迭代使左、上、下边界影响传入全区域，但右边界影响仍仅传入一个网格。
ADI ：一轮迭代包括一次逐行、一次逐列的扫描；所以在每一轮迭代后所有边界的影响均传入计算区域内部，从而加快了收敛速度。
收敛速度的比较，正方形区域，1B 、C ，Laplace 方程五点格式，均匀网格步长为h 。
迭代方法 点迭代 线迭代 Jacobi h 2/2 h 2 Gouss-Seidel
h 2 2h 2 SOR
2h
22h
（2）与边界条件的性质有关
tw
x
t
tw
x
t
tf tw 定
向点
λ/l
x t
1B 、C 规定了边界节点的温度，影响直接传入计算区域内部； 3B 、C 规定了环境温度及定向点位置，
∞-=??t t x
t
l λ，对边界温度的限定程度比1B 、C 时
弱，所以对内部的影响也较弱；或将f t 视为外部温度，其对计算区域内部的影响被外部换热热阻削弱，而1B 、C 可视为∞→α或外部热阻0→的极限情况，故3B 、C 的影响比1B 、C 时弱。
2B 、C 仅规定了壁面的钭率，壁温完全不确定，对内部节点温度值的改进提供的信息最少，收敛最慢。
可见，为了提高代数方程迭代解法的收敛速度，应力求使边界条件的影响迅速传入计算区域内部，措施：
①增加迭代解法中直接解法的成分，从点迭代→线迭代→ADI ；
②适当选择扫描的始边，多以1B 、C 或3B 、C 的边界为始边，少以2B 、C （尤其是绝热边界）为扫描始边。
5-4 不规则区域的处理—网格生成技术
如何对不规则区域进行有效的处理，以便于进行传热与流动过程的数值模拟，是近年来计算传热研究中的一个重要课题。
以上讨论的传热过程大都发生在规则而简单的区域中，但许多实际的热传递现象是在不规则的区域中进行的，例如：
①套片管中肋片的传热 ②渐扩通道中的流动与换热
③环形空间中的自然对流 ④流体外掠管束
以上四种情形中的流动与换热不是直角坐标、圆柱轴对称坐标或极坐标所能方便地予以描述。 虽然有限元法在处理不规则边界方面显示了极大的优越性，但就流动与传热而言，在计算技巧与方法方面，有限差分法都比有限元法成熟。
用有限差分法处理这类问题的方法可以归纳为以下几种。 1、采用阶梯形边界（网格）
用阶梯形边界近似代替四分之一圆弧边界。阶梯形边界（网格）是采用有限差分法计算不规则区域的最普通的方法。
缺点：程序缺少通用性；曲面边界上网格必须划得比较细密（否
则会引起较大误差）。
2、采用区域扩充法
当计算区域的边界不规则程度不很严重时，可以采用区域扩充法，把计算区域扩充为直角坐标，圆柱坐标等常规正交坐标系中易于描述的形状。条件：保证原计算区域的情况不变！’
q
α,t f
② ①
v=实际值
v=1030
流动：把计算区域扩充到图中虚线所示的整个圆形通道，从而可以应用圆柱轴对称坐标系中的控制方程加以描述。
如何保证原来的情况不变？孔板区视为粘性无限大的“流体”，而其余区域的流体粘性值就等于真实值，边界上流速赋为零，计算中零边值将迅速传到孔板区域内，有效地模拟了孔板的存在。
传热：对三种不同的边界条件，具体的处理方式不同
（1）均匀壁温边界条件
令扩充区域中的导热系数为无限大，而扩充后的区域边界温度则等于已知值。 （2）绝热的边界条件 令扩充区域中的材料导热系数为零即可实现此条件 （3）均匀热流边界条件
可应用附加源项法来实现，真实边界上均匀热流可以附加源项的形式置于与真实边界相邻的控制容积中去，而扩充区域则处于绝热状态。
周期性二维渐扩、渐缩通道中的换热，倾斜的边界上作用有均匀的热流。采用左下所示阶梯形网格，并把计算区域扩充到一个长方形，以便利用直角坐标系求解。 看一个网格单元的放大后的情形，P 控制容积的附加源项为
abcd of ad V qL S ?=/
L ef —实际边界与控制容积P 的两条边界相交部分的长度；
abcd V ?—控制容积P 的体积
扩充区域0=扩λ，则控容P 中的附加源项ad S 不会向扩充区域传递，从而实现了实附边界上的均匀热流加热条件。
（4）外部对流换热边界条件，f t l -
根据附加源项法，此时P 控制容积的两个附加源项为
V
L S V
L t S e ad p of
f
ad c ??+-
=??
+=
λδαλδα//11
//1,,
δ—网格节点P 到实际边界的距离0=扩λ
区域扩充法的优点：可以用按规则区域编制的通用程序来计算非规则区域的问题；易于实施。缺点：浪费一些计算机内存及计算时间。 3、采用三角形网格 外节点法对于不规则区域中的导热问题，采用三角形网格可以得
到比较满意的结果。
从不规则区域的三角形网格中划出围绕节点P 的多边形来分析。
确定节点P 所代表的控制容积：①三角形外心作为控容的顶点，要求：三角形为锐角三角形，以保证外心一定在三角形内；
②以三角形重心作为控容的顶点，对三角形形状无限制；外心为控制
容积顶点，P 控制容积如图所示，各部分平均导热系数分别为f c b a λλλλ ,,,，用控制容积能量平衡法建立离散方程。
q q
实际边界 d
c p b
δ e
a
λ扩≡0 计算 边界 q
4 5
6
1
2
3 e d c
b f
P-1之间的热导1'11'11///p b b p a a p L L L L c λλ+=
21'111'1cot 21
cot 21ββp b p a L L L L ==
所以)cot cot (2
1
cot 21cot 2121211βλβλβλβλb a b a p C +=+=
常物性时，)cot (cot 2
1
211ββλ+=p C
其它各节点（2-6）与P 节点之间的热导亦按上式计算，而只需把21,ββ改换成相应顶角即可，P 控制容积的热容：
变物性时，图示各部阴影部分的热容量之和∑?i
i
i p V
c )(,ρ
由控制容积P 的热平衡可得
∑∑?+-=?-?i
p p i pi o p p i
V S t t C t t V c )()()
(τ
ρ
常物性时，上式变为
∑?+-=?-?i
p p i pi o p p p
V S t t C t t V c )()(τ
ρ
以上二式右端温度值的时层未标出，它取决于所采用的格式。 边界条件的处理：
第一类边界条件，无需专门处理；
第二、三类B 、C 边界上的节点，需要仿上述方法建立补充方程和图第二、三类边界上的节点P 的控制容积
3311cot 21
,cot 21βλβλc p b p C C ==
)cot cot (2
1
422
βλβλc b p C += 边界加入热量则计入到P 控制容积的源项中，对第三类边界条件，传入热量ad p f L t t )(-α，其中p t 的时层取决于所采用的格式。 Winslow 指出，当导热系数为常数时，三角形网格所建立的离散方程的系数矩阵具有正定性，可用SOR 求解；导热系数与温度有关时，为求迭代收敛，不宜采用SOR ，且两次迭代之间的pi C 可用欠松弛。
根据实际问题的需要，可以采用三角形网格与矩形网格的组合结构，例如二维纵向肋片管。
三角形网格的缺点：节点位置的确定、编号，节点间距的计算比较费时程序比较复杂。 内节点法：以每个三角形的外心为节点，其余同上。 4、组合坐标法
2
3
P 1 1,t f
b
c
a
d
β2 β4
二维平直一弯曲段组成的通道，可采用坐标系组合的方法，例如直角坐标与极坐标的组合，把该计算系统分成三个区域，在区域I 、III 中采用直角坐标系，在区域II 中采用极坐标系。
三个区域中垂直于主流方向的坐标增量是一致
的，dr dy =，但在主流方向上的坐标不一致，为求得一致，在I 、III 区域中定义i R x x /~=，则I 、
III 区域中的x ~也是弧度，与II 区中一致，从而可以
从入口截面开始就统一地连续计量x ~坐标。 在各区域中分别按相应的控制方程建立其离散方程，然后在整个区域内统一地求解代数方程。
三个区域的两个分界线作为主控制容积的界面线，所以当采用交错网格时，主流方向速度的
控制容积将跨越两个区域，如图所示。该速控制容积的w E a a ,分别按级坐标和直角坐标系计算，而S N a a ,可按并联处理：先分别按两区域中的半控制容积计算，后相加，源项亦仿此。
5、采用特殊的正交坐标系
三维正交曲线坐标系中，稳态对流一扩散问题的通用控制方程为
)(1)(1)(112133
3213122321321321φρφρφρh h u x h h h h h u x h h h h h u x h h h ??
+??+??
=
S x h h h x h h h x h h h x h h h +??Γ??
+??Γ??)(1)(13
321332122311321φφ
控制方程写成这种形式时，一切由曲线坐标，滞流模型引起的附加项都包括在S 中去了，因
此，S 的表达式一般比较复杂。
曲面边界平面通道中层流充分发展对流换热，设x 3主为主流方向，且为直线坐标，再利用充分发展条件，可得
C h Sh x h h x x h h x =-=??Γ??+??Γ??212
2121121)()(φφ 式中，C ,Γ例于表中
φ
Γ
C
u 3 u 321/dx dp h h t
u/p r
3321/)(x t u h h ??ρ
y
x
Ri r
θ
P u E
速度u 的控制容积
在该流道截面的二维正交曲线坐标系中取出一个控制容积P 来分析，x 3方向控制容积厚度为
3x ?，将上面的控制方程在该控制容积上作积分，等式左端第一项为
?
???Γ???=n
s
e
dx dx x h h x x I 2111213)(φω
=21
121123
])()[(dx x h h x h h x n
s
e ?
??Γ-??Γ
?ωφφ 利用度规系数的性质，上式进一步变为
])()([2212
132?
?-Γ--Γ??n
s
n
s p dx e
PE
p E e
dx h wp s h s x I ω
ωω
δφφδφφ =])()()()()()
([21123SP
e
w p PE e P E e s s s s x δφφδφφω?-Γ-?-Γ?
类似地，左端第二项积分为：
])()()()()()
([21213SP
s
S p s NP n p N n s s s s x II δφφδφφ?-Γ-?-Γ?≈
等号右端的积分为：21)(h h S S C p p c φ+-=
?
?
???+-=+-=?=n
s
p n
s
p p c e
p p c e
V S S dx
dx h h S S dx Cdx x III )((2
1
2
1213φφωω
其中控容体积p V ?为：
311222)()(2)()(x s s s s V s n e p ???
?
????+????????+???ω 将I 、II 、III 代入积分式，整理得到
b a a a a a S S N N w w E E p p ++++=φφφφφ
其中
SP
s
s
S NP n n N w PE e e E s s a s s a wp s s a s s a )()(,)()(,)()(),)()((
21211212δδδδω?Γ=?Γ=?=?Γ=
p C p p S N w E p V S b V S a a a a a ?=?-+++=,
前面导出的三种典型坐标系中的二维导热的离散方程只不过是上式的特殊形式，从上式出
N
a
n d
e E P W
ω s
b c
x 2
x 1
S
发，可以编制正交坐标系中二维导热的通用程序，不同坐标系中只是规定21,h h 的不同计算式即可。
（对流动问题，采用原始变量法时，亦需应用交错网格）
对于扩散（导热）问题，当计算区域的边界正好与正交曲线坐标系的坐标相吻合时，采用相当导热体的概念可以使这类问题的计算得到简化，无内热源时，能量方程为
0))(()(32332132232121213211=??????+????x t h h h h x t h h h h x x t h h h h x λλλ 或
0)(3
12
321=????
∑=i i i i
x t h h h h x λ 如果把21,x x 及3x 视为一新直角坐标系中的三个坐标，把2
1
3
21h h h h λ
作为相当导热系数，记着*i λ，则上式化为
∑==????3
1*0)(i i i i
x t
x λ 这样，在原直角坐标系中的一个曲面物体就被转换成了新直角坐标系中的规则长方体，导热系数由λ变成了各向异性的相当导热系数*i λ，原直角坐标系所代表的空间称为物理空间，而新直角坐标系所代表的则称为计算空间，如图所示：
Z x
x 1 x 2
x 1 y
x 3
x 1
x 2
计算空间
物理空间
于是，该导热问题可先在计算空间中求解，由于它在计算空间中是规则区域，数值计算易于
进行，然后把计算空间中的计算结果传送到物理空间中去，两个空间中点的对应关系由曲线坐标的定义所规定。
在由物理空间向计算空间转换时，边界条件亦应作相应的转换而成为计算空间长方体的边界条件。
第一类边界条件：把给定温度赋给计算空间中相应的点即可。 对第二、第三类边界条件，相应的模式列于表中
类别 空间
第二类边界条件
第三类边界条件
物理空间
n q n
t
=??-λ
B
B f n t
t t ??=-λ
α)( 计算空间
n i i i
q h h h h x t 3
21*=??-λ B
i
i
B f i x t
t t h h h h ??-=-*321)(λα
把物理空间中的曲线坐标系当作计算空间中的直解坐标系，从而把物理空间中不规则形状的求解区域变换成计算空间中规则的长方体（三维问题）或矩形（二维问题），这种做不
仅适用于导热问题，也可用于对流换热问题，
于是，有限差分法的一个主要弱点—不易处理
不规则边界问题，可以说原则上已根本解决了。 现在的问题是，现成的曲线坐标系的数目
毕竟是有限的，而不规则形状的物体则是千变万化的，不可能满足上述变换的需要。要使上述变换思想付诸实现，必须发展一套方法，对于那些没有现成曲线坐标系可以利用的复杂形状，可以利用这套方法来建立一个与该物体相适应的坐标系，这就是适体坐标系的思想。
5-5 适体坐标系（Body-Fitted coordinates ）
一、适体坐标系的基本概念
如上所述，最理想的坐标系是其各坐标轴与所计算物体的边界一一相符合的坐标系，称为适体坐标系，又称贴体坐标系，附体坐标系，无现成的坐标系可资利用时，希望通过计算来构造这样的坐标系，如图所示。
y x
η D ξ1
A
η1
B
ξ ξ2
C
η2 η
η2 η1 ξ1
A
ξ2 ξ B C
Δξ
Δη
x-y 坐标系中的不规则区域ABCD ，把其相交的两个边界作为曲线坐标系的两个轴，记为ηξ,，且：
（1）一条边界上，只能一个坐标单值地发生为化，另一个坐标保持为常数；
（2）在两条对应边界上，另一组曲线坐标的最大值与最小值对应相等，以便在计算平面上得出矩形区域。
问题：物理平面上，点（x ，y ）与计算平面上（ηξ,）的对应关系？求解结果的传送？
z y x x 3 x 2 x 1
f t ,*
2α const
x =2
n n q h h q h h h h q 312321*
2== α31*
2h h a = 2
x
n α,t f q n
x 2=const
如果把ηξ,视为物理平面上的两个未知函数，那么上述确定ηξ,的问题就是物理平面上的一个边值问题，因此，从物理平面上来说，所谓要生成一个适体坐标系，实际上相当于要求解物理平面上的一个边值问题。
反之，首先把物理平面上的区域ABCD 按已规定的边界上的ηξ,值，画成计算平面上的ηξ,坐标系中对应的矩形，然后以均匀网格离散化，于是问题变为：已知计算区域边界上各节点（ηξ,）相应的（ηξ,），问在计算区域内部任意一点（ηξ,）对应的（ηξ,）值是多少？这样，如把x 、y 视为计算平面的未知函数，则生成适体坐标系的问题即为计算平面中的一个边值问题。
从数值计算的观点，对生成的适体坐标系有以下几点要求：
（1）物理平面上的点与计算平面上的点一一对应，同一簇中的曲线不能相变，不同簇中的两曲线只能相交一次；
（2）在适体坐标系中，每一个节点应当是一系列曲线坐标轴的交点，而不是一群三角形元素的顶点或无序的点群，以便设计有效，经济的算法及程序，采用矩形网格即可；
（3）物理平面求解区域内网格疏密程度易于控制；
（4）在物理平面的适体坐标的边界上，网格线最好与边界线正交或接近正交，以便于边界条件的离散化。
二、生成适体坐标系的方法分类 大致可分为三类 1、复变函数法
利用复变函数的映射，可以把相当一批二维不规则区域变换成矩形区域，而且可以得出解极或部分解析的变换关系式。
y
r x
r 1
r 2
θ
θ 2π
r 2
r 1
r
η 1
1 ξ
λr =? λη
=?
λr =λr λB =Nr
)/()(121t r r r --=ξ
πη2/0=
相应地： )(ξηx
=)2cos(])([121πηξr r r -+
),(ηξy
=)2(])([121πηξxin r r r -+
复变函数映射 θθi r l re l z l W n i n n +===)(
ηξi + 所以??
?==θ
ηξr
l n
控制方程的变化：0122222=??+??+??θt
r r t r
t 因为θη=，∴2222//y t t ??=??θ
nr l =ξ，∴0/2
2=??+????=??ηξξt
r t r t …仍然是各向同性！
代入原控制方程，得
（条件0≠r ）0222
2=??+??y
t
t ξ…仍然是各向同性。 可见，正交坐标系法与复变函数法并不相同！
2、代数变换法—利有一些代数关系式来进行区域变换，具体实施方法颇多，最简单的是边界规范化的方法。
3、解微分方程的方法—通过求解边值问题的微分方程建立物理平面与计算平面上各点间的对应关系。
对该边值问题的控制方程的类型，物理问题本身无任何限定，可以比较自由地按照对所生成网格的要求来选择控制方程。
三、计算平面上求解区域形状的选择
1、物理平面上的区域由四条两两相交曲线构成的单连通域→计算平面上取为正方形或矩形；
2、物理平面的L 型区域：有两种选择
6
5 4
3
2
1
x
y 6
5 4
3 2
1
ξ η
① η
ξ
6 5 4
3
② 1
2
3、物理平面上的型区域：两种选择
特别注意：节点间一一对应关系不要受破坏，计算平面中求解节点方程时，对重迭线上的节点，每一个变量应有两套数组，以分别存放从左边和右边计算而得之值。
η 2π
λ 0
l n r 1 l n r 2
ξ
注意：计算平面上的尖角点未必对应于物理平面上的尖角点，反之亦然，例如物理
平面上一条连续封闭的曲线所围成的区域可以转换成计算平面上的一个正方形。 选择计算平面上求解区域形状的一条
件：生成的网格要适用于计算的问题，通常，
在计算平面上都采用正方形均匀网格。 四、用适体坐标系求解流动与换热问题的总体步骤
1、网格生成：在计算平面上选择与物
理平面上复杂区域相应的求解区域，并找出两个区域内部节点之间的对应关系，
),(),,(ηξηξy y x x ==或其反函数。这些关系可以是解析的，也可以是数值的。
2、控制方程的生成与离散，把物理平面上的控制方程及边界条件转换成计算平面上的相应形式，并利用控制容积积人法律立离散方程。
3、离散方程的求解及解的传递。在计算平面上获得收敛的解后，根据节点间的对应关系，传送到物理平面上去。
5-6 控制方程的转换及离散化
一、方程的转换
讨论如何把物理平面上的稳态对流一扩散方程的通用形式转换成一般曲线坐标系ηξ-中的形式，物理平面上的通用形式为
),()()()()(y x S y y x x y u x u +??Γ??+??Γ??=??+??φ
φφρφρ
y
y y x
x x ????+????=??????+
????=??ηηφξξφφηηφξξφφ 设已有函数),()
,(ηξηξy y x x ==要求找出反函数ηξ,的上述导数（用已知函数
),(),,(ηξηξy y x x ==的导数来表示），这涉及到多元函数的微分关系：
,1
1,11ξηξηηξy J y J x y J y J x -=??-=??=??=??
ξηξηηξx J
x J y x J x J y 1
1,11=??=??-=??-=??
8 7 7 8 3 4 5 6 1 2 y x η
8 7
① 3 4 5 6
1 2 3 9 4 5 2 6 1 2 5 6
7 8 1 η ξ ξ y x 3 9 4
式中，J 代有Jacobi 行列式：ξηξηξ
ηξy x y x y
y x y x J -=????-????=
∴])()([1)(1][1η
φξφφφξηφηξφφξηξηηξ??-??=-=????-????=??y y J y y J y y J x =
])()[(1
ηξξηφφy y J
-
])()([1
)(1ηξξηξηηξφφφφφx x J
x x J y +--=+-=?? 同理，对函数φρφρu u ,可以写出：
])()([1)(ηφρξφρφρξη??-??=??y u y u J x u
])
()([1)(ξ
φρηφρφρηξ??-??=??x u x u J y u 扩散项的转换： 由：)(1ξηη
φξφφy y J x ??-??=?? 有
：
?
?????-Γ??--Γ??=-Γ??=??Γ??])([])([1)](1[)(ξξηηξηξηηξξηηξφφηφφξφφφy y y J y y y J J y y J x x x
=
?
?????-Γ??--Γ??)]([)]([12
2ξηξξηηξηηξφφηφφξx x x J y y y J J 同理有：
?
?????-Γ??--Γ??-=??Γ??)]([)]([1)(22ξηξξηηξηξηφφηφφξφx x x J x x x J J y y 源项：),()],(),,([),(ηξηξηξS y x S y x S ==
两对流项相加：得])
()([1])()([1ξ
φρηφρηφρξφρηξξη??-??+??-??x u x u J y u y u J
=
])()([1})]([)]([{1η
φρξφρηρφξρφξξηη??+??=?-?+?-?v u J uy ux ux u J
式中，ξξηηuy ux V ux uy U -=-=,，可分别视为计算平面上ξ及η方向上的速度分量。 两扩散项相加：得
]})([])([{1]})([])([{12222ξηξηξηξξηηξηξξηηφηφηφξφξx x y y J
y x J J x x y y J y x J J +Γ??-+Γ??++Γ??-+Γ??
=
)]([1)]([1ηξηξφβφηβφφξx J
J l J J +-Γ
??+-Γ?? 式中222
2
,,ξξηξηξηηη
ξηξβαy x r y
y x x y y x x y x +=????+????=+=+=
计算平面上的守恒方程为：
),()]([1)]([1)(1)(1ηξφβφβφαφξφρηφρξηξηξS r J
J J J v J u J ++-Γ
+-Γ??=??+?? 可见，式中各项仍然保持了直角坐标中相应各项的意义，其中源项S 完全由直角坐标中的源项转换而来，其它各项在转换过程中并不给源项增添新的成分。
生成的曲线坐标在物理平面上是否正交？
0=β的点处，曲线坐标正交，且α和r 即为相应度规系数的平方；
0≠β的点处，曲线坐标斜交，β偏离零的多少反应了局部曲线坐标偏离正交的程度。
从物理平面到计算平面区域的简化（规则化）是以控制方程的复杂化为代价的。 二、方程的离散化
热物理过程的数值模拟-计算传热学3汇总
四、非线笥问题迭代式解法的收敛性 每一层次上满足迭代法求解的收敛条件+相邻次间代数方程的系数变化不太大（亦即未知量 的变化不太大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 利用迎风格式、混合格式、乘方格式求解φ的分布.
热物理模拟设备的发展
物理模拟设备的发展综述 摘要：物理模拟技术，作为材料成形工艺的简单实验，可以对复杂成形技术提供可靠的支持，在材料的加工领域里面有不可取代的作用。早期使用橡皮泥，铅块，石蜡等塑性较好的材料来进行复杂成形过程的模拟，以提供合理的设计参数，这种方法浪费大，时间长，效率较低，随着计算机技术的发展，目前更多的模拟同在在电脑上进行，先在热物理模拟机上进行的简单的模拟，得到材料的性能参数，然后在电脑上利用专门的商业软件进行模拟，这样不仅花费小，开发周期短，而且可以使材料的数据得到最大的用途。因此，热物理模拟设备的发展对物理模拟的进步有着举足轻重的作用。 关键词：物理模拟，热物理模拟机，Gleeble
前言 “物理模拟”是一个内涵十分丰富的广义概念，也是一种重要的科学方法和工程手段。通常，“物理模拟”是指缩小或放大比例，或简化条件，或待用材料，用实验的模型来代替原型的研究。对材料和热加工工艺来说，物理模拟通常指利用小试样，借助某种实验装置在线材料制备或热加工过程中受热火受力的物理过程，充分而准确的揭示材料或工件在制备和热加工过程中的组织和性能变化规律，用这些来评定或预测材料制备或加工过程中可能出现的问题，为制定合理的加工工艺和参数，以及研制新材料提供理论指导和技术支持。物理实验可以分为以下两种，一种是在模拟过程中进行的实验，另一种是模拟完成后进行的实验。 以往我们在进行科学研究或者工件的生产过程，为评价工艺方案对材料性能或产品质量的影响，多采用实验的方法，这种简单直接的实验不仅仅要消耗大量的时间，材料和金钱，而且得到结果仅仅能够表示在该工艺下的结果，并不能对其他工艺有太多的指导意义，因此我们必须在实验工艺和方法上进行有一定的创新和改造。 近些年来，随着计算机技术和工程检测技术的迅速发展，物理模拟，数值模拟以及与模拟相关的专业软件都有了长足的进步，相关软件在材料科学和工程领域的运用都取得了非常好的效果，材料学科的研究开始从“经验”走向“科学”。新模拟技术的应用使得人们不仅可以对变形过程有了更加直观的认识，对模具的设计参数好坏有了更加直观的评价，为工艺的制定和工艺参数的设计提供了更加可靠的依据，从而大大减少了新产品和新材料的开发周期和开发费用，降低了企业的成本，提高企业的竞争力。
传热学数值计算大作业2014011673
数值计算大作业 一、用数值方法求解尺度为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 二维导热物体温度场的数值模拟 一、物理问题 有一个用砖砌成的长方形截面的冷空气通道， 于纸面方向上用冷空气及砖墙的温度变化很小， 可以近似地予以忽略。 在下列两种情况下试计算： 砖墙横截面上的温度分布；垂直于纸面方向的每 米长度上通过砖墙的导热量。 第一种情况：内外壁分别均匀维持在 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 取步长δ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、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、时间项差分的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. 传热学的发展概述 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 为物体周围流体的温度。一个 物体或两个物体之间通过电磁波形式进行的能量传递交换称为热辐射，通常由斯 粘土的热物理参数和冻结过程中的温度场演变研究冻土是一种由固、液、气、冰四相物质组成的复合材料。冻土中的冰相对其物理、力学和工程性质具有决定影响。 因此,研究冰相的存在形式和形成过程,是揭示冻土复杂性质的基础和条件。而认识冻土的热参数及其导热特性,则是揭示冻土形成规律必不可少的条件。 本论文从冻土的材料组成入手,研究了冻结过程中的热量构成,建立了冻土 热物理参数的计算新模型,研发了室内模型试验装置,开展室内模型试验,通过模型试验结合数值模拟分析,验证了本文提出的冻土热物理参数计算新模型的有效性。为寒冷地区的工程建设和冻结法施工,提供了设计理论基础。 获得的成果主要包括以下几个方面。(1)揭示了土在冻结过程中的热量构成。 将未冻土当作颗粒、孔隙水和孔隙气组成的复合体,则冻结过程实质上是孔隙水的降温过程、孔隙水的结冰过程、孔隙冰的降温过程、土颗粒的降温过程共四个物理过程的组合。因此,土在冻结过程中的热量消耗主要用于颗粒的温度变化、水(冰)的温度变化和冰水相变。 (2)揭示了土在冻结过程中,比热和导热系数的演变特点。由于冰和水的含量在冻结阶段是不断变化的,而冰的比热和导热系数完全不同于水的比热和导热系数,因此土体的比热和导热系数是一个动态变化过程。 因此,研究土冻结阶段的导热系数和比热,必须与孔隙水的相变过程和相应 潜热相联系。(3)根据是否存在冰水相变或潜热,提出了一种冻结过程的划分方法。 随温度降低,若土中水一直以液态形式存在则称之为未冻阶段;若存在固态 冰且液态水随温度降低逐渐减少则称之为冻结阶段;若存在固态冰但液态水不再随温度降低而减少则称之为冻实阶段。(4)建立了完整的冻结过程三阶段比热计 热物理过程的数值模拟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型钢淬火模型的编程计算与数据处理等。 以下是我们小组各组员的任务完成情况： 以下是我们小组各成员的排序情况： 传热学经典计算题 热传导 1. 用热电偶测量气罐中气体的温度。热电偶的初始温度为20℃，与气体的表面传热系数为()210/W m K ?。热电偶近似为球形，直径为0.2mm 。试计算插入10s 后，热电偶的过余温度为初始过余温度的百分之几？要使温度计过余温度不大于初始过余温度的1％，至少需要多长时间？已知热电偶焊锡丝的()67/W m K λ=?，7310ρ= 3/kg m ，()228/c J kg K =?。 解： 先判断本题能否利用集总参数法。 3 5100.110 1.491067hR Bi λ--??===?＜0.1 可用集总参数法。 时间常数 3 73102280.110 5.563103c cV c R hA h ρρτ-??===?= s 则10 s 的相对过余温度 0θθ=exp c ττ??-= ???exp 1016.65.56??-= ???％ 热电偶过余温度不大于初始过余温度1％所需的时间，由题意 0θθ=exp c ττ??- ??? ≤0.01 exp 5.56τ?? - ???≤0.01 解得 τ≥25.6 s 1、空气以10m/s 速度外掠0.8m 长的平板，C t f 080=,C t w 030=，计算 该平板在临界雷诺数c e R 下的c h 、全板平均表面传热系数以及换热量。 （层流时平板表面局部努塞尔数 3/12/1332.0r e x P R Nu =，紊流时平板表面局部努塞尔数3/15/40296.0r e x P R Nu =，板宽为1m ，已知5105?=c e R ，定性 温度C t m 055=时的物性参数为： )/(1087.22K m W ??=-λ，s m /1046.1826-?=ν，697.0=r P ） 解：(1)根据临界雷诺数求解由层流转变到紊流时的临界长度 C t t t w f m 055)(21=+=,此时空气得物性参数为： )/(1087.22K m W ??=-λ，s m /1046.1826-?=ν，697.0=r P )(92.0101046.1810565m u R X ul R c c e c e =???==?=-ν ν 由于板长是0.8m ，所以，整个平板表面的边界层的流态皆为层流 ? ==3/12/1332.0r e x P R hl Nu λ)/(41.7697.0)105(8.01087.2332.0332.023/12/1523/12 /1C m W P R l h r e c c ?=????==-λ (2)板长为0.8m 时，整个平板表面的边界层的雷诺数为： 561033.41046.188.010?=??==-νul R e 全板平均表面传热系数: )/(9.13697.0)1033.4(8.01087.2664.0664.023/12/1523/12 /1C m W P R l h r e c ?=????==-λ 全板平均表面换热量W t t hA w f 9.557)3080(18.09.13)(=-???=-=Φ西安交通大学传热学大作业二维温度场热电比拟实验1
计算传热学中国石油大学(华东)第四章大作业
计算传热学数值模拟
热物理过程数值模拟考试题(2009冉)
传热学大作业报告 二维稳态导热
计算传热学-传热基本原理及其有限元应用
粘土的热物理参数和冻结过程中的温度场演变研究
热物理过程的数值模拟-计算传热学1
传热学计算例题
热物理过程的数值模拟-计算传热学2
西安交通大学传热学大作业
传热学数值模拟
传热学经典计算题