專利名稱：一種有限元求解軋制過程溫度場的集中熱容矩陣方法
技術(shù)領(lǐng)域：
本發(fā)明屬于軋制技術(shù)領(lǐng)域，特別涉及一種有限元求解軋制過程溫度場的集中熱容 矩陣方法。
背景技術(shù)：
在現(xiàn)代工程技術(shù)領(lǐng)域中，有限元法作為一種有效的數(shù)值分析方法已廣泛應(yīng)用于結(jié)構(gòu)、力 和熱等的分析過程中；近年來，許多研究人員采用有限元法分析熱軋和淬火等過程中瞬態(tài)溫 度場。在用有限元法分析瞬態(tài)溫度場時，通常采用的方法是在空間上進(jìn)行結(jié)構(gòu)離散、在時間 上采用有限差分格式進(jìn)行求解。該方法從初始時刻的溫度場開始，每隔一個時間步長進(jìn)行一 次迭代計算，進(jìn)而求出各個時刻的溫度場；其優(yōu)點(diǎn)是節(jié)省計算機(jī)內(nèi)存，但是往往產(chǎn)生時間和
空間上的振蕩現(xiàn)象，影響了計算的穩(wěn)定性和精度。
對于上述振蕩現(xiàn)象，研究人員從差分格式、網(wǎng)格劃分方式、加權(quán)變系數(shù)、迭代變步長、
在實數(shù)子空間上求解和用Norsette法代替差分格式計算等多方面進(jìn)行了深入的研究，得出了 一些降低或消除振蕩的方法，但效果都不太理想。

發(fā)明內(nèi)容
本發(fā)明的目的就是通過采用集中熱容矩陣的方法克服上述采用有限元法求解瞬態(tài)溫度場 時常常產(chǎn)生的時間和空間上的振蕩現(xiàn)象，保證計算的穩(wěn)定性，提高的計算精度。
實現(xiàn)本發(fā)明巨的的技術(shù)解決方案是
① 采集軋制過程數(shù)據(jù)，包括軋制參數(shù)，材料熱物性參數(shù)，單元劃分信息。 軋制參數(shù)軋制時間，軋件寬度，軋件厚度，初始時間，初始溫
度，軋件周圍介質(zhì)溫度，時間步長。 材料熱物性參數(shù)熱傳導(dǎo)系數(shù)，黑度，比熱，密度。 單元劃分信息寬度單元數(shù)和厚度單元數(shù)
② 根據(jù)單元劃分?jǐn)?shù)據(jù)、軋件寬度和厚度尺寸建立有限元分析模型(見圖1)，然后進(jìn)行單 元節(jié)點(diǎn)編號、確定換熱邊界和計算節(jié)點(diǎn)坐標(biāo)。
單元和節(jié)點(diǎn)編號沿厚度方向和寬度方向逐漸增加，圖1中，i為單元編號，j為節(jié)點(diǎn)編號，
H為厚度，W為寬度。確定邊界條件AB和AD邊界絕熱，BC和CD邊界換熱。以A點(diǎn)坐
標(biāo)為零，計算各節(jié)點(diǎn)坐標(biāo)，在寬度方向和厚度方向上單元均勻劃分。③根據(jù)不同軋制過程，確定邊界換熱系數(shù)力。
在板帶軋制過程中，不同軋制階段，具有不同換熱系數(shù)計算模型。熱軋板帶 在空冷過程中，其表面換熱方式主要為輻射和自然對流，輻射系數(shù)表述為
<formula>formula see original document page 5</formula> (i)
其中/B 為輻射系數(shù)，cr為Stefan-Boltzman常數(shù)，cr = 5,67xl(T8r/(m2 .〖4); s為黑度系數(shù)，￡與溫度的關(guān)系式為e = 0.125(771000)2 - 0.38(771000) + 1.1 。
熱軋板帶在高壓水除磷過程中，主要換熱方式為強(qiáng)迫對流和側(cè)面輻射，輻射計算方法同
上，對流系數(shù)表達(dá)式為
<formula>formula see original document page 5</formula> (2)
其中w (i:/min.附"為水流密度；r(《)板帶表面溫度。
在軋制過程中，板帶與軋輥之間接觸換熱是主要熱損失方式。接觸換熱系數(shù)與氧化鐵皮 厚度、氧化鐵皮導(dǎo)熱率和軋制壓力有關(guān)。接觸換熱系數(shù)表達(dá)式為
<formula>formula see original document page 5</formula> (3)
式中凡(M尸")-軋制壓力
利用有限元基本原理，計算四邊形等參單元的形函數(shù)iV、 B矩陣和雅克 比矩陣J。
以二維熱傳導(dǎo)基本方程為基礎(chǔ)，利用歐拉方程建立等效泛函，確定溫度場求解的系統(tǒng) 方程。
以熱力學(xué)第一定律為依據(jù)建立無內(nèi)熱源強(qiáng)度的二維熱傳導(dǎo)的微分方程為
<formula>formula see original document page 5</formula>
式中[&]-溫度剛度矩陣，fc]=i(k(e)]+[At)D; [A]-變溫矩陣，[&]=i;k3w];
常數(shù)項列式，(p} = it(e)}; {T^-溫度列式；五-單元總數(shù)；上標(biāo)e表示每個單元。
<formula>formula see original document page 5</formula>對每個單元來說，剛度矩陣、變溫矩陣和常數(shù)項可以通過下式求解:
'<formula>formula see original document page 6</formula>r-瞬時溫度(a:) p-材料密度(&/w3)
c-材料比熱 f-時間(s)
h熱傳導(dǎo)系數(shù)
利用歐拉方程將二維熱傳導(dǎo)問題方程(4)變?yōu)榈刃Х汉?
<formula>formula see original document page 6</formula>(5)
根據(jù)熱傳導(dǎo)問題的變分原理，對泛函式(5)求一階偏導(dǎo)數(shù)并置零，得到溫度求解的系統(tǒng)方
程:
(6)
式中
其中:
<formula>formula see original document page 6</formula>
yo-材料密度(紐/w3);
c-材料比熱(々(％.尺))；
h-換熱系數(shù)； 7V-形函數(shù)；
,'，y'-節(jié)點(diǎn)編號；
◎利用二點(diǎn)向后差分格式，將系統(tǒng)方程轉(zhuǎn)化為瞬態(tài)溫度場求解的線性方程組。將系統(tǒng)方 程(6)中的溫度對時間偏導(dǎo)數(shù)表示為二點(diǎn)向后差分格式
<formula>formula see original document page 7</formula> (7) 將時間向后差分格式(4)帶入系統(tǒng)方程式(3)得到溫度場求解的線性方程組
<formula>formula see original document page 7</formula>(8)
⑦集中熱容矩陣
式(8)中，[K」是由各單元質(zhì)量陣和比熱的乘積組裝成的矩陣，稱為協(xié)調(diào)質(zhì)量熱容矩陣(簡 稱熱容矩陣)，該矩陣為n (n為節(jié)點(diǎn)總數(shù))階對稱帶狀矩陣。由于并非只有對角元素為零， 由(1)式中各方程可以看出流入某節(jié)點(diǎn)的熱量不僅與該節(jié)點(diǎn)的溫度變化有關(guān)還與其周圍節(jié)點(diǎn) 的溫度變化有關(guān)，這是不合理的。
對于該問題一個合理的解決辦法是將單元的質(zhì)量集中到各節(jié)點(diǎn)上去。 一種可行的做法是 將熱容矩陣的同行或同列元素相加代替對角線元素，新的熱容矩陣只有對角線元素有值，其 余元素均為零，如(3)式所示。
<formula>formula see original document page 7</formula>
集中熱容矩陣后(1)式中各方程的物理意義是十分明確且合理的，即流入某節(jié)點(diǎn)的熱量 與該節(jié)點(diǎn)溫度的變化相平衡。
⑧ 求解線性方程組，獲得溫度場
⑨ 根據(jù)迭代次數(shù)判斷程序計算是否結(jié)束，如果迭代次數(shù)大于設(shè)定迭代次數(shù)則退出程序， 否則繼續(xù)迭代計算。
計算流程如圖2所示。
本發(fā)明在傳統(tǒng)算法的基礎(chǔ)上通過采用集中熱容矩陣的方法，在不影響計算效率的情況下 有效地克服了振蕩現(xiàn)象，能夠保證了計算的穩(wěn)定性，提高計算精度。


圖1本發(fā)明方法有限元分析模型圖，
圖2本發(fā)明方法單元劃分圖，
圖3本發(fā)明方法計算流程圖，
圖4實施本發(fā)明方法之前溫度場的空間分布圖，
圖5實施本發(fā)明方法之前溫度隨時間變化圖，
圖6實施本發(fā)明方法之后溫度場的空間分布圖，
圖7實施本發(fā)明方法之后溫度隨時間變化圖。
圖中i為單元編號，j為節(jié)點(diǎn)編號，H為厚度，W為寬度，l為換熱邊界，P，Q，R為有 限元分析模型上任意選擇的測溫點(diǎn)，H/2為半厚度，B/2為半寬度。
具體實施例方式
以中厚板出加熱爐后空冷過程中橫斷面二維瞬態(tài)溫度場的計算為例加以說明。
(1) 計算條件
鋼板寬B-2.0m，厚H:0.22m，密度P =7800.0 kg/m3，比熱c-670.0 J/(kg.K)、導(dǎo)熱系數(shù) k=30.0 W/(m.K)，與空氣之間的換熱系數(shù)視為溫度的函數(shù)，鋼板的初始溫度均勻分布(丁0=1150 °C)，環(huán)境溫度丁00=25 °C。
(2) 有限元網(wǎng)格
假設(shè)軋件對稱，取斷面四分之一作為研究對象,以其對稱中心為原點(diǎn)坐標(biāo)，采用四邊形等 參單元，寬度方向單元數(shù)n^30，厚度方向單元數(shù)11=10，單元劃分如圖3所示。 集中熱容矩陣前的計算結(jié)果
取迭代時間步長為0.5 s，采用向后差分格式迭代四步后，鋼板橫斷面溫度的計算結(jié)果 如圖4所示，邊界附近節(jié)點(diǎn)的溫度值超出了合理的范圍(TO)，產(chǎn)生了空間上的振蕩現(xiàn)象。
圖5為振蕩區(qū)域中點(diǎn)P (66， 99)、 Q (967， 99)和R (967， 22)的溫度值隨時間的變化 情況；由圖4可以看出溫度值在迭代開始階段都有一個違背傳熱規(guī)律的上升的過程，產(chǎn)生了 時間上的振蕩現(xiàn)象?？臻g和時間上的振蕩影響了計算的穩(wěn)定性和精度，必須加以消除。
集中熱容矩陣后的計算結(jié)果
圖6和圖7為對上述算例采用集中熱容矩陣法后的計算結(jié)果。可見集中熱容矩陣后，空 間上的溫度振蕩得到很好的抑制；邊界振蕩較劇烈區(qū)中的P、 Q禾卩R點(diǎn)的溫度值在迭代開始階 段的上升現(xiàn)象也基本消失。
經(jīng)計算可知，采用集中熱容矩陣后，用有限元法計算瞬態(tài)溫度場可以獲得較為穩(wěn)定精確 的計算結(jié)果。
權(quán)利要求
1、一種有限元求解軋制過程溫度場的集中熱容矩陣方法，其特征在于該方法包括以下步驟(1)采集軋制過程數(shù)據(jù)，包括軋制參數(shù)，材料熱物性參數(shù)，單元劃分信息軋制參數(shù)初始時間，軋制時間，軋件寬度，軋件厚度，初始溫度，軋件周圍介質(zhì)溫度，時間步長材料熱物性參數(shù)熱傳導(dǎo)系數(shù)，黑度，比熱，密度單元劃分信息寬度單元數(shù)和厚度單元數(shù)(2)建立有限元分析模型，進(jìn)行單元節(jié)點(diǎn)編號、確定換熱邊界和計算節(jié)點(diǎn)坐標(biāo)(3)根據(jù)不同軋制過程，確定邊界換熱系數(shù)h熱軋板帶，在空冷過程中，輻射系數(shù)表述為HR＝σ·ε·(T+Tair)(T2+Tair2)式中HR為輻射系數(shù)，σ＝5.67×10-8W/(m2·K4)ε＝0.125(T/1000)2-0.38(T/1000)+1.1熱軋板帶在高壓水除磷過程中，輻射系數(shù)HR，對流系數(shù)表達(dá)式為HCW＝124.7×w0.663×10-0.00147(T-273.16)式中w為水流密度，T板帶表面溫度在軋制過程中，接觸換熱系數(shù)表達(dá)式為IHTC＝695pm-34400(W/m2K)式中pm-軋制壓力(4)利用有限元基本原理，計算四邊形等參單元的形函數(shù)N、B矩陣和雅克比矩陣J(5)以二維熱傳導(dǎo)基本方程為基礎(chǔ)，利用歐拉方程建立等效泛函，確定溫度場求解的系統(tǒng)方程以熱力學(xué)第一定律為依據(jù)建立無內(nèi)熱源強(qiáng)度的二維熱傳導(dǎo)的微分方程為
全文摘要
一種有限元求解軋制過程溫度場的集中熱容矩陣方法，通過將有限元線性方程組中熱容矩陣的同行或同列元素相加代替對角線元素進(jìn)行計算，在不影響計算效率的情況下有效地克服了傳統(tǒng)有限元方法計算瞬態(tài)溫度場時產(chǎn)生的振蕩現(xiàn)象，保證了計算的穩(wěn)定性；該方法計算的中厚板軋制過程中的鋼板表面溫度值與實測溫度的結(jié)果對比表明，該方法保證了計算瞬態(tài)溫度場時具有較高的計算精度。
文檔編號G06F17/50GK101178748SQ20071015898
公開日2008年5月14日 申請日期2007年12月18日 優(yōu)先權(quán)日2007年12月18日
發(fā)明者剛 劉, 劉相華, 李長生 申請人:東北大學(xué)