在線算法演示虛擬仿真平臺開發
——以熱傳導為例

宋小鵬， 李志偉， 謝云云

（桂林航天工業學院能源與建筑環境學院，廣西 桂林 541004）

0 引言

隨著信息技術的發展和教學改革工作的推進，當實驗操作困難或可視條件不佳時，虛擬仿真技術為教學提供了諸多便利［1-2］。教育部在2020 年發布了虛擬仿真實驗教學課程建設指南，旨在更好地應用虛擬仿真技術服務教學工作［3-4］。文獻［5-6］中分別利用虛擬仿真技術開展發電廠實習、實踐教學，并取得了較好的效果。

傳熱學中一些常規實驗不夠直觀，無法直接觀察到溫度場、等溫線、熱流矢量分布等。田昌［7］針對“傳熱學”實驗開展虛擬仿真教學，探討和發揮其優勢。侴愛輝等［8］探討傳熱學實驗虛擬仿真和實驗結合的教學方式，加強學生對課程內容的直觀理解。大量教學工作者，基于Matlab開發傳熱學實驗虛擬仿真［9-11］。Matlab功能齊全，需使用許可，部分高校使用受限。梁秀俊等［12］基于Java 語言，王輝等［13］基于Fluent，李大鵬等［14］基于Energy2D，分別開發了傳熱學實驗虛擬仿真。Java語言、Fluent和Energy2D 等軟件具有較高的學習成本，需要較長學習時間，不利于快速入門。

鑒于此，使用Web 前端技術（HTML5／JavaScript）開發算法演示平臺，以便對仿真結果進行可視化前處理或后處理。

1 仿真平臺的設計與實現

熱傳導數值求解實驗虛擬仿真的常規流程為搭建開發環境（如安裝C／C ++ ／Fortran 編譯器和IDE），程序編制，編譯鏈接和運行，對運行結果后處理、結果分析與程序調試，如圖1 所示。為縮短算法演示流程，將實驗集中在算法編程的實現上。使用JavaScript 作為編程語言，借助瀏覽器便能對仿真結果實時可視化。

圖1 短流程算法演示

基于Vue-repl、Element plus playground、MathJax和ECharts等開源項目，搭建交互式的實驗虛擬仿真平臺。Element plus能快速在網頁中構建用戶界面，提高交互性；MathJax 支持在網頁中展示Latex 數學公式；ECharts為數據可視化圖表庫，可在網頁中對仿真結果可視化，如繪制溫度曲線或者溫度場的Contour 圖。仿真平臺界面如圖2 所示，主要由代碼編輯器和計算程序執行結果預覽區構成。使用基于vue3.2 的模板（通常包含了用戶界面、程序邏輯腳本和樣式）編寫仿真程序，可對代碼關鍵字進行高亮。當程序修改時，系統可自動編譯vue 文件并運行（基于Vue-repl），程序無誤時在預覽區實時顯示仿真結果。

圖2 仿真平臺界面

用戶編輯程序過程中，系統自動檢查語法錯誤，如有錯誤則顯示之，如無語法錯誤，則對程序進行預覽運行，如圖3 所示。通常仿真程序運行后，顯示用戶界面，當用戶與頁面交互時（如仿真程序的運行參數調整），根據用戶輸入進行網格剖分、邊界條件設置、代數方程的迭代計算以及對結果的可視化展示。仿真程序的運行和數據可視化均在瀏覽器端進行，服務器端負荷很低。

圖3 仿真平臺程序邏輯

2 熱傳導實驗數值解虛擬仿真算法

傳熱是典型的擴散傳輸現象之一，單純從偏微分方程理解溫度場分布較為困難和抽象，將溫度場可視化，能直觀理解和掌握“場”的概念。基于前述的仿真平臺，開發熱傳導數值仿真程序，將溫度場可視化。

2.1 計算域內部節點迭代格式說明

二維矩形計算域的非穩態溫度場控制偏微分方程

式中：ρ、CP、λ和T分別為控制體的密度、等壓比熱、導熱系數和溫度；為單位體積內熱源功率。考慮到控制體內可能有內熱源，若內熱源功率隨節點溫度線性變化，則

以規則矩形計算域的正交網格為例，如圖4 所示，與節點P直接相鄰的有東、南、西、北4 個節點，分別為E、S、W和N，節點P所在控制體大小為Δx·Δy·1（設控制體另一維度為1），節點P所在控制體（形狀為6 面體）在東、南、西、北4 個面的面積分別為Se、Ss、Sw和Sn，顯然Se＝Sw＝Δy·1，Sn＝Ss＝Δx·1。

圖4 網格及節點控制體示意圖

為使迭代計算有明確的物理意義，使用能量平衡法推導面心結構單元的溫度場迭代公式。以矩形網格為例，假設計算域內部控制體P僅和與其直接相鄰的控制體E、S、W和N 有熱量交換，控制體P與四周交換的熱量

式中：nb下標為與控制體P直接相鄰的4 個控制體；TP為控制體節點P的溫度；Tnb為相鄰節點的溫度；Snb為與控制體P相鄰控制體界面的面積。控制體界面上的凈傳熱量與內熱源發熱量共同決定了控制體內所包含的熱量，故控制體溫度

式中：V和Δt分別為控制體的體積和時間步長；上標0為上一時刻的值；knb為相鄰控制體界面間的傳熱系數，兩相鄰控制體換熱方式為導熱時knb＝λnb／dnb，λnb為界面處的導熱系數，dnb為節點P到相鄰控制體中心點的距離，dnb在東南西北方向上的長度分別為圖中的dxPE、dxPW、dyPS和dyPN。特別地，當東西方向網格均勻剖分時Δx＝dxPE＝dxPW；當南北方向網格均勻剖分時Δy＝dyPS＝dyPN。如為穩態導熱問題，則左側為0。若各物性參數為常數，網格均勻，穩態導熱且無內熱源時，可得節點P的顯式

2.2 邊界節點迭代格式

邊界節點的溫度求解直接影響整個溫度場的求解。第1 類邊界條件時，即計算域邊界上的節點為已知溫度，無須求解。第2 類和第3 類邊界條件都可將其等效為源項（內熱源），使用“附加源項法”［15］將邊界節點作內部節點處理。如圖5 所示的邊界節點P所在的控制體（陰影部分），有熱量輸入，并給定邊界熱流強度為qbound，則對于節點P的熱源源項

圖5 邊界上的控制體有熱量傳入

顯然SC＝qbound·Se／V，SP＝0。

程序實現時，令P節點東側傳熱系數kE＝0，P節點內熱源相關的系數SC＝qbound·Se／V，SP＝0 即可。

若邊界為第3 類邊界條件，給定邊界處的對流換熱系數α 和環境換熱溫度Tenv。則對于P節點的熱源源項

顯然

程序實現時，令P節點的aE＝0，內熱源相關系數SC＝αTenv·Se／V，SP＝-αSe／V即可。

對于不規則邊界單元，以圖6 中邊界上的控制體P為例，其東側和南側處于邊界節點上，東側為第3 類邊界條件（對流換熱系數為α，環境溫度Tenv），南側為熱流邊界條件（熱流強度為qs），其中假設Pe 段和Ps段均為絕熱邊界條件。為簡化計算，使計算域邊界上節點溫度計算仍然可使用內部節點的計算公式，增加實際上并不存在的虛擬節點S 和E，如此邊界上的節點可以等效為“內部節點”。處理與控制體S換熱的第2 類邊界條件時，令該邊界上傳熱系數kS＝0，內熱源為qsSs／VP，顯然此處的換熱面積Ss與控制體積VP

均與內部節點不同，需要做特殊處理。處理與控制體E換熱的第3 類邊界條件時，令該邊界上傳熱系數kE＝α，相鄰（虛擬）控制體節點溫度TE為Tenv。

通過調整節點的傳熱系數，增加虛擬節點和等效熱源的方法，可將邊界節點按照內部節點來處理，簡化了程序迭代計算。

2.3 仿真結果分析

在矩形計算域0 ＜x＜a，0 ＜y＜b內，求解熱物性參數為常數時的二維穩態溫度場。設邊界條件為＝0，根據數學物理方法中的分離變量法，可得其解析解

式中，θ ＝（2k+1）π／b。若a＝50、b＝30、A＝0.05、B＝10，用以上算法和解析解得到絕對誤差分布如圖7所示。可見計算結果誤差較小，表明算法能夠求解得到合理數值解。

圖7 溫度場數值解的絕對誤差分布

3 熱傳導實驗仿真過程設計

對于算法類虛擬仿真，算法的實現與驗證是關鍵。對首次使用平臺或接觸JavaScript編程的學生，提供無內熱源時的一維穩態導熱求解案例。長度為15 個單位的一維物體，兩端溫度分別為100 ℃和0 ℃。案例設計7 個關鍵步驟，引導學生了解導熱溫度場數值求解的基本流程。步驟有：JavaScript 基本語法；用戶界面及交互；物理模型與網格剖分［見圖8（a）］；溫度場初始化；邊界條件處理［見圖8（b）］；代數方程迭代求解［見圖8（c）］；殘差計算；計算結果（可參考圖2）與殘差數據可視化。圖8（c）中可見計算結果呈線性分布，說明了算法的合理性。

圖8 一維無內熱源穩態導熱數值解求解過程

經學習一維熱傳導溫度場數值解的案例后，學生已初步掌握了JavaScript編程、數值求解過程及數據可視化等基礎知識，基于此，可完成二維穩態溫度場數值解求解。其計算條件同2.3，運行界面如圖9 所示，左側為溫度場數值求解程序源代碼，右側為溫度場Contour圖和殘差趨勢圖，其驗證已在圖7 中完成。核心的穩態溫度場計算迭代程序實現僅12 行。殘差在100 次迭代時迅速收斂，后期收斂速度較慢。演示頁面同樣支持手機客戶端使用，使實驗環境可適用范圍更為廣泛。

圖9 二維無內熱源穩態導熱的數值求解

4 結語

（1）仿真平臺不僅適用于熱傳導虛擬實驗，也適用于常規的算法演示。二維熱傳導溫度場數值解虛擬仿真結果表明JavaScript可用于算法演示。

（2）基于Web 的短流程仿真平臺免去了軟件安裝與開發環境配置環節，能實現腳本程序運行和計算結果可視化，為算法演示類實驗虛擬仿真提供了便利。

（3）算法演示類實驗虛擬仿真，程序算法的編寫和調試至關重要，是理解實驗原理和目的關鍵。

（4）平臺基于開源項目構建，支持JavaScript腳本運行、用戶界面展示及數據圖表的顯示，在較簡單應用場景中可實現對Matlab一類軟件的替代，而不依賴商業軟件。