陳小星 趙林建

（第七一五研究所，杭州，310023）

隨著主動拖曳聲吶的廣泛使用，對拖體的水動力數(shù)值模擬方法研究也逐漸在設(shè)計中廣泛采用，相對以往的工程近似估算，該方法精度有明顯提高。同時與水池試驗相比，具有成本低、周期短等優(yōu)點。文獻[1，2]對粘性水動力系數(shù)進行了數(shù)值計算方法研究，采取Fluent軟件進行計算，通過二次編程增加動量源項，將旋臂試驗轉(zhuǎn)換為定常計算研究。本文以Star-CCM+為計算軟件，采用場函數(shù)直接添加動量源項，避免了復(fù)雜的編程。通過數(shù)值計算仿真并與實測數(shù)據(jù)進行了比對，表明該方法可以用于方案設(shè)計階段，計算結(jié)果可以達到工程使用要求。

1 數(shù)值模型

1.1 源項計算

在慣性坐標(biāo)系下，流體的連續(xù)性方程為

動量RANS（Reynolds-Averaged Navier-Stokes）方程為

式中，ρ表示流場密度，U為流場速度，f為體積力，為流體應(yīng)力張量，在角速度為常數(shù)Ω的旋轉(zhuǎn)參考系下，連續(xù)性方程沒有速度的導(dǎo)數(shù)項，保持不變，RANS方程變?yōu)?/p>

式中，Ur為相對速度，r為旋轉(zhuǎn)半徑。相對慣性坐標(biāo)系下的增加項MS為

額外增加的體積力MS，可以作為動量源項添加進流場中，從而實現(xiàn)穩(wěn)態(tài)計算，避免了Star-CCM+中采用旋轉(zhuǎn)參考系的復(fù)雜做法。

1.2 k-e湍流模型

湍動能k、湍流耗散率ε的傳輸方程為

式中，為平均速度，μ為動力粘度，ρ為流體密度，μt為湍流渦粘度，σk、σε、Cε1、Cε2為模型系數(shù)，Sk、Sε為用戶指定的源項。

1.3 k-ω湍流模型

k-ω湍流模型是雙方程模型，它可對湍動能k和單位耗散率ω的傳輸方程進行求解，以確定湍流渦粘度。一般適用于強旋湍流，本文中懸臂計算采用k-ω湍流模型。其傳輸方程為

式中，β*、β為模型系數(shù)，Pk、Pω為結(jié)果項。

2 直航試驗

2.1 流體域及網(wǎng)格劃分

為消除計算域壁面對流動的影響，流域范圍應(yīng)越大越好，但這必然會導(dǎo)致網(wǎng)格數(shù)量過大影響計算效率。本文選擇矩形計算流域，為保證模型周圍流場得到充分發(fā)展，長寬高為13 m×6 m×6 m，采用六面體網(wǎng)格。

采用k-e湍流模型進行仿真計算，確保Wall Y+都在30～300之間。首層邊界層厚度公式如下[3]：

式中，Y+為期望的取值，這里取100，L為拖體特征長度，Re為雷諾數(shù)。邊界層數(shù)量為6層，邊界層總厚度為0.02 m。計算收斂后，查看計算的Y+分布如圖1所示，大部分拖體表面的Y+值都在30～300之間。

圖1 拖體表面Wall Y+分布圖

對模型進行了網(wǎng)格無關(guān)性驗證，建立了四套網(wǎng)格分別進行計算，按照ITTC（International Towing Tank Conference）推薦設(shè)定，BaseSize按2比率縮小，計算結(jié)果見表1。選擇網(wǎng)格4的大小，模型的阻力變化不大，本文計算選擇BaseSize為0.025 m。

航標(biāo)是引導(dǎo)船舶安全航行的重要設(shè)施，是我國海事安全保障體系的重要組成部分，航標(biāo)可靠的導(dǎo)助航效能是船舶安全通航的保障。伴隨著科學(xué)技術(shù)的發(fā)展，航標(biāo)科技的發(fā)展也越來越快，我們有理由相信，未來航標(biāo)維護管理工作肯定會比現(xiàn)在更加的科學(xué)、安全和可靠。

表1 網(wǎng)格數(shù)量與阻力關(guān)系表

2.2 邊界條件和求解器設(shè)置

本文一共使用了4種邊界條件：（1）上游壁面為速度入口，給定來流方向及大小；（2）四周壁面為滑移壁面；（3）下游壁面為壓力出口；（4）拖體表面為無滑移壁面。使用分離求解器，離散方式采用二階迎風(fēng)格式，收斂標(biāo)準(zhǔn)選取為10-5，迭代步數(shù)為1 000步。

2.3 試驗結(jié)果

直航試驗主要計算各個拖曳速度下的縱向阻力，用來確定縱向阻力系數(shù)速度下的中縱剖面的速度分布及壓力分布見圖2、圖3。計算阻力和實驗阻力擬合曲線見圖4。

圖2 中縱剖面速度分布圖

圖3 中縱剖面壓力分布圖

圖4 各拖曳速度下的阻力

對計算獲得的數(shù)據(jù)進行整體線性回歸分析，取得拖體的縱向阻力系數(shù)，無因次化后結(jié)果為與水池試驗的無因次結(jié)果進行比較，誤差為6.7%。

3 斜航試驗

斜航試驗用來確定由速度引起的水動力系數(shù)，主要為位置導(dǎo)數(shù) (水平面及垂直面位置導(dǎo)數(shù)通過擬合Y-v、N-v、Z-w、M-w曲線獲得。斜航試驗的計算工況見表2，通過調(diào)整拖體的偏轉(zhuǎn)角度計算各個流體力及流體力矩。所有工況的來流速度為2m/s。

表2 斜航試驗數(shù)值模擬計算工況

斜航試驗采用與阻力試驗相同的網(wǎng)格尺寸、邊界條件和求解器設(shè)置，只是模型進行偏轉(zhuǎn)。其Y-v、N-v、Z-w、M-w曲線的對比如圖5～8所示。

對仿真結(jié)果采用整體回歸分析方法，得到水平面和垂直面的位置導(dǎo)數(shù)，進行無因次化后的結(jié)果見表3。除了結(jié)果偏差略大，其余的水動力系數(shù)誤差均在10%以內(nèi)。

圖5 俯仰力矩M-垂向速度w曲線

圖6 垂向力Z-垂向速度w曲線

圖7 偏航力矩N-橫向速度v曲線

圖8 側(cè)向力Y-橫向速度v曲線

表3 斜航試驗數(shù)值模擬計算結(jié)果

4 旋臂試驗

旋臂試驗主要用來確定由角速度引起的水動力系數(shù)，主要是旋轉(zhuǎn)導(dǎo)數(shù) (水平面以及垂直面旋轉(zhuǎn)導(dǎo)數(shù)通過擬合Y-r、N-r、Z-q、M-q曲線獲得。旋臂試驗的計算工況見表4。通過調(diào)整拖體的回轉(zhuǎn)半徑計算各個流體力及流體力矩。所有工況的來流速度為2m/s。

表4 旋臂試驗數(shù)值模擬計算工況

旋臂試驗采用回轉(zhuǎn)體流域，邊界條件和求解器設(shè)置與斜航試驗相同，為了提高收斂效率，采用多面體網(wǎng)格k-ω湍流模型進行仿真計算。

對增加的源項MS，在求解區(qū)域中按式（4）增加了動量源項場函數(shù)。同時入口速度指定為Ω×r場函數(shù)。

邊界層厚度確定按照式（9）進行計算，取Y+=4，邊界層數(shù)量20層，邊界層總厚度為0.03 m。計算收斂后，Y+分布見圖9，圖中大部分拖體表面的Y+值都在5以下，其Y-r、N-r、Z-q、M-q計算曲線見圖10～11，其回轉(zhuǎn)半徑R=30 m下的速度分布見圖12。

圖9 拖體表面Wall Y+分布圖

圖10 側(cè)向力Y、偏航力矩N-偏航角速度r曲線

圖11 垂向力Z、俯仰力矩N-縱傾角速度q曲線

圖12 2m/s速度、R=30 m時回轉(zhuǎn)速度分布圖

對仿真結(jié)果采用整體回歸分析方法，得到水平面和垂直面的旋轉(zhuǎn)導(dǎo)數(shù)，進行無因次化后的結(jié)果見表5。從表中數(shù)據(jù)可知，誤差都在20%以內(nèi)。

表5 旋臂試驗數(shù)值模擬計算結(jié)果

5 結(jié)論

本文通過仿照拖體的運動形式，采用兩種湍流模型，對拖體水動力系數(shù)進行了數(shù)值模擬計算，通過使用附加動量源項方法，將瞬態(tài)計算轉(zhuǎn)換為穩(wěn)態(tài)計算，大大縮短了計算時間。同時將計算結(jié)果與試驗數(shù)據(jù)進行了比對，驗證了該方法的可行性，計算結(jié)果可以滿足工程使用精度要求。本文的研究可以應(yīng)用于拖體方案設(shè)計階段的穩(wěn)定性和操縱性預(yù)報，從而縮短開發(fā)周期，降低開發(fā)成本，對拖體的水動力系數(shù)計算有一定的參考價值。