Lua腳本在有限元分析中的應用

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

摘要：

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

關鍵詞：

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

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

文獻標志碼：? 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腳本程序具有跨平臺性能好、執行速度快、支持熱更新、可擴展性能好等優勢，常用于程序配置。[12]陸喬[3]在不更新應用程序的前提下，基于Lua腳本實現動態增加程序功能的目標。程君[4]以Lua腳本作為移動互聯網中間件，解決不同平臺的編程接口差異導致信息傳遞不暢的問題。宋宏江等[5]基于Lua腳本編制驗證程序，明顯提升衛星綜合測試過程中指令生成的效率。

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

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

1 蓄熱磚數學物理模型

1.1 物理模型

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

1.2 數學模型

1.2.1 控制偏微分方程

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

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

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

1.2.2 定解條件

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

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

2 溫度場有限元求解過程

2.1 求解過程

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

2.2 Lua腳本封裝

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

3 計算結果與討論

3.1 程序驗證

為評估有限元程序計算結果的準確性，在矩形

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

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

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

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

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

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

5. addMaterial（34，7200，680，""） 定義材料熱物性：導熱系數，密度和比熱容

6. addNode（0，0）計算域左下節點

7. addNode（a，0） 計算域右下節點

8. addNode（a，b） 計算域右上節點

9. addNode（0，b） 計算域左上節點

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） 設定計算域材料為#01，最大三角形單元面積為10

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

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

17. Export（"luaResult.dat"，0）導出計算結果

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

C++程序也可以通過Lua/C API調用Lua腳本中預定義的函數，圖4溫度場左側的邊界條件定義如下。

1. function LuaFunc_BC_left（x， y） 左側邊的值

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

3. end

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

通過調用Lua腳本中邊界條件設定的函數，完成邊界條件設置，但應當避免每次迭代過程都調用Lua腳本，以保證計算速度。

3.2 蓄熱磚溫度場分析

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

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： 定義節點"）

7. addNode（0，0）Node #00 左下節點

8. addNode（0.4，0）Node #01 右下節點

9. addNode（0.4，0.2）Node #02 右上節點

10. addNode（0，0.2）Node #03 左上節點

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）孔道直徑節點1

16.? ? addNode（x0+i*dis+r，y0+j*dis）孔道直徑節點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： 定義計算域內的孔道"）

21.? ? addHole（x0+i*dis，y0+j*dis）設定孔洞區域

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）右側邊，剖分30段，邊界條件#02 絕熱

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

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

30. print（"Step5： 設定計算域及其材料，網格最大面積"）

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迭代次數，時間步長，初始值

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

36. print（"—Step7： 導出計算結果"）

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

38. print（"Soulution Done！"）

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

使用電阻絲加熱蓄熱磚10 h后，其溫度場的計算結果見圖6。蓄熱磚左右對稱，邊界條件也對稱，因此溫度場也呈現左右對稱特點。下底面（熱面）為熱流邊界條件，其溫度分布最高，上底面（冷面）溫度分布最低。蓄熱磚溫度隨高度增加呈非線性減小;越靠近熱面，溫度場變化越劇烈。孔道附近溫度場變化相對劇烈，但孔道附近的等溫線都與孔道正交，這是因為孔道為絕熱邊界條件，所以無熱流通過。

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

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

4 結束語

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

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

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

參考文獻：

[1]

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

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

[3] 陸喬. 基于Lua的iOS動態化系統的設計與實現[D]. 武漢： 華中科技大學， 2019.

[4] 程君. 基于Lua的移動互聯網中間件系統的研究與實現[D]. 南京： 東南大學， 2017.

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

[6] 趙頔， 王啟民. 基于ANSYS分析的蓄熱磚蓄熱特性數值模擬及實驗研究[J]. 沈陽工程學院學報（自然科學版）， 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]. 上海電力學院學報， 2019， 35（6）： 525530. DOI： 10.3969/j.issn.10064729.2019.06.003.

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

（編輯 武曉英）