三維有限元非線性導(dǎo)熱計(jì)算程序FemHC

劉瑜， 朱可可， 萬(wàn)秋里， 王成恩

(1.上海交通大學(xué) 機(jī)械與動(dòng)力工程學(xué)院，上海 200240； 2.中國(guó)航空發(fā)動(dòng)機(jī)集團(tuán)有限公司沈陽(yáng)發(fā)動(dòng)機(jī)研究所，沈陽(yáng) 110015)

0 引 言

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

1 非線性結(jié)構(gòu)導(dǎo)熱問(wèn)題的有限元法

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

(1)

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

(2)

式中：為邊界上點(diǎn)的坐標(biāo)；為應(yīng)用通量，一般給定具體的值；為對(duì)流換熱通量，

=(,,)(-)

(3)

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

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

方程的初始條件為

(,0)=()

(4)

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

1.1 半離散有限元模型

將對(duì)流傳熱區(qū)域離散為適合有限單元的集合,在方程兩側(cè)乘以權(quán)函數(shù)()，在單元上進(jìn)行積分，對(duì)高階導(dǎo)數(shù)項(xiàng)進(jìn)行分部積分，在邊界積分中應(yīng)用邊界條件，得到離散元弱形式為

(5)

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

(6)

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

將其表示為矩陣形式為

(7)

(8)

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

(9)

(10)

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

(11)

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

(12)

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

(13)

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

(14)

對(duì)邊界面通量進(jìn)行計(jì)算，最終得到

(15)

其中：

(16)

式中：為源項(xiàng)；為邊界上指定的熱流通量在邊界上的積分；為對(duì)流換熱邊界對(duì)剛度矩陣的貢獻(xiàn)；為對(duì)流換熱邊界對(duì)右端項(xiàng)的貢獻(xiàn)。和的表達(dá)式為

(17)

(18)

式中：為對(duì)流換熱邊界。

1.2 非線性方程組的求解

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

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

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

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

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

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

(19)

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

=-,,=

(20)

式中：為溫度；為熱通量；為熱傳導(dǎo)系數(shù)。控制方程的伽遼金有限元離散公式為

(21)

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

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

以三維問(wèn)題為例，對(duì)于體積分，有

(22)

對(duì)于面積分，有

(23)

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

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

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

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

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

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

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

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

2.2.2 邊界條件文件

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

圖 4 邊界條件文件

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

2.3 程序框架的核心部分

(a)節(jié)點(diǎn)類Node。Node類存儲(chǔ)網(wǎng)格節(jié)點(diǎn)的坐標(biāo)。

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

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

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

(e)有限元節(jié)點(diǎn)上存儲(chǔ)變量編號(hào)類FemIndex。FemIndex類建立單元節(jié)點(diǎn)上所存儲(chǔ)變量的編號(hào)與單元節(jié)點(diǎn)編號(hào)的關(guān)系。為靈活處理多物理場(chǎng)耦合問(wèn)題的變量存儲(chǔ)和調(diào)用，可以對(duì)所有的變量統(tǒng)一進(jìn)行編號(hào)，也可以對(duì)某一變量單獨(dú)進(jìn)行編號(hào)。

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

圖 5 exprtk字符串解析實(shí)例

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

圖 6 核心類之間的關(guān)系

3 程序驗(yàn)證和應(yīng)用

為驗(yàn)證所采用算法及其程序?qū)崿F(xiàn)的準(zhǔn)確性，對(duì)若干算例進(jìn)行計(jì)算，并與分析解或參考解進(jìn)行比較，然后應(yīng)用該程序?qū)θS渦輪葉片的導(dǎo)熱進(jìn)行數(shù)值模擬。

3.1 二維矩形條傳熱問(wèn)題

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

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

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

圖 7 2階四邊形單元計(jì)算得到的溫度場(chǎng)，K

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

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

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

圖 8 一維和二維計(jì)算測(cè)點(diǎn)溫度變化曲線

3.3 三維功能梯度材料模擬

研究具有指數(shù)函數(shù)物性系數(shù)的功能梯度材料的導(dǎo)熱問(wèn)題，計(jì)算域?yàn)閇0,]的立方體。材料的導(dǎo)熱系數(shù)和熱容系數(shù)沿方向變化，其變化方程為

(,,)=52

(24)

(,,)=2

(25)

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

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

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

(26)

該問(wèn)題的精確解為

(27)

(28)

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

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

3.4 具有多種材料的物體導(dǎo)熱計(jì)算

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

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

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

(29)

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

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

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

表 3 有限元計(jì)算的測(cè)點(diǎn)溫度與參考解對(duì)比

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

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

圖 12 三維葉片導(dǎo)熱邊界條件

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

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

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

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

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

(30)

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

葉片葉身以及每個(gè)冷卻流道指定對(duì)流換熱邊界條件，換熱系數(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)，在其余邊界上指定絕熱邊界條件。

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

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

4 結(jié)束語(yǔ)

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