GPU加速有限元接觸傳熱過程

摘 要：在鋼板的沖印過程中，模具的淬火效果在板料的轉變中起著重要的作用。因此，本文對鋼板的淬火過程建立了一個有限元模型，該模型的關鍵點主要是接觸熱傳導過程，利用有限元計算機程序來計算淬火過程溫度場的變化，同時通過溫度測量裝置來得到模具和板料的溫度-時間曲線來驗證計算結果。最后，采用基于GPU（圖形處理單元）技術的并行計算來加快計算。

關 鍵 詞：接觸傳熱；有限元；GPU

引言

接觸傳熱工程技術普遍存在，如板材熱沖壓成形技術、電子器件的冷卻與加熱、高效絕熱結構安裝等。該問題的研究涉及到材料表面特性、溫度、受載環境等多種因素的影響，問題的分析呈現出明顯的耦合特性。使問題求解難度較一般的傳熱計算難度提高，體現出明顯的雙重非線性特征——傳熱非線性與接觸非線性。尋求高效穩定的問題求解方法，同時具有理論和工程應用價值。一般說來，由于問題的復雜性，除了對有限的簡單問題可以得到解析解外，對于多數問題，一般采用有限元方法進行求解。

有限元進行數值計算存在的一大問題是計算時間過長，尋找高效的計算平臺及相適應的算法一直以來都是研究的熱點。近年來利用GPU進行科學運算成為高性能計算的發展方向之一，在圖像處理、數據分析和有限元模擬領域已經取得了諸多成果。利用GPU進行科學運算的優點在于運算速度快、成本低廉，成本和功耗相當于計算機的1/10左右，缺點在于顯存小、計算精度低，作為一種新型的科學計算平臺，傳統的串行算法和適用于編程模式的并行算法均難以充分發揮的GPU運算能力，開發適用于GPU的并行算法是推廣科學運算的迫切需求[1]?；谏鲜龇治?，本研究開發了二維熱傳導有限元程序，并使用GPU改進實現計算的加速。

接觸傳熱計算模型

1.本構方程

接觸傳熱有限元程序設計

根據上述理論分析，按照如圖1所示流程圖設計有限元程序，輸入記錄有限元模型網格節點坐標、單元信息的文件以及相關材料性能參數的文件，遍歷模型中所有單元，根據單元坐標，材料性能參數計算單元剛度矩陣以及單元熱容矩陣，然后根據單元拓撲信息組裝得到總體剛度矩陣K以及總體熱容矩陣C；在每個時間步長內，根據式（3）及結點間的接觸關系計算每個接觸節點的載荷向量，并組裝成總體載荷向量P[3]。接觸傳熱過程屬于瞬態傳熱，在時間域內用有限元單元網格劃分；在時間域內則用有限差分法，將連續的時間變量離散為若干時間步。常用的差分方式有向前差分、向后差分、Crank-Nicolson格式及Galerkin格式，由于Crank-Nicolson格式具有較高的計算精度，并且無條件穩定，故本文采用此種格式，如式（5）所示：

數值算例

以鋼板淬火過程作為算例，如圖2所示，上下為低溫模具，初始溫度為20℃，中間為高溫板料，初始溫度為700℃。計算模型的導熱系數和比熱容如表1和表2所示：接觸導熱系數為4694 W /m2。

計算結果分析

采用的計算平臺為CPU Intel雙核主頻2.6GHz、內存4GB、顯卡GTX260、顯卡核心575MHz、顯存896MB。

為縮短接觸傳熱有限元計算時間，先分析當前程序中各步驟所占時間比例，結果如圖3所示。記錄求解線性方程式（5）在整個有限元程序中所占時間比例，雖然剛度矩陣因為材料的非線性特性，每個時間步長都需要重新計算并組裝，但求解線性方程組仍占據計算的大部分時間，并且占據百分比隨節點數的增加而增加[4]。在目前的情況下，總計算時間占求解方程式（5）的比例可能高達95%。因此，找到一個高效的算法來求解方程式（5）成為加快計算速度的關鍵因素。

根據有限元的理論可知，方程式（5）中的總體剛度矩陣以及總體熱容矩陣均為稀疏矩陣。對于大規模稀疏線性方程組，直接解法的計算量和存儲空間要求都很大。近年來雖然人們發展了如PARDISO和GSS等直接解法求解器，但是對于超大規模的稀疏矩陣求解，迭代解法仍是目前主要的方法，共軛梯度迭代法由于其所需存儲空間和計算量較小且容易進行并行計算，在大型對稱正定稀疏矩陣求解中越來越受到青睞。在共軛梯度法中迭代過程的每步之間是串行的，但向量內積和向量更新、稀疏矩陣與向量乘法都是數據級的，并行操作可以交付給GPU。因此，負責迭代前預置步驟的并行計算以及每次迭代過程中的向量內積和向量更新以及稀疏矩陣與向量乘法，并行計算負責迭代循環和收斂條件的控制，CPU每次迭代過程中的標量除法操作顯存和內存的開辟和清理工作。為考察GPU加速效果，采用標準的線性系統軟件包LINPACK作為標準。最終得到的結果如圖4所示，從圖4中可以看出GPU并行算法加速比隨結點數增加而增加最高可達4.21。其原因在于隨著結點數目的增加、矩陣階數增加、矩陣數目增加、矩陣乘法增加，計算時利用了更多的計算核心算法，并行度更好，更多的節省了計算時間。

參考文獻

[1] 王順緒，戴華.廣義特征值問題的并行塊Jacobi-Davidson方法及應用.計算力學學報[J].2008，25（4）：428-433.

[2] 李洪林，張海蕾，王希誠.一種量子化學有限元并行計算方法.大連理工大學學報[J].2005，25（4）：469-472.

[3] 喻永聲，林家浩.超大型結構特征值問題求解的多重子結構子空間迭代.工程力學[J].2003，20（6）：149-154.

[4] 王勖成.有限單元法[M].北京：清華大學出版社，2003.