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

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

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

0 引 言

隨著CFD和計算機的飛速發展及廣泛應用，許多流體力學的復雜問題都通過計算機數值模擬進行解決，有限元法就是眾多數值計算方法的一種［1］。

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

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

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

1 問題的提出

1.1 有限元法簡介

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

1.2 問題的物理模型

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

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

坐標系按照圖中所示建立，無窮遠處來流速度和壓強分別為V∞=1，P∞=0，分別求解網格密度為1，0.5，0.1，0.05 下的流場流函數和壓強分布。

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

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

1.4 邊界條件

區域Ω和相關邊界條件值如下:

2 有限元數值求解

2.1 積分表達式

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

應用Green－Gauss 公式可將上述積分改寫為:

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

2.2 單元剖分

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

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

2.3 確定單元基函數

單元及結點數確定后，單元基函數也確定了。本算例中采用三角形三結點的單元基函數，單元基函數的線性插值函數為:

三結點分別對應3 個基函數。記e 單元的三角形的3 個結點坐標為;則式(4)所滿足的插值條件為:

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

其中下標i，j，k 按1，2，3 順序循環取值。求解得到:

其中:

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

2.4 單元分析

注意到有如下變形式:

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

其中:

2.5 方程求解

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

各個單元速度的求解:

求出單元的速度后，進而可以得到各個單元節點的速度，通過Bernouli 方程計算出各節點的壓力值，假設求解區域位于同一水平面內，介質密度ρ=1 來流壓力P=0，那么結點壓力:

3 計算結果及分析

3.1 流線圖和壓力云圖

通過編制程序，得到點間隔1，0.5，0.1，0.05 的網格相關數據節點個數和單元個數，如表1所示。

表1 各個網格密度節點和單元數Tab.1 Node and element number of each mesh density

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

圖2 網格密度流線圖Fig.2 Flow pattern when the mesh density equals

圖3 網格密度壓力圖Fig.3 Pressure profile when the mesh density equals

3.2 不同網格密度對壓力的影響

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

3.3 流線圖和壓力云圖

從結果輸出圖中可以得到以下分析結果:

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

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

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

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

4 結 語

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

［1］章本照．流體力學中的有限元方法［M］．北京:機械工業出版社，1986．

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

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

［4］岳蕾，張志國，蔣奉兼．圓球繞流的大渦模擬分析研究［C］．第十一屆全國水動力學學術會議暨第二十四屆全國水動力學研討會，無錫，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］章本照，印建安，張宏基．流體力學數值方法［M］．北京:機械工業出版社，2003．