一種艦船流場特性計算方法研究*

張永坤 王樹樂

(91439部隊 大連 116041)

?

一種艦船流場特性計算方法研究*

張永坤 王樹樂

(91439部隊 大連 116041)

以DTMB5415為研究對象,采用RANS方法結合VOF模型,研究了網格因素以及湍流模型對阻力、波形和流場預報精度的影響。首先設計了粗糙、中等和精細共三套網格,分別進行興波流場計算,通過比較試驗測試的阻力和波形數據,確定中等網格劃分方式既能保證計算的準確性又能保持較高的計算效率。其次選擇三種湍流模型,對總阻力及舷側波高進行數值計算,SST k-ω湍流模型在保證計算效率的情況下可滿足計算精度要求。最后采用中等網格劃分以及SST k-ω模型進行計算,并同試驗測量結果進行比對,結果表明船側表面附近數值計算波形分布與試驗結果基本一致。

流場; 計算網格; 流體體積函數; 數值計算

Class Number U661.2

1 引言

船體是一個極其復雜的曲面,由于流體的粘性作用,船舶航行時產生邊界層,形成非常不均勻的三向伴流場。船舶伴流場的準確獲取是螺旋槳設計、船體水動力和噪聲性能預報等研究的基礎條件。鑒于此,開展水面艦船伴流場數值計算研究,分析艦船阻力、特征截面伴流場分布,對于提高水面艦船水動力性能的預報精度具有重要意義。使用CFD方法進行流場計算時,首先要將計算域進行離散,數值計算是在離散的網格節點上滿足流體力學方程,因此網格分辨率對數值計算結果有較大的影響。網格數太少,太稀疏難以得到精確的結果[1],網格數太大、過密又會增加計算量。

本文以DTMB5415為研究對象,采用RANS方法[2～7]結合VOF模型[9～10]研究了網格和湍流模型等因素對阻力、波形和流場預報精度的影響。首先進行網格依賴性研究,采取全結構化網格劃分方法,固定某種湍流模型,設計了粗糙、中等和精細共三套網格,分別進行興波流場計算,通過比較試驗測試的阻力和波形數據,確定中等網格劃分既能保證計算的準確性又能保持較高的計算效率。在此基礎上,選取理論發展比較完善、船舶粘性流場數值計算常用的幾種湍流模型,RNGK-ε和SST k-ω湍流模型和雷諾應力模型進行湍流模型的影響研究,通過與試驗數據的對比,三種湍流模型均能較準確地預報該船的阻力,而在船體舷側波形預報上,SST k-ω湍流模型和雷諾應力模型更具有優勢。

2 計算對象

本文采用DTMB5415為研究對象,進行數值計算方法的驗證。DTMB5415為美國海軍驅逐艦(DDG51)的模型,自1996年后,DTMB5415成了ITTC用于船舶CFD水動力計算驗證的表準模型之一,意大利和美國對它做了一系列試驗,且公布了完整的試驗報告,有詳細的阻力、波形和流場試驗數據。其外形輪廓如圖1所示,主尺度見表1。

圖1 DTMB5415的外形輪廓

項目實尺度模型尺度λ24．824Lpp(m)1425．72B(m)19．060．76D(m)6．150．248Δ(m3)8424．40．549Sw(m2)2972．64．786

根據哥德堡2010會議公布的測試案例,選取的計算狀態為:固定升沉1.82×10-3,縱傾(-0.108°),航速2.097m/s,雷諾數為1.19×107,長度傅汝德數為0.28。計算平臺為ANSYS FLUENT14.0。

3 計算方法

計算考慮興波,即帶自由液面,自由液面由于要考慮空氣和水的相互作用,所以屬于多項流問題。采用VOF方法處理自由液面問題,以流體占據網格單元體積份額的途徑來跟蹤自由面演化,具體算法使用固定網格系統捕捉多種相參混流體的交界面,所有流體滿足同一動量方程,在整個計算區域上跟蹤每一種流體在每個計算單元中的體積分數。

假設第q種流體在單元中的體積分數為αq,則存在以下三種可能: 1)αq=0,單元內不存在第q種流體; 2)αq=1,單元內充滿第q種流體; 3) 0<αq<1,單元內存在不同流體的交界面。

在VOF方法中,αq定義一個單元中流體所占有的體積分數,用αq還可以確定交界面的位置,因為在交界面法向上αq變化最快。在確定交界面法向以及體積分數αq的值后,單元中就可以確定一條分割線,用來近似表達兩種液體的交界面。

為了跟蹤流體之間界面的位置,需要求解體積分數的連續性方程,對q種流體有:

(1)

各種流體體積分數滿足式(2):

(2)

每個控制體中,流體密度由體積分數平均得到:

ρ=∑αqρq

(3)

其它流體性質如粘性也由類似的方法計算。

在整個計算區域只要計算單個動量方程,動量方程通過密度ρ和粘性系數μ依賴于各相流體的體積分數,動量方程如式(4):

(4)

對于其它輸運量,如湍動能、耗散率或者雷諾應力等,所有流體也只要求解一組輸運方程。

微分方程采用有限體積法離散時,數值計算必須控制各個面上的對流和擴散通量,以便和單元內源相平衡,這里對單元表面通量的計算采用的是精度較高的幾何重構法。幾何重構法使用分段線性的方法描繪流體之間的界面,它假設兩種流體之間的界面在每個單元內有個線性斜面,并使用這個線性形狀來計算穿過單元面的流體對流。

4 湍流模型及邊界條件設置

4.1 船體網格劃分

采用全結構化網格離散計算域[11],網格采用H-O拓撲形式,在球鼻艏和艉部流動變化較大的區域進行加密,保證船中近壁面第一層網格單元y+=80左右,船體表面網格如圖2所示。

圖2 DTMB5415表面網格劃分

自由液面興波流場的計算區域從靜水面處分為兩個部分,水域和空氣層。上游入口離船艏1倍船長,下游出口離船尾2.5倍船長,空氣層的厚度為1.5倍吃水高度,外側面離船中1倍船長。

4.2 邊界條件設置

邊界條件的設置如下:上游入口采用速度入口,給定均勻來流的速度值;下游出口采用FLUENT提供的用戶接口,編寫用戶UDF函數,給定靜壓出口,忽略擾動;上邊界、側邊和中間面均采用對稱面邊界條件,根據計算經驗,上邊界和側邊采取對稱面與無擾動的速度來流條件效果是一樣的;在船體表面定義無滑移、不可穿透的邊界條件。

4.3 湍流模型選取

湍流模型對船舶粘性流場的計算精度有重要的影響。選取理論發展比較完善、船舶粘性流場數值計算常用的幾種湍流模型,RNGK-ε和SST k-ω湍流模型和雷諾應力模型。

1) RNG k-ε模型

為了準確描述強旋流或帶有彎曲壁面的流動,Yakhot等提出了重整化群k-ε(即RNG k-ε)模型,該模型在時均連續方程和動量守恒方程的基礎上,引入湍動能k和湍動耗散率ε的方程,再建立湍動粘度μt與k和ε的關系式,使方程組封閉。

2) SST k-ω模型

改進二：熱源改酒精燈為恒溫調奶器（溫度可在37～100℃之間任意恒定）；改水浴加熱為80℃恒溫直接加熱；改分組脫色為全班集中脫色。恒溫調奶器為內熱熱源，沒有明火，除去了明火熱源可能點燃酒精引發火災的安全隱患。加熱溫度恒定在酒精沸點溫度之下，因此不會出現因酒精沸騰飛濺到學生身上造成的意外傷害。溫度恒定在酒精沸點之下，可以減慢酒精氣化速度，節約酒精的同時也加快了脫色速度。

為了準確求解分離流問題,Menter提出了SST k-ω湍流模型,該模型對壁面處流場應用k-ω模型而對總體流動則應用k-ε模型求解,兩模型間的過渡通過混合函數完成。

3) 雷諾應力模型

渦粘模型采用各向同性的湍動粘度來計算湍流應力,這種假設往往導致此類模型難以準確考慮旋轉流動和流動方向表面曲率變化的影響,當流場中存在高曲率流動、逆壓梯度、分離流及強旋轉流動時,湍動粘度的各向異性影響較大,導致應用渦粘模型模擬的流場分布與實驗結果明顯不符。為了克服這個缺點,需要對RANS方程中的Reynolds應力項直接建立微分方程式進行求解,以此得到Reynolds應力的各個分量值。

采用有限體積法離散控制方程和湍流模式,壓力項采用body-forced-weighted格式離散,其他項均采用二階迎風格式離散,體積分數采用幾何重構方案,壓力速度耦合迭代采用SIMPLEC方法。

5 結果分析

5.1 網格依賴性研究

為了研究網格依賴性,共設計了粗糙、中等、精細共三套網格,三套網格的拓撲形式完全一致,僅僅是長度方向、展向和腰身方向的節點數量不一致,三套網格的加細率rG=1.4142,三套網格的具體設置如表2所示。

表2 三套網格節點布置情況

圖3為中等網格總阻力系數收斂的時間歷程曲線,從圖中可以看出計算時間達到10s后,總阻力系數成比較穩定的周期性變化,取穩定周期內總阻力的時歷平均值作為計算結果。

圖3 總阻力系數收斂曲線

三套網格數量下,Fn=0.28時DTMB5415裸船體的總阻力系數數值計算結果與試驗數據的對比如表3所示。

表3 不同網格數量的總阻力系數計算結果與試驗數據的對比

從表3中可以看出:隨著網格數量的增加,阻力計算精度也相應增加。在此傅汝德數下,中等網格的計算精度已足夠工程應用范圍。

采用26屆ITTC推薦規程對阻力計算的網格收斂性進行判斷。粗糙網格計算結果記為sGC;中等網格計算結果記為sGM;精細網格計算結果記為sGF;計算粗糙網格與中等網格、中等網格與精細網格結果之差為

εGMC=SGM-SGC=4.257-4.352=-0.095

εGFM=SGF-SGM=4.246-4.257=-0.011

網格收斂因子:

由于網格收斂因子〈RG〉小于1,說明計算網格的收斂性得到滿足,表明計算網格數越大數值計算結果越精確,不確定度越小。

進一步比較波形對網格數量的依賴性。圖4和圖5分別為三種網格數量下y/LPP=0.082,y/LPP=0.172船側波高的計算值與試驗值的對比。

圖4 y/LPP=0.082時船體舷側波高沿船長分布曲線

圖5 y/LPP=0.172時船體舷側波高沿船長分布曲線

從圖中可以看出:當非常靠近船體表面時,即y/LPP=0.082,三套網格均能很好地模擬出船體舷側波高沿船長的分布,波峰和波谷的位置捕捉得非常精確,具體幅值稍有差異,整體來說,精細網格分布效果最好。而稍微遠離船體表面,即y/LPP=0.172時,它們之間的差異變大,粗糙網格模擬的船體舷側波高在船中部分已不能捕捉二次小波谷,在艉部捕捉的波谷值也過小,而中等和精細網格依然能較好地模擬船體舷側波高沿船長的分布,波峰和波谷的位置非常精確,幅值大小也相差無幾。

綜合阻力和興波波形的計算結果,同時考慮到計算資源,湍流模型研究可采取中等網格劃分形式。

5.2 湍流模型影響研究

1) 總阻力

以DTMB5415為對象,采取中等網格劃分,分別結合以上三種湍流模型對傅汝德數Fr=0.28時的艦船興波流場進行求解。計算結果如表4所示。

可以看出:RNGK-ε模型比試驗值稍低而雷諾應力模型則比試驗值稍高,總體來說,三種湍流模型預報的總阻力精度均很高。

表4 不同湍流模型的總阻力系數計算結果與試驗數據的對比

2) 船體舷側波高

進一步比較這幾種湍流模型計算得到的船體舷側波高。圖6和圖7分別為三種湍流模型y/LPP=0.082,y/LPP=0.172船側波高的計算值與試驗值的對比。

圖6 y/LPP=0.082時船體舷側波高沿船長分布曲線

圖7 型y/LPP=0.172時船體舷側波高沿船長分布曲線

從圖中可以看出:這三種湍流模型均能較好地預報船體舷側波形,波峰和波谷的位置精確。總體而言,雷諾應力模型效果最好,SST k-ω模型次之,RNGK-ε模型稍微差點。綜合計算精度和計算效率,可選擇SST k-ω模型進行計算。

5.3 SST k-ω模型計算結果與試驗值比對

由圖8可知,在船側表面附近,數值計算波形分布與試驗數據吻合得相當好。而遠離船體表面區域特別是艉部,由于網格的稀疏而導致波形消弱。

圖8 SST k-ω模型計算波高等值線與試驗數據的對比

觀察圖9可以發現,由于球鼻艏的存在而產生的漩渦向后發展,在此截面依然很強,這與試驗觀測是一致的;同時CFD還顯示了此截面另外存在舭渦,而試驗并未捕捉到這個現象。

圖9 x/LPP=0.2截面速度場的計算值與試驗值的比較

圖10 x/LPP=1.1截面速度場的計算值與試驗值的比較

圖10為方艉正后方的某個截面。數值計算結果成功地捕捉到了方艉的流動現象,從圖中可以看出,流動產生的主渦被正確地模擬。與試驗結果相比,數值計算結果有兩點差異。一是主渦中心位置更加靠近船體中縱剖面;二是船體底部的橫向流動亦服從主渦的流動行為,而試驗結果則顯示該部分橫向流動在渦流之外。產生這種現象的原因可能是CFD計算的流線更易遭受漩渦的影響。

6 結語

通過本文研究得到一些有意義的結論:

1) 計算網格數越大數值計算結果越精確,不確定度越小。

2) 稍微遠離船體表面,粗糙網格模擬的船體舷側波高在船中部分已不能捕捉二次小波谷,在艉部捕捉的波谷值也過小。

3) 通過比較試驗測試的阻力和波形數據,中等網格劃分既能保證計算的準確性又能保持較高的計算效率。在滿足計算精度的條件下,采用中等數量計算網格可節省計算資源。

4) 就總阻力預報而言,RNGK-ε和SST k-ω湍流模型和雷諾應力模型預報精度都很高,滿足計算要求。

5) 從預報船體舷側波形來講,三種湍流模型均能較好地預報,波峰和波谷的位置精確。總體而言,雷諾應力模型效果最好,SST k-ω模型次之,RNGK-ε模型稍微差點。綜合考慮計算精度及計算效率,SST k-ω模型最優。

6) 選取DTMB5415計算對象,采用中等網格劃分以及SST k-ω模型進行計算,計算結果同試驗測量結果進行比對。結果表明船側表面附近數值計算波形分布與試驗結果一致,同時捕捉到球鼻艏導致的漩渦發展、舭渦以及方艉的流動等現象。

[1] 劉志華,熊鷹.潛艇流場數值計算網格與湍流模型選取[J].華中科技大學學報(自然科學版),2009,37(7):98-101.

[2] 黃振宇,周連第.帶全附體潛艇尾流場的數值預報研究[J].中國造船,2001,42(2):6-11.

[3] 陳明周.基于CFD的潛艇尾流場波數預報研究[D].武漢:海軍工程大學船舶與動力學院,2007.

[4] 吳召華,陳作鋼.基于體積力法的船體自航性能數值預報[J].上海交通大學學報,2013,47(6):943-949.

[5] 沈興榮,馮學梅,蔡榮泉.大型集裝箱船槳舵干擾黏性流場的數值計算[J].船舶力學,2009,13(4):540-550.

[6] 邱遼原,侯國祥.軸對稱回轉體繞流場數值模擬[J].華中科技大學學報,自然科學版,2004,32(10):46-48.

[7] Stern F, Kim H T, Patel V C, et al. A viscous flow approach to the computation of propeller-hull interaction[J]. Journal of Ship Research,1988,32(4):246-262.

[8] Hirt C W, Nichol B D. Volume of fluid (VOF) method for the dynamics of free boundaries[J]. J. Comp. Phys.,1981,39:201-221.

[9] 李誼樂,劉應中.二階精度的VOF自由面追蹤方法及其應用[J].船舶力學,1999,3(1):44-52.

[10] 萬鵬程,傅慧萍.艦船氣泡尾流數值模擬[J].上海交通大學學報,2013,47(2):193-202.

[11] 張振江,陳作鋼.船舶CFD結構化網格自動生成技術的開發[J].上海交通大學學報,2012,46(2):317-322.

A Method of Ship Flow Field Numerical Calculation

ZHANG Yongkun WANG Shule

(No. 91439 Troops of PLA, Dalian 116041)

The influence of grid division and turbulent model to resistance, wave ship and precision of flow prediction has been investigated by selecting ship model DTMB5415 as calculated object and using RANS method and VOF model. First, coarseness, middle and huge grid division have been designed to calculate wave flow field. By comparing resistance and wave ship date of the test with calculating data, the results show that middle count grid divisions can obtain veracity of calculate results and higher calculate efficiency. Second, three turbulent flow models have been used to compute total resistance and ship wave. The SST k-ω turbulent flow model can satisfy calculated results accuracy under high calculate rate. Finally, middle grid division and SST k-ω turbulent flow model have been selected to calculated ship wave. The results showed that the calculated results have been agreed with the experiment data.

flow field, computational grid, VOF method, numerical calculation

2015年4月12日,

2015年5月27日

張永坤,男,博士,工程師,研究方向:水下爆炸。

U661.2

10.3969/j.issn.1672-9730.2015.10.040