鋼鐵冶金過程計算機仿真ppt課件

1、計算機數(shù)值模擬的最終目的是處理工藝優(yōu)化設計算機數(shù)值模擬的最終目的是處理工藝優(yōu)化設計問題。一旦全面實現(xiàn)了冶金加工過程的計算機數(shù)計問題。一旦全面實現(xiàn)了冶金加工過程的計算機數(shù)值模擬，資料的加工消費將會產(chǎn)生深化的變革。值模擬，資料的加工消費將會產(chǎn)生深化的變革。 CAD/CAM/CAE的開展，波音的開展，波音767飛機的設計飛機的設計和加工。和加工。據(jù)、曲線、圖形輸出。據(jù)、曲線、圖形輸出。優(yōu)化工藝優(yōu)化工藝分析計算機模擬的結果，假設結果不理想，調(diào)分析計算機模擬的結果，假設結果不理想，調(diào)整工藝參數(shù)，再進展計算機模擬，直到模擬結果為整工藝參數(shù)，再進展計算機模擬，直到模擬結果為最正確結果。最正確結果。 資料的導

2、熱系數(shù)，資料的導熱系數(shù)，W/(m.K) 溫度梯度，溫度梯度，K/mn負號表示導熱的熱流量向溫度低的方向傳送，即與溫度梯負號表示導熱的熱流量向溫度低的方向傳送，即與溫度梯度的方向相反。度的方向相反。比熱流量是個向量，即它有大小和方向。比熱流量是個向量，即它有大小和方向。假設比熱流量的分量和X、Y、Z坐標系相聯(lián)絡，那么X、Y、Z方向的熱流量分量應是： qx= - , qy= - , qz= - 1-1-2 X Y Z比熱流量 q=iqx+jqy+kqz 1-1-3導熱系數(shù) 物理意義：沿導熱方向的單位長度上，溫度減低,物質所允許經(jīng)過的熱流量。方向性：大多數(shù)液體和固體屬于各向同性的物質。各向異性資料的

3、導熱系數(shù)具有方向性，如石墨。溫度函數(shù)：值還隨溫度而變化。大多數(shù)金屬的導熱系數(shù)隨溫度的升高而降低。大多數(shù)液體水和甘油除外其導熱系數(shù)隨溫度的升高而降低。 氣體的導熱系數(shù)隨溫度的升高而添加。 對流對流是流體氣體和液體中溫度不同的各部分相互混合的宏觀運動引起熱量傳送的景象。工程上最具有實踐意義的是：相對運動著的流體與所接觸的固體壁面之間的熱量交換過程，普通稱為對流換熱。工程上在研討固體壁面和流體之間的對流換熱時，除了高度稀薄的氣體外，人們不去留意流體的單個質點，而把流體看成是延續(xù)介質。實踐的流體總有粘性，流動時，受粘性和壁面摩擦的影響，在接近壁面附近的流體將降低流速，在壁面上完全被滯止不動，即X=0時

4、，=0，如圖1-1所示。因此，熱量從壁面?zhèn)鹘o貼壁的那部分流體，將依托導熱。T(K) V(m/s) Tw qwc Tf V X(m)圖1-1臨近壁面的流體速度分布和溫度分布流體和固體壁接觸面上的“相際熱流密度為： qwc= - f x=+0 1-1-4 X式中：qwc熱流密度，W/m2 f流體的導熱系數(shù)，W/(m.K) T流體的溫度，K液體的導熱系數(shù)較小，氣體的導熱系數(shù)更小，所以受熱時，在貼壁處的流體溫度勢必沿著軸的反方向急劇升高。隨著離壁面的間隔的添加，流體的流動將帶走更多的熱量，使軸方向的溫度梯度延續(xù)下降，溫度分布趨向平坦化。正是經(jīng)過這種導熱和對流的共同作用，使熱量在流體內(nèi)部得到傳播，越臨近

5、壁面，導熱越起主導作用。圖1-1所示，假設厚度為的貼壁靜止膜，膜內(nèi)溫度線性地從壁面溫度TW降到遠離壁面，尚未被加熱的流體溫度Tf，那么 Tf-TW qwc= - f 1-1-5 無界對流時壁面與流體的換熱，鋼錠與周圍空氣的對流換熱屬于這種情況。流體在管和槽道內(nèi)部的流動，稱為有界對流，這時不存在遠離壁面，尚未被加熱的流體溫度，那么采用截面平均溫度作為流體的特征溫度Tf，那么qWC=aC(TW-Tf), W/m2 1-1-6 這就是所謂的“牛頓冷卻定律。式中：aC為放熱系數(shù)，W/(m2.K)其實“牛頓冷卻定律并不是表達對流換熱景象本質的物理定律。凡能影響流體流動的各種要素，包括流體的種類和形狀、運

6、動的速度、流道的外形和大小，以及固體壁面粗糙度等，都會對aC值產(chǎn)生影響。式1-1-6只不過方式地把放熱過程的一切復雜性和計算上的困難，都轉移到并且集中在放熱系數(shù)這樣一個物理量上罷了。aC代表流體和所接觸的固體外表之間溫度每相差，該流體與外表之間“相際熱流量的大小。運用式1-1-6可以進展對流換熱的計算。但由于對流換熱的復雜性，該式中的放熱系數(shù)aC需從相應的準那么方程式求出。準那么方程式是針對不同對流換熱情況，在綜合實驗結果的根底之上，運用類似實際將表征某景象的物理量整理成一些類似準那么，用因次分析法得到的各種類型的表達式。輻射只需溫度高于絕對零度，任何物體都隨時向周圍空間輻射才干。輻射用斯蒂芬

7、玻耳茲曼定律表達： E=，w/m2 1-1-7 式中：物體的輻射率或黑度,(0-1)；斯蒂芬玻耳茲曼常數(shù)或絕對黑度的輻射常數(shù)； 5.6710-8 W/m2 K4溫度,K。實踐上，輻射往往涉及二個物體間輻射熱交換。假設二個物體的外表溫度不同，中間被空氣所隔開，T1T2時，那么相互輻射的結果，外表溫度為T1的物體發(fā)射出去的輻射熱超越了吸收來自外表溫度為T2的物體輻射熱，引起輻射換熱的熱流量那么為： Q12= 12(1-2)F112 或q12= 12(1-2)12 1-1-8 式中：12 物體1與2綜合黑度； F1物體1的外表積；12 物體的外表向外發(fā)射出去的輻射熱量中，能投射到物體外表上的份額，稱

8、為角系數(shù),(0-1)。壁面在氣體中冷卻，存在輻射換熱和對流換熱。思索到壁面與氣體之間存在著輻射換熱，其熱流密度為Qr=arF(Tw-Tf) , w 或qwr=ar(Tw-Tf) , w/m2 1-1-9式中：qwr單位面積的輻射量，w/m2 ar輻射放熱系數(shù)，w/(m2.k)Tw輻射物體外表溫度，k Tf透明的氣體介質的特征溫度，k思索到壁面與氣體之間還存在著對流換熱，其熱流密度為qwc=ac(Tw-Tf) , w/m2 1-1-10 由壁面?zhèn)髯叩目偀崃髅芏萹w應是qwr和qwc二者之和qw=ac(Tw-Tf)+ar(Tw-Tf) 1-1-11 令a0=ac+ar 那么qw=a0(Tw-Tf)

9、 ,w/m2 1-1-12用輻射放熱系數(shù)ar，可以方式地把輻射換熱折合成對流換熱，用總放熱系 數(shù)a0兼顧輻射換熱的影響，從而有利于簡化復雜傳熱的分析和計算。如遠離外表的外界外表溫度趨于環(huán)境溫度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.k) 1-1-16從此可見，(w-

10、f)降為零時，ar并非零值，而是以4 12 m3。隨著溫差的擴展和平均溫度m的升高，ar值將迅速添加。由于ac隨溫差的變化較小，在高溫范圍和大溫差情況下，ar有能夠成為a0的主要組成部分。ar與氣體的運動情況無關，而ac與氣體速度的降低而減小，在氣體自然對流的情況下，ac 5,w/(m2.k)，即使m只需300K，ar就能夠占a0的一半左右。 鋼錠普通屬于各向同性的物質，其加熱或冷卻過程數(shù)學模擬計算根據(jù)的根本數(shù)學模型，是不穩(wěn)定導熱偏微分方程。下面討論各向同性資料導熱方程式的建立。直角坐標系導熱偏微分方程導熱偏微分方程的建立，是經(jīng)過調(diào)查處于導熱過程中的物質的微元體積xy z)的能量平衡來建立。如

11、圖1-2所示。在時間內(nèi)，經(jīng)過六個面的導熱所獲得的能量，加上微元體內(nèi)產(chǎn)生的內(nèi)熱源熱量，要等于微元體積內(nèi)物質積存熱量的改動，即溫度升高或降低。xzydQx+xdQxdQzdQydQz+zdQy+y圖1-2直角坐標系導熱方程式的微元體圖1-2中，微元控制體尺寸x、y、 z，按照傅立葉導熱定律，在X方向流入微元體左外表的熱流可表示為： TdQx=-(y z)，W 1-1-17 X 式中： 導熱系數(shù)，W/(m.k) y z X方向微元體外表積，m2T X方向的溫度梯度，k/mX在X方向流出微元體右外表的熱流，可以運用泰勒級數(shù)展開，保管級數(shù)的第一項和第二項而得到： TdQx+x= dQx + (dQx)x

12、 = 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 用同樣的方法，可以得出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假設單位時間、單位空間所產(chǎn)生的熱量為Q(x,y,z,

13、)，那么微元體的發(fā)熱量為：Q(x,y,z,)x y z 1-1-23由于導熱傳進微元體的凈熱流式1-22和微元體內(nèi)產(chǎn)生的熱量式1-23一同用于增大微元體的內(nèi)能。微元體的內(nèi)能的添加表如今微元體能量存儲隨時間的變化率，即 T Cp x y z 1-1-24 式中： 密度，kg/m3 Cp比熱，J/(kg.k) 時間，s 對微元體進展能量平衡，使能量存儲的時間變化率與由導熱引起的流入微元體的凈熱流和微元體內(nèi)產(chǎn)生的熱量之和相等，得下式： T T T T Cp = ( )+ ( ) + ( )+Q 1-1-25 X X Y Y Z Z 圓柱坐標系導熱偏微分方程實踐問題經(jīng)常涉及柱面對稱問題，且邊境條件給定

14、在一個外表上，因此外表具有一個坐標堅持不變的性質。在這種情況下，采用圓柱坐標系是適宜的。 Z Y X dZ r dr d z O 圖1-3圓柱坐標系導熱方程式的微元體圖1-3所示圓柱坐標系，直接按內(nèi)外兩個圓弧面和其它四個平面組成的微元體積，在導熱過程中熱量必需按收支平衡的原那么導出，此時，微元體積為：dv=(r d)dzdr 1-1-26 沿內(nèi)圓弧面流入微元體積的熱流：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 = d

15、Qr- (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 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 根據(jù)直角坐標系導熱微分方程推導的思緒，可得到： 1

16、T 1 T T (r) + () + () + Q = Cp 1-1-33 r r r r2 z z 球坐標系導熱偏微分方程對于球坐標系，如圖1-4所示。 Z Y X O r rSin rSin Z Y X r dr d O d rSin d rSin 圖1-4 球坐標系導熱方程式的微元體由內(nèi)、外兩個球面、兩個圓弧面和兩個平面所組成的微元體積為： dv=(r d)rSind)dr 1-1-34沿半徑方向流入微元體積的熱流： TdQr=-(r Sind.rd) 1-1-35 r沿半徑方向流出微元體積的熱流： TdQr+dr= dQr + (dQr)dr = dQr + -(rSind.rd)dr

17、 r r r T 1 T = dQr+(-r2)Sind.ddr=dQr- (r2)dv 1-1-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 沿方向流入微元體積的熱流：TdQ=-(dr.rd) 1-1-39 (rSin)沿方向流出微元體積的熱流： dQ+d= dQ + (dQ).(rSin)d (rSin) T = dQ

18、+ (-dr.rd)(rSin)d (rSin) (rSin) T = dQ - ()dv 1-1-40 (r2Sin2)同理可整理得： T T T(r2)+(Sin)+()+Q = Cp r2r r r2Sin r2Sin 1-1-41 求解不穩(wěn)定導熱偏微分方程的數(shù)值解法，主要有：有限差分解法、有限元法、邊境元法三類。邊境元法正在研討和完善之中，目前常用的是有限差分解法和有限元法。我們專門討論有限差分解法的數(shù)學根底，數(shù)值方程的建立，差分方程的穩(wěn)定性和收斂性等問題。有限差分解法是用差分方程近似地替代微分方程，建立差分方程有直接法和能量平衡法兩種。1.2.1直接代換法直接代換法是從微分方式出發(fā)，

19、用差商直接代換微商的方法建立差分方程。1.2.1.1數(shù)學根底微商和差商的定義假設T(x)是延續(xù)函數(shù)，那么其導數(shù)為：dT (x+x)-T(x) = lim = lim 1-2-1dx x0 x x0 x1-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 表示差商替代微商的偏向。誤差多大？需求做誤差分析，才干大膽地運用。誤差分析假設函數(shù)T(x)在x時的值為T(

20、x)，在x+ x時的值為T(x+ x)，假設函數(shù)T(x)在x處的各階導數(shù)存在，那么按照泰勒級數(shù)展開，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等于函數(shù)T(x)在x處的各階導數(shù)的線性組合，只能是近似地等于差商。兩者之間也必然有偏向。圖1-2-1表示了一階差商與一階微商之間的關系

21、。用不同方法得到的差商去替代微商，它們帶來的誤差是不同的。即向前差商：dT/dxT(x+ x)-T(x)/ x 1-2-4向后差商：dT/dxT(x)-T(x- x)/ x 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 x- x x x+ x X T T(x+ x) T(x) T(x- x) dT/dx T(x+ x)-T(x)/ x T(x+ x)-T(x- x)/2 x T(x)-T(x- x)/ x 圖1-2-1差商與微商按照泰勒級數(shù)展開，T(x)與T(x+ x)的關系如下式所示：

22、 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 即向前差商的偏向是截去了泰勒級數(shù)展開式中的高階項而引起的，常稱“截斷誤差，其截斷誤差為與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!

23、 dx2 即向后差商的截斷誤差為與x同級的小量O(x) 。由式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的函數(shù)，對它們還可以求差商。這種一階差商的差商稱為二階差商，它是二階微商的近似，常用向前和向后差商來二階微商，即d2T T(x+ x)-T(x) T(x)-T(x - x) - / x dx2 xx T(x+ x)-2T(x)+T(x- x)

24、= 1-2-12 (x)2由式1-2-7和式1-2-9相加，經(jīng)簡化后得： T(x+ x)-2T(x)+T(x- x)d2T - = O(x)2 1-2-13 (x)2 dx2即二階差商的截斷誤差為與(x)2同級的小量O(x)2 。一維系統(tǒng)假定有一高寬無限(即高寬方向上無熱流)，厚度為L的平板，T=f(x,)即溫度是x方向位置和時間的函數(shù)，所謂一維系統(tǒng)是指幾何空間為一維。初始時辰=0，T=T0，為了簡化，思索無內(nèi)熱源，、Cp均為常數(shù)。選取網(wǎng)格點間距x和時間步長將研討對象離散化。顯式差分格式 T T 一維不穩(wěn)定導熱方程為 Cp = ( ) X X 該方程在區(qū)域0，0 xL內(nèi)全部點都成立，如圖1-2

25、-2所示。將方程運用于內(nèi)節(jié)點 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圖1-2-2一維顯式差分將式1-2-14和式1-2-15代入式1-2-13得： Tp+1i Tpi Tpi+1-2Tpi+Tpi-1 Cp+O( )= + O(x)2 ) 1-2-16 (x)2假設在上

26、式中去掉O( )和O(x)2),整理得： Tp+1i= Tpi+F0(Tpi+1-2Tpi+Tpi-1) 1-2-17 式中F0= /Cp(x)2=導熱速率/蓄熱速率，稱F0為傅立葉數(shù)。 1 2 i-1 i i+1 N p+1 p 1 0 x F0的數(shù)值小意味著加熱或冷卻此物體所需求的時間長，反之，所需求的時間短， F0是一個無因次數(shù)字。當初始條件和邊境條件知時，用式1-2-17就可模擬區(qū)域內(nèi)各節(jié)點隨時間增長的溫度值Tpi(i=2,3,，N-1；p=1,2,3,隱式差分格式 一維隱式差分如圖1-2-3所示，將一維不穩(wěn)定導熱微分方程運用于內(nèi)節(jié)點i，那么： T p 2T p+1 Cp () = (

27、 ) 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= Tp

28、i+F0(Tp+1i+1-2Tp+1i+Tp+1i-1) 1-2-20 式中F0= /Cp(x)2。比較式1-2-17和1-2-20，可以看出顯式差分格式的突出優(yōu)點是每個節(jié)點方程都可以獨立求解，整個計算過程非常簡便。但是，的選取要遭到限制，有時為了滿足差分格式穩(wěn)定性條件，能夠選得很小，使計算任務量加大。隱式差分格式的最大優(yōu)點是，它對恣意F0值都是穩(wěn)定的。這種穩(wěn)定是絕對的，即不受邊境條件、x、熱物參數(shù)的影響，于是可以選的較大，計算速度加快。但是，對于節(jié)點i，只從式1-2-20不能獨立求解，必需涉及Tp+1i+1、Tp+1i、Tp+1i-1的聯(lián)立線性代數(shù)方程組才干求解，也就是說，它含有三個未知數(shù)。

29、時間步長每前進一步，從坐標0 xL，網(wǎng)格點1iN整個區(qū)域的每個點，上述方程都要列出一次見圖1-2-3)。因此，向一個新的時間步長每挪動一步就必需解一個方程組。當按順序列出這些方程時，除了要的第一點和最后一點的方程只需兩個未知數(shù)外，其他每一個方程都含有三個未知數(shù)，于是方程是三對角的。對于這種情況，已有適用的求解程序，后面將討論。關于差分方程穩(wěn)定性將在后面討論。圖1-2-3一維隱式差分二維系統(tǒng)顯式和隱式差分格式建立的方法和兩種差分格式的特點在前面討論過，對二維系統(tǒng)同樣適用。為簡單起見，在此只討論二維系統(tǒng)顯式差分方程的建立。為使問題簡化，依然假定熱物性值為常數(shù)，思索無內(nèi)熱源。首先把所討論的區(qū)域離散化

30、，如圖1-2-4所示。 1 2 i-1 i i+1 N p+1 p 1 0 0 x 網(wǎng)線：用平行于X、Y軸的 直線網(wǎng)線，進 行空間離散化節(jié)點：網(wǎng)線與網(wǎng)線的交點步長：節(jié)點間的間隔圖1-2-4 二維均質物體的分割x與y可以是不變的常量，即等步長，也可以是變量即在區(qū)域內(nèi)的不同處是不同的，即變步長。假設區(qū)域內(nèi)各點處的溫度梯度相差很大，那么在溫度變化猛烈處，網(wǎng)格布得密些，在溫度變化不猛烈處，網(wǎng)格布得疏些。至于網(wǎng)格多少，步長取多少為宜，要根據(jù)計算精度與計算任務量等要素而定。對于0 xL1和0y1/2，由上式可見，當p時辰i節(jié)點溫度Tpi越大，由于(1-2F0)為0時，全部內(nèi)節(jié)點的溫度Tp+1i的值總是處于

31、Tpi+1，Tpi，Tpi-1 三個值的最大值與最小值之間的某個中間數(shù)值。這一現(xiàn)實顯然是符合熱力學第二定律的。為什么隱式差分方程無條件穩(wěn)定？將隱式差分方程稍加整理，可得：Tpi= Tp+1i+F0(2Tp+1i-Tp+1i+1-Tp+1i-1) 1-2-75假定p+1時辰，在區(qū)域內(nèi)某一節(jié)點i處獲得區(qū)域內(nèi)最大溫度值，即Tp+1i為這一時辰區(qū)域內(nèi)的最大溫度,Tp+1i Tp+1i+1 ；Tp+1iTp+1i-1,式1-2-75等號右邊括號內(nèi)三項的代數(shù)和必然大于零，從而必然有TpiTp+1i。也就是說，在p時辰，區(qū)域內(nèi)的最大溫度必然大于p+1時辰區(qū)域內(nèi)的最大溫度。依此類推，必然將最大溫度值推到初始條

32、件。假定p+1時辰，在區(qū)域內(nèi)某一節(jié)點i處獲得區(qū)域內(nèi)最小溫度值，即Tp+1i為這一時辰區(qū)域內(nèi)的最小溫度,Tp+1i Tp+1i+1 ；Tp+1iTp+1i-1,式1-2-75等號右邊括號內(nèi)三項的代數(shù)和必然小于零，從而必然有TpiTp+1i。按照上面的分析方法可知，整個過程的溫度最小值，必然出如今初始條件中?？傊?-2-75不論F0為何值，它的運算邏輯都是符合熱力學原理的。即隱式差分格式是無條件穩(wěn)定的。由上作出結論：隱式差分格式是無條件穩(wěn)定的，顯式差分格式是有條件穩(wěn)定的。這一結論，原那么上對二維、三維系統(tǒng)也是適用的。顯式差分格式的穩(wěn)定性條件是要求Tpi項的系數(shù)0。由式1-2-17一維、式1-2-

33、25二維、式1-2-28三維和式1-2-52二維對流邊境等差分方程可導出詳細的穩(wěn)定性條件如下：內(nèi)部節(jié)點：一維直角坐標：F01/2 二維直角坐標：F01/4 三維直角坐標：F01/6邊境節(jié)點設為對流邊境 1 1 一維直角坐標：F0 二維直角坐標：F02(1+acx/) 2(2+acx/) 1 三維直角坐標：F02(3+acx/)分析 ac=0得到內(nèi)部節(jié)點判據(jù)； 當其它條件不變時，ac增大，收斂條件要求更小，即更嚴厲。根據(jù)這一原那么，其它邊境條件下邊境節(jié)點差分方程的穩(wěn)定性也可以相應導出。對于圓柱坐標系、球坐標系的顯式差分方程，也要按上述原那么分別導出詳細的穩(wěn)定性條件?？傊?，用顯式差分格式數(shù)值求解導

34、熱問題，離散化處置時，時間步長和空間步長x的選取是遭到穩(wěn)定性條件制約的。普通的習慣是先定空間步長，再計算時間步長。顯式差分格式的優(yōu)點：每個節(jié)點方程能獨立求解，整個計算過程非常簡單。缺陷：為滿足穩(wěn)定性條件，空間步長和時間步長的選取必需遵守以下關系： (x)2Cp 一維系統(tǒng)： (F0 1/2) 2 (x)2Cp 二維系統(tǒng)： (F0 1/4) 4 (x)2Cp 三維系統(tǒng)： (F0 1/6) 6例如， 一維系統(tǒng)， x=1mm，=7840kg/m3，Cp=610J/(kg.)， =30J/(s.m.)，求得：0.079,s。即1h約5 萬次。這就闡明：為了滿足精度要求，當網(wǎng)格劃分較細時，時間步長大大減小

35、，當維數(shù)愈高，在同樣的空間步長條件下，時間步長也愈小。因此，為了得到穩(wěn)定而又比較準確的答案，就要進展大量的計算。為了消除穩(wěn)定性時間步長的限制，可以采用隱式格式求解。但是，在每 一個時間層上，對二維系統(tǒng)就要解一個五對角的N*M個未知量的代數(shù)方程組(N為x方向的網(wǎng)格數(shù)，M為y方向的網(wǎng)格數(shù))，即挪動一個時間步長，就必需解整個方程組，每一個方程都有五個未知數(shù)和一個對角占優(yōu)元素。這樣一來，就沒有解三對角方程組那樣方便。因此，需求建立新的差分格式，希望這一差分格式能滿足如下要求：無條件穩(wěn)定；有合理精度的解；所產(chǎn)生的代數(shù)方程組易于求解。交替方向隱式方法正是滿足這些要求的一種差分格式。簡稱IAD1.2.6.1

36、差分方程的建立交替方向隱式方法需求對給定時間步長推導出二組有限差分方程。這些方程是顯式和隱式項的混合，對于前半個時間步長，x方向各項是隱式方式，y方向各項是顯式方式；而在下半個時間步長時那么相反。假設把二維導熱微分方程式運用于節(jié)點(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)2 第二個

37、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層是一個過渡時間層，經(jīng)過式1-2-77計算過渡時間層的值，然后用式1-2-79計算第p+1時辰的值。由于如今有一種迅速解三對角線矩陣的算法，所以交替方向隱式法是有效的。1.2.7三對角線矩陣的解法前面引見了顯式、隱式和交替方

38、向隱式法等差分格式，從數(shù)學方面來看，它們都是一個線性代數(shù)方程組。線性代數(shù)方程組的求解普通有兩類方法：一類是直接方法，如三對角線矩陣的解法、高斯消去法、矩陣分解法；另一類是迭代方法，如高斯塞德爾迭代、超松弛迭代、強隱式迭代。普通說來，直接方法只需進展有限此運算就可以得到方程的解。它的缺陷是必需有較大的計算機存儲量和破費較多的計算時間，而且在運用時， 還必需留意抑制舍入誤差的影響，可以采用添加計算數(shù)字的位數(shù)，或利用主元消去等手段。迭代法沒有舍入誤差積累的問題，普通需求較少的計算機存儲量，計算時間也比較節(jié)省。但是，如何能確保迭代收斂和收斂得比較快，這在運用中必需留意處理。當離散化方程是三節(jié)點關系式時

39、，高斯消去法就轉換成人們所常用的三對角線矩陣法，這種方法也稱為托馬斯(Thomas)算法或TDMA、追逐法。它之所以稱為三對角線矩陣算法，是由于用矩陣方式寫出這些三節(jié)點關系式時，一切系數(shù)都是沿矩陣的三對角線陳列的。設離散化方程是如下方式的三節(jié)點關系式： - aiTi-1 + biTi - ciTi+1 = di 1-2-80i=1,2,N，其中i=1和i=N分別為邊境節(jié)點。由于T0和TN+1沒有定義，故取a1=0和CN=0。因此，當i=1時，式1-2-80變成T1和T2之間的關系式，即可以用T2表示T1，T1=f1(T2)。一個方程二個未知數(shù)當i=2時，式1-2-80變成T1、T2和T3之間的

40、關系式，可以用T2表示T3，T2=f2(T3)。以上二個方程三個未知數(shù)這種交換關系，可以不斷進展到TN 用TN+1來表示，即TN=fN(TN+1)，由于CN=0(TN+1的系數(shù)) ，于是可以求得TN的值。以上添加一個方程，沒有添加未知數(shù)，所以可求出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

41、 ) + biTi - ciTi+1 = di ，經(jīng)整理得： 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的。根據(jù)以上分析，其求解步驟是：取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得T

42、N=QN，即求出點的溫度；對于i=N-1,N-2,2,1，用式1-2-81分別求出TN-1，T2，T1 。1.2.7.1三對角線矩陣解法程序Turbo BASIC言語設A(I)為ai系數(shù)，B(I)為bi系數(shù)，C(I)為Ci系數(shù)，D(I)為di系數(shù)，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(

43、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 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

44、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 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-

45、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 SUBQuick BASIC言語DECLARE SUB TDMA ()REM simul N = 11: M = 19: s = .35: lamda = 30!: cp = 610!: rou = 7840!: dt = 120! dx = s / (

46、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)

47、 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 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 B

48、AP(1 TO N), DAQ(1 TO 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 SUB程序解釋：BAP(

49、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 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) TM

50、ID(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.2將數(shù)值模擬問題變成三對角線矩陣解法方式一維隱式差分格式式1-2-20 Tp+1i= Tpi+F0(Tp+1i+1-2Tp+1i+Tp+1i-1),整理成：-F0 Tp+1i-1 +

51、(1+2F0)Tp+1i F0 Tp+1i+1 = Tpi 所以 ai=F0 bi=1+2F0 ci=F0 di= Tpi 二維交替隱式格式第一個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，F(xiàn)0=/2Cp(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=

52、F0 di=Tpi,j+F0(Tpi,j+1-2Tpi,j+Tpi,j-1) 第二個1/2 Tp+1i,jTp+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)2設x=y，F(xiàn)0=/2Cp(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(

53、Tp+1/2i+1,j-2Tp+1/2i,j+Tp+1/2i-1,j)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(-);prin

54、tf(n); for(i=1;i=N;i+)Ai=F0;Bi=2*F0+1;Ci=F0;Di=15.0; A1=0;B1=1.0;C1=0;D1=100.0;AN=0;BN=1;CN=0;DN=100.0; for(p=1;p=M;p+) TDMA(); printf(% 5.2f ,(dt*p/60); for(i=1;i=N;i+) printf(%5.2f ,TMIDi); for(i=1;i=N;i+) Di=TMIDi; printf(n);getch();void TDMA (void)float BAPN+1,DAQN+1; int i; BAP1=B1;DAQ1=D1; for

55、 (i=1;i=1;i-) TMIDi=(DAQi+Ci*TMIDi+1)/BAPi;運轉結果 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.0

56、0 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

57、.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

58、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-為了使工件內(nèi)的溫度場被獨一地確定下來，必需把一些復雜的要素思索到數(shù)值方程的求解條件中去。如何確定初始條件、邊境條件、熱物參數(shù)隨各種要素的變化、潛熱等。1.3.1初始條件A-B兩物質組成的一維系統(tǒng)，其離散化處置如圖1-3-1所示。圖1-3-1A-B兩物質系統(tǒng)離散化處置所謂初始條件，即計算時間的初始起點，=0時，模擬系統(tǒng)的條件，如物質

59、內(nèi)部點溫度處處為TA0，B物質內(nèi)部點溫度處處為TB0，TA0和TB0也可以是x、y、z坐標的函數(shù)。問題是：在=0時，A-B兩物質交界面溫度Tb0確定為何值？關于Tb0確實定，有以下兩種方法。第一種是實測法，進展實測統(tǒng)計獲得Tb0的普遍規(guī)律。第二種是解析法，即以為-B兩物質接觸的瞬間分界面上溫度堅持不變， 1 2 i-1 i i+1 n-1 n x 物質 物質 且只需與分界面垂直方向有熱交換。根據(jù)一維不穩(wěn)定導熱原理，假設物質為均勻的，且其熱物性值不隨時間溫度而變化，于是可得A-B兩物質交界面的溫度T0inf為： (T0A- T0inf)bA=(T0inf- T0B)bB ，整理得 1-3-1bA

60、 T0A+ bB T0B T0inf = 1-3-2 bA+ bB 式中bA、bB 為、物質的蓄熱系數(shù)，b=Cp假設分界面不在節(jié)點上，如圖1-3-2所示。如圖1-3-2分界面不在節(jié)點上 1 2 i-1 i i+1 n-1 n x 物質 物質 T0b(B 側) T0inf T0b(側) 先求出T0inf，然后分別求出分界面、物質一側的溫度，即： x T0b (A/B)=T0inf+(T0(A/B)-T0inf)erf( ) 1-3-3 2t式中，x為界面兩側節(jié)點離界面的間隔，普通為x/2；為或物質的導溫系數(shù)， =/(Cp)； t為構成A-B界面的時間；erf()為誤差函數(shù)，有誤差函數(shù)表可查得。例