Lua腳本在有限元分析中的應(yīng)用

宋小鵬 何秀錦 古小敏 吳國(guó)珊 王斌武

摘要：

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

關(guān)鍵詞：

蓄熱磚; 溫度場(chǎng); 前處理; 有限元

中圖分類號(hào)：? TU522.07; TB115.1

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

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

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

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

1.1 物理模型

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

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

1.2.1 控制偏微分方程

蓄熱磚沿孔道方向較長(zhǎng)，在蓄熱過(guò)程中熱量的主要傳遞方向?yàn)橛上碌酌娴缴系酌妗闇p少計(jì)算量，可只計(jì)算蓄熱磚橫截面上的溫度分布，故其計(jì)算域見(jiàn)圖2。

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

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

1.2.2 定解條件

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

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

2 溫度場(chǎng)有限元求解過(guò)程

2.1 求解過(guò)程

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

2.2 Lua腳本封裝

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

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

3.1 程序驗(yàn)證

為評(píng)估有限元程序計(jì)算結(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]，可得到其解析解為

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

1. a，b=50，30 矩形計(jì)算域長(zhǎng)度和寬度

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）計(jì)算域左下節(jié)點(diǎn)

7. addNode（a，0） 計(jì)算域右下節(jié)點(diǎn)

8. addNode（a，b） 計(jì)算域右上節(jié)點(diǎn)

9. addNode（0，b） 計(jì)算域左上節(jié)點(diǎn)

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

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

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

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

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

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

16. Solve（iterations，timsStep，initialValue，coordnate） 計(jì)算

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

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

C++程序也可以通過(guò)Lua/C API調(diào)用Lua腳本中預(yù)定義的函數(shù)，圖4溫度場(chǎng)左側(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

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

3.2 蓄熱磚溫度場(chǎng)分析

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

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

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

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

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

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é)點(diǎn)1

16.? ? addNode（x0+i*dis+r，y0+j*dis）孔道直徑節(jié)點(diǎn)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： 定義計(jì)算域內(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è)定計(jì)算域及其材料，網(wǎng)格最大面積"）

31. materialID，maxAreaSize=1，0.00005

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

33. print（"Step6： 求解"）

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

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

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

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

38. print（"Soulution Done！"）

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

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

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

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

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

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

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

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

參考文獻(xiàn)：

[1]

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

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

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

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

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

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

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

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

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

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

（編輯 武曉英）