許煬塏,梁玉珂,沈林維
(浙江大學(xué)海洋學(xué)院,浙江 舟山 316021)
隨著海洋事業(yè)的發(fā)展,越來(lái)越多的艦船對(duì)低航速甚至零航速下作業(yè)呈現(xiàn)需求的趨勢(shì)。普通減搖鰭因其升力機(jī)制影響無(wú)法在零航速下實(shí)現(xiàn)減搖,而不受船舶航速影響的減搖水艙存在占用船舶內(nèi)部空間和排水量的劣勢(shì)[1],尤其是它的減搖效率存在不確定性。因此,如何減輕低航速下艦船的橫搖運(yùn)動(dòng)是個(gè)亟待解決的問(wèn)題。
2002年,Naiad和MARIN 首次聯(lián)合提出了零航速減搖鰭,并通過(guò)試驗(yàn)發(fā)現(xiàn)該減搖系統(tǒng)的減搖效率達(dá)到63%~95%[2-3]?;谂拇蛞淼挠^察和研究,MARIN[4]提出了一種以弦長(zhǎng)為轉(zhuǎn)軸進(jìn)行拍打的零航速減搖鰭,也稱(chēng)拍打鰭,并在與常規(guī)零航速減搖鰭模型的對(duì)比試驗(yàn)中發(fā)現(xiàn),在恒定轉(zhuǎn)矩下拍打鰭產(chǎn)生的沖量會(huì)高于常規(guī)鰭,所需的瞬時(shí)功率則低于常規(guī)鰭。Gaillarde[5]分析了零航速減搖鰭慣性力和粘性力的作用,以及鰭相對(duì)于船舶橫搖的響應(yīng)時(shí)間,發(fā)現(xiàn)減小鰭的展弦比有利于船舶作業(yè)。國(guó)內(nèi)研究人員通過(guò)分析零航速減搖鰭的展弦比、翼梢形狀等特征,針對(duì)其升力機(jī)制,提出了一種估計(jì)鰭基本尺度的方法[6]。
金鴻章等[7]對(duì)零航速減搖鰭進(jìn)行了受力分析,將升力分解成形狀力、旋渦力和附加質(zhì)量力,發(fā)現(xiàn)鰭的升力與轉(zhuǎn)動(dòng)角速度、角加速度等因素相關(guān)。宋吉廣[8]確定了零航速減搖鰭升力公式中的未知參數(shù),其計(jì)算結(jié)果與試驗(yàn)值接近,為進(jìn)一步研究零航速減搖鰭的升力機(jī)制提供了基礎(chǔ)。葉青云[9]提出了計(jì)算零航速減搖鰭升力的經(jīng)驗(yàn)公式,與CFD 方法得到的模擬值比較符合。Liang 等[10]基于徑向基函數(shù)和回歸神經(jīng)網(wǎng)絡(luò)進(jìn)行了船模橫搖試驗(yàn)和模擬,利用與波浪相關(guān)的相位匹配實(shí)現(xiàn)船舶減橫搖,為零航速減搖鰭減搖方法提供了理論思路。
總之,國(guó)內(nèi)外對(duì)零航速減搖鰭的研究已近20年,零航速減搖鰭幾何模型建立方法已經(jīng)基本成型;但在減搖機(jī)制上大部分是基于實(shí)際應(yīng)用的控制策略研究[11]。本文將利用CFD 方法,在數(shù)值模擬減搖鰭簡(jiǎn)單運(yùn)動(dòng)的基礎(chǔ)上,分析得到減搖鰭運(yùn)動(dòng)與船舶橫搖運(yùn)動(dòng)之間的關(guān)系,最后利用減搖鰭的主動(dòng)簡(jiǎn)諧運(yùn)動(dòng),對(duì)零航速船模在4種規(guī)則波中的橫搖運(yùn)動(dòng)進(jìn)行直接數(shù)值模擬,獲得減搖鰭的減搖效率。
由于工程實(shí)際問(wèn)題所涉及的流動(dòng)復(fù)雜性,考慮到工程精度要求和計(jì)算發(fā)展現(xiàn)狀,流動(dòng)控制方程多會(huì)采用如下雷諾時(shí)均后的連續(xù)性方程和納維-斯托克斯(RANS)方程:
本文采用SSTk-ω湍流模型來(lái)計(jì)算雷諾應(yīng)力。SSTk-ω湍流模型在處理近壁和遠(yuǎn)場(chǎng)流動(dòng)時(shí),分別采用類(lèi)似標(biāo)準(zhǔn)的k-ω和k-ε湍流模型計(jì)算方式,從而較好地實(shí)現(xiàn)由內(nèi)部低雷諾數(shù)到外部高雷諾數(shù)的過(guò)渡[12]。
零航速減搖鰭的運(yùn)動(dòng)可通過(guò)重疊網(wǎng)格技術(shù)實(shí)現(xiàn)。網(wǎng)格由背景區(qū)域和減搖鰭運(yùn)動(dòng)區(qū)域兩部分組成,前者靜止,后者隨鰭作剛性運(yùn)動(dòng),區(qū)域之間進(jìn)行信息傳遞。為確保信息傳遞的連續(xù)和高效,鰭運(yùn)動(dòng)區(qū)域要求始終在背景區(qū)域內(nèi)。本章對(duì)有航速和無(wú)航速兩種情況下的鰭運(yùn)動(dòng)進(jìn)行數(shù)值模擬和驗(yàn)證。
減搖鰭的具體模型及其運(yùn)動(dòng)形式如圖1所示[13]:
圖1 翼型幾何模型運(yùn)動(dòng)形式[13]Fig.1 Motion pattern[13]of airfoil model
均勻來(lái)流速度U=1.5 m/s,鰭圍繞半圓形前緣的圓心進(jìn)行簡(jiǎn)諧運(yùn)動(dòng),旋轉(zhuǎn)角度φ表示為
式中,φ0為鰭旋轉(zhuǎn)的最大角度,sin(φ0)=12A0/11C,f為鰭旋轉(zhuǎn)的頻率。考察鰭模型的弦長(zhǎng)C=1.2×10-2m,斯特勞哈爾數(shù)St=Df/U=0.22,AD=2A0/D=1.81,Re=UC/ν=440,ν為運(yùn)動(dòng)粘度系數(shù),不考慮湍流的影響,采用二維模型。
圖2 為計(jì)算區(qū)域和邊界條件設(shè)置。左邊為速度入口,右邊為壓力出口,頂部、底部及減搖鰭表面均設(shè)置為無(wú)滑移壁面。
圖2 計(jì)算域大小和邊界條件Fig.2 Computational domain and boundary conditions
圖3 為網(wǎng)格設(shè)置示意圖,在鰭的運(yùn)動(dòng)區(qū)域附近設(shè)置兩層網(wǎng)格加密區(qū)。第一層加密網(wǎng)格的尺寸為3.125×10-4m,第二層加密網(wǎng)格的尺寸為6.25×10-4m,最外層背景網(wǎng)格尺寸為1.25×10-3m,總網(wǎng)格數(shù)量為17萬(wàn)。取時(shí)間步長(zhǎng)Δt=2.5×10-5s。
圖3 網(wǎng)格加密情況Fig.3 Mesh encryption
將數(shù)值模擬獲得的升阻力系數(shù)與Andersen等[13]文獻(xiàn)中的試驗(yàn)值進(jìn)行對(duì)比,結(jié)果如圖4所示。二維模擬中,升力系數(shù)CL和阻力系數(shù)Cd的表達(dá)式分別為
圖4 升力系數(shù)和阻力系數(shù)模擬值和試驗(yàn)值的比較Fig.4 Comparison of simulated and tested values of lift and drag coefficients
式中,F(xiàn)L為鰭受到的升力,F(xiàn)d為鰭受到的阻力。
可以發(fā)現(xiàn),升力系數(shù)的模擬值與試驗(yàn)值非常接近,但阻力系數(shù)存在一定差異,這主要有兩方面的原因。一方面,本計(jì)算模型是二維,而試驗(yàn)中采用的是三維模型;另一方面,升力系數(shù)主要由翼型上、下表面壓力差所產(chǎn)生,而阻力系數(shù)受剪切應(yīng)力的影響較大,后者對(duì)流動(dòng)的分離及尾部流動(dòng)更加敏感。
采用茹可夫斯基翼型結(jié)構(gòu)進(jìn)行二維減搖鰭模擬,如圖5所示,其中,d1是鰭軸和弦長(zhǎng)中心之間的距離,C為弦長(zhǎng),S為展長(zhǎng),φ(t)為轉(zhuǎn)動(dòng)角度。本節(jié)考慮鰭弦長(zhǎng)C=2.838 m,展長(zhǎng)S為單位長(zhǎng)度,d1=0.852 m,鰭作簡(jiǎn)諧運(yùn)動(dòng)如公式(3),φ0=π/6,f=1/8。
圖5 翼型幾何模型[8]Fig.5 Geometric model of airfoil
計(jì)算域及網(wǎng)格結(jié)構(gòu)設(shè)置與有航速鰭的類(lèi)似,在鰭運(yùn)動(dòng)區(qū)域周?chē)袃蓪泳W(wǎng)格加密,第一層加密區(qū)的網(wǎng)格尺寸為0.05 m,第二層加密區(qū)的網(wǎng)格尺寸為0.1 m,背景區(qū)域的網(wǎng)格尺寸為0.2 m,總網(wǎng)格數(shù)量合計(jì)為10.6萬(wàn)。
利用尾鰭處最大速度Umax=0.934 m/s 來(lái)計(jì)算升力系數(shù)CL、力矩系數(shù)CM=3Fd/(4ρU2C)和Re,則Re=3×106,因此模擬采用SSTk-ω湍流模型。滿(mǎn)足庫(kù)朗數(shù)CFL=(U·Δt/Δx)<1[14-15],速度變量U取值Umax,Δx為最小網(wǎng)格尺寸0.05 m,計(jì)算合理的時(shí)間步長(zhǎng)范圍Δt<0.54 s,確定時(shí)間步長(zhǎng)為0.05 s。
對(duì)減搖鰭進(jìn)行運(yùn)動(dòng)模擬,并將模擬結(jié)果與Wang 等[16]文獻(xiàn)中的結(jié)果進(jìn)行對(duì)比,如圖6所示。
圖6 模擬值和經(jīng)驗(yàn)值的對(duì)比Fig.6 Comparison of simulated values and empirical values
觀察圖6可知,數(shù)值模擬的水動(dòng)力系數(shù)與Wang等[16]文獻(xiàn)中的模擬值和經(jīng)驗(yàn)值非常接近,力矩系數(shù)則在總體趨勢(shì)上比較一致,但在極值附近存在較大差異。主要原因是當(dāng)減搖鰭在極值位置突然反向運(yùn)動(dòng)時(shí),與周?chē)黧w形成巨大的速度差,產(chǎn)生流動(dòng)分離和旋渦脫落等復(fù)雜現(xiàn)象,從而使鰭的形狀力、旋渦力和附加質(zhì)量力的計(jì)算變得比較困難[16]。另外,高雷諾數(shù)下的湍流模型給計(jì)算帶來(lái)許多不確定因素[17]。
對(duì)某船模的橫搖運(yùn)動(dòng)進(jìn)行模擬,其基本尺度如表1所示。
表1 船模的基本尺度Tab.1 Principal dimensions of ship model
計(jì)算域中,船模的船長(zhǎng)方向與大地坐標(biāo)系的Y軸平行,但仍與其隨體坐標(biāo)系中的X軸平行(如圖7所示)。計(jì)算區(qū)域根據(jù)模型船的船長(zhǎng)來(lái)設(shè)置,分別為-3<(x/Lpp)<1.5、-1<(y/Lpp)和-2<(z/Lpp)<1,坐標(biāo)原點(diǎn)設(shè)置在水面船體中心處。采用重疊網(wǎng)格對(duì)船舶橫搖運(yùn)動(dòng)進(jìn)行模擬,需要設(shè)置一個(gè)包括船體的小長(zhǎng)方體,長(zhǎng)度為1.5倍船長(zhǎng),寬度為0.6倍船長(zhǎng),高度為1倍船長(zhǎng)。
圖7 計(jì)算域及船體運(yùn)動(dòng)區(qū)域Fig.7 Computational domain and ship motion region
計(jì)算域的邊界條件設(shè)置如圖8 所示:入口處對(duì)應(yīng)入射波0.5H·sin(2πt/T1)的速度,其中,H為波高,T1為波浪周期;頂部速度為0;側(cè)部為對(duì)稱(chēng)面;出口處設(shè)為壓力出口,并采用阻尼消波的方式來(lái)防止波浪壁面反射;底部為非滑移壁面。
圖8 橫搖模擬計(jì)算域及邊界條件Fig.8 Computational domain and boundary conditions in ship roll simulation
對(duì)船體近壁面、船體運(yùn)動(dòng)區(qū)域和自由液面附近的流場(chǎng)區(qū)域進(jìn)行網(wǎng)格加密。圖9為在中橫剖面處,未安裝減搖鰭的船舶網(wǎng)格分布圖。背景區(qū)域的網(wǎng)格數(shù)量為70 萬(wàn)左右,船體運(yùn)動(dòng)區(qū)域的網(wǎng)格數(shù)量為117萬(wàn)左右,總網(wǎng)格數(shù)為187萬(wàn)左右。
圖9 中橫剖面處的網(wǎng)格分布圖Fig.9 Grid distribution at the middle section
根據(jù)ITTC[18]對(duì)船舶耐波性計(jì)算案例的建議,本節(jié)對(duì)時(shí)間步長(zhǎng)為0.005 s、0.002 5 s和0.00 125 s的橫搖運(yùn)動(dòng)進(jìn)行獨(dú)立性驗(yàn)證。圖10 顯示的是船舶橫搖穩(wěn)定后,不同時(shí)間步長(zhǎng)下船舶橫搖角隨時(shí)間變化的曲線(xiàn)圖。可以發(fā)現(xiàn),不同時(shí)間步長(zhǎng)下橫搖角、角速度和角加速度的幅值接近,說(shuō)明所選擇的時(shí)間步長(zhǎng)對(duì)結(jié)果影響不大,計(jì)算結(jié)果收斂。綜合考慮計(jì)算效率及波浪模擬的精度,后續(xù)計(jì)算的時(shí)間步長(zhǎng)取為0.002 5 s。
圖10 船舶橫搖運(yùn)動(dòng)時(shí)間步長(zhǎng)獨(dú)立性驗(yàn)證Fig.10 Independent verification of time step in ship rolling motion
根據(jù)Dallinga 和Rapuc[4]的研究,當(dāng)零航速減搖鰭的力矩和船舶橫搖角速度存在相位差π 時(shí),鰭的減搖效果最明顯。國(guó)外研究人員曾做過(guò)有關(guān)鰭-船體適配性研究,試圖探索減搖鰭和作六自由度運(yùn)動(dòng)船舶之間水動(dòng)力的影響,發(fā)現(xiàn)減搖鰭性能的損失主要來(lái)自船體邊界層的影響[9];另外本文主要考慮零航速下的減搖鰭運(yùn)動(dòng),靜止船體對(duì)作簡(jiǎn)諧運(yùn)動(dòng)減搖鰭的水動(dòng)力影響甚小。因此,假設(shè)零航速船體對(duì)減搖鰭運(yùn)動(dòng)水動(dòng)力的影響可以忽略,那么確定鰭運(yùn)動(dòng)相位角的步驟為:對(duì)敞水條件下,作正弦運(yùn)動(dòng)的零航速減搖鰭進(jìn)行數(shù)值模擬,獲得鰭升力與鰭運(yùn)動(dòng)的相位差θ0,如圖11所示;接著,數(shù)值模擬船舶在規(guī)則波作用下的橫搖運(yùn)動(dòng),使得鰭運(yùn)動(dòng)與船舶橫搖角速度之間的相位差為θ0+π。
圖11 減搖鰭轉(zhuǎn)動(dòng)角度與升力系數(shù)之間相位差Fig.11 Phase angle between fin stabilizer turning angle and lift coefficient
將一對(duì)零航速減搖鰭安裝在船中舭部位置,鰭剖面為NACA 0020,弦長(zhǎng)C=0.2 m,展長(zhǎng)S=0.1 m,展弦比為0.5,鰭軸設(shè)置在距鰭前緣1/4位置處。減搖鰭的幾何模型如圖12所示:
圖12 鰭的幾何模型Fig.12 Geometric model of fin
用無(wú)量綱的面積系數(shù)Cs來(lái)預(yù)估鰭和船體之間的幾何大小關(guān)系,
式中,A1為安裝在船體上所有減搖鰭的最大面積C×S之和,Lpp為船長(zhǎng),B為船寬。綦志剛[19]、Huang[20]和Ghassemi 等[21]的Cs分別為0.012、0.014 和0.01,本文采用的減搖鰭Cs=0.016,與其他文獻(xiàn)的結(jié)果比較接近。根據(jù)船模對(duì)應(yīng)實(shí)船航行時(shí)的要求,本計(jì)算取5 級(jí)海況下對(duì)應(yīng)的模擬波高為0.08 m,并在4 個(gè)不同周期的規(guī)則波下進(jìn)行減搖效率驗(yàn)證。
圖13 為船舶橫搖穩(wěn)定后,CFD 模擬得到的無(wú)減搖鰭和安裝減搖鰭的零航速船舶橫搖角變化曲線(xiàn)比較。定義減搖效率η為
圖13 無(wú)減搖鰭和安裝減搖鰭的零航速船舶橫搖角曲線(xiàn)Fig.13 Roll angle of ship at zero speed with and without fin stabilizer
式中,φ0和φ1分別為船舶橫搖穩(wěn)定時(shí),未安裝減搖鰭和安裝零航速減搖鰭的船舶最大橫搖角。則4個(gè)入射波條件下對(duì)應(yīng)的鰭的減搖效率分別為61%、62%、63%和81%。
本文利用CFD軟件,在數(shù)值模擬減搖鰭簡(jiǎn)單運(yùn)動(dòng)的基礎(chǔ)上,分析得到減搖鰭運(yùn)動(dòng)與船舶橫搖運(yùn)動(dòng)之間的關(guān)系,最后利用減搖鰭的主動(dòng)簡(jiǎn)諧運(yùn)動(dòng),對(duì)零航速船模在規(guī)則波中的橫搖運(yùn)動(dòng)進(jìn)行直接數(shù)值模擬,得出如下結(jié)論:
(1)在敞水條件下,對(duì)減搖鰭在有、無(wú)航速兩種狀態(tài)下作簡(jiǎn)諧運(yùn)動(dòng)時(shí)所受的流體動(dòng)力進(jìn)行數(shù)值計(jì)算,得到的結(jié)果與文獻(xiàn)中的試驗(yàn)和數(shù)值計(jì)算結(jié)果比較吻合,為在復(fù)雜運(yùn)動(dòng)模式下的流體動(dòng)力計(jì)算做好準(zhǔn)備。
(2)利用零航速減搖鰭運(yùn)動(dòng)與船舶橫搖運(yùn)動(dòng)之間的關(guān)系,對(duì)零航速船模在相同波高、4 個(gè)不同周期規(guī)則波下的橫搖運(yùn)動(dòng)進(jìn)行數(shù)值預(yù)報(bào),發(fā)現(xiàn)減搖鰭的減搖效率在61%和81%之間,為進(jìn)一步研究零航速船舶在不規(guī)則波下的減搖運(yùn)動(dòng)打下基礎(chǔ)。