金属凝固计算机模拟b.ppt

1,金属凝固计算机模拟,前言温度场数值模拟浓度场数值模拟,北京科技大学材料学院吴春京,尼加拉瓜瀑布2000.05.18摄,2,前言,、材料加工过程铸造、锻造、焊接、挤、压、拉、拔、热处理。专业设置。（大学、研究生、专业学位）2、计算机在材料加工中应用领域有五个方面：（1）科学计算数值模拟优化计算如最低成本配料优化计算，其是应用线性规划解法，求满足化学成分要求范围内，在已知化学成分和价格的几种或几百种原料中，同时还可考虑某种原料最多，最少或必须用多少的条件，求出最低成本的配料单来。也可以用此方法求出最高成本的配料单。传统手算配料单，由于其不考虑成本，也无法考虑成本，因此，传统手算的配料成本落在最高和最低成本之间。,3,4,5,6,7,8,9,前言,（2）信息处理信息处理的特点是计算简单、数据量大。现在7585的计算机用于信息处理。数据库数据库（DATABASE）的概念由美国军方于1963年提出。生产、经营、人事、工资等管理都可应用数据库管理。材料加工的生产、经营、人事、工资等管理都可应用数据库管理。,10,前言,专家系统第一个专家系统是1965年由美国斯坦福大学开发的，用质谱仪得到的数据来决定一个未知化合物分子结构的程序。专家系统是将专家的知识和经验按一定规则存在知识库或数据库里，用户则可通过专家系统，查寻检索到存在专家系统内的专家的知识和经验。材料加工的计算机专家系统正在研究之中，如“铸造缺陷分析专家系统”，“辊型设计的专家系统”。,11,12,前言,斯坦福大学（StanfordUniversity）创建于1885年，位于加州圣克拉拉县。总面积超过8000英亩，校园约5700英亩，住宅约500英亩，剩下大半个校园，除了一些缓冲地段之外，就是植物园、高尔夫球场和若干个科学实验场。校园构建誉为美国城市规划和庭院设计的典范。紧凑、集中、和谐。主要建筑：土黄色石墙环绕下的红屋顶建筑，拱廊相接，棕榈成行。最有名的建筑是斯坦福纪念教堂。校园内公路长达46英里，三个水坝和湖。1英亩=6.070市亩，1英里=1.609公里,13,前言,斯坦福的工业园是闻名遐迩的“硅谷”诞生地，在64名硅谷最重要的先驱中，有28人是斯坦福大学的校友、教授或研究人员，占42%。7个学院（地球科学、工、文理、商、教育、法、医）专职教师1714人，从建校到现在已有25位诺贝尔奖得主。现有教师中，有17位诺贝尔奖得主，4位普利策奖得主，23位麦克阿瑟奖得主，21位国家科学奖得主，3位国家技术奖得主，220位美国文理科学院院士，128位国家科学院院士，83位国家工程院院士，25位国家教育科学院院士，42位美国哲学社会科学学会成员，6位Wolf基金奖获得者，6位Koret基金奖获得者，2位总统自由奖获得者。,14,前言,（3）计算机辅助设计和制造。计算机辅助设计，简称CAD。计算机辅助制造，简称CAM。近年来由于在计算机上进行设计，就连设计图的概念也有所改变，如在波音767飞机设计中，图纸定义为点列的集合，然后从美国把录有这些点列数据的磁盘寄到日本等国即可生产。输出设计图只是为了便于人的理解，而并非照图加工。机械加工时全部使用数控加工机械。材料加工工艺设计的CAD，材料加工工装制造的CAM，材料加工车间设计的CAD系统正在研究发展之中。由于材料加工的复杂性，材料加工车间的多样性和研究发展水平限制，现在材料加工的CADCAM只局限于某一类型的CADCAM。,15,前言,现代英国铸造的CADCAM技术发展到每周能够出图18000多幅。日本日立公司研究的设计金属压铸金属型CAD系统，能够在获得产品数据之后，启动画出型腔图，并对模型进行计算，输出用于CAM的数据磁盘。（4）计算机测试和控制系统计算机和检测仪表、控制部件结合，可形成计算机测试、控制系统。材料加工生产中的测试和控制问题，一般地说，都可以采用计算机测试、控制系统。（5）机械手和机器人,16,前言,、数值模拟两个不同领域的现象，能用同一数学形式来描述，称这两个现象彼此是可模拟。模拟的方法是把一个领域内求解的问题过渡到另一个领域中去解决。随着计算机和计算数学的发展，用计算机数值计算法求解问题的计算精度已经达到很接近解析解，此法称为计算机数值模拟。仿真的定义，与模拟的差异。、工艺优化计算机数值模拟的最终目的是解决工艺优化设计问题。一旦全面实现了材料加工过程的计算机数值模拟，材料的加工生产将会产生深刻的变革。,17,前言,、计算机数值模拟的步骤给定研究对象几何条件、物理条件、初始条件、边界条件离散化处理将过程所涉及的区域在空间和时间上进行离散化处理。建立数值方程建立内节点和边界节点的数值方程。选择计算方法选择适当的计算方法求解线性代数方程组。编制计算机程序编制计算机程序，由计算机算出结果，并用数据、曲线、图形输出。优化工艺分析计算机模拟的结果，如果结果不理想，调整工艺参数，再进行计算机模拟，直到模拟结果为最佳结果。,18,前言,、讲解的主要内容通过一、二个例子，让大家亲自经历计算机数值模拟的一般步骤，从中掌握其基本原理和方法，为更好应用计算机数值模拟，优化工艺参数打下基础。材料加工过程中，比较常用、典型的数值模拟有：温度场数值模拟浓度场数值模拟、教学安排和考核讲课18学时自己上机练习4学时考试（笔试开卷）2学时成绩评定：笔试开卷（70%）；交上机的程序及运行结果（30%）,19,20,1温度场计算机数值模拟,1.1传热的基本知识1.1.1传热的基本方式导热导热属于接触传热，是连续介质就地传递热量，没有各部部分物质之间宏观的相对位移。在不透明固体实体内部，由于各部分物质之间无法作宏观的相对位移，不透明无法传递辐射能，实体保证接触，所以只能依靠导热方式传递热量。导热的基本定律是傅立叶定律，即导热的热流密度q与温度梯度成正比，即qgradT=-gradT=-(W/m2)1.1.1n式中：q传热方向上单位时间、单位面积上所通过的热量，J/(s.m2)=W/m2；材料的导热系数，W/(m.K)；温度梯度，K/m。n负号表示导热的热流量向温度低的方向传递，即与温度梯度的方向相反。热流密度是个向量，即它有大小和方向。,21,1.1.1传热的基本方式,如果热流密度的分量和（X、Y、Z）坐标系相联系，那么X、Y、Z方向的热流密度分量应是：qx=-,qy=-,qz=-1.1.2XYZ热流密度q=iqx+jqy+kqz1.1.3导热系数物理意义：沿导热方向的单位长度上，温度减低,物质所容许通过的热流量。方向性：大多数液体和固体属于各向同性的物质。各向异性材料的导热系数具有方向性，如石墨。温度函数：值还随温度而变化。大多数金属的导热系数随温度的升高而降低。大多数液体（水和甘油除外）其导热系数随温度的升高而降低。气体的导热系数随温度的升高而增加。,22,1.1.1传热的基本方式,对流对流是流体（气体和液体）中温度不同的各部分相互混合的宏观运动引起热量传递的现象。工程上最具有实际意义的是：相对运动着的流体与所接触的固体壁面之间的热量交换过程，一般称为对流换热。工程上在研究固体壁面和流体之间的对流换热时，除了高度稀薄的气体外，人们不去注意流体的单个质点，而把流体看成是连续介质。实际的流体总有粘性，流动时，受粘性和壁面摩擦的影响，在靠近壁面附近的流体将降低流速，在壁面上完全被滞止不动，即X=0时，=0，如图1.1所示。因此，热量从壁面传给贴壁的那部分流体，将依靠导热。,T(K)V(m/s)TwqwcTfVX(m)图1.1邻近壁面的流体速度分布和温度分布,23,1.1.1传热的基本方式,流体和固体壁接触面上的“相”际热流密度为：qwc=-fx=+01.1.4X式中：qwc热流密度，W/m2；f流体的导热系数，W/(m.K)；T流体的温度，K。液体的导热系数较小，气体的导热系数更小，所以受热时，在贴壁处的流体温度势必沿着轴的反方向急剧升高。随着离壁面的距离的增加，流体的流动将带走更多的热量，使轴方向的温度梯度连续下降，温度分布趋向平坦化。正是通过这种导热和对流的共同作用，使热量在流体内部得到传播，越临近壁面，导热越起主导作用。图1.1所示，假如厚度为的贴壁静止膜，膜内温度线性地从壁面温度TW降到远离壁面，尚未被加热的流体温度Tf，则Tf-TWqwc=-f1.1.5无界对流时壁面与流体的换热，钢锭与周围空气的对流换热属于这种情况。,24,1.1.1传热的基本方式,流体在管和槽道内部的流动，称为有界对流，这时不存在远离壁面，尚未被加热的流体温度，则采用截面平均温度作为流体的特征温度Tf，则qWC=aC(TW-Tf),W/m21.1.6这就是所谓的“牛顿冷却定律”。式中：aC为放热系数，W/(m2.K)。其实“牛顿冷却定律”并不是表达对流换热现象本质的物理定律。凡能影响流体流动的各种因素，包括流体的种类和状态、运动的速度、流道的形状和大小，以及固体壁面粗糙度等，都会对aC值产生影响。式1.1.6只不过形式地把放热过程的一切复杂性和计算上的困难，都转移到并且集中在放热系数这样一个物理量上罢了。aC代表流体和所接触的固体表面之间温度每相差，该流体与表面之间“相”际热流量的大小。运用式1.1.6可以进行对流换热的计算。但由于对流换热的复杂性，该式中的放热系数aC需从相应的准则方程式求出。准则方程式是针对不同对流换热情况，在综合实验结果的基础之上，运用相似理论将表征某现象的物理量整理成一些相似准则，用因次分析法得到的各种类型的表达式。连铸钢坯二冷区,25,辐射只要温度高于绝对零度，任何物体都随时向周围空间辐射能量。辐射热流密度用斯蒂芬玻耳兹曼定律表达：E=，w/m21.1.7式中：物体的辐射率或黑度,(0-1)；斯蒂芬玻耳兹曼常数或绝对黑度的辐射常数；5.6710-8W/m2K4温度,K。实际上，辐射往往涉及二个物体间辐射热交换。如果二个物体的表面温度不同，中间被空气所隔开，T1T2时，则相互辐射的结果，表面温度为T1的物体发射出去的辐射热超过了吸收来自表面温度为T2的物体辐射热，引起辐射换热的热流量则为：Q12=12(1-2)F112或q12=12(1-2)121.1.8式中：12物体1与2综合黑度；F1物体1的表面积；12物体的表面向外发射出去的辐射热量中，能投射到物体表面上的份额，称为角系数,(0-1)。,1.1.1传热的基本方式,26,轧辊喷淬时的工作环境,27,喷淬机的简易俯视图,28,红外测温仪+吹扫装置1：红外测温仪；2：固定支架；3：瞄准管；4：进气管：5：瞄准管前端,29,辊身不同深度处的冷却曲线,30,喷雾过程辊身处表面综合换热系数与时间的关系,31,1.1.1传热的基本方式,壁面在气体中冷却，存在辐射换热和对流换热。考虑到壁面与气体之间存在着辐射换热，其热流量（或热流密度）为Qr=arF(Tw-Tf),w或qwr=ar(Tw-Tf),w/m21.1.9式中：qwr单位时间、单位面积辐射的热量（辐射热流密度），w/m2；ar辐射放热系数，w/(m2.k)；Tw辐射物体表面温度，k；Tf透明的气体介质的特征温度，k。考虑到壁面与气体之间还存在着对流换热，其热流密度为qwc=ac(Tw-Tf),w/m21.1.10由壁面传走的总热流密度qw应是qwr和qwc二者之和qw=ac(Tw-Tf)+ar(Tw-Tf)1.1.11令a0=ac+ar则qw=a0(Tw-Tf),w/m21.1.12用辐射放热系数ar，可以形式地把辐射换热折合成对流换热，用总放热系,32,1.1.1传热的基本方式,数a0兼顾辐射换热的影响，从而有利于简化复杂传热的分析和计算。如远离表面的外界表面温度趋于环境温度Tf，并且12=1时，由式1.1.8得qwr=12(w-f)1.1.13由式1.1.9得：qwr=ar(w-f)1.1.14由式1.1.13和式1.1.14得：12(w-f)=ar(w-f)1.1.15因w-f=(w-f)(w3+w2f+wf2+f3)设m=1/2(w+f)，(w3+w2f+wf2+f3)4m3ar412m3,w/(m2.k)1.1.16从此可见，(w-f)降为零时，ar并非零值，而是以412m3。随着温差的扩大和平均温度m的升高，ar值将迅速增加。由于ac随温差的变化较小，在高温范围和大温差情况下，ar有可能成为a0的主要组成部分。ar与气体的运动状况无关，而ac与气体速度的降低而减小，在气体自然对流的情况下，ac5,w/(m2.k)，即使m只有300K，ar就可能占a0的一半左右。,33,1.1.2导热的偏微分方程式,钢锭一般属于各向同性的物质，其加热或冷却过程数学模拟计算依据的基本数学模型，是不稳定导热偏微分方程。下面讨论各向同性材料导热方程式的建立。直角坐标系导热偏微分方程导热偏微分方程的建立，是通过考察处于导热过程中的物质的微元体积（xyz)的能量平衡来建立。如图1.2所示。在时间内，通过六个面的导热所获得的能量，加上微元体内产生的内热源热量，要等于微元体积内物质积蓄热量的改变，即温度升高或降低。,x,z,y,dQx+x,dQx,dQz,dQy,dQz+z,dQy+y,图1.2直角坐标系导热方程式的微元体,34,1.1.2导热的偏微分方程式,图1.2中，微元控制体尺寸x、y、z，按照傅立叶导热定律，在X方向流入微元体左表面的热流可表示为：TdQx=-(yz)，W1.1.17X式中：导热系数，W/(m.k)；yzX方向微元体表面积，m2；TX方向的温度梯度，k/m。X在X方向流出微元体右表面的热流可表示为：TdQx+x=-(yz)，WX在X方向导热的净热流为：?dQx-dQx+=0,35,1.1.2导热的偏微分方程式,图1.2中，微元控制体尺寸x、y、z，按照傅立叶导热定律，在X方向流入微元体左表面的热流可表示为：TdQx=-(yz)，W1.1.17X式中：导热系数，W/(m.k)；yzX方向微元体表面积，m2；TX方向的温度梯度，k/m。X在X方向流出微元体右表面的热流，可以应用泰勒级数展开，保留级数的第一项和第二项而得到：TdQx+x=dQx+(dQx)x=dQx+-(yz)xXXXT=dQx-()xyz1.1.18XX在X方向导热的净热流为TdQx-dQx+=()xyz1.1.19XX,36,1.1.2导热的偏微分方程式,用同样的方法，可以得出Y,Z方向与式1.19相类似的导热的净热流方程式，即TdQy-dQy+y=()xyz1.1.20YYTdQz-dQz+z=()xyz1.1.21ZZ在三个坐标方向净热流的总和为：TTT()+()+()xyz1.1.22XXYYZZ如果单位时间、单位空间所产生的热量为Q(x,y,z,)，那么微元体单位时间的发热量（热流）为：Q(x,y,z,)xyz1.1.23由于导热传进微元体的净热流式1.1.22和微元体内产生的热量式1.1.23一起用于增大微元体的内能。微元体的内能的增加表现在微元体能量存储随时间的变化率，即TCpxyz1.1.24,37,1.1.2导热的偏微分方程式,式中：密度，kg/m3；Cp比热，J/(kg.k)；时间，s。对微元体进行能量平衡，使能量存储的时间变化率与由导热引起的流入微元体的净热流和微元体内产生的热流之和相等，得下式：TTTTCp=()+()+()+Q1.1.25XXYYZZ圆柱坐标系导热偏微分方程实际问题常常涉及柱面对称问题，且边界条件给定在一个表面上，因此表面具有一个坐标保持不变的性质。在这种情况下，采用圆柱坐标系是合适的。,图1.3圆柱坐标系导热方程式的微元体,38,1.1.2导热的偏微分方程式,图1.3所示圆柱坐标系，直接按内外两个圆弧面和其它四个平面组成的微元体积，在导热过程中热量必须按收支平衡的原则导出，此时，微元体积为：dv=(rd)dzdr1.1.26沿内圆弧面流入微元体积的热流：TdQr=-(rddz)1.1.27r沿外圆弧面流出微元体积的热流：TdQr+dr=dQr+(dQr)dr=dQr+-(rddz)drrrrT1T=dQr+(-r)ddzdr=dQr-(r)dv1.1.28rrrrr沿平面流入微元体积的热流：TdQ=-(drdz)1.1.29r沿+d平面流出微元体积的热流：TdQ+d=dQ+(dQ)rd=dQ+-(drdz)rdrrr,39,1.1.2导热的偏微分方程式,1=dQ-()dv1.1.30r2沿z面流入微元体积的热流：TdQz=-(rddr)1.1.31z沿z+dz面流出微元体积的热流：TdQz+dz=dQz+(dQz)dz=dQz+-(rddr)dzzzzT=dQz-()dv1.1.32zz根据直角坐标系导热微分方程推导的思路，可得到：1T1TT(r)+()+()+Q=Cp1.1.33rrrr2zz,40,1.1.2导热的偏微分方程式,球坐标系导热偏微分方程对于球坐标系，如图1.4所示。,图1.4球坐标系导热方程式的微元体,由内、外两个球面、两个圆弧面和两个平面所组成的微元体积为：dv=(rd)（rSind)dr1.1.34,41,1.1.2导热的偏微分方程式,沿半径方向流入微元体积的热流：TdQr=-(rSind.rd)1.1.35r沿半径方向流出微元体积的热流：TdQr+dr=dQr+(dQr)dr=dQr+-(rSind.rd)drrrrT1T=dQr+(-r2)Sind.ddr=dQr-(r2)dv1.1.36rrr2rr沿方向流入微元体积的热流：TdQ=-(rSind.dr)1.1.37r沿方向流出微元体积的热流：TdQ+d=dQ+(dQ)rd=dQ+-(rSinddr)rdrrrTT=dQ-(Sin)rddr.rd=dQ-(Sin)dvr2Sinr21.1.38,42,1.1.2导热的偏微分方程式,沿方向流入微元体积的热流：TdQ=-(dr.rd)1.1.39(rSin)沿方向流出微元体积的热流：dQ+d=dQ+(dQ).(rSin)d(rSin)T=dQ+(-dr.rd)(rSin)d(rSin)(rSin)T=dQ-()dv1.1.40(r2Sin2)同理可整理得：TTT(r2)+(Sin)+()+Q=Cpr2rrr2Sinr2Sin1.1.41,43,1.2导热方程的有限差分解法,求解不稳定导热偏微分方程的数值解法，主要有：有限差分解法、有限元法、边界元法三类。边界元法正在研究和完善之中，目前常用的是有限差分解法和有限元法。我们这门课专门讨论有限差分解法的数学基础，数值方程的建立，差分方程的稳定性和收敛性等问题。有限差分解法是用差分方程近似地代替微分方程，建立差分方程有直接法和能量平衡法两种。1.2.1直接代换法直接代换法是从微分形式出发，用差商直接代换微商的办法建立差分方程。1.2.1.1数学基础微商和差商的定义若T(x)是连续函数，则其导数为：dT(x+x)-T(x)=lim=lim1.2.1dxx0xx0x,44,1.2.1.1数学基础,式1.2.1右边/x是有限的差商，x和都不为零。式左边d/dx是T/x当x0时的差商，称之为微商。在x没有达到零之前，T/x只是dT/dx的近似。实际应用x0。如果把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)ndnTT(x+x)=T(x)+x+1.2.2dx2!dx2n!dxn整理后得：TT(x+x)-T(x)dTxd2T(x)n-1dnT=+1.2.3xxdx2!dx2n!dxn从上式可知，T(x)在x处的差商T/x等于函数T(x)在x处的各阶导数的线性组合，只能是近似地等于差商。两者之间也必然有偏差。,45,1.2.1.1数学基础,图1.2.1表示了一阶差商与一阶微商之间的关系。用不同方法得到的差商去代替微商，它们带来的误差是不同的。即向前差商：dT/dxT(x+x)-T(x)/x1.2.4向后差商：dT/dxT(x)-T(x-x)/x1.2.5中心差商：dT/dxT(x+x)-T(x-x)/2x1.2.6,图1.2.1差商与微商,46,1.2.1.1数学基础,按照泰勒级数展开，T(x)与T(x+x)的关系如下式所示：dT(x)2d2T(x)ndnTT(x+x)=T(x)+x+1.2.7dx2!dx2n!dxn整理后得：T(x+x)-T(x)dTxd2T-=+=O(x)1.2.8xdx2!dx2即向前差商的偏差是截去了泰勒级数展开式中的高阶项而引起的，常称“截断误差”，其截断误差为与x同级的小量O(x)。同理dT(x)2d2T(x)3d3TT(x-x)=T(x)-x+-+1.2.9dx2!dx23!dx3整理后得：T(x)-T(x-x)dTxd2T-=-+=O(x)1.2.10xdx2!dx2即向后差商的截断误差为与x同级的小量O(x)。,47,1.2.1.1数学基础,由式1.2.7减式1.2.9，将2xdT/dx移至等式左边，两边再除以2x，得：T(x+x)-T(x-x)dT(x)2d3T-=+=O(x)21.2.112xdx3!dx3即中心差商的截断误差为与(x)2同级的小量O(x)2。当x固定时，上述一阶差商一般仍是x的函数，对它们还可以求差商。这种一阶差商的差商称为二阶差商，它是二阶微商的近似，常用向前和向后差商来二阶微商，即d2TT(x+x)-T(x)T(x)-T(x-x)-/xdx2xxT(x+x)-2T(x)+T(x-x)=1.2.12(x)2由式1.2.7和式1.2.9相加，经简化后得：T(x+x)-2T(x)+T(x-x)d2T-=O(x)21.2.13(x)2dx2即二阶差商的截断误差为与(x)2同级的小量O(x)2。,48,1.2.1.2建立内节点差分方程,一维系统假定有一高宽无限(即高宽方向上无热流)，厚度为L的平板，T=f(x,)即温度是x方向位置和时间的函数，所谓一维系统是指几何空间为一维。初始时刻=0，T=T0，为了简化，考虑无内热源，、Cp均为常数。选取网格点间距x和时间步长将研究对象离散化。显式差分格式TT一维不稳定导热方程为Cp=()XX该方程在区域0，0xL内全部点都成立，如图1.2.2所示。将方程应用于内节点i可写成：Tp2TpCp()=()p?1.2.13iX2i上式等号左端的微分式用温度对时间的一阶向前差商来近似，即：TpTp+1i-Tpi()=+O()1.2.14i上式等号右端的微分式用温度对空间的二阶差商来近似，即：TpTpi+1-2Tpi+Tpi-1()=+O(x)2)1.2.15x2i(x)2,49,1.2.1.2建立内节点差分方程,图1.2.2一维显式差分将式1.2.14和式1.2.15代入式1.2.13得：Tp+1iTpiTpi+1-2Tpi+Tpi-1Cp+O()=+O(x)2)1.2.16(x)2若在上式中去掉O()和O(x)2),整理得：Tp+1i=Tpi+F0(Tpi+1-2Tpi+Tpi-1)1.2.17式中:F0=/Cp(x)2=导热速率/蓄热速率，称F0为傅立叶数。,50,1.2.1.2建立内节点差分方程,F0的数值小意味着加热或冷却此物体所需要的时间长，反之，所需要的时间短，F0是一个无因次数字。当初始条件和边界条件已知时，用式1.2.17就可模拟区域内各节点随时间增长的温度值Tpi(i=2,3,，N-1；p=1,2,3,）隐式差分格式一维隐式差分如图1.2.3所示，将一维不稳定导热微分方程应用于内节点i，则：Tp2Tp+1Cp()=()1.2.18iX2i式1.2.18和1.2.13相比，式的左边完全一样，温度对时间的一阶偏微商，仍用一阶向前差商来近似，而式1.2.18和1.2.13右边有所不同。式1.2.13中温度对距离的二阶偏微商是对应时刻p的，在用差商近似微商时，用p时刻的T值；而式1.2.18中，温度对距离的二阶偏微商是对应时刻p+1的，用差商近似微商时，用p+1时刻的T值。即式1.2.18相应的差分方程为：Tp+1iTpiTp+1i+1-2Tp+1i+Tp+1i-1Cp+O()=+O(x)2)1.2.19(x)2若在上式中去掉O()和O(x)2),整理得：Tp+1i=Tpi+F0(Tp+1i+1-2Tp+1i+Tp+1i-1)1.2.20,51,1.2.1.2建立内节点差分方程,式中:F0=/Cp(x)2。比较式1.2.17和1.2.20，可以看出显式差分格式的突出优点是每个节点方程都可以独立求解，整个计算过程十分简便。但是，的选取要受到限制，有时为了满足差分格式稳定性条件，可能选得很小，使计算工作量加大。隐式差分格式的最大优点是，它对任意F0值都是稳定的。这种稳定是绝对的，即不受边界条件、x、热物参数的影响，于是可以选的较大，计算速度加快。但是，对于节点i，只从式1.2.20不能独立求解，必须涉及Tp+1i+1、Tp+1i、Tp+1i-1的联立线性代数方程组才能求解，也就是说，它含有三个未知数。时间步长每前进一步，从坐标0xL，网格点1iN整个区域的每个点，上述方程都要列出一次（见图1.2.3)。因此，向一个新的时间步长每移动一步就必须解一个方程组。当按顺序列出这些方程时，除了要的第一点和最后一点的方程只有两个未知数外，其余每一个方程都含有三个未知数，于是方程是三对角的。对于这种情况，已有适用的求解程序，后面将讨论。关于差分方程稳定性将在后面讨论。,52,1.2.1.2建立内节点差分方程,图1.2.3一维隐式差分二维系统显式和隐式差分格式建立的方法和两种差分格式的特点在前面讨论过，对二维系统同样适用。为简略起见，在此只讨论二维系统显式差分方程的建立。为使问题简化，仍然假定热物性值为常数，考虑无内热源。首先把所讨论的区域离散化，如图1.2.4所示。,53,1.2.1.2建立内节点差分方程,网线：用平行于X、Y轴的直线（网线），进行空间离散化节点：网线与网线的交点步长：节点间的距离图1.2.4二维均质物体的分割x与y可以是不变的常量，即等步长，也可以是变量（即在区域内的不同处是不同的），即变步长。如果区域内各点处的温度梯度相差很大，则在温度变化剧烈处，网格布得密些，在温度变化不剧烈处，网格布得疏些。至于网格多少，步长取多少为宜，要根据计算精度与计算工作量等因素而定。对于0xL1和0Tp+1i+1；Tp+1iTp+1i-1,式1.2.75等号右边括号内三项的代数和必然大于零，从而必然有TpiTp+1i。也就是说，在p时刻，区域内的最大温度必然大于p+1时刻区域内的最大温度。依此类推，必然将最大温度值推到初始条件。假定p+1时刻，在区域内某一节点i处取得区域内最小温度值，即Tp+1i为这一时刻区域内的最小温度,Tp+1iTp+1i+1；Tp+1iTp+1i-1,式1.2.75等号右边括号内三项的代数和必然小于零，从而必然有TpiTp+1i。按照上面的分析方法可知，整个过程的温度最小值，必然出现在初始条件中。总之式1.2.75不管F0为何值，它的运算逻辑都是符合热力学原理的。即隐式差分格式是无条件稳定的。由上作出结论：隐式差分格式是无条件稳定的，显式差分格式是无条件稳定的。这一结论，原则上对二维、三维系统也是适用的。显式差分格式的稳定性条件是要求Tpi项的系数0。由式1.2.17（一维）、式1.2.25（二维）、式1.2.28（三维）和式1.2.52（二维对流边界）等差分方程可导出具体的稳定性条件如下：内部节点：,77,1.2.5差分方程的稳定性和收敛性,一维直角坐标：F01/2二维直角坐标：F01/4三维直角坐标：F01/6边界节点（设为对流边界）11一维直角坐标：F0二维直角坐标：F02(1+acx/)2(2+acx/)1三维直角坐标：F02(3+acx/)分析ac=0得到内部节点判据；当其它条件不变时，ac增大，收敛条件要求更小，即更严格。根据这一原则，其它边界条件下边界节点差分方程的稳定性也可以相应导出。对于圆柱坐标系、球坐标系的显式差分方程，也要按上述原则分别导出具体的稳定性条件。总之，用显式差分格式数值求解导热问题，离散化处理时，时间步长和空间步长x的选取是受到稳定性条件制约的。一般的习惯是先定空间步长，再计算时间步长。,78,1.2.6二维交替隐式格式,显式差分格式的优点：每个节点方程能独立求解，整个计算过程十分简单。缺点：为满足稳定性条件，空间步长和时间步长的选取必须遵守下列关系：(x)2Cp一维系统：(F01/2)2(x)2Cp二维系统：(F01/4)4(x)2Cp三维系统：(F01/6)6例如，一维系统，x=1mm，=7840kg/m3，Cp=610J/(kg.)，=30J/(s.m.)，求得：0.079,s。即1h约5万次。这就表明：为了满足精度要求，当网格划分较细时，时间步长大大减小，当维数愈高，在同样的空间步长条件下，时间步长也愈小。因此，为了得到稳定而又比较精确的答案，就要进行大量的计算。为了消除稳定性时间步长的限制，可以采用隐式格式求解。但是，在每,79,1.2.6二维交替隐式格式,一个时间层上，对二维系统就要解一个五对角的N*M个未知量的代数方程组(N为x方向的网格数，M为y方向的网格数)，即移动一个时间步长，就必须解整个方程组，每一个方程都有五个未知数和一个对角占优元素。这样一来，就没有解三对角方程组那样方便。因此，需要建立新的差分格式，希望这一差分格式能满足如下要求：无条件稳定；有合理精度的解；所产生的代数方程组易于求解。交替方向隐式方法正是满足这些要求的一种差分格式。（简称IAD）1.2.6.1差分方程的建立交替方向隐式方法需要对给定时间步长推导出二组有限差分方程。这些方程是显式和隐式项的混合，对于前半个时间步长，x方向各项是隐式形式，y方向各项是显式形式；而在下半个时间步长时则相反。若把二维导热微分方程式应用于节点(i,j)，可以写成如下形式：第一个1/2Tp2Tp+1/22TpCp()=()+()1.2.76i,jx2i,jy2i,jTp+1/2i,jTpi,jTp+1/2i+1,j-2Tp+1/2i,j+Tp+1/2i-1,jCp()=+/2(x)2Tpi,j+1-2Tpi,j+Tpi,j-11.2.77(y)2,80,1.2.6二维交替隐式格式,第二个1/2Tp+1/22Tp+1/22Tp+1Cp()=()+()1.2.78i,jx2i,jy2i,jTp+1i,jTp+1/2i,jTp+1/2i+1,j-2Tp+1/2i,j+Tp+1/2i-1,jCp()=+/2(x)2Tp+1i,j+1-2Tp+1i,j+Tp+1i,j-11.2.79(y)2从上可见：第p+1/2层是一个过渡时间层，通过式1.2.77计算过渡时间层的值，然后用式1.2.79计算第p+1时刻的值。由于现在有一种迅速解三对角线矩阵的算法，所以交替方向隐式法是有效的。1.2.7三对角线矩阵的解法前面介绍了显式、隐式和交替方向隐式法等差分格式，从数学方面来看，它们都是一个线性代数方程组。线性代数方程组的求解一般有两类方法：一类是直接方法，如三对角线矩阵的解法、高斯消去法、矩阵分解法；另一类是迭代方法，如高斯塞德尔迭代、超松弛迭代、强隐式迭代。一般说来，直接方法只要进行有限此运算就可以得到方程的解。它的缺点是必须有较大的计算机存储量和花费较多的计算时间，而且在使用时，,81,1.2.7三对角线矩阵的解法,还必须注意抑制舍入误差的影响，可以采用增加计算数字的位数，或利用主元消去等手段。迭代法没有舍入误差积累的问题，一般需要较少的计算机存储量，计算时间也比较节省。但是，如何能确保迭代收敛和收敛得比较快，这在使用中必须注意解决。当离散化方程是三节点关系式时，高斯消去法就转换成人们所常用的三对角线矩阵法，这种方法也称为托马斯(Thomas)算法或TDMA、追赶法。它之所以称为三对角线矩阵算法，是由于用矩阵形式写出这些三节点关系式时，所有系数都是沿矩阵的三对角线排列的。设离散化方程是如下形式的三节点关系式：-aiTi-1+biTi-ciTi+1=di1.2.80i=1,2,N，其中i=1和i=N分别为边界节点。由于T0和TN+1没有定义，故取a1=0和CN=0。因此，当i=1时，式1.2.80变成T1和T2之间的关系式，即可以用T2表示T1，T1=f1(T2)。（一个方程二个未知数）当i=2时，式1.2.80变成T1、T2和T3之间的关系式，可以用T2表示T3，T2=f2(T3)。（以上二个方程三个未知数）,82,1.2.7三对角线矩阵的解法,这种替换关系，可以一直进行到TN用TN+1来表示，即TN=fN(TN+1)，由于CN=0(TN+1的系数)，于是可以求得TN的值。（以上增加一个方程，没有增加未知数，所以可求出TN）而TN-1可从TN求得，TN-2可从TN-1求得，T可从T3求得，T1可从T2求得，这就是三对角线矩阵解法的实质。具体公式推导如下：设待求的递推替换关系式为：Ti=PiTi+1+Qi1.2.81则Ti-1=Pi-1Ti+Qi-11.2.82将式1.2.82代入式1.2.80，得-ai(Pi-1Ti+Qi-1)+biTi-ciTi+1=di，经整理得：cidi+aiQi-1Ti=Ti+1+1.2.83bi-aiPi-1bi-aiPi-1比较式1.2.83和式1.2.81得cidi+aiQi-1Pi=，Qi=1.2.84bi-aiPi-1bi-aiPi-1式1.2.84是个递推关系式，因为它们是用Pi-1，和Qi-1来表示Pi，和Qi的。,83,1.2.7三对角线矩阵的解法,根据以上分析，其求解步骤是：取i=1，因为a1=0，所以P1=c1/b1，Q1=d1/b1；对于i=2,3,N-1，用式1.2.84分别求出Pi和Qi；当i=N时，因为CN=0，所以PN=0，由式1.2.81得TN=QN，即求出点的温度；对于i=N-1,N-2,2,1，用式1.2.81分别求出TN-1，T2，T1。1.2.7.1三对角线矩阵解法程序TurboBASIC语言设A(I)为ai系数，B(I)为bi系数，C(I)为Ci系数，D(I)为di系数，N为N个温度变量，TMID(I)为Ti温度，M为计算M个时间步长。REM主程序N=11:M=12:S=0.35:lamda=30!:cp=610!:rou=7840!:dt=120!dx=S/(N-1):f0=lamda*dt/(rou*cp*dx*dx)DIMA(1:N),B(1:N),C(1:N),D(1:N),TMID(1:N)A(1)=0:B(1)=1:C(1)=0:D(1)=100!FORI=2TON-1A(I)=f0:B(I)=1+2*f0:C(I)=f0:D(I)=15!NEXTI,84,1.2.7三对角线矩阵的解法,A(N)=0:B(N)=1:C(N)=0:D(N)=100!CLS:PRINTdx=;dx*1000;(mm)PRINTdt=;dt/60;minPRINTTemperatureResultsPRINTmin1P2P3P4P5P6P7P8P9P10P11PFORJ=1TOMCALLTDMAPRINTJ*dt/60;FORI=1TONPRINTINT(TMID(I)*100)/100;NEXTIPRINTFORI=2TON-1D(I)=TMID(I)NEXTINEXTJEND,85,1.2.7三对角线矩阵的解法,REM解三对角线矩阵子程序SUBTDMASHAREDTMID(),A(),B(),C(),D(),NLOCALBAP(),DAQ()DIMBAP(N),DAQ(N)BAP(1)=B(1):DAQ(1)=D(1)FORI=1TON-1BAP(I+1)=B(I+1)-A(I+1)*C(I)/BAP(I)DAQ(I+1)=D(I+1)+A(I+1)*DAQ(I)/BAP(I)NEXTITMID(N)=DAQ(N)/BAP(N)FORI=N-1TO1STEP1TMID(I)=(DAQ(I)+C(I)*TMID(I+1)/BAP(I)NEXTIENDSUB,86,1.2.7三对角线矩阵的解法,QuickBASIC语言DECLARESUBTDMA()REMsimulN=11:M=19:s=.35:lamda=30!:cp=610!:rou=7840!:dt=120!dx=s/(N-1):f0=lamda*dt/(rou*cp*dx*dx)DIMA(1TON),B(1TON),C(1TON),D(1TON),TMID(1TON)A(1)=0:B(1)=1:C(1)=0:D(1)=100!FORi=2TON-1A(i)=f0:B(i)=1+2*f0:C(i)=f0:D(i)=15!NEXTiA(N)=0:B(N)=1:C(N)=0:D(N)=100!CLSPRINTdx=;dx*1000;(mm)PRINTdt=;dt/60;minPRINTTemperatureResultsPRINTmin1P2P3P4P5P6P7P8P9P10P11P,87,1.2.7三对角线矩阵的解法,FORj=1TOMCALLTDMAPRINTj*dt/60;FORi=1TONPRINTINT(TMID(i)*100)/100;NEXTiPRINTFORi=2TON-1D(i)=TMID(i)NEXTiNEXTjENDSUBTDMASHAREDTMID(),A(),B(),C(),D(),NDIMBAP(1TON),DAQ(1TON),88,1.2.7三对角线矩阵的解法,BAP(1)=B(1):DAQ(1)=D(1)FORi=1TON-1BAP(i+1)=B(i+1)-A(i+1)*C(i)/BAP(i)DAQ(i+1)=D(i+1)+A(i+1)*DAQ(i)/BAP(i)NEXTiTMID(N)=DAQ(N)/BAP(N)FORi=N-1TO1STEP-1TMID(i)=(DAQ(i)+C(i)*TMID(i+1)/BAP(i)NEXTiENDSUB程序解释：BAP(I)=B(I)-A(I)*P(I-1)，其是式1.2.84中bi-aiPi-1DAQ(I)=D(I)+A(I)*Q(I-1)，其是式1.2.84中di+aiQi-1BAP(1)=B(1)-A(1)*P(0)=B(1)，因为P0=0DAQ(1)=D(1)+A(1)*Q(0)=D(1)，因为Q0=0BAP(2)=B(2)-A(2)*C(1)/BAP(1)=B(2)-A(2)*P(1)，推广到I=1TON-1,89,1.2.7三对角线矩阵的解法,DAQ(2)=D(2)+A(2)*DAQ(1)/BAP(1)=D(2)+A(2)*Q(1)，推广到I=1TON-1D(N)+A(N)*DAQ(N-1)/BAP(N-1)TMID(N)=DAQ(N)/BAP(