焊接热传导的有限差分计算

1、焊接热传导的有限差分计算在焊接工程中经常遇到的一些问题，如焊接热传导、焊接热应力以及焊接过程中的氢扩散等问题，可以归结为解某一个特定的微分方程。然而只有在十分简单的情况下并且作许多简化的假定，才有可能求得这些微分方程的解析解。事实上，由于实际问题多种多样，边界条件十分复杂，用解析法来求解微分方程是十分困难的。为了满足生产和工程上的需要，必须应用近似计算方法。有限差分法和有限单元法都是正在被广泛应用的两种近似计算方法.有限差分法从微分方程出发，将区域经过离散处理后，近似地用差分、差商来代替微分、微商，微分方程和边界条件的求解就可归结为求解一个线性代数方程组，得到的是数值解。有限差分法的长处是：对

2、于具有规则的几何特性和均匀的材料特性问题，差分法的程序设计比较简单，收敛性好，计算过程简单。有限差分法的缺点是局限于规则的差分网格(正方形网格、矩形网格或正三角形网格)，显得死板僵硬。本章介绍有限差分法在焊接传热计算方面的应用。4.1均匀网格中函数的导数偏微分方程组差分解法的基本原理是根据在几个相邻点的函数值来表示函数f (x)在某点的导数，这个点可以在邻点的一侧，也可以位于邻点之间。这些邻点之间的距离可以相等，都等于ax ,也可以各不相等。在本节中节点的距离是均匀的，并等于ax ,下节讨论不均匀间距的问题。fl -工1 x。 11图4-1等网格间距时相邻点的函数值4.1.1直观法考察图4-1

3、所示的函数f(x), x,x0和x1是等间距的节点，x1 - x0 = x0 - x=ax ,亦即xi = x0 +mx , f (x )和fi = f (x | j = f (xi )表示函数值，中间值由赋予x的相应的非整数来表示,例如f2=fx工 2 f,|x0_ . f在x =x-2处的直观估计值为dffo - f 二dx|x - 12 x(4-1)_ . f在x = x12处时的直观估计值为f|f1 - fodx：12.x(4-2)df .如果我们推断这两个式子都是在x=x0处的近似值，则称它们为单侧差分，(4-1)式为后向dx差分，(4-2)式为前向差分。显然，只要利用x0两侧的两点

4、就可以得到在x = x0处较好的估计值竺 _ 3dx|x0- 2 x(4-3)此式称为基于f1和f两值的中心差分，实际上它等于单侧差分的平均值。中心一阶导数公式为dff1 2 - f2dx|x 3x中心二阶导数为df dff - fofo - 3u=f-_dxj-dx_ ax axdx2 x0 dxdxjxf&xxxf1一2 fo f 4x 2其精度随ax变小而提高。(4-4)4.1.2泰勒级数法尽管上述直观方法得到的在 x = x0处的一阶和二阶导数的公式是合理的，是以其相邻点的函数值为基础的，在xo的左边是x1，右边是xi,但不能估计误差，用这种方法也不能建立更精确的公式。而运用泰勒级数却

5、能做到。首先考察如图4-1所示的等距节点的情况。并将函数f(% )= f(% + ax)ax = x1 xo按泰勒级数展开234f f. “.以 f.以 f .以fl_ f。-xfof。f。2!3!4!fo(4)式中fo（4）代表在（xo , x1 ）内某点的四阶导数。同样，将函数f （x）= f（x。 ax）ax = xo - x按泰勒级数展开f x 2：x 3二f。丁。-“ 4f (4) 丁 3，。式中ff0代表在（x，）内某点的四阶导数。当截去（ax）2后面诸项时，从这些展开式可得出f xo ）的前向差分和后向差分公式f _ fi - f。 dx |x=xx,1图4-2不等距节点的函数值

6、df _ f。- 3 , x | _ f,。dx|x=x。x 2式中具有ax阶的误差项。另一方面，将上述两个泰勒级数展开式相减可得df fi - f_i 3x 2 t,、|= a 一1f 4,。+ fo,1 ）dx|x %2 x 3!此式具有（ax f阶的截断误差项，虽然此中心差分公式确实具有二阶截断误差，但并不意味着应用它一定得出比具有一阶截断误差的单侧差分公式更为精确的结果。为了得到x=x。处的二阶导数，将上述两个泰勒级数展开式相加d2f dx2. 2式中误差是（ax）级的。4.2非均匀网格中函数的导数4.2.1 直观法在如图4-2所示不等距节点x,x0和x1的情况下，我们仍然可以做f （

7、x比x = x。处的前向差1 一 ,，“dfdx分与后向差分估计。在 x0 +seax处（是在x0与x1之间的中点）的前向差分估计可用 f1和f0表示 2(4-5)，1 一 ，在x0 324*处（是x和x0的中点）的后向差分用f和f0表不为(4-6)f 1fo - 3dx |x 1x0 4swx swlx中心差分为(4-7)1 一一.1 一由于se asw ,点x =x二十一（se十x = xo十一（se sw x是在图中小右侧。这时, 22此中心差分公式是最好的估计值。取两个单侧导数的加权平均更好一些，因为加权更有利于得到更接近的估计值。如一阶导数fsedf 1swdf 1dx|xb - s

8、e +sw dx|x/swa se +sw dx|xbse(4-8)sef0 - f 二swf1se swswx sw se se =-x2 se se swflsw - sese swf0sesw sw - sefj而对于二阶导数f史d 2 fdx lx x 2 se . x dx r 3 _2sw . xdxf 屋二(se sw) x2f1 f0f0 f(4-9)se xsw x=se sw t21fif1一 x i se se swse swsw sw se4.2.2泰勒级数法对于如图4-2所示的不等距节点的情况，写出f(x1 )= f(x0+seax)和f(x,)=f(x0 -swax

9、)的泰勒级数展开式，并加以处理后，同样也可得到式(4-8)和式(4-9)。4.3稳定态热传导问题的有限差分方程4.3.1偏微分方程替代法对于多变量函数 t =t(x, y,z),把y, z看成常数时，函数t对于x的偏导数就是t对x的普通导数，反之亦然，因而 4.1节和4.2节中所有导数的概念和公式便可以直接应用了。当然，在应用由4.1节和4.2节所得到的公式对 x求偏导时，必须保持 y = y0,z = z0,反之亦然。1 、均匀网格讨论在求解区域上划分均匀网格，并具有均匀导热系数九的三维导热区域。首先只考察内部节点，如果所有边界位于其一个坐标为一常数，例如 x = a的表面上，而且网格系统直

10、接联着边界， 那末，只要已知边界温度场问题，就可以求解内部温度场问题。图 4-3表示了一个典型网格点 p及 其六个相邻点，它们分别记为n, s, e, w, i, o。令a =取=凶=az为均匀网格间距。稳定态热传导的基本方程为2 22t t ftqvd一 2 一 2 一 2二0x 二y二 z式中，qvd为内热源，为单位体积内热量产生的速率。利用 4.1节的结果，可以得到近似式to - 2tptite-2tptwtn- 2tptsqvdn 0 0 0 = 0彳：22式中，如果qvd随位置或温度而变化，则在p点取值。上式可简化为to ti te tw tn ts一 6tpqvd.：2 = 0(4

11、-10)2 .非均匀网格对于非均匀网格，从 p点到n, s, ew, i和o各点的距离写为 snsssea, soa ,那么，4.2节的公式仍可应用。这样，对于非均匀网格的三维问题，可以认为类似于均匀网格的问题，其稳定态热传导的方程式可以表示为211f tn -tp：sn sn -sssnssts +-2 - tess(sn +ss ) _ se(se +sw )1-1 t_ _ tp-. twse swsw (sw 十 se )2111二0.二1totp 1t.-2go&ssos- s-s-so(4-11)式中，如果qvd随位置或温度而变化，则 qvd点在p点取值。注意，当假设所有s都等于1

12、时，便回到均匀网格公式。当没有内热源时，最后一项可以消掉，一般将（4-11 ）式改写如下(4-12)hntnhstsheteawtwa tat。aptpac =0此式等价于（4-12）和（4-10）。例如均匀网格,=a s = a e = aw - a , a o 1= apacqvd .2 &对于非均匀网格,hn2 sn sssnapk sn ssse swsosj；acas, he ,aw, ai,ao 与 an 不中以。式（4-12）是焊接温度场有限差分方程的一般形式。将待求的工件区域划分为网格（不管网格间距是否均匀），有个n节点，对每个节点都有一个（4-12）所示的差分方程。n个节点，

13、对应n个 差分方程（代数方程），联立求解，得出 n个节点的温度值。4.3.2能量平衡法求某一节点温度的相应的有限差分方程也可用能量平衡法推导出来。每个节点用来表示该系统内的一个称为单元控制体积的区域，考察从相邻区域通过这一区域的表面的热流，并推导出一个用相邻节点温度表示的求这一节点温度的公式。这个方法能导出与前面所得到的同样结果，并完全一致，而且十分方便和容易应用，特别是在不均匀网格、对流边界条件和异形区域等的情况下。1.均匀网格为了清楚起见，我们考察三维直角坐标中具有均匀网格间距x = ax = z = az情况下、围绕 p点的单元控制体积，p点及其六个邻近点如图 4-4所示。所有这些单元面

14、（图中为节点e所画的阴影 部分）的面积为a2,这六个点的任何一点到 p的距离都为ao如p点不是边界点。而且假定 w位于 西部边界，并且tw是已知的，则所有内部节点都可彳如下考虑：计算出节点i与其相邻节点间的传图4-4围绕p点的单元控制体积导系数kj ,然后可写出从j到i的热流qj =kj（tj -1）（4-13）在稳定态条件下，对节点i的能量平衡方程为q qj +qvd vi =2 kj% -1）+qvd m =0（4-14）式中工 是对所有单元面求和。此方程是利用了通过六个面进入单元控制体制m的热流加上内热源qvd（w/m3）而得到的。因为所有格距都是均匀的且等于,因此对节点i的单元控制体积

15、vi = axayaz = r内部节点的传导系数用下式计算aijkj =（4-15）j lj式中，九为导热系数， aj为垂直于j和i之间热流方向的平均表面积，lj为两个节点i和j之间的距离。对于三维直角坐标问题，研究图4-4所示的具有六个相邻点的p点的能量平衡，应用公式(4-13)和(4-15)可得到22.2tn -tp ts -tp te - tp aa.-2.-：2.-23tw -tp t -tp to -tp qvd 3 =0& tn -2tptsjte-2tptw t-2tptoqvd=0&aa可见，这与直接的偏微分方程替代是相等的，进而这方程与前面方程(4-10)的形式是一样的q -

16、 2tn ts 丁亚 t . % , 丁 丁。vd-= h2.非均匀网格讨论一下与图 4-4相似，但格距变化的三维问题.从 p点到n, s, e, i和o点的距离用sn,ssa,se4so 等表示，其单元控制体积为1 11vv = sn . ss - se sw- s - - so2 22.-3二8 snss sesw s so北侧和南侧的面积(an , as)为11.：2an = as= 2 & sw ：a ss。心=4 seswss。vv_ a-sn ss2对于其余各节点也可推导出类似的表达式。在此单元控制体积上的热平衡为vv2snsssntn -tp v2 sn ss1_一（ts tp

17、）+四个对 e、w、i 和 o se点的类似项+qvd vv =0式中的单元控制体积 w已在上面给出。此式等于tn tp +1ts!+两个类似项 +q也=0+ ss)snssss(sn+ss)九此式与偏微分方程的直接替代所得到的公式(4-11)完全相同。从上面的讨论中可以看出，建立热传导有限差分方程的两个方法一偏微分方程替代法与能量平 衡法是一样的。因此建议，这两个方法中，那一个对解决当前问题方便就采用那一个。4.4瞬态热传导问题的有限差分方程瞬态热传导与稳定态热传导的主要差别之点是增加了时间的因素。因而不仅需要几何上把区城离散化，确定内节点、边界节点，而且还需要把热过程经历的时间区域离散化。

18、在构成瞬态热传导 的差分方程时，必须特别留意它的稳定性。因为，用不稳定的差分方程进行求解是没有意义的。在 4.3节中，我们已经发现，用偏微分方程替代法和能量平衡法所导出的差分方程是完全一样的，对 于瞬态热传导问题，也是如此。因此，本节主要用偏微分方程替代法建立均匀网格中的差分方程。 参照上节中的方法，可以进行必要的修正以适用于不规则网格。考察在直角坐标中，在整个求解区域布置均匀网格系统的、并具有均匀导热系数的三维导热区域。取 =取=3=&z,并把注意力集中于如图 4-3所示的内部节点 p,在该图中还表示了它的 6个邻点，先暂不考虑内热源，则基本偏微分方程为222(4-16)；titit1 盯=

19、xcy 二za ct式中，a为导温系数，a=一。4.4.1 显式差分格式在式(4-16)中，方程的左边取时间t时的值，方程的右端用前向差分表示，它只包含在p点的t(t +&月t(t)。应用4.1节的结果，则对 p点的有限差分方程为(4-17)to -2tp t te -2tp tw tn -2tp ts1 tp -tp.2r个 一 a ：t式中，tp =tpg+& ) , tp =tpg ),其余不带()的项都在时间t取值。上式整理后，得to ti te tw tn ts-6 tpfo(4-18)式中，fo为傅立叶数，f0=a2t =上。:cp，式(4-18)中tp =tp(t +zkt )可

20、以根据它本身及其相邻六个点在时刻t时的温度来计算，而它们在时刻t的温度是已知的。这样，根据前一个时刻各节点的温度值，可以用(4-18)式直接得出下一个时刻各节点的温度值；一个一个时间步长地推进下去，就可以得出任意时刻各节点的温度值。关于这些公式的稳定性、精确性和收敛性有一些需要考虑的事项，由于这些考虑，对于给定的空间步长限制了 &的极大值，这些将在 4.6节中讨论。如果有内热源的话，则上式可以把 p点在时刻t取值的 q1&2项补充进去。4.4.2 隐式差分格式关于tp的隐式公式包含邻点在新的时间层t +&的其他温度，因此，从前面一个已知的时间层t的值不能简单地计算出 tp,而且不得不为每个时间

21、步长解一组联立方程以求出新的温度，这组方程的数目等于待求温度的节点的总数。各项都在时间层t十工取值时，对p点的有限差分方程为iiiiiiii/ i(4-19)to -2tp +t.te -2tp +tw 3-2tp +ts _ 1 tp-tp 飞z2z2- a m式中的记号同前，此式可化简为to +ti +te +tw +tn +ts 6 +1 t, tpfo(4-20)不论采用什么时间步长 at ,上式对于舍入误差的传播是稳定的，但是不够精确，而且推进一个步长所需要的工作量比更为简单的显式公式要大(因为每个时间步长都必须解一组联立方程)4.4.3 c-n 格式为了进一步减小截断误差，在显式差

22、分和完全隐式差分的基础上又提出crank-nicolson格式，简称c n格式：1 to -2tp +ti +te -2tp +tw .tn -2tp +ts l2 j 及及 一1 to -2tpt-te-2tptwtn-2tpts1tp-tp2_ 飞一%t由此式可得1 一_ _ _.11! to+tw +te +tn+ts、 3 十 tp(4-21)2 foj1 ，丁，丁 ，丁 ，丁 ，丁二 1 =(to +t +tw +te +tn +ts)4 3 tp2 1jc-n格式的最大优点是截断误差小。而且也是稳定的格式，因此常被采用。4.4.2 加权差分格式如果通过下式来应用加权参数6s( 00

23、s 1)：所取之值=(在t+西的值)+(1-)(在t的值)(4-22)则三维问题可表示为 stotitetwtnts- 6tp 厂!1- ?stot i - tetwtn t- tp _tpf0由此可得v f _1入 to +ti +te +tw +tn +ts )- 68s +tp =fo j- 1 - to t tw te tn ts6 1-tp_f0久=0,式(4-24)与式(4-18)相同，对应的差分格式是显式格式。久=1,式(4-24)与式(4-20)相同，对应的差分格式是完全隐式格式.6s =0.5,式(4-24)与式(4-21)相同，对应的差分格式是c-n格式。在19s从。到1的

24、变化过程中，可以得到不同的差分格式。除 a=0的情况外, 分格式都是隐式的。在0aw1条件下，对于给定的 餐与,随着名的增长，计算的精确度下降, 能得到保证。其中当19s=2/3时，称为加列金格式，也是常用的一种差分格式.-6tp(4-23)(4-24)其他所构成的差而稳定性却越4.5 边界、界面、复合传热工件和非均匀物性的工件4.5.1 边界点到目前为止，我们讨论的都是内部节点。现在考察这样一个问题：它有一个覆盖求解区域的网 格系统，区域的边界与网格线相重合。为了使得这个问题能够求解，必须对整个这组节点的各个节点列出有限差分方程。这组节点一般来说由全部内部节点和全部边界节点所组成。各个内部节

25、点的差分方程可以用前面介绍的方法推导出来，现在的任务是要对各个边界点提供一个方程式。就传热问题而言，各边界点常出现两种情况，如图 4-5的b点为例：(1)边界表面温度已知，此时tb=tg(4-25)式中，tg是在规定的问题中所给定的已知温度。因此公式(4-25)为b点提供了一个简单的差分方程。如果没有这个方程，那么 tb就是未知数，于是问题得不到解决。n_.f b e 1图4-5边界节点b(2)边界表面温度是未知的。在这种情况下，需要补充关于外法线方向温度梯度工所取的形二 nft式。可以有下面两种形式：:n：t口 _ =q(qo是已知的)(4-26).n:t=(tb -tf )(对流换热)(4

26、-27).n ,式中，豆为表面放热系数，儿是导热系数，tb和tf分别是边界温度和环境温度。给出这两种一般情况之后，现在的问题就是讨论怎样把公式(4-25)至(4-27)与边界点的有限差分方程结合起来。事实上，式 (4-25)并不存在什么问题，因为如果边界温度是已知的，那么它们就 可以直接写进有关的有限差分方程中去。而其它的边界条件一表面放热一对于稳定态和瞬态，将采 用能量平衡法分别予以考虑。4.5.2表面放热边界条件（稳定态）图4-6说明边界节点b的单元控制体积（面积）是如何建立起来的，该图还表明环境温度 tf是怎样通过传热系数 k而影响tb的。就现在我们这里所讨论的情况而言， 图中的b点是位

27、于表面上的点, 而其四个邻点就是 n、s、e和a（周围环境）。对图4-6虚线所示的b点的单元控制体积，运用能量 平衡法，得到：tn-tb ：ts-tb .te-tb+入+必+口 altf -tb )=0(4-28)2.：2;整理后得出1-tn ts teb3tb = 2 (4-29)2 bir t _ a -式中，bi =是毕欧数。图4-6边界节点b在边上图4-7边界节点b在角上如果边界节点b在角上，如图4-7所示边界点在西南角上那样，公式 （4-29）就不适用了。运用能量平衡法，对于边界点b二u 二二 tf -tb =0（4-30）2:2:22得出上式的过程中，还利用了这样一个事实：（1、，

28、十a ）的表面积所处的环境温度为2ta。（x a经整理，并用bi =毕欧数代替得到1(tnte)bitftb =-1 bi(4-31)4.5.3表面放热边界条件（瞬态）如图4-8所示，在具有边长为 的正方形网格的二维矩形区域里出现两类边界节点。一类是在边上（例如节点1）,另一类是在角上（例如节点o）。对节点0,应用能量平衡法to -t0t3 - tot1 - to=/u十九22 tf-to整理后得到显式形式to =2f。t1t32bitfto2 -2bi 2fo(4-32)也可以隐式形式表示to= 2f0 -t1 -t3 -2bitfto2 2bi 2fo(4-33)式中，fo ；：t。；2c

29、t abi =对节点1,显式形式:t1 =fo to 2t2 t4 2bitf-4 - 2 bi |(4-34 )隐式形式t1 =fo -to -2t2 -丁4 -2bitf t；!fo+ 4 + 2b1(4-35)注意，隐式公式包含几个新的时间层的未知数，因而增加了求解的复杂性。而显式公式则可直接得出新的温度值。图4-8边界节点1和0图4-9材料的界面示意图上面是以二维情况为例推导边界点的差分方程，至于三维问题，用同样的方法也能建立起边界 节点的差分方程。4.5.3 界面传热现象经常发生在复合结构上。这种复合结构由若干不同的材料组成，各种材料之间由一个薄的界面连接。图 4-9描绘了通过这样的

30、界面进行的稳定态导热情况。在界面处有一个明显的温度间断，假想界面的温度差 &tm可由间断处每种材料的局部线性温度分布图的外插法求出，其稳定热流为(4-36 )(4-37)_ 汀t qab a|a - b一 |b:x 1a. x ib单位界面传导系数为qabai 二tnt当两种材料接触很理想时（例如熔化的或焊接的材料），就不存在温度差atint ,于是a t g ,这时界面温度是连续的，即ta =tb。但是一般来说，由于ka丰kb ,所以温度梯度往往是不连续的。用能量平衡法也可以精确地导出有限差分方程。4.5.4 非均匀物性当温度变化不大时，可认为导热系数是一常数。但是如果预先估计到温度会有很大

31、变化，则必须考虑导热系数和温度之间的依赖关系。例如可写成t t 九（丁）（ 4-38 ）如果入与t是线性关系，则可采用下式，仃+、mt)+九3)2 j 2在实际解题的时候，为了计算温度场，可以用 行迭代，直至获得满意的收敛为止。(4-39)九的平均值。根据新的温度场，连续对 八和t进4.6 热传导差分解法的精确性、稳定性和收敛性4.6.1 误差限差分格式的选择以及在计算中傅立叶数所谓误差，指的是有限差分方程的近似解和偏微分方程精确解之间的差值。尽管我们可能不知 道这种精确解，甚至也可能根本不存在这种精确解的表达式。误差有两类，第一类叫截断误差，它 是由于用有限差分代替导数的计算所造成的。这类误

32、差取决于初始给定的温度分布、边界条件、有fo的选择。第二类误差是数值误差，即舍入误差。这类误差是在任何一种计算中对有效数字的限制所引起的，而且对于与傅立叶数相关的时间空间网格来说,这类误差在某些有限差分格式里是固有的。这里指的是那些随着计算解的进行其摆动值的传播和增 长方面不稳定的有限差分格式。收敛性指的是将时间空间网格逐步细分以使得近似计算解接近与精确解的意思。并不是所有的偏微分方程的有限差分法模拟都具有这种收敛性。4.6.2 稳定性所谓差分格式的稳定性，就是计算误差随着步数的增加是否会积累到超出所允许的范围，或者说，最后计算结果对初始条件和边界条件的数据误差及计算中的舍入误差是否敏感的问题

33、。保证差分格式的稳定性在实际运算中是很重要的。这是因为初始条件与边界条件在很多情况下 是测量数据，这些数据不可避免地存在误差，在数值计算中也有舍入误差。如果这些误差在计算过 程中不断地被放大，导致解的不稳定，那么计算结果就会失去了意义。1 .显式差分格式的稳定性现在考察一般的显式瞬态有限差分方程_1l 1 _(4-40)to +ti +te +tw +tn +ts - 6- tp = -tp fofot . ：t_ ,.,一式中，傅立叶数fo =一 =二。事实上，对于a = zkx = ay = az的均匀规则网格，o :2 dtj(j =o, i,n,s,e,w )的系数全是正的。 这里所说

34、的 tj指的是围绕节点p所有六个节点的温度,假如tp的系数均为负值的话(这里所说的tp是指时间层为t时p点的温度)，那么该式表明计算时间 为t时的tp越低，则p点在新的时间层(t + at)时的温度就越高，但是这样就违背了热力学原理。由此得出结论：为了不违背热力学原理所必须满足的条件是：对于给定的空间网格划分应当这样来选择为，以使得tp的系数不为负值。事实证明这就是显式有限差分方程稳定性的充分条件。将这一简单原则用于各种显式差分方程，我们可以推演出下列对于傅立叶数a bi =的限制：九内部节点：，八，一一1一维直角坐标：f0 2，八，一一1二维直角坐标：f0 4，八，一一1三维直角坐标：f06

35、foa. ：t和毕欧数(4-41a )(4-41b )(4-41c )边界节点(表面放热):，八一1一维直角坐标：f0 -2(1 bi )，八，一一1二维直角坐标：f0 -2(2 bi)，八，一一1三维直角坐标：f0 _12(3 bi)(4-42a)(4-42b )(4-42c )二维直角坐标(4-42d )1(角上下点) ：f0 4(1 bi)根据这一简单的原则，即使加上定常的内热源项也不应影响其稳定性。很明显，公式(4-42)的条件对于具有较大的毕欧数的传热系统是不利的，即对于表面传热系数与内部导热系数之比具有较大比值的系统是不利的。例如在给定口，儿和之后，在以显式方法求解 三维不稳定态问

36、题时，时间步长受下式限制:2 2(4-43),i2 l -2-t =f0 _a 2a 3 bi因而当毕欧数bi较大时，为了稳定性的原因，担就只得取小一些。这样就引起工作量的大量增加而 精确度并无显著的提高。解决这个问题的简单办法是不用隐式方法，而是重新安排节点，即将节点 间距增加，这样对稳定性判据就放宽了，但精确度较差。2.隐式格式的稳定性我们也可以对加权平均隐式公式建立稳定性判据，虽然它不完全象刚才讨论的显式形式那样简 单。特别应当指出的是，如果1fo (4-44 )21 -2吸则六点差分公式，即(4-24)式，在内部节点上是稳定的。公式 (4-44)中仇是隐式加权因子。如果1-9s 1 ,

37、式(4-44)显然成立。所以不论fo取何值，差分格式总是稳定的。 完全隐式格式(3=1)、21 .2 cn格式(/=一)、加列金格式(/=一)都属于这个范围，因而都是无条件稳定的格式。2 34.6.3精确度如果差分格式是稳定的，则舍入误差很少或根本不会影响精确度，截断误差对精确度影响最为显著。在内部节点上六点差分公式的误差为：显式：et| e a,& +b1(ax)2(4-45a)隐式：et| a2at +b2(ax) 2(4-45b)c-n格式：et| m a3(at)2 +b3(ax)2(4-45c)可以看出，c-n格式的截断误差与(at 2同阶，c-n格式的精度较高。无论那种差分格式，随

38、着步数的增加，定步长计算的精度都是不断提高的，比较好的方法是采用变步长计算。在开始阶段选小步长，经短时间后，逐步加大步长，达到既保证精度又节约计算机机时的目的。4.7非直角坐标系的热传导有限差分方程非直角坐标更容易表示某些物体的形状。例如，用柱坐标就能更方便地描述轴对称问题。柱坐 标的一般导热方程，也可以用有限差分方程的形式来表示，采用的方法与直角坐标系相同。柱坐标的一般导热方程式如下:222(4-46):t 1 ft1 ; t ftqvd1 ；t -3-=:r r 才 r ；z a ft图4-10和图4-11绘出了柱坐标的几何图形。z, r和8分别具有相等的增量。图4-10柱坐标示意图图4-

39、11柱坐标系的单元控制体积直接替代偏微分方程式(4-46)之后，便得到节点0的有限差分方程如下-2t0 +丁4 +丁2-t4 +丁1 -2t0 +t3 .丁5 -2t0 +丁6 +qvdl3r 22r0mr。2 (日 f3 f 九t0 -t0at(4-47)上式是完全的显式表示。我们也可以利用能量平衡法建立差分方程。参考图4-10和4-11 ,根据由节点0相邻各点传到0点的导热热流速率和 0点所表示的控制体积内的热源放热速率，以及控 制体积内的温度上升速率，建立0点的能量平衡：从l点到0点的导热热流速率=z(azar )tlzt-从2点到0点的导热热流速度=m9az r r0 + it2 -t

40、02t 一 t从3点到0点的导热热流速率=m3 -一-r 衿从4点到。点的导热热流速率=z(a0az r0j4 丁。从5点到0点的导热热流速率=z(r0a6ar t二t0z从6点到0点的导热热流速率=?jr0er t6 -t0 ) z传到控制体积的热量=：cp ro,-.-：z.-：r t0 一70.：t控制体积内热源放热速率 =qvd(r0a0araz )72-2to+t4j_丁4j _2to+t3+ t5- 2to+丁6i 22roarro2(a0 2(iz)九于是，可得出0点的完全能量平衡如下：to -toat(4-48 )公式(4-47)与(4-48)是完全一样的。除了对称轴上的一点外

41、，该式对其它任一节点都是适用的。因为在轴上 =0,因此-i为不定值(=)。但是，利用罗比塔法则，可以证明：当rt 0彳r cr )0时2(4-49)1 ;t二 2t一- -t r cr 二 r这个公式是有定值的。对于圆柱形轴对称且 qvd = 0的有限差分方程，可以取最近邻点的温度作为r = 0处的温度。而t对于非轴对称(,#0 )的稳定态问题，则把其周围各网格点温度的平均值作为r = 0处的温度。c6如果公式(4-48)左边各项都按上一个时间层进行计算(计算不带撇的t),那么，经过适当整理，允许直接用前一个时间层的 0点及其周围六个相邻点的温度来计算t0,这就是显式表示法。如果左边各项都以新

42、的时间层进行计算(计算带撇的t),则得到完全的隐式表达式:i _iiiii .iii .iit2 -2t0+t4t3-t4j -2t,tt2j -2t07丁6=t0-t。(“。) 22r：rr；2：z2- a t般说来，人们可以自由选择哪个时间层来计算左端项，从而分别得出空间导数项的显式的、完全隐式的以及加权平均的差分方程表达式。在前面已经讨论过,to不可能简单地计算出来，为了同时计算各点在新的时间层的 t ,必须解一组代数方程。4.8有限差分方程的计算机解法稳定态和瞬态导热问题的有限差分模拟，将热传导的偏微分方程转换为一组节点的代数方程,使得我们必须同时求解有许多代数方程的方程组。首先讨论含

43、有于有n个节点的传热问题）：n个未知数、n个方程的方程组（对aiiti a12t2 a21t1 a22t2 anlti - an2t2- alntn = bi, a2ntn - b2-anntn 二 bn(4-51 )对于几种不同的传热情况，已经导出了这类有限差分方程。因为每个未知节点温度只与它周围邻点的温度有关，传热差分方程构成的线性代数方程组（4-51）的系数矩阵是稀疏矩阵，采用迭代法求解更为合适。可将方程（4-51）简写成n.二 ai jtj = bi j 1i =1,2, ,n(4-52 )迭代法的基本思想是，构造一个由丁12,，组成的矢量序列，使其收敛于某个极限矢量丁1*万；工*,而

44、且1*万；，匚就是方程组（4-51）的精确解。根据构造矢量序列的方法不同, 有简单迭代法、高斯一赛德尔迭代法与超松驰迭代法等。4.8.1简单迭代法简单迭代法，也称为同步迭代法。迭代的最终目的是求解方程组（4-51）中的t1,t2,，tn。若系数矩阵a对角线上元素不为零，即aii #0 （i =1,2，n）,则可将（4-51 ）式改写成:工=a1； h - a12t2 -a1ntnt2 =a22 b2211- a23t3 - a2ntn(4-53 )tn =an：0 - anx - an2t2 - ann 二二)其中任一方程均可写成:ti = ainb- a/j ijdi =1,2, ,n(4-

45、54 )任意名定tif） （i=1,2;，n作为解的第零次近似，把它们代入（4-54）式的右端，由此可算得:ti1 =anbi j aijtj0j 4jdi =1,2, ,n作为解的第一次近似。把第一次近似得到的解再代入式（4-54）的右端，得到解的第二次近似。一般地，将已得到的解的第 k次近似t（k x弋入式（4-54）的右端，得:nt尸laj bi -z aijtj“）（i =1,2,，n）（4-55）jm 1母）作为解的第 k+1次近似。这样所得到序列tk ,t2k , ,tnk k =0,1,2； ）为方程组的近似解。只要方程组（4-51 ）存在唯一解，则不论零次近似如何选取，当kt

46、g时，此序列 k1,丁2,，tn）必然收敛，且收敛于万程组的精确解t1 ,t2，tn 。实际计算中，k不可能取必，但是可以认为 k充分大时，序列1,丁2,，tn 已足够精确地接近方程组的解。通常，对于充分大的 k值，如果其相邻两次迭代解ti（k书和ti（k ） （i=1,2,n）之间的偏差小于预先给定的适当小量式ea0）,即满足ti 尸）-ti（k|8（i =1,2，,n）（4-56）就可以结束迭代过程，而取 ti（kht） （i = 1,2，n ）作为方程组（4-51）的近似解。4.8.2高斯-塞德尔迭代法简单迭代法虽然计算程序简单，但它计算每一个ti（k+）（i=1,2，n）都要用到全部t

47、（k加值。因此在计算机上，必须有两套工作单元来存放全部未知量的旧值和新值。为了节省工作单元并加快收敛速度，对上述计算过程作如下修改：方程组（4-53）的第一式仍写成t1k 1 = a111bl - a12t2k - - -amtnk,f n、即工+）=&； b1 -s aij#）i j三）方程组（4-53）中第二式中，t（k牍成第一次算得的t1ck卡），即t2kb2 - a2iti(k*n-、a2jtjkj=3按此规律，可得t（k*）的一般关系式:fi ant）=ab _ ajjc） z 近6（4-57）、j4j二4j这种计算方法称为高斯一塞德尔迭代法，也称异步迭代法。用该迭代法时，必须将未知

48、量按顺序排 列，并按顺序逐个进行迭代。由式（4-57 ）可以看出，用高斯塞德尔迭代法进行迭代求解时，只需用一套工作单元存放近似值t。或ti（k41）,这样节省了工作单元，同时在迭代过程中，有一半用迭代的新值，加速了收敛 速度。因此，该法是一种常用的方法。4.8.3 超松驰迭代法实际计算表明，当未知量个数很多时，高斯塞德尔迭代法的收敛速度仍然较慢，作为改进， 人们又提出了超松驰迭代法。对于方程组（4-51 ）,先假定第零次近似的矢量序列为斤1c ,t2（0卜，tn（0现在分两步来求第次近似 打工1,11 .:。 第一步用高斯一塞德尔迭代法：i二nzq=abi -z aijtjo-z 小#0）（i

49、=1,2，n）（4-58）ijtj注书求出第一个近似解zq,第二步按下述公式（4-59）来改善z（）,并得第一次的近似解 t）：(4-59)ti1 = 1 - tj-zj1i =1,2, n式中，切是一个适当选择的参数。用这个来代替zf i由式（4-58）与（4-59）中消去t。）可以彳出由工（后接求出tf牺公式。一般来说，由 t（k）算出tk + xk =0,1,2,，n）的公式为n”t）=（1-sti（k）+sa/ bi -z 小丁产）- j。）（4-60）、 j4j=t+/这就是超松弛迭代法的计算公式，0称为松弛因子。由式（4-60）可见，当0=1时，式（4-60）化简为高斯一塞德尔迭代

50、法的计算公式。因此，超松弛法是更为一般的迭代公式。对于不同的缶值，超松弛法收敛的快慢也不同。可以找出一个最佳*松弛因子0 ,使超松弛法收敛得比别的 切值更快。值范围为10 2。对于00 1的情况称为欠松弛；对于 o a 1的情况称为超松弛。4.9焊接瞬态温度场的有限差分计算实例 瞬态焊接热传导微分方程的一般形式为st 口 f rrr 口 f 口 /(4-61 )八c cl c cl 1 g j . cl ,。 . t 1 pcp=九下十九 十八icte 士dy 工cy ）zz.z.）式中，（之y, z）是固定坐标系，原点是 o;动坐标系为（x,y,z）,原点是o。焊接电弧沿 自轴（x） 以速度v0移动，如图4-12所示。固定坐标系与移动坐标系的关系如下x =x +v0t, x =占一vt因此，有, t _ 工 卫 _ .t .x .t .t,:x. t. x .t . t . t式中，t是动坐标系中的时间。实际上，动坐标系和固定坐标系只是对空间位置而言的；对时间来说，不存在动坐标与固定坐标的区别,冬=1。因为x = xv0t,