基于LBM與Matlab GUI的流體力學(xué)可視化教學(xué)與虛擬實(shí)驗(yàn)

劉向東, 戴含暉, 陳永平,,3

(1.揚(yáng)州大學(xué) 電氣與能源動(dòng)力工程學(xué)院，江蘇 揚(yáng)州 225127；2.東南大學(xué) 能源與環(huán)境學(xué)院，南京 210096；3.蘇州科技大學(xué) 環(huán)境科學(xué)與工程學(xué)院，江蘇 蘇州 215009)

0 引 言

流體力學(xué)是力學(xué)的一個(gè)重要分支，主要研究流體內(nèi)的力學(xué)規(guī)律，是水利工程、化學(xué)工程、動(dòng)力工程等諸多學(xué)科的理論基礎(chǔ)，在航空航天、船舶制造、能源電力等領(lǐng)域有著廣泛應(yīng)用[1-2]。流體力學(xué)有著理論性強(qiáng)、抽象難懂的特點(diǎn)。特別是，在多數(shù)高校的流體力學(xué)教學(xué)中，各類理論性極強(qiáng)的概念、案例往往只能通過陳述/推導(dǎo)公式的方式傳授給學(xué)生，學(xué)生在學(xué)習(xí)過程中難以獲得直觀體驗(yàn)，這大大削弱了他們的接受度和學(xué)習(xí)興趣，進(jìn)而明顯影響到了教學(xué)效果[3-5]。

隨著計(jì)算機(jī)技術(shù)和計(jì)算流體力學(xué)理論的發(fā)展，對(duì)流體力學(xué)教學(xué)及實(shí)驗(yàn)案例進(jìn)行數(shù)值虛擬仿真能夠直觀清晰地展現(xiàn)流體力學(xué)相關(guān)規(guī)律，從而為改觀傳統(tǒng)教學(xué)方法形式單一、教學(xué)成效有限的現(xiàn)狀提供了一種有效手段[6-8]。傳統(tǒng)計(jì)算流體力學(xué)方法通過數(shù)值化求解以非線性偏微分納維爾-斯托克斯方程組(Navier-Stokes equations)為核心的宏觀連續(xù)性流體力學(xué)控制方程來實(shí)現(xiàn)對(duì)流體力學(xué)問題的虛擬仿真，其數(shù)值化編程和計(jì)算求解過程都相當(dāng)復(fù)雜，并且求解計(jì)算量大、計(jì)算效率一般，從而限制了其在流體力學(xué)教學(xué)與實(shí)驗(yàn)領(lǐng)域的普及。為解決這一問題，本文將使用有別于傳統(tǒng)計(jì)算流體力學(xué)方法的格子玻爾茲曼方法(Lattice Boltzmann method，LBM)[9-10]來進(jìn)行典型流體力學(xué)問題的數(shù)值模擬仿真。與傳統(tǒng)計(jì)算流體力學(xué)模擬方法不同，LBM基于分子動(dòng)理論，通過對(duì)微觀粒子行為的統(tǒng)計(jì)平均來反應(yīng)流體宏觀動(dòng)力學(xué)特征，是一種典型的介觀模擬方法。相較于傳統(tǒng)計(jì)算流體力學(xué)方法，LBM具備先天的并行特性，且物理背景清晰、邊界條件處理簡(jiǎn)單、程序易于實(shí)現(xiàn)、計(jì)算效率高[11]，已經(jīng)在許多傳統(tǒng)計(jì)算流體力學(xué)模擬方法難以勝任的領(lǐng)域都獲得了成功應(yīng)用[12-13]，近年來逐漸受到了國(guó)內(nèi)外研究者的關(guān)注。

值得注意的是，Matlab軟件中的圖形用戶界面(Graphic User Interface，GUI)工具可通過添加控件設(shè)計(jì)人機(jī)交互界面，便于搭建可視化仿真實(shí)驗(yàn)平臺(tái)，已應(yīng)用于數(shù)學(xué)、物理、力學(xué)等課程教學(xué)中[14-16]。鑒于LBM方法在流體力學(xué)高效模擬計(jì)算和Matlab GUI工具在用戶操控與可視化方面所具有的優(yōu)勢(shì)，結(jié)合LBM方法和Matlab GUI工具設(shè)計(jì)搭建流體力學(xué)虛擬實(shí)驗(yàn)仿真平臺(tái)，用于高校流體力學(xué)教學(xué)實(shí)驗(yàn)中，為進(jìn)一步豐富流體力學(xué)課堂與實(shí)驗(yàn)教學(xué)手段、提高課堂與實(shí)驗(yàn)的教學(xué)效果提供一種嶄新工具。

1 LBM簡(jiǎn)介

LBM基于分子動(dòng)理論，將流體抽象為大量只有質(zhì)量沒有體積的微觀粒子，這些粒子可以向空間的若干方向任意移動(dòng)。LBM將粒子的演化過程分為碰撞和遷移兩個(gè)階段，因此格子玻爾茲曼方程可以分解為兩個(gè)部分。

(1)

在數(shù)值模擬過程中根據(jù)流場(chǎng)空間維數(shù)的不同，有不同的離散速度集和與之對(duì)應(yīng)的平衡態(tài)分布函數(shù)。針對(duì)二維的計(jì)算區(qū)域，選擇D2Q9模型(2維空間，9個(gè)離散速度，見圖1)，其流場(chǎng)被一系列正方形劃分。其速度、局部平衡分布函數(shù)、權(quán)函數(shù)及聲速分別如下：

(2)

圖1 D2Q9模型示意圖

在獲得速度分布函數(shù)后，流體的宏觀密度與速度通過以下公式計(jì)算得到：

(5)

計(jì)算過程中固體界面處的邊界節(jié)點(diǎn)使用非平衡外推的方法，該節(jié)點(diǎn)上分布函數(shù)計(jì)算方法如下

(6)

式中,x+為壁面法向的相鄰節(jié)點(diǎn)的位置。

LBM的演化過程可分解為碰撞和遷移兩個(gè)過程，其計(jì)算流程圖如圖2所示，具體的計(jì)算步驟如下：

步驟1給定初始密度ρ和速度場(chǎng)u；

步驟2計(jì)算平衡態(tài)分布函數(shù)feq；

步驟3在時(shí)刻t執(zhí)行碰撞過程，發(fā)生改變的節(jié)點(diǎn)x上的粒子分布函數(shù)為原有粒子分布和碰撞項(xiàng)之和；

步驟4執(zhí)行遷移，節(jié)點(diǎn)x上的粒子以速度ei運(yùn)動(dòng)到與節(jié)點(diǎn)x相鄰的x+ei，粒子分布函數(shù)也隨之遷移；

步驟5邊界處理；

步驟6計(jì)算時(shí)刻t+1的密度ρ與速度u；

步驟7重復(fù)步驟3～6直到滿足終止條件。

圖2 LBM方法計(jì)算流程圖

2 功能模塊

2.1 整體設(shè)計(jì)思路

采用Matlab GUI工具對(duì)LBM模擬計(jì)算方法進(jìn)行封裝，便可建立起流體力學(xué)虛擬實(shí)驗(yàn)仿真平臺(tái)，其用戶界面如圖3所示。由圖可知，本文所建立的仿真平臺(tái)主要包含頂蓋驅(qū)動(dòng)流、泊肅葉流和卡門渦街模塊。3個(gè)模塊對(duì)應(yīng)著流體力學(xué)中的幾個(gè)經(jīng)典問題。該仿真平臺(tái)的邏輯結(jié)構(gòu)如圖4所示。

圖3 流體力學(xué)虛擬實(shí)驗(yàn)仿真平臺(tái)用戶界面

圖4 流體力學(xué)虛擬實(shí)驗(yàn)仿真平臺(tái)模塊結(jié)構(gòu)圖

2.2 頂蓋驅(qū)動(dòng)流模塊

頂蓋驅(qū)動(dòng)流是通過腔體頂蓋水平滑移運(yùn)動(dòng)驅(qū)動(dòng)腔體內(nèi)流體產(chǎn)生流動(dòng)的典型流體力學(xué)問題。不同工況條件下的頂蓋驅(qū)動(dòng)流流場(chǎng)內(nèi)蘊(yùn)含著復(fù)雜流體力學(xué)現(xiàn)象(諸如不穩(wěn)定層流、多級(jí)渦及湍流等)，從而使之成為流體力學(xué)教學(xué)中的經(jīng)典物理模型之一。為此，本文建立的流體力學(xué)虛擬實(shí)驗(yàn)仿真平臺(tái)特包含了頂蓋驅(qū)動(dòng)流模塊，其用戶界面如圖5(a)所示。該模塊為用戶提供了自定義計(jì)算區(qū)域尺寸與流動(dòng)雷諾數(shù)的功能。用戶可通過文本框分別鍵入計(jì)算區(qū)域水平方向網(wǎng)格數(shù)、豎直方向網(wǎng)格數(shù)、流動(dòng)雷諾數(shù)，然后點(diǎn)擊“計(jì)算”按鈕完成模擬操作。頂蓋驅(qū)動(dòng)流的不同流動(dòng)狀態(tài)主要受流動(dòng)雷諾數(shù)影響。在輸入計(jì)算區(qū)域幾何尺寸后，在顯示計(jì)算結(jié)果的區(qū)域內(nèi)會(huì)給出網(wǎng)格劃分示意圖(圖5(b))。圖6顯示了不同雷諾數(shù)條件下，同一個(gè)方腔內(nèi)的流體流動(dòng)的LBM模擬計(jì)算結(jié)果。

由圖6可見，整個(gè)流場(chǎng)中的流線由藍(lán)色線條表示，計(jì)算區(qū)域中的速度矢量場(chǎng)用紅色箭頭表示。值得一提的是，為了讓圖中表示速度的矢量箭頭顯示得更加清晰，在繪制過程中對(duì)整個(gè)速度場(chǎng)做了標(biāo)準(zhǔn)化處理，即所有速度矢量的模都設(shè)為相同。這樣，圖中紅色箭頭便可直觀展示流場(chǎng)內(nèi)各處的速度方向信息。從圖中可以發(fā)現(xiàn)，隨著雷諾數(shù)Re的增大，在中央1級(jí)渦的右下部、左下部和左上部會(huì)依次逐漸出現(xiàn)2級(jí)渦。Re=200時(shí)，方腔右下部的2級(jí)渦已有形成的趨勢(shì)，到Re=1 000時(shí)右下部的2級(jí)渦已完全成型。當(dāng)Re提升到4 000時(shí)方腔左下部會(huì)出現(xiàn)2級(jí)渦。當(dāng)Re=700時(shí)，方腔左上部的2級(jí)渦也已初步形成。基于該模塊中流線與流體速度矢量的可視化展示，學(xué)生可以在圖中清晰的獲得方腔內(nèi)流體流動(dòng)的詳細(xì)流場(chǎng)分布，并能夠直觀地感受到流動(dòng)狀態(tài)在很大程度上與流體流動(dòng)的雷諾數(shù)直接相關(guān)。

2.3 泊肅葉流模塊

深入了解泊肅葉流的流動(dòng)規(guī)律是研究更為復(fù)雜流動(dòng)問題的基礎(chǔ)。本文建立的流體力學(xué)虛擬實(shí)驗(yàn)仿真平臺(tái)也納入了泊肅葉流模塊，其用戶界面如圖7所示，計(jì)算區(qū)域的網(wǎng)格尺寸與入口流速可通過對(duì)應(yīng)的文本框鍵入，再點(diǎn)擊“計(jì)算”按鈕便可獲得模擬結(jié)果。圖中顯示的是計(jì)算區(qū)域內(nèi)水平速度的分布云圖。從圖右側(cè)的顏色圖例可以看出，顏色由淺入深表示的是速度由大到小的變化過程。水平方向的速度分布總體呈現(xiàn)從流道中央向兩側(cè)遞減的規(guī)律。圖中還顯示了中央某界面上的速度分布圖，用紅色箭頭與輪廓線標(biāo)出，這使得計(jì)算結(jié)果更加直觀清晰。值得注意的是，在計(jì)算區(qū)域左側(cè)流體流入部分存在一定的過渡段。這是由于在給定入口流速時(shí)是平均分布的，由速度平均分布演變?yōu)槌浞职l(fā)展段內(nèi)拋物線狀的速度分布需要一定的過渡段。同時(shí)，考慮到不同入口流速下泊肅葉流速度云圖分布相似，該模塊還設(shè)置了計(jì)算模擬數(shù)據(jù)導(dǎo)出功能，為用戶獲得不同工況下通道內(nèi)充分發(fā)展段主流速度分布的具體計(jì)算模擬結(jié)果提供了有效接口。利用該功能，圖8給出了不同入口流速條件下充分發(fā)展段的水平方向速度分布的模擬結(jié)果。如圖所示，學(xué)生可以觀察到兩平板間的流動(dòng)大致呈現(xiàn)二次方函數(shù)的拋物線分布，垂直于流體流動(dòng)方向的流動(dòng)速度由流道中央向通道兩側(cè)壁面遞減，直至壁面上速度降為0。該模塊可以為學(xué)生直觀形象地認(rèn)識(shí)流體力學(xué)中典型的泊肅葉流主流速度分布特征提供一種有效手段。

(a) 用戶界面圖

(b) 計(jì)算網(wǎng)格圖

圖5 頂蓋驅(qū)動(dòng)流模塊用戶界面及計(jì)算網(wǎng)格圖

(a) Re=200

(b)Re=1 000

(c)Re=4 000

(d)Re=7 000

圖6 不同雷諾數(shù)條件下頂蓋驅(qū)動(dòng)流的計(jì)算結(jié)果

圖7 泊肅葉流模塊用戶界面

圖8 不同入口流速條件下泊肅葉流主流速度分布

2.4 卡門渦街模塊

卡門渦街是當(dāng)定常來流繞過流場(chǎng)中的障礙物時(shí)，障礙物兩側(cè)周期性脫落形成旋轉(zhuǎn)方向相反、排列規(guī)則的雙列線渦的流體力學(xué)現(xiàn)象。在不同的來流速度與障礙物尺寸條件下，障礙物之后的卡門渦街特征各異。本文流體力學(xué)虛擬實(shí)驗(yàn)仿真平臺(tái)所包含的卡門渦街模塊用戶界面如圖9所示，演示的是二維圓柱繞流問題。該模塊允許用戶調(diào)整入口流速u、障礙物直徑與流道寬度比pr，可以實(shí)現(xiàn)對(duì)卡門渦街特征的模擬與直觀展示。從模擬計(jì)算結(jié)果可見，該模塊可以直觀地展示出不同工況下計(jì)算區(qū)域內(nèi)流體速度大小連續(xù)分布的流場(chǎng)云圖，計(jì)算結(jié)果界面還給出了每個(gè)模擬計(jì)算工況所對(duì)應(yīng)的不同尺寸圓柱障礙物的輪廓示意。圖9(a)～(c)顯示的是圓柱直徑一定、入口流速不同的條件下，各個(gè)流場(chǎng)計(jì)算結(jié)果的對(duì)比。圖9(a)中的入口流速為30 m/s，該條件下整個(gè)計(jì)算區(qū)域內(nèi)的流動(dòng)較為平穩(wěn)，在圓柱的后方基本沒有出現(xiàn)線渦脫落的現(xiàn)象。當(dāng)入口流速增大到250 m/s(見圖9(b))，圓柱之后出現(xiàn)了明顯的擾動(dòng)，但隨著流動(dòng)向下游發(fā)展，這一擾動(dòng)也逐漸趨于緩和。而繼續(xù)增加入口流速到500 m/s(見圖9(c))，可以發(fā)現(xiàn)在圓柱之后的整個(gè)計(jì)算區(qū)域內(nèi)都出現(xiàn)了明顯的線渦脫落現(xiàn)象。當(dāng)在圖9(c)的基礎(chǔ)上將障礙物半徑占流道寬度的比例從0.3減小到0.1時(shí)(見圖9(d))，圓柱后的流場(chǎng)波動(dòng)仍然顯著，但波動(dòng)幅度相較于圖9(c)有所減小而波動(dòng)頻率則有所增加，這就說明隨著障礙物半徑占流道寬度比例的減小，卡門渦街所產(chǎn)生的線渦強(qiáng)度有所減小，同時(shí)其脫落頻率則有所提升。通過以上算例，學(xué)生可以直觀地認(rèn)識(shí)到，在流體繞障礙物流動(dòng)的過程中，流體流速較高時(shí)更易在障礙物后產(chǎn)生卡門渦街，并且渦街產(chǎn)生的頻率和強(qiáng)度與障礙物與流道間的尺寸比例直接相關(guān)。

3 虛擬仿真實(shí)驗(yàn)平臺(tái)在教學(xué)活動(dòng)中的應(yīng)用方式與優(yōu)勢(shì)所在

本文建立的流體力學(xué)虛擬實(shí)驗(yàn)仿真平臺(tái)為學(xué)生提供了一個(gè)直觀、清晰認(rèn)知經(jīng)典流體力學(xué)的有效窗口，為改觀傳統(tǒng)陳述/推導(dǎo)公式型的教學(xué)方法提供了契機(jī)。任課教師在課堂教學(xué)過程中可以通過引導(dǎo)學(xué)生在計(jì)算機(jī)上進(jìn)行仿真平臺(tái)的實(shí)際操作來加深對(duì)經(jīng)典流體力學(xué)的理解認(rèn)知，并直觀體會(huì)到關(guān)鍵工況參數(shù)對(duì)流體流動(dòng)特性的影響。同時(shí)，教師還可以就此向?qū)W生引申介紹后續(xù)《計(jì)算流體力學(xué)》《數(shù)值傳熱學(xué)》等相關(guān)課程的基本思想，進(jìn)一步拓展學(xué)生的知識(shí)面，并為其未來專業(yè)課的學(xué)習(xí)打下良好基礎(chǔ)。

相對(duì)于傳統(tǒng)流體力學(xué)教學(xué)方法，使用Matlab GUI設(shè)計(jì)的流體力學(xué)仿真模擬實(shí)驗(yàn)平臺(tái)進(jìn)行輔助教學(xué)的優(yōu)勢(shì)明顯：

(1) 具有良好的可視化效果。虛擬實(shí)驗(yàn)仿真平臺(tái)最顯著的優(yōu)點(diǎn)就是具備良好的可視化效果。通過矢量場(chǎng)圖、流線、速度云圖將各個(gè)算例清晰直觀地展現(xiàn)在學(xué)生面前，有利于強(qiáng)化學(xué)生對(duì)該方面知識(shí)的理解，為后續(xù)內(nèi)容的學(xué)習(xí)打好基礎(chǔ)。

(a) u=30 m/s，pr=0.3

(b) u=250 m/s，pr=0.3

(c) u=500 m/s，pr=0.3

(d) u=500 m/s，pr=0.1

(2) 操作簡(jiǎn)單，運(yùn)行速度快。本文所建立的虛擬仿真實(shí)驗(yàn)平臺(tái)界面清晰，操作簡(jiǎn)便。同時(shí)，該平臺(tái)基于LBM方法進(jìn)行數(shù)值模擬仿真，且本文已根據(jù)Matlab GUI的運(yùn)行特點(diǎn)對(duì)程序內(nèi)的運(yùn)算邏輯做了優(yōu)化。因此，相較于基于傳統(tǒng)計(jì)算流體力學(xué)方法的數(shù)值仿真平臺(tái)，該平臺(tái)運(yùn)行速度得到了明顯提升，可在幾分鐘內(nèi)獲得合理計(jì)算結(jié)果。

(3) 其他優(yōu)勢(shì)。該平臺(tái)所基于的LBM方法，邊界條件處理簡(jiǎn)單，程序易于實(shí)現(xiàn)且開源，學(xué)有余力的學(xué)生可以通過自行編譯程序，對(duì)幾個(gè)流體力學(xué)問題進(jìn)行自主模擬。這充分體現(xiàn)了“因材施教”的教育理念，并有助于進(jìn)一步提升學(xué)生的計(jì)算機(jī)程序編譯水平。虛擬仿真實(shí)驗(yàn)平臺(tái)節(jié)約了實(shí)驗(yàn)過程中的耗材耗能，成本低廉、節(jié)能環(huán)保，從而實(shí)現(xiàn)了“綠色教學(xué)”。

4 結(jié) 語

采用LBM對(duì)流體力學(xué)教學(xué)中的經(jīng)典案例進(jìn)行了數(shù)值模擬，并基于Matlab GUI工具對(duì)其進(jìn)行封裝，開發(fā)了流體力學(xué)虛擬實(shí)驗(yàn)仿真平臺(tái)。該仿真平臺(tái)主要涵蓋了流體力學(xué)知識(shí)體系中頂蓋驅(qū)動(dòng)流、泊肅葉流和卡門渦街3個(gè)經(jīng)典案例模塊。學(xué)生可以在平臺(tái)用戶界面上輸入可變工況參數(shù)進(jìn)行計(jì)算，并可以獲得流場(chǎng)中的速度矢量場(chǎng)、流線等重要信息的直觀展示。同時(shí)，學(xué)有余力的學(xué)生還可以通過自行編譯程序，對(duì)幾個(gè)流體力學(xué)問題進(jìn)行自主模擬。

相對(duì)于流體力學(xué)傳統(tǒng)課堂與實(shí)驗(yàn)教學(xué)方法，本虛擬仿真實(shí)驗(yàn)平臺(tái)具有可視化效果好，計(jì)算速度快，界面清晰操作簡(jiǎn)單，教學(xué)成本低廉，方便交流等優(yōu)點(diǎn)，有助于激發(fā)學(xué)生的學(xué)習(xí)興趣，加深學(xué)生對(duì)抽象概念與復(fù)雜流動(dòng)狀態(tài)的理解認(rèn)識(shí)，課堂與實(shí)驗(yàn)教學(xué)效果的提高，值得在流體力學(xué)教學(xué)中進(jìn)行推廣。