基于CFD的方柱繞流水動(dòng)力數(shù)值模擬

劉志江，郭建廷，孟小峰，劉浩然

(1. 華中科技大學(xué) 船舶與海洋工程學(xué)院，湖北 武漢 430074；2. 江蘇科技大學(xué) 船舶與海洋工程學(xué)院，江蘇 鎮(zhèn)江 212000)

0 引 言

出于減噪的目的，一些船舶設(shè)置有自循環(huán)系統(tǒng)，利用船舶航行時(shí)迎流面所產(chǎn)生的動(dòng)壓頭抽吸海水，降低冷凝器的用水量減少水泵的運(yùn)行次數(shù)[1-2]。然而在系統(tǒng)與海水連通處，為防止各種雜物進(jìn)入系統(tǒng)管道，通常在入口處設(shè)置格柵，當(dāng)自循環(huán)系統(tǒng)工作時(shí)，海水流經(jīng)格柵會(huì)產(chǎn)生剪切層分離等現(xiàn)象，從而在格柵表面及附近流場(chǎng)產(chǎn)生強(qiáng)大的壓力脈動(dòng)場(chǎng)，分析格柵的壓力脈動(dòng)特性，對(duì)船舶降噪具有重大意義。

本文基于計(jì)算流體力學(xué)（CFD）分析方法，在不同雷諾數(shù)Re、不同長(zhǎng)寬比L/B下進(jìn)行數(shù)值模擬，分析單方柱[3-4]的主要繞流參數(shù)和壓力脈動(dòng)的變化規(guī)律。在此基礎(chǔ)上，取定截面長(zhǎng)寬比,分析雷諾數(shù)以及間距比對(duì)多方柱[5-6]繞流參數(shù)以及壓力脈動(dòng)的影響。

1 數(shù)學(xué)模型

計(jì)算流體力學(xué)（CFD）是流體力學(xué)的重要分支，在計(jì)算機(jī)科學(xué)發(fā)展下產(chǎn)生的介于數(shù)學(xué)、流體力學(xué)和計(jì)算機(jī)之間的交叉學(xué)科，主要研究?jī)?nèi)容是通過(guò)計(jì)算機(jī)和數(shù)值方法來(lái)求解流體力學(xué)的控制方程，對(duì)流體力學(xué)問(wèn)題進(jìn)行模擬和分析，主要遵循質(zhì)量守恒、動(dòng)量守恒和能量守恒三大定律，使用三大基礎(chǔ)方程對(duì)其進(jìn)行描述，文章中未考慮溫度變化，因此不考慮能量守恒方程，剩余兩項(xiàng)的表述如下：

1）連續(xù)性方程

連續(xù)性方程是質(zhì)量守恒定律在流體力學(xué)中具體表達(dá)，對(duì)流體采用連續(xù)介質(zhì)模型，密度和速度在時(shí)間以及空間上都是連續(xù)、可微函數(shù)，數(shù)學(xué)表達(dá)式如下：

2）運(yùn)動(dòng)方程

運(yùn)動(dòng)方程是牛頓第二定律即動(dòng)量守恒定律在流體力學(xué)中的具體表達(dá)，假設(shè)流體不可壓縮時(shí)，數(shù)學(xué)表達(dá)式如下：

本文計(jì)算主要采用SSTk-ω，集成k-ε與k-ω兩種模型的優(yōu)勢(shì)，近壁面處使用k-ω模型，邊界層使用kε模型，考慮湍流剪應(yīng)力的輸運(yùn)效應(yīng)，對(duì)有壓力梯度引起的分離現(xiàn)象模擬更加精確，數(shù)學(xué)表達(dá)式如下：

式中：Gk為層流速度梯度產(chǎn)生的湍流動(dòng)能；Gω由ω方程產(chǎn)生；Γk和Γω是k和ω的擴(kuò)散率；Yk和Yω由擴(kuò)散產(chǎn)生的湍流；Dω為正交發(fā)散項(xiàng)。

2 計(jì)算模型

2.1 計(jì)算域及邊界條件選擇

計(jì)算域大小的選擇至關(guān)重要，從理論上來(lái)說(shuō)計(jì)算域越接近真實(shí)流場(chǎng)其計(jì)算結(jié)果越真實(shí)，但是對(duì)于CFD計(jì)算來(lái)說(shuō)，計(jì)算資源也是需要考慮的因素，受限于電腦的運(yùn)行速度，計(jì)算域越小計(jì)算越快，因此需要在保證計(jì)算結(jié)果的前提下選擇較小的計(jì)算域。參考時(shí)忠民等[7]對(duì)繞流邊界影響的研究，取數(shù)值模擬計(jì)算域的進(jìn)口距方柱取25B，上下兩側(cè)邊界距方柱取25B，出口距方柱取40B，其中方柱長(zhǎng)為L(zhǎng)，寬為B。計(jì)算域邊界左側(cè)為速度入口，右側(cè)為壓力出口，上下兩側(cè)為對(duì)稱邊界，方柱表面取為無(wú)滑移壁面，如圖1所示。

圖1 方柱繞流幾何模型示意圖Fig. 1 Geometric model diagram of flow around rectangular cylinder

2.2 網(wǎng)格獨(dú)立性及數(shù)值方法驗(yàn)證

在數(shù)值模擬計(jì)算中，首先要驗(yàn)證網(wǎng)格獨(dú)立性。網(wǎng)格劃分是前期最重要的工作且耗費(fèi)時(shí)間長(zhǎng)，網(wǎng)格的疏密及質(zhì)量直接影響著計(jì)算規(guī)模的大小和計(jì)算結(jié)果的精確程度。通常認(rèn)為當(dāng)網(wǎng)格密度達(dá)到一定程度后，繼續(xù)增加網(wǎng)格量對(duì)于計(jì)算結(jié)果的影響非常小，此時(shí)可以認(rèn)為網(wǎng)格疏密對(duì)于計(jì)算結(jié)果的影響可以忽略。網(wǎng)格劃分方式時(shí)均采用結(jié)構(gòu)化網(wǎng)格，質(zhì)量均可保證在0.8以上，可忽略質(zhì)量差異對(duì)結(jié)果的影響。圖2給出了本文驗(yàn)證網(wǎng)格獨(dú)立性所采用4套網(wǎng)格。

圖2 不同密度網(wǎng)格示意圖Fig. 2 Different density grid diagram

對(duì)于非穩(wěn)態(tài)計(jì)算，時(shí)間步長(zhǎng)的選取也將影響數(shù)值模擬結(jié)果的準(zhǔn)確性。采用庫(kù)朗數(shù)（CFL）來(lái)調(diào)節(jié)計(jì)算的穩(wěn)定性與收斂性。CFL越大，收斂速度加快但穩(wěn)定性逐漸降低。表1給出雷諾數(shù)Re=5 000時(shí)，不同網(wǎng)格和時(shí)間步長(zhǎng)下方柱繞流的升力系數(shù)的均方根值Cl.rms，阻力系數(shù)的平均值Cd.mean，斯特勞哈爾數(shù)St及各結(jié)果間的相對(duì)誤差。

表1 不同網(wǎng)格和時(shí)間步長(zhǎng)下方柱繞流計(jì)算結(jié)果Tab. 1 Calculation results of flow around square column with different grid and time step

從表1可以看出，網(wǎng)格密度和CFL的變化對(duì)平均阻力系數(shù)Cd.mean和斯特勞哈爾數(shù)St的影響較小，最大相對(duì)誤差僅為3.87%。升力系數(shù)的均方根Cl.rms在網(wǎng)格數(shù)量較低時(shí)計(jì)算結(jié)果偏小，當(dāng)網(wǎng)格數(shù)量到達(dá)11萬(wàn)以上后計(jì)算結(jié)果也趨于穩(wěn)定。因此考慮到計(jì)算消耗，本文劃分網(wǎng)格時(shí)采用114 619網(wǎng)格數(shù)量的劃分方式，時(shí)間步長(zhǎng)上CFL值設(shè)置為5。

針對(duì)所選用的網(wǎng)格以及時(shí)間步長(zhǎng)，并與文獻(xiàn)進(jìn)行對(duì)比，驗(yàn)證數(shù)值方法的正確性，結(jié)果如表2所示。

表2 數(shù)值文獻(xiàn)結(jié)果對(duì)比Tab. 2 Comparison of numerical and literature results

2.3 監(jiān)測(cè)點(diǎn)及多方柱位置布置

探究方柱后流場(chǎng)壓力脈動(dòng)分布規(guī)律，在后方均布25個(gè)監(jiān)測(cè)點(diǎn)，左上角為1號(hào)監(jiān)測(cè)點(diǎn)，從左向右依次編號(hào)，具體情況如圖3所示。多方柱繞流時(shí)，模擬實(shí)船格柵布置情況，以5個(gè)單方柱并列，在保持截面長(zhǎng)寬比L/B=2不變的情況下，改變兩方柱中心距離D的尺寸從而改變方柱間距比D/B，本文取D/B=（5,4,3,2），多方柱時(shí)監(jiān)測(cè)點(diǎn)與單方柱相同，具體尺寸布置及網(wǎng)格情況如圖4和圖5所示。

圖3 監(jiān)測(cè)點(diǎn)分布示意圖Fig. 3 Distribution diagram of monitoring points

圖4 方柱位置布置示意圖Fig. 4 Schematic diagram of the position layout of the rectangular cylinder

圖5 網(wǎng)格示意圖Fig. 5 Grid diagram

3 結(jié)果分析

3.1 單方柱繞流分析

3.1.1 雷諾數(shù)影響分析

給定截面長(zhǎng)寬比L/B為2，不同的雷諾數(shù)Re下繞流參數(shù)（升力系數(shù)的均方根值Cl.rms、阻力系數(shù)的平均值Cd、渦脫頻率F、斯特勞哈爾數(shù)St）如圖6所示，脈動(dòng)壓力變化如圖7所示。

圖6 繞流系數(shù)變化曲線Fig. 6 Curve of variation of flow coefficient

圖7 脈動(dòng)均方根值曲線Fig. 7 Pulsating root-mean-square curve

如圖中所示，升力系數(shù)、阻力系數(shù)，斯特勞哈爾數(shù)變化很小，可以忽略不計(jì)，渦脫頻率隨雷諾數(shù)線性增加。各個(gè)監(jiān)測(cè)點(diǎn)的脈動(dòng)均方根值隨著雷諾數(shù)的增加而增加，總體來(lái)看，不同雷諾數(shù)下，壓力脈動(dòng)均方根值變化趨勢(shì)趨于吻合。從監(jiān)測(cè)點(diǎn)的位置來(lái)看，上面2行監(jiān)測(cè)點(diǎn)均方根值遠(yuǎn)離方柱呈現(xiàn)下降的趨勢(shì)，但變化幅度較小。位于流場(chǎng)中間的監(jiān)測(cè)點(diǎn)壓力脈動(dòng)均方根值遠(yuǎn)離方柱也呈現(xiàn)下降趨勢(shì)，但變化劇烈。下方2行監(jiān)測(cè)點(diǎn)壓力脈動(dòng)均方根值隨距離增加而逐漸增大。

3.1.2 方柱截面長(zhǎng)寬比影響分析

方柱截面的長(zhǎng)寬比會(huì)對(duì)方柱后的流動(dòng)產(chǎn)生顯著影響，本文分別取長(zhǎng)寬比L/B為0.5，1，1.5，2，3的方柱進(jìn)行數(shù)值模擬。各雷諾數(shù)下，不同長(zhǎng)寬比時(shí)繞流參數(shù)如圖8所示。選取Re=15 385時(shí)，針對(duì)不同截面長(zhǎng)寬比，各個(gè)監(jiān)測(cè)點(diǎn)壓力脈動(dòng)均方根值曲線如圖9所示。

由圖可知，不同雷諾數(shù)下，各繞流參數(shù)隨截面長(zhǎng)寬比變化趨勢(shì)相同。當(dāng)L/B=1，2時(shí)，阻力系數(shù)、升力系數(shù)以及斯特勞哈爾數(shù)重合較好；L/B=0.5，1.5，3時(shí)，阻力系數(shù)、升力系數(shù)以及斯特勞哈爾數(shù)在各自雷諾數(shù)下存在差異。阻力系數(shù)曲線與升力系數(shù)曲線先下降再上升，阻力系數(shù)曲線在L/B=1.5取得最小值，升力系數(shù)曲線在L/B=2.0時(shí)取得最小值。渦脫頻率曲線的變化規(guī)律與斯特勞哈爾數(shù)曲線相同，在監(jiān)測(cè)范圍內(nèi)先增大再減小，在L/B=1.5取得最小值時(shí)，逐漸增大。長(zhǎng)寬比L/B=1.5，2時(shí)，壓力脈動(dòng)均方根值明顯小于其他長(zhǎng)寬比下結(jié)果，其中L/B=2時(shí)結(jié)果最小。長(zhǎng)寬比L/B=0.5，1，3時(shí)，壓力脈動(dòng)的均方根值較接近，交替出現(xiàn)最大值。

圖8 繞流系數(shù)變化曲線Fig. 8 Curve of variation of flow coefficient

圖9 脈動(dòng)均方根值曲線Fig. 9 Pulsating root-mean-square curve

3.2 二維方柱群分析

3.2.1 雷諾數(shù)影響分析

方柱的間距比D/B=5，不同雷諾數(shù)下，方柱群繞流參數(shù)如圖10所示，方柱群后壓力脈動(dòng)壓力如圖11所示。

圖10 繞流系數(shù)變化曲線Fig. 10 Curve of variation of flow coefficient

從圖中可以看出，與單方柱相比，柱群的升力系數(shù)與阻力系數(shù)都有明顯的增加，其中位于中間位置的方柱3改變最小，中間兩側(cè)的方柱2與方柱4改變最大。方柱3的渦脫頻率與斯特勞哈爾數(shù)較單方柱有所下降，其余方柱明顯增加。從總體結(jié)果來(lái)看，對(duì)稱布置的方柱1與方柱5，方柱2與方柱4繞流參數(shù)變化趨勢(shì)與幅度接近一致。與單方柱相比，方柱群的壓力脈動(dòng)均方根值均顯著變大，從均方根值變化趨勢(shì)看，除了中間一行監(jiān)測(cè)點(diǎn)外,其余監(jiān)測(cè)點(diǎn)變化趨勢(shì)與單方柱相同。

3.2.2 方柱間距影響分析

雷諾數(shù)為15 385時(shí)，不同間距比時(shí)方柱群繞流參數(shù)如圖12所示，不同方柱間距比下各監(jiān)測(cè)點(diǎn)壓力脈動(dòng)脈動(dòng)壓力如圖13所示。

從圖中可以看出，繞流參數(shù)曲線（升力系數(shù)的均方根值Cl.rms，阻力系數(shù)的平均值Cd，渦脫頻率F、斯特勞哈爾數(shù)St）均具有很好的對(duì)稱性，方柱1和方柱5，方柱2和方柱4主要參數(shù)基本吻合。隨著間距比的增大，升力系數(shù)與阻力系數(shù)呈現(xiàn)減小的趨勢(shì)，方柱3變化最為劇烈。內(nèi)部三方柱渦脫頻率隨間距比的增大而減小，但最外側(cè)兩方柱渦脫頻率隨間距比的增大而增大，當(dāng)D/B=2,3時(shí)，中間方柱與其周圍方柱渦脫頻率差幅較小，隨著間距增大，呈明顯下降趨勢(shì)。斯特勞哈爾數(shù)曲線變化規(guī)律說(shuō)渦脫頻率變化曲線大致相同。方柱間距比對(duì)流場(chǎng)壓力脈動(dòng)的影響很明顯，位于流場(chǎng)前方的監(jiān)測(cè)點(diǎn)壓力脈動(dòng)隨間距比變化混亂，無(wú)規(guī)律可循。流場(chǎng)后方的監(jiān)測(cè)點(diǎn)可以明顯發(fā)現(xiàn)，間距比D/B=2時(shí)壓力脈動(dòng)的均方根值大于其它間距比結(jié)果。

圖11 脈動(dòng)均方根值曲線Fig. 11 Pulsating root-mean-square curve

4 結(jié) 語(yǔ)

本文分析了方柱的主要繞流參數(shù)和壓力脈動(dòng)隨雷諾數(shù)以及方柱截面長(zhǎng)寬比L/B的變化規(guī)律。在此基礎(chǔ)上，取定L/B=2,分析了雷諾數(shù)以及間距比對(duì)繞流參數(shù)以及壓力脈動(dòng)的影響，現(xiàn)得出以下結(jié)論：

單個(gè)二維方柱下，雷諾數(shù)主要影響脈動(dòng)壓力，各監(jiān)測(cè)參數(shù)均隨截面長(zhǎng)寬比發(fā)生較大改變；多方柱并列時(shí)，監(jiān)測(cè)點(diǎn)參數(shù)較單方柱均有所提高，方柱間距比對(duì)脈動(dòng)壓力影響較大。

圖12 繞流系數(shù)變化曲線Fig. 12 Curve of variation of flow coefficient

圖13 脈動(dòng)均方根值曲線Fig. 13 Pulsating root-mean-square curve