基于MATLAB GUI 的工程流體力學(xué)可視化教學(xué)與上機(jī)實驗

鄧梓龍，徐 澤，張程賓

（東南大學(xué) 能源與環(huán)境學(xué)院，江蘇 南京 210096）

“工程流體力學(xué)”是能源動力工程、機(jī)械工程、土木工程等專業(yè)的重要專業(yè)基礎(chǔ)課程[1-3]。在工程流體力學(xué)教學(xué)中遇到復(fù)雜的公式，單純進(jìn)行理論推導(dǎo)的效果往往不佳[4-5]。本文通過 MATLAB(graphical user interface, GUI)設(shè)計出一款可視化工程流體力學(xué)實驗教學(xué)軟件，對工程流體力學(xué)中頂蓋驅(qū)動流、泊肅葉流、科特流這3 種典型算例進(jìn)行了流場的數(shù)值求解，并上機(jī)進(jìn)行實驗演示。該實驗教學(xué)軟件的應(yīng)用，豐富了課堂教學(xué)內(nèi)容，可以幫助學(xué)生更好地理解流體動力學(xué)原理，取得了較好的課堂教學(xué)效果。

1 SIMPLE 算法

1.1 SIMPLE 算法計算流程

計算流體力學(xué)仿真需要采用數(shù)值方法對納維-斯托克斯方程組進(jìn)行求解，從而預(yù)測流場中的流動。流體力學(xué)中不可壓縮流動的數(shù)值解法分為聯(lián)立求解法和分離式求解法[6]。聯(lián)立求解法需聯(lián)立求解所有或部分變量的代數(shù)方程，計算效率低，應(yīng)用較少。分離式求解法是按順序求解各變量，又可分為以速度、壓力作為變量的原始變量法和以渦量、流函數(shù)作為變量的非原始變量法。原始變量法中的壓力修正法是目前被使用最廣泛的一種計算方法[7]。求解壓力耦合方程組的半隱式方法(semi-implicit method for pressure-linked equations，SIMPLE)是一種常用于求解不可壓縮流場的壓力修正法[6,8]。SIMPLE 算法由Patankar 和Spalding在 1972 年提出，后期經(jīng)過多次改進(jìn)，衍生出SIMPLER、SIMPLEC、SIMPLEX、SIMPLE Date、PISO等多種方案。以SIMPLE 算法為基礎(chǔ)，可以模擬求解工程流體力學(xué)中的典型流動問題，該算法成為上機(jī)實驗演示的理論基礎(chǔ)。

商業(yè)數(shù)學(xué)軟件MATLAB 具有計算效率高、圖形功能豐富、交互性好、擴(kuò)充能力好等特點[9-10]，其GUI由圖形控件、窗口、事件響應(yīng)等元素構(gòu)成，用戶在圖形界面輸入信息，程序則會通過回調(diào)函數(shù)響應(yīng)并輸出結(jié)果。MATLAB GUI 具有強大的可視化功能，控件豐富，用戶可以簡單、便捷地設(shè)計出友好、美觀的人機(jī)交互界面[11-13]。

SIMPLE 算法解決了不可壓縮流場數(shù)值計算中壓力與速度耦合的問題，并且數(shù)值特性好、表達(dá)概念清晰、易于實施，因而被廣泛應(yīng)用于計算流體力學(xué)、計算傳熱學(xué)等領(lǐng)域[14]。SIMPLE 算法的主要思想是“假設(shè)流場-改進(jìn)流場”[15]。

SIMPLE 算法通常包含4 個步驟：

步驟1假設(shè)壓力場p*，速度場u*，v*；

步驟2利用p*求解u*、v*的動量離散方程，SIMPLE 算法近似省略了鄰點速度的影響，即∑anbu*nb=0，∑anbv*nb=0，求解得到u*，v* ：

步驟3上一步計算出的u*，v*未能滿足質(zhì)量守恒，利用質(zhì)量守恒方程改進(jìn)壓力場，求解得到壓力修正量p′：

式中，ρ 為流體密度，ρ0表示初始密度，b 值代表一個控制體不滿足連續(xù)性方程的剩余質(zhì)量大小。當(dāng)各控制體的b 值的代數(shù)和小于容許誤差，可判定速度場收斂。

步驟4修正壓力p、速度u 與v。若速度場收斂，計算結(jié)束，若未收斂，使用本層修正后的p、u、v 作為下一輪迭代初值，進(jìn)行下一輪迭代，直到收斂為止[16]。

1.2 交錯網(wǎng)格

若將速度與壓力存儲在同一個控制容積節(jié)點，則在P 點沿u 方向的離散方程中出現(xiàn)pE與pW的壓差，而沒有pp，結(jié)果可能會帶來波動的流場。SIMPLE 算法常采用交錯網(wǎng)格，將標(biāo)量存儲在主控制容積，u 與v的控制容積分別與主控制容積在x 方向和y 方向有半個網(wǎng)格步長的錯位[6]，如圖1 所示。

圖1 交錯網(wǎng)格

2 軟件介紹

2.1 軟件界面

工程流體力學(xué)可視化教學(xué)軟件包含選擇案例、選擇工況、物理參數(shù)、模型尺寸、計算繪圖5 個部分。在“計算案例”選單中，可選擇頂蓋驅(qū)動流、泊肅葉流、科特流這3 個基本算例，然后可在“選擇工況”選單中設(shè)置流動工況。軟件設(shè)置了Re=1 000、Re=2 000、Re=5 000 這3 個默認(rèn)工況選項和1 個自定義選項。

設(shè)置流動工況后，軟件會給出默認(rèn)物性參數(shù)和模型尺寸。用戶也可點自定義選項，自行設(shè)置密度、黏度、流速等物性參數(shù)及模型尺寸參數(shù)。點擊計算繪圖區(qū)域中“計算雷諾數(shù)”按鈕，即可計算當(dāng)前工況下的雷諾數(shù)Re。參數(shù)設(shè)置完成后，點擊“計算案例”按鈕便可開始計算，計算結(jié)束后在繪圖區(qū)顯示計算結(jié)果云圖。學(xué)生可自由選擇不同算例、設(shè)置各項參數(shù)，通過軟件計算得到不同案例在不同工況下的速度分布云圖，并分析計算結(jié)果。

2.2 算例1：頂蓋驅(qū)動流

頂蓋驅(qū)動流是通過方腔頂蓋運動驅(qū)動方腔內(nèi)流體運動的經(jīng)典流體力學(xué)問題，常用于驗證理論計算的正確性。頂蓋驅(qū)動流內(nèi)部流場的形態(tài)會隨著Re 的變化而形成顯著的差異。直觀清晰地了解頂蓋驅(qū)動流中的流場變化規(guī)律有助于學(xué)生更深入認(rèn)識流體力學(xué)中的無量綱數(shù)與流動狀態(tài)的關(guān)系。因此，本工程流體力學(xué)教學(xué)實驗平臺將頂蓋驅(qū)動流算例納入其中。

頂蓋驅(qū)動流方腔大小為1 m×1 m，流體密度ρ＝1 kg/m3，黏度μ＝0.001 Pa?s，方腔左右及下壁面因固壁不滑移，u 與v 皆為0。當(dāng)頂蓋流體流速u 分別為1 m/s、2 m/s、5 m/s 時，Re 分別為1 000、2 000、5 000，計算選取X 方向長度為特征長度。

在頂蓋驅(qū)動流計算案例中，分別選擇Re=1 000、Re=2 000、Re=5 000 這3 種默認(rèn)工況，物性參數(shù)及模型尺寸亦自動給出默認(rèn)值，點擊“計算案例”按鈕進(jìn)行計算，計算結(jié)果如圖2 示。

當(dāng)Re＝1 000 時，方腔內(nèi)部存在3 個旋渦，最大的一級渦位于中間偏右上角位置，另外2 個二級渦分別位于左下角及右下角，且右下角的渦大于左下角的渦，如圖2(a)所示。

當(dāng)Re＝2 000 時，一級渦位置向左下偏移，左下角的二級渦變大，位置向右上角移動，右下角的二級渦略微變小且位置向下偏移，如圖2(b)所示。

當(dāng)Re＝5 000 時，左下角及右下角的二級渦均變大，在左上角又出現(xiàn)一個新的渦，如圖2(c)所示。

學(xué)生通過該算例的結(jié)果分析可直觀了解Re 對頂蓋驅(qū)動流流場的影響。

圖2 頂蓋驅(qū)動流模型示意圖及不同Re 下速度分布

2.3 算例2：泊肅葉流

泊肅葉流是不可壓黏性流體在管內(nèi)或兩平板間流動，是流體力學(xué)中的一種基本的流動現(xiàn)象。研究泊肅葉流的規(guī)律，對于醫(yī)藥學(xué)、工業(yè)物料運輸、生物科學(xué)等領(lǐng)域具有重要意義。泊肅葉流算例可以幫助學(xué)生直觀地了解泊肅葉流在充分發(fā)展段的速度分布以及雷諾數(shù)對流動狀況的影響。

泊肅葉流模型長10 m，寬1 m，上下壁面u 與v皆為0，假設(shè)右側(cè)出口充分發(fā)展。流體密度ρ＝1 kg/m3，黏度μ＝0.001 Pa?s。設(shè)置左側(cè)入口流體流速u 分別為1 m/s、2 m/s、5 m/s 時，Re 分別為1 000、2 000、5 000，計算選取入口寬度為特征長度。

在泊肅葉流算例中，分別選擇Re＝1 000、Re＝2 000、Re＝5 000 這3 種默認(rèn)工況，物性參數(shù)及模型尺寸亦自動給出默認(rèn)值。點擊“計算案例”按鈕進(jìn)行計算，計算結(jié)果如圖3 所示。

在圖3 中，黃色區(qū)域為高流速區(qū)域，藍(lán)色區(qū)域為低流速區(qū)域，近壁面處流速為0，泊肅葉流流速呈現(xiàn)中間高、兩邊低的拋物線分布。在充分發(fā)展段截面上，流體最高流速高于入口速度。隨著Re 增大，充分發(fā)展的截面上最高流速也隨之增大。學(xué)生可通過該案例直觀地理解泊肅葉流的截面速度分布以及影響流動狀況的因素。

2.4 算例3：科特流

科特流是黏性流體在兩個相對運動平板間的流動，是流體力學(xué)中基本且重要的流動現(xiàn)象，在流體黏度測量、汽輪機(jī)運行、大氣科學(xué)、洋流研究等場景均有應(yīng)用?？铺亓魉憷箤W(xué)生了解科特流中流體剪切力對流動狀況的影響以及不同Re 下科特流流場差異。

科特流模型長為10 m，寬為1 m，下層流體流速為u2=－2 m/s，左右出口皆為充分發(fā)展。流體密度ρ＝1 kg/m3，黏度μ＝0.001 Pa?s。設(shè)置上層流體流速u1分別為1 m/s、2 m/s、5 m/s 時，Re 分別為1 000、2 000、5 000，計算選取X 方向長度為特征長度。

在科特流算例中，分別選擇Re=1 000、Re=2 000、Re=5 000 這3 種默認(rèn)工況，物性參數(shù)及模型尺寸亦自動給出默認(rèn)值。點擊“計算案例”按鈕進(jìn)行計算，計算結(jié)果如圖4 所示。

科特流在Y 方向存在速度梯度，流速呈現(xiàn)中間低、兩邊高的分布，且在中間區(qū)域存在旋渦。當(dāng)Re＝1 000時，旋渦位置接近中間；當(dāng)Re＝2 000 時，旋渦位置向下偏移，當(dāng)Re＝5 000 時，上層流體流速大于下層流體流速，旋渦位置又向上偏移。學(xué)生通過該案例的計算結(jié)果，能夠更直觀地體會剪切流的流動狀況。

圖3 泊肅葉流模型示意圖及及不同Re 下速度分布

3 工程流體力學(xué)可視化教學(xué)軟件優(yōu)勢

基于MATLAB GUI 的工程流體力學(xué)可視化教學(xué)軟件的交互界面友好、清晰，操作步驟簡單易學(xué)，能夠強化學(xué)生對流體流動過程的理解，促進(jìn)學(xué)生對關(guān)鍵知識點的掌握，從而達(dá)到提升教學(xué)質(zhì)量的效果。此外，教師還可以以此為基礎(chǔ)，向?qū)W生介紹計算流體力學(xué)與數(shù)值傳熱學(xué)的相關(guān)知識，為學(xué)生后續(xù)課程的學(xué)習(xí)打下基礎(chǔ)。在教學(xué)安排方面，教師不僅可以在課堂上親自示范帶領(lǐng)學(xué)生開展流體力學(xué)相關(guān)內(nèi)容的學(xué)習(xí)，還可以在講解完相關(guān)的理論基礎(chǔ)后將可視化軟件的操作作為課后作業(yè)，留給學(xué)生自行探索與完成。這樣能夠充分調(diào)動學(xué)生的主觀能動性，培養(yǎng)學(xué)生的自主學(xué)習(xí)能力，并且還能使學(xué)生在課后對所學(xué)知識進(jìn)行鞏固，加深學(xué)生對所學(xué)知識的理解。

相比傳統(tǒng)工程流體力學(xué)課堂教學(xué)中采用的講授模式，使用工程流體力學(xué)可視化教學(xué)軟件輔助教學(xué)具有以下優(yōu)勢：

圖4 科特流模型示意圖及不同Re 下速度分布

（1）可視化效果好。工程流體力學(xué)可視化教學(xué)軟件圖形界面易于操控，軟件可對頂蓋驅(qū)動流、泊肅葉流、科特流進(jìn)行多種工況下的計算，計算所得結(jié)果能夠幫助學(xué)生直觀地了解各流動特性，增強學(xué)生對理論知識的理解。

（2）計算效率高。工程流體力學(xué)可視化教學(xué)軟件中算例的流場采用SIMPLE 算法求解，該算法概念清晰、計算效率高、數(shù)值特性好，學(xué)生可結(jié)合軟件源碼中SIMPLE 算法實現(xiàn)過程加深對該算法的理解。

（3）軟件實行開源，可供學(xué)生學(xué)習(xí)使用，易于推廣到課堂輔助教學(xué)，降低教學(xué)成本。學(xué)生可以通過學(xué)習(xí)源碼加深對流動的數(shù)值解法及MTALAB GUI 的認(rèn)識與理解，也可以在軟件現(xiàn)有內(nèi)容的基礎(chǔ)上添加新算例，提高對流動的數(shù)值解法編程計算的能力。

4 結(jié)論

相比較傳統(tǒng)工程流體力學(xué)實驗教學(xué)方式，基于MATLAB GUI 的工程流體力學(xué)可視化軟件教學(xué)具有可視化效果好、成本低、易于推廣等特點。算例采用SIMPLE 算法求解，學(xué)生可自由設(shè)定物性參數(shù)、模型尺寸等參數(shù)，計算得到速度場云圖。學(xué)生可將計算結(jié)果與理論相結(jié)合，強化對課堂理論知識的理解，從而提升課堂教學(xué)效果。