金属凝固计算机模拟b

1前言１、材料加工过程铸造、锻造、焊接、挤、压、拉、拔、热处理。专业设置。（大学、研究生、专业学位）2、计算机在材料加工中应用领域有五个方面：（1）科学计算①数值模拟②优化计算如最低成本配料优化计算，其是应用线性规划解法，求满足化学成分要求范围内，在已知化学成分和价格的几种或几百种原料中，同时还可考虑某种原料最多，最少或必须用多少的条件，求出最低成本的配料单来。也可以用此方法求出最高成本的配料单。传统手算配料单，由于其不考虑成本，也无法考虑成本，因此，传统手算的配料成本落在最高和最低成本之间。234567前言（2）信息处理信息处理的特点是计算简单、数据量大。现在75～85％的计算机用于信息处理。①数据库数据库（DATABASE）的概念由美国军方于1963年提出。生产、经营、人事、工资等管理都可应用数据库管理。材料加工的生产、经营、人事、工资等管理都可应用数据库管理。8前言②专家系统第一个专家系统是1965年由美国斯坦福大学开发的，用质谱仪得到的数据来决定一个未知化合物分子结构的程序。专家系统是将专家的知识和经验按一定规则存在知识库或数据库里，用户则可通过专家系统，查寻检索到存在专家系统内的专家的知识和经验。材料加工的计算机专家系统正在研究之中，如“铸造缺陷分析专家系统”，“辊型设计的专家系统”。910前言斯坦福大学（StanfordUniversity）创建于1885年，位于加州圣克拉拉县。总面积超过8000英亩，校园约5700英亩，住宅约500英亩，剩下大半个校园，除了一些缓冲地段之外，就是植物园、高尔夫球场和若干个科学实验场。校园构建誉为美国城市规划和庭院设计的典范。紧凑、集中、和谐。主要建筑：土黄色石墙环绕下的红屋顶建筑，拱廊相接，棕榈成行。最有名的建筑是斯坦福纪念教堂。校园内公路长达46英里，三个水坝和湖。1英亩=6.070市亩，1英里=1.609公里11前言斯坦福的工业园是闻名遐迩的“硅谷”诞生地，在64名硅谷最重要的先驱中，有28人是斯坦福大学的校友、教授或研究人员，占42%。7个学院（地球科学、工、文理、商、教育、法、医）专职教师1714人，从建校到现在已有25位诺贝尔奖得主。现有教师中，有17位诺贝尔奖得主，4位普利策奖得主，23位麦克阿瑟奖得主，21位国家科学奖得主，3位国家技术奖得主，220位美国文理科学院院士，128位国家科学院院士，83位国家工程院院士，25位国家教育科学院院士，42位美国哲学社会科学学会成员，6位Wolf基金奖获得者，6位Koret基金奖获得者，2位总统自由奖获得者。12前言（3）计算机辅助设计和制造。计算机辅助设计，简称CAD。计算机辅助制造，简称CAM。近年来由于在计算机上进行设计，就连设计图的概念也有所改变，如在波音767飞机设计中，图纸定义为点列的集合，然后从美国把录有这些点列数据的磁盘寄到日本等国即可生产。输出设计图只是为了便于人的理解，而并非照图加工。机械加工时全部使用数控加工机械。材料加工工艺设计的CAD，材料加工工装制造的CAM，材料加工车间设计的CAD系统正在研究发展之中。由于材料加工的复杂性，材料加工车间的多样性和研究发展水平限制，现在材料加工的CAD／CAM只局限于某一类型的CAD／CAM。13前言现代英国铸造的CAD／CAM技术发展到每周能够出图18000多幅。日本日立公司研究的设计金属压铸金属型CAD系统，能够在获得产品数据之后，启动画出型腔图，并对模型进行计算，输出用于CAM的数据磁盘。

（4）计算机测试和控制系统计算机和检测仪表、控制部件结合，可形成计算机测试、控制系统。材料加工生产中的测试和控制问题，一般地说，都可以采用计算机测试、控制系统。（5）机械手和机器人14前言２、数值模拟两个不同领域的现象，能用同一数学形式来描述，称这两个现象彼此是可模拟。模拟的方法是把一个领域内求解的问题过渡到另一个领域中去解决。随着计算机和计算数学的发展，用计算机数值计算法求解问题的计算精度已经达到很接近解析解，此法称为计算机数值模拟。仿真的定义，与模拟的差异。３、工艺优化计算机数值模拟的最终目的是解决工艺优化设计问题。一旦全面实现了材料加工过程的计算机数值模拟，材料的加工生产将会产生深刻的变革。15前言４、计算机数值模拟的步骤

①给定研究对象几何条件、物理条件、初始条件、边界条件

②离散化处理将过程所涉及的区域在空间和时间上进行离散化处理。

③建立数值方程建立内节点和边界节点的数值方程。

④选择计算方法选择适当的计算方法求解线性代数方程组。

⑤编制计算机程序编制计算机程序，由计算机算出结果，并用数据、曲线、图形输出。

⑥优化工艺分析计算机模拟的结果，如果结果不理想，调整工艺参数，再进行计算机模拟，直到模拟结果为最佳结果。16前言５、讲解的主要内容通过一、二个例子，让大家亲自经历计算机数值模拟的一般步骤，从中掌握其基本原理和方法，为更好应用计算机数值模拟，优化工艺参数打下基础。材料加工过程中，比较常用、典型的数值模拟有：

①温度场数值模拟

②浓度场数值模拟

６、教学安排和考核

①讲课18学时

②自己上机练习4学时

③考试（笔试开卷）2学时

④成绩评定：Ⅰ笔试开卷（70%）；

Ⅱ交上机的程序及运行结果（30%）17181温度场计算机数值模拟1.1传热的基本知识1.1.1传热的基本方式

①导热导热属于接触传热，是连续介质就地传递热量，没有各部部分物质之间宏观的相对位移。在不透明固体实体内部，由于各部分物质之间无法作宏观的相对位移，不透明无法传递辐射能，实体保证接触，所以只能依靠导热方式传递热量。导热的基本定律是傅立叶定律，即导热的热流密度q与温度梯度成正比，即

Ｔ

q∝gradT=-λgradT=-λ——

(W/m2

n

式中：q—传热方向上单位时间、单位面积上所通过的热量，[J/(s.m2)]=W/m2；

λ—材料的导热系数，W/(m.K)；

Ｔ

———温度梯度，K/m。

n负号表示导热的热流量向温度低的方向传递，即与温度梯度的方向相反。

热流密度是个向量，即它有大小和方向。

191.1.1传热的基本方式如果热流密度的分量和（X、Y、Z）坐标系相联系，那么X、Y、Z方向的热流密度分量应是：

Ｔ

Ｔ

Ｔ

qx=-λ——,qy=-λ——,qz=-λ——

XYZ

热流密度q=iqx+jqy+kqz导热系数物理意义：沿导热方向的单位长度上，温度减低１℃,物质所容许通过的热流量。方向性：大多数液体和固体属于各向同性的物质。各向异性材料的导热系数具有方向性，如石墨。温度函数：λ值还随温度而变化。大多数金属的导热系数随温度的升高而降低。大多数液体（水和甘油除外）其导热系数随温度的升高而降低。

气体的导热系数随温度的升高而增加。

201.1.1传热的基本方式②对流对流是流体（气体和液体）中温度不同的各部分相互混合的宏观运动引起热量传递的现象。工程上最具有实际意义的是：相对运动着的流体与所接触的固体壁面之间的热量交换过程，一般称为对流换热。工程上在研究固体壁面和流体之间的对流换热时，除了高度稀薄的气体外，人们不去注意流体的单个质点，而把流体看成是连续介质。实际的流体总有粘性，流动时，受粘性和壁面摩擦的影响，在靠近壁面附近的流体将降低流速，在壁面上完全被滞止不动，即X=0时，Ｖ=0，如图1.1所示。因此，热量从壁面传给贴壁的那部分流体，将依靠导热。T(K)V(m/s)TwqwcTf

VX(m)δ图1.1邻近壁面的流体速度分布和温度分布211.1.1传热的基本方式流体和固体壁接触面上的“相”际热流密度为：

Ｔ

qwc=-λf——

X式中：qwc—热流密度，W/m2；

λf—流体的导热系数，W/(m.K)；

T—流体的温度，K。液体的导热系数较小，气体的导热系数更小，所以受热时，在贴壁处的流体温度势必沿着Ｘ轴的反方向急剧升高。随着离壁面的距离Ｘ的增加，流体的流动将带走更多的热量，使Ｘ轴方向的温度梯度连续下降，温度分布趋向平坦化。正是通过这种导热和对流的共同作用，使热量在流体内部得到传播，越临近壁面，导热越起主导作用。图1.1所示，假如厚度为δ的贴壁静止膜，膜内温度线性地从壁面温度TW降到远离壁面，尚未被加热的流体温度Tf，则

Tf-TW

qwc=-λf——

δ

无界对流时壁面与流体的换热，钢锭与周围空气的对流换热属于这种情况。221.1.1传热的基本方式流体在管和槽道内部的流动，称为有界对流，这时不存在远离壁面，尚未被加热的流体温度，则采用截面平均温度作为流体的特征温度Tf，则

qWC=aC(TW-Tf),W/m2

这就是所谓的“牛顿冷却定律”。式中：aC为放热系数，W/(m2.K)。其实“牛顿冷却定律”并不是表达对流换热现象本质的物理定律。凡能影响流体流动的各种因素，包括流体的种类和状态、运动的速度、流道的形状和大小，以及固体壁面粗糙度等，都会对aC值产生影响。式只不过形式地把放热过程的一切复杂性和计算上的困难，都转移到并且集中在放热系数这样一个物理量上罢了。aC代表流体和所接触的固体表面之间温度每相差１℃，该流体与表面之间“相”际热流量的大小。运用式可以进行对流换热的计算。但由于对流换热的复杂性，该式中的放热系数aC需从相应的准则方程式求出。准则方程式是针对不同对流换热情况，在综合实验结果的基础之上，运用相似理论将表征某现象的物理量整理成一些相似准则，用因次分析法得到的各种类型的表达式。连铸钢坯二冷区23

③辐射只要温度高于绝对零度，任何物体都随时向周围空间辐射能量。辐射热流密度用斯蒂芬—玻耳兹曼定律表达：

E=εσ０Ｔ４，w/m2

式中：ε—物体的辐射率或黑度,(0-1)；

σ０—斯蒂芬—玻耳兹曼常数或绝对黑度的辐射常数；

5.67×10-8W/m2K4

Ｔ—温度,K。实际上，辐射往往涉及二个物体间辐射热交换。如果二个物体的表面温度不同，中间被空气所隔开，T1〉T2时，则相互辐射的结果，表面温度为T1的物体发射出去的辐射热超过了吸收来自表面温度为T2的物体辐射热，引起辐射换热的热流量则为：

Q1→2=ε12σ０(Ｔ1４-Ｔ2４)F1φ12或

q1→2=ε12σ０(Ｔ1４-Ｔ2４)φ12

式中：ε12

—物体1与2综合黑度；

F1—物体1的表面积；

φ12—物体１的表面向外发射出去的辐射热量中，能投射到物体２表面上的份额，称为角系数,(0-1)。1.1.1传热的基本方式24轧辊喷淬时的工作环境

25喷淬机的简易俯视图

26红外测温仪+吹扫装置

1：红外测温仪；2：固定支架；3：瞄准管；4：进气管：5：瞄准管前端27辊身不同深度处的冷却曲线28喷雾过程辊身处表面综合换热系数与时间的关系

291.1.1传热的基本方式壁面在气体中冷却，存在辐射换热和对流换热。考虑到壁面与气体之间存在着辐射换热，其热流量（或热流密度）为

Qr=arF(Tw-Tf),w或

qwr=ar(Tw-Tf),w/m2式中：qwr—单位时间、单位面积辐射的热量（辐射热流密度），w/m2；

ar—辐射放热系数，w/(m2.k)；

Tw—辐射物体表面温度，k；

Tf—透明的气体介质的特征温度，k。考虑到壁面与气体之间还存在着对流换热，其热流密度为

qwc=ac(Tw-Tf),w/m2

由壁面传走的总热流密度qw应是qwr和qwc二者之和

qw=ac(Tw-Tf)+

ar(Tw-Tf

令a0=ac+ar

则qw=a0(Tw-Tf),w/m2

用辐射放热系数ar，可以形式地把辐射换热折合成对流换热，用总放热系

301.1.1传热的基本方式

数a0兼顾辐射换热的影响，从而有利于简化复杂传热的分析和计算。如远离表面的外界表面温度趋于环境温度Tf，并且φ12=1时，由式得

qwr=ε12σ０(Ｔw４-Ｔf４)

由式得：

qwr=ar(Ｔw-Ｔf

由式和式得：

ε12σ０(Ｔw４-Ｔf４)=ar(Ｔw-Ｔf

因Ｔw４-Ｔf４=(Ｔw-Ｔf)(Ｔw3+Ｔw2Ｔf

+ＴwＴf2+Ｔf3)

设Ｔm=1/2(Ｔw+Ｔf)，(Ｔw3+Ｔw2Ｔf

+ＴwＴf2+Ｔf3)≈4Ｔm3

ar≈4ε12σ０Ｔm3

,w/(m2从此可见，(Ｔw-Ｔf)降为零时，ar并非零值，而是以4ε12σ０Ｔm3。随着温差的扩大和平均温度Ｔm的升高，ar值将迅速增加。由于ac随温差的变化较小，在高温范围和大温差情况下，ar有可能成为a0的主要组成部分。ar与气体的运动状况无关，而ac与气体速度的降低而减小，在气体自然对流的情况下，ac≈5,w/(m2.k)，即使Ｔm只有300K，ar就可能占a0的一半左右。

311.1.2导热的偏微分方程式钢锭一般属于各向同性的物质，其加热或冷却过程数学模拟计算依据的基本数学模型，是不稳定导热偏微分方程。下面讨论各向同性材料导热方程式的建立。

①直角坐标系导热偏微分方程

导热偏微分方程的建立，是通过考察处于导热过程中的物质的微元体积（δxδyδz)的能量平衡来建立。如图1.2所示。在时间δ内，通过六个面的导热所获得的能量，加上微元体内产生的内热源热量，要等于微元体积内物质积蓄热量的改变，即温度升高或降低。xzydQx+δxdQxdQzdQydQz+δzdQy+δy图1.2直角坐标系导热方程式的微元体321.1.2导热的偏微分方程式图1.2中，微元控制体尺寸δx、δy、δz，按照傅立叶导热定律，在X方向流入微元体左表面的热流可表示为：

T

dQx=-λ(δyδz)—，

X

式中：λ—

导热系数，W/(m.k)；

δyδz—X方向微元体表面积，m2；

T

—

—

X方向的温度梯度，k/m。

X在X方向流出微元体右表面的热流可表示为：

T

dQx+δx=-λ(δyδz)—，W

X

在X方向导热的净热流为：?

dQx-dQx+δ=0331.1.2导热的偏微分方程式图1.2中，微元控制体尺寸δx、δy、δz，按照傅立叶导热定律，在X方向流入微元体左表面的热流可表示为：

T

dQx=-λ(δyδz)—，

X

式中：λ—

导热系数，W/(m.k)；

δyδz—X方向微元体表面积，m2；

T

—

—

X方向的温度梯度，k/m。

X在X方向流出微元体右表面的热流，可以应用泰勒级数展开，保留级数的第一项和第二项而得到：

T

dQx+δx=dQx+—(dQx)δx=dQx+—[-λ(δyδz)—]δx

XXX

T=dQx-—(λ—)δxδyδz

XX在X方向导热的净热流为

TdQx-dQx+δ=—(λ—)δxδyδz

XX

341.1.2导热的偏微分方程式用同样的方法，可以得出Y,Z方向与式1.１.19相类似的导热的净热流方程式，即

TdQy-dQy+δy=—(λ—)δxδyδz

YY

TdQz-dQz+δz=—(λ—)δxδyδz1.

1.21

ZZ在三个坐标方向净热流的总和为：

TTT

[—(λ—)+—(λ—)+—(λ—)]δxδyδz1.1

.22

XXYYZZ如果单位时间、单位空间所产生的热量为Q’(x,y,z,τ)，那么微元体单位时间的发热量（热流）为：

Q’(x,y,z,τ)δxδyδz

由于导热传进微元体的净热流式和微元体内产生的热量式一起用于增大微元体的内能。微元体的内能的增加表现在微元体能量存储随时间的变化率，即

TρCpδxδyδz—

351.1.2导热的偏微分方程式

式中：ρ—密度，kg/m3；Cp—比热，J/(kg.k)；—时间，s。

对微元体进行能量平衡，使能量存储的时间变化率与由导热引起的流入微元体的净热流和微元体内产生的热流之和相等，得下式：

TTTT

ρCp—=—(λ—)+—(λ—)+—(λ—)+Q’

XXYYZZ②圆柱坐标系导热偏微分方程实际问题常常涉及柱面对称问题，且边界条件给定在一个表面上，因此表面具有一个坐标保持不变的性质。在这种情况下，采用圆柱坐标系是合适的。图1.3圆柱坐标系导热方程式的微元体361.1.2导热的偏微分方程式图1.3所示圆柱坐标系，直接按内外两个圆弧面和其它四个平面组成的微元体积，在导热过程中热量必须按收支平衡的原则导出，此时，微元体积为：

dv=(rdφ)dzdr1.1.26沿内圆弧面流入微元体积的热流：

T

dQr=-λ—(rdφdz)

r沿外圆弧面流出微元体积的热流：

T

dQr+dr=dQr+—(dQr)dr=dQr+—[-λ—(rdφdz)]dr

rrr

T

1T

=dQr+—(-λr—)dφdzdr=dQr-—

—(λr—)dv

rrrrr沿φ平面流入微元体积的热流：

T

dQφ=-λ——

(drdz

rφ沿φ+dφ平面流出微元体积的热流：

T

dQφ+dφ=dQφ+—(dQφ)rdφ=dQφ+—[-λ—(drdz)]rdφ

rφ

rφrφ

371.1.2导热的偏微分方程式

1

Ｔ

=dQφ-—

—(λ—)dv

r2φφ沿z面流入微元体积的热流：

T

dQz=-λ—(rdφdr

z沿z+dz面流出微元体积的热流：

T

dQz+dz=dQz+—(dQz)dz=dQz+—[-λ—(rdφdr)]dz

zzz

T

=dQz-—(λ—)dv

zz根据直角坐标系导热微分方程推导的思路，可得到：

1T1ＴTT

—

—(λr—)+—

—(λ—)+—(λ—)+Q’=ρCp—

rrrr2φφzz381.1.2导热的偏微分方程式

③球坐标系导热偏微分方程对于球坐标系，如图1.4所示。图1.4球坐标系导热方程式的微元体由内、外两个球面、两个圆弧面和两个平面所组成的微元体积为：

dv=(rdθ)（rSinθdφ)dr391.1.2导热的偏微分方程式沿半径方向流入微元体积的热流：

T

dQr=-λ—

r沿半径方向流出微元体积的热流：

T

dQr+dr=dQr+—(dQr)dr=dQr+—[-λ—(rSinθdφ.rdθ)]dr

rrr

T

1T

=dQr+—(-λr2—)Sinθdφ.dθdr=dQr-—

—(λr2—)dv1.1.36

rrr2rr沿θ方向流入微元体积的热流：

T

dQθ=-λ——

(rSinθdφ.dr)1.1.37

rθ沿θ方向流出微元体积的热流：

T

dQθ+dθ=dQθ+—(dQθ)rdθ=dQθ+—[-λ—(rSinθdφdr)]rdθ

rθ

rθrθ

T

T

=dQθ-—(λ—Sinθ)rdφdr.rdθ=dQθ-————(λ—Sinθ)dv

r2θθSinθr2θ

θ

401.1.2导热的偏微分方程式沿φ方向流入微元体积的热流：

T

dQφ=-λ—————

(dr.rdθ)1.1.39

(rSinθ)φ沿φ方向流出微元体积的热流：

dQφ+dφ=dQφ+——————(dQφ).(rSinθ)dφ

(rSinθ)φ

T

=dQφ+——————(-λ—————dr.rdθ)(rSinθ)dφ

(rSinθ)φ

(rSinθ)φ

T

=dQφ-——————(λ——)dv

(r2Sin2θ)φ

φ同理可整理得：

TＴ

TT

——(λr2—)+—————(λSinθ—)+—————(λ—)+Q’=ρCp—

r2rrr2Sinθθ

θr2Sin２θφφ

1.1.41

411.2导热方程的有限差分解法求解不稳定导热偏微分方程的数值解法，主要有：有限差分解法、有限元法、边界元法三类。边界元法正在研究和完善之中，目前常用的是有限差分解法和有限元法。我们这门课专门讨论有限差分解法的数学基础，数值方程的建立，差分方程的稳定性和收敛性等问题。有限差分解法是用差分方程近似地代替微分方程，建立差分方程有直接法和能量平衡法两种。直接代换法直接代换法是从微分形式出发，用差商直接代换微商的办法建立差分方程。数学基础

①微商和差商的定义若T(x)是连续函数，则其导数为：

dTＴ(x+Δx)-T(x)ΔＴ

—=lim——————=lim——

dxΔx→0Δx

Δx→0Δx42数学基础式右边ΔＴ/Δx是有限的差商，Δx和ΔＴ都不为零。式左边dＴ/dx是ΔT/Δx当Δx→0时的差商，称之为微商。在Δx没有达到零之前，ΔT/Δx只是dT/dx的近似。实际应用Δx≠0。如果把ΔT/Δx趋于dT/dx的过程认为是近似向精确过渡，那么，用ΔT/Δx代替dT/dx，则两者的差值│ΔT/Δx-dT/dx│表示差商代替微商的偏差。误差多大？需要做误差分析，才能大胆地应用。

②误差分析假设函数T(x)在x时的值为T(x)，在x+Δx时的值为T(x+Δx)，如果函数T(x)在x处的各阶导数存在，则按照泰勒级数展开，T(x)与T(x+Δx)的关系如下式所示：

dT(Δx)2d2T(Δx)ndnT

dx

2!dx2n!dxn整理后得：

ΔTT(x+Δx)-T(x)dTΔxd2T(Δx)n-1dnT

ΔxΔxdx

2!dx2n!dxn从上式可知，T(x)在x处的差商ΔT/Δx等于函数T(x)在x处的各阶导数的线性组合，只能是近似地等于差商。两者之间也必然有偏差。43数学基础图表示了一阶差商与一阶微商之间的关系。用不同方法得到的差商去代替微商，它们带来的误差是不同的。即向前差商：dT/dx≈[T(x+Δx)-T(x)]/Δx

向后差商：dT/dx≈[T(x)-T(x-Δx)]/Δx中心差商：dT/dx≈[T(x+Δx)-T(x-Δx)]/2Δx图差商与微商44数学基础按照泰勒级数展开，T(x)与T(x+Δx)的关系如下式所示：

dT(Δx)2d2T(Δx)ndnT

dx

2!dx2n!dxn整理后得：

T(x+Δx)-T(x)dTΔxd2T

Δxdx

2!dx2即向前差商的偏差是截去了泰勒级数展开式中的高阶项而引起的，常称“截断误差”，其截断误差为与Δx同级的小量O(Δx)。同理

dT(Δx)2d2T(Δx)3d3T

dx

2!dx23!dx3整理后得：

T(x)-T(x-Δx)dTΔxd2T

——————-—=-————+……=O(Δx)1.2.10

Δxdx

2!dx2即向后差商的截断误差为与Δx同级的小量O(Δx)。45数学基础由式减式，将2ΔxdT/dx移至等式左边，两边再除以2Δx，得：

T(x+Δx)-T(x-Δx)dT(Δx)2d3T

—————————-—=————+……=O[(Δx)2

2Δxdx

3!dx3即中心差商的截断误差为与(Δx)2同级的小量O[(Δx)2]。当Δx固定时，上述一阶差商一般仍是x的函数，对它们还可以求差商。这种一阶差商的差商称为二阶差商，它是二阶微商的近似，常用向前和向后差商来二阶微商，即

d2T

T(x+Δx)-T(x)

T(x)-T(x-Δx)

——≈[——————-——————]/Δx

dx2

Δx

ΔxT(x+Δx)-2T(x)+T(x-Δx)

(Δx)2由式和式相加，经简化后得：

T(x+Δx)-2T(x)+T(x-Δx)

d2T

———————————-——=O[(Δx)2

(Δx)2dx2即二阶差商的截断误差为与(Δx)2同级的小量O[(Δx)2]。461.2.1.2建立内节点差分方程

①一维系统假定有一高宽无限(即高宽方向上无热流)，厚度为L的平板，T=f(x,τ)即温度是x方向位置和时间的函数，所谓一维系统是指几何空间为一维。初始时刻τ=0，T=T0，为了简化，考虑无内热源，λ、ρ、Cp均为常数。选取网格点间距Δx和时间步长Δτ将研究对象离散化。显式差分格式

TT一维不稳定导热方程为

ρCp—=—(λ—)

XX

该方程在区域τ>0，0<x<L内全部点都成立，如图所示。将方程应用于内节点i可写成：

Tp2Tp

ρCp(—)=λ(—iX2i上式等号左端的微分式用温度对时间的一阶向前差商来近似，即：

TpTp+1i-Tpi

(—)=————+O(Δτ)i

Δτ上式等号右端的微分式用温度对空间的二阶差商来近似，即：

２TpTpi+1-2Tpi+Tpi-1

(——)=——————+O((Δx)2)

x2

i(Δx)2471.2.1.2建立内节点差分方程图一维显式差分将式和式代入式得：

Tp+1i–TpiTpi+1-2Tpi+Tpi-1

ρCp[————+O(Δτ)]=λ[——————+O((Δx)2)]

Δτ(Δx)2若在上式中去掉O(Δτ)和O((Δx)2),整理得：

Tp+1i=Tpi+F0(Tpi+1-2Tpi+Tpi-1)

式中:F0=

λΔτ/[ρCp(Δx)2]=导热速率/蓄热速率，称F0为傅立叶数。481.2.1.2建立内节点差分方程F0的数值小意味着加热或冷却此物体所需要的时间长，反之，所需要的时间短，F0是一个无因次数字。当初始条件和边界条件已知时，用式就可模拟区域内各节点随时间τ增长的温度值Tpi(i=2,3,…，N-1；p=1,2,3,…）隐式差分格式

一维隐式差分如图所示，将一维不稳定导热微分方程应用于内节点i，则：

Tp2Tp+1

ρCp(—)=λ(—iX2i式和相比，式的左边完全一样，温度对时间的一阶偏微商，仍用一阶向前差商来近似，而式和右边有所不同。式中温度对距离的二阶偏微商是对应时刻p的，在用差商近似微商时，用p时刻的T值；而式中，温度对距离的二阶偏微商是对应时刻p+1的，用差商近似微商时，用p+1时刻的T值。即式相应的差分方程为：

Tp+1i–TpiTp+1i+1-2Tp+1i+Tp+1i-1

ρCp[————+O(Δτ)]=λ[———————+O((Δx)2)]

Δτ(Δx)2若在上式中去掉O(Δτ)和O((Δx)2),整理得：

Tp+1i=Tpi+F0(Tp+1i+1-2Tp+1i+Tp+1i-1)

491.2.1.2建立内节点差分方程

式中:F0=

λΔτ/[ρCp(Δx)2]。比较式和，可以看出显式差分格式的突出优点是每个节点方程都可以独立求解，整个计算过程十分简便。但是，Δτ的选取要受到限制，有时为了满足差分格式稳定性条件，Δτ可能选得很小，使计算工作量加大。隐式差分格式的最大优点是，它对任意F0值都是稳定的。这种稳定是绝对的，即不受边界条件、Δτ、Δx、热物参数的影响，于是Δτ可以选的较大，计算速度加快。但是，对于节点i，只从式不能独立求解，必须涉及Tp+1i+1、Tp+1i、Tp+1i-1的联立线性代数方程组才能求解，也就是说，它含有三个未知数。时间步长每前进一步，从坐标0≤x≤L，网格点1≤i≤N整个区域的每个点，上述方程都要列出一次（见图1.2.3)。因此，向一个新的时间步长每移动一步就必须解一个方程组。当按顺序列出这些方程时，除了要的第一点和最后一点的方程只有两个未知数外，其余每一个方程都含有三个未知数，于是方程是三对角的。对于这种情况，已有适用的求解程序，后面将讨论。关于差分方程稳定性将在后面讨论。501.2.1.2建立内节点差分方程图一维隐式差分

②二维系统显式和隐式差分格式建立的方法和两种差分格式的特点在前面讨论过，对二维系统同样适用。为简略起见，在此只讨论二维系统显式差分方程的建立。为使问题简化，仍然假定热物性值为常数，考虑无内热源。首先把所讨论的区域离散化，如图所示。×××××

×××

×××

×××

×××

×××

×××

××

511.2.1.2建立内节点差分方程网线：用平行于X、Y轴的直线（网线），进行空间离散化节点：网线与网线的交点步长：节点间的距离图1.2.4二维均质物体的分割Δx与Δy可以是不变的常量，即等步长，也可以是变量（即在区域内的不同处是不同的），即变步长。如果区域内各点处的温度梯度相差很大，则在温度变化剧烈处，网格布得密些，在温度变化不剧烈处，网格布得疏些。至于网格多少，步长取多少为宜，要根据计算精度与计算工作量等因素而定。对于0<x<L1和0<y<L2的矩形区域内，将二维不稳定导热方程式应用于节点（i,j)可写成：

Tp2T

2Tp

ρCp(—)=λ(——+——i,jx2

y2i,j521.2.1.2建立内节点差分方程式左边一阶偏微商用温度对时间一阶向前差商近似：

TpTp+1i,j–Tpi,j

(—)=—————+O(Δτ)i,j

Δτ式右边的二阶偏微商所对应的差商可表示为：

２TpTpi+1,j-2Tpi,j+Tpi-1,j

(——)=————————+O((Δx)2)

x2

i,j

(Δx)2２TpTpi,j+1-2Tpi,j+Tpi,j-1

(——)=————————+O((Δy)2)

y2

i,j(Δy)2当Δτ、Δx、Δy较小时，忽略O(Δτ)、O((Δx)2)、O((Δy)2)项，当Δx=Δy时，即x、y方向网格划分步长相等。将上三式代入式则得到节点(i,j)的差分方程：

Tp+1i,j=Tpi,j+F0(Tpi+1,j+Tpi-1,j+Tpi,j+1+

Tpi,j-1-4Tpi,j

式中:F0=

λΔτ/[ρCp(Δx)2]。上述条件，隐式差分格式，则除Tpi,j项外，其它项的Tp改为Tp+1。531.2.1.2建立内节点差分方程

③三维系统和二维系统的假定相同，热物性值为常数，无内热源。设Δx=Δy=Δz，即均匀网格间距，将三维不稳定导热方程式应用于节点(i,j,k)，可以写成：

Tp2T

2T2Tp

ρCp(—)=λ(——+——+——i,j,kx2

y2z2i,j,k即可得到节点(i,j,k)的近似式

Tp+1i,j,k–Tpi,j,kTpi+1,j,k-2Tpi,j,k+Tpi-1,j,k

ρCp(———————)=λ[———————————+

Δτ(Δx)2Tpi,j+1,k-2Tpi,j,k+Tpi,j-1,kTpi,j,k+1-2Tpi,j,k+Tpi,j,k-1

———————————+———————————(Δy)2(Δz)2上式经简化可得节点(i,j,k)的近似式：

Tp+1i,j,k=Tpi,j,k+F0(Tpi+1,j,k+Tpi-1,j,k+Tpi,j+1,k+Tpi,j-1,k+Tpi,j,k+1+Tpi,j,k-1-6Tpi,j,k)

式中:F0=

λΔτ/[ρCp(Δx)2]。如果写成隐式差分格式，则除Tpi,j,k项外，其它项的Tp改为Tp+1。541.2.2能量平衡法能量平衡法能量平衡法是不再从微分方程出发用差商代替微商的办法建立差分方程。而是将导热的基本定律直接近似，局部地应用于围绕每个节点的各单元控制体积，用热量平衡法建立差分方程。由于它的基本思想是把能量守恒原则应用到每个单元体，因此常称为能量守恒原则的有限差分法。建立内节点差分方程①一维系统高宽无限(即高宽方向上无热流)，一定厚度的平板，T=f(x,τ)，初始时刻τ=0，T=T0，为了简化，考虑无内热源，λ、ρ、Cp均为常数。沿板厚方向（即热流方向）把均质物体分割为若干单元体，各单元体的端面积为F，几何步长为Δx，时间步长为Δτ（用P上角标表示第几个时间步长）。单元体中心为节点，节点温度代表该单元体的温度。二相邻节点间的导热面积是F。如图所示。图一维均质物体分割55建立内节点差分方程考虑内节点i，从时间等于τ开始的Δτ时间间隔内，通过截面F从(i-1)单元体流入到i单元体的热量为：

Tpi-1–Tpi

Qi-1→i=λF(————

Δx从(i+1)单元体流入到i单元体的热量为：

Tpi+1–Tpi

Qi+1→i=λF(————)

Δx由于热量流入i单元体，引起i单元体的内能变化。在Δτ时间内，i单元体的内能变化为：

ΔUTp+1i–Tpi

——=ρCp(FΔx)————

ΔΔ

因为能量平衡或守恒，所以∑Q=ΔU/Δτ，即：

Tp+1i–TpiTpi-1–TpiTpi+1–Tpi

ρCp(FΔx)————=λF(————)+λF(————

Δ

ΔxΔx整理可得：

Tp+1i=Tpi+F0(Tpi+1-2Tpi+Tpi-1)

式中:F0=

λΔτ/[ρCp(Δx)2]。上式和直接代换法式完全一样。56建立内节点差分方程②二维系统如果物体的形状使之只在X、Y方向上有热流，Z方向上无热流，或相对于X、Y方向可忽略不计，这种情况就可以当作二维系统处理。对于如图所示的二维矩形区域，划分成有限个小单元体。图中虚线所示的(i,j)单元体及周围四个相邻的单元体，用于推导内单元体差分方程。每个小单元体的几何步长为Δx和Δy；时间步长Δτ。为使问题简化，仍然假定热物性值为常数，考虑无内热源。图二维系统的分割57建立内节点差分方程对于(i,j)单元体，Δτ时间间隔内的内能变化为：

ΔUTp+1i,j–Tpi,j

——=ρCp(ΔxΔy.1)—————

ΔΔ在Δτ时间内从周围四个相邻单元体流入(i,j)单元体的热量分别为：

Tpi-1,j–Tpi,j

Qi-1,j→i,j=λ(Δy.1)(——————

ΔxTpi+1,j–Tpi,j

Qi+1,j→i,j=λ(Δy.1)(——————

ΔxTpi,j-1–Tpi,j

Qi,j-1→i,j=λ(Δx.1)(——————

ΔyTpi,j+1–Tpi,j

Qi,j+1→i,j=λ(Δx.1)(——————

Δy因为∑Q=ΔU/Δτ，若Δx=Δy，经整理后得：

Tp+1i,j=Tpi,j+F0(Tpi+1,j+Tpi-1,j+Tpi,j+1+

Tpi,j-1-4Tpi,j

式中:F0=

λΔτ/[ρCp(Δx)2]。上式和直接代换法式完全一样。③三维系统同理可得和式(1.2.28)一样的节点差分方程。58能量平衡法与直接代换法的比较能量平衡法与直接代换法的比较尽管用能量平衡法与直接代换法所得到的差分方程在形式上完全一样，但它们所依赖的基础有所不同。用直接代换法得出差分方程时，要求区域内温度分布连续而光滑，因为按泰勒级数展开，要求温度T(x)在x处的各阶导数都存在，连续而光滑即可导。而能量平衡法只要求每个单元体边界上热流可积，即可以算出每个单元体边界上的热量，它比温度分布连续光滑的要求显然放宽了。如果区域由导热系数截然不同的两种物质构成，并在交界面上存在接触热阻，或在交界面上存在间隙，那么在那些交界面上或间隙处，温度分布就不连续或不光滑。在这种情况下，可用能量平衡法建立差分方程，只需在单元划分时，使那些交界面正好落在单元的边界即可，如图所示。图异种物质接触时单元的划分图界面温度分布示意图59能量平衡法与直接代换法的比较现假定二维矩形分割，交界面两侧的物体导热系数不同，讨论相邻两节点(i+1,j)节点流出(i,j)节点的热量公式的建立。当异种物质接触时，除传导传热外，还要考虑界面热阻Rb，则

Tpi+1,j–Tp2

Qi+1,j→T2=λi+1,j(Δy.1)(—————

Δx/2Tp2–Tp1

QT2→T1=(Δy.1)(————Rb

Tp1–Tpi,j

QT1→i,j=λi,j(Δy.1)(————

Δx/2因为三个Q相似，令其相等，由式和式1.2.42得：

QΔx/2

———

———

=Tpi+1,j–Tp2

Δy.1λi+1,j

QΔx/2

———

———

=Tp1-

Tpi,j

Δy.1λi,j

QΔx/2Δx/2

———(———+———)

=(Tpi+1,j–Tpi,j)-(Tp2-Tp1)

Δy.1λi+1,jλi,j

将式代入上式得：

60能量平衡法与直接代换法的比较

Δy.1

Q=—————————————(Tpi+1,j–Tpi,j)

Rb+Δx/(2λi+1,j)+Δx/(2λi,j)

与直接代换法相比，能量平衡法的主要优点是：在二维系统中能使用三角形、四角形，在三维系统中能使用四、五、六面体等多种分割单元，因此能解复杂形状的问题；对于每个节点单元能指定热物性值，易于求解由多种物质组成的系统的问题。能量平衡法的缺点是：输入数据量多，程序复杂，计算时间较长。所以，对于矩形或圆柱形等简单形状，能进行有规律分割的系统，采用直接代换法是有利的。可是，能有规律分割时，能量平衡法的程序与直接代换法的程序相同。一般说，能量平衡法更实用些。1.2.3边界节点差分方程的建立物体内部的温度场必然受到物体表面条件的影响，如表面散热条件差，表面温度高，物体内部的温度场平缓。反之，物体内部温度场的变化也影响着表面上的条件，如物体内部的导热系数大，物体内部温度场平缓，使表面温度升高，表面散热量也大。为了数值模拟，必须建立边界节点的差分方程。在此，先讨论一般换热条件的边界节点差分方程的建立。611.2.3边界节点差分方程的建立图所示了绝热、给定热流密度qr、对流、定温、辐射换热等五种边界表面的二维矩形区域网格，现用能量平衡法建立各种换热边界条件下的节点方程。图各种换热边界表面图绝热边界单元绝热边界把图中AB面上任一节点(i,j)及其相邻节点取出分析，如图所示。

Tpi+1,j–Tpi,jTpi,j-1–Tpi,j

λ(Δy.1)(—————)+λ(Δx/2.1)(—————)+

ΔxΔyTpi,j+1–Tpi,jTp+1i,j–Tpi,j

λ(Δx/2.1)(—————)+0=ρCp(Δx/2Δy.1)—————

ΔyΔ

62绝热边界若Δx=Δy，经整理后得：

Tp+1i,j=Tpi,j+F0(2Tpi+1,j+Tpi,j-1+Tpi,j+1-4Tpi,j

式中:F0=

λΔτ/[ρCp(Δx)2]。给定热流密度qr边界把图中BC面上任一节点(i,j)及其相邻节点取出分析，如图所示。用能量平衡法建立给定热流密度qr边界节点方程。

Tpi-1,j–Tpi,jλ(Δy/2.1)(—————)+

ΔxTpi+1,j–Tpi,jλ(Δy/2.1)(—————)+ΔxTpi,j+1–Tpi,j

λ(Δx.1)(—————)+qr(Δx.1)

ΔyTp+1i,j–Tpi,j=ρCp(ΔxΔy/2.1)—————Δ

若Δx=Δy，经整理后得：图给定热流密度边界单元