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

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

0 引 言

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

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

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

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

1 問題的提出

1.1 有限元法簡(jiǎn)介

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

1.2 問題的物理模型

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

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

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

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

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

1.4 邊界條件

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

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

2.1 積分表達(dá)式

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

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

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

2.2 單元剖分

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

除此之外，單元剖分還要建立本質(zhì)邊界條件和自然邊界條件的節(jié)點(diǎn)表。自然邊界條件對(duì)解此題并無貢獻(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)分別對(duì)應(yīng)3 個(gè)基函數(shù)。記e 單元的三角形的3 個(gè)結(jié)點(diǎn)坐標(biāo)為;則式(4)所滿足的插值條件為:

將式(5)代入式(4)，即得到含有9 個(gè)待定系數(shù)ai，bi，ci(i=1，2，3)的9 個(gè)代數(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 方程求解

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

各個(gè)單元速度的求解:

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

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

3.1 流線圖和壓力云圖

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

表1 各個(gè)網(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)格密度對(duì)壓力的影響

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

3.3 流線圖和壓力云圖

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

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

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

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

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

4 結(jié) 語(yǔ)

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

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

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

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

［4］岳蕾，張志國(guó)，蔣奉兼．圓球繞流的大渦模擬分析研究［C］．第十一屆全國(guó)水動(dòng)力學(xué)學(xué)術(shù)會(huì)議暨第二十四屆全國(guó)水動(dòng)力學(xué)研討會(huì)，無錫，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．