《钢铁冶金过程计算机仿真.ppt》由会员分享，可在线阅读，更多相关《钢铁冶金过程计算机仿真.ppt（100页珍藏版）》请在三一办公上搜索。
1、1,钢铁冶金过程计算机仿真,温度场数值模拟,北京科技大学冶金学院刘青,2,前言,、钢铁冶金加工过程炼铁炼钢（铸造、锻造、焊接、挤、压、拉、拔、热处理）。专业设置。、数值模拟两个不同领域的现象，能用同一数学形式来描述，称这两个现象彼此是可模拟。模拟的方法是把一个领域内求解的问题过渡到另一个领域中去解决。随着计算机和计算数学的发展，用计算机数值计算法求解问题的计算精度已经达到很接近解析解，此法称为计算机数值模拟。仿真的定义，与模拟的差异。、工艺优化 计算机数值模拟的最终目的是解决工艺优化设计问题。一旦全面实现了冶金加工过程的计算机数值模拟，材料的加工生产将会产生深刻的变革。CAD/CAM/CAE的
2、发展，波音767飞机的设计和加工。,3,前言,、计算机数值模拟的步骤给定研究对象几何条件、物理条件、初始条件、边界条件离散化处理将过程所涉及的区域在空间和时间上进行离散化处理。建立数值方程建立内节点和边界节点的数值方程。选择计算方法选择适当的计算方法求解线性代数方程组。编制计算机程序编制计算机程序，由计算机算出结果，并用数据、曲线、图形输出。优化工艺分析计算机模拟的结果，如果结果不理想，调整工艺参数，再进行计算机模拟，直到模拟结果为最佳结果。,4,前言,、讲解的主要内容通过一、二个例子，让大家亲自经历计算机数值模拟的一般步骤，从中掌握其基本原理和方法，为更好应用计算机数值模拟，优化工艺参数打下
3、基础。冶金过程中，比较常用、典型的数值模拟有：温度场数值模拟浓度场数值模拟组织与性能模拟,5,6,1 温度场计算机数值模拟,1.1 传热的基本知识1.1.1 传热的基本方式导热导热属于接触传热，是连续介质就地传递热量，没有各部部分物质之间宏观的相对位移。在不透明固体实体内部，由于各部分物质之间无法作宏观的相对位移，不透明无法传递辐射能，实体保证接触，所以只能依靠导热方式传递热量。导热的基本定律是傅立叶定律，即导热的比热流量q与温度梯度成正比，即：qgrad T=-grad T=-(W/m2)（1-1-1）n式中：q 传热方向上单位时间、单位面积上所通过的热量，J/(s.m2)=W/m2 材料的
4、导热系数，W/(m.K)温度梯度，K/m n负号表示导热的热流量向温度低的方向传递，即与温度梯度的方向相反。比热流量是个向量，即它有大小和方向。,7,1.1.1 传热的基本方式,如果比热流量的分量和（X、Y、Z）坐标系相联系，那么X、Y、Z方向的热流量分量应是：qx=-,qy=-,qz=-（1-1-2）X Y Z比热流量 q=iqx+jqy+kqz（1-1-3）导热系数 物理意义：沿导热方向的单位长度上，温度减低,物质所容许通过的热流量。方向性：大多数液体和固体属于各向同性的物质。各向异性材料的导热系数具有方向性，如石墨。温度函数：值还随温度而变化。大多数金属的导热系数随温度的升高而降低。大多
5、数液体（水和甘油除外）其导热系数随温度的升高而降低。气体的导热系数随温度的升高而增加。,8,1.1.1 传热的基本方式,对流对流是流体（气体和液体）中温度不同的各部分相互混合的宏观运动引起热量传递的现象。工程上最具有实际意义的是：相对运动着的流体与所接触的固体壁面之间的热量交换过程，一般称为对流换热。工程上在研究固体壁面和流体之间的对流换热时，除了高度稀薄的气体外，人们不去注意流体的单个质点，而把流体看成是连续介质。实际的流体总有粘性，流动时，受粘性和壁面摩擦的影响，在靠近壁面附近的流体将降低流速，在壁面上完全被滞止不动，即X=0时，=0，如图1-1所示。因此，热量从壁面传给贴壁的那部分流体，
6、将依靠导热。,T(K)V(m/s)Tw qwc Tf V X(m)图1-1邻近壁面的流体速度分布和温度分布,9,1.1.1 传热的基本方式,流体和固体壁接触面上的“相”际热流密度为：qwc=-f x=+0（1-1-4）X式中：qwc热流密度，W/m2 f流体的导热系数，W/(m.K)T流体的温度，K液体的导热系数较小，气体的导热系数更小，所以受热时，在贴壁处的流体温度势必沿着轴的反方向急剧升高。随着离壁面的距离的增加，流体的流动将带走更多的热量，使轴方向的温度梯度连续下降，温度分布趋向平坦化。正是通过这种导热和对流的共同作用，使热量在流体内部得到传播，越临近壁面，导热越起主导作用。图1-1所示
7、，假如厚度为的贴壁静止膜，膜内温度线性地从壁面温度TW降到远离壁面，尚未被加热的流体温度Tf，则 Tf-TW qwc=-f（1-1-5）无界对流时壁面与流体的换热，钢锭与周围空气的对流换热属于这种情况。,10,1.1.1 传热的基本方式,流体在管和槽道内部的流动，称为有界对流，这时不存在远离壁面，尚未被加热的流体温度，则采用截面平均温度作为流体的特征温度Tf，则qWC=aC(TW-Tf),W/m2（1-1-6）这就是所谓的“牛顿冷却定律”。式中：aC为放热系数，W/(m2.K)其实“牛顿冷却定律”并不是表达对流换热现象本质的物理定律。凡能影响流体流动的各种因素，包括流体的种类和状态、运动的速度
8、、流道的形状和大小，以及固体壁面粗糙度等，都会对aC值产生影响。式1-1-6只不过形式地把放热过程的一切复杂性和计算上的困难，都转移到并且集中在放热系数这样一个物理量上罢了。aC代表流体和所接触的固体表面之间温度每相差，该流体与表面之间“相”际热流量的大小。运用式1-1-6可以进行对流换热的计算。但由于对流换热的复杂性，该式中的放热系数aC需从相应的准则方程式求出。准则方程式是针对不同对流换热情况，在综合实验结果的基础之上，运用相似理论将表征某现象的物理量整理成一些相似准则，用因次分析法得到的各种类型的表达式。,11,辐射只要温度高于绝对零度，任何物体都随时向周围空间辐射能力。辐射用斯蒂芬玻耳
9、兹曼定律表达：E=，w/m2 1-1-7 式中：物体的辐射率或黑度,(0-1)；斯蒂芬玻耳兹曼常数或绝对黑度的辐射常数；5.6710-8 W/m2 K4温度,K。实际上，辐射往往涉及二个物体间辐射热交换。如果二个物体的表面温度不同，中间被空气所隔开，T1T2时，则相互辐射的结果，表面温度为T1的物体发射出去的辐射热超过了吸收来自表面温度为T2的物体辐射热，引起辐射换热的热流量则为：Q12=12(1-2)F112 或q12=12(1-2)12 1-1-8 式中：12 物体1与2综合黑度；F1物体1的表面积；12 物体的表面向外发射出去的辐射热量中，能投射到物体表面上的份额，称为角系数,(0-1)
10、。,1.1.1 传热的基本方式,12,1.1.1 传热的基本方式,壁面在气体中冷却，存在辐射换热和对流换热。考虑到壁面与气体之间存在着辐射换热，其热流密度为Qr=arF(Tw-Tf),w 或qwr=ar(Tw-Tf),w/m2 1-1-9式中：qwr单位面积的辐射量，w/m2 ar辐射放热系数，w/(m2.k)Tw辐射物体表面温度，k Tf透明的气体介质的特征温度，k考虑到壁面与气体之间还存在着对流换热，其热流密度为qwc=ac(Tw-Tf),w/m2 1-1-10 由壁面传走的总热流密度qw应是qwr和qwc二者之和qw=ac(Tw-Tf)+ar(Tw-Tf)1-1-11 令a0=ac+ar
11、 则qw=a0(Tw-Tf),w/m2 1-1-12用辐射放热系数ar，可以形式地把辐射换热折合成对流换热，用总放热系,13,1.1.1 传热的基本方式,数a0兼顾辐射换热的影响，从而有利于简化复杂传热的分析和计算。如远离表面的外界表面温度趋于环境温度Tf，并且12=1时，由式1-1-8得 qwr=12(w-f)1-1-13由式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)4m3 ar 4 12 m3,w/(m2
12、.k)1-1-16从此可见，(w-f)降为零时，ar并非零值，而是以4 12 m3。随着温差的扩大和平均温度m的升高，ar值将迅速增加。由于ac随温差的变化较小，在高温范围和大温差情况下，ar有可能成为a0的主要组成部分。ar与气体的运动状况无关，而ac与气体速度的降低而减小，在气体自然对流的情况下，ac 5,w/(m2.k)，即使m只有300K，ar就可能占a0的一半左右。,14,1.1.2 导热的偏微分方程式,钢锭一般属于各向同性的物质，其加热或冷却过程数学模拟计算依据的基本数学模型，是不稳定导热偏微分方程。下面讨论各向同性材料导热方程式的建立。直角坐标系导热偏微分方程导热偏微分方程的建立
13、，是通过考察处于导热过程中的物质的微元体积（xy z)的能量平衡来建立。如图1-2所示。在时间内，通过六个面的导热所获得的能量，加上微元体内产生的内热源热量，要等于微元体积内物质积蓄热量的改变，即温度升高或降低。,图1-2直角坐标系导热方程式的微元体,15,1.1.2 导热的偏微分方程式,图1-2中，微元控制体尺寸x、y、z，按照傅立叶导热定律，在X方向流入微元体左表面的热流可表示为：TdQx=-(y z)，W 1-1-17 X 式中：导热系数，W/(m.k)y z X方向微元体表面积，m2T X方向的温度梯度，k/mX在X方向流出微元体右表面的热流，可以应用泰勒级数展开，保留级数的第一项和第
14、二项而得到：TdQx+x=dQx+(dQx)x=dQx+-(y z)x X X X T=dQx-()x y z 1-1-18 X X 在X方向导热的净热流为 T dQx-dQx+=()x y z 1-1-19 X X,16,1.1.2 导热的偏微分方程式,用同样的方法，可以得出Y,Z 方向与式-19相类似的导热的净热流方程式，即 T dQy-dQy+y=()x y z 1-1-20 Y Y T dQz-dQz+z=()x y z 1-1-21 Z Z 在三个坐标方向净热流的总和为：T T T()+()+()x y z 1-1-22 X X Y Y Z Z如果单位时间、单位空间所产生的热量为Q(
15、x,y,z,)，那么微元体的发热量为：Q(x,y,z,)x y z 1-1-23由于导热传进微元体的净热流式1-22和微元体内产生的热量式1-23一起用于增大微元体的内能。微元体的内能的增加表现在微元体能量存储随时间的变化率，即 T Cp x y z 1-1-24,17,1.1.2 导热的偏微分方程式,式中：密度，kg/m3 Cp比热，J/(kg.k)时间，s 对微元体进行能量平衡，使能量存储的时间变化率与由导热引起的流入微元体的净热流和微元体内产生的热量之和相等，得下式：T T T T Cp=()+()+()+Q 1-1-25 X X Y Y Z Z 圆柱坐标系导热偏微分方程实际问题常常涉及
16、柱面对称问题，且边界条件给定在一个表面上，因此表面具有一个坐标保持不变的性质。在这种情况下，采用圆柱坐标系是合适的。,图1-3圆柱坐标系导热方程式的微元体,18,1.1.2 导热的偏微分方程式,图1-3所示圆柱坐标系，直接按内外两个圆弧面和其它四个平面组成的微元体积，在导热过程中热量必须按收支平衡的原则导出，此时，微元体积为：dv=(r d)dzdr 1-1-26 沿内圆弧面流入微元体积的热流：TdQr=-(r ddz)1-1-27 r沿外圆弧面流出微元体积的热流：TdQr+dr=dQr+(dQr)dr=dQr+-(r ddz)dr r r r T 1 T=dQr+(-r)ddzdr=dQr-
17、(r)dv 1-1-28 r r r r r 沿平面流入微元体积的热流：TdQ=-(dr dz)1-1-29r沿+d平面流出微元体积的热流：TdQ+d=dQ+(dQ)rd=dQ+-(drdz)rd r r r,19,1.1.2 导热的偏微分方程式,1=dQ-()dv 1-1-30 r2 沿z面流入微元体积的热流：TdQz=-(r ddr)1-1-31 z沿z+dz面流出微元体积的热流：TdQz+dz=dQz+(dQz)dz=dQz+-(r ddr)dz z z z T=dQz-()dv 1-1-32 z z 根据直角坐标系导热微分方程推导的思路，可得到：1 T 1 T T(r)+()+()+Q
18、=Cp 1-1-33 r r r r2 z z,20,1.1.2 导热的偏微分方程式,球坐标系导热偏微分方程对于球坐标系，如图1-4所示。,图1-4 球坐标系导热方程式的微元体,由内、外两个球面、两个圆弧面和两个平面所组成的微元体积为：dv=(r d)（rSind)dr 1-1-34,21,1.1.2 导热的偏微分方程式,沿半径方向流入微元体积的热流：TdQr=-(r Sind.rd)1-1-35 r沿半径方向流出微元体积的热流：TdQr+dr=dQr+(dQr)dr=dQr+-(rSind.rd)dr r r r T 1 T=dQr+(-r2)Sind.ddr=dQr-(r2)dv 1-1-
19、36 r r r2 r r 沿方向流入微元体积的热流：TdQ=-(rSind.dr)1-1-37 r沿方向流出微元体积的热流：TdQ+d=dQ+(dQ)rd=dQ+-(rSinddr)rd r r r T T=dQ-(Sin)rddr.rd=dQ-(Sin)dv r2 Sinr2 1-1-38,22,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-()dv 1-1-40(r2Sin2)同
20、理可整理得：T T T(r2)+(Sin)+()+Q=Cp r2r r r2Sin r2Sin 1-1-41,23,1.2 导热方程的有限差分解法,求解不稳定导热偏微分方程的数值解法，主要有：有限差分解法、有限元法、边界元法三类。边界元法正在研究和完善之中，目前常用的是有限差分解法和有限元法。我们专门讨论有限差分解法的数学基础，数值方程的建立，差分方程的稳定性和收敛性等问题。有限差分解法是用差分方程近似地代替微分方程，建立差分方程有直接法和能量平衡法两种。直接代换法直接代换法是从微分形式出发，用差商直接代换微商的办法建立差分方程。数学基础微商和差商的定义若T(x)是连续函数，则其导数为：dT(
21、x+x)-T(x)=lim=lim 1-2-1dx x0 x x0 x,24,数学基础,1-2-1式右边/x是有限的差商，x和都不为零。式左边d/dx是T/x当x0时的差商，称之为微商。在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)的关
22、系如下式所示：dT(x)2 d2T(x)n dnT T(x+x)=T(x)+x+1-2-2 dx 2!dx2 n!dxn整理后得：T T(x+x)-T(x)dT x d2T(x)n-1 dnT=+1-2-3 x x dx 2!dx2 n!dxn从上式可知，T(x)在x处的差商T/x等于函数T(x)在x处的各阶导数的线性组合，只能是近似地等于差商。两者之间也必然有偏差。,25,数学基础,图1-2-1表示了一阶差商与一阶微商之间的关系。用不同方法得到的差商去代替微商，它们带来的误差是不同的。即向前差商：dT/dxT(x+x)-T(x)/x 1-2-4向后差商：dT/dxT(x)-T(x-x)/x
23、1-2-5中心差商：dT/dx1/2T(x+x)+T(x)/x+T(x)-T(x-x)/x=T(x+x)-T(x-x)/2x 1-2-6,图1-2-1差商与微商,26,数学基础,按照泰勒级数展开，T(x)与T(x+x)的关系如下式所示：dT(x)2 d2T(x)n dnT T(x+x)=T(x)+x+1-2-7 dx 2!dx2 n!dxn整理后得：T(x+x)-T(x)dT x d2T-=+=O(x)1-2-8 x dx 2!dx2 即向前差商的偏差是截去了泰勒级数展开式中的高阶项而引起的，常称“截断误差”，其截断误差为与x同级的小量O(x)。同理dT(x)2 d2T(x)3 d3T T(x
24、-x)=T(x)-x+-+1-2-9 dx 2!dx2 3!dx3整理后得：T(x)-T(x-x)dT x d2T-=-+=O(x)1-2-10 x dx 2!dx2 即向后差商的截断误差为与x同级的小量O(x)。,27,数学基础,由式1-2-7减式1-2-9，将2 x dT/dx移至等式左边，两边再除以2 x，得：T(x+x)-T(x-x)dT(x)2 d3T-=+=O(x)2 1-2-11 2x dx 3!dx3 即中心差商的截断误差为与(x)2同级的小量O(x)2。当x固定时，上述一阶差商一般仍是x的函数，对它们还可以求差商。这种一阶差商的差商称为二阶差商，它是二阶微商的近似，常用向前和
25、向后差商来二阶微商，即d2T T(x+x)-T(x)T(x)-T(x-x)-/x dx2 xx T(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)2 1-2-13(x)2 dx2即二阶差商的截断误差为与(x)2同级的小量O(x)2。,28,1.2.1.2 建立内节点差分方程,一维系统假定有一高宽无限(即高宽方向上无热流)，厚度为L的平板，T=f(x,)即温度是x方向位置和时间的函数，所谓一维系统是指几何空间为一维。初始时刻=0，T=T0，为了简化，考虑无内热源，、Cp均为常数。选取网
26、格点间距x和时间步长将研究对象离散化。显式差分格式 T T 一维不稳定导热方程为 Cp=()X X 该方程在区域0，0 xL内全部点都成立，如图1-2-2所示。将方程应用于内节点 i 可写成：T p 2T p Cp()=()1-2-13 i X2 i上式等号左端的微分式用温度对时间的一阶向前差商来近似，即：T p Tp+1i-Tpi()=+O()1-2-14 i 上式等号右端的微分式用温度对空间的二阶差商来近似，即：T p Tpi+1-2Tpi+Tpi-1()=+O(x)2)1-2-15 x2 i(x)2,29,1.2.1.2 建立内节点差分方程,图1-2-2一维显式差分将式1-2-14和式1
27、-2-15代入式1-2-13得：Tp+1i Tpi Tpi+1-2Tpi+Tpi-1 Cp+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为傅立叶数。,30,1.2.1.2 建立内节点差分方程,F0的数值小意味着加热或冷却此物体所需要的时间长，反之，所需要的时间短，F0是一个无因次数字。当初始条件和边界条件已知时，用式1-2-17就可模拟区域内各节点随时间增长的温度值Tpi(i=2,3,，N-1；p=1,2,3,）隐式差分格式 一
28、维隐式差分如图1-2-3所示，将一维不稳定导热微分方程应用于内节点i，则：T p 2T p+1 Cp()=()1-2-18 i X2 i 式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+1i Tpi Tp+1i+1-2Tp+1i+Tp+1i-1 Cp+O()=+O(x)2)
29、1-2-19(x)2若在上式中去掉O()和O(x)2),整理得：Tp+1i=Tpi+F0(Tp+1i+1-2Tp+1i+Tp+1i-1)1-2-20,31,1.2.1.2 建立内节点差分方程,式中F0=/Cp(x)2。比较式1-2-17和1-2-20，可以看出显式差分格式的突出优点是每个节点方程都可以独立求解，整个计算过程十分简便。但是，的选取要受到限制，有时为了满足差分格式稳定性条件，可能选得很小，使计算工作量加大。隐式差分格式的最大优点是，它对任意F0值都是稳定的。这种稳定是绝对的，即不受边界条件、x、热物参数的影响，于是可以选的较大，计算速度加快。但是，对于节点i，只从式1-2-20不能
30、独立求解，必须涉及Tp+1i+1、Tp+1i、Tp+1i-1的联立线性代数方程组才能求解，也就是说，它含有三个未知数。时间步长每前进一步，从坐标0 xL，网格点1iN整个区域的每个点，上述方程都要列出一次（见图1-2-3)。因此，向一个新的时间步长每移动一步就必须解一个方程组。当按顺序列出这些方程时，除了要的第一点和最后一点的方程只有两个未知数外，其余每一个方程都含有三个未知数，于是方程是三对角的。对于这种情况，已有适用的求解程序，后面将讨论。关于差分方程稳定性将在后面讨论。,32,1.2.1.2 建立内节点差分方程,图1-2-3一维隐式差分二维系统显式和隐式差分格式建立的方法和两种差分格式的
31、特点在前面讨论过，对二维系统同样适用。为简略起见，在此只讨论二维系统显式差分方程的建立。为使问题简化，仍然假定热物性值为常数，考虑无内热源。首先把所讨论的区域离散化，如图1-2-4所示。,33,1.2.1.2 建立内节点差分方程,网线：用平行于X、Y轴的 直线（网线），进 行空间离散化节点：网线与网线的交点步长：节点间的距离图1-2-4 二维均质物体的分割x与y可以是不变的常量，即等步长，也可以是变量（即在区域内的不同处是不同的），即变步长。如果区域内各点处的温度梯度相差很大，则在温度变化剧烈处，网格布得密些，在温度变化不剧烈处，网格布得疏些。至于网格多少，步长取多少为宜，要根据计算精度与计算
32、工作量等因素而定。对于0 xL1和0yL2的矩形区域内，将二维不稳定导热方程式应用于节点（i,j)可写成：T p 2T 2T p Cp()=(+)1-2-21 i,j x2 y2 i,j,34,1.2.1.2 建立内节点差分方程,式1-2-21左边一阶偏微商用温度对时间一阶向前差商近似：T p Tp+1i,j Tpi,j()=+O()1-2-22 i,j 式1-2-21右边的二阶偏微商所对应的差商可表示为：T p Tpi+1,j-2Tpi,j+Tpi-1,j()=+O(x)2)1-2-23 x2 i,j(x)2 T p Tpi,j+1-2Tpi,j+Tpi,j-1()=+O(y)2)1-2-2
33、4 y2 i,j(y)2当、x、y较小时，忽略O()、O(x)2)、O(y)2)项，当x=y时，即x、y方向网格划分步长相等。将上三式代入式1-2-21则得到节点(i,j)的差分方程：Tp+1i,j=Tpi,j+F0(Tpi+1,j+Tpi-1,j+Tpi,j+1+Tpi,j-1-4Tpi,j)1-2-25 式中F0=/Cp(x)2。上述条件，隐式差分格式，则除Tpi,j项外，其它项的Tp 改为Tp+1。,35,1.2.1.2 建立内节点差分方程,三维系统和二维系统的假定相同，热物性值为常数，无内热源。设x=y=z，即均匀网格间距，将三维不稳定导热方程式应用于节点(i,j,k)，可以写成：T
34、p 2T 2T 2T p Cp()=(+)1-2-26 i,j,k x2 y2 z2 i,j,k即可得到节点(i,j,k)的近似式 Tp+1i,j,k Tpi,j,k Tpi+1,j,k-2Tpi,j,k+Tpi-1,j,k Cp()=+(x)2 Tpi,j+1,k-2Tpi,j,k+Tpi,j-1,k Tpi,j,k+1-2Tpi,j,k+Tpi,j,k-1+1-2-27(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-6T
35、pi,j,k)1-2-28 式中F0=/Cp(x)2。如果写成隐式差分格式，则除Tpi,j,k项外，其它项的Tp 改为Tp+1。,36,1.2.2 能量平衡法,1.2.2 能量平衡法能量平衡法是不再从微分方程出发用差商代替微商的办法建立差分方程。而是将导热的基本定律直接近似，局部地应用于围绕每个节点的各单元控制体积，用热量平衡法建立差分方程。由于它的基本思想是把能量守恒原则应用到每个单元体，因此常称为能量守恒原则的有限差分法。1.2.2.1 建立内节点差分方程一维系统高宽无限(即高宽方向上无热流)，一定厚度的平板，T=f(x,)，初始时刻=0，T=T0，为了简化，考虑无内热源，、Cp均为常数。
36、沿板厚方向（即热流方向）把均质物体分割为若干单元体，各单元体的端面积为F，几何步长为x，时间步长为（用P上角标表示第几个时间步长）。单元体中心为节点，节点温度代表该单元体的温度。二相邻节点间的导热面积是F。如图1-2-5所示。图1-2-5 一维均质物体分割,37,1.2.2.1 建立内节点差分方程,考虑内节点 i，从时间等于开始的时间间隔内，通过截面F从(i-1)单元体流入到i单元体的热量为：Tpi-1 Tpi Qi-1i=F()1-2-29 x 从(i+1)单元体流入到i单元体的热量为：Tpi+1 Tpi Qi+1i=F()1-2-30 x 由于热量流入 i 单元体，引起 i 单元体的内能变
37、化。在时间内，i 单元体的内能变化为：U Tp+1i Tpi=Cp(Fx)1-2-31 因为能量平衡或守恒，所以Q=U/，即：Tp+1i Tpi Tpi-1 Tpi Tpi+1 Tpi Cp(Fx)=F()+F()1-2-32 x x整理可得：Tp+1i=Tpi+F0(Tpi+1-2Tpi+Tpi-1)1-2-33 式中F0=/Cp(x)2。上式和直接代换法式1-2-17完全一样。,38,1.2.2.1 建立内节点差分方程,二维系统如果物体的形状使之只在X、Y方向上有热流，Z方向上无热流，或相对于X、Y方向可忽略不计，这种情况就可以当作二维系统处理。对于如图1-2-6所示的二维矩形区域，划分成
38、有限个小单元体。图中虚线所示的(i,j)单元体及周围四个相邻的单元体，用于推导内单元体差分方程。每个小单元体的几何步长为x和y；时间步长。为使问题简化，仍然假定热物性值为常数，考虑无内热源。图1-2-6 二维系统的分割,39,1.2.2.1 建立内节点差分方程,对于(i,j)单元体，时间间隔内的内能变化为：U Tp+1i,j Tpi,j=Cp(xy.1)1-2-34 在时间内从周围四个相邻单元体流入(i,j)单元体的热量分别为：Tpi-1,j Tpi,j Qi-1,ji,j=(y.1)()1-2-35 x Tpi+1,j Tpi,j Qi+1,ji,j=(y.1)()1-2-36 x Tpi,
39、j-1 Tpi,j Qi,j-1i,j=(x.1)()1-2-37 y Tpi,j+1 Tpi,j Qi,j+1i,j=(x.1)()1-2-38 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)1-2-39 式中F0=/Cp(x)2。上式和直接代换法式1-2-25完全一样。三维系统同理可得和式(1-2-28)一样的节点差分方程。,40,1.2.2.2 能量平衡法与直接代换法的比较,1.2.2.2 能量平衡法与直接代换法的比较尽管用能量平衡法与直接代换法所得到的差分方程在形式上完全一样，但它
40、们所依赖的基础有所不同。用直接代换法得出差分方程时，要求区域内温度分布连续而光滑，因为按泰勒级数展开，要求温度T(x)在x处的各阶导数都存在，连续而光滑即可导。而能量平衡法只要求每个单元体边界上热流可积，即可以算出每个单元体边界上的热量，它比温度分布连续光滑的要求显然放宽了。如果区域由导热系数截然不同的两种物质构成，并在交界面上存在接触热阻，或在交界面上存在间隙，那么在那些交界面上或间隙处，温度分布就不连续或不光滑。在这种情况下，可用能量平衡法建立差分方程，只需在单元划分时，使那些交界面正好落在单元的边界即可，如图1-2-7所示。图1-2-7 异种物质接触时单元的划分图1-2-8 界面温度分布
41、示意图,41,1.2.2.2 能量平衡法与直接代换法的比较,现假定二维矩形分割，交界面两侧的物体导热系数不同，讨论相邻两节点(i+1,j)节点流出(i,j)节点的热量公式的建立。当异种物质接触时，除传导传热外，还要考虑界面热阻Rb，则 Tpi+1,j Tp2 Qi+1,jT2=i+1,j(y.1)()1-2-40 x/2 Tp2 Tp1 QT2T1=(y.1)()1-2-41 Rb Tp1 Tpi,j QT1i,j=i,j(y.1)()1-2-42 x/2因为三个Q相似，令其相等，由式1-2-40和式1-2-42 得：Q x/2=Tpi+1,j Tp2 1-2-43 y.1 i+1,j Q x
42、/2=Tp1-Tpi,j 1-2-44 y.1 i,j Q x/2 x/2(+)=(Tpi+1,j Tpi,j)-(Tp2-Tp1)1-2-45 y.1 i+1,j i,j 将式1-2-41代入上式得：,42,1.2.2.2 能量平衡法与直接代换法的比较,y.1 Q=(Tpi+1,j Tpi,j)1-2-46 Rb+x/(2i+1,j)+x/(2i,j)与直接代换法相比，能量平衡法的主要优点是：在二维系统中能使用三角形、四角形，在三维系统中能使用四、五、六面体等多种分割单元，因此能解复杂形状的问题；对于每个节点单元能指定热物性值，易于求解由多种物质组成的系统的问题。能量平衡法的缺点是：输入数据
43、量多，程序复杂，计算时间较长。所以，对于矩形或圆柱形等简单形状，能进行有规律分割的系统，采用直接代换法是有利的。可是，能有规律分割时，能量平衡法的程序与直接代换法的程序相同。一般说，能量平衡法更实用些。1.2.3 边界节点差分方程的建立物体内部的温度场必然受到物体表面条件的影响，如表面散热条件差，表面温度高，物体内部的温度场平缓。反之，物体内部温度场的变化也影响着表面上的条件，如物体内部的导热系数大，物体内部温度场平缓，使表面温度升高，表面散热量也大。为了数值模拟，必须建立边界节点的差分方程。在此，先讨论一般换热条件的边界节点差分方程的建立。,43,1.2.3 边界节点差分方程的建立,图1-2
44、-9所示了绝热、给定热流密度qr、对流、定温、辐射换热等五种边界表面的二维矩形区域网格，现用能量平衡法建立各种换热边界条件下的节点方程。图1-2-9 各种换热边界表面图1-2-10 绝热边界单元绝热边界把图1-2-9中AB面上任一节点(i,j)及其相邻节点取出分析，如图1-2-10所示。Tpi+1,jTpi,j Tpi,j-1 Tpi,j(y.1)()+(x/2.1)()+x y Tpi,j+1Tpi,j Tp+1i,j Tpi,j(x/2.1)()+0=Cp(x/2y.1)1-2-47 y,44,绝热边界,若x=y，经整理后得：Tp+1i,j=Tpi,j+F0(2Tpi+1,j+Tpi,j-
45、1+Tpi,j+1-4Tpi,j)1-2-48 式中F0=/Cp(x)2。给定热流密度qr边界把图1-2-9中BC面上任一节点(i,j)及其相邻节点取出分析，如图1-2-11所示。用能量平衡法建立给定热流密度qr边界节点方程。Tpi-1,jTpi,j(y/2.1)()+x Tpi+1,jTpi,j(y/2.1)()+x Tpi,j+1Tpi,j(x.1)()+qr(x.1)y Tp+1i,j Tpi,j=Cp(xy/2.1)1-2-49 若x=y，经整理后得：图1-2-11 给定热流密度边界单元 Tp+1i,j=Tpi,j+F0(Tpi-1,j+Tpi+1,j+2Tpi,j+1-4Tpi,j)
46、+2qr/(Cpx)1-2-50 式中F0=/Cp(x)2。,45,1.2.3.3 对流边界,1.2.3.3 对流边界把图1-2-9中CD面上任一节点(i,j)及其相邻节点取出分析，如图1-2-12所示。在CD面上，ac为对流换热系数，Tf为物质周围介质温度。用能量平衡法建立对流边界节点方程。Tpi-1,jTpi,j(y.1)()+x Tpi,j-1Tpi,j(x/2.1)()+y Tpi,j+1Tpi,j(x/2.1)()+ac(y.1)(Tf-Tpi,j)y Tp+1i,j Tpi,j=Cp(x/2y.1)1-2-51 若x=y，经整理后得：图1-2-12对流边界单元 Tp+1i,j=Tp
47、i,j+F0(2Tpi-1,j+Tpi,j-1+Tpi,j+1-4Tpi,j)+2ac(Tf-Tpi,j)/(Cpx)1-2-52 式中F0=/Cp(x)2。1.2.3.4 给定温度边界图1-2-9中AD面上给定温度Tw，把AD面上任一节点(i,j)及其相邻节点取出分析，如图1-2-13所示。对这些单元直接写出结果：Tp+1i,j=Tw 1-2-53,46,1.2.3.5 辐射边界,把图1-2-9中AD面上受辐射热流的作用，把AD面上任一节点(i,j)及其相邻节点取出分析，如图1-2-13所示。为黑度1，为斯蒂芬玻耳兹曼常数，Tf为物质周围介质温度。Tpi-1,jTpi,j(y/2.1)()+
48、x Tpi+1,jTpi,j(y/2.1)()+x Tpi,j-1Tpi,j(x.1)()+y(x.1)Tf4(Tpi,j)4 Tp+1i,j Tpi,j=Cp(xy/2.1)1-2-54 若x=y，经整理后得：图1-2-13定温或辐射边界单元 Tp+1i,j=Tpi,j+F0(Tpi-1,j+Tpi+1,j+2Tpi,j-1-4Tpi,j)+2Tf4(Tpi,j)4/(Cpx)1-2-55 式中F0=/Cp(x)2。,47,1.2.3.6 混合边界,图1-2-9中四个角上的边界单元，除AD边定温边界条件时，A、D 两个拐角节点的温度为一定值（Tw）的情况以外，在其余换热边界条件下，它们的两个
49、外边界从属于两种边界条件，如C拐角节点，从图1-9可知，它的一个边界处于对流换热，另一个边界处于给定热流密度qr。在这种情况下，将C拐角节点及其相邻节点取出分析，如图1-14所示，按能量平衡法可得：Tpi-1,jTpi,j(y/2.1)()+x Tpi,j+1Tpi,j(x/2.1)()+y qr(x/2.1)+ac(y/2.1)(Tf Tpi,j)Tp+1i,j Tpi,j=Cp(x/2y/2.1)1-2-56 图1-2-14混合边界单元 若x=y，经整理后得：Tp+1i,j=Tpi,j+F0(2Tpi-1,j+2Tpi,j+1-4Tpi,j)+2qrF0 x/+2acF0 x(TfTpi,
50、j)/式中F0=/Cp(x)2。1-2-57,48,1.2.3.6 混合边界,对图1-2-9中A、B、D等拐角节点，在绝热、给定热流密度qr、对流和辐射换热等混合边界条件下，都可以按上述类似方法，分别建立它们的边界单元差分方程。对于物体出现内拐点（凹角）的情况，也可按能量守恒定律得出凹角节点的差分方程。如图1-2-15所示的凹角节点（i,j），若凹角两边处于对流换热边界条件时，把（i,j）节点及其相邻节点取出分析，则可得：Tpi-1,jTpi,j(y.1)()+x Tpi+1,jTpi,j(y/2.1)()+x Tpi,j-1Tpi,j(x/2.1)()+y Tpi,j+1Tpi,j(x.1)
