孫現(xiàn)有，劉顯龍

(1.海軍駐昆明七五〇試驗場軍事代表室；2.中國船舶重工集團公司七五〇試驗場，云南 昆明 650051)

0 引言

拖纜作為水面及水下武器的承載設備發(fā)揮著巨大的作用。常見的拖纜拖曳形式有水面艦拖曳航行體，水下航行體拖曳聲學線列陣等。為了計算整個拖纜的阻力性能，需要將其參數(shù)化，計算各個的阻力性能，并將其整體矢量求和。

以往常采用經(jīng)驗公式法來計算拖纜微元的阻力性能[1]。經(jīng)驗公式是將來流速度按拖纜的切向、法向2個方向進行分解，再由公式計算切向力和法向力[2]。采用這種方法由于未考慮切向與法向速度相互干擾因素，結(jié)果不夠準確，連璉等人[3]通過試驗的方法嚴格說明了此種方法帶來的誤差，本文不再贅述。隨著計算流體力學(CFD)技術的發(fā)展，采用CFD方法能夠快速地計算流體性能。另外，CFD方法僅能計算有限個攻角時拖纜的阻力性能，為了得到實際不同攻角的速度，需要采用三次樣條插值的方法來插值求得。

1 流體力學基礎理論

1.1 控制方程

不可壓縮粘性流體的連續(xù)性方程和動量方程為

(1)

(2)

1.2 湍流模型

要使上述方程封閉，須引入新的湍流模型方程，把應力項中的脈動值與時均值聯(lián)系起來。RNGk-ε方程中湍流動能k和湍流脈動強度ε的輸運方程為

(3)

(4)

常數(shù)為：G1ε=1.42；G2ε=1.68；Cμ=0.0845；

σk=1.0；σε=1.3；η0=4.38；β=0.012。

2 CFD模擬

2.1 流體計算域

拖纜微元與迎流速度的夾角(即攻角)從0°到90°范圍變化。0°及90°情況較為簡單，可以較快地得到結(jié)構化網(wǎng)格，較為復雜的是0°至90°之間的情況。對于此種情況，將速度入口由平面改為弧形面，可以較好地解決這一問題[4]，其流體計算域及邊界設置情況見圖1。

需要說明的是，拖纜微元在不同攻角時阻力性能是不同的，而計算時又不可能由CFD的方法計算每種攻角情況，為了解決這一問題，采用三次樣條插值的方法[5]。

2.2 網(wǎng)格劃分

為了得到高質(zhì)量的網(wǎng)格，對整個計算域采取分塊化處理，再劃分為結(jié)構化網(wǎng)格。在拖纜周圍設置邊界層網(wǎng)格，第一層控制近壁面的網(wǎng)格高度和質(zhì)量，近壁面高度按估算公式(5)確定：

(5)

式中：Δy為第一層邊層界厚度；L為長度值；Re為雷諾數(shù)。實際取y+為30，整體網(wǎng)格數(shù)量為90萬個。拖纜流場網(wǎng)格見圖2。

2.3 邊界條件

邊界條件的設定是保證CFD實現(xiàn)的必要條件之一。本文計算域的邊界條件設定見圖1，包括：速度入口、壓力出口、壁面和對稱面。

1)速度入口：在進流面設定迎流速度的大小和方向，Vin=V0；

2)壓力出口：在出流面設定相對于參考壓力點的流體靜壓力，Pout=0；

3)壁面條件：潛艇表面及圓柱表面設定為不可滑移邊界條件，u=v=w=0；

2.4 數(shù)值計算

采用有限體積法離散控制方程和湍流模式，并用耦合隱式求解法計算。在差分格式中，壓力項使用標準格式，動量項、湍流動能項以及耗散率項先使用一階迎風格式迭代500次，再采用二階迎風格式[6]。

采用RNGk-ε湍流模型，并采用有限體積法對控制方程進行離散。壓力-速度耦合迭代求解采用協(xié)調(diào)一致的SIMPLEC算法，壓力項采用標準離散格式外，動量、湍流動能和湍流耗散率采用二階迎風格式，亞松弛因子均采用默認值。近壁處區(qū)域采用壁面函數(shù)法[7]。

3 三次樣條插值

在實際計算過程中，通過CFD方法計算拖纜微元在攻角分別為0°、10°、20°、…、90°的阻力情況。實際計算時若碰到不同屬于上述攻角情況時，采用三次樣條插值的方法進行估算[8]。

曲線插值一般有指數(shù)函數(shù)插值、拉格朗日插值、三次樣條插值、B樣條曲線、NURBS曲線等方法。這里采用三次樣條進行插值。

在區(qū)間[a,b]內(nèi)給定一個分劃Δ：a=x0

1)每個子區(qū)間[xi-1,xi],i=1,2,…,n上為三次多項式；

2)y(x)在整個區(qū)間[a,b]上二階連續(xù)可導。

則稱y(x)為區(qū)間[a,b]上關于分劃Δ的三次樣條函數(shù)，其中xi(i=0,1,…,n)稱為節(jié)點。插值點如表1所示。

表1 三次樣條曲線插值點

則三次樣條插值函數(shù)為

(6)

式(6)中：

hk=xk+1-xk(k=0,1,…,n-1)

(7)

λkmk-1+2mk+μkmk+1=gk(k=1,2,…,n-1)

(8)

其中：

4 拖纜計算方法簡介

由于拖纜是柔性的，受力將發(fā)生變形[8]。因此，各處拖纜的來流攻角是不同的。采用有限差分法，將拖纜處理成有限個小段，對各小段分別受力分析，由初始點的力、轉(zhuǎn)角條件求出未知位置的力、轉(zhuǎn)角及坐標點。

為了滿足工程中的需要，且基本滿足實際情況，做出以下幾點假設：

1)拖纜在張力作用下的伸縮忽略不計；

2)拖纜為柔性的，不傳遞彎矩；

3)整個流動為二維流動；

4)作用于拖纜微元上的水動力系數(shù)等同于作用于無限長的拖纜。

將一小段拖纜微元進行受力分析，見圖3。

按照切向力與法向力方向進行分解，由力的平衡可得：

(T+dT)sin(dφ)+pdscos(φ)=Dds

(9)

(T+dT)cos(dφ)=pdssin(φ)+Fds+T

(10)

式 (9)、(10)中：T為拉力；D為單位長度拖纜法線方向的力(由CFD及三次樣條插值計算)；F為單位長度拖纜切線方向的力(由CFD及三次樣條插值計算)；p為單位長度拖纜在水中重力。

忽略二階小量，由式(10)得：

(11)

將式(11)代入式(9)得：

(12)

由式(11)、(12)可求出不同位置拖纜拉力、轉(zhuǎn)角等值。

另外：

dx=ds×cosφ

(13)

dy=ds×sinφ

(14)

由式(13)、(14)可求出相應位置點的坐標值。

5 計算結(jié)果分析

在實際計算時，可能會遇到不同直徑的拖纜。眾所周知，在CFD計算中，最為復雜的是采用ICEM來劃分結(jié)構化網(wǎng)格，耗費了大量的人力、物力。目前，Workbench(15.0版)后已經(jīng)將ICEM集成進去，但是當修改拖纜直徑后，ICEM的拓撲分塊處理并不能自動化解決，還是需要大量、重復性人工干預。為了較好地解決這一問題，采用ICEM自帶的腳本語言，將劃分網(wǎng)格中的建組、分塊、劃網(wǎng)格、導出網(wǎng)格全過程監(jiān)控、錄制并保存，當模型的直徑發(fā)生變化時，重新讀入一次腳本語言，即可一鍵化劃分網(wǎng)格，將原本需要2 h左右的劃分網(wǎng)格時間縮減至1 min以內(nèi)。需要說明的是，建立模型的直徑需要在Workbenchk中的DM模塊下進行，建立好模型后，僅修改其直徑等參數(shù)，以便于腳本語言無縫識別，否則在其它建模環(huán)境中，如Pro/E、Rhino等，腳本語言對點、線、面、體的識別極易出錯。

以直徑為10 mm的拖纜(長度為1 m)為例，其在迎流速度為5 kn，其切向力及法向力見表2。

表2 10 mm拖纜在不同攻角時切向力與法向力(v=5 kn)

以攻角為0°、45°和90°為例，其壓力云圖見圖4。需要說明的是，表2中攻角沒有計算接近0°和90°的角度，因為此時網(wǎng)格斜角過大，劃分網(wǎng)格質(zhì)量不夠好，影響了計算的精度，采用三次樣條插值進行插值計算。攻角在5°時，計算得到切向力0.48 N，法向力0.32 N。

6 結(jié)束語

上述計算表明，本文采用的方法有以下優(yōu)勢：

1)采用CFD的方法計算拖纜微元的阻力性能，考慮了切向速度與法向速度間干擾，優(yōu)于以往經(jīng)驗公式；

2)采用參數(shù)化的方法建模、劃分結(jié)構化網(wǎng)格，大大節(jié)省了前處理時間；

3)采用三次樣條插值的方法，可以快速計算實際中任一攻角的阻力性能。

綜合前三者優(yōu)勢，結(jié)合編程的方法，可以將整個拖曳系統(tǒng)(含多節(jié)點)在不同速度拖曳時的總阻力及姿態(tài)情況計算出來。