Matlab求解庫埃特流在計算流體力學課程中的應用1)

劉 益 付小莉

(同濟大學土木工程學院，上海 200092)

課程內容設計應該注重開放性、改變學生的知識觀和學習觀，并構建以教師為主導、以學生為中心的教學模式[1]。計算流體力學作為機械、航天、車輛乃至土木工程等專業的必修課，兼具有流體力學和工程數學的特點。但對于計算流體力學在內的相關課程，公式多、符號抽象、計算復雜是影響學生學習的絆腳石。另外，由于對以往理論力學等課程的基礎概念已經生疏，學生在學習理論基礎時常常感到困難[2]。為了擺脫傳統教學模式，讓學生更好地參與進課程教學中，鍛煉學生的編程求解能力，同時拓寬流體力學的教學范圍，開展了不可壓縮庫埃特流的數值解實驗，并對不同方法的求解精度與求解效率等方面進行討論。庫埃特流動是一種存在理論解析解的經典流動，非常適合剛剛入門計算流體力學的同學上手實踐。

1 問題的提出

庫埃特流動是最簡單的由剪切力驅動的黏性流動，與復雜的邊界層流動具有很多相同的物理特征。其具體問題描述為：設有兩個平行的平板，距離為D，上平板以恒定的速度ue運動，下平板則保持靜止，試考慮這兩個平板之間的黏性流動[3-6]。

1.1 解析解

由納維?斯托克斯方程進行簡化，得到不可壓恒溫庫埃特流動的控制方程

求解該微分方程較為簡單，得到解析解如式（2），速度分布為圖1所示。

圖1 庫埃特流動

1.2 數值解

在數值求解時，傳統方法采用顯式和隱式方法，為了加深學生對數值方法的理解，擬增加第三種迭代算法，即壓力修正法進行計算。為了方便后續的討論，將參數進行無量綱化，定義為

經過化簡，得到數值求解的控制方程為

式中，Re為兩平板之間距離計算得到的雷諾數。

為了簡化記號，本文之后所有討論的變量上帶有的“′”號都將省略。在計算之前，設定初始邊界條件為

（1）利用顯式方法計算時，對時間進行一階向前差分，對坐標進行二階中心差分，得到數值求解控制方程為

（2）利用隱式方法計算時，采用克蘭克?尼克爾森方法求解，該格式無條件穩定，控制方程為

（3）壓力修正法的求解基本思路

先假定一個初始條件，在計算區內形成二維流動，然后在迭代過程中觀察流場是否收斂到精確解，具體求解流程見圖2。

圖2 壓力修正法求解流程圖

2 Matlab求解控制方程

第1節通過差分法得到了顯式方法與隱式方法的控制方程，并提出了基于壓力修正的第三種方法。本節將采用Matlab編程進行問題求解。

對于顯式方法與隱式方法，需要設定邊界條件與第一迭代步的初始條件，即為式（7），通過前一迭代步與本迭代步的關系進行求解，顯式方法直接將參數代入即可，但是隱式方法需要求解方程組。而壓力修正法則與上述二者不同，在給定邊界條件時，需要施加一個實際擾動以形成流場，而并非直接進行迭代計算。這里采用兩種方式，一種是在中心網格點插入豎直脈沖速度，另一種是在中心網格點插入水平脈沖速度，雖然二者的方向不同，但是在計算收斂性和精確度上是相近的。

圖3分別為顯式方法、隱式方法、壓力修正法（引入豎直脈沖）和壓力修正法（引入水平脈沖）的迭代求解收斂圖，可以看到，每個關鍵網格點的速度由一開始設定的初值逐漸向理論解析解靠近。

圖3 不同算法收斂圖

3 數值計算誤差與求解效率分析

計算流體力學乃至整個數值計算的研究領域內，最關心的就是求解精度與效率。為了探究以上問題，本節將集中討論四個自變量：推進步數、等分數、時間間隔、雷諾數對四種求解方法的誤差影響。由于所有網格點中最上方和最下方點的速度已由給定初始條件確定，擬采用中點的速度值作誤差對比，其解析解的值為0.5。

3.1 推進步數的影響

控制等分數為10，時間間隔為0.001 s，雷諾數為63.6，而推進步數在50～10 000中變動，其中點速度結果值如表1。

表1 不同推進步數下中點速度值

經過計算可以發現，壓力修正法的收斂速度較快，在很少的步數內達到求解收斂，而顯式方法與隱式方法求解速度相近，相同情況下需要迭代的步數較多。

3.2 等分數的影響

控制推進步數為10 000，時間間隔為0.001 s，雷諾數為63.6，而等分數在4～16中變動，其中點速度結果值如表2。

表2 不同等分數下中點速度值

表2表明，隨著等分數增大，網格變密，計算速度顯著變慢。而壓力修正法在不同等分數下仍有很高的精度。尤其在等分數較多時，10 000步內仍然可以達到收斂，也體現其算法的高效。

3.3 時間間隔的影響

控制推進步數為10 000，等分數為10，雷諾數為63.6，而時間間隔在0.001～0.636 s中變動，其中點速度結果值如表3。

表3 不同時間間隔下中點速度值

計算后，伴隨著時間間隔變大，算法開始不收斂，而只有隱式算法是無條件穩定的，但有可能失去時間精度，事實上，存在一個最優的時間間隔使得隱式方法效率最高。而另外三個算法的要求則較高，顯式方法一定需要保證穩定參數E≤0.5才能保證收斂，而壓力修正法需要在一定限度的時間間隔內才能得到正確答案。

3.4 雷諾數的影響

控制等分數為10，時間間隔為0.001 s，推進步數為10 000，而雷諾數在10～2000中變動，其中點速度結果值如表4。

表4 不同雷諾數下中點速度值

從表4可以看出，雷諾數對壓力修正法影響較大，由于壓力修正法是經過連續性方程和動量方程耦合聯立推導得到的，在數值解的情況下更具有實際意義，即需保證存在真實的流場滿足參數設定，否則有可能在極端條件下不收斂。另一方面，雷諾數變大也會影響顯式與隱式方法的收斂速度，但是對于收斂性的影響不大，隱式方法在不同雷諾數情況下均能收斂。

從以上四個變量的計算中得到結論，顯式算法與隱式算法對時間間隔的大小與等分數比較敏感，即編程者需要確定一個合適的網格劃分尺寸與計算精度；而壓力修正法對雷諾數等較為敏感，在雷諾數較高或不符合實際的情況下，有可能導致數值振蕩，影響求解的精度與效率。

4 結語

大學教師往往更重視對學生專門知識與多項技能的培養[7]，而大學課堂上的教學除了最基本的教授知識外，更重要的是教給學生獲取知識的方法和激發學生的學習興趣，由此才能讓學生主動去學習對自己有幫助的課程與知識[8]。本文通過一個非常經典的流動模型，來啟發學生通過多種方式進行數值求解，并探究基本參數與誤差之間的影響關系。經過練習，學生能很好地入門計算流體力學領域，了解求解的一般方式與步驟，鍛煉編程能力和作圖能力，為后續的學習科研打下基礎。