一種有限元求解軋制過程溫度場的集中熱容矩陣方法

2023-05-31 22:05:46

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

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


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