李 超,宋 健,劉金勇,艾艷輝,佘湖清
(中國船舶集團(tuán)有限公司第七一〇研究所,湖北 宜昌 443000)
半潛航行器作為一種近水面航行器,通過自身舵板調(diào)節(jié)姿態(tài),借助舵板和主體產(chǎn)生定深力控制航行深度,降低波浪干擾,提高海況適應(yīng)性。研究航行器流體動力性能并求解相關(guān)系數(shù),有助于預(yù)報半潛航行器運(yùn)動特性,掌握控制規(guī)律。
針對半潛航行器流體動力性能研究方法主要有模型試驗、數(shù)值計算、半經(jīng)驗估算等。Mori 等[1]提出一種帶翼板的高速半潛式航行器,借助翼板產(chǎn)生的向下升力增加潛深,翼板也有助于消減行波,降低興波阻力,之后Mori 等[2]采用數(shù)值模擬方法研究了半潛航行器的運(yùn)動特性,分析了不同航行深度、翼板位置和攻角下的流場。談果戈[3]通過CFD 方法模擬平面運(yùn)動機(jī)構(gòu),計算了不同航行深度下半潛航行器慣性類和速度水動力系數(shù)。Ueno 等[4]基于Hirayama[5]提出的半潛航行器模型通過理論和半經(jīng)驗公式估算了航行器線性流體動力系數(shù),闡明了該航行器在垂直面和水平面運(yùn)動中的動力特性,同時指出通過試驗確定流體動力系數(shù)的必要性。后續(xù)通過水池試驗,確定了不同潛深下速度和舵角流體動力系數(shù),并指出通過估算的方法對舵某些系數(shù)的準(zhǔn)確度有待提高[6]。為了驗證半潛航行器多狀態(tài)轉(zhuǎn)換的可行性,Cong[7]對帶有壓載水艙的半潛式航行器開展了自航試驗,研究了不同前后翼角組合對航行器狀態(tài)轉(zhuǎn)變的影響,測量了航行過程中翼和主體的阻力、升力和縱搖力矩的變化。Cong 指出航行器運(yùn)動的非線性、非線性水動力和剩余浮力等因素增加了狀態(tài)轉(zhuǎn)換的難度,準(zhǔn)確測量機(jī)翼上的流體動力有助于更好地理解多狀態(tài)轉(zhuǎn)換過程的機(jī)理。對于該種類型的半潛式航行器,董文才等[8]也做過相關(guān)研究。大多數(shù)研究集中在半潛航行器的速度系數(shù)、舵角系數(shù)、加速度系數(shù)的求解,并對半潛航行器航態(tài)轉(zhuǎn)換性能進(jìn)行了分析,對關(guān)系到半潛航行器大幅度運(yùn)動性能的速度和舵角的耦合系數(shù)研究較少。
針對目前研究的不足,通過數(shù)值方法對半潛航行器的縱向速度和舵角系數(shù)進(jìn)行仿真計算,并與風(fēng)洞試驗結(jié)果進(jìn)行對比驗證,通過最小二乘法回歸了速度、舵角耦合系數(shù)。對于耦合系數(shù)中舵效降低的情況結(jié)合局部流場進(jìn)行了深入分析,確定了適用于該航行器的最佳舵角范圍。
圖1 所示的半潛航行器由回轉(zhuǎn)體、桅桿、下附體、前水平舵、后水平舵及后垂直舵組成。前舵左右對稱分布,后舵為十字分布。航行器長L為6.8 m,回轉(zhuǎn)體部分直徑 φ為0.92 m,總高度H為3.6 m,寬度B為1.84 m,前后舵板采用NACA0012 翼型。計算域劃分和邊界條件設(shè)置見圖2,流體域入口到首部為1 倍航行器長度,出口到尾部為3 倍航行器長度,左右到中垂線各2 倍航行器長度。圖中①號面設(shè)為對稱平面,②號面為壓力出口,其余面為速度入口。
圖1 航行器外型Fig.1 Vehicle appearance
圖2 流體域劃分Fig.2 Division of fluid domain
數(shù)值模擬主要通過求解RANS 方程得到仿真結(jié)果,其控制方程如下:
Reynolds 應(yīng)力實際上有6 個不同的應(yīng)力項,為了構(gòu)成封閉方程組,需要引入其他條件對方程進(jìn)行求解。剪切應(yīng)力輸運(yùn)k-ω模型[9](簡稱SSTk-ω模型)在計算航行器繞流方面應(yīng)用廣泛[10-12]。其方程表示如下:
式中:
式中:F1,F(xiàn)2為混合函數(shù)。
式中:Gk,Gω分別為湍動能和 ω方程;Yk,Yω分別為k和k-ω的發(fā)散項;Gω為正交發(fā)散項。
半潛航行器航行深度超過0.35L,該深度下水面文丘里效應(yīng)對航行器影響不再明顯[13],航行器通過流線型桅桿與水面接觸,興波對升阻力影響暫不考慮,因此在計算流體動力系數(shù)的過程中不考慮自由液面。影響航行器計算精度的要素主要有湍流模型、網(wǎng)格密度、時間步長等,其中網(wǎng)格密度是網(wǎng)格無關(guān)性問題談?wù)摰闹攸c[14-15]。本文航行器基礎(chǔ)網(wǎng)格分別為0.8 m,0.65 m,0.5 m,如圖3 和圖4 所示。為了更好地契合舵面物理形狀,需對舵板表面、端面進(jìn)行網(wǎng)格加密,網(wǎng)格大小為航行器表面網(wǎng)格的一半,航行器周圍整體網(wǎng)格進(jìn)行多層加密,網(wǎng)格大小為表面網(wǎng)格的2 倍、4 倍、8 倍。按此規(guī)律得到網(wǎng)格分別為258 萬、325 萬、384 萬,3 種網(wǎng)格能夠貼合舵面形狀,繼續(xù)加密會消耗更多的計算資源,使得計算周期變長。
圖3 流體域網(wǎng)格Fig.3 Grid of fluid domain
圖4 不同表面網(wǎng)格密度Fig.4 Different surface mesh density
對正浮態(tài)下不同網(wǎng)格密度航行器進(jìn)行仿真計算,得到不同航速下升、阻力結(jié)果,如表1 所示。阻力和升力隨網(wǎng)格數(shù)量的增加有變小的趨勢,以最密網(wǎng)格為基準(zhǔn),剩余2 種網(wǎng)格阻力結(jié)果與之平均誤差分別為2.2%和5.15%,升力結(jié)果平均誤差2.1%和5.25%,低速情況下因基礎(chǔ)數(shù)值較小,誤差較大。綜合考慮計算成本與結(jié)果準(zhǔn)確度,網(wǎng)格基礎(chǔ)為0.65 m。
表1 主要參數(shù)Tab.1 Primary parameters
試驗?zāi)P椭黧w為玻璃鋼,內(nèi)設(shè)鋼骨架,表面噴漆。風(fēng)洞模型試驗在西北工業(yè)大學(xué)風(fēng)洞試驗室進(jìn)行。采用PSI8400 電子掃描閥數(shù)據(jù)系統(tǒng)采集天平輸出的電壓信號,采集速度為5 萬點/秒,測量精度為0.1%。模型支撐方式為單支桿腹部支撐,為減少撐桿對試驗結(jié)果的影響,撐桿進(jìn)行導(dǎo)流處理,測力天平采用內(nèi)測式六分量應(yīng)變天平。模型縮比為1∶2,試驗平均氣溫10℃。
圖5 試驗?zāi)P虵ig.5 Test model
由于風(fēng)洞試驗雷諾數(shù)小于航行器實際雷諾數(shù),首先驗證雷諾數(shù)對試驗結(jié)果的影響,風(fēng)速40 m/s,50 m/s,60 m/s;然后開展航行器速度和舵角系數(shù)試驗,固定航速,通過改變航行器攻角來改變其垂向速度,攻角變化范圍為-9°~9°,中間間隔3°,保持航行器正浮態(tài),依次改變前、后舵角,舵角變化范圍-15°~15°。為了消除腹部支撐桿對結(jié)果的影響,還需進(jìn)行支撐桿干擾消除試驗。
航行器最低航速下雷諾數(shù)Re=1.6×107,該雷諾數(shù)下需要風(fēng)洞風(fēng)速為76 m/s,現(xiàn)有試驗風(fēng)速最大到60 m/s,該風(fēng)速下模型雷諾數(shù)超過了阻力系數(shù)變化較大的過渡區(qū),但仍無法滿足雷諾相似,因此需要研究雷諾數(shù)對結(jié)果的影響規(guī)律。舵角保持0°,風(fēng)速60 m/s,50 m/s,40 m/s,30 m/s,航行器攻角α 角度范圍-9°~9°。圖6中曲線斜率即為航行器的流體動力系數(shù),可以看出雷諾數(shù)對阻力系數(shù)曲線影響較大,對升力系數(shù)和縱搖力矩系數(shù)影響較小。風(fēng)速從60m/s 到30m/s,對應(yīng)的阻力系數(shù)分別為1.08×10-4,1.07×10-4,1.04×10-4,9.8×10-5,正浮狀態(tài)下阻力系數(shù)CX0分別為-0.018 8,-0.018 5,-0.017 9,-0.017 3。根據(jù)結(jié)果趨勢,76 m/s 風(fēng)速下各力和力矩系數(shù)可采用60 m/s 風(fēng)速的結(jié)果,初始阻力系數(shù)需按式(10)進(jìn)行擬合修正。
圖6 不同風(fēng)速結(jié)果曲線Fig.6 Result curves of different wind speeds
本文水動力系數(shù)采用潛艇方程表達(dá)形式[13],將試驗數(shù)據(jù)按式(11)~式(18)方式進(jìn)行處理。X,Z,M分別表示局部坐標(biāo)下的阻力、升力和縱搖力矩,按右手坐標(biāo)系建立局部坐標(biāo)o-xyz,x軸指向航行器首部,z軸豎直向下。下標(biāo)m和s分別表示模型數(shù)據(jù)和仿真數(shù)據(jù),*efd為風(fēng)洞試驗結(jié)果,*EFD為試驗換算后結(jié)果,為換算后的無因次化結(jié)果。
換算公式如下:
3.3.1 速度系數(shù)
航行器初始航速10 kn,對應(yīng)的不同攻角下縱向速度分別為-0.805 m/s,-0.538 m/s,-0.269 m/s,0,0.269 m/s,0.538 m/s,0.805 m/s,正、負(fù)號表示方向,航行器抬首,攻角為正,垂向速度為正。圖(7)給出了局部坐標(biāo)系下航行器阻力、升力、縱搖力矩CFD 仿真和風(fēng)洞試驗曲線。由圖7(a)可以看出,阻力隨航行器垂向速度變化呈較強(qiáng)的偶函數(shù)性質(zhì),垂向速度較小時,CFD 結(jié)果和試驗結(jié)果誤差較小,垂向速度-0.805 m/s時,二者誤差達(dá)到9%。根據(jù)函數(shù)性質(zhì)和式(14)求得無因次化后的流體動力系數(shù)和,具體結(jié)果如表2 所示。
圖7(b)與圖7(c)中CFD 仿真結(jié)果和風(fēng)洞試驗結(jié)果誤差較小,最大誤差分別為4.9%和6.3%。升力和縱搖力矩系數(shù)隨航行器垂向速度變化呈奇函數(shù)特征,半潛航行器上下結(jié)構(gòu)不對稱導(dǎo)致曲線未呈現(xiàn)完全對稱,抬首產(chǎn)生的升力和力矩都略大于埋首。垂向速度大于0.269m/s(對應(yīng)航行器抬首或埋首3°以上),升力和縱搖力矩呈現(xiàn)非線性。航行器抬首產(chǎn)生會抬首力矩,埋首自身產(chǎn)生埋首力矩,表明航行器本身的流體動力并未提供相應(yīng)的回復(fù)力矩。
圖7 速度系數(shù)曲線Fig.7 Velocity coefficient curve
根據(jù)函數(shù)性質(zhì)和式(15)~式(18)求得無因次化后的流體動力系數(shù)具體結(jié)果見表2。
表2 速度系數(shù)Tab.2 Velocity coefficient
3.3.2 舵角系數(shù)
半潛航行器由前水平舵、后水平舵、后垂直舵組成,航行過程中通過前、后水平舵保持航行器的航行深度和調(diào)節(jié)運(yùn)動姿態(tài)。前、后舵角變化范圍從-15°~15°,間隔3°,對應(yīng)的舵角從-0.209 rad 到0.209 rad,正負(fù)號代表方向,負(fù)號表示打下舵角。
由圖8 可以看出,前、后舵舵角變化產(chǎn)生的阻力曲線接近偶函數(shù)性質(zhì),前舵最大誤差3.8%,后舵最大誤差6.1%,求得無因次舵角系數(shù),如表3 所示。
表3 舵角流體動力系數(shù)Tab.3 Hydrodynamic coefficient of rudder angle
圖8 舵阻力系數(shù)曲線Fig.8 Rudder resistance coefficient curve
如圖9 和圖10 所示,前舵升力和力矩具有良好的線性,最大舵角15°時未發(fā)生舵效降低的現(xiàn)象或趨勢。后舵的舵角9°以內(nèi)變化趨勢與前舵相近,升力數(shù)值大于前舵,打上舵角9°時開始出現(xiàn)舵效降低的現(xiàn)象,繼續(xù)增加舵角升力增加不明顯,該種現(xiàn)象主要原因在于前舵稍渦對后舵來流產(chǎn)生了影響,導(dǎo)致后舵失速角變小。對縱搖力矩起主要作用的是升力,后舵升力比前舵要大,但大舵下的縱搖力矩要小于前舵,表明在該種情況下后舵壓力中心位置發(fā)生了偏移。在實際操縱過程中前舵在-15°~15°舵角范圍內(nèi)進(jìn)行操舵,后舵在-9°~9°舵角范圍內(nèi)進(jìn)行操舵,得到前、后舵無因次流體動力系數(shù)見表3。
圖9 舵升力系數(shù)曲線Fig.9 Rudder lift coefficient curve
圖10 舵縱搖力矩系數(shù)曲線Fig.10 Rudder pitching moment coefficient curve
3.3.3 耦合系數(shù)
航行器受波浪影響會在埋首或抬首姿態(tài)情況下進(jìn)行操舵,舵周期流場會受到航行器干擾,該種干擾表現(xiàn)為縱向耦合水動力。耦合水動力受航行器垂向速度和舵角的共同影響,升力可表示為:
通過改變航行器垂向速度和舵角,得到多組仿真數(shù)據(jù),可表示為:
ZE的最小二估計量為(X′X)-1X′Z,Z的最小二乘估計量為X(X′X)-1X′Z,力矩的處理方式和力相同,可得到耦合系數(shù)Zwδ和Mwδ。圖11 為升力和縱搖力矩擬合曲面,無因次化系數(shù)和如表4 所示。
表4 耦合流體動力系數(shù)Tab.4 Coupled hydrodynamic coefficient
圖11 速度舵角升力擬合曲面Fig.11 Fitting surface of velocity and lift generated rudder
可知,垂向速度的變化范圍為-0.538~0.538 m/s,前、后舵角變化范圍-0.157~0.157rad。耦合曲面與仿真數(shù)據(jù)擬合較好,滿足線性擬合趨勢。航行器垂向速度0.538 m/s,舵角0.157rad,前舵和后舵產(chǎn)生的升力皆降低,該工況對應(yīng)的航行器攻角為6°,舵角為9°。上文前舵15°舵角未出現(xiàn)舵面升力降低的情況,表明前舵升力降低的主要原因為航行器攻角改變了周圍流場的分布,對前舵產(chǎn)生影響。后舵升力降低的原因除了受前舵影響外可能和航行器攻角也有關(guān),需結(jié)合流場進(jìn)行進(jìn)一步分析。
如圖12 所示,縱搖力矩擬合曲面與仿真數(shù)據(jù)擬合較好。與升力耦合曲面相似,航行器攻角為6°,舵角為9°,縱搖力矩數(shù)據(jù)點未落到曲面上,且誤差較大,表明此時舵的調(diào)節(jié)效果降低,開始出現(xiàn)非線性變化規(guī)律,增加舵角甚至還會起到反作用。
圖12 速度舵角縱搖力矩擬合曲面Fig.12 Fitting surface of velocity and pitch moment about rudder
根據(jù)仿真結(jié)果,前、后舵在航行器有攻角情況下更易發(fā)生舵效降低的情況。為了研究該種現(xiàn)象產(chǎn)生的原因,對前、后舵周圍流場進(jìn)行分析。圖13 為航行器攻角0°和6°,前、后舵9°情況下流線分布圖。為了更好描述,將靠近航行器主體的舵面稱為貼體端,另一面為自由端,帶·的圖航行器攻角為0°。
圖13 前、后舵流線Fig.13 Streamlines of front and rear rudder
如圖13(a)和圖13(b)所示,自由端稍渦由壓力面和吸力面兩部分渦系組成,在舵板后緣匯聚。舵板前緣流線能夠快速成型并向下遷移,靠近中后緣的流線會被稍渦捕獲。航行器攻角增加了吸力面流線和稍渦脫落速度。
如圖13(c)和圖13(d)所示,貼體端稍渦主要由吸力面邊界層分離產(chǎn)生,航行器的存在阻礙了壓力面稍渦的形成,此時流體在壓力的作用下流向上表面。隨著航行器攻角的增加,前、后舵在舵面1/4 以后位置形成直徑與舵板弦長相當(dāng)?shù)谋砻鏈u,且后舵的渦更劇烈,這也是造成升力損失的主要原因。
結(jié)合前、后舵板表面極限流線進(jìn)行進(jìn)一步分析,舵面較小的是前舵,較大的是后舵。如圖14 所示,同一片舵中左側(cè)為舵板壓力面,右側(cè)為舵板吸力面。兩端稍渦使得極限流線發(fā)生偏轉(zhuǎn),且吸力面偏轉(zhuǎn)更嚴(yán)重,表明吸力面受稍渦的影響更明顯。由圖15 可以看出,在兩端稍渦的影響下,極限流線偏轉(zhuǎn)更嚴(yán)重,在舵板吸力面中間形成表面渦,表面渦與兩端稍渦共同作用增加了吸力面壓力。表明航行器攻角改變了周圍流場環(huán)境,使得舵板失速角變小。
圖14 航行器0°攻角舵板極限流線Fig.14 Limit streamline of rudder plate at 0° attack angle of vehicle
圖15 航行器6°攻角舵板極限流線Fig.15 Limit streamline of rudder plate at 6° attack angle of vehicle
針對半潛航行器縱向水動力性能預(yù)報問題,通過數(shù)值方法開展了半潛航行器的縱向速度和舵角系數(shù)仿真計算,并與風(fēng)洞試驗結(jié)果進(jìn)行對比驗證;通過最小二乘法回歸了速度、舵角耦合系數(shù),分析了影響舵效的主要原因。得出主要結(jié)論如下:
1)CFD 仿真結(jié)果與試驗結(jié)果吻合較好,升力和縱搖力矩最大誤差不超過6.5%,可滿足工程應(yīng)用的需要;
2)航行器無攻角情況下,后舵會受前舵影響,導(dǎo)致有效舵角范圍變小,可通過增加前、后舵垂直方向距離進(jìn)一步優(yōu)化性能;
3)增加航行器攻角會影響前后舵板貼體端稍渦的形成,迫使稍渦向上偏移,兩端稍渦迫使舵面流線發(fā)生偏移,在吸力面形成一個直徑與舵板弦長相當(dāng)?shù)臏u。