钢铁冶金过程计算机仿真.ppt

钢铁冶金过程计算机仿真钢铁冶金过程计算机仿真温度场数值模拟温度场数值模拟北京科技大学冶金学院　刘青北京科技大学冶金学院　刘青北京科技大学冶金学院　刘青北京科技大学冶金学院　刘青1前　言前　言１、钢铁冶金加工过程１、钢铁冶金加工过程１、钢铁冶金加工过程１、钢铁冶金加工过程　　　炼铁－炼钢（铸造、锻造、焊接、挤、压、拉、拔、热处理）            专业设置２、数值模拟２、数值模拟２、数值模拟２、数值模拟　　　两个不同领域的现象，能用同一数学形式来描述，称这两个现象彼此是可模拟　　　模拟的方法是把一个领域内求解的问题过渡到另一个领域中去解决　　　随着计算机和计算数学的发展，用计算机数值计算法求解问题的计算精度已经达到很接近解析解，此法称为计算机数值模拟           仿真的定义，与模拟的差异３、工艺优化３、工艺优化３、工艺优化３、工艺优化 　　　计算机数值模拟的最终目的是解决工艺优化设计问题一旦全面实现了冶金加工过程的计算机数值模拟，材料的加工生产将会产生深刻的变革            CAD/CAM/CAE的发展，波音767飞机的设计和加工            2前　言前　言４、计算机数值模拟的步骤４、计算机数值模拟的步骤４、计算机数值模拟的步骤４、计算机数值模拟的步骤　　①给定研究对象　　　几何条件、物理条件、初始条件、边界条件　　②离散化处理　　　将过程所涉及的区域在空间和时间上进行离散化处理。

　　③建立数值方程　　　建立内节点和边界节点的数值方程　　④选择计算方法　　　选择适当的计算方法求解线性代数方程组　　⑤编制计算机程序　　　编制计算机程序，由计算机算出结果，并用数据、曲线、图形输出　　⑥优化工艺　　　分析计算机模拟的结果，如果结果不理想，调整工艺参数，再进行计算机模拟，直到模拟结果为最佳结果3前　言前　言５、讲解的主要内容５、讲解的主要内容　　　　　　　　通过一、二个例子，让大家亲自经历计算机数值模拟的一般步骤，从中掌握其基本原理和方法，为更好应用计算机数值模拟，优化工艺参数打下基础　　冶金过程中，比较常用、典型的数值模拟有：　　①温度场数值模拟　　②浓度场数值模拟　　③组织与性能模拟451 温度场计算机数值模拟温度场计算机数值模拟1.1 1.1 传热的基本知识传热的基本知识传热的基本知识传热的基本知识1.1.1 1.1.1 传热的基本方式传热的基本方式传热的基本方式传热的基本方式　①导热　　　导热属于接触传热，是连续介质就地传递热量，没有各部部分物质之间宏观的相对位移　　　在不透明固体实体内部，由于各部分物质之间无法作宏观的相对位移，不透明无法传递辐射能，实体保证接触，所以只能依靠导热方式传递热量。

　　　导热的基本定律是傅立叶定律，即导热的比热流量q与温度梯度成正比，即：　　　　　　　　　　　　 　 Ｔ q∝grad T=-λgrad T= - λ——　(W/m2) （1-1-1）　　　　　　　　　　　　 　 n　　　　式中：q  — 传热方向上单位时间、单位面积上所通过的热量，[J/(s.m2)]=W/m2                   λ — 材料的导热系数，W/(m.K)　　　　   Ｔ                                  —— —　温度梯度，K/m　　　　    n　　负号表示导热的热流量向温度低的方向传递，即与温度梯度的方向相反　　比热流量是个向量，即它有大小和方向　　61.1.1 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）•导热系数 　　物理意义：沿导热方向的单位长度上，温度减低１℃,物质所容许通过的热流量。

　　方向性：大多数液体和固体属于各向同性的物质各向异性材料的导热系数具有方向性，如石墨　　温度函数：λ值还随温度而变化大多数金属的导热系数随温度的升高而降低大多数液体（水和甘油除外）其导热系数随温度的升高而降低 气体的导热系数随温度的升高而增加 71.1.1 1.1.1 传热的基本方式传热的基本方式②对流•对流是流体（气体和液体）中温度不同的各部分相互混合的宏观运动引起热量传递的现象•工程上最具有实际意义的是：相对运动着的流体与所接触的固体壁面之间的热量交换过程，一般称为对流换热•工程上在研究固体壁面和流体之间的对流换热时，除了高度稀薄的气体外，人们不去注意流体的单个质点，而把流体看成是连续介质•实际的流体总有粘性，流动时，受粘性和壁面摩擦的影响，在靠近壁面附近的流体将降低流速，在壁面上完全被滞止不动，即X=0时，Ｖ=0，如图1-1所示因此，热量从壁面传给贴壁的那部分流体，将依靠导热T(K)                              V(m/s) Tw         qwc                       Tf                                            V                                        X(m)δ图1-1邻近壁面的流体速度分布和温度分布　　81.1.1 1.1.1 传热的基本方式传热的基本方式•流体和固体壁接触面上的“相”际热流密度为：　 Ｔ qwc= - λf—— x=+0 （1-1-4）　　　　　 X式中：qwc—热流密度，W/m2 　       λf—流体的导热系数，W/(m.K) T—流体的温度，K•液体的导热系数较小，气体的导热系数更小，所以受热时，在贴壁处的流体温度势必沿着Ｘ轴的反方向急剧升高。

随着离壁面的距离Ｘ的增加，流体的流动将带走更多的热量，使Ｘ轴方向的温度梯度连续下降，温度分布趋向平坦化正是通过这种导热和对流的共同作用，使热量在流体内部得到传播，越临近壁面，导热越起主导作用•图1-1所示，假如厚度为δ的贴壁静止膜，膜内温度线性地从壁面温度TW降到远离壁面，尚未被加热的流体温度Tf，则　 Tf-TW qwc= - λf—— 　　（1-1-5）　　　　 δ •无界对流时壁面与流体的换热，钢锭与周围空气的对流换热属于这种情况91.1.1 1.1.1 传热的基本方式传热的基本方式•流体在管和槽道内部的流动，称为有界对流，这时不存在远离壁面，尚未被加热的流体温度，则采用截面平均温度作为流体的特征温度Tf，则　　　qWC=aC(TW-Tf), W/m2 （1-1-6） 这就是所谓的“牛顿冷却定律”　　式中：aC为放热系数，W/(m2.K)•其实“牛顿冷却定律”并不是表达对流换热现象本质的物理定律。

凡能影响流体流动的各种因素，包括流体的种类和状态、运动的速度、流道的形状和大小，以及固体壁面粗糙度等，都会对aC值产生影响•式1-1-6只不过形式地把放热过程的一切复杂性和计算上的困难，都转移到并且集中在放热系数这样一个物理量上罢了aC代表流体和所接触的固体表面之间温度每相差１℃，该流体与表面之间“相”际热流量的大小•运用式1-1-6可以进行对流换热的计算但由于对流换热的复杂性，该式中的放热系数aC需从相应的准则方程式求出准则方程式是针对不同对流换热情况，在综合实验结果的基础之上，运用相似理论将表征某现象的物理量整理成一些相似准则，用因次分析法得到的各种类型的表达式10　③辐射•只要温度高于绝对零度，任何物体都随时向周围空间辐射能力辐射用斯蒂芬—玻耳兹曼定律表达： E=εσ０Ｔ４，w/m2 1-1-7 式中：ε—物体的辐射率或黑度,(0-1)；　　　　σ０—斯蒂芬—玻耳兹曼常数或绝对黑度的辐射常数； 5.67×10-8 W/m2 K4　　　　　　　Ｔ—温度,K。

•实际上，辐射往往涉及二个物体间辐射热交换如果二个物体的表面温度不同，中间被空气所隔开，T1〉T2时，则相互辐射的结果，表面温度为T1的物体发射出去的辐射热超过了吸收来自表面温度为T2的物体辐射热，引起辐射换热的热流量则为：　　 Q1→2= ε12σ０(Ｔ1４-Ｔ2４)F1φ12 　或　　　　　　　　　　　　　q1→2= ε12σ０(Ｔ1４-Ｔ2４)φ12 　　　　　　　　　　 　　1-1-8 式中：ε12 —物体1与2综合黑度；　　　　 F1—物体1的表面积；　　　　φ12 —物体１的表面向外发射出去的辐射热量中，能投射到物体　　　　　２表面上的份额，称为角系数,(0-1)1.1.1 1.1.1 传热的基本方式传热的基本方式111.1.1 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 　则qw=a0(Tw-Tf) ,w/m2 1-1-12　•用辐射放热系数ar，可以形式地把辐射换热折合成对流换热，用总放热系 121.1.1 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+Ｔw2Ｔf +ＴwＴf2 +Ｔf3) 设Ｔm=1/2(Ｔw+Ｔf)，(Ｔw3+Ｔw2Ｔf +ＴwＴf2 +Ｔf3)≈ 4Ｔm3 ar ≈ 4 ε12σ０ Ｔm3　,w/(m2.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的一半左右 131.1.2 1.1.2 导热的偏微分方程式导热的偏微分方程式•钢锭一般属于各向同性的物质，其加热或冷却过程数学模拟计算依据的基本数学模型，是不稳定导热偏微分方程下面讨论各向同性材料导热方程式的建立　①①①①直角坐标系直角坐标系直角坐标系直角坐标系导热导热导热导热偏微分方程偏微分方程偏微分方程偏微分方程　　导热偏微分方程的建立，是通过考察处于导热过程中的物质的微元体积（δxδy δz)的能量平衡来建立如图1-2所示在时间δ内，通过六个面的导热所获得的能量，加上微元体内产生的内热源热量，要等于微元体积内物质积蓄热量的改变，即温度升高或降低xzydQx+δxdQxdQzdQydQz+δzdQy+δy图1-2　直角坐标系导热方程式的微元体141.1.2 1.1.2 导热的偏微分方程式导热的偏微分方程式•图1-2中，微元控制体尺寸δx、δy、 δz，按照傅立叶导热定律，在X方向流入微元体左表面的热流可表示为：　　　　　　　　　 T　　dQx=-λ(δy δz)—　，W 1-1-17　　　　　　　　　 X 式中： λ — 导热系数，W/(m.k)                   δy δz — X方向微元体表面积，m2　　　　　T　　 — — X方向的温度梯度，k/m　　　　　X•在X方向流出微元体右表面的热流，可以应用泰勒级数展开，保留级数的第一项和第二项而得到：　　　　　　　   T　　dQx+δ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 151.1.2 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’(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  　　　　　　161.1.2 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 ②②②②圆柱坐标系圆柱坐标系圆柱坐标系圆柱坐标系导热导热导热导热偏微分方程偏微分方程偏微分方程偏微分方程•实际问题常常涉及柱面对称问题，且边界条件给定在一个表面上，因此表面具有一个坐标保持不变的性质。

在这种情况下，采用圆柱坐标系是合适的图1-3　圆柱坐标系导热方程式的微元体171.1.2 1.1.2 导热的偏微分方程式导热的偏微分方程式•图1-3所示圆柱坐标系，直接按内外两个圆弧面和其它四个平面组成的微元体积，在导热过程中热量必须按收支平衡的原则导出，此时，微元体积为：　　dv=(r dφ)dzdr 1-1-26 •沿内圆弧面流入微元体积的热流：　　　　　T　　dQr=-λ—(r dφdz)　 1-1-27　　　　　 r•沿外圆弧面流出微元体积的热流：　　　　　　　　   T　　dQr+dr= dQr + —(dQr)dr = dQr + —[-λ—(r dφdz)]dr 　　　　　　　　　 r r r 　　　　　　  T 　　　 1  T　　 = dQr + —(-λr—) dφdzdr = dQr- — —(λr—) dv 　1-1-28　　　　　　　 r r r r r •沿φ平面流入微元体积的热流：　　　　　　T　　dQφ=-λ——　(dr dz) 1-1-29　　　　　　rφ•沿φ+dφ平面流出微元体积的热流：　　　　　　　　  　　  T　　dQφ+dφ= dQφ + — (dQφ)rdφ = dQφ + — [-λ— (drdz)] rdφ 　　　　　　　　　 rφ 　 rφ rφ　　181.1.2 1.1.2 导热的偏微分方程式导热的偏微分方程式　　 1  　 Ｔ 　　 = dQφ - — —(λ—) dv 　　　　　 　　　　　　　 1-1-30 　　　　　 r2 φ φ•沿z面流入微元体积的热流：　　　 　T　　dQz=-λ—(r dφdr) 1-1-31　　　　　 z•沿z+dz面流出微元体积的热流：　　　　　　　　   T　　dQz+dz= dQz + —(dQz)dz = dQz + —[-λ—(r dφdr)]dz 　　　　　　　　　 z z z 　　　　　　  T 　　 = dQz - —(λ—) dv 　 1-1-32　　　　　　　 z z •根据直角坐标系导热微分方程推导的思路，可得到：　 1  T 1  Ｔ  T T　 — —(λr—) + — —(λ—) + —(λ—) + Q’ = ρCp  — 1-1-33　 r r r r2 φ φ z z  191.1.2 1.1.2 导热的偏微分方程式导热的偏微分方程式　③③③③球坐标系球坐标系球坐标系球坐标系导热导热导热导热偏微分方程偏微分方程偏微分方程偏微分方程•对于球坐标系，如图1-4所示。

图1-4  球坐标系导热方程式的微元体•由内、外两个球面、两个圆弧面和两个平面所组成的微元体积为： dv=(r dθ)（rSinθdφ)dr 1-1-34201.1.2 1.1.2 导热的偏微分方程式导热的偏微分方程式•沿半径方向流入微元体积的热流：　　　　 　T　　dQr=-λ—(r Sinθdφ.rdθ) 1-1-35　　　　　 r•沿半径方向流出微元体积的热流：　　　　　　　　   T　　dQr+dr= dQr + —(dQr)dr = dQr + —[-λ—(rSinθdφ.rdθ)]dr 　　　　　　　　　 r r r 　　　　  T 　 1  T　　 = dQr+—(-λr2—)Sinθdφ.dθdr=dQr- — —(λr2—)dv 1-1-36 　　　　　 r r r2 r r •沿θ方向流入微元体积的热流：　　　　　　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θ　 θ　　　　　　　　　　　　　　　　　　　　　　　　　　　　　 1-1-38　　 211.1.2 1.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 1-1-40　　　　 　　　　　(r2Sin2θ)φ　φ•同理可整理得：  T  Ｔ  　　T T　——(λr2—)+—————(λSinθ—)+—————(λ—)+Q’ = ρCp —　 r2r r r2Sinθθ　 θ r2Sin２θφ φ 　 　 　  1-1-41 　　221.2 1.2 导热方程的有限差分解法导热方程的有限差分解法•求解不稳定导热偏微分方程的数值解法，主要有：有限差分解法、有限元法、边界元法三类。

边界元法正在研究和完善之中，目前常用的是有限差分解法和有限元法•我们专门讨论有限差分解法的数学基础，数值方程的建立，差分方程的稳定性和收敛性等问题•有限差分解法是用差分方程近似地代替微分方程，建立差分方程有直接法和能量平衡法两种1.2.11.2.1直接代换法直接代换法直接代换法直接代换法•直接代换法是从微分形式出发，用差商直接代换微商的办法建立差分方程1.2.1.1数学基础　①微商和差商的定义•若T(x)是连续函数，则其导数为：　dT　       Ｔ(x+Δx)-T(x)              ΔＴ　— = lim   ——————  =  lim   ——                                                      1-2-1　dx    Δx→0        Δx　　　 Δx→0  Δx231.2.1.11.2.1.1数学基础数学基础•1-2-1式右边ΔＴ/Δ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)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处的各阶导数的线性组合，只能是近似地等于差商。

两者之间也必然有偏差241.2.1.11.2.1.1数学基础数学基础•图1-2-1表示了一阶差商与一阶微商之间的关系用不同方法得到的差商去代替微商，它们带来的误差是不同的即　　向前差商：dT/dx≈[T(x+ Δx)-T(x)]/ Δx                                                1-2-4　　向后差商：dT/dx≈[T(x)-T(x- Δx)]/ Δx                                                      1-2-5　　中心差商：dT/dx≈1/2{[T(x+ Δx)+T(x)]/ Δx+[T(x)-T(x- Δx)]/ Δx}                                       =[T(x+ Δx)-T(x- Δx)]/ 2Δx                                              1-2-6图1-2-1　差商与微商251.2.1.11.2.1.1数学基础数学基础•按照泰勒级数展开，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- Δ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) 261.2.1.11.2.1.1数学基础数学基础•由式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 　            2Δx                     dx　    3!     dx3 •即中心差商的截断误差为与(Δx)2同级的小量O[(Δx)2] 。

•当Δx固定时，上述一阶差商一般仍是x的函数，对它们还可以求差商这种一阶差商的差商称为二阶差商，它是二阶微商的近似，常用向前和向后差商来二阶微商，即　　d2T  　 T(x+ Δx)-T(x) 　T(x)-T(x - Δx)            　   —— ≈[  —————— - —————— ]/ Δx                                　    dx2       　　　Δx　　　　　Δx                     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] 。

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

291.2.1.2 1.2.1.2 建立内节点差分方程建立内节点差分方程•F0的数值小意味着加热或冷却此物体所需要的时间长，反之，所需要的时间短， F0是一个无因次数字•当初始条件和边界条件已知时，用式1-2-17就可模拟区域内各节点随时间τ增长的温度值Tpi(i=2,3,…，N-1；p=1,2,3,…）•隐式差分格式             •一维隐式差分如图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 )]       1-2-19                            Δτ                                      (Δx)2•若在上式中去掉O(Δτ )和O((Δx)2),整理得：      Tp+1i= Tpi+F0(Tp+1i+1-2Tp+1i+Tp+1i-1)                           1-2-20   301.2.1.2 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的联立线性代数方程组才能求解，也就是说，它含有三个未知数时间步长每前进一步，从坐标0≤x≤L，网格点1≤i≤N整个区域的每个点，上述方程都要列出一次（见图1-2-3)因此，向一个新的时间步长每移动一步就必须解一个方程组当按顺序列出这些方程时，除了要的第一点和最后一点的方程只有两个未知数外，其余每一个方程都含有三个未知数，于是方程是三对角的对于这种情况，已有适用的求解程序，后面将讨论关于差分方程稳定性将在后面讨论311.2.1.2 1.2.1.2 建立内节点差分方程建立内节点差分方程图1-2-3　一维隐式差分　②②二维系统二维系统•显式和隐式差分格式建立的方法和两种差分格式的特点在前面讨论过，对二维系统同样适用为简略起见，在此只讨论二维系统显式差分方程在此只讨论二维系统显式差分方程的建立的建立为使问题简化，仍然假定热物性值为常数，考虑无内热源•首先把所讨论的区域离散化，如图1-2-4所示×××××　×××　　×××　　　×××　　　　×××　　　　　×××　　　　　　×××　　　　　　　××　　　　　　321.2.1.2 1.2.1.2 建立内节点差分方程建立内节点差分方程　网线：用平行于X、Y轴的               直线（网线），进               行空间离散化　节点：网线与网线的交点　步长：节点间的距离　　　　　　　　　　　　　　图1-2-4   二维均质物体的分割•Δx与Δy可以是不变的常量，即等步长，也可以是变量（即在区域内的不同处是不同的），即变步长。

如果区域内各点处的温度梯度相差很大，则在温度变化剧烈处，网格布得密些，在温度变化不剧烈处，网格布得疏些至于网格多少，步长取多少为宜，要根据计算精度与计算工作量等因素而定•对于0

将上三式代入式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+1341.2.1.2 1.2.1.2 建立内节点差分方程建立内节点差分方程　③③三维系统三维系统•和二维系统的假定相同，热物性值为常数，无内热源设Δx=Δy=Δz，即均匀网格间距，将三维不稳定导热方程式应用于节点(i,j,k)，可以写成： T 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-6Tpi,j,k) 1-2-28 式中F0= λΔτ/[ρCp(Δx)2]。

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

在Δτ时间内，i 单元体的内能变化为：　 ΔU Tp+1i– Tpi —— = ρCp(FΔx)———— 1-2-31 　 Δ Δ 　　　•因为能量平衡或守恒，所以∑Q=ΔU/Δτ，即： Tp+1i– Tpi Tpi-1– Tpi Tpi+1– Tpi 　ρCp(FΔx)———— = λ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完全一样 371.2.2.1 1.2.2.1 建立内节点差分方程建立内节点差分方程②二维系统•如果物体的形状使之只在X、Y方向上有热流，Z方向上无热流，或相对于X、Y方向可忽略不计，这种情况就可以当作二维系统处理。

对于如图1-2-6所示的二维矩形区域，划分成有限个小单元体图中虚线所示的(i,j)(i,j)单元单元体体及周围四个相邻的单元体，用于推导内单元体差分方程每个小单元体的几何步长为Δx和Δy；时间步长Δτ为使问题简化，仍然假定热物性值为常数，考虑无内热源图1-2-6 二维系统的分割381.2.2.1 1.2.2.1 建立内节点差分方程建立内节点差分方程•对于(i,j)单元体，Δτ时间间隔内的内能变化为：　 ΔU Tp+1i,j– Tpi,j       —— = ρCp(ΔxΔy.1)————— 1-2-34                      　 Δ Δ•在Δτ时间内从周围四个相邻单元体流入(i,j)单元体的热量分别为：　　　　　 Tpi-1,j– Tpi,j  　Qi-1,j→i,j = λ(Δy.1)(——————) 1-2-35       　 Δx  Tpi+1,j– Tpi,j  　Qi+1,j→i,j = λ(Δy.1)(——————) 1-2-36      　 Δx  Tpi,j-1– Tpi,j  　Qi,j-1→i,j = λ(Δx.1)(——————) 1-2-37     　 Δy  Tpi,j+1– Tpi,j  　Qi,j+1→i,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)一样的节点差分方程　　391.2.2.2 1.2.2.2 能量平衡法与直接代换法的比较能量平衡法与直接代换法的比较1.2.2.2 1.2.2.2 1.2.2.2 1.2.2.2 能量平衡法与直接代换法的比较能量平衡法与直接代换法的比较能量平衡法与直接代换法的比较能量平衡法与直接代换法的比较•尽管用能量平衡法与直接代换法所得到的差分方程在形式上完全一样，但它们所依赖的基础有所不同用直接代换法直接代换法得出差分方程时，要求区域内温度分布连续而光滑，因为按泰勒级数展开，要求温度T(x)在x处的各阶导数都存在，连续而光滑即可导•而能量平衡法能量平衡法只要求每个单元体边界上热流可积，即可以算出每个单元体边界上的热量，它比温度分布连续光滑的要求显然放宽了•如果区域由导热系数截然不同的两种物质构成，并在交界面上存在接触热阻，或在交界面上存在间隙，那么在那些交界面上或间隙处，温度分布就不连续或不光滑在这种情况下，可用能量平衡法建立差分方程，只需在单元划分时，使那些交界面正好落在单元的边界即可，如图1-2-7所示。

　　图1-2-7 异种物质接触时单元的划分　图1-2-8 界面温度分布示意图401.2.2.2 1.2.2.2 能量平衡法与直接代换法的比较能量平衡法与直接代换法的比较•现假定二维矩形分割，交界面两侧的物体导热系数不同，讨论相邻两节点(i+1,j)节点流出(i,j)节点的热量公式的建立当异种物质接触时，除传导传热外，还要考虑界面热阻Rb，则　　 Tpi+1,j– Tp2  　Qi+1,j→T2 = λi+1,j(Δy.1)(—————) 1-2-40       Δx/2  Tp2– Tp1  　QT2→T1 = (Δy.1)(————) 1-2-41      Rb  Tp1– Tpi,j  　QT1→i,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/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代入上式得： 411.2.2.2 1.2.2.2 能量平衡法与直接代换法的比较能量平衡法与直接代换法的比较　      　　　　　 Δy.1 　　　　　　　　　　　 　Q= ————————————— (Tpi+1,j– Tpi,j ) 　 　　　 1-2-46 Rb+Δx/(2λi+1,j)+Δx/(2λi,j) •与直接代换法相比，能量平衡法的主要优点是：在二维系统中能使用三角形、四角形，在三维系统中能使用四、五、六面体等多种分割单元，因此能解复杂形状的问题；对于每个节点单元能指定热物性值，易于求解由多种物质组成的系统的问题。

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

Tpi+1,j–Tpi,j Tpi,j-1– Tpi,j λ(Δy .1)(—————)+λ(Δx/2 .1)(—————)+       Δx                                           Δy Tpi,j+1–Tpi,j Tp+1i,j– Tpi,j λ(Δx/2 .1)(—————)+0 =ρCp(Δx/2Δy.1)————— 1-2-47   Δy                                                     Δ　 431.2.3.11.2.3.1绝热边界绝热边界•若Δx=Δy，经整理后得：   Tp+1i,j= Tpi,j+F0(2Tpi+1,j + Tpi,j-1 + Tpi,j+1 - 4Tpi,j) 1-2-48 式中F0= λΔτ/[ρCp(Δx)2]。

1.2.3.21.2.3.21.2.3.21.2.3.2给定热流密度给定热流密度给定热流密度给定热流密度qrqrqrqr边界边界边界边界•把图1-2-9中BC面上任一节点(i,j)及其相邻节点取出分析，如图1-2-11所示用能量平衡法建立给定热流密度qr边界节点方程 Tpi-1,j–Tpi,j λ(Δy/2 .1)(—————)+ Δx Tpi+1,j–Tpi,j λ(Δy/2 .1)(—————)+                                   Δx Tpi,j+1–Tpi,j λ(Δx .1)(—————)+ qr(Δx.1) Δy Tp+1i,j– Tpi,j    =ρCp(ΔxΔy/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)+2qrΔτ/(ρCpΔx) 1-2-50 式中F0= λΔτ/[ρCp(Δx)2]。

　441.2.3.3 1.2.3.3 对流边界对流边界1.2.3.3 1.2.3.3 1.2.3.3 1.2.3.3 对流边界对流边界对流边界对流边界•把图1-2-9中CD面上任一节点(i,j)及其相邻节点取出分析，如图1-2-12所示在CD面上，ac为对流换热系数，Tf为物质周围介质温度用能量平衡法建立对流边界节点方程 Tpi-1,j–Tpi,j λ(Δy .1)(—————)+ Δx Tpi,j-1–Tpi,j λ(Δx/2 .1)(—————)+                                   Δy Tpi,j+1–Tpi,j λ(Δx/2 .1)(—————)+ac(Δy.1)(Tf- Tpi,j) Δy Tp+1i,j– Tpi,j    =ρCp(Δx/2Δy .1)————— 1-2-51  Δ　　　　　　　　　　　　　　　　　　　　　　•若Δx=Δy，经整理后得： 　　　　　　　　　　 图1-2-12对流边界单元   Tp+1i,j= Tpi,j+F0(2Tpi-1,j+Tpi,j-1+Tpi,j+1-4Tpi,j)+2acΔτ(Tf -Tpi,j )/(ρCpΔx) 　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　 1-2-52 式中F0= λΔτ/[ρCp(Δx)2]。

　1.2.3.4 1.2.3.4 1.2.3.4 1.2.3.4 给定温度边界给定温度边界给定温度边界给定温度边界•图1-2-9中AD面上给定温度Tw，把AD面上任一节点(i,j)及其相邻节点取出分析，如图1-2-13所示对这些单元直接写出结果：　　Tp+1i,j = Tw 1-2-53451.2.3.5 1.2.3.5 辐射边界辐射边界•把图1-2-9中AD面上受辐射热流的作用，把AD面上任一节点(i,j)及其相邻节点取出分析，如图1-2-13所示ε为黑度ε≤1，σ０为斯蒂芬—玻耳兹曼常数，Tf为物质周围介质温度 Tpi-1,j–Tpi,j λ(Δy/2 .1)(—————)+ Δx Tpi+1,j–Tpi,j λ(Δy/2 .1)(—————)+ Δx Tpi,j-1–Tpi,j λ(Δx .1)(—————)+ Δy εσ０(Δx.1)[Tf4 – (Tpi,j)4] Tp+1i,j– Tpi,j =ρCp(ΔxΔy/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)+ 2εσ０Δτ[Tf4 – (Tpi,j)4]/(ρCpΔx) 1-2-55 式中F0= λΔτ/[ρCp(Δx)2]。

　461.2.3.6 1.2.3.6 混合边界混合边界•图1-2-9中四个角上的边界单元，除AD边定温边界条件时，∠A、∠D 两个拐角节点的温度为一定值（Tw）的情况以外，在其余换热边界条件下，它们的两个外边界从属于两种边界条件，如∠C拐角节点，从图1-２-9可知，它的一个边界处于对流换热，另一个边界处于给定热流密度qr在这种情况下，将∠C拐角节点及其相邻节点取出分析，如图1-２-14所示，按能量平衡法可得：　　　　　　　 Tpi-1,j–Tpi,j λ(Δy/2 .1)(—————)+ Δx Tpi,j+1–Tpi,j λ(Δx/2 .1)(—————)+ Δy qr(Δx/2 .1)+ ac(Δy/2.1)(Tf – Tpi,j) Tp+1i,j– Tpi,j =ρCp(Δx/2Δy/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(Tf–Tpi,j)/λ 式中F0= λΔτ/[ρCp(Δx)2]。

1-2-57471.2.3.6 1.2.3.6 混合边界混合边界•对图1-2-9中∠A、 ∠ B、∠D等拐角节点，在绝热、给定热流密度qr、对流和辐射换热等混合边界条件下，都可以按上述类似方法，分别建立它们的边界单元差分方程•对于物体出现内拐点（凹角）的情况，也可按能量守恒定律得出凹角节点的差分方程如图1-2-15所示的凹角节点（i,j），若凹角两边处于对流换热边界条件时，把（i,j）节点及其相邻节点取出分析，则可得：　　　　　　 Tpi-1,j–Tpi,j λ(Δy .1)(—————)+ Δx Tpi+1,j–Tpi,j λ(Δy/2 .1)(—————)+ Δx Tpi,j-1–Tpi,j λ(Δx/2 .1)(—————)+ Δy Tpi,j+1–Tpi,j λ(Δx .1)(—————)+ Δy ac[(Δx/2+ Δy/2).1](Tf – Tpi,j) Tp+1i,j– Tpi,j =ρCp(3/4ΔxΔy .1)————— 1-2-58 图1-2-15 凹角边界单元 Δ　　　481.2.41.2.4圆柱坐标系的有限差分方程圆柱坐标系的有限差分方程•圆柱坐标系更容易表示圆柱体工件的形状，如果物体无内热源、热物性值为常数，圆柱坐标系的一般导热方程式如下： T 1  T 1  Ｔ  T　 ρCp — =λ[— —(r — ) + — —( — ) + —( — )]  　 r r r r2 φ φ z z 1 r T 1  T 1  Ｔ  T　 =λ{—(— — )+—[—(— )r]+— —( — ) + —( — )} 　 r r r r r r r2 φ φ z z 1 T 2T 1 2Ｔ ２T　 =λ( — — + — + — —— + —— ) 1-2-59 　 r r r2 r2 φ2 z２•上式的有限差分方程，可以用直接代换法和能量平衡法导出。

①①直接代换法直接代换法•图1-2-16表示圆柱坐标的几何图形，z、r、φ分别具有相等的增量，对于三维系统来说，直接代替偏微分方程之后，则可以得到节点０的显式和隐式的有限差分方程　　 Tp+10-Tp0 　 Tp2-Tp4 Tp2-2Tp0+Tp4 Tp1-2Tp0+Tp3 Tp5-2Tp0+Tp6 　ρCp ———=λ[——— +————— +————— +—————] 1-2-60 　Δ 　 r0２Δr (Δr)2 r02(Δφ)2 　(Δz)２•如果上式右端各项都以新的时间层p+1时刻计算，则得到隐式的差分方程491.2.4 1.2.4 圆柱坐标系的有限差分方程圆柱坐标系的有限差分方程　　　　　　　　　　　　图1-2-16　圆柱坐标②能量平衡法•在Δτ时间内，从周围相邻各点流入０点的热量为：　　　　　　 Tp1–Tp0 　Q1→0=λ(ΔzΔr)(———) 1-2-61 　 r0Δφ　　　　　　 Tp3–Tp0 　Q3→0=λ(ΔzΔr)(———) 1-2-62 　 r0Δφ　　　　　　 Tp2–Tp0 　Q2→0=λ[(r0+Δr/2)ΔφΔz](———) 1-2-63 　 Δr501.2.4 1.2.4 圆柱坐标系的有限差分方程圆柱坐标系的有限差分方程　　　　 Tp4–Tp0 Q4→0=λ[(r0-Δr/2)ΔφΔz](———) 1-2-64 　 Δr　　　　 Tp5–Tp0 Q5→0=λ[(r0ΔφΔr](———) 1-2-65 　 Δz　　　　 Tp6–Tp0 Q6→0=λ[(r0ΔφΔr](———) 1-2-66 　 Δz ΔU Tp+10– Tp0       ∑Q= —— = ρCp(r0ΔφΔrΔz)————— 1-2-67                      　 Δ Δ•经整理后，可以得出节点0的显式有限差分方程：　　 Tp+10-Tp0 　 Tp2-Tp4 Tp2-2Tp0+Tp4 Tp1-2Tp0+Tp3 Tp5-2Tp0+Tp6 　ρCp ———=λ[——— +————— +————— +—————] 1-2-68 　Δ 　 r0２Δr (Δr)2 r02(Δφ)2 　(Δz)２•式1-2-68和式1-2-60完全一样。

除了对称轴上的一点外，该式对其它任一节点都是适用的因为在轴上r0=0,分母为零•用能量平衡法导出对称轴上的点的有限差分方程，如图1-2-17所示•从２点传入到０点的热量　　　　 Tp2–Tp0 Q2→0=λ2(Δr/2)πΔz/N(———) 1-2-69 　 Δr•从５、６点传入到０的热量511.2.4 1.2.4 圆柱坐标系的有限差分方程圆柱坐标系的有限差分方程注：半径Δr/2        上下点为5、6点，节距Δz       ０点周围温度为T2  （或平均温度）　　设φ角方向分成N份　　　　　　　　　　　　　　　　　　图1-2-17对称轴上点　　　 Tp5–Tp0 Q5→0=λ(Δr/2)2π/N(———) 1-2-70 　 Δz　　　 Tp6–Tp0 Q6→0=λ(Δr/2)2π/N(———) 1-2-71 　 Δz ΔU Tp+10– Tp0       ∑Q= —— = ρCp(Δr/2)2 π Δz/N ————— 1-2-72                      　 Δ Δ•整理后得：　　 Tp+10-Tp0 　 Tp2-Tp0 Tp5-2Tp0+Tp6 　ρCp ———=λ[4 ——— + ————— ] 1-2-73 　Δ 　 (Δr)2 (Δz)2 521.2.5 1.2.5 差分方程的稳定性和收敛性差分方程的稳定性和收敛性•前面介绍了数值模拟中常用的两种差分格式，显式和隐式两种差分格式，显式和隐式，并讨论了不同情况下差分方程的建立，这些差分方程从形式上看都是可以进行计算的，但还有两个问题尚需要研究，那就是：①当空间步长Δx→０时和时间步长Δτ→０时，差分方程的解是否逼近偏微分方程的解，即差分方程是否收差分方程是否收敛？敛？②差分方程在计算过程中所产生的误差是不断增大，还是可以控制，即差分方程是否稳定？差分方程是否稳定？•数学上可以证明二者在适当的条件下，从稳定性可以推出收敛性，即稳定稳定是收敛的充分必要条件是收敛的充分必要条件，因此，只要讨论稳定性问题。

•保证解的稳定性在实际计算中是十分重要的因为误差存在于：初始条件、边界条件和数值计算的舍入误差初始条件与边界条件数据一般是实际测量的数据，其存在测量误差存在的误差如果在计算过程中不断被放大就会导致解的不稳定，则计算出的结果变的无意义•下面针对一维系统的两种差分格式，从物理意义方面讨论它们稳定性的问题从中得出的有关结论，原则上对二维、三维系统的两种差分格式也是适用的•一维显式差分方程格式的差分方程为：　　Tp+1i= Tpi+F0(Tpi+1-2Tpi+Tpi-1)=F0Tpi+1+(1-2F0)Tpi+F0Tpi-1 1-2-74　　由上式可以看到：531.2.51.2.5差分方程的稳定性和收敛性差分方程的稳定性和收敛性•①i节点在p+1时刻的温度Tp+1i只受p时刻i-1,i,i+1三个节点温度Tpi+1，Tpi，Tpi-1的影响•② Tpi+1，Tpi，Tpi-1三者的系数之和为１•从上可说明，Tp+1i是Tpi+1，Tpi，Tpi-1三者的加权平均值，其中F0或(1-2F0)是一种加权平均的系数•由此可知，要使显式差分格式的计算结果符合物理意义，Tpi+1，Tpi，Tpi-1三者的系数均应不小于零。

由显式差分方程可见，Tpi+1，Tpi-1的系数显然大于零，而要求Tpi的系数不小于零，必然的结果是F0≤1/2•如果F0>1/2，由上式可见，当p时刻i节点温度Tpi越大，因为(1-2F0)为<0则，p+1时刻i节点的温度值Tp+1i越小，这显然是违反热力学第二定律的•热力学第二定律指出，一切自发过程，如高温到低温传热，气体由高压区扩散到低压区，都是不可逆的，一直进行到平衡为止如果F0≤1/2，则当p>0时，全部内节点的温度Tp+1i的值总是处于Tpi+1，Tpi，Tpi-1 三个值的最大值与最小值之间的某个中间数值这一事实显然是符合热力学第二定律的•为什么隐式差分方程无条件稳定？•将隐式差分方程稍加整理，可得：541.2.51.2.5差分方程的稳定性和收敛性差分方程的稳定性和收敛性　　Tpi= Tp+1i+F0(2Tp+1i-Tp+1i+1-Tp+1i-1) 1-2-75•假定p+1时刻，在区域内某一节点i处取得区域内最大温度值，即Tp+1i为这一时刻区域内的最大温度,Tp+1i> Tp+1i+1 ；Tp+1i>Tp+1i-1,式1-2-75等号右边括号内三项的代数和必然大于零，从而必然有Tpi>Tp+1i。

也就是说，在在p p时时刻，区域内的最大温度必然大于刻，区域内的最大温度必然大于p+1p+1时刻区域内的最大温度时刻区域内的最大温度依此类推，必然将最大温度值推到初始条件•假定p+1时刻，在区域内某一节点i处取得区域内最小温度值，即Tp+1i为这一时刻区域内的最小温度,Tp+1i< Tp+1i+1 ；Tp+1i

•根据这一原则，其它边界条件下边界节点差分方程的稳定性也可以相应导出对于圆柱坐标系、球坐标系的显式差分方程，也要按上述原则分别导出具体的稳定性条件•总之，用显式差分格式数值求解导热问题，离散化处理时，时间步长总之，用显式差分格式数值求解导热问题，离散化处理时，时间步长Δτ和和空间步长空间步长Δx的选取是受到稳定性条件制约的的选取是受到稳定性条件制约的一般的习惯是先定空间步长，再计算时间步长561.2.6 1.2.6 二维交替隐式格式二维交替隐式格式•显式差分格式的优点：每个节点方程能独立求解，整个计算过程十分简单•缺点：为满足稳定性条件，空间步长和时间步长的选取必须遵守下列关系：                                      (Δx)2ρCp     一维系统：Δτ≤—————     (F0 ≤1/2)　         　　　　　　　　2λ                                      (Δx)2ρCp     二维系统：Δτ≤—————     (F0 ≤1/4)　         　　　　　　　　4λ                                      (Δx)2ρCp      三维系统：Δτ≤—————     (F0 ≤1/6)　         　　　　　　　　6λ•例如， 一维系统， Δx=1mm，ρ=7840kg/m3，Cp=610J/(kg.℃)，    λ=30J/(s.m.℃)，求得：Δτ≤0.079,s。

即1h约5 万次•这就表明：为了满足精度要求，当网格划分较细时，时间步长大大减小，当维数愈高，在同样的空间步长条件下，时间步长也愈小因此，为了为了得到稳定而又比较精确的答案，就要进行大量的计算得到稳定而又比较精确的答案，就要进行大量的计算•为了消除稳定性时间步长的限制，可以采用隐式格式求解为了消除稳定性时间步长的限制，可以采用隐式格式求解但是，在每571.2.61.2.6二维交替隐式格式二维交替隐式格式     一个时间层上，对二维系统就要解一个五对角的N*M个未知量的代数方程组(N为x方向的网格数，M为y方向的网格数)，即移动一个时间步长，就必须解整个方程组，每一个方程都有五个未知数和一个对角占优元素这样一来，就没有解三对角方程组那样方便•因此，需要建立新的差分格式，希望这一差分格式能满足如下要求：①无条件稳定；②有合理精度的解；③所产生的代数方程组易于求解•交替方向隐式方法正是满足这些要求的一种差分格式交替方向隐式方法正是满足这些要求的一种差分格式简称IAD）1.2.6.11.2.6.1差分方程的建立差分方程的建立差分方程的建立差分方程的建立•交替方向隐式方法需要对给定时间步长推导出二组有限差分方程。

这些方程是显式和隐式项的混合，对于前半个时间步长，x方向各项是隐式形式，y方向各项是显式形式；而在下半个时间步长时则相反•若把二维导热微分方程式应用于节点(i,j)，可以写成如下形式：　 第一个1/2Δτ T p 2T p+1/2 2T p       ρCp (—) = λ[(——) +(——) ] 1-2-76                     i,j x2 i,j　 y2 i,j　 　 Tp+1/2i,j– Tpi,j Tp+1/2i+1,j-2Tp+1/2i,j+Tp+1/2i-1,j      ρCp(——————)=λ[—————————————+         Δτ/2                                       (Δx)2 Tpi,j+1-2Tpi,j+Tpi,j-1 ————————] 　　　　　　　　　　　　　　　　 1-2-77  (Δy)2581.2.6 1.2.6 二维交替隐式格式二维交替隐式格式      第二个1/2Δτ T p+1/2 2T p+1/2　 2T p+1       ρCp (——) = λ[(——) + (——) ] 1-2-78                    i,j x2 i,j　 y2 i,j　 　 Tp+1i,j– Tp+1/2i,j Tp+1/2i+1,j-2Tp+1/2i,j+Tp+1/2i-1,j      ρCp(———————)=λ[—————————————+         Δτ/2                                (Δx)2 Tp+1i,j+1-2Tp+1i,j+Tp+1i,j-1 ——————————] 　　　　　　　　　　　　　　　　 1-2-79 　 (Δy)2•从上可见：第p+1/2层是一个过渡时间层，通过式1-2-77计算过渡时间层的值，然后用式1-2-79计算第p+1时刻的值。

由于现在有一种迅速解三对角线矩阵的算法，所以交替方向隐式法是有效的1.2.71.2.7三对角线矩阵的解法三对角线矩阵的解法三对角线矩阵的解法三对角线矩阵的解法•前面介绍了显式、隐式和交替方向隐式法等差分格式，从数学方面来看，它们都是一个线性代数方程组线性代数方程组的求解一般有两类方法线性代数方程组的求解一般有两类方法：一类是直接方法直接方法，如三对角线矩阵的解法、高斯消去法、矩阵分解法；另一类是迭代方法迭代方法，如高斯—塞德尔迭代、超松弛迭代、强隐式迭代•一般说来，直接方法只要进行有限此运算就可以得到方程的解它的缺点是必须有较大的计算机存储量和花费较多的计算时间，而且在使用时，591.2.7 1.2.7 三对角线矩阵的解法三对角线矩阵的解法　 还必须注意抑制舍入误差的影响，可以采用增加计算数字的位数，或利用主元消去等手段•迭代法没有舍入误差积累的问题，一般需要较少的计算机存储量，计算时间也比较节省但是，如何能确保迭代收敛和收敛得比较快，这在使用中必须注意解决•当离散化方程是三节点关系式时，高斯消去法就转换成人们所常用的三对角线矩阵法，这种方法也称为托马斯(Thomas)算法或TDMA、追赶法。

它之所以称为三对角线矩阵算法，是由于用矩阵形式写出这些三节点关系式时，所有系数都是沿矩阵的三对角线排列的•设离散化方程是如下形式的三节点关系式：     -  aiTi-1  +  biTi  -  ciTi+1  =  di                                                                     1-2-80•i=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)以上二个方程三个未知数）601.2.7 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求得，这就是三对角线矩阵解法的实质。

•具体公式推导如下：　　设待求的递推替换关系式为：        T i= Pi  Ti+1  +  Qi                                                                                    1-2-81　　则Ti-1= Pi-1  Ti  +  Qi-1                                                                             1-2-82       将式1-2-82代入式1-2-80，得　　- ai( Pi-1 Ti  + Qi-1 ) + biTi  -  ciTi+1  =  di ，经整理得：                                                                    　         Ci                          di+aiQi-1            Ti= ———  Ti+1 + ———                                                                     1-2-83              bi-aiPi-1                    bi-aiPi-1 •比较式1-2-83和式1-2-81得　　　 　Ci                                di+aiQi-1            Pi= ———  ， Qi = ———                                                                  1-2-84              bi-aiPi-1                         bi-aiPi-1 •式1-2-84是个递推关系式，因为它们是用Pi-1，和Qi-1来表示Pi，和Qi的。

611.2.71.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.11.2.7.1三对角线矩阵解法程序三对角线矩阵解法程序三对角线矩阵解法程序三对角线矩阵解法程序①Turbo BASIC语言•设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)     DIM A(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!     FOR I= 2 TO N - 1       A(I) = f0: B(I) = 1 + 2 * f0: C(I) = f0: D(I) = 15!     NEXT I621.2.71.2.7三对角线矩阵的解法三对角线矩阵的解法    A(N) = 0: B(N) = 1: C(N) = 0: D(N) = 100!    CLS: PRINT " dx="; dx * 1000; "(mm) "    PRINT " dt="; dt / 60; "min"    PRINT " Temperature Results"    PRINT "min  1P    2P     3P     4P     5P     6P     7P     8P     9P    10P   11P"      FOR J = 1 TO M         CALL TDMA          PRINT J * dt / 60;       FOR I = 1 TO N           PRINT INT(TMID(I) * 100) / 100;        NEXT I          PRINT       FOR I = 2 TO N - 1          D(I) = TMID(I)        NEXT I     NEXT J     END     631.2.71.2.7三对角线矩阵的解法三对角线矩阵的解法•REM　解三对角线矩阵子程序SUB TDMA     SHARED TMID(),A(),B(),C(),D(),N     LOCAL BAP(),DAQ()     DIM BAP(N),DAQ(N)     BAP(1)=B(1):DAQ(1)=D(1)     FOR I=1 TO N-1        BAP(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)     NEXT I        TMID(N)=DAQ(N)/BAP(N)     FOR I=N-1 TO 1 STEP –1        TMID(I)=(DAQ(I)+C(I)*TMID(I+1))/BAP(I)     NEXT IEND SUB641.2.71.2.7三对角线矩阵的解法三对角线矩阵的解法②Quick BASIC语言DECLARE SUB TDMA ()REM  simul     N = 11: M = 19: s = .35: lamda = 30!: cp = 610!: rou = 7840!: dt = 120!     dx = s / (N - 1): f0 = lamda * dt / (rou * cp * dx * dx)     DIM A(1 TO N), B(1 TO N), C(1 TO N), D(1 TO N), TMID(1 TO N)     A(1) = 0: B(1) = 1: C(1) = 0: D(1) = 100!     FOR i = 2 TO N - 1     A(i) = f0: B(i) = 1 + 2 * f0: C(i) = f0: D(i) = 15!     NEXT i     A(N) = 0: B(N) = 1: C(N) = 0: D(N) = 100!     CLS     PRINT " dx="; dx * 1000; "(mm) "     PRINT " dt="; dt / 60; "min"     PRINT " Temperature Results"     PRINT "min  1P    2P     3P     4P     5P     6P     7P     8P     9P    10P   11P"     651.2.71.2.7三对角线矩阵的解法三对角线矩阵的解法FOR j = 1 TO MCALL TDMA     PRINT j * dt / 60;     FOR i = 1 TO N     PRINT INT(TMID(i) * 100) / 100;     NEXT i     PRINT     FOR i = 2 TO N - 1     D(i) = TMID(i)     NEXT iNEXT jEND SUB TDMA     SHARED TMID(), A(), B(), C(), D(), N     DIM BAP(1 TO N), DAQ(1 TO N)     661.2.71.2.7三对角线矩阵的解法三对角线矩阵的解法BAP(1) = B(1): DAQ(1) = D(1)FOR i = 1 TO N - 1     BAP(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)NEXT i     TMID(N) = DAQ(N) / BAP(N) FOR i = N - 1 TO 1 STEP -1     TMID(i) = (DAQ(i) + C(i) * TMID(i + 1)) / BAP(i) NEXT iEND SUB程序解释：　①BAP(I)=B(I)-A(I)*P(I-1)，其是式1-2-84中bi-aiPi-1     ②DAQ(I)=D(I)+A(I)*Q(I-1)，其是式1-2-84中di+aiQi-1     ③BAP(1)=B(1)-A(1)*P(0)=B(1)，因为P0=0    ④DAQ(1)=D(1)+A(1)*Q(0)=D(1)，因为Q0=0    ⑤BAP(2)=B(2)-A(2)*C(1)/BAP(1)=B(2)-A(2)*P(1)，推广到I=1 TO N-1    671.2.71.2.7三对角线矩阵的解法三对角线矩阵的解法   ⑥DAQ(2)=D(2)+A(2)*DAQ(1)/BAP(1)=D(2)+A(2)*Q(1)，推广到I=1 TO N-1                                                          D(N)+A(N)*DAQ(N-1)/BAP(N-1)    ⑦TMID(N)=DAQ(N)/BAP(N)=——————————————=Q(N)                                                         　B(N)-A(N)*C(N-1)/BAP(N-1)        (TN=QN)    ⑧TMID(I)=(DAQ(I)+C(I)*TMID(I+1))/BAP(I)    　　　　　D(I)+A(I)*Q(I-1)+ C(I)*TMID(I+1)                      =———————————————                                     B(I)-A(I)*P(I-1)                  Ci                            di+aiQi-1            (Ti= ———  Ti+1 + ———    1-2-83)               bi-aiPi-1                     bi-aiPi-1 1.2.7.21.2.7.2将数值模拟问题变成三对角线矩阵解法形式将数值模拟问题变成三对角线矩阵解法形式将数值模拟问题变成三对角线矩阵解法形式将数值模拟问题变成三对角线矩阵解法形式①一维隐式差分格式　式1-2-20 Tp+1i= Tpi+F0(Tp+1i+1-2Tp+1i+Tp+1i-1),整理成：　-F0 Tp+1i-1 + (1+2F0)Tp+1i – F0 Tp+1i+1 = Tpi 所以 ai=F0 bi=1+2F0 ci=F0 di= Tpi     681.2.71.2.7三对角线矩阵的解法三对角线矩阵的解法②二维交替隐式格式•第一个1/2Δτ 　　　　 Tp+1/2i,j– Tpi,j Tp+1/2i+1,j-2Tp+1/2i,j+Tp+1/2i-1,j       式1-2-77 ρCp(——————)=λ[—————————————+         Δτ/2                                       (Δx)2 Tpi,j+1-2Tpi,j+Tpi,j-1 ————————] (Δy)2•设Δx=Δy，F0=λΔτ/[2ρCp(Δx)2], 整理成: -F0Tp+1/2i-1,j+(1+2F0)Tp+1/2i,j-F0Tp+1/2i+1,j=Tpi,j+F0(Tpi,j+1-2Tpi,j+Tpi,j-1) 所以ai=F0 bi=1+2F0 ci=F0 di=Tpi,j+F0(Tpi,j+1-2Tpi,j+Tpi,j-1)     •第二个1/2Δτ 　　　　 Tp+1i,j–Tp+1/2i,j Tp+1/2i+1,j-2Tp+1/2i,j+Tp+1/2i-1,j       式1-2-79 ρCp(——————)=λ[—————————————+         Δτ/2 Tp+1i,j+1-2Tp+1i,j+Tp+1i,j-1 ——————————] 　 (Δy)2691.2.71.2.7三对角线矩阵的解法三对角线矩阵的解法设Δx=Δy，F0=λΔτ/[2ρCp(Δx)2], 整理成: -F0Tp+1i,j-1+(1+2F0)Tp+1i,j-F0Tp+1i,j+1=Tp+1/2i,j+F0(Tp+1/2i+1,j-2Tp+1/2i,j+Tp+1/2i-1,j) 所以aj=F0 bj=1+2F0 cj=F0 dj=Tp+1/2i,j+F0(Tp+1/2i+1,j-2Tp+1/2i,j+Tp+1/2i-1,j)701.2.71.2.7三对角线矩阵的解法三对角线矩阵的解法1.2.7.31.2.7.3三对角线矩阵解法三对角线矩阵解法三对角线矩阵解法三对角线矩阵解法Ｃ语言Ｃ语言程序程序程序程序设A[i]为ai系数，B[i]为bi系数，C[i]为Ci系数，D[i]为di系数，N为N个温度变量，TMID[i]为Ti温度，M为计算M个时间步长 。

include#define N 11#define M 12#define S 0.35void TDMA(void);float A[N+1],B[N+1],C[N+1],D[N+1],TMID[N+1];float lamda=30.0,cp=610.0,rou=7840.0;float dx=S/(N-1),dt;float F0;/*main function*/main(){ int i,p;  float F0;  clrscr();  711.2.71.2.7三对角线矩阵的解法三对角线矩阵的解法dt=120.0;  F0=lamda*dt/(rou*cp*dx*dx);  printf(" dx= %5.2f (mm)\n",dx*1000);  printf(" dt= %5.2f (min)\n",(dt/60));  printf(" Temperature Results(C)\n");  for(i=0;i<=N;i++) printf("------");printf("\n");  printf("  min  ");  for(i=1;i<=N;i++) printf("%d P   ",i);printf("\n");  for(i=0;i<=N;i++) printf("------");printf("\n");  for(i=1;i<=N;i++){A[i]=F0;B[i]=2*F0+1;C[i]=F0;D[i]=15.0;}  A[1]=0;B[1]=1.0;C[1]=0;D[1]=100.0;A[N]=0;B[N]=1;C[N]=0;D[N]=100.0;  for(p=1;p<=M;p++){ TDMA();  printf("% 5.2f ",(dt*p/60));  for(i=1;i<=N;i++) {printf("%5.2f ",TMID[i]);}  for(i=1;i<=N;i++) D[i]=TMID[i]; printf("\n");}721.2.71.2.7三对角线矩阵的解法三对角线矩阵的解法getch();}void TDMA (void){float BAP[N+1],DAQ[N+1];  int i;  BAP[1]=B[1];DAQ[1]=D[1];  for (i=1;i<=N-1;i++){ BAP[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];}  TMID[N]=DAQ[N]/BAP[N];  for (i=N-1;i>=1;i--)  TMID[i]=(DAQ[i]+C[i]*TMID[i+1])/BAP[i];}731.2.71.2.7三对角线矩阵的解法三对角线矩阵的解法•运行结果 dx=  35.0 (mm) dt=   2.0   (min) Temperature Results(℃) -------------------------------------------------------------------------------------------------------------  min       1 P        2 P       3 P       4 P      5 P       6 P       7 P       8 P       9 P      10 P      11 P-------------------------------------------------------------------------------------------------------------- 2.00      100.00   40.55   22.69   17.33   15.76   15.42   15.76   17.33   22.69   40.55    100.00 4.00      100.00   54.30   30.97   21.12   17.45   16.54   17.45   21.12   30.97   54.30    100.00  6.00      100.00   62.44   38.13   25.48   19.93   18.41   19.93   25.48   38.13   62.44    100.00  8.00      100.00   67.72   44.02   29.89   22.95   20.91   22.95   29.89   44.02   67.72    100.00 10.00     100.00   71.41   48.84   34.12   26.28   23.87   26.28   34.12   48.84   71.41    100.00 12.00     100.00   74.18   52.86   38.09   29.77   27.12   29.77   38.09   52.86   74.18    100.00 14.00     100.00   76.37   56.29   41.78   33.29   30.52   33.29   41.78   56.29   76.37    100.00 16.00     100.00   78.17   59.27   45.22   36.76   33.96   36.76   45.22   59.27   78.17    100.00 18.00     100.00   79.71   61.91   48.42   40.15   37.37   40.15   48.42   61.91   79.71    100.00 20.00     100.00   81.05   64.30   51.42   43.41   40.70   43.41   51.42   64.30   81.05    100.00 22.00     100.00   82.26   66.47   54.22   46.53   43.91   46.53   54.22   66.47   82.26    100.0024.00     100.00   83.35   68.47   56.85   49.50   46.99   49.50   56.85   68.47   83.35    100.00--------------------------------------------------------------------------------------------------------------741.3 1.3 数值模拟的求解条件数值模拟的求解条件•为了使工件内的温度场被唯一地确定下来，必须把一些复杂的因素考虑到数值方程的求解条件中去。

如何确定初始条件、边界条件、热物参数随各种因素的变化、潜热等1.3.11.3.1初始条件初始条件初始条件初始条件•A-B两物质组成的一维系统，其离散化处理如图1-3-1所示图1-3-1　A-B两物质系统离散化处理•所谓初始条件，即计算时间的初始起点，τ=0时，模拟系统的条件，如Ａ物质内部点温度处处为TA0，B物质内部点温度处处为TB0，TA0和TB0也可以是x、y、z坐标的函数•问题是：在τ=0时，A-B两物质交界面温度Tb0确定为何值？•关于Tb0的确定，有以下两种方法•第一种是实测法，进行实测统计获得Tb0的普遍规律•第二种是解析法，即认为Ａ-B两物质接触的瞬间分界面上温度保持不变，751.3 1.3 数值模拟的求解条件数值模拟的求解条件　 且只有与分界面垂直方向有热交换根据一维不稳定导热原理，若物质为均匀的，且其热物性值不随时间（温度）而变化，于是可得A-B两物质交界面的温度T0inf为：     (T0A- T0inf)bA=(T0inf- T0B)bB   　，整理得 　　　　　　　　　　　　　　　　　　　1-3-1　　　　　　bA T0A+ bB T0B      T0inf = ——————　　　　　　　　　　　　　　　　　　　　1-3-2                    bA+ bB       式中bA、bB 为Ａ、Ｂ物质的蓄热系数，b=√λρCp•若分界面不在节点上，如图1-3-2所示。

如图1-3-2　分界面不在节点上761.3 1.3 数值模拟的求解条件数值模拟的求解条件•先求出T0inf　，然后分别求出分界面Ａ、Ｂ物质一侧的温度，即：                                                               x     T0b (A/B)=T0inf+(T0(A/B)-T0inf)erf( ——— )                                                1-3-3　　　　　　　　　　　　　　   2√αt•式中，x为界面两侧节点离界面的距离，一般为Δx/2；　　　　α为Ａ或Ｂ物质的导温系数， α=λ/(ρCp)；　　　　 t为形成A-B界面的时间；　　　　erf()为误差函数，有误差函数表可查得•例如，对350×350×35mm的铝合金板件砂型铸造，浇注时间为4s，沿板厚方向一维数值模拟，选取Δx=2.5mm，设T0A=640℃，λA=192.6w/(m ℃)，ρA=2680kg/m3，(Cp)A=1.34kJ/(kg ℃) ，T0B=15℃，λB=0.565w/(m ℃)，ρB=1520kg/m3，(Cp)A=0.95kJ/(kg ℃) 。

•根据式1-3-2可计算出：　 T0inf  ≈619 ℃•根据式1-3-3可计算出：     T0b (A) ≈621 ℃       T0b (B) ≈304 ℃ 771.3 1.3 数值模拟的求解条件数值模拟的求解条件1.3.21.3.2边界条件边界条件边界条件边界条件①A-B两物质界面•A-B两物质界面分为理想接触和非理想接触两种情况所谓理想接触，就是A-B两物质之间相互紧贴，相互接触的两表面具有相同的温度，如图1-3-3所示图1-3-3理想接触•在这种情况下，若导热系数不随温度而变化，如图1-3-2所示的网格分割（界面无气隙）若节点之间没有热阻，T0b(A侧)→ T0b(B侧)可按串联热阻叠加特性求出：      Δx        Δx/2     Δx/2                           2λAλB　 ——　= ——　+ ——            λA→B   = ———                                  1-3-4    λA→B            λA              λB 　　　　　　　　　　λA+ λB781.3 1.3 数值模拟的求解条件数值模拟的求解条件•实际上，Ａ-B两物质界面常是属于非理想接触。

如钢液在钢锭模内冷却凝固时，钢锭的体积一般要逐渐收缩，而钢锭模受热要逐渐膨胀此时，在钢锭和钢锭模之间形成一层很薄的气体层，这一气体层和钢锭模内表面的涂料层合在一起称为气隙，如图1-3-4所示这样，在钢锭和钢锭模之间就要产生界面温度差ΔTm，如图1-3-5所示　　　　图1-3-4界面气隙　　　　　　　图1-3-5因气隙产生的界面温度差•在气隙的涂料层中，热量传递是靠导热，而在气隙的气体层内同时存在着导热、辐射和对流换热现象气隙中的物质在受热时能蓄热，在冷却时能放热因此，气隙层中的传热过程是很复杂的气隙的厚度一般是零至几791.3 1.3 数值模拟的求解条件数值模拟的求解条件     毫米，所以气隙的厚度同钢锭和钢锭模的尺寸相比，通常是很小的这意味着复杂形状的气隙可允许看作为传热学意义上的平壁气隙热阻Rb如图1-3-4可用下式计算：　　　　　　　　  δg　    δco　　Rb = Rg + Rco = ——  + ——                                                                  1-3-5                                   λg　    λco•式中，Rg 为气体层热阻；Rco为涂料层热阻；　　　　δg、δco分别为气体层和涂料层的厚度；　　　　λg、λco分别为气体层和涂料层的导热系数；•气隙中的传热系数h=1/Rb[w/(m2℃）]，钢锭和钢锭模间产生了气隙，通过传热系数h必然要影响钢锭的凝固过程。

为此，有关学者对此问题作了不少的研究，但目前还没有找到一种令人满意的方法来确定界面热交换系数h值的大小现将普遍用于研究和处理该问题的方法介绍如下，供参考•第一种方法是试算法第一种方法是试算法该法假定界面传热系数，在气隙形成前为h1，在气隙形成后为h2，以及形成气隙时钢锭表面温度为T２通过多次数值模拟计算，得出钢锭和钢锭模之间采用不同的h1,h2和T2时的计算凝固时间，并与实测的凝固时间相比较选择计算结果同实测结果相符者，计算时采用的h1和h2即为所求的传热系数这是一种较简化的方法801.3 1.3 数值模拟的求解条件数值模拟的求解条件•第二种方法是求解界面热流法第二种方法是求解界面热流法该法通过测量钢锭内部和钢锭模内部的温度场，并将其推算到界面，求出界面温度降（ΔT），并从界面的温度曲线斜率中计算出通过界面的热流通过界面的热流与界面温度降的比值，即h=q/ΔT便是界面热交换系数用此法能求出h-τ曲线，但应注意到，在推算界面温度曲线及从曲线斜率中求热流时，误差易放大•第三种方法是试凑法第三种方法是试凑法该法是在测量了钢锭和钢锭模内部几点温度后，通过控制│Tε-Tε’│<ε（ε为人为规定的误差），来求解界面热阻Rb(h=1/Rb)。

 如果Tε和Tε’分别为钢锭模内距界面10mm处的计算和测量温度值也就是说，设在一定的时间间隔内Rb是某一常数，将不同的Rb分别代入计算方程，直到满足│Tε-Tε’│<ε为止，然后再转入下一个时间间隔的计算这种方法可求出h-τ曲线，但在钢锭模内距界面10mm处的计算温度虽然较准确，而在其它点的计算温度误差较大•第四种方法是界面温度函数法第四种方法是界面温度函数法通过大量实测，获得一系列的界面温度差随时间变化值，再经过回归分析，得到界面温差函数ΔT=f(τ)将此函数代入差分方程式导出一运算温度，这样就可以绕过求解h值的问题，而且避免了钢锭模温度场的计算，直接计算出钢锭的温度场•第五种方法用界面温度法建立通用的钢锭第五种方法用界面温度法建立通用的钢锭-钢锭模边界方程钢锭模边界方程这种方法的基本思想是：建立无量纲关系，即相对温度和相对时间的关系假定浇注工艺的影响用界面上钢锭一侧的最高温度Tp来衡量，钢锭和钢锭模热物性811.3 1.3 数值模拟的求解条件数值模拟的求解条件　 值及其尺寸大小的影响表现为界面上钢锭一侧温度Ti的降低速度，而且这一影响可以某一估计时间参数来表示，称此估计的时间参数为计算凝固时间，记为tf　。

令界面最高温度Tp为基准温度，按当量厚度法得出的计算凝固时间tf为基准时间，则无量纲温度θ定义为：　　θ＝Ti / Tp                                                                                           1-3-6•无量纲时间τ定义为从浇注完开始到现时刻的时间t与基准时间tf的比，即       τ＝  t  /  tf                                                                                                                                                1-3-7•在不同工艺条件下进行浇注试验，实测界面上钢锭一侧温度Ti随时间 t 的变化根据实测数据，进行数学处理和无量纲转换，可以得到θ-τ关系式　 θ＝A0+A1τ+A2τ2+A3τ3+A4τ4+A5τ5                                            1-3-8     式中： A0、A1、A2、A3、A4、A5为常数。

•数值模拟时，将t=PΔt代入式1-3-7得到τ，再将τ代入式1-3-8得θ，将θ代入式1-3-6即可算得Ti，即钢锭一侧温度这样一来，便可抛开钢锭模，只对钢锭进行计算•界面温度法看来是一种比较理想的方法，但由于铸造过程的复杂性，这种方法所得结果实用的可靠性还有待进一步实验检验针对不同工艺条件、合金种类和工件不同部位，分别展开试验研究，得出各自条件下的821.31.3数值模拟的求解条件数值模拟的求解条件　 具有一定通用性的边界条件的方程，是今后努力的方向企图通过具体实验，得出普遍适用的通用边界方程，显然是困难的②与大气的边界•如果数值模拟对象的表面温度较低，其与大气之间主要发生自然对流热交换这时，边界条件容易处理假定环境温度为Tf，自然对流换热系数为ac，网格点间距为Δx，模拟对象的导热系数为λ，如图1-3-8所示的边界节点的差分方程将如何建立a)边界节点为半个控制单元　　　　　 (b)边界节点为整个控制单元 图1-3-8　与大气的边界•边界节点为半个控制单元　　　　　　　 Tp+1i– Tpi Tpi-1– Tpi    ρCp(ＦΔx/2)———— =λF(————)+acF(Tf- Tpi) 1-3-9   Δ Δx •边界节点为整个控制单元831.3 1.3 数值模拟的求解条件数值模拟的求解条件　　　　　　 Tp+1i– Tpi Tpi-1– Tpi    ρCp(ＦΔx)———— =λF(————) + hi F(Tf- Tpi) 1-3-10   Δ Δx                             1     式中：hi = ————                       １    Δx/2                       — + ——                        ac λ•如果数值模拟对象的表面温度较高，其与大气之间发生自然对流和辐射热交换。

这时，模拟对象与大气的传热系数应采用下式计算     h = ac + ar                                                                                                                                                          1-3-11•此处的ac是指仅仅为自然对流换热时的换热系数[w/(m2k)]•对于垂直平板或者圆柱      ac=1.42×(ΔT/L)1/4                      (104〈Re〈109)层流　  ac=0.95×(ΔT)1/3                            (109〈Re )        紊流•对于水平圆柱      ac=1.32×(ΔT/d)1/4                       (104〈Re〈109)层流　  ac=1.24×(ΔT)1/3                            (109〈Re )        紊流•对于水平平板（上面）      ac=1.32×(ΔT/L)1/4                       (104〈Re〈109)层流841.3 1.3 数值模拟的求解条件数值模拟的求解条件      ac=1.43×(ΔT)1/3                            (109〈Re )        紊流•对于水平平板（下面）      ac=0.61×(ΔT/L)1/3                       (层流、 紊流)•上列式中ΔT为物体表面温度与大气温度之差，即ΔT=(Ti-Tf)，K      L为物体的长度，m     d为物体的直径，m　  Re为雷诺数•另外，ar是热辐射传热系数，可以用下式表示：     ar=热辐射比热流量/（Ti-Tf）= εσ0φ1,,2（Ti4-Tf4）/（ Ti-Tf）       = εσ0φ1,,2（Ti2+Tf2）（Ti2-Tf2） /（ Ti-Tf）       = εσ0φ1,,2（Ti2+Tf2）（Ti+Tf）                                                        1-3-12　式中：ε为物体的辐射率或黑度；　　　　σ０为斯蒂芬—玻耳兹曼常数；5.67×10-8 w/m2k4　　　　　　Ｔi为外表面温度；　　　　φ1,2 为角系数，因大空间小物体，这时φ1,2 =1。

851.3 1.3 数值模拟的求解条件数值模拟的求解条件1.3.3 1.3.3 热物性值热物性值热物性值热物性值•模拟对象的热物性值主要包括比热Cp，导热系数λ，密度ρ，潜热L其中除ρ可大致看作常数外， Cp、λ、L一般随温度而变化这种变化反映到数值模拟方程中，每一时间步长都可能包含很多个热物性值，其对模拟精度有影响实际上影响热物性值的因素很多，各家所得数据又很不一致因此，如何精确地掌握热物性值是数值模拟研究中一项很重要的工作•解决热物性值问题的途径主要有三个方面：　①查阅手册和资料；　②进行专门的实测；　③实验和数学处理相结合•热加工是在高温下进行的，因此需要高温下的热物性参数，但实际上高温的热物性参数值很难查到专门实测当然最为理想，可需要应用高温比热仪、导热仪和量热计，但能测1200-1500℃以上温度的这类设备，如日本生产的用激光法测量热物性值的TC-3000-HNC型（测试温度范围，室温—1500℃）和TC-3000-UVH型（最高温度达2200℃），这类设备，制造成本高，价格昂贵在可能条件下，采用实验与数学处理相结合的方法，但其精度如何有待进一步验证•对Cp、λ常采用以下方法处理：861.3 1.3 数值模拟的求解条件数值模拟的求解条件　①常数法     此种方法基于这样的假定，如果数值模拟的温度范围不大，数值模拟时取其温度范围内的某一温度下的Cp和λ值作为其平均值，其误差不大。

　②线性函数法　 假定Cp和λ以线性方式随温度变化，    Cp=Cp1+a1(T-T1)，λ=λ1+b1(T-T1) ，                                             1-3-13    式中Cp和λ为Ｔ温度下的比热和导热系数；　　   Cp1和λ1为T1温度下的比热和导热系数；　　　a1和b1为比例常数   ③插值法　用已知某几点温度下的Cp和λ值，根据数学上的插值公式，计算出每一个不同温度下的Cp和λ值•不同的含碳量的普碳钢，其Cp和λ不同•相同的含碳量的铁，如果其凝固成白口铁、灰口铁、球铁，其Cp和λ不同•Fe-C合金熔化潜热Ｌ随含碳量变化871.3 1.3 数值模拟的求解条件数值模拟的求解条件1.3.4 1.3.4 潜热处理潜热处理潜热处理潜热处理•当合金由液态转变为固态时，由于液态的内能EL大于固态的内能ES，必须产生ΔE= EL-ES的内能变化这个内能变化ΔE（通常用Ｌ表示）称为凝固潜热，或称为熔化潜热 1.3.4.1 1.3.4.1 考虑了析出潜热的热能守恒式考虑了析出潜热的热能守恒式考虑了析出潜热的热能守恒式考虑了析出潜热的热能守恒式•假定单位体积、单位时间内固相率的增加率为gS/τ，因此，潜热放出的热量为：ρＬgS/τ，考虑此发热量以后的热能守恒式，对于一维问题为：　　　 T  T gS ρCp — = — (λ— )+ ρＬ—                1-3-14　  X X τ•此处如果作如下变换：　 gS gS T ρＬ— = ρＬ— — 1-3-15　 τ T τ•把潜热项移到左端，则变成：　　　 gS T  T ρ（Cp-L — ）— = — (λ— )   1-3-16 T  X X •因此，如果固相率gS和温度T的关系已知，则式1-3-16就能够很容易地进行数值模拟。

881.3 1.3 数值模拟的求解条件数值模拟的求解条件•严格地说，质量固相率fS和体积固相率gS是不同的，但是以下也可以近似地认为： fS＝gS 1.3.4.2 1.3.4.2 固相率和温度的关系固相率和温度的关系固相率和温度的关系固相率和温度的关系•从相图可知道固相率（质量固相率fS）和温度的关系对于恒温下凝固的纯金属、共晶凝固和包晶凝固，其固相率不能根据温度来确定对于具有一定结晶温度范围的合金，固相结晶析出的固-液共存区中，液相线温度是与液相浓度相对应的对于一般合金的凝固，若假定液相线为直线，平衡分配系数k为常数的话，对于k<1的合金，则液相线温度和液相中的溶质浓度CL呈线性关系，即　 T = T0 – m CL                                                                                                                                           1-3-17     式中：T0 为溶剂的纯金属的熔点； m 为液相线斜率•对于含碳量0.18-0.59％的碳钢、低合金钢和中合金钢，当每种合金元素含量的极限为：5%Si、9%Mn、5%Cr、3%Ni、6%Cu、7%Mo、5%V、12%W、10%AL、15%Co、3%Ti、0.7%P、0.08%S、0.03%O2、0.03%N2、0.1%B和5%Zr的情况下，可利用下式计算液相线和固相线温度。

 •TL=1539-(78C+12Si+5Mn+1.5Cr+4Ni+5Cu+12Mo+2V+1.5W+2AL+2Co            +10Ti+150B+6Zr+30P+25S+80O2+90N2+1300H2)                         1-3-18•TS=TL-[13+140C+4.5Si+2.5Mn+Cr+4.5Ni+3(Cu+Mo+V+W+AL)+Co+5Ti]   1-3-19891.3 1.3 数值模拟的求解条件数值模拟的求解条件•因此，为了求液相线温度，必须知道固液共存区的溶质浓度，而溶质浓度区随固相率而变化，也就是说：要了解固相率和温度的关系，无非是了解固相率和溶质浓度的关系•对于实用的多元合金，往往不清楚其正确的状态图，所以分配系数也无法知道在这种情况下可以先采用热分析求出凝固开始温度TL和结束温度TS，然后假定如下：•假定其为线性分布，T = TL - ( TL – TS ) fS ，则　 fS/T = -1 / ( TL – TS )                                                                     1-3-20•假定其为二次分布，T = TL - ( TL – TS ) fS２ ，则 TL – T fS = ( ————）1/2　　　 TL – TS 　 fS　　1 TL – T              -1　　　                 -1　　         　 — = —(———)-1/2(———)= —————————— 1-3-21　 T 2 TL – TS                TL – TS           2  (TL – TS )1/2 (TL – TS)1/2  1.3.4.3 1.3.4.3 潜热的实际处理方法潜热的实际处理方法潜热的实际处理方法潜热的实际处理方法•知道了固相率和温度的关系后，可以采用以下的方法求解考虑了潜热放出后的热能守恒式。

901.3 1.3 数值模拟的求解条件数值模拟的求解条件①等价比热法•将式1-3-16中的（Cp-L gS/T)用CpE置换，即 　　　 gS T T  T ρ（Cp-L — ）— = ρCpE — = — (λ— ) 1-3-22 T   X X•根据前面所述的固相率随温度变化的不同情况，就很容易作数值模拟计算•对于凝固温度区间较小的合金，按此法处理，当通过液相线温度和固相线温度时会产生显著的误差 ②热焓法•凝固相变时物质的热焓Ｈ为：     Ｈ＝Ｈ0+∫TT0 Cp dT +（1-fS）L                                                          1-3-23      式中，Ｈ0为基准温度T0时的热焓•如果式1-3-23对温度求导，则得：　 Ｈ 　 fS — ＝ Cp - L — 1-3-24 T 　　　 T•假定质量固相率fS等于体积固相率gS，由式1-3-16得出：911.3 1.3 数值模拟的求解条件数值模拟的求解条件　　 gS T fS T H T H 　 　 T ρ(Cp-L— )— =ρ(Cp-L—)— =ρ— — =ρ— = —(λ—) 1-3-25 T  T  T   X X•由式1-3-25求得Ｈ随时间τ的变化，再根据Ｈ和Ｔ的关系，即可求出Ｔ与τ的变化关系。

•也可以进行下列数学推导，直接求出温度变化例如对二维问题，若热物性值为常数，用热焓法处理潜热，考虑了潜热的热能守恒式为：　　 H 2T 2T ρ — = λ( — + — ) 1-3-26  X2 Y2 •用差商代替微商，下式成立：　 H HP+1 – HP (HP+1 – H0)-(HP – H0) ΔHP+1 -ΔHP — = ———— = ———————— = —————— 1-3-27  Δτ 　　　　　Δτ Δτ•假定质量固相率fS与温度Ｔ近似成线性分布，则： 　　fS = （TL- T）/（TL-TS）       　[ T = TL - ( TL – TS ) fS ]•由式1-3-23得 TL-T     ΔＨ＝Ｈ-Ｈ0=∫TT0 Cp dT +（1-fS）L= Cp（T-T0）+（1- ———）L                                                                                                        TL-TS•所以921.31.3数值模拟的求解条件数值模拟的求解条件 TL-TP+1   ΔＨP+1-ΔＨP =[ Cp（TP+1 -T0）+（1- ———）L ]                                                                    TL-TS TL-TP                           - [ Cp（TP - T0 ）+（1- ——— ）L ]                                                                    TL-TS L                           =  Cp（TP+1 - TP ）（1+ ———— ）                               1-3-28                                                                    (TL-TS)Cp•将式1-3-27式1-3-28代入式1-3-26得　 　 L Tp+1i,j–Tpi,j ρCp [ 1 + ———— ] —————         (TL-TS)Cp Δτ    Tpi+1,j - 2Tpi,j + Tpi-1,j Tpi,j+1 - 2Tpi,j + Tpi,j-1 ＝ —————————— + —————————— ] 1-3-29  (Δx)2 (Δy)2•从式1-3-29可从差分方程直接获得温度场的数值解。

③温度回升法•此法与上述方法相比，其物理意义更明确，适合于数值计算•现在，假定某个领域（体积）中固相率增加ΔgS，其放出的潜热QS可以用下式表示：     QS＝ρV ΔgS L                                                                                      1-3-30931.3 1.3 数值模拟的求解条件数值模拟的求解条件•温度回升法中，先不考虑潜热的放出进行温度计算，求出微小时间Δτ内的温度降低量，ΔT=TL-T如果ΔT>0，就产生凝固，由于放出潜热，温度应该回升到TL（假定没有过冷），因此下式成立：　 QS＝ρCpV ΔT                                                                                      1-3-31•式1-3-31表示：原应降温ΔT，但由于潜热放出使温度回升了ΔT，即温度不下降，其需要的热量正好等于潜热放出量所以： 　ρV ΔgS L = ρCpV ΔT      ΔgS = CpΔT / L                                                                                    1-3-32•此法采用固相率的增加来代替潜热的放出。

如果固相率∑ΔgS＝１，则表示领域Ｖ的凝固结束即使是在等温凝固场合，此法也易于求出固相率•编程步骤：　　 ①计算ΔT=TL-T，如果ΔT>0，则　 ② 计算ΔgS = CpΔT / L                    n                                                                        n-1     ③如果 ∑ΔgS≤１ 则温度回升到TL，否则ΔT’=(1- ∑ΔgS)L/Cp，温度回                            i=1                                                                                                            i=1 升到T+ ΔT’　 941.4 1.4 产品质量预测和优化工艺产品质量预测和优化工艺•计算机数值模拟的最终目的是解决工艺优化设计问题计算机数值模拟的最终目的是解决工艺优化设计问题分析计算机模拟的结果，如果结果不理想，调整工艺参数，再进行计算机模拟，直到模拟结果为最佳结果。

1.4.1 1.4.1 等温曲线法等温曲线法等温曲线法等温曲线法•将某一时刻工件上相同的温度联成线，即得到等温曲线实际数值模拟中常用将温度值和颜色连续起来，将工件上不同的温度用不同的颜色显示，相同的颜色即是相同的温度•颜色常用RGB(0-255,0-255,0-255)函数表示，如果一种颜色代表一个温度值范围，则RGB()函数可以表示256×256×256＝16,777,216个温度范围•用等温曲线法可以预测铸件的缩孔位置，如图1-4-1所示图1-4-1　等温曲线法预测缩孔的位置951.4 1.4 产品质量预测和优化工艺产品质量预测和优化工艺•可以将铸件的临界固相率作为宏观停止流动和补缩的界限，描绘出该界限下铸件内时间变化曲线，即等温（或等固相率）曲线，如果等温（或等固相率）曲线形成了封闭回路，则认为在那个封闭回路内将产生缩孔如图1-4-1所示，图中60s的等固相率曲线形成封闭区间，即区间外为全部是固相，区间内为有液相，在后面的凝固过程中，由于液相凝固收缩得不到补缩，所以会产生缩孔如果在冒口以外的铸件部分不存在封闭的等温曲线，则在铸件凝固过程中补缩通道一直畅通，冒口就能补给液态金属，铸件内就不会产生缩孔。

•钢锭在加热时，奥氏体的形成、晶粒的尺寸等与温度和时间密切相关，用等温曲线法，可以预测完成奥氏体转变的分界线，相同初始奥氏体晶粒的尺寸的区域•等温曲线法在预测缩孔的发生部位上具有一定的可靠性，此法简单而实用但是这种方法只是求出了补充液体金属的流路被“截断”，流动压力损失变的非常大，确实是产生缩孔缺陷的条件•如果液体金属补充流路没有被“截断”，例如当离开冒口距离较远、补缩困难的地方，虽然也属于顺序凝固，但顺序凝固程度较差，相等于等温曲线呈深“Ｕ”字形，仍会在铸件中心线上形成缩孔，这时用等温曲线来作为判断标准，就无法预计在这种情况下缺陷的产生961.41.4产品质量预测和优化工艺产品质量预测和优化工艺1.4.2 1.4.2 温度梯度法温度梯度法温度梯度法温度梯度法•铸件中心线上的缩孔受凝固时温度梯度的支配，如图1-4-2所示a)                  　　　　                                  (b)图1-4-2凝固时温度梯度不同造成补缩难易不同(a)温度梯度大　　　　　　　　(b)温度梯度小•当温度梯度大时，铸件在整个凝固过程中朝向冒口张开的扩张角越大，补缩通道越宽，冒口能充分发挥补缩作用，铸件内就不会产生缩孔与缩松。

温度梯度小到一定程度时，补缩通道将在凝固结束之前的某个时刻“截断”，铸件内将出现缩孔与（或）缩松•温度梯度法就是根据凝固末期的温度梯度的大小来预计缩孔的产生的971.41.4产品质量预测和优化工艺产品质量预测和优化工艺①温度梯度的定义•当计算截面上的某一单元在某一时刻τ的温度为T0时，以二维数值模拟为例，见图1-4-3(a)，经过Δτ时刻后当其温度T0降至评价温度（固相线温度或临界固相率对应的温度）或评价温度以下T’0时，在τ+Δτ时刻计算出从中心点到周围８个相邻点之间的温度梯度G，见图1-4-3(b)，（三维为26个相邻点3×3×3－1），以其中的最大正值作为该单元凝固末期的温度梯度，即：G = max [ ( T’-T’0 ) / ΔL ]                                                   1-4-1     式中T’为与该单元节点相邻的某一单元节点在τ+Δτ时刻的温度；　　　ΔL为与该单元节点相邻的某一单元节点间的距离                           （a）                                                        (b)图1-4-3 温度梯度G的定义(a)计算的单元　　  (b)计算单元的温度梯度981.41.4产品质量预测和优化工艺产品质量预测和优化工艺•若这个温度梯度小于某个临界温度梯度值，则该单元将产生缩孔或缩松。

将同一G值的点连接起来，便得到铸件断面上的等G线，在低于某个临界温度梯度值的等G线内，将产生缩孔或缩松②预测缩孔的临界温度梯度•要使铸件在凝固过程中建立良好的补缩条件而获得致密的铸件，在工艺上应使铸件上远离冒口或浇口的部分到冒口或浇口之间建立一个梯增的温度梯度而温度梯度是随时间和铸件的部位发生变化，因此对一定形状，尺寸的每种合金铸件应有一个不产生缩孔的临界温度梯度值•临界温度梯度值的获得：实际浇成铸件，用超声波或Ｘ射线探伤检验，找到缩孔或缩松的部位，出现缩孔或缩松部位在小于某个温度梯度的等温度梯度线所包围的区域内，此温度梯度则为临界温度梯度判据•铬钼钒铸钢、Cr13Ni5铸钢，其临界温度梯度值GSC2-3℃/cm，壁厚大于100mm取下限值，小于100mm取上限值•当材料进入弹性变形阶段，温度梯度越大，其产生的热应力越大，因此可以用其预测产品的热应力，和产生裂纹的倾向99谢　谢谢　谢 ！！100083 北京科技大学冶金与生态工程学院北京科技大学冶金与生态工程学院 刘刘 青青：：010-62332358(O)：：13331116466 E-Mail: qliu@尼加拉瓜瀑布2000.05.18摄100。