Matlab求解庫(kù)埃特流在計(jì)算流體力學(xué)課程中的應(yīng)用1)

劉 益 付小莉

(同濟(jì)大學(xué)土木工程學(xué)院，上海 200092)

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

1 問(wèn)題的提出

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

1.1 解析解

由納維?斯托克斯方程進(jìn)行簡(jiǎn)化，得到不可壓恒溫庫(kù)埃特流動(dòng)的控制方程

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

圖1 庫(kù)埃特流動(dòng)

1.2 數(shù)值解

在數(shù)值求解時(shí)，傳統(tǒng)方法采用顯式和隱式方法，為了加深學(xué)生對(duì)數(shù)值方法的理解，擬增加第三種迭代算法，即壓力修正法進(jìn)行計(jì)算。為了方便后續(xù)的討論，將參數(shù)進(jìn)行無(wú)量綱化，定義為

經(jīng)過(guò)化簡(jiǎn)，得到數(shù)值求解的控制方程為

式中，Re為兩平板之間距離計(jì)算得到的雷諾數(shù)。

為了簡(jiǎn)化記號(hào)，本文之后所有討論的變量上帶有的“′”號(hào)都將省略。在計(jì)算之前，設(shè)定初始邊界條件為

（1）利用顯式方法計(jì)算時(shí)，對(duì)時(shí)間進(jìn)行一階向前差分，對(duì)坐標(biāo)進(jìn)行二階中心差分，得到數(shù)值求解控制方程為

（2）利用隱式方法計(jì)算時(shí)，采用克蘭克?尼克爾森方法求解，該格式無(wú)條件穩(wěn)定，控制方程為

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

先假定一個(gè)初始條件，在計(jì)算區(qū)內(nèi)形成二維流動(dòng)，然后在迭代過(guò)程中觀察流場(chǎng)是否收斂到精確解，具體求解流程見(jiàn)圖2。

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

2 Matlab求解控制方程

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

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

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

圖3 不同算法收斂圖

3 數(shù)值計(jì)算誤差與求解效率分析

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

3.1 推進(jìn)步數(shù)的影響

控制等分?jǐn)?shù)為10，時(shí)間間隔為0.001 s，雷諾數(shù)為63.6，而推進(jìn)步數(shù)在50～10 000中變動(dòng)，其中點(diǎn)速度結(jié)果值如表1。

表1 不同推進(jìn)步數(shù)下中點(diǎn)速度值

經(jīng)過(guò)計(jì)算可以發(fā)現(xiàn)，壓力修正法的收斂速度較快，在很少的步數(shù)內(nèi)達(dá)到求解收斂，而顯式方法與隱式方法求解速度相近，相同情況下需要迭代的步數(shù)較多。

3.2 等分?jǐn)?shù)的影響

控制推進(jìn)步數(shù)為10 000，時(shí)間間隔為0.001 s，雷諾數(shù)為63.6，而等分?jǐn)?shù)在4～16中變動(dòng)，其中點(diǎn)速度結(jié)果值如表2。

表2 不同等分?jǐn)?shù)下中點(diǎn)速度值

表2表明，隨著等分?jǐn)?shù)增大，網(wǎng)格變密，計(jì)算速度顯著變慢。而壓力修正法在不同等分?jǐn)?shù)下仍有很高的精度。尤其在等分?jǐn)?shù)較多時(shí)，10 000步內(nèi)仍然可以達(dá)到收斂，也體現(xiàn)其算法的高效。

3.3 時(shí)間間隔的影響

控制推進(jìn)步數(shù)為10 000，等分?jǐn)?shù)為10，雷諾數(shù)為63.6，而時(shí)間間隔在0.001～0.636 s中變動(dòng)，其中點(diǎn)速度結(jié)果值如表3。

表3 不同時(shí)間間隔下中點(diǎn)速度值

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

3.4 雷諾數(shù)的影響

控制等分?jǐn)?shù)為10，時(shí)間間隔為0.001 s，推進(jìn)步數(shù)為10 000，而雷諾數(shù)在10～2000中變動(dòng)，其中點(diǎn)速度結(jié)果值如表4。

表4 不同雷諾數(shù)下中點(diǎn)速度值

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

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

4 結(jié)語(yǔ)

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