100mw汽輪機(jī)滑參數(shù)停機(jī)畢業(yè)論文

1、<p>　　1000MW汽輪機(jī)滑參數(shù)停機(jī)優(yōu)化研究之汽輪機(jī)轉(zhuǎn)子溫度場及應(yīng)力場的有限元計(jì)算</p><p>　　1000MW sliding-parameter shutdown of optimization of steam turbine rotor of steam turbine on the study of finite element calculation of temperature f

2、ield and stress field</p><p><b>　　中文摘要</b></p><p>　　汽輪機(jī)在啟停、變負(fù)荷運(yùn)行時(shí)，由于轉(zhuǎn)子表面蒸汽參數(shù)的不斷變化，其內(nèi)部處于非穩(wěn)態(tài)溫度場，將產(chǎn)生較大的溫度梯度及熱應(yīng)力，汽輪機(jī)轉(zhuǎn)子是汽輪機(jī)的最重要部件，汽輪機(jī)轉(zhuǎn)子的熱應(yīng)力是限制機(jī)組啟停速度及安全運(yùn)行的關(guān)鍵因素。因此為了保證機(jī)組運(yùn)行的安全性與經(jīng)濟(jì)性，對汽輪機(jī)轉(zhuǎn)子進(jìn)行溫

3、度場與應(yīng)力場的計(jì)算分析，實(shí)現(xiàn)單元機(jī)組啟動、停運(yùn)優(yōu)化，是電力系統(tǒng)關(guān)注的熱點(diǎn)問題之一。</p><p>　　依據(jù)軸對稱彈性理論，運(yùn)用工具軟件Ansys對1000MW汽輪機(jī)轉(zhuǎn)子進(jìn)行建模，該模型在空間域應(yīng)用有限元網(wǎng)格劃分，在時(shí)間域用有限差分網(wǎng)格劃分，對轉(zhuǎn)子上應(yīng)力集中的部位進(jìn)行網(wǎng)格加密，依據(jù)轉(zhuǎn)子的溫度場與應(yīng)力場的數(shù)學(xué)模型。本文用有限元法對1000MW機(jī)組滑參數(shù)停機(jī)過程中汽輪機(jī)高壓、中壓轉(zhuǎn)子建立二維模型進(jìn)行溫度場、應(yīng)力場的分

4、析計(jì)算，得到機(jī)組滑參數(shù)停機(jī)運(yùn)行下的轉(zhuǎn)子主要部位的溫度場與應(yīng)力場，以此得到不同時(shí)刻的溫度場與應(yīng)力場分布及熱應(yīng)力集中部位應(yīng)力值隨時(shí)間變化的情況，找出高壓轉(zhuǎn)子溫度場、應(yīng)力場的變化規(guī)律，使該汽輪機(jī)轉(zhuǎn)子表面最大合成實(shí)際應(yīng)力控制在280MPa以下，對計(jì)算結(jié)果進(jìn)行綜合分析，對1000MW汽輪機(jī)滑參數(shù)停機(jī)提出合理化建議。</p><p>　　關(guān)鍵詞：汽輪機(jī)，轉(zhuǎn)子，溫度場，應(yīng)力場，有限元計(jì)算， </p><p&

5、gt;<b>　　Abstract</b></p><p>　　Steam turbine in start-stop， variable load operation， because the rotor surface steam parameters， the internal changes in the steady-state temperature field， will pr

6、oduce the large temperature gradient and thermal stress，steam turbine rotor is the most important components， steam turbine rotor thermal stress is restricted unit start-stop speed and safe operation of the key factors.

7、So in order to ensure the security of the unit operation， and economy of steam turbine rotor of temperature field a</p><p>　　According to axisymmetric elasticity theory， using tools software Ansys modeling 1

8、000MW steam turbine rotor in spatial domain， this model finite element mesh division， in the time domain meshing with finite difference of rotor position of stress concentration grid encryption， upon the rotor temperatur

9、e field and stress field of the mathematical model. In this paper， using the finite element method to stop 1000MW units slippery parameters during high voltage and medium-pressure steam turbine roto</p><p>　

10、　Key words: steam turbine， rotor， temperature field， thermal stress field ，Life monitoring </p><p><b>　　目錄</b></p><p><b>　　中文摘要I</b></p><p>　　AbstractII<

11、;/p><p><b>　　目錄III</b></p><p>　　1.引言- 1 -</p><p>　　1.1 課題的目的和意義- 1 -</p><p>　　1.2 國內(nèi)外研究發(fā)展?fàn)顩r- 2 -</p><p>　　1.2.1 國外發(fā)展現(xiàn)狀- 2 -</p><p&

12、gt;　　1.2.2 國內(nèi)發(fā)展現(xiàn)狀- 3 -</p><p>　　1.3 本文研究內(nèi)容- 5 -</p><p>　　2. 溫度場、應(yīng)力場模型的建立- 7 -</p><p>　　2.1 概述- 7 -</p><p>　　2.2 二維離散溫度場模型- 7 -</p><p>　　2.2.1 溫度場模型的建立

13、- 7 -</p><p>　　2.2.2 各部分放熱系數(shù)的計(jì)算- 10 -</p><p>　　2.3 溫度場的變量說明及運(yùn)行程序- 11 -</p><p>　　2.4 應(yīng)力場數(shù)學(xué)模型的建立- 19 -</p><p>　　2.5 汽輪機(jī)轉(zhuǎn)子的合應(yīng)力與強(qiáng)度指標(biāo)- 21 -</p><p>　　2.6 滑參

14、數(shù)停機(jī)應(yīng)力場分析- 22 -</p><p>　　2.7 應(yīng)力場變量說明及運(yùn)行程序- 22 -</p><p>　　2.8 本章小結(jié)- 32 -</p><p>　　3 轉(zhuǎn)子溫度場與應(yīng)力場的有限元計(jì)算- 33 -</p><p>　　3.1 概述- 33 -</p><p>　　3.2 轉(zhuǎn)子有限元分析網(wǎng)格劃分

15、- 33 -</p><p>　　3.3 幾何模型的建立- 35 -</p><p>　　3.4 有限元計(jì)算分析方案- 35 -</p><p>　　3.5 轉(zhuǎn)子溫度場、應(yīng)力場計(jì)算分析- 35 -</p><p>　　3.6 本章小結(jié)- 40 -</p><p>　　4 轉(zhuǎn)子壽命分析- 41 -</p

16、><p>　　4.1 轉(zhuǎn)子壽命損耗評價(jià)- 41 -</p><p>　　4.1.1 轉(zhuǎn)子壽命損耗- 41 -</p><p>　　4.1.2 影響轉(zhuǎn)子壽命的因素- 41 -</p><p>　　4.1.2.1 材料因素- 41 -</p><p>　　4.1.2.2 應(yīng)力因素- 42 -</p>&

17、lt;p>　　4.1.2.3 結(jié)構(gòu)尺寸因素- 42 -</p><p>　　4.2 滑參數(shù)停機(jī)方案的比較- 42 -</p><p>　　4.3 本章小結(jié)- 43 -</p><p>　　5 結(jié)論與展望- 44 -</p><p>　　5.1 本文主要工作- 44 -</p><p>　　5.2 后續(xù)工

18、作展望- 44 -</p><p>　　致 謝- 45 -</p><p>　　參考文獻(xiàn)- 46 -</p><p><b>　　1.引言</b></p><p>　　1.1 課題的目的和意義</p><p>　　能源是國家經(jīng)濟(jì)發(fā)展的重要支柱，電力在能源生產(chǎn)中起著重要的作用，電力生產(chǎn)是現(xiàn)代工

19、業(yè)的的基礎(chǔ)之一，其產(chǎn)量直接或間接的影響一個(gè)國家發(fā)展的現(xiàn)代化水平。很多地方由于受到資源環(huán)境制約，水利資源缺乏，形成以火電為主的電源結(jié)構(gòu)。文獻(xiàn)指出，在2003年，全國發(fā)電量18462億千瓦，其中火電發(fā)電量14193億千瓦，占全國發(fā)電總量的76.8%，所以工業(yè)和生活用電大部分依靠火力電站和核電站。絕大多數(shù)火力和核電站發(fā)電是以加熱蒸汽推動安裝于轉(zhuǎn)子上的葉片旋轉(zhuǎn)，給轉(zhuǎn)子施以動力帶動連接在轉(zhuǎn)子上的發(fā)電機(jī)發(fā)電生產(chǎn)。隨著現(xiàn)代電力資源日趨貧乏 ，電力生產(chǎn)

20、中將逐步加大對新能源的開發(fā)利用，如太陽能發(fā)電。有資料一顯示，國外在太陽能發(fā)電技術(shù)應(yīng)用上，德國Schlaich教授1978年提出利用太陽能加熱氣流形成的動力驅(qū)動渦輪機(jī)帶動發(fā)電機(jī)發(fā)電。另一方面汽輪機(jī)作為原動機(jī)械具有功率大、熱經(jīng)濟(jì)性高、單位功率和制造成本低等一系列優(yōu)點(diǎn)，所以在冶金、化工等部門廣泛應(yīng)用以作為驅(qū)動各種泵、風(fēng)機(jī)等的動力機(jī)械。轉(zhuǎn)子是汽輪機(jī)中是帶動發(fā)電機(jī)旋轉(zhuǎn)，間接把熱能轉(zhuǎn)化為機(jī)械能的部件，是汽輪機(jī)中的重要部件之一。大型汽輪機(jī)機(jī)組轉(zhuǎn)子在運(yùn)

21、轉(zhuǎn)過程中，工作條件惡劣月復(fù)雜，運(yùn)行中不僅要</p><p>　　因此對大容量機(jī)組轉(zhuǎn)子熱應(yīng)力和壽命進(jìn)行研究。合理控制機(jī)組的啟停不僅關(guān)系到機(jī)組</p><p>　　的經(jīng)濟(jì)性，更涉及到機(jī)組的安全性。對轉(zhuǎn)子進(jìn)行熱應(yīng)力和壽命的有限元計(jì)算正是根據(jù)傳熱學(xué)，彈性力學(xué)和疲勞強(qiáng)度等基本原理，通過溫度場和應(yīng)力場的計(jì)算，采用局部應(yīng)力計(jì)算法，結(jié)合金屬材料具體的疲勞試驗(yàn)曲線，而求得轉(zhuǎn)子的低周疲勞損傷，不僅對轉(zhuǎn)子的壽命

22、進(jìn)行正確的預(yù)測，更對在線壽命監(jiān)測系統(tǒng)可靠性做出正確的評價(jià)。</p><p>　　本學(xué)位論文主要針對汽輪機(jī)高中壓轉(zhuǎn)子壽命問題中采用有限元方法進(jìn)行分析根據(jù)其熱應(yīng)力等參數(shù)的分布特性開展研究。結(jié)論，有其深遠(yuǎn)的理論意義和極大的實(shí)際應(yīng)用價(jià)值：</p><p>　　(1) 通過對汽輪機(jī)高中壓轉(zhuǎn)子工況的數(shù)值模擬，可以進(jìn)一步加深對汽輪機(jī)高中壓轉(zhuǎn)子運(yùn)行中的特點(diǎn)認(rèn)識。</p><p>　

23、　(2) 通過對汽輪機(jī)高中壓轉(zhuǎn)子運(yùn)行中的特性認(rèn)識，可以對進(jìn)行汽輪機(jī)高中壓轉(zhuǎn)子熱應(yīng)力分析有一定的理論指導(dǎo)作用。</p><p>　　(3) 通過對汽輪機(jī)高中壓轉(zhuǎn)子的數(shù)值模擬研究，可以清楚地了解其存在現(xiàn)有產(chǎn)生不足問題的主要原因，找到容易產(chǎn)生裂紋或者疲勞的部位，并可以節(jié)省大量的實(shí)驗(yàn)過程和實(shí)驗(yàn)時(shí)間及花費(fèi)；對提高汽輪機(jī)高中壓轉(zhuǎn)子在實(shí)際應(yīng)用中的綜合性能，減少疲勞損傷以增加壽命有著非常重要的價(jià)值。</p><

24、;p>　　1.2 國內(nèi)外研究發(fā)展?fàn)顩r</p><p>　　1.2.1 國外發(fā)展現(xiàn)狀</p><p>　　由于目前轉(zhuǎn)子的溫度和應(yīng)力尚不能直接進(jìn)行測量，只能通過間接方法，建立相應(yīng)的數(shù)學(xué)模型，測量相關(guān)參數(shù)，求出轉(zhuǎn)子金屬溫度和應(yīng)力的變化及壽命損耗?，F(xiàn)在轉(zhuǎn)子應(yīng)力的數(shù)學(xué)模型大多數(shù)是采用一維溫度場理論解的簡化式，其計(jì)算精度較低，只能反映應(yīng)力的變化趨勢，而不能得到應(yīng)力的精確值。若在此基礎(chǔ)上計(jì)算轉(zhuǎn)子在

25、啟停和變負(fù)荷過程中的壽命損耗，將會產(chǎn)生較大的誤差。</p><p>　　國外機(jī)組壽命管理的應(yīng)用在日本、美國和歐洲較為普遍。美國自60年代Gollin電站汽輪機(jī)失事以后，一些大的公司和研究機(jī)構(gòu)如GE、WESTINGHOUSE、EPRI等對轉(zhuǎn)子的安全性更為重視；進(jìn)行了深入的研究。他們將有限元等先進(jìn)數(shù)值方法由于汽輪機(jī)轉(zhuǎn)子的分析計(jì)算，對轉(zhuǎn)子材料的低周疲勞、高溫蠕變、低溫脆性和裂紋擴(kuò)展規(guī)律等諸多方面的問題進(jìn)行了大量的研究，

26、并在汽輪發(fā)電機(jī)組上安裝了應(yīng)力及壽命損耗指示器以指導(dǎo)機(jī)組運(yùn)行。</p><p>　　日本在汽輪機(jī)壽命管理方面也做了很多工作，除了預(yù)測可能出現(xiàn)裂紋的壽命外，還對轉(zhuǎn)子剩余壽命做出計(jì)算。日本的Kagawa University的Ebara等對汽輪機(jī)動葉片采用的12Cr鋼和Ti-6AI-4V合金的疲勞特性進(jìn)行研究，F(xiàn)ujiyama，Kazunari；Takaki，Keisuke；Nakatani，Yujiro等根據(jù)統(tǒng)計(jì)損傷

27、和隨機(jī)損傷仿真研究，對汽輪機(jī)設(shè)備進(jìn)行壽命評估，采用先進(jìn)技術(shù)設(shè)計(jì)汽輪機(jī)流通部分，以提高機(jī)組的性能和設(shè)備的可靠性。另外，日本在無損探傷的研究方面處于世界先進(jìn)水平，日立、三菱重工、東芝、富士機(jī)電等著名大公司相繼提出脆化一腐蝕法、硬度法、金屬組織法、電極化法等無損探傷方法作為改進(jìn)轉(zhuǎn)子壽命評估的手段。</p><p>　　德國的Wichtmann，Andreas研究了高溫對汽輪機(jī)部件的蠕變損傷；Zaviska，Reiche

28、l研究了汽輪機(jī)冷態(tài)啟動過程中的轉(zhuǎn)子溫度變化，在此基礎(chǔ)上建立了冷態(tài)啟動仿真模型； Scheefer，M；Knodler；Scarlin，B等對電廠抗高溫、高壓材料進(jìn)行了探討，一方面是發(fā)展新的材料，一方面是在已有的材料表面噴涂抗氧化性能強(qiáng)的圖層；以及關(guān)于機(jī)組安全經(jīng)濟(jì)運(yùn)行方面也進(jìn)行了大量的研究.</p><p>　　由美國薩金倫迪公司負(fù)責(zé)人總體設(shè)計(jì)，并由美國燃料公司和瑞士蘇而壽公司提供鍋爐，由瑞士ABB提供汽輪機(jī)發(fā)電機(jī)

29、組的超臨界600MW機(jī)組配有與自動啟動序接口的熱應(yīng)力控制裝置，它安裝在高、中壓轉(zhuǎn)子進(jìn)汽端的轉(zhuǎn)子溫度探針和一套計(jì)算回路組成。在啟動升速和加負(fù)荷過程中，該裝置都參與控制。運(yùn)行部門反映該機(jī)組的最大優(yōu)點(diǎn)是啟動速度快，但也提出不能以應(yīng)力作為控制的唯一依據(jù)，還應(yīng)考慮到汽缸的膨脹是否均勻。</p><p>　　英國GEC汽輪發(fā)電機(jī)公司在對大型汽輪機(jī)設(shè)計(jì)的最新進(jìn)展進(jìn)行推論時(shí)認(rèn)為，在不遠(yuǎn)的將來提高目前的蒸汽溫度和采用兩次再循環(huán)的可

30、能性都很小，因而循環(huán)熱效率不會有很大提高?，F(xiàn)代葉型效率已達(dá)到如此水平，將來的改進(jìn)僅是有限的?？深A(yù)見的主要發(fā)展是繼續(xù)增大機(jī)組容量來提高經(jīng)濟(jì)性。研制更大的轉(zhuǎn)子、更長的末級葉片和更大的排汽環(huán)狀面積會有利于提高單機(jī)容量。然而，如果提高單機(jī)功率，采用大轉(zhuǎn)子長葉片的同時(shí)，可靠性達(dá)不到較小機(jī)組的同樣水平，對大型汽輪機(jī)的考慮就會降低，因?yàn)榘l(fā)電廠大型汽輪機(jī)設(shè)備的啟停費(fèi)用是極大的。目前國外先進(jìn)大機(jī)組十分重視根據(jù)熱應(yīng)力狀態(tài)和疲勞壽命自動控制機(jī)組的動態(tài)運(yùn)行，提

31、高機(jī)組動態(tài)運(yùn)行的經(jīng)濟(jì)性。美國EPRI早在十年前就提出其進(jìn)一步的研究方向，其中前兩項(xiàng)是：</p><p>　　(1) 對機(jī)組動態(tài)特性進(jìn)行在線監(jiān)視；</p><p>　　(2) 機(jī)組起停過程中的經(jīng)濟(jì)性能分析，調(diào)峰機(jī)組的性能監(jiān)測和機(jī)組運(yùn)行壽命分析。</p><p>　　1.2.2 國內(nèi)發(fā)展現(xiàn)狀</p><p>　　目前，我國大部分電力由火力發(fā)電機(jī)組

32、承擔(dān)，其中將近三分之一的火力發(fā)電機(jī)組，運(yùn)行時(shí)間已接近或超過其設(shè)計(jì)年限，屬于老化設(shè)備，從目前我國用電負(fù)荷現(xiàn)狀來看，仍需這機(jī)組繼續(xù)發(fā)揮效益，即使在發(fā)電量相當(dāng)大的發(fā)達(dá)國家如美、蘇、英等國家，也都在努力探討如何對這些設(shè)備進(jìn)行壽命管理，進(jìn)一步實(shí)現(xiàn)根據(jù)設(shè)備狀態(tài)進(jìn)行的維修管理。從運(yùn)行方面看，這些火電機(jī)組原先承擔(dān)的基本負(fù)荷已讓位給大容量火電機(jī)組和核電機(jī)組，它們要參與調(diào)峰運(yùn)行。本來由于多年運(yùn)行，轉(zhuǎn)子因熱疲勞、蠕變等因素產(chǎn)生了損傷積累和材質(zhì)退化現(xiàn)象，那么現(xiàn)

33、在的工作條件更為苛刻，勢必又要加速其時(shí)效老化過程。為確定這些老化設(shè)備能否安全的繼續(xù)運(yùn)行，必須對老化設(shè)備進(jìn)行剩余壽命監(jiān)測和預(yù)測，制訂更新失效老化部件的計(jì)劃，實(shí)施適當(dāng)?shù)木S修管理，這對改善這些火電機(jī)組的老化狀況，延長運(yùn)行壽命，是非常重要的.目前有關(guān)機(jī)組調(diào)峰運(yùn)行過程中的熱應(yīng)力變化和壽命管理方面還有若干問題沒有徹底解決。例如在進(jìn)行機(jī)組非穩(wěn)態(tài)溫度狀態(tài)和熱應(yīng)力計(jì)算中蒸汽參數(shù)和換熱系數(shù)的確定，壽命預(yù)測中我國轉(zhuǎn)子用鋼高溫疲勞曲線的確定，都有待進(jìn)一步的研究

34、和完善。汽輪機(jī)在高溫、高壓和高轉(zhuǎn)速的條件下工作，實(shí)際運(yùn)行中參數(shù)的變動、負(fù)荷的波動與設(shè)計(jì)</p><p>　　我國從上世紀(jì)80年代初開始進(jìn)行轉(zhuǎn)子壽命損耗預(yù)測和壽命分配研究。多年來，我國有關(guān)研究機(jī)構(gòu)、高等院校以及制造部門、電廠針對機(jī)組調(diào)峰的需要，以國產(chǎn)機(jī)組為對象，研究了汽輪機(jī)主要零部件在非穩(wěn)態(tài)下的溫度及熱應(yīng)力分布、變化規(guī)律、金屬材料的疲勞特性以及部件的壽命損耗。對國產(chǎn)大容量機(jī)組參與調(diào)峰運(yùn)行的可行性、調(diào)峰運(yùn)行的安全性和

35、經(jīng)濟(jì)性、調(diào)峰幅度進(jìn)行了深入的探討，對低負(fù)荷、少汽無功和兩班制等不同的調(diào)峰方式在經(jīng)濟(jì)性和安全性方面進(jìn)行了理論分析和實(shí)驗(yàn)研究，很多單位都相繼開展了汽輪機(jī)轉(zhuǎn)子應(yīng)力監(jiān)測和壽命損耗計(jì)算的研究工作。</p><p>　　清華大學(xué)的黃仙、倪維斗等提出汽輪機(jī)轉(zhuǎn)子熱應(yīng)力的復(fù)頻法建模，根據(jù)轉(zhuǎn)子熱應(yīng)力的形成機(jī)理，通過建立“蒸汽溫度變化率—轉(zhuǎn)子體全溫差"與“轉(zhuǎn)子體全溫差—轉(zhuǎn)子體內(nèi)外表面熱應(yīng)力"的復(fù)頻域數(shù)學(xué)模型，然后由

36、控制理論方便的得到“蒸汽溫度變化率—轉(zhuǎn)子體內(nèi)外表面熱應(yīng)力”的數(shù)學(xué)模型。倪維斗、黃以騰等在原有汽輪機(jī)轉(zhuǎn)子熱應(yīng)力解析計(jì)算式的基礎(chǔ)上，給出了遞推求解熱應(yīng)力值的算法。利用遞推式，不僅加快了計(jì)算速度，還使得在監(jiān)控過程中用來保存歷史信息所需的存儲空間維持不變。支小牧、寇可新建立了有中心孔轉(zhuǎn)子熱應(yīng)力解析遞推計(jì)算法。該方法計(jì)入了換熱系數(shù)的變化，對汽溫非線性變化采用分段線性疊加的分析方法，這個(gè)對于邊界介質(zhì)溫度變化和換熱情況復(fù)雜的有中心孔的受熱部件的應(yīng)力計(jì)

37、算和研究具有普遍意義，并應(yīng)用雨流法實(shí)時(shí)監(jiān)測和管理轉(zhuǎn)子的壽命損耗。</p><p>　　華中科技大學(xué)的王坤、黃樹紅等對125MW汽輪機(jī)轉(zhuǎn)子的啟停調(diào)峰實(shí)驗(yàn)進(jìn)行了研究，計(jì)算了啟停調(diào)峰過程中的應(yīng)力譜和低周疲勞壽命損耗，探討了如何對轉(zhuǎn)子壽命損耗進(jìn)行優(yōu)化控制的問題。劉華堂、李樹人等對國產(chǎn)200MW汽輪機(jī)參加調(diào)峰運(yùn)行的壽命進(jìn)行了分析計(jì)算。</p><p>　　浙江大學(xué)的陳堅(jiān)紅、阮偉等采用大型結(jié)構(gòu)有限元分

38、析軟件，對冷態(tài)啟動溫升率對轉(zhuǎn)子壽命的影響進(jìn)行了研究。郭乙木、魯祖統(tǒng)在工業(yè)汽輪機(jī)轉(zhuǎn)子材料低周疲勞實(shí)驗(yàn)及對實(shí)際工況下轉(zhuǎn)子瞬態(tài)溫度場、應(yīng)力場分析的基礎(chǔ)上，對轉(zhuǎn)子的有效壽命進(jìn)行了預(yù)測舊。</p><p>　　華北電力大學(xué)的張光、張寶衡等提出一種計(jì)算大型汽輪機(jī)轉(zhuǎn)子熱應(yīng)力的二維離散模型，這種二維離散模型也可以進(jìn)一步被用來計(jì)算轉(zhuǎn)子軸向熱膨脹和監(jiān)測機(jī)組脹差。黃仙、楊昆等對汽輪機(jī)熱應(yīng)力自適應(yīng)模型進(jìn)行了研究。崔健、安利強(qiáng)等對轉(zhuǎn)子熱應(yīng)

39、力計(jì)算中的有限元方法的精度問題進(jìn)行了探討。苑素玲、葛永慶等應(yīng)用MATLAB軟件中SIMULINK軟件包對汽輪機(jī)轉(zhuǎn)子熱應(yīng)力進(jìn)行了實(shí)時(shí)仿真，得到實(shí)時(shí)監(jiān)控下汽輪機(jī)轉(zhuǎn)子熱應(yīng)力對溫升率的仿真曲線。安利強(qiáng)、王璋奇提出了有限元熱應(yīng)力分析的前后處理可視化的實(shí)現(xiàn)方法。劉靜靜、楊昆通過直接測出的高中壓缸內(nèi)缸壁的溫度來模擬高中壓轉(zhuǎn)子表面溫度，作為第一類邊界條件，導(dǎo)出了汽輪機(jī)轉(zhuǎn)子溫度和熱應(yīng)力的計(jì)算公式。安江英、楊昆、張保衡對汽輪機(jī)轉(zhuǎn)子長期在高溫高壓狀態(tài)下材質(zhì)軟

40、化現(xiàn)象進(jìn)行了研究，并建立了硬度與蠕變和應(yīng)力松弛性能之間的解析關(guān)系式。</p><p>　　西安交大對汽輪機(jī)轉(zhuǎn)子抗疲勞的模糊可靠性理論進(jìn)行了分析，主要對汽輪機(jī)葉片的疲勞壽命預(yù)測方法進(jìn)行了研究。國家電力公司熱工研究院也對汽輪機(jī)在變負(fù)荷運(yùn)行的動葉片和疲勞壽命特性進(jìn)行了相關(guān)的研究工作。2002年云南電力實(shí)驗(yàn)研究院與上海成套設(shè)備研究所合作研究了“汽輪機(jī)高壓轉(zhuǎn)子在蠕變和疲勞作用下壽命的設(shè)計(jì)和評定體系，為汽輪機(jī)高壓轉(zhuǎn)子的壽命設(shè)

41、計(jì)、壽命評定和壽命延長奠定了堅(jiān)實(shí)的基礎(chǔ)[2]。</p><p>　　近年來，我國電廠節(jié)能工作受到高度重視。增加高參數(shù)大容量機(jī)組的裝機(jī)容量，汽輪機(jī)通流部分的改造，利用各種節(jié)能理論對火電廠熱力系統(tǒng)進(jìn)行節(jié)能診斷等等，都取得顯著的成效。與此同時(shí)，在火電大機(jī)組參與電網(wǎng)兩班制調(diào)峰和改善大機(jī)組動態(tài)運(yùn)行的安全性與經(jīng)濟(jì)性方面，也做了大量工作。有關(guān)專家學(xué)者在該領(lǐng)域的理論研究中，以汽輪機(jī)轉(zhuǎn)子疲勞壽命為重點(diǎn)，對國產(chǎn)汽輪機(jī)轉(zhuǎn)子冷、熱態(tài)啟動

42、及停機(jī)工況熱應(yīng)力及疲勞壽命的計(jì)算方法與壽命合理分配原則，國產(chǎn)汽輪機(jī)轉(zhuǎn)子不同啟停工況下溫度場、應(yīng)力場、壽命損耗率分析，國產(chǎn)機(jī)組調(diào)峰性能評價(jià)，不同調(diào)峰運(yùn)行方式安全經(jīng)濟(jì)性分析等等，進(jìn)行了多方面的開拓。在汽輪機(jī)轉(zhuǎn)子熱應(yīng)力和疲勞壽命在線監(jiān)測方面，我國自行研究的監(jiān)測裝置己在國產(chǎn)200MW汽輪機(jī)組上投入運(yùn)行，它使運(yùn)行人員直觀地看到各種參數(shù)的連續(xù)變化趨勢，迅速地得出正確地判斷及時(shí)指導(dǎo)運(yùn)行操作，并兼有運(yùn)行參數(shù)巡檢、記錄打印等功能，受到運(yùn)行人員的歡迎。我國

43、力行業(yè)中，在努力實(shí)現(xiàn)機(jī)組最佳動態(tài)運(yùn)行，充分發(fā)揮其經(jīng)濟(jì)效益上己開始取得共識，并從來的起步狀態(tài)逐步走向更加深入的探索。我國按轉(zhuǎn)讓技術(shù)制造和配套的機(jī)組，具有上世紀(jì)80年代初的世界水平。其中300MW和600MW機(jī)組均可</p><p>　　1.3 本文研究內(nèi)容</p><p>　　壽命管理是指以機(jī)組經(jīng)濟(jì)地實(shí)現(xiàn)其服役全壽命為目標(biāo)，在對設(shè)備狀態(tài)進(jìn)行監(jiān)測和評估的基礎(chǔ)上優(yōu)化設(shè)備運(yùn)行與維修管理的新技術(shù):汽

44、輪機(jī)轉(zhuǎn)子是發(fā)電機(jī)組的關(guān)鍵部位，汽輪機(jī)轉(zhuǎn)子的壽命問題是機(jī)組壽命管理的主要問題[3]。汽輪機(jī)在啟?；蜇?fù)荷變動時(shí)，轉(zhuǎn)子金屬內(nèi)部將產(chǎn)生較大的溫度梯度并由此而產(chǎn)生熱應(yīng)力，這種過度狀況下的熱應(yīng)力是影響機(jī)組對負(fù)荷變動的適應(yīng)能力和機(jī)組壽命損耗的重要因素。盡管大型汽輪機(jī)轉(zhuǎn)子的壽命管理技術(shù)已在世界范圍內(nèi)取得了很大的進(jìn)展，技術(shù)日見成熟，但仍存在許多問題需要從理論實(shí)驗(yàn)方面加以研究。這些問題主要集中于轉(zhuǎn)子壽命消耗的計(jì)算方法和壽命預(yù)測的方法，比較焦點(diǎn)的問題包括:&

45、lt;/p><p>　　(1)有限單元法雖然被廣泛采用于溫度場和應(yīng)力場的計(jì)算，但如何更精確更快速的進(jìn)行計(jì)算和分析以建立更有效的計(jì)算模型仍有待研究，包括其他數(shù)值計(jì)算方法的運(yùn)用研究和實(shí)際軟件應(yīng)用研究；</p><p>　　(2)材料的低周疲勞曲線的實(shí)驗(yàn)研究和合理運(yùn)用；</p><p>　　(3)材料老化(包括軟化、脆化)對轉(zhuǎn)子低周疲勞損耗、蠕變損耗和裂縫演變的影響，從而對轉(zhuǎn)

46、子壽命計(jì)算和斷裂預(yù)測的影響；</p><p>　　(4)低周疲勞損傷和蠕變損傷的交互影響和他們之間耦合計(jì)算方法，以及對壽命預(yù)測的影響；</p><p>　　(5)損傷力學(xué)的運(yùn)用，可靠性研究以及斷裂力學(xué)的發(fā)展。</p><p>　　針對這些潛在的研究發(fā)展方向和問題，我的課題的主要研究內(nèi)容是：</p><p>　　(1)滑參數(shù)停機(jī)轉(zhuǎn)子熱應(yīng)力數(shù)學(xué)模

47、型的建立。</p><p>　　1)有中心孔轉(zhuǎn)子外表面、中心部位應(yīng)力計(jì)算公式的推導(dǎo)(無中心孔轉(zhuǎn)子應(yīng)力理論計(jì)算公式的推導(dǎo))。</p><p>　　2)轉(zhuǎn)子有限元建模。</p><p>　　3)轉(zhuǎn)子溫度場的計(jì)算，轉(zhuǎn)子應(yīng)力場的計(jì)算。</p><p>　　(2)考慮轉(zhuǎn)子復(fù)雜應(yīng)力狀態(tài)下低周疲勞壽命模型的建立并編寫程序。</p><p

48、>　　(3)開發(fā)以轉(zhuǎn)子熱應(yīng)力約束條件的超超臨界汽輪機(jī)滑停優(yōu)化的計(jì)算程序。</p><p>　　(4)于汽輪機(jī)轉(zhuǎn)子在實(shí)際運(yùn)行中總是出現(xiàn)疲勞破壞，最后導(dǎo)致斷裂，這將影響整個(gè)機(jī)組的運(yùn)行。根據(jù)上述計(jì)算和分析結(jié)果，推斷出汽輪機(jī)在實(shí)際運(yùn)轉(zhuǎn)過程中由于溫度變化引起的應(yīng)力集中的主要部件，提出一些有建設(shè)性的改進(jìn)，延長汽輪機(jī)的運(yùn)行壽命。</p><p>　　(5)于汽輪機(jī)轉(zhuǎn)子在實(shí)際運(yùn)行中總是出現(xiàn)疲勞破壞，

49、最后導(dǎo)致斷裂，這將影響整個(gè)機(jī)組的運(yùn)行。根據(jù)上述計(jì)算和分析結(jié)果，推斷出汽輪機(jī)在實(shí)際運(yùn)轉(zhuǎn)過程中由于溫度變化引起的應(yīng)力集中的主要部件，提出一些有建設(shè)性的改進(jìn)，延長汽輪機(jī)的運(yùn)行壽命。</p><p>　　2. 溫度場、應(yīng)力場模型的建立</p><p><b>　　2.1 概述</b></p><p>　　汽輪機(jī)在啟動、停機(jī)或負(fù)荷變動時(shí)，轉(zhuǎn)子金屬內(nèi)部將產(chǎn)

50、生較大的溫度梯度并由此而產(chǎn)生熱應(yīng)力，這種過度狀況下的熱應(yīng)力是影響機(jī)組對負(fù)荷變動的適應(yīng)能力和機(jī)組壽命損耗的重要因素。</p><p>　　由熱彈性理論可知，部件中的熱應(yīng)力與部件金屬梯度成正比。溫度梯度是由于汽輪機(jī)啟動、停機(jī)或負(fù)荷變動過程中，部件被連續(xù)加熱或冷卻，使部件處于熱不穩(wěn)定狀態(tài)而造成的，汽輪機(jī)在穩(wěn)態(tài)工況運(yùn)行時(shí)，不存在或者只有輕微徑向溫度梯度和熱應(yīng)力。汽輪機(jī)轉(zhuǎn)子的前軸封段和前幾級，在啟停級負(fù)荷變動過程中，氣溫的

51、變化最為劇烈，導(dǎo)致該部位溫度梯度和熱應(yīng)力最大，成為最危險(xiǎn)的部位，由于設(shè)計(jì)時(shí)間有限，本文只研究了汽輪機(jī)高壓轉(zhuǎn)子的調(diào)節(jié)級葉輪根部的溫度場及應(yīng)力場的分布情況，由于汽輪機(jī)轉(zhuǎn)子是高速旋轉(zhuǎn)的部件，目前尚無直接測量其金屬溫度和熱應(yīng)力的有效手段，所以需要理論計(jì)算來解決。</p><p>　　然而，對于汽輪機(jī)轉(zhuǎn)子這樣的幾何形狀不規(guī)則、邊界條件變化較大的軸對稱體，要獲得精確的計(jì)算結(jié)果，需要應(yīng)用數(shù)值解法。本文所應(yīng)用的有限元法是數(shù)值分析

52、方法的一種，它將連續(xù)彈性體加以離散化，把受載物體看成是由各個(gè)單元所組成的網(wǎng)絡(luò)。單元的形狀可以是各種各樣，本論文所介紹的是三角形單元，各單元在他們的頂點(diǎn)處相連接，這些連接點(diǎn)稱為結(jié)點(diǎn)。有限元法在處理具有復(fù)雜形狀和條件的物體時(shí)有較好的效果，而且單元形狀和稀疏程度可以是任意變化的，它可以用較少的節(jié)點(diǎn)而使區(qū)域達(dá)到更好的近似，因此是較好的分析轉(zhuǎn)子的溫度場和應(yīng)力場的方法。</p><p>　　對于轉(zhuǎn)子來說，可以將轉(zhuǎn)子考慮為軸對

53、稱二維計(jì)算模型，將復(fù)雜的幾何形狀連續(xù)體離散化，并將介質(zhì)對轉(zhuǎn)子表面的放熱系數(shù)和材料物性作為隨啟動時(shí)間的變數(shù)來處理。有限元法是最常用于熱應(yīng)力計(jì)算的數(shù)值解法，它是在能量法和差分法的基礎(chǔ)上發(fā)展而來的，兼具了兩者的優(yōu)點(diǎn)。有限元法的實(shí)質(zhì)是通過單元內(nèi)的變分計(jì)算并總體合成代替微分方程求解，只要劃分單元足夠小，就能用簡單的線性插值函數(shù)得到足夠精度的結(jié)果，而且網(wǎng)格劃分靈活，對不規(guī)則邊界的適應(yīng)性好。</p><p>　　2.2 二維離

54、散溫度場模型</p><p>　　2.2.1 溫度場模型的建立</p><p>　　計(jì)算轉(zhuǎn)子不穩(wěn)定溫度場時(shí)，將轉(zhuǎn)子視為一個(gè)均勻、各向同性、無內(nèi)熱源的物體，屬于軸對稱非定常溫度函數(shù)，溫度T(z，r，r)在區(qū)域D中滿足下列偏微分方程式：</p><p><b>　　(2.1)</b></p><p>　　式中 —— 轉(zhuǎn)子

55、鋼材的導(dǎo)熱系數(shù)，W/(m·℃)；</p><p>　　—— 材料的密度，kg/m3；</p><p>　　Cp—— 材料的比熱，kJ/(kg·℃)。</p><p>　　確定上面微分方程式(2.1)的解，需要滿足的初始條件。，此外還應(yīng)滿足一定的邊界條件。對于汽輪機(jī)轉(zhuǎn)子來說，外表面的邊界條件由蒸汽對轉(zhuǎn)子表面的換熱速度來確定，屬于傳熱學(xué)中的第三類邊界

56、條件，即邊界與介質(zhì)的熱交換條件為已知: </p><p><b>　　(2.2)</b></p><p>　　式中 Tf —— 與轉(zhuǎn)子表面接觸的氣溫；</p><p>　　—— 蒸汽與轉(zhuǎn)子表面的放熱系數(shù)。</p><p>　　由變分原理，公式(2-1)的第三類邊值問題可以等價(jià)轉(zhuǎn)化為求一泛函數(shù)J的極值問題即：。采用有限單

57、元法對區(qū)域D進(jìn)行離散化，在分成的小單元上以上泛函數(shù)式是成立的。這樣盡管T在D中的分布很復(fù)雜，但在每個(gè)小單元可以近似看作是線性分布。由于D是全部單元的總體，所以：</p><p><b>　　(2.3)</b></p><p>　　式中 Je—— 對單元e進(jìn)行計(jì)算；</p><p>　　E—— 在區(qū)域D內(nèi)劃分單元總數(shù)；</p>&

58、lt;p>　　n—— 待求溫度值的節(jié)點(diǎn)個(gè)數(shù)。</p><p>　　邊界單元的泛函數(shù)形式：</p><p>　　+ (2.4)</p><p>　　內(nèi)部單元泛函數(shù)形式：</p><p>　　若假設(shè)溫度T(z，r，)在單元e中呈線性分布，在單元作變微分計(jì)算，并經(jīng)過一系列</p><p><b>　　

59、的處理，最后得：</b></p><p><b>　　(2.5)</b></p><p>　　式中 i、j、k—— 三角形單元三個(gè)頂點(diǎn)編號，如圖2.1所示。</p><p><b>　　i</b></p><p><b>　　j </b></p>

60、;<p>　　j k</p><p>　　圖 2.1 i、j、k三角形單元三個(gè)頂點(diǎn)編號</p><p>　　將式(2.5)帶入式(2.4)中，可得到n階線性代數(shù)方程組，以求出n個(gè)節(jié)點(diǎn)的溫度，其形式為：</p><p><b>　　(2.6)</b></p><p>　　

61、式中 ——剛度矩陣(nn階正定對稱方陣)；</p><p>　　——非穩(wěn)態(tài)溫度矩陣(nn階正定對稱方陣)；</p><p>　　、、和——分別為溫度、熱載荷和溫度變化率列向量(n1階方陣)。</p><p>　　其中 項(xiàng)用有限元差分法展開。</p><p>　　本文通過對各種差分格式的比較，選用伽遼金格式，則式(2.6)變?yōu)椋?lt;/p

62、><p><b>　　(2.7)</b></p><p>　　式(2.7)穩(wěn)定且不振蕩的條件是：</p><p><b>　　0</b></p><p>　　對于軸對稱有限元分割，可近似地得出第一類邊界條件時(shí)時(shí)間步長 限制條件：</p><p>　　式中x—— 單元的平均邊長。&

63、lt;/p><p>　　伽遼金格式的特點(diǎn)是在超過上式的范圍時(shí)發(fā)生的震蕩比較輕微。詳細(xì)的推導(dǎo)在文獻(xiàn)上都有詳細(xì)的介紹，由于篇幅所限，不再詳述。二維對稱溫度場的計(jì)算精度關(guān)鍵在于正確處理邊界條件。溫度場計(jì)算過程中，根據(jù)邊界位置及當(dāng)時(shí)的蒸汽溫度和壓力，實(shí)時(shí)地計(jì)算該處蒸汽的各項(xiàng)物性參數(shù)及換熱系數(shù)。溫度場計(jì)算分為穩(wěn)態(tài)溫度場計(jì)算和瞬態(tài)溫度場計(jì)算[4]。</p><p>　　2.2.2 各部分放熱系數(shù)的計(jì)算&l

64、t;/p><p>　　在啟停過程中，轉(zhuǎn)子外表面上蒸汽溫度和放熱系數(shù)隨不同的軸向位置和時(shí)間而變化，亦即是時(shí)間和空間的函數(shù)。有關(guān)汽輪機(jī)轉(zhuǎn)子表面放熱系數(shù)的計(jì)算公式，國內(nèi)外尚未統(tǒng)一定論，放熱系數(shù)對熱應(yīng)力的計(jì)算影響不大，當(dāng)啟動未進(jìn)入準(zhǔn)穩(wěn)態(tài)時(shí)，不同計(jì)算公式所得熱應(yīng)力有一定差值，進(jìn)入準(zhǔn)穩(wěn)態(tài)之后，放熱系數(shù)的大小對熱應(yīng)力值無明顯影響。本文采用哈汽一南工公式，轉(zhuǎn)子各部位放熱系數(shù)公式如下：</p><p>　　(1

65、)葉輪兩側(cè)的放熱系數(shù)：</p><p>　　W/(m2·℃) (2.8)</p><p>　　式中 Nu—— 努賽爾特征數(shù)，可按下式計(jì)算：</p><p>　　當(dāng)Re<2.4時(shí)，Nu=0.675</p><p><b>　　當(dāng)Re時(shí)，</b></p><p

66、>　　式中 —— 氣流的運(yùn)動粘度系數(shù)，m2/s；</p><p>　　u —— 葉輪外圓Rb處得圓周速度，m/s；</p><p>　　Rb —— 葉輪的外圓半徑，m；</p><p>　　Re —— 介質(zhì)雷諾數(shù)；</p><p>　　—— 蒸汽導(dǎo)熱系數(shù)，W/(m·℃)。</p><p>　　(2)

67、汽封處的放熱系數(shù)(曲齒)：</p><p>　　W/(m·℃) (2.9)</p><p>　　式中 d—— 汽封間隙，mm；</p><p>　　S—— 汽封齒距，mm；</p><p>　　H—— 轉(zhuǎn)子表面到汽封環(huán)的高度，m；</p><p>　　wi —— 縫隙

68、中氣流的平均速度，m/s。</p><p>　　(3)光軸處的放熱系數(shù)：</p><p>　　W/(m·℃) (2.10)</p><p>　　式中 Ra—— 光軸外圓半徑，m。</p><p>　　(4)輪緣處的放熱系數(shù)：</p><p>　　W/(m·℃)

69、 (2.11)</p><p>　　式中 —— 葉片材料的導(dǎo)熱率 ；</p><p>　　r0—— 葉輪頂部的半徑，m。</p><p>　　也可以近似處理，高、中壓轉(zhuǎn)子葉輪外緣處的蒸汽放熱系數(shù)349 W/(m·℃)、160 W/(m·℃)【5】。</p><p>　　(5)調(diào)節(jié)級兩

70、側(cè)的放熱系數(shù)：</p><p><b>　　當(dāng)Re<2.4時(shí)，</b></p><p><b>　　當(dāng)Re時(shí)，</b></p><p>　　式中 r0 —— 葉輪半徑，m；</p><p>　　W —— 葉輪旋轉(zhuǎn)角速度，rad/s； </p><p>

71、　　—— 葉輪切向速度，m/s；</p><p>　　r —— 葉輪的平均半徑，m； </p><p>　　u —— 蒸汽的粘度，kg/(m·s)；</p><p>　　Cp —— 蒸汽的平均比熱，kJ/(kg·℃)； </p><p>　　1 —— 蒸汽的密度，kg/m3；</p>&l

72、t;p>　　—— 運(yùn)動粘度，kg/(m·s)。</p><p>　　2.3 溫度場的變量說明及運(yùn)行程序</p><p>　　溫度場程序框圖如圖2.1所示，采用Basic語言對程序進(jìn)行編制。</p><p><b>　　標(biāo)識符說明：</b></p><p><b>　　(1)簡單變量</b&

73、gt;</p><p>　　A8表示單元邊界上的對流換熱系數(shù)(W/m2·℃)；</p><p>　　A9表示單元邊界上對流換熱系數(shù)和介質(zhì)溫度Tf(℃)的乘積；</p><p>　　BI，BJ，BM分別表示公式(5.19)中的bi、bj、bm；</p><p>　　B6，B7分別第二、三類邊界單元的分組數(shù)；</p>&l

74、t;p>　　CI，CJ，CM分別表示公式(5.19)中的ci、cj、cm；</p><p>　　C4，C5分別表示圖紙尺寸放大倍數(shù)和時(shí)間變步長倍數(shù)；</p><p>　　圖2.1 溫度場程序圖</p><p>　　D1表示三角形單元的面積；</p><p>　　分別表示內(nèi)部、第一類和第二類邊界單元的最后一個(gè)編號； </p>

75、<p>　　表示單元總數(shù)，也是第三類邊界單元的最后一個(gè)編號；</p><p><b>　　表示公式中的；</b></p><p>　　判別量，和2分別作穩(wěn)態(tài)和非穩(wěn)態(tài)溫度場計(jì)算；</p><p>　　判別量，只作單元三結(jié)點(diǎn)編號信息的分解；</p><p>　　判別量，和1分別表示軸對稱問題和平面問題；</

76、p><p>　　分別表示單元頂點(diǎn)i，j，m的編號；</p><p><b>　　表示結(jié)點(diǎn)總數(shù)；</b></p><p>　　表示第一類邊界單元的邊界結(jié)點(diǎn)數(shù)；</p><p>　　分別矩陣的行號和列號；</p><p>　　分別表示組成單元材料的總數(shù)和該單元的材料類別；</p><p

77、>　　判別量，分別表示向后差分格式，C-N格式和伽遼金格式；</p><p><b>　　表示公式中的乘積；</b></p><p>　　分別表示第二類和第三類邊界單元的數(shù)目；</p><p>　　RI，RJ，RM分別表示單元三頂點(diǎn)的縱坐標(biāo)；</p><p>　　表示系數(shù)矩陣下三角陣中非零帶內(nèi)元素的總數(shù)；</

78、p><p>　　SI表示邊界單元的邊長；</p><p>　　表示物體的均勻初溫(℃)；</p><p>　　表示初始時(shí)間步長t0(s)；</p><p>　　表示非穩(wěn)定導(dǎo)熱過程所進(jìn)行的時(shí)間(s)；</p><p>　　表示變時(shí)間步長，若作定時(shí)間步長計(jì)算，則；</p><p>　　表示規(guī)定時(shí)間推進(jìn)積

79、分的總次數(shù)；</p><p>　　表示時(shí)間推進(jìn)積分的累計(jì)次數(shù)，為初始時(shí)刻；</p><p>　　XI，XJ，XM分別表示單元三頂點(diǎn)的橫坐標(biāo)。</p><p><b>　　(2)組數(shù)</b></p><p>　　用來保存P數(shù)組在非穩(wěn)態(tài)導(dǎo)熱計(jì)算中的初值；</p><p>　　用來保存E數(shù)組在非穩(wěn)態(tài)導(dǎo)熱

80、計(jì)算中的初值；</p><p>　　表示第二類邊界單元各分組中的最后一個(gè)單元號，并記</p><p>　　用來保存K數(shù)組在非穩(wěn)態(tài)導(dǎo)熱計(jì)算中的初值；</p><p>　　表示式(5.63)中N矩陣下三角元素一維壓縮排列的數(shù)組；</p><p>　　表示第三類邊界單元的對流換熱系數(shù)；</p><p>　　表示第三類邊界單元

81、的周圍介質(zhì)溫度；</p><p>　　表示第三類單元信息的數(shù)組；</p><p>　　表示第三類邊界單元各分組的周圍介質(zhì)溫度Tf(℃)；</p><p>　　表示第三類邊界單元各分組的對流換熱系數(shù)(W/m2·℃)；</p><p>　　表示公式(5.63)中K矩陣下三角元素一維壓縮排列的數(shù)組；</p><p>

82、;　　表示第一類邊界單元的邊界結(jié)點(diǎn)溫度(℃)；</p><p>　　表示第二類邊界單元的熱流密度；</p><p>　　表示主對角元素在一維壓縮排列中的地址；</p><p>　　表示第一類邊界單元的邊界的結(jié)點(diǎn)號；</p><p>　　表示線性代數(shù)方程組的右端列向量，最后存放解向量；</p><p>　　表示單元的單位

83、體積熱產(chǎn)生率；</p><p>　　表示結(jié)點(diǎn)的縱坐標(biāo)y或r，輸入單位用；</p><p>　　作為臨時(shí)存放非零帶內(nèi)元素的一個(gè)數(shù)組；</p><p>　　存放非穩(wěn)態(tài)導(dǎo)熱前一時(shí)刻的溫度場；</p><p>　　表示第二類邊界單元各分組的熱流密度，流出取正值，流入取負(fù)值；</p><p>　　表示第三類邊界單元各分組中的最后

84、一個(gè)單元，并記；</p><p>　　表示結(jié)點(diǎn)的x坐標(biāo)(橫坐標(biāo))，輸入單位用；</p><p>　　表示單元的導(dǎo)熱系數(shù)(W/m2·℃)；</p><p>　　表示單元的容積比熱。</p><p><b>　　(3)子程序</b></p><p>　　DIV用來分解單元信息；</p&

85、gt;<p>　　ADDRESS用來確定主對角元素在一維壓縮排列中的地址；</p><p>　　JUDGE矩陣消元或相乘時(shí)為查明零元素而用的判據(jù)；</p><p>　　SOALG用消元法求解線性代數(shù)方程組；</p><p>　　TH11用來計(jì)算系數(shù)矩陣元素及內(nèi)熱源對右端列向量的貢獻(xiàn)；</p><p>　　BOUND用來計(jì)算邊界條

86、件對右端列向量的貢獻(xiàn)；</p><p>　　DIFF表示非穩(wěn)態(tài)導(dǎo)熱計(jì)算中的三種時(shí)間差分格式；</p><p>　　NSTEA作式(5.84)中非穩(wěn)態(tài)導(dǎo)熱右端列向量的加法和懲罰運(yùn)算。</p><p>　　10 REM TEMPERATURE FIELD SOLVED BY FEM</p><p>　　12 READ L0，L1，E0，E1，E2

87、，E3，M0，C4，B6，B7，H0，T0，T1，C5，U0 </p><p>　　20 Q2=E2-E1:Q3=E3-E2 </p><p>　　22 DIM A(L0)，N(L0，)P(L0)，R(L0)，S(L0)，T(L0)，X(L0)</p><p>　　24 DIM F(Q3)，G(Q3).H#(E3)，I

88、(B7)，J(B7)，W(B7) </p><p>　　26 DIM Q(M0)，Y(M0)，Z(MO)</p><p>　　40 FOR A0=1 TO L0:READ X(A0)， R(A0)</p><p>　　42 PRINT X(A0);R(AO):NEXT A0</p><p>　　44 FOR A0=1 TO E3:READ H

89、#(A0) </p><p>　　46 PRINT H#(A0):NEXT A0</p><p>　　48 FOR A0=1 TO M0:READ Y(A0)，Z(A0)，Q(A0)</p><p>　　50 PRINT Y(A0)，Z(A0)，Q(A0)</p>

90、<p>　　52 NEXT A0 </p><p>　　54 IF Ll=0 THEN 62 </p><p>　　56 PRINT''ST B.C.NODE NUMBER AND TEMPERATURE</p><p

91、>　　58 FOR A0=1 TO Ll:READ，O(A0)，L(A0) </p><p>　　60 PRINT O(A0).L(A0):NEXT A0 </p><p>　　62 IF Q2=0 THEN 76 &

92、lt;/p><p>　　64 PRINT''2ND B.C.GROUP LAST ELEMENT NUMBEP.& HEAT FLUX''</p><p>　　66 FOR A0=0 TO B6:READ C(AO)，V(AO) </p><p>　　68 PRINT C(A0);V(A0):NEXT A0 <

93、/p><p>　　70 FOR A0=1 TO B6:V(A0)=V(A0)*(.000001)</p><p>　　72 FOR E9=C(A0-l)+l TO C(A0) ，</p><p>　　74 M(E9-El)=V(A0):NEXT E9:NEXT A0</p><p>　　76 IF Q3=0 THEN 94

94、 </p><p>　　78 PRINT"3RD B.C.GROUP LAST NO.FILM COEF. & MEDIUM TEMP''</p><p>　　80 FOR A0=0 TO B7</p><p>　　82 READ W(A0)，J

95、(A0)，I(A07)</p><p>　　84 PRINT W(A0);J(A0); I(A0 ) :NEXT A0 </p><p>　　86 FOR A0=1 TO B7:J(A0):J(A0)*.00000l </p><p>　　88 FOR E9=W(A0-l)十1 TO W(A0)</p><p>　　90 F(E9

96、-E2)=J(A0):G(E9-E2)=I(A0)</p><p>　　92 NEXT E9:NEXT， A0</p><p>　　94 FOR A0=1 TO L0 </p><p>　　96T(A0)=T0:NEXT A0 </p>&l

97、t;p>　　98 IF HO=0 THEN 104</p><p>　　100 PRTNT"PLANE TEMPERATURE PROBLEM'' </p><p>　　102 GOTO 106</p><p>　　104 PRINT"AXIS SYMME'I'RY TEMPRATURE PROBLEM

98、9;'</p><p>　　106IF ABS(C4-1)<.001 THEN 112 </p><p>　　108 FOR A0=1 TO L0:X(AO)=X(AO)/C4</p><p>　　110 R(A0)=R(A0)/C4:NEXT A0 <

99、/p><p>　　112 GOSUB l030</p><p>　　1l4 S0=N(L0):PRINT''S0='' ;S0</p><p>　　116 DIM K(S0)，B(S0).D(S0)，E(S0)</p><p>　　118 FOR A0=1 TO M0:Y(A0)=Y(A0)*A (.001)<

100、;/p><p>　　120 Q(A0)=Q(A0)*(1E-09) </p><p>　　122 IF T0>99999! THEN 126</p><p>　　124 Z(A0)*(1E-09)</p><p>　　126 NEXT A0</p><p>　　128 T2=0:T7=T1:U9=1

101、 </p><p>　　130 FOR A0=1 TO S0:K(A0)=0:E(A0)=0:NEXT A0 </p><p>　　132 FOR A0=l TO LO:P(A0)=0:NEXT A0</p><p>　　134FOR E9=l TO E3:G1=2

102、:GO SUB 1000 </p><p>　　13b IF Dl>0 THEN l44 </p><p>　　138 PRINT''FAIL IN ELEMENT NUMBER E9='' ;E9</p><p>　　140 GOTO 260<

103、;/p><p>　　144 IF H0=0 THEN 148</p><p>　　146 RI=I:RJ=1:RM=1，</p><p>　　148 F0=Y(M2)*(RI+RJ+RM/(12*D1)</p><p>　　152 N2=D1*(M2)</p><p>　　153 GOSUB 1250</p>

104、<p>　　154 I7=I1:J7=J1:B8=BI:B9=BJ</p><p>　　156 C8=C1:C9=CJ:R7=RI:R8=BJ</p><p>　　158 R9=RM:GOSUB 1200</p><p>　　160 I7=J1:J7=M1:B8=BJ:B9=BM</p><p>　　162 C8=CJ:C9=CM

105、:R7=RJ:R8=RM</p><p>　　164 R9=RI:COSUB 1200</p><p>　　166 I7=M1:J7=I1:B8=BM:B9=BI</p><p>　　168 C8=CM:C9=CI:R7=RM:R8=RI</p><p>　　170 R9=RJ:GOSUB 1200</p><p>　

106、　200 NEXT E9</p><p>　　202 IF T1>9999! THEN208 </p><p>　　204 G0=2:N0=1:GOSUB1300</p><p>　　206 T2=T2+T7:GOSUB1350

107、 </p><p>　　208 IF L1=0 THEN 216</p><p>　　210 FOR A0=1 TO L1:I0=O(A0)</p><p>　　212 K(N(I0))=K(N(I0))*1E+08</p><p>　　214 P(I0)=K(N(I0))*L(

108、A0):NEXT A0</p><p>　　216 GOSUB 1100:IF T1<99999! THEN 224</p><p>　　218 PRINT"STEADY TEMPRATURE FIELD"</p><p>　　220 FOR I0=1 TO L0:PRINT P(I0);</p><p>　　224

109、 PRINT"NONSTEADY TEMPRATURE FIELD"</p><p>　　226 PRINT"TIME=";T2;'STEP SIZE=";T7;"C5=';C5</p><p>　　228 FOR I0=1 TO L0</p><p>　　230 PRINT"T(

110、";A0;")=";P(A0)</p><p>　　232 T(A0)=P(A0):P(A0)</p><p>　　234 NEXT A0:U9=U9+1:PRINT</p><p>　　236 FOR I0=1 TO S0:K(I0)=D(I0)</p><p>　　238 E-(I0)=B(I0):NEXT

111、I0</p><p>　　240 IF C5<1.2 THEN 244</p><p>　　242C5=C5-.05 </p><p>　　244 T7=T7*C5:IF U9<U0+.1 THEN 204</p><

112、;p><b>　　260 END </b></p><p>　　300 DATA 5，0，3，3，3，4，1，1，0，1，1，20，1，1，3</p><p>　　302 DATA 0，0，10，0，5，5，0，10，10，10</p><p>　　304 DATA 0.0010030041，0.0010020031，0.00300500

113、41</p><p>　　306 DATA 0.0030020051，16，1E6，0</p><p>　　308 DATA 3，0，0，4，1600，100</p><p>　　1000 REM SUBROUTIME OF DIV</p><p>　　1002 Z1#=H#(E9)*1000:I1=INT(Z1#)</p>&

114、lt;p>　　1004 Z1#=(Z1#-I1)*1000:J1=INT(Z1#)</p><p>　　1006 Z1#=(Z1#-J1)*1000:M1=INT(Z1#)</p><p>　　1008 IF G1<1.5 THEN 1024</p><p>　　1010 Z1#=(Z1#-M1)*10:M2=INT(Z1#+.5)</p>

115、<p>　　1012 XI=X(I1):XJ=X(J1):XM=X(M1)</p><p>　　1014 RI=R(I1):RJ=R(J1):RM=R(M1)</p><p>　　1016 BI=RJ=RJ-RM:BJ=RM-RI:BM=RI-RJ</p><p>　　1018 CI=XM-XJ:CJ=XI-XM:CM=XJ-XI</p>