圓球勢流繞流的有限元方法

黃 璐，陳 立，邱遼原，許 輝

(中國艦船研究設(shè)計中心，湖北 武漢430064)

0 引 言

隨著CFD和計算機(jī)的飛速發(fā)展及廣泛應(yīng)用，許多流體力學(xué)的復(fù)雜問題都通過計算機(jī)數(shù)值模擬進(jìn)行解決，有限元法就是眾多數(shù)值計算方法的一種［1］。

在工程實(shí)踐中，鈍體繞流是一個十分典型的問題。船舶與海洋工程中就有很多例子可以簡化為鈍體繞流的模型進(jìn)行研究，如推進(jìn)器，海洋平臺的立管和柱體繞流等。圓球繞流是鈍體繞流的典型例子。

圓球繞流的問題已有許多研究，吳望一等［2］利用有限元法和傳統(tǒng)的分析方法相結(jié)合，研究了無界域上理想流體的圓球繞流問題。王吉飛和萬德成［3］利用有限元法和并行計算相結(jié)合，研究了理想不可壓流體的二維和三維鈍體勢流繞流問題，給出了二維無限域圓柱繞流和三維無限域圓球繞流的數(shù)值模擬結(jié)果，并做了并行效率的分析。岳蕾、張志國等［4］利用大渦模擬等湍流模型，分析了圓球繞流的尾部流場及圓球表面的壓力變化。Hanazaki［5］通過差分方法求解三維非穩(wěn)態(tài)N－S 方程，得到一些圓球繞流結(jié)果(Re ＜200)。Sungsu Lee［6］用有限元方法模擬Re 數(shù)在100 ～500 之間的圓球繞流。

本文利用Fortran 語言自編程序，采用有限元法求解圓管內(nèi)圓球在理想不可壓縮條件下的流動問題，將求解區(qū)域劃分成離散的網(wǎng)格，采用1，0.5，0.1，0.05 等不同的網(wǎng)格密度，比較不同的網(wǎng)格密度對有限元法求解圓球繞流問題的影響，為有限元求解圓球繞流問題時采用合適的網(wǎng)格密度提供參考。

1 問題的提出

1.1 有限元法簡介

有限元法利用泛函變分原理或方程余量與權(quán)函數(shù)正交化原理。利用“分塊逼近”進(jìn)行離散求解［7］，將流場的求解區(qū)域劃分成互不重疊的單元;在每個單元內(nèi)，選擇若干點(diǎn)作為求解函數(shù)的插值點(diǎn)。單元中的函數(shù)將由線性化的基函數(shù)取代，然后通過單元體的積分就可得到單元有限元方程，經(jīng)過累加獲得總體有限元方程，通過求解總體有限元方程得到節(jié)點(diǎn)上的函數(shù)值，從而求解得到需要的物理量。

1.2 問題的物理模型

考慮位于圓管內(nèi)無限長的圓球體的繞流問題，幾何尺寸如圖1所示。無窮遠(yuǎn)處來流為Vz=1，Vy=0。由于流場具有上下左右的軸對稱性，取圓球中心處的縱向剖面，且只考慮左上角1/4 的計算區(qū)域a－b－c－d－e，把它作為有限元的求解區(qū)域Ω。

圖1 求解問題的物理模型Fig.1 Physical model of the problem

坐標(biāo)系按照圖中所示建立，無窮遠(yuǎn)處來流速度和壓強(qiáng)分別為V∞=1，P∞=0，分別求解網(wǎng)格密度為1，0.5，0.1，0.05 下的流場流函數(shù)和壓強(qiáng)分布。

1.3 物理問題滿足的微分方程式

由于圓管內(nèi)流動具有上下、前后對稱特點(diǎn)，計算區(qū)域Ω 取流場的1/4 區(qū)域;Γ1是Ω 的本質(zhì)邊界，由進(jìn)口邊界a－e，固壁邊界e－d，圓球面固壁邊界a－b－c 組成;Γ2是自然邊界，它是前后對稱軸線cd。

1.4 邊界條件

區(qū)域Ω和相關(guān)邊界條件值如下:

2 有限元數(shù)值求解

2.1 積分表達(dá)式

對流函數(shù)ψ 滿足的Stokes 方程采用Galerkin 方法，可得到:

應(yīng)用Green－Gauss 公式可將上述積分改寫為:

其中:Ω 為(r，z)平面上的求解區(qū)域;Γ2為自然邊界。

2.2 單元剖分

根據(jù)物理問題的特點(diǎn)以及區(qū)域的形狀，把計算區(qū)域分成許多幾何形狀規(guī)則但大小可以不同的單元，確定單元節(jié)點(diǎn)的數(shù)目和位置，建立表示網(wǎng)格的數(shù)據(jù)結(jié)構(gòu)。采用的單元形狀和節(jié)點(diǎn)的分布，以及插值函數(shù)的選取還應(yīng)考慮到計算精度和可微性的要求。本題中根據(jù)給定的網(wǎng)格密度劃分相應(yīng)的網(wǎng)格，采用三角形結(jié)構(gòu)化網(wǎng)格。為使形成的總體有限元方程組系數(shù)矩陣具有較小的帶寬，區(qū)域節(jié)點(diǎn)序號沿著區(qū)域短邊方向，自下而上逐列編號。每個單元節(jié)點(diǎn)序號采用逆時針的編號方法，最后列出總體序號與單元序號之間的表格。

除此之外，單元剖分還要建立本質(zhì)邊界條件和自然邊界條件的節(jié)點(diǎn)表。自然邊界條件對解此題并無貢獻(xiàn)，由網(wǎng)格表特征可知，本題的本質(zhì)邊界條件為a－b－c 段上，ψ 為0;d－e 段上，ψ 為2;e－a段上，ψ 為r2/2。

2.3 確定單元基函數(shù)

單元及結(jié)點(diǎn)數(shù)確定后，單元基函數(shù)也確定了。本算例中采用三角形三結(jié)點(diǎn)的單元基函數(shù)，單元基函數(shù)的線性插值函數(shù)為:

三結(jié)點(diǎn)分別對應(yīng)3 個基函數(shù)。記e 單元的三角形的3 個結(jié)點(diǎn)坐標(biāo)為;則式(4)所滿足的插值條件為:

將式(5)代入式(4)，即得到含有9 個待定系數(shù)ai，bi，ci(i=1，2，3)的9 個代數(shù)方程組:

其中下標(biāo)i，j，k 按1，2，3 順序循環(huán)取值。求解得到:

其中:

式中A(e)為e 單元三角形的面積。

2.4 單元分析

注意到有如下變形式:

其中(zc，rc)為單元三角形的形心坐標(biāo)。將式(4)和式(10)代入式(9)，即可推導(dǎo)出有限元方程:

其中:

2.5 方程求解

將各個單元的方程按照順序合成后形成總體方程，在形成的總體方程中利用消行修正法處理本質(zhì)邊界條件，得到修正方程。利用Gauss 列主元素消元法直接求解方程組。求出所有的待求量后，便得到了近似函數(shù)的表達(dá)式，并可以計算出相關(guān)的物理量。

各個單元速度的求解:

求出單元的速度后，進(jìn)而可以得到各個單元節(jié)點(diǎn)的速度，通過Bernouli 方程計算出各節(jié)點(diǎn)的壓力值，假設(shè)求解區(qū)域位于同一水平面內(nèi)，介質(zhì)密度ρ=1 來流壓力P=0，那么結(jié)點(diǎn)壓力:

3 計算結(jié)果及分析

3.1 流線圖和壓力云圖

通過編制程序，得到點(diǎn)間隔1，0.5，0.1，0.05 的網(wǎng)格相關(guān)數(shù)據(jù)節(jié)點(diǎn)個數(shù)和單元個數(shù)，如表1所示。

表1 各個網(wǎng)格密度節(jié)點(diǎn)和單元數(shù)Tab.1 Node and element number of each mesh density

不同網(wǎng)格密度得到的流線圖形如圖2所示，不同網(wǎng)格密度得到的壓力云圖如圖3所示。

圖2 網(wǎng)格密度流線圖Fig.2 Flow pattern when the mesh density equals

圖3 網(wǎng)格密度壓力圖Fig.3 Pressure profile when the mesh density equals

3.2 不同網(wǎng)格密度對壓力的影響

為了比較不同的網(wǎng)格密度對計算結(jié)果的影響，選取壓力P 在不同網(wǎng)格密度下，分別在y=0，y=1，y=2 處計算結(jié)果的比較，比較結(jié)果如圖4所示(圖中t 代表網(wǎng)格密度)。

3.3 流線圖和壓力云圖

從結(jié)果輸出圖中可以得到以下分析結(jié)果:

1)從流線圖中可以看到，隨著網(wǎng)格密度的增加，流線變得更加光滑，更加均勻，說明網(wǎng)格密度越密得到的流場流線圖越接近真實(shí)的流場。從圖2中可以看到，y=0 時是0 流線，隨著y 值的增加，流線的值也在增加，在y=2 時達(dá)到最大。

2)從壓力云圖中可以看到，在圓球前端，壓力最大。隨著圓球的壁面向上，壓力值慢慢變小，直達(dá)圓球頂端，壓力達(dá)到最小。這與根據(jù)伯努利定理分析速度矢量得到的結(jié)果基本符合。

3)由圖4 可以看出，在網(wǎng)格密度為1和0.5時，由于網(wǎng)格較稀疏，相關(guān)物理量的變化還比較大;在網(wǎng)格密度為0.05和0.1 時，相關(guān)物理量的變化比較小，說明隨著網(wǎng)格的加密計算趨于穩(wěn)定。

圖4 y=0，1，2 時不同網(wǎng)格密度壓力P 變化圖Fig.4 Pressure profile of different mesh density when y=0，1，2

4 結(jié) 語

以上結(jié)果與文獻(xiàn)［3］中計算結(jié)果比較，符合度較好，說明文中所用的有限元法適合于解決相關(guān)的流體力學(xué)問題，為不同密度網(wǎng)格的有限元法在解決流體力學(xué)問題研究時提供參考。以后可以進(jìn)一步研究粘性流動下有限元解法，也可同CFD 模擬作進(jìn)一步的比較。

［1］章本照．流體力學(xué)中的有限元方法［M］．北京:機(jī)械工業(yè)出版社，1986．

［2］吳望一，梁鐳，溫功碧，等．無界區(qū)域上理想繞流問題的混合元方程［J］．空氣動力學(xué)學(xué)報，1989，7(4):449－454．

［3］王吉飛，萬德成．鈍體勢流繞流的有限元并行計算［C］．第二十一屆全國水動力學(xué)研討會暨全國第八屆水動力學(xué)學(xué)術(shù)會議，成都，2008:227－233．

［4］岳蕾，張志國，蔣奉兼．圓球繞流的大渦模擬分析研究［C］．第十一屆全國水動力學(xué)學(xué)術(shù)會議暨第二十四屆全國水動力學(xué)研討會，無錫，2012:237－244．YUE lei，ZHANG Zhi-guo，JIANG Feng-jian．Investigation of flow cross sphere using large eddy simulation［C］．The 11st Session of the National Seminar on Hydrodynamics and National Eighth Hydrodynamics Academic Conference，Wuxi，2012:237－244．

［5］HANAZAKIH A．A numerical study of three-dimensional stratified flow past a sphere［J］．Journal of Fluid Mechanics，1998，192:393－419．

［6］SUNGSU L．A numerical study of the unsteady wake behind a sphere in uniform flow at moderate Reynolds numbers［J］．Computer and Fluids，2000，29:639－667．

［7］章本照，印建安，張宏基．流體力學(xué)數(shù)值方法［M］．北京:機(jī)械工業(yè)出版社，2003．