宋小鵬 何秀錦 古小敏 吳國珊 王斌武

摘要：

針對復(fù)雜計算域有限元分析前處理過程中手動建模和邊界條件設(shè)定煩瑣和耗時的問題，以蓄熱磚溫度場分析為例，討論Lua腳本在有限元前處理過程中的應(yīng)用;開發(fā)溫度場有限元分析程序，基于Lua/C API定義若干接口，使用特例的解析解驗證溫度場計算結(jié)果，基于Lua腳本對蓄熱磚溫度場進行計算。仿真結(jié)果表明，Lua腳本可簡化有限元前處理過程，特別是幾何建模過程，也便于對計算結(jié)果進行參數(shù)化分析。

關(guān)鍵詞：

蓄熱磚; 溫度場; 前處理; 有限元

中圖分類號：? TU522.07; TB115.1

文獻標(biāo)志碼：? B

Application of Lua script in preprocessing

of finite element analysis

SONG Xiaopeng， HE Xiujin， GU Xiaoming， WU Guoshan， WANG Binwu

（

School of Energy and Building Environment， Guilin University of Aerospace Technology，

Guilin 541004， Guangxi， China）

Abstract：

As to the issue that the manual modeling and boundary condition setting are tedious and timeconsuming in the preprocessing of finite element analysis in complex computing domain，? the application of Lua script in the preprocessing of finite element is studied

by taking temperature field analysis of thermal storage brick as an example. The finite element analysis program for temperature field of thermal storage brick is developed. The interfaces based on Lua/C API is defined， and the calculation results of temperature field are verified by the analytic solution of a special case. The temperature field of thermal storage brick is calculated based on Lua script. The simulation results show that the preprocessing of finite element， especially the geometric modeling process is simplified using the Lua script， and the parametric analysis on the computing results is easy to be completed.

Key words：

thermal storage brick; temperature field; preprocessing; finite element

0 引 言

Lua腳本程序具有跨平臺性能好、執(zhí)行速度快、支持熱更新、可擴展性能好等優(yōu)勢，常用于程序配置。[12]陸喬[3]在不更新應(yīng)用程序的前提下，基于Lua腳本實現(xiàn)動態(tài)增加程序功能的目標(biāo)。程君[4]以Lua腳本作為移動互聯(lián)網(wǎng)中間件，解決不同平臺的編程接口差異導(dǎo)致信息傳遞不暢的問題。宋宏江等[5]基于Lua腳本編制驗證程序，明顯提升衛(wèi)星綜合測試過程中指令生成的效率。

大量科研工作者已對蓄熱磚的蓄熱特性進行分析，如：趙頔等[6]使用ANSYS軟件分析蓄熱磚蓄熱特性，分析結(jié)果與實驗研究結(jié)果較為一致;徐德璽等[7]基于有限元法計算蓄熱磚孔道直徑與其熱效率之間的關(guān)系。蓄熱磚的孔道形狀和布置對其蓄熱特性有顯著影響：胡思科等[8]研究圓形與方形孔道固體蓄熱磚的蓄熱和放熱特性;李晶晶等[9]研究固體蓄熱裝置的自然對流放熱特性受圓形和橢圓形孔道布置的影響。

有限元分析是研究蓄熱材料溫度場的重要手段。在蓄熱磚有限元分析過程中，當(dāng)孔道數(shù)量較多時，孔道的建模和邊界條件設(shè)置較為煩瑣。本文以蓄熱磚溫度場分析為例，探討Lua腳本在有限元分析過程中的應(yīng)用，旨在提高有限元建模和計算效率。

1 蓄熱磚數(shù)學(xué)物理模型

1.1 物理模型

計算采用的蓄熱磚結(jié)構(gòu)示意見圖1，其底面寬為0.4 m，底面長為1.0 m，高為0.2 m。在蓄熱磚長度方向上布置若干孔道，用于蓄熱后釋放顯熱，其中孔道半徑為14 mm，孔道中心間距為45 mm。蓄熱磚底面通過電阻絲加熱，以增加蓄熱磚顯熱。

1.2 數(shù)學(xué)模型

1.2.1 控制偏微分方程

蓄熱磚沿孔道方向較長，在蓄熱過程中熱量的主要傳遞方向為由下底面到上底面。為減少計算量，可只計算蓄熱磚橫截面上的溫度分布，故其計算域見圖2。

在蓄熱磚蓄熱過程中，孔道內(nèi)無流體的強迫流動，可僅考慮熱傳導(dǎo)，因此其控制偏微分方程為

式中：T為蓄熱磚溫度;t為蓄熱時間;ρ、Cp和λ分別為蓄熱磚的密度、定壓熱容和熱導(dǎo)率。

1.2.2 定解條件

蓄熱磚初始溫度為20 ℃，下底面由電阻絲提供

恒定熱流加熱，熱流密度設(shè)定為2 500 W/m2，連續(xù)加熱10 h，假設(shè)蓄熱磚的其他邊和孔道內(nèi)部均為絕熱邊界條件。蓄熱磚主要材料為氧化鎂，假設(shè)加熱過程物性參數(shù)恒定，其密度為3 000 kg/m3，定壓熱容為1 kJ/（kg·K），熱導(dǎo)率為3 W/（m·K）。

2 溫度場有限元求解過程

2.1 求解過程

在有限元計算過程中，先結(jié)合邊界條件對計算域進行網(wǎng)格劃分，將計算域離散為有限數(shù)量的三角形單元，根據(jù)單元形狀函數(shù)、邊界條件和物性參數(shù)等計算各個三角形單元的單元系數(shù)矩陣。使用伽遼金加權(quán)余量法計算穩(wěn)態(tài)溫度場并合成全局矩陣，得到關(guān)于溫度場的代數(shù)方程。代數(shù)方程求解完成后，對計算結(jié)果進行后處理。

2.2 Lua腳本封裝

使用C+編寫二維溫度場有限元計算程序，編寫5類關(guān)鍵接口，實現(xiàn)幾何模型計算域的設(shè)定、定解條件（邊界條件、初始條件、材料屬性等）的定義、網(wǎng)格劃分、代數(shù)方程求解和結(jié)果輸出等。將C++程序編譯為可執(zhí)行程序，使用配置文件設(shè)定程序執(zhí)行過程中的參數(shù)。Lua腳本解釋器體積小、執(zhí)行速度快，可用于程序配置文件。與常規(guī)配置文件相比，Lua腳本配置文件可以執(zhí)行各類復(fù)雜數(shù)學(xué)運算，支持循環(huán)語句、條件轉(zhuǎn)移等常規(guī)語言具有的特性。本文使用Lua配置有限元程序，在C++中通過Lua/C API定義9個關(guān)鍵內(nèi)置程序接口（見圖3），用于Lua腳本調(diào)用。通過編寫Lua腳本，實現(xiàn)通用二維溫度場有限元分析程序的幾何模型計算域定義、定解條件設(shè)定、網(wǎng)格劃分、計算和結(jié)果輸出。

3 計算結(jié)果與討論

3.1 程序驗證

為評估有限元程序計算結(jié)果的準(zhǔn)確性，在矩形

假設(shè)邊界條件為T|x=0=Ay（b-y）、T|x=a=0、T|y=0=Bsin（πx/a）和T|y=b=0，根據(jù)數(shù)學(xué)物理方法中的分離變量法[10]，可得到其解析解為

編寫Lua腳本以求解該溫度場。定義邊界條件、材料熱物性參數(shù)和幾何模型計算域的Lua腳本如下。

1. a，b=50，30 矩形計算域長度和寬度

2. addBC（1，0，0，"10*sin（3.1415926*x/50）"） 左側(cè)邊的邊界條件 #01

3. addBC（1，0，0，"0.05*y*（30y）"） 下側(cè)邊的邊界條件 #02

4. addBC（1，0，0，""） 上側(cè)和右側(cè)都是絕熱邊界條件 #03

5. addMaterial（34，7200，680，""） 定義材料熱物性：導(dǎo)熱系數(shù)，密度和比熱容

6. addNode（0，0）計算域左下節(jié)點

7. addNode（a，0） 計算域右下節(jié)點

8. addNode（a，b） 計算域右上節(jié)點

9. addNode（0，b） 計算域左上節(jié)點

10. addLineSgmt（0，1，1，1） 計算域下底邊邊界條件為#01

11. addLineSgmt（1，2，1，3） 計算域下底邊邊界條件為#02

12. addLineSgmt（2，3，1，3） 計算域下底邊邊界條件為#03

13. addLineSgmt（3，0，1，2） 計算域下底邊邊界條件為#02

14. addDomain（25，15，1，10） 設(shè)定計算域材料為#01，最大三角形單元面積為10

15. iterations，timsStep，initialValue，coordnate=1，0，0，0? 迭代次數(shù)、步長和初始值

16. Solve（iterations，timsStep，initialValue，coordnate） 計算

17. Export（"luaResult.dat"，0）導(dǎo)出計算結(jié)果

溫度場計算的解析解和使用上述Lua腳本求解得到結(jié)果對比見圖4。由此可見，計算結(jié)果誤差較小，說明本文有限元計算程序能夠求解得到合理的數(shù)值解。

C++程序也可以通過Lua/C API調(diào)用Lua腳本中預(yù)定義的函數(shù)，圖4溫度場左側(cè)的邊界條件定義如下。

1. function LuaFunc_BC_left（x， y） 左側(cè)邊的值

2.? ?return 10*sin（3.1415926*x/50）

3. end

4. addBC（1，0，0，"LuaFunc_BCleft"） 左側(cè)邊的邊界條件 #01

通過調(diào)用Lua腳本中邊界條件設(shè)定的函數(shù)，完成邊界條件設(shè)置，但應(yīng)當(dāng)避免每次迭代過程都調(diào)用Lua腳本，以保證計算速度。

3.2 蓄熱磚溫度場分析

結(jié)合在C++中設(shè)定的且可在Lua腳本中被調(diào)用的函數(shù)，編寫蓄熱磚溫度場求解腳本，求解步驟如下：（1）定義下底邊的熱流邊界條件和其他邊的絕熱邊界條件;（2）定義蓄熱磚材料的熱物性參數(shù);（3）定義蓄熱磚的4個頂點;（4）定義蓄熱磚計算域和4個邊的邊界條件，定義蓄熱磚內(nèi)的孔道（這些孔道要在計算域中去除，每個孔道需要定義節(jié)點、弧線段并設(shè)置絕熱邊界條件）;（5）定義計算域及其材料，設(shè)定網(wǎng)格劃分最大面積閾值;（6）求解代數(shù)方程;（7）輸出計算結(jié)果。具體代碼如下。

1. print（"Step1： 定義邊界條件"）

2. addBC（2，2500，0，""）熱流邊界條件 #01

3. addBC（2，0，0，""）絕熱邊界條件 #02

4. print（"Step2： 定義材料"）

5. addMaterial（5，3000，1000，""）氧化鎂熱物性

6. print（"Step3： 定義節(jié)點"）

7. addNode（0，0）Node #00 左下節(jié)點

8. addNode（0.4，0）Node #01 右下節(jié)點

9. addNode（0.4，0.2）Node #02 右上節(jié)點

10. addNode（0，0.2）Node #03 左上節(jié)點

11. 孔道建模

12. nodeID，r，dis，x0，y0=4，0.014，0.045，0.04，0.035孔道半徑、間距等

13. for j=0，3，1 do垂直方向4行孔道

14.? ?for i=0，7，1 do水平方向8列孔道

15.? ?addNode（x0+i*disr，y0+j*dis）孔道直徑節(jié)點1

16.? ? addNode（x0+i*dis+r，y0+j*dis）孔道直徑節(jié)點2

17.? ? print（"Step4.1：添加弧線段"）

18.? ? addArcSgmt（nodeID，nodeID+1，180，15，2）孔道上弧段，剖分15段，采用#2絕熱條件

19.? ? addArcSgmt（nodeID+1，nodeID，180，15，2）孔道下弧段，剖分15段，采用#2絕熱條件

20.? ? print（"Step4.2： 定義計算域內(nèi)的孔道"）

21.? ? addHole（x0+i*dis，y0+j*dis）設(shè)定孔洞區(qū)域

22.? ? nodeID=nodeID+2

23.? ?end

24. end

25. print（"Step4.3：Add 蓄熱磚4條邊"）

26. addLineSgmt（0，1，60，1）;下底邊，剖分60段，邊界條件#01 熱流邊界條件

27. addLineSgmt（1，2，30，2）右側(cè)邊，剖分30段，邊界條件#02 絕熱

28. addLineSgmt（2，3，60，2）上底邊，剖分60段，邊界條件#02 絕熱

29. addLineSgmt（3，0，30，2）左側(cè)邊，剖分30段，邊界條件#02 絕熱

30. print（"Step5： 設(shè)定計算域及其材料，網(wǎng)格最大面積"）

31. materialID，maxAreaSize=1，0.00005

32. addDomain（0.001，001，materialID，maxAreaSize）計算域材料為前述定義#01

33. print（"Step6： 求解"）

34. iterations，timsStep，initialValue，coordnate=60*60*10，1，20，0迭代次數(shù)，時間步長，初始值

35. Solve（iterations，timsStep，initialValue，coordnate）

36. print（"—Step7： 導(dǎo)出計算結(jié)果"）

37. Export（"luaResult.dat"，0）

38. print（"Soulution Done！"）

使用上述Lua腳本生成的蓄熱磚網(wǎng)格見圖5，共計2 931個節(jié)點、4 784個三角形單元。蓄熱磚寬面被剖分為60段，窄面被剖分為30段，孔道被剖分為30段。通過Lua語句的兩重循環(huán)，實現(xiàn)孔道幾何尺寸設(shè)定、孔道布置、孔道網(wǎng)格分段數(shù)設(shè)定和孔道邊界條件的設(shè)定，可極大簡化有限元分析的前處理過程。如果使用常規(guī)商業(yè)有限元分析軟件建模，每個孔洞設(shè)置一次邊界條件，則需要設(shè)定32次邊界條件，而通過Lua腳本批量建模并設(shè)定邊界條件，前處理操作時間明顯縮短。

使用電阻絲加熱蓄熱磚10 h后，其溫度場的計算結(jié)果見圖6。蓄熱磚左右對稱，邊界條件也對稱，因此溫度場也呈現(xiàn)左右對稱特點。下底面（熱面）為熱流邊界條件，其溫度分布最高，上底面（冷面）溫度分布最低。蓄熱磚溫度隨高度增加呈非線性減小;越靠近熱面，溫度場變化越劇烈?？椎栏浇鼫囟葓鲎兓鄬×?，但孔道附近的等溫線都與孔道正交，這是因為孔道為絕熱邊界條件，所以無熱流通過。

如果需要改變蓄熱磚孔道直徑、排列間距等幾何特征參數(shù)，只需要修改Lua腳本中的參數(shù)即可，無須重新手工建模和設(shè)定邊界條件，便于對計算進行參數(shù)化分析。

分別繪制圖6中A點（蓄熱磚最底部兩個相鄰孔道中心）、B點（蓄熱磚幾何中心）和C點（蓄熱磚上部兩個相鄰孔道中心）的溫度變化曲線，見圖7。加熱2 h后，3個點的溫度均呈線性增加。由于C點距電阻絲最遠，其溫度在加熱早期（2 h內(nèi)）升溫緩慢;A點距電阻絲較近，能快速升溫;B點到A點和到C點的距離相近，但B點溫度并非A點與C點溫度的平均值，說明蓄熱磚加熱過程是非穩(wěn)態(tài)導(dǎo)熱。

4 結(jié)束語

探究使用Lua腳本簡化有限元前處理建模過程，得到如下結(jié)論。

（1）使用C++開發(fā)二維溫度場有限元分析程序，計算結(jié)果與解析解吻合較好，說明建立的有限元分析程序算法準(zhǔn)確。

（2）在蓄熱磚溫度場有限元分析過程中，通過Lua/C API編寫Lua腳本可簡化有限元前處理過程，便于對有限元計算結(jié)果進行參數(shù)化分析。

參考文獻：

[1]

陸艮峰， 張蓉， 張賽橋. Lua在軌道交通綜合監(jiān)控系統(tǒng)模擬器中的應(yīng)用[J]. 工業(yè)控制計算機， 2017， 30（3）： 1315. DOI： 10.3969/j.issn.1001182X.2017.03.006.

[2] 陳佳銘， 王風(fēng)立， 鄧君湘， 等. Lua腳本與C++交互流程及在其HSTPN仿真軟件的應(yīng)用[J]. 計算機應(yīng)用與軟件， 2018， 35（10）：1316. DOI： 10.3969/j.issn.1000386x.2018.10.003.

[3] 陸喬. 基于Lua的iOS動態(tài)化系統(tǒng)的設(shè)計與實現(xiàn)[D]. 武漢： 華中科技大學(xué)， 2019.

[4] 程君. 基于Lua的移動互聯(lián)網(wǎng)中間件系統(tǒng)的研究與實現(xiàn)[D]. 南京： 東南大學(xué)， 2017.

[5] 宋宏江， 高何， 盧成志， 等. 一種基于Lua腳本的航天器遙控快速測試驗證技術(shù)[J]. 航天器工程， 2020， 29（3）： 182186. DOI： 10.3969/j.issn.16738748.2020.03.028.

[6] 趙頔， 王啟民. 基于ANSYS分析的蓄熱磚蓄熱特性數(shù)值模擬及實驗研究[J]. 沈陽工程學(xué)院學(xué)報（自然科學(xué)版）， 2020， 16（2）： 3438. DOI： 10.13888/j.cnki.jsie（ns）.2020.02.008.

[7] 徐德璽， 金映麗， 邢作霞， 等. 基于有限元的固體電蓄熱裝置蓄熱模擬及實驗[J]. 機械工程與自動化， 2016（4）： 2325. DOI： 10.3969/j.issn.16726413.2016.04.009.

[8] 胡思科， 劉建宇， 邢姣嬌. 具有圓、方孔道的固體蓄、放熱特性的分析與比較[J]. 流體機械， 2015， 43（9）： 7378. DOI： 10.3969/j.issn.10050329.2015.09.015.

[9] 李晶晶， 李永光， 王治源， 等. 固體蓄熱裝置圓形和橢圓形孔自然對流放熱特性的研究[J]. 上海電力學(xué)院學(xué)報， 2019， 35（6）： 525530. DOI： 10.3969/j.issn.10064729.2019.06.003.

[10] 梁昆淼. 數(shù)學(xué)物理方法[M]. 4版. 北京： 高等教育出版社， 2009： 143162.

（編輯 武曉英）