劉瑜， 朱可可， 萬秋里， 王成恩

(1.上海交通大學 機械與動力工程學院，上海 200240； 2.中國航空發(fā)動機集團有限公司沈陽發(fā)動機研究所，沈陽 110015)

0 引 言

當在更嚴格和更復雜的條件下進行結構導熱研究時，非線性導熱問題分析變得更加重要，如渦輪葉片導熱計算、功能梯度材料導熱計算等。導熱問題中存在各種非線性，例如，非線性熱源和匯、非線性邊界條件(如邊界上的輻射和非線性對流換熱)、具有相變的導熱和與溫度相關的材料特性等。為解決這些非線性問題，通常使用有限差分法、有限元法、有限體積法和邊界元方法等。本文介紹一種新開發(fā)的基于有限元法求解溫度相關各向同性材料熱傳導的程序FemHC。

1 非線性結構導熱問題的有限元法

描述區(qū)域中固體導熱的控制方程為

(1)

式中：為溫度；為密度；為比熱容；是熱傳導張量分量；為單位體積的內部熱生成源項；下標和滿足愛因斯坦求和規(guī)則。所有的變量都可能是空間坐標=(,,)和時間的函數(shù)。方程的求解需要指定合適的初邊值條件，即

(2)

式中：為邊界上點的坐標；為應用通量，一般給定具體的值；為對流換熱通量，

=(,,)(-)

(3)

式中：為對流換熱系數(shù)。

以上邊界條件包含固體與環(huán)境之間的熱傳導和對流傳熱。物性系數(shù)、、和可以是溫度的函數(shù)，因而各方程可以是非線性的。

方程的初始條件為

(,0)=()

(4)

初邊值問題的有限元求解分為2步：(1)空間離散，將控制方程的弱形式在典型的有限單元上進行離散，得到關于溫度的常微分方程，即得到溫度節(jié)點值的常微分方程組；(2)采用合適的方法，如有限差分法，對第一步得到的常微分方程組進行時間離散，得到關于+1時刻節(jié)點值的代數(shù)方程組。定常問題可以不考慮時間項，根據(jù)問題的非線性特點選擇合適的迭代方法求解即可。

1.1 半離散有限元模型

將對流傳熱區(qū)域離散為適合有限單元的集合,在方程兩側乘以權函數(shù)()，在單元上進行積分，對高階導數(shù)項進行分部積分，在邊界積分中應用邊界條件，得到離散元弱形式為

(5)

將溫度的有限元近似代入到弱形式中，得到半離散的有限元模型。在選擇的近似時，假設時間變化和空間變化可以分離，即

(6)

式中：為節(jié)點的溫度向量；為單元的節(jié)點數(shù)目。

將其表示為矩陣形式為

(7)

(8)

令權函數(shù)()=e,()(=1,2,…,)，并將溫度的插值函數(shù)代入到弱形式中，得

(9)

(10)

采用有限元法對控制方程和邊界條件離散并進行裝配，得到非線性常微分方程組

(11)

單元系數(shù)矩陣和向量可以表示為向量形式

(12)

本文設計軟件采用的方程形式是最一般的情形，材料物性、邊界條件和體積源項都是溫度的函數(shù)。裝配矩陣可以表示為

(13)

在式(13)中，求和在網(wǎng)格所有的單元上進行。一旦插值函數(shù)確定，單元幾何就確定，可以得到全局裝配方程為

(14)

對邊界面通量進行計算，最終得到

(15)

其中：

(16)

式中：為源項；為邊界上指定的熱流通量在邊界上的積分；為對流換熱邊界對剛度矩陣的貢獻；為對流換熱邊界對右端項的貢獻。和的表達式為

(17)

(18)

式中：為對流換熱邊界。

1.2 非線性方程組的求解

根據(jù)是否為時間相關問題，可將非線性方程組的求解分為2類：穩(wěn)態(tài)問題和瞬態(tài)問題。穩(wěn)態(tài)問題可采用Picard迭代方法求解。瞬態(tài)問題要先對時間導數(shù)項進行離散，然后對每一時間步的非線性方程組采用Picard迭代。為增加穩(wěn)定性，還可以采用松弛算法，具體細節(jié)見文獻[6]。在程序中，線性方程組采用BiCGStab算法，并采用開源線性方程組求解器Eigen求解。

2 有限元程序的實現(xiàn)

為開發(fā)有限元計算渦輪葉片結構導熱的程序，首先建立有限元計算框架。該框架具有如下特點：(1)可以統(tǒng)一處理一維、二維和三維問題；(2)能夠根據(jù)具體問題靈活指定邊界條件，將網(wǎng)格與求解過程完全解耦；(3)不限制有限單元類型，支持添加新的有限單元；(4)支持混合單元；(5)能夠方便處理多區(qū)域問題(如果計算域不同區(qū)域具有不同的物性參數(shù))和多物理場耦合問題。

2.1 統(tǒng)一處理一維、二維和三維問題

當前流行的CAE軟件，如FLUENT和OpenFOAM等，不具有全維度模擬能力，F(xiàn)LUENT只能處理二維和三維網(wǎng)格，OpenFOAM只能處理三維網(wǎng)格。OpenFOAM要計算一維和二維問題時，只能在三維網(wǎng)格上指定合適的邊界條件模擬一維和二維問題。在程序開發(fā)和問題求解的初始階段，往往要先從簡單的一維和二維問題著手，在一維和二維問題取得滿意的效果后再解決復雜的三維問題。如果算法能同時處理一維、二維和三維問題，那么可以為程序的開發(fā)和問題的求解提供很大便利。

以經(jīng)典熱傳導方程為例，說明有限元法如何統(tǒng)一處理一維、二維和三維問題。其控制方程為

(19)

根據(jù)傅里葉定律，

=-,,=

(20)

式中：為溫度；為熱通量；為熱傳導系數(shù)?？刂品匠痰馁み|金有限元離散公式為

(21)

式中：為單元形狀函數(shù)。

根據(jù)問題維度的不同，式(21)的積分對象不同：對于三維問題，其為體積分和邊界上的面積分；對于二維問題，其為面積分和邊界上的線積分;對于一維問題，其為線積分和邊界上的點積分。因此，可以根據(jù)問題維度將體網(wǎng)格單元和組成邊界的面，或者面網(wǎng)格單元和組成邊界的線，或者線網(wǎng)格單元和組成邊界的點都視為有限單元存儲。將點統(tǒng)一以三維坐標存儲，采用高斯積分公式對這些積分項進行近似計算。

以三維問題為例，對于體積分，有

(22)

對于面積分，有

(23)

在程序中建立有限元空間時，除建立體單元的有限元，還需要建立邊界上面單元的有限元；如果偏微分方程的求解算法需要考慮網(wǎng)格內面的面積分，還可以建立內面有限元。如此處理后，可以采用同樣的方法計算體積分和面積分，因而可以統(tǒng)一處理一維、二維和三維網(wǎng)格，增加程序的靈活性和適用范圍。

2.2 將網(wǎng)格與求解過程完全解耦

根據(jù)具體問題，靈活指定邊界條件，將網(wǎng)格與求解過程完全解耦。數(shù)值模擬往往需要考慮多種邊界條件類型，如果在網(wǎng)格文件中包含所求解問題的具體邊界條件，當需要改變邊界條件類型時，需要在網(wǎng)格生成軟件中更改，較為繁瑣。為將網(wǎng)格與求解過程完全解耦，本軟件設計單純網(wǎng)格文件和邊界條件文件。

2.2.1 單純計算網(wǎng)格文件

網(wǎng)格輸入依賴于3個文件，分別是網(wǎng)格節(jié)點文件node.txt、網(wǎng)格單元文件element.txt和網(wǎng)格邊界文件boundary.txt。將某二維計算域(1,5)×(1,1)均勻剖分為9個線性四邊形單元，網(wǎng)格文件截圖見圖1～3。網(wǎng)格邊界文件是在計算域劃分網(wǎng)格時得到的，與具體的問題無關，因而可將網(wǎng)格文件與求解過程解耦。

圖 1 網(wǎng)格節(jié)點文件

圖 2 網(wǎng)格單元文件

圖 3 網(wǎng)格邊界文件

2.2.2 邊界條件文件

程序求解需要根據(jù)具體問題指定相應的邊界條件，由boundaryType.txt指定，相應的文件示例截圖見圖4。對于具有多個變量的偏微分方程，不同的變量在同一邊界上一般會有不同的邊界條件。為此，每個變量都需要存儲邊界條件。存儲的邊界條件數(shù)據(jù)分為2類：(a)Dirichlet邊界，存儲的基本數(shù)據(jù)包含Dirichlet邊界的體單元編號和該體單元在Dirichlet邊界上的節(jié)點編號，可以用C++的容器map存儲；(b)需要進行面積分的邊界，存儲的數(shù)據(jù)是邊界面元對應的有限元空間的編號。

圖 4 邊界條件文件

這樣處理的優(yōu)點是可以將網(wǎng)格生成與方程求解分開，同時可以對不同變量靈活指定邊界條件，容易增加求解變量，添加新的邊界條件也容易。當然，就導熱計算而言，只需指定溫度邊界條件即可。

2.3 程序框架的核心部分

(a)節(jié)點類Node。Node類存儲網(wǎng)格節(jié)點的坐標。

(b)有限元空間類FESpace。為靈活添加有限元空間，有限元空間類設計為抽象基類。該類用于計算有限元形狀函數(shù)及其導數(shù)，存儲高斯積分點坐標和積分權重等?；惖呐缮惿删唧w的有限元空間。

(c)單元類Element及其集合ElementSet。Element類存儲單元類型、單元節(jié)點編號等。Element類還存儲指向有限元空間基類的智能指針，在建立具體單元時對有限元空間進行構造。程序不限制有限單元的類型，支持混合單元，可以根據(jù)需要添加合適的單元。ElementSet是Element的集合，包含對Element進行操作的成員函數(shù)。

(d)邊界類Boundary。Boundary類存儲邊界單元、邊界單元所在的體單元的編號和邊界條件類型。

(e)有限元節(jié)點上存儲變量編號類FemIndex。FemIndex類建立單元節(jié)點上所存儲變量的編號與單元節(jié)點編號的關系。為靈活處理多物理場耦合問題的變量存儲和調用，可以對所有的變量統(tǒng)一進行編號，也可以對某一變量單獨進行編號。

(f)方程裝配類。方程裝配類包括裝配有限元離散得到的左端項矩陣BiLinearForm和右端項向量LinearForm。這2個類被設計為抽象基類，方程中具體的左端項矩陣和右端項向量為相應基類的派生類。稀疏矩陣以壓縮行形式或壓縮列形式存儲。

圖 5 exprtk字符串解析實例

由此可見，exprtk與常規(guī)的函數(shù)表達形式十分接近，exprtk讀入這些字符串后將其解析為數(shù)學函數(shù)。核心類之間的關系見圖6。

圖 6 核心類之間的關系

3 程序驗證和應用

為驗證所采用算法及其程序實現(xiàn)的準確性，對若干算例進行計算，并與分析解或參考解進行比較，然后應用該程序對三維渦輪葉片的導熱進行數(shù)值模擬。

3.1 二維矩形條傳熱問題

以長寬比為2∶1的矩形條的傳熱問題為例，建立傳熱系統(tǒng)的二維模型。矩形內有加熱源項，左側有熱流流進，右側溫度恒定，上側施以對流換熱冷卻，下側為絕熱壁面。利用格林函數(shù)可以得到該問題的定常精確解。用于測試的相關參數(shù)如下：矩形長0.10 m、寬0.05 m，材料的熱傳導系數(shù)為0.4 W/(m·K)，體積加熱源項為=1.353×10W/m。邊界條件設定為：左側熱流通量為3 500 W/m，右側壁面給定溫度25 ℃。計算時對網(wǎng)格進行逐次細化，采用1階和2階四邊形單元進行計算，并估計有限元的精度階。左下角點、上側壁面中點和左上角點等3個測點有限元計算網(wǎng)格收斂性見表1，估計的精度階見表2，可見計算結果符合有限元的理論精度。2階四邊形單元計算得到的溫度場見圖7。

表 1 有限元計算網(wǎng)格收斂性計算結果

表 2 4節(jié)點單元有限元解的精度階估計結果

圖 7 2階四邊形單元計算得到的溫度場，K

3.2 直條的瞬態(tài)導熱

直條長度5 m，AISI 304無縫鋼管。初始溫度為300 K。直條左邊界指定溫度900 K，右側指定Neumann條件。熱擴散系數(shù)=為溫度的函數(shù)()=+，其中，=2.0×10，=0.003 7。

該問題本質上是一維問題，首先采用一維網(wǎng)格計算。時間離散采用隱式Euler格式，時間步長取1。分別采用自由度為200的1階單元和2階單元計算，得到=0.4測點處的溫度變化情況。采用二維網(wǎng)格進行計算，取直條的寬度為1，采用自由度為800的1階和2階四邊形網(wǎng)格，計算測點=04、=0.5處的溫度情況，并與一維結果進行對比，見圖8。一維計算和二維計算都與文獻[9]參考解吻合很好。

圖 8 一維和二維計算測點溫度變化曲線

3.3 三維功能梯度材料模擬

研究具有指數(shù)函數(shù)物性系數(shù)的功能梯度材料的導熱問題，計算域為[0,]的立方體。材料的導熱系數(shù)和熱容系數(shù)沿方向變化，其變化方程為

(,,)=52

(24)

(,,)=2

(25)

(0,,,)=0;(1,,,)=0;

(,0,,)=0;(,1,,)=0;

(,,0,)=0;(,,1,)=0

(26)

該問題的精確解為

(27)

(28)

采用1階六面體網(wǎng)格計算，時間離散采用隱式Euler格式，步長取0.001。不同時刻，直線(0.5,0.5,)上的溫度分布曲線見圖9，有限元計算結果與精確解吻合很好。

圖 9 直線(0.5,0.5,z)上的溫度分布曲線

3.4 具有多種材料的物體導熱計算

如果導熱物體的不同部分采用不同的導熱材料，而且不同材料的導熱系數(shù)差別很大，那么溫度計算就很困難。以建筑行業(yè)的標準算例為例，該算例計算墻橫截面的導熱，截面的長500.0 mm、寬47.5 mm。墻體由4種不同的材料組成，最大和最小導熱系數(shù)分別為230和0.029 W/(m·K)。

參照文獻[10]，多種材料物體的邊界條件和不同材料的導熱系數(shù)見圖10。上表面環(huán)境溫度為0 ℃，表面熱阻為0.06 m·K/W；下表面環(huán)境溫度為20 ℃，表面熱阻為0.11 m·K/W。表面熱阻與換熱系數(shù)的關系為

圖 10 計算域和不同區(qū)域的導熱系數(shù)和邊界條件示意,mm

(29)

式中：為表面熱阻；和分別為對流換熱系數(shù)和輻射換熱系數(shù)。忽略輻射的影響，近似可得物體上、下表面換熱系數(shù)分別為16.667和9.090 9W/(m·K)。

計算網(wǎng)格包含11 636個三角形單元、4 575個四邊形單元以及53 208個節(jié)點。溫度場分布計算結果見圖11。不同測點的溫度計算結果與參考解的對比見表3。本文計算結果與參考解吻合很好，驗證程序計算多區(qū)域問題的能力。

圖 11 具有多種材料的物體的溫度云圖，℃

表 3 有限元計算的測點溫度與參考解對比

3.5 不含冷卻流道的三維葉片的穩(wěn)態(tài)導熱計算

某不含冷卻流道的三維葉片及其不同面上的邊界條件設定見圖12，其計算網(wǎng)格見圖13，含有67 617個四面體單元、15 802個節(jié)點。進行線性導熱計算，葉片的導熱系數(shù)為12 W/(m·K)。計算得到的葉片表面溫度云圖見圖14，葉片內部切片溫度云圖見圖15。

圖 12 三維葉片導熱邊界條件

圖 13 三維葉片導熱計算網(wǎng)格

圖 14 三維葉片溫度云圖，K

圖 15 葉片切片上的溫度云圖，K

3.6 含冷卻流道的三維葉片的穩(wěn)態(tài)導熱計算

選取某含9個冷卻流道的三維葉片，計算網(wǎng)格見圖16，具有436 952個四面體單元、90 021個節(jié)點。葉片的導熱系數(shù)是關于溫度的分段線性函數(shù)，作為算例測試的導熱系數(shù)函數(shù)為

(30)

圖 16 帶冷卻流道的渦輪葉片計算網(wǎng)格

葉片葉身以及每個冷卻流道指定對流換熱邊界條件，換熱系數(shù)和換熱溫度各不相同。葉身上的換熱溫度和換熱系數(shù)分別為1 700 K和200 W/(m·K)，冷卻流道1～10表面的換熱溫度依次為400、500、600、700、800、900、950、1 000、1 050和1 100 K，冷卻流道1～10的換熱系數(shù)依次為1 000、1 100、1 200、1 300、1 400、1 500、1 600、1 700、1 800和1 900 W/(m·K)，在其余邊界上指定絕熱邊界條件。

計算得到葉片的溫度云圖見圖17。計算在GNU/Linux系統(tǒng)上進行，處理器為AMD Ryzen5 3550H，其主頻為2.10 GHz，當殘值小于1×10時計算終止，所用時間為112 s。本例說明程序指定復雜物性參數(shù)和處理任意多個邊界條件的能力。

圖 17 帶冷卻流道的渦輪葉片溫度云圖，K

4 結束語

介紹有限元非線性導熱計算程序FemHC所采用的算法和具體的實現(xiàn)細節(jié)。FemHC通過將網(wǎng)格與具體的問題解耦，可以靈活指定邊界條件，通過字符串解析庫實現(xiàn)物性參數(shù)、初始邊值條件、源項函數(shù)和計算控制參數(shù)從外部文本文件中讀取，大大增加程序的靈活性。通過若干算例驗證FemHC在計算線性/非線性，穩(wěn)態(tài)/瞬態(tài)計算中的準確性。通過三維渦輪葉片的導熱計算，表明FemHC可以用于實際渦輪葉片的計算。