半潛航行器縱平面流體動力數值計算與試驗研究

2023-05-04 14:01:36劉金勇艾艷輝佘湖清

艦船科學技術 2023年5期

李超，宋健，劉金勇，艾艷輝，佘湖清

(中國船舶集團有限公司第七一〇研究所，湖北宜昌 443000)

0 引言

半潛航行器作為一種近水面航行器，通過自身舵板調節姿態，借助舵板和主體產生定深力控制航行深度，降低波浪干擾，提高海況適應性。研究航行器流體動力性能并求解相關系數，有助于預報半潛航行器運動特性，掌握控制規律。

針對半潛航行器流體動力性能研究方法主要有模型試驗、數值計算、半經驗估算等。Mori 等[1]提出一種帶翼板的高速半潛式航行器，借助翼板產生的向下升力增加潛深，翼板也有助于消減行波，降低興波阻力，之后Mori 等[2]采用數值模擬方法研究了半潛航行器的運動特性，分析了不同航行深度、翼板位置和攻角下的流場。談果戈[3]通過CFD 方法模擬平面運動機構，計算了不同航行深度下半潛航行器慣性類和速度水動力系數。Ueno 等[4]基于Hirayama[5]提出的半潛航行器模型通過理論和半經驗公式估算了航行器線性流體動力系數，闡明了該航行器在垂直面和水平面運動中的動力特性，同時指出通過試驗確定流體動力系數的必要性。后續通過水池試驗，確定了不同潛深下速度和舵角流體動力系數，并指出通過估算的方法對舵某些系數的準確度有待提高[6]。為了驗證半潛航行器多狀態轉換的可行性，Cong[7]對帶有壓載水艙的半潛式航行器開展了自航試驗，研究了不同前后翼角組合對航行器狀態轉變的影響，測量了航行過程中翼和主體的阻力、升力和縱搖力矩的變化。Cong 指出航行器運動的非線性、非線性水動力和剩余浮力等因素增加了狀態轉換的難度，準確測量機翼上的流體動力有助于更好地理解多狀態轉換過程的機理。對于該種類型的半潛式航行器，董文才等[8]也做過相關研究。大多數研究集中在半潛航行器的速度系數、舵角系數、加速度系數的求解，并對半潛航行器航態轉換性能進行了分析，對關系到半潛航行器大幅度運動性能的速度和舵角的耦合系數研究較少。

針對目前研究的不足，通過數值方法對半潛航行器的縱向速度和舵角系數進行仿真計算，并與風洞試驗結果進行對比驗證，通過最小二乘法回歸了速度、舵角耦合系數。對于耦合系數中舵效降低的情況結合局部流場進行了深入分析，確定了適用于該航行器的最佳舵角范圍。

1 數值計算模型及空間離散

1.1 幾何模型

圖1 所示的半潛航行器由回轉體、桅桿、下附體、前水平舵、后水平舵及后垂直舵組成。前舵左右對稱分布，后舵為十字分布。航行器長L為6.8 m，回轉體部分直徑 φ為0.92 m，總高度H為3.6 m，寬度B為1.84 m，前后舵板采用NACA0012 翼型。計算域劃分和邊界條件設置見圖2，流體域入口到首部為1 倍航行器長度，出口到尾部為3 倍航行器長度，左右到中垂線各2 倍航行器長度。圖中①號面設為對稱平面，②號面為壓力出口，其余面為速度入口。

圖1 航行器外型Fig.1 Vehicle appearance

圖2 流體域劃分Fig.2 Division of fluid domain

1.2 數值方法

數值模擬主要通過求解RANS 方程得到仿真結果，其控制方程如下：

Reynolds 應力實際上有6 個不同的應力項，為了構成封閉方程組，需要引入其他條件對方程進行求解。剪切應力輸運k-ω模型[9](簡稱SSTk-ω模型)在計算航行器繞流方面應用廣泛[10-12]。其方程表示如下：

式中：

式中：F1，F2為混合函數。

式中：Gk，Gω分別為湍動能和 ω方程；Yk，Yω分別為k和k-ω的發散項；Gω為正交發散項。

1.3 網格劃分及無關性討論

半潛航行器航行深度超過0.35L，該深度下水面文丘里效應對航行器影響不再明顯[13]，航行器通過流線型桅桿與水面接觸，興波對升阻力影響暫不考慮，因此在計算流體動力系數的過程中不考慮自由液面。影響航行器計算精度的要素主要有湍流模型、網格密度、時間步長等，其中網格密度是網格無關性問題談論的重點[14-15]。本文航行器基礎網格分別為0.8 m，0.65 m，0.5 m，如圖3 和圖4 所示。為了更好地契合舵面物理形狀，需對舵板表面、端面進行網格加密，網格大小為航行器表面網格的一半，航行器周圍整體網格進行多層加密，網格大小為表面網格的2 倍、4 倍、8 倍。按此規律得到網格分別為258 萬、325 萬、384 萬，3 種網格能夠貼合舵面形狀，繼續加密會消耗更多的計算資源，使得計算周期變長。

圖3 流體域網格Fig.3 Grid of fluid domain

圖4 不同表面網格密度Fig.4 Different surface mesh density

對正浮態下不同網格密度航行器進行仿真計算，得到不同航速下升、阻力結果，如表1 所示。阻力和升力隨網格數量的增加有變小的趨勢，以最密網格為基準，剩余2 種網格阻力結果與之平均誤差分別為2.2%和5.15%，升力結果平均誤差2.1%和5.25%，低速情況下因基礎數值較小，誤差較大。綜合考慮計算成本與結果準確度，網格基礎為0.65 m。

表1 主要參數Tab.1 Primary parameters

2 風洞試驗模型及內容

2.1 網格劃分及無關性討論

試驗模型主體為玻璃鋼，內設鋼骨架，表面噴漆。風洞模型試驗在西北工業大學風洞試驗室進行。采用PSI8400 電子掃描閥數據系統采集天平輸出的電壓信號，采集速度為5 萬點/秒，測量精度為0.1%。模型支撐方式為單支桿腹部支撐，為減少撐桿對試驗結果的影響，撐桿進行導流處理，測力天平采用內測式六分量應變天平。模型縮比為1∶2，試驗平均氣溫10℃。

圖5 試驗模型Fig.5 Test model

2.2 試驗內容

由于風洞試驗雷諾數小于航行器實際雷諾數，首先驗證雷諾數對試驗結果的影響，風速40 m/s，50 m/s，60 m/s；然后開展航行器速度和舵角系數試驗，固定航速，通過改變航行器攻角來改變其垂向速度，攻角變化范圍為-9°～9°，中間間隔3°，保持航行器正浮態，依次改變前、后舵角，舵角變化范圍-15°～15°。為了消除腹部支撐桿對結果的影響，還需進行支撐桿干擾消除試驗。

3 結果分析

3.1 雷諾數影響分析

航行器最低航速下雷諾數Re=1.6×107，該雷諾數下需要風洞風速為76 m/s，現有試驗風速最大到60 m/s，該風速下模型雷諾數超過了阻力系數變化較大的過渡區，但仍無法滿足雷諾相似，因此需要研究雷諾數對結果的影響規律。舵角保持0°，風速60 m/s，50 m/s，40 m/s，30 m/s，航行器攻角α 角度范圍-9°～9°。圖6中曲線斜率即為航行器的流體動力系數，可以看出雷諾數對阻力系數曲線影響較大，對升力系數和縱搖力矩系數影響較小。風速從60m/s 到30m/s，對應的阻力系數分別為1.08×10-4，1.07×10-4，1.04×10-4，9.8×10-5，正浮狀態下阻力系數CX0分別為-0.018 8，-0.018 5，-0.017 9，-0.017 3。根據結果趨勢，76 m/s 風速下各力和力矩系數可采用60 m/s 風速的結果，初始阻力系數需按式(10)進行擬合修正。

圖6 不同風速結果曲線Fig.6 Result curves of different wind speeds

3.2 結果換算

本文水動力系數采用潛艇方程表達形式[13]，將試驗數據按式(11)～式(18)方式進行處理。X，Z，M分別表示局部坐標下的阻力、升力和縱搖力矩，按右手坐標系建立局部坐標o-xyz，x軸指向航行器首部，z軸豎直向下。下標m和s分別表示模型數據和仿真數據，*efd為風洞試驗結果，*EFD為試驗換算后結果，為換算后的無因次化結果。

換算公式如下：

3.3 結果對比分析

3.3.1 速度系數

航行器初始航速10 kn，對應的不同攻角下縱向速度分別為-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，正、負號表示方向，航行器抬首，攻角為正，垂向速度為正。圖(7)給出了局部坐標系下航行器阻力、升力、縱搖力矩CFD 仿真和風洞試驗曲線。由圖7(a)可以看出，阻力隨航行器垂向速度變化呈較強的偶函數性質，垂向速度較小時，CFD 結果和試驗結果誤差較小，垂向速度-0.805 m/s時，二者誤差達到9%。根據函數性質和式(14)求得無因次化后的流體動力系數和，具體結果如表2 所示。

圖7(b)與圖7(c)中CFD 仿真結果和風洞試驗結果誤差較小，最大誤差分別為4.9%和6.3%。升力和縱搖力矩系數隨航行器垂向速度變化呈奇函數特征，半潛航行器上下結構不對稱導致曲線未呈現完全對稱，抬首產生的升力和力矩都略大于埋首。垂向速度大于0.269m/s(對應航行器抬首或埋首3°以上)，升力和縱搖力矩呈現非線性。航行器抬首產生會抬首力矩，埋首自身產生埋首力矩，表明航行器本身的流體動力并未提供相應的回復力矩。

圖7 速度系數曲線Fig.7 Velocity coefficient curve

根據函數性質和式(15)～式(18)求得無因次化后的流體動力系數具體結果見表2。

表2 速度系數Tab.2 Velocity coefficient

3.3.2 舵角系數

半潛航行器由前水平舵、后水平舵、后垂直舵組成，航行過程中通過前、后水平舵保持航行器的航行深度和調節運動姿態。前、后舵角變化范圍從-15°～15°，間隔3°，對應的舵角從-0.209 rad 到0.209 rad，正負號代表方向，負號表示打下舵角。

由圖8 可以看出，前、后舵舵角變化產生的阻力曲線接近偶函數性質，前舵最大誤差3.8%，后舵最大誤差6.1%，求得無因次舵角系數，如表3 所示。

表3 舵角流體動力系數Tab.3 Hydrodynamic coefficient of rudder angle

圖8 舵阻力系數曲線Fig.8 Rudder resistance coefficient curve

如圖9 和圖10 所示，前舵升力和力矩具有良好的線性，最大舵角15°時未發生舵效降低的現象或趨勢。后舵的舵角9°以內變化趨勢與前舵相近，升力數值大于前舵，打上舵角9°時開始出現舵效降低的現象，繼續增加舵角升力增加不明顯，該種現象主要原因在于前舵稍渦對后舵來流產生了影響，導致后舵失速角變小。對縱搖力矩起主要作用的是升力，后舵升力比前舵要大，但大舵下的縱搖力矩要小于前舵，表明在該種情況下后舵壓力中心位置發生了偏移。在實際操縱過程中前舵在-15°～15°舵角范圍內進行操舵，后舵在-9°～9°舵角范圍內進行操舵，得到前、后舵無因次流體動力系數見表3。

圖9 舵升力系數曲線Fig.9 Rudder lift coefficient curve

圖10 舵縱搖力矩系數曲線Fig.10 Rudder pitching moment coefficient curve

3.3.3 耦合系數

航行器受波浪影響會在埋首或抬首姿態情況下進行操舵，舵周期流場會受到航行器干擾，該種干擾表現為縱向耦合水動力。耦合水動力受航行器垂向速度和舵角的共同影響，升力可表示為：

通過改變航行器垂向速度和舵角，得到多組仿真數據，可表示為：

ZE的最小二估計量為(X′X)-1X′Z，Z的最小二乘估計量為X(X′X)-1X′Z，力矩的處理方式和力相同，可得到耦合系數Zwδ和Mwδ。圖11 為升力和縱搖力矩擬合曲面，無因次化系數和如表4 所示。

表4 耦合流體動力系數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。耦合曲面與仿真數據擬合較好，滿足線性擬合趨勢。航行器垂向速度0.538 m/s，舵角0.157rad，前舵和后舵產生的升力皆降低，該工況對應的航行器攻角為6°，舵角為9°。上文前舵15°舵角未出現舵面升力降低的情況，表明前舵升力降低的主要原因為航行器攻角改變了周圍流場的分布，對前舵產生影響。后舵升力降低的原因除了受前舵影響外可能和航行器攻角也有關，需結合流場進行進一步分析。

如圖12 所示，縱搖力矩擬合曲面與仿真數據擬合較好。與升力耦合曲面相似，航行器攻角為6°，舵角為9°，縱搖力矩數據點未落到曲面上，且誤差較大，表明此時舵的調節效果降低，開始出現非線性變化規律，增加舵角甚至還會起到反作用。

圖12 速度舵角縱搖力矩擬合曲面Fig.12 Fitting surface of velocity and pitch moment about rudder

4 舵效降低原因分析

根據仿真結果，前、后舵在航行器有攻角情況下更易發生舵效降低的情況。為了研究該種現象產生的原因，對前、后舵周圍流場進行分析。圖13 為航行器攻角0°和6°，前、后舵9°情況下流線分布圖。為了更好描述，將靠近航行器主體的舵面稱為貼體端，另一面為自由端，帶·的圖航行器攻角為0°。

圖13 前、后舵流線Fig.13 Streamlines of front and rear rudder

如圖13(a)和圖13(b)所示，自由端稍渦由壓力面和吸力面兩部分渦系組成，在舵板后緣匯聚。舵板前緣流線能夠快速成型并向下遷移，靠近中后緣的流線會被稍渦捕獲。航行器攻角增加了吸力面流線和稍渦脫落速度。

如圖13(c)和圖13(d)所示，貼體端稍渦主要由吸力面邊界層分離產生，航行器的存在阻礙了壓力面稍渦的形成，此時流體在壓力的作用下流向上表面。隨著航行器攻角的增加，前、后舵在舵面1/4 以后位置形成直徑與舵板弦長相當的表面渦，且后舵的渦更劇烈，這也是造成升力損失的主要原因。

結合前、后舵板表面極限流線進行進一步分析，舵面較小的是前舵，較大的是后舵。如圖14 所示，同一片舵中左側為舵板壓力面，右側為舵板吸力面。兩端稍渦使得極限流線發生偏轉，且吸力面偏轉更嚴重，表明吸力面受稍渦的影響更明顯。由圖15 可以看出，在兩端稍渦的影響下，極限流線偏轉更嚴重，在舵板吸力面中間形成表面渦，表面渦與兩端稍渦共同作用增加了吸力面壓力。表明航行器攻角改變了周圍流場環境，使得舵板失速角變小。