水利機(jī)械計(jì)算水動(dòng)力學(xué)------提交報(bào)告

1、一、求解N-S方程的數(shù)值方法求解 N-S 方程的數(shù)值方法可分為有網(wǎng)格和無(wú)網(wǎng)格兩類。有網(wǎng)格方法的求解思路是將連續(xù)空間劃分為離散網(wǎng)格 (或單元)，并在網(wǎng)格點(diǎn) (或單元) 上離散控制方程，從而尋求控制方程的數(shù)值解；無(wú)網(wǎng)格方法則將著眼點(diǎn)回歸到流體運(yùn)動(dòng)的基本單元質(zhì)點(diǎn)上，通過追蹤質(zhì)點(diǎn)運(yùn)動(dòng)性質(zhì)的變化過程來(lái)獲得流體解答。網(wǎng)格數(shù)值模擬方法發(fā)展較早，求解方法較成熟，應(yīng)用也較廣；而無(wú)網(wǎng)格方法則剛剛起步，具有很大的發(fā)展空間。下面分別就兩類中較有代表性的方法予以介紹。1.1 網(wǎng)格數(shù)值模擬方法數(shù)值計(jì)算網(wǎng)格可分為不隨流體運(yùn)動(dòng)而改變形狀的歐拉網(wǎng)格和隨流體運(yùn)動(dòng)改變形狀的拉格朗日網(wǎng)格兩類。歐拉網(wǎng)格劃分較簡(jiǎn)單，但在用于具有自由表面

2、流體運(yùn)動(dòng)數(shù)值計(jì)算中需要對(duì)自由表面位置進(jìn)行適時(shí)追蹤；拉格朗日網(wǎng)格雖然無(wú)需進(jìn)行自由表面追蹤，但需要在計(jì)算過程中隨時(shí)調(diào)整計(jì)算網(wǎng)格，以適應(yīng)流體自由表面形狀，額外增加了計(jì)算量，并且在流體變形劇烈處還會(huì)出現(xiàn)網(wǎng)格失效現(xiàn)象。有限差分法 (Finite Difference Method, FDM)和邊界元法 (Boundary Element Method, BEM) 都是歐拉網(wǎng)格方法中常用的方法。在拉格朗日網(wǎng)格法中，有限元法和邊界元法可直接應(yīng)用，而有限差分法一般需對(duì)坐標(biāo)做適體變換。用歐拉網(wǎng)格法求解含有自由表面的水波動(dòng)力問題難點(diǎn)在于追蹤自由表面。目前比較有代表性的自由表面追蹤方法有 PIC (Particle

3、-In-Cell) 法、MAC (Marker-And-Cell) 法和VOF (VolumeOfFluid) 法等。以下分別予以介紹。1.1.1 質(zhì)點(diǎn)網(wǎng)格法 (Particle-in-Cell，簡(jiǎn)稱 PIC 法) PIC法是 Harlow提出的一種用歐拉矩形網(wǎng)格計(jì)算多種介質(zhì)流體運(yùn)動(dòng)的方法。該方法把流體既視為連續(xù)介質(zhì)，又視為帶有一定質(zhì)量的質(zhì)點(diǎn)，然后研究質(zhì)點(diǎn)在經(jīng)過固定歐拉網(wǎng)格上的運(yùn)動(dòng)性質(zhì)。PIC 法具有計(jì)算多相流和處理三維自由表面的能力，曾經(jīng)成功運(yùn)用于模擬二維流體發(fā)生劇烈變形的情況。其缺點(diǎn)是占用較多計(jì)算機(jī)內(nèi)存；模式只能給出自由表面單元位置，而不能給出自由表面精確位置；另外在流體變形劇烈處易出現(xiàn)較

4、大插值誤差。1.1.2 MAC 法 (Marker-And-Cell) 針對(duì) PIC 法的上述缺點(diǎn)，20 世紀(jì) 60 年代中期，Harlow, Welch提出了一種改進(jìn)方法，稱為 MAC 法。最初的 MAC 法在所有單元內(nèi)部都布滿無(wú)質(zhì)量的標(biāo)記點(diǎn)，通過跟蹤這些點(diǎn)即可判斷自由表面的位置。MAC 法在求解時(shí)以速度分量和壓力為因變量，在固定網(wǎng)格上采用有限差分法離散控制方程。在 MAC 法之后，人們又先后提出了 SMAC 法、ABAMC 法、SUMAC 法和 TUMMAC 法11系列。新版本模型改為只在自由表面處設(shè)標(biāo)記點(diǎn)，對(duì)這些點(diǎn)的連續(xù)跟蹤可得到自由表面的準(zhǔn)確位置，且所耗機(jī)時(shí)大為減少。Hirt 和 Ni

5、chlos在總結(jié)當(dāng)時(shí)的自由表面跟蹤方法的基礎(chǔ)上還提出了“線段法”(Line Segment Method) 的概念，可近似確定自由表面的位置，但這種方法推廣到三維的情況較困難。MAC 法于 20世紀(jì) 60 年代提出后先是用于涌潮的傳播研究，后自 80 年代始大量用于波浪研究。如，Miyata運(yùn)用一種二維 MAC 格式研究了波浪破碎；Sakai等也提出一種基于早期 MAC 法的數(shù)值模型，來(lái)模擬水波潑濺以及波浪破碎過程中二階乃至三階漩渦的生成。Gao 和 Zhao運(yùn)用二維 MAC 模型，研究了波浪、建筑物及沙灘的交互作用。在 MAC 法之后，學(xué)者們開發(fā)了更具通用性和計(jì)算效率的模型。如，Chorin

6、提出的投影法 (Projection Method) 就是以其數(shù)學(xué)上的精確、簡(jiǎn)明而著稱的一個(gè)。以 MAC 法為基礎(chǔ)，Hirt 等和 Nichols 等開發(fā)了易于操作的 SOLA (solution-algorithm) 模型。Kothe等利用不完全 Cholesky共軛梯度法，將表面張力模擬為一種體積力，開發(fā)出了 RIPPLE模型。1.1.3 VOF 法（Volume-Of-Fluid）1.1.3.1 算法思想Nichols 等是最早用利用VOF法的研究群體。VOF 法假設(shè)在整個(gè)流動(dòng)區(qū)域內(nèi)流體密度為一常數(shù) ，空氣密度設(shè)為 0，引進(jìn)流體體積參數(shù) (或者是 VOF 函數(shù)) ,從連續(xù)方程出發(fā)建立關(guān)于

7、F的對(duì)流方程：(1-1)F 在流體單元內(nèi)為 1，在空單元中則為 0，而在自由表面單元中， F 值介于 0和 1 之間。給定速度場(chǎng)就可以通過上述對(duì)流方程確定F 值，進(jìn)而確定任一時(shí)刻的自由表面單元位置。1.1.3.2 方法介紹 供體受體法 (Donor-Acceptor Method)：由于 VOF 函數(shù)在自由表面附近有很大梯度，用一般的激波捕捉方法雖可求解上述 VOF 對(duì)流方程，但由于存在嚴(yán)重?cái)?shù)值耗散，難于精確確定自由表面位置。為此，Hirt 和 Nichols提出了供體受體法 (Donor-Acceptor Method)。該法首先依據(jù)流速在 x , y 方向的梯度大小確定自由水面的走向，并根

8、據(jù)流速的方向判斷自由水面附近的單元格為供體格還是受體格，并進(jìn)一步在運(yùn)動(dòng)方程求解中決定采用逆風(fēng)或順風(fēng)格式。這就有效地控制了數(shù)值耗散。有關(guān)數(shù)值實(shí)現(xiàn)的詳細(xì)信息可參閱 Hirt和 Nichlos、Kothe 、Liu 和 Lin等相關(guān)文獻(xiàn)。傾斜界面法 (Sloping Interface Method)：上述供體受體方法對(duì)自由表面幾何形狀進(jìn)行了簡(jiǎn)化處理，導(dǎo)致供體受體法精度不高。Youngs引入了傾斜界面法，較之于供體受體法，此類方法可更精確確定自由表面的位置。但目前在算法實(shí)現(xiàn)以及向三維擴(kuò)展時(shí)有一定困難，也就妨礙了傾斜界面法的實(shí)際應(yīng)用。層面設(shè)置法 (Level Set Method)：該法由 Sussm

9、an 等于 1994 年提出，是早期借助移動(dòng)網(wǎng)格跟蹤界面方法的改進(jìn)版本。該方法很容易確定分界面并且數(shù)值耗散較小。但由于這種格式?jīng)]有明確利用質(zhì)量守恒方程，因此復(fù)雜情況下容易造成質(zhì)量虧損。1.1.3.3 應(yīng)用情況 Delft 水力學(xué)研究小組自 20世紀(jì) 90年代早期就利用 VOF 方法進(jìn)行波浪動(dòng)力學(xué)研究。Van der Meer 等提出一種模型，研究作用在沿岸建筑物上的卷破波。Iwata 等研究了由于水下建筑物所致的波浪破碎變形。Lin 和 Liu k - e 利用紊流模型封閉 N-S方程，開發(fā)了一種基于 VOF 法的模型，對(duì)傾斜海灘上孤立波破碎等諸多水波動(dòng)問題進(jìn)行了研究。Sussman等耦合 V

10、OF 法和層面設(shè)置法，用于三維情況的不可壓縮流體的研究。上述方法雖然解決了自由表面的追蹤問題，但當(dāng)自由表面有劇烈變形時(shí)，各類方法在處理對(duì)流項(xiàng)時(shí)，數(shù)值耗散問題仍然都是不可避免的。1.2 無(wú)網(wǎng)格數(shù)值模擬方法無(wú)網(wǎng)格數(shù)值方法也稱為粒子方法，該方法完全摒棄了數(shù)值網(wǎng)格的概念，由于運(yùn)動(dòng)方程中不直接出現(xiàn)對(duì)流項(xiàng)，自然也就不會(huì)出現(xiàn)網(wǎng)格方法中的數(shù)值耗散現(xiàn)象。目前文獻(xiàn)上報(bào)道的粒子方法有 PAF法 (Particle-And-Force) 、無(wú)網(wǎng)格 Euler/N-S 方程解法、自由元 (Free Element) Galekin 法、離散元 (Discrete Element) 法、SPH(Smoothed- Par

11、ticle-Hydrodynamics) 法以及 MPS 法 (Moving Particle Semi-implicit Method)等。在水動(dòng)力學(xué)研究中開展較多的方法主要有 MPS 和 SPH 兩種，以下對(duì)其予以介紹。1.2.1 SPH 方法 SPH 方法的基礎(chǔ)是插值理論。SPH 通過某個(gè)插值函數(shù)在質(zhì)點(diǎn)上給出變量的近似值，從而將以偏微分方程形式給出的連續(xù)流體的守恒方程轉(zhuǎn)換成積分方程?？刂品匠讨械膲毫μ荻软?xiàng)、擴(kuò)散項(xiàng)等空間導(dǎo)數(shù)項(xiàng)也通過核子函數(shù)表示。Cummins 和 Rudman 提出了一種 SPH投影法 (SPH projection method)，通過求解壓力泊松方程確定壓力，從而確立

12、了具有嚴(yán)格算法的不可壓縮模型。SPH 方法源于天體物理學(xué)，后擴(kuò)展至水動(dòng)力學(xué)研究。隨后，Edmond 等用不可壓縮 SPH模型求解 Lagrange形式的 N-S 方程，并將模型與大渦數(shù)值模擬相耦合，以模擬近岸孤立波的作用機(jī)理。繼而，Shao 等將不可壓縮 SPH 模型用于模擬帶有自由表面的牛頓和非牛頓流體的研究。1.2.2 MPS方法 Koshizuka 等提出了一種粒子方法，用于計(jì)算不可壓縮黏性流體，并將之命名為移動(dòng)粒子半隱式方法 (Moving-Particle Semi-implicit method)，簡(jiǎn)稱為 MPS 法。該方法中，粒子間相互作用通過一個(gè)核子函數(shù)來(lái)表征，梯度和擴(kuò)散算子運(yùn)

13、算都通過核子函數(shù)表達(dá)。流體的不可壓縮性是通過控制粒子密度數(shù)不變實(shí)現(xiàn)的。計(jì)算過程分為顯式 (Explicit Stage) 和半隱式兩個(gè)階段。顯式階段僅考慮重力和黏滯項(xiàng)，半隱式階段用隱式方法求解壓力泊松方程，但速度的二次校正仍然為顯式。目前，該方法已用于波浪破碎機(jī)理和兩相流的研究中。另外，Gotoh 還將亞格子尺度紊流模型引入 MPS 模型，獲得了 N-S方程的 Lagrange 解?？傮w來(lái)看，水波動(dòng)力學(xué)中的無(wú)網(wǎng)格數(shù)值方法還處于初期探索階段，運(yùn)用于實(shí)際工程中的例子還少見報(bào)道。但由于其獨(dú)有的可很方便地追蹤自由表面的特性，決定了其在某些特定的工程實(shí)際問題中有網(wǎng)格方法無(wú)可比擬的優(yōu)越性，因而也具有很好的

14、發(fā)展前景。二、湍流及其數(shù)值模擬方法2.1 湍流湍流是流體的一種流動(dòng)狀態(tài)。當(dāng)流速很小時(shí)，流體分層流動(dòng)，互不混合，稱為層流，也稱為穩(wěn)流或片流；逐漸增加流速，流體的流線開始出現(xiàn)波浪狀的擺動(dòng)，擺動(dòng)的頻率及振幅隨流速的增加而增加，此種流況稱為過渡流；當(dāng)流速增加到很大時(shí)，流線不再清楚可辨，流場(chǎng)中有許多小漩渦，層流被破壞，相鄰流層間不但有滑動(dòng)，還有混合。這時(shí)的流體作不規(guī)則運(yùn)動(dòng)，有垂直于流管軸線方向的分速度產(chǎn)生，這種運(yùn)動(dòng)稱為湍流，又稱為亂流、擾流或紊流。2.2 湍流的數(shù)值模擬工程中絕大多數(shù)流體都是湍流。深層次地揭示流體運(yùn)動(dòng)的紊動(dòng)結(jié)構(gòu)及能量耗散過程是計(jì)算流體力學(xué)應(yīng)用研究的極其重要的課題。出于不同的研究目的，湍流

15、數(shù)值模擬有以下三個(gè)不同的層次：直接數(shù)值模擬 ( Direct Numerical Simulation，簡(jiǎn)稱 DNS)、雷諾平均數(shù)值模擬(Reynolds Averaged Navier-Stokes，簡(jiǎn)稱 RANS) 和大渦數(shù)值模擬 (Large Eddy Simulation，簡(jiǎn)稱 LES)。以下分別簡(jiǎn)要介紹。2.2.1 直接數(shù)值模擬由于 N-S 方程本身是封閉的，故從原則上講可以求解所有湍流問題。用 DNS 法直接求解 N-S 方程，能獲得最精細(xì)的流場(chǎng)信息。但在目前計(jì)算機(jī)發(fā)展水平下，進(jìn)行 DNS 的應(yīng)用研究是不現(xiàn)實(shí)的，目前只限于小規(guī)模的低雷諾數(shù)簡(jiǎn)單湍流物理機(jī)制的研究。DNS 常用的數(shù)值方

16、法是譜方法或偽譜法。Orszag 等最早用 DNS 計(jì)算了各向同性湍流。2.2.2 雷諾平均數(shù)值模擬由于雷諾平均方程是不封閉的，故 RANS 的核心思想是建立雷諾應(yīng)力封閉模型，使得平均運(yùn)動(dòng)方程可解。目前比較常用的模型有零方程模型、一方程模型 (k方程模型)、二方程模型 (渦黏性模型，k -e 模型)、代數(shù)應(yīng)力模型 (k -e - A模型)。其中k -e 模型是目前應(yīng)用最廣泛的湍流模型。在可以預(yù)見的將來(lái)，即使 DNS求解復(fù)雜湍流成為現(xiàn)實(shí) (耗時(shí)很長(zhǎng))，能夠快速地算出滿足一定工程精度要求的湍流統(tǒng)計(jì)模式仍會(huì)受到工程師們的青睞。事實(shí)上，RANS 是目前工程界處理湍流問題的唯一方法。2.2.3 大渦數(shù)值

17、模擬LES 的基本思想是：把湍流瞬時(shí)運(yùn)動(dòng)通過某種濾波方法分解為大尺度運(yùn)動(dòng)和小尺度運(yùn)動(dòng)兩部分。大渦運(yùn)動(dòng)通過直接求解 N-S 方程計(jì)算；小渦運(yùn)動(dòng)的影響概化為亞格子雷諾應(yīng)力，需通過建立模型求解，模型稱為亞格子尺度模型 (Subgrid Scale Model)。氣象學(xué)家 Deardoff 首次把大渦模擬用于有工程意義的槽道中的流體運(yùn)動(dòng)數(shù)值計(jì)算中。1972 年起，Stanford 大學(xué)的 Ferziger和 Reynolds 領(lǐng)導(dǎo)的集體開始對(duì) LES 做深入系統(tǒng)的研究。蘇銘德 1982 年以來(lái)提出了一種代數(shù)應(yīng)力模型，計(jì)算了槽道中的流體流動(dòng)。Watanabe 用大渦模擬計(jì)算了三維波浪破碎。LES 至今在

18、氣體動(dòng)力學(xué)中的研究開展較多，在水動(dòng)力學(xué)中的研究還十分有限，最主要的困難還不在于計(jì)算機(jī)的限制，而在其方法本身，如現(xiàn)有的亞格子尺度模型仍很不完善以及近壁模型的入流出流邊界等問題。三、水力機(jī)械空化數(shù)值模擬問題描述：高速離心泵在工作時(shí)，容易在葉輪葉片處產(chǎn)生汽蝕，會(huì)嚴(yán)重影響離心泵的性能。在葉輪入口處配置誘導(dǎo)輪,可以提高離心泵的抗汽蝕性能。3.1 兩相流基本理論在低比轉(zhuǎn)數(shù)離心泵發(fā)生汽蝕時(shí)，從氣泡的產(chǎn)生到消失時(shí)間極短，是一個(gè)復(fù)雜的動(dòng)態(tài)氣液兩相流動(dòng)過程。氣液兩相流流場(chǎng)的組成，從宏觀看，包括氣相、液相和氣液界面3部分。氣液界面的存在對(duì)兩相流的流動(dòng)特性有很大影響，流場(chǎng)中的界面在數(shù)量上和分布上經(jīng)常發(fā)生變化。 在等溫

19、兩相流中，隨流動(dòng)的進(jìn)行，壓力呈現(xiàn)變化，兩相界面的大小和分布狀況即發(fā)生變化。當(dāng)壓力變化較小時(shí)，沿流道壓降與介質(zhì)壓強(qiáng)相比，數(shù)值很小，界面及其分布可以看成是不變的。但在壓力變化很大的兩相流中，介質(zhì)壓力大小和氣液兩相界面的變化均很顯著，因此流動(dòng)特性也會(huì)隨之改變。氣液兩相流的數(shù)值模擬包括氣相場(chǎng)和氣液間的相互干擾計(jì)算。相間干擾即氣液兩相間的動(dòng)量、能量、質(zhì)量的交換過程。常見算法可以分為歐拉一歐拉型算法和歐拉一拉格朗日算法兩種類型。3.2 兩相流控制方程連續(xù)方程：(3-1)式中：液滴蒸發(fā)進(jìn)入氣相的質(zhì)量動(dòng)量方程：(3-2)式中：F離散相粒子對(duì)氣流的作用力3.3 氣泡項(xiàng)體積分?jǐn)?shù)方程(3-3)其中： 3.4 變螺距誘導(dǎo)輪的設(shè)計(jì)3.4.1 誘導(dǎo)輪提高泵抗汽蝕性能的原因高抗汽蝕性能的誘導(dǎo)輪在一定的程度上是以犧牲能量指標(biāo)換來(lái)的，其效率并不高，但誘導(dǎo)輪本身負(fù)荷不大，可以改變?nèi)~輪進(jìn)口流動(dòng)情況，不明顯降低水泵的效率。(a)無(wú)誘導(dǎo)輪 (b)有誘導(dǎo)輪圖1 葉輪進(jìn)口速度三角形圖2 誘導(dǎo)輪模型由上圖可知，加誘導(dǎo)輪之后，主葉輪葉片進(jìn)口前相對(duì)速度和無(wú)誘導(dǎo)輪相比減小了。則主葉輪的汽蝕余量減小，泵的抗汽蝕性能提高。同時(shí)加誘導(dǎo)輪之后主葉輪葉片進(jìn)口之前液體的絕對(duì)速度和未加誘導(dǎo)輪時(shí)相比也增加了，此的增加由誘導(dǎo)輪作用產(chǎn)生，并不引起壓力下降。因此可以認(rèn)為不影響主葉輪內(nèi)最低