朱信堯，韓 月，朱仁慶，王樹齊

(江蘇科技大學(xué) 船舶與海洋工程學(xué)院，鎮(zhèn)江 212100)

海洋中蘊(yùn)含著豐富的資源，水下潛器及船舶是海洋開發(fā)的重要手段.微小型水下潛器絕大部分時間在深水中工作，但當(dāng)其從水面搭載平臺布放和回收、航行一段時間后上浮進(jìn)行導(dǎo)航定位及傳輸數(shù)據(jù)時，潛體受到自由液面影響.小水線面雙體船的水下潛體一般為回轉(zhuǎn)體，而其離自由液面很近[1].當(dāng)潛體在近水面附近運動時，潛體會與靜水面之間有顯著的相互作用，水面會發(fā)生明顯的變形，不再為靜水面，甚至可能會在水面上出現(xiàn)波浪破碎等現(xiàn)象，自由液面對潛體受力造成明顯影響.

文獻(xiàn)[2-6]利用STAR-CCM+對DARPPA SUBOFF AFF-1(裸艇體)模型軸線距自由液面高度(潛深)h從1.1倍最大直徑D到∞、無量綱速度Fn從0.205到0.512、側(cè)滑角β從0到18.11°時的流場進(jìn)行了模擬，分析了自由液面及粘性對潛體受力、操縱性的影響；文獻(xiàn)[7]對SUBOFF有/無附體潛艇模型無量綱速度Fn從0.128～0.84、無量綱深度h/D為1.3、2.2、3.3、4.4時的流場進(jìn)行了模擬并分析了附體對近自由液面潛體受力的影響；文獻(xiàn)[8]研究了SUBOFF模型無量綱速度Fn從0.3～0.5、無量綱深度h/D為1.1、2.2、3.3、4.4時與自由液面的相互作用；文獻(xiàn)[9-10]利用ANSYS CFX及STAR-CCM+對Afterbody-1模型潛深h/D從0.75到4.0、速度從0.4 m/s到1.4 m/s時的流場進(jìn)行了模擬并分析了自由液面對潛體受力的影響；文獻(xiàn)[11]利用CFD方法對回轉(zhuǎn)體Afterbody-1模型的流場進(jìn)行了模擬，研究了潛深、速度等對潛器近水面受力的影響；文獻(xiàn)[12]利用CFD方法對一魚雷狀水下潛器在潛深h/D為1.5、2.5、3.5三種工況下的流場進(jìn)行了模擬，研究了近水面航行器的水動力性能；文獻(xiàn)[13]利用泰勒展開邊界元法對近水面潛艇迎浪狀態(tài)下的垂向二階波浪力進(jìn)行了計算.

國內(nèi)外學(xué)者研究潛體與自由液面耦合作用時潛深h/D一般取值較大，對h/D為0.5～0.75的研究很少.文中基于CFD理論，利用重疊網(wǎng)格方法對DARPPA SUBOFF裸艇體模型h/D從0.5到3.3時的流場進(jìn)行模擬，研究潛深對近水面潛體與自由液面耦合作用的影響.

1 數(shù)值方法

采用雷諾時均N-S方程方法(RANS，Reynolds Averaged Navier-Stokes equation)求解潛體近水面粘性流動問題.主要求解模型設(shè)置為:

(1) 湍流模型采用SSTk-ω湍流模型；

(2) 自由液面捕捉采用VOF方法；

(3) 采用瞬態(tài)分析，時間步長設(shè)置為0.01 s；

(4) 耦合方式：壓力速度耦合方式采用Coupled耦合求解方式；

(5) 離散方式：壓力采用PRESTO！離散方式；體積分?jǐn)?shù)采用compressive離散方式；動量采用Second Order Upwind離散方式；

(6) 收斂監(jiān)視：采用RMS殘差監(jiān)視方式，殘差值取10-6；以迭代時所監(jiān)測力恒定或成周期性的變化作為收斂與否的判定標(biāo)準(zhǔn).

2 幾何模型、計算工況及網(wǎng)格劃分

2.1 幾何模型

以DARPPA SUBOFF AFF-1為研究對象[16]，模型的主要參數(shù)如表1，模型三維圖如圖1.

潛深h與艇體最大直徑D之間的比值定義為h*，即h*=h/D.

表1 SUBOFF模型主要參數(shù)

圖1 DARPPA SUBOFF AFF-1模型

2.2 計算工況

為了研究潛深對近水面潛體與自由液面耦合作用的影響，文中對SUBOFF模型在傅汝德數(shù)Fn為0.205、0.359、0.512，潛深h*為0.50、0.55、0.60、0.65、0.70、0.75、0.80、0.90、1.00、1.10、2.20、3.30時共36種工況下的流場進(jìn)行模擬.傅汝德數(shù)Fn為0.205、0.359、0.512時速度分別為1.248 621、2.186 61、3.118 508 m/s.

2.3 網(wǎng)格劃分

對于潛體近水面水動力計算這類繞流問題，一般通過建立一個適當(dāng)?shù)奶摂M邊界，將繞流問題轉(zhuǎn)化為內(nèi)流問題，在虛擬邊界與潛體形成的空間區(qū)域內(nèi)求解RANS方程.外圍流場采用長方體區(qū)域的形式，其計算域邊界條件及尺寸如圖2.由于流場關(guān)于縱剖面對稱，為了提高計算精度、減小計算時間，所選計算域為關(guān)于縱剖面xoz對稱的一半?yún)^(qū)域.

圖2 計算域Fig.2 Computational domain

由于結(jié)構(gòu)化的網(wǎng)格便于調(diào)整網(wǎng)格疏密、便于處理邊界層的網(wǎng)格并能夠提高數(shù)值仿真結(jié)果的精度，文中采用結(jié)構(gòu)化的六面體網(wǎng)格形式.

當(dāng)潛深h*為0.5時，潛體上表面與自由液面相切，為了保證此時的網(wǎng)格質(zhì)量，采用重疊網(wǎng)格的方法進(jìn)行網(wǎng)格劃分.

(1) 部件網(wǎng)格

部件網(wǎng)格為SUBOFF周圍區(qū)域的網(wǎng)格，其區(qū)域從為11D，頭部前面為1D，直徑為1.6D.對其流域劃分結(jié)構(gòu)化網(wǎng)格，當(dāng)速度為3.118 508 m/s時，為了保證第一層邊界層y+在30～120，第一層網(wǎng)格尺寸取0.7 mm.SUBOFF周圍重疊網(wǎng)格如圖3.網(wǎng)格數(shù)量為78萬，網(wǎng)格質(zhì)量都大于0.36，質(zhì)量較高.

圖3 重疊網(wǎng)格

(2) 背景網(wǎng)格

背景網(wǎng)格區(qū)域取長45D、高13D、半寬6D.對其流域劃分結(jié)構(gòu)化網(wǎng)格，在自由液面上下各0.2 m加密，垂向網(wǎng)格尺寸為0.01 m.背景網(wǎng)格如圖4.網(wǎng)格數(shù)量為102萬，網(wǎng)格質(zhì)量都大于0.97，質(zhì)量很高.

圖4 背景網(wǎng)格Fig.4 Background mesh

3 計算準(zhǔn)確性研究與驗證

3.1 網(wǎng)格收斂性研究

以劃分的網(wǎng)格參數(shù)為基準(zhǔn)，對SUBOFF AFF-1模型無限水深0攻角及潛深h*為1.1時網(wǎng)格的收斂性進(jìn)行研究，計算工況及網(wǎng)格參數(shù)如表2.

表2 網(wǎng)格收斂性研究網(wǎng)格參數(shù)

阻力系數(shù)CT、升力系數(shù)CL及俯仰力矩系數(shù)my計算公式分別為：

(1)

(2)

(3)

式中：Fx,Fz,My分別為阻力、升力及俯仰力矩；ρ為水的密度；v為潛器速度；Ω及L為潛器的濕表面積及長度.

為了與文獻(xiàn)[17]中試驗結(jié)果對比，無限水深0攻角計算速度取3.045 248 m/s，計算結(jié)果如表3，結(jié)果隨網(wǎng)格變化如表4.

表3 無限水深計算結(jié)果

表4 無限水深計算結(jié)果隨網(wǎng)格變化百分比

為了與文獻(xiàn)[2]中試驗結(jié)果對比，潛深h*=1.1時計算速度取3.118 508 m/s，阻力收斂曲線如圖5，從圖5可以看出，120 s后阻力波動很小，說明此時已經(jīng)收斂，結(jié)果處理時取120 s后的平均值.計算結(jié)果如表5，結(jié)果隨網(wǎng)格變化如表6.

圖5 h/D=1.1時阻力收斂過程曲線Fig.5 Drag convergence curve of h/D=1.1

表5 潛深h*=1.1時計算結(jié)果

表6 潛深h*=1.1時計算結(jié)果隨網(wǎng)格變化百分比

從表4可以看出，無限水深時，工況2-3時結(jié)果變化僅為-0.16%；再增加網(wǎng)格密度，對計算機(jī)硬件要求太高，計算時間更長，且對結(jié)果影響很小，所以網(wǎng)格尺度達(dá)到工況2時已經(jīng)滿足精度要求；

從表6可以看出，潛深h*=1.1時，除去俯仰力矩系數(shù)外，工況5-6時結(jié)果變化均小于0.5%；俯仰力矩系數(shù)變化較大的主要原因是俯仰力矩系數(shù)本身就是一個很小的量，導(dǎo)致求解時相對變化較大.

綜合考慮計算時間及計算結(jié)果的穩(wěn)健性，文中計算時取比例為100%的網(wǎng)格進(jìn)行計算分析.

3.2 計算準(zhǔn)確性驗證

為了驗證求解的精度及可信性，采用文中計算方法對DARPPA SUBOFF AFF-1模型無限水域及近自由液面的流場進(jìn)行了模擬，并將計算阻力與試驗結(jié)果進(jìn)行了比較.

(1) 無限水域

利用2.3劃分的網(wǎng)格，對SUBOFF AFF-1模型0攻角情況下的流場進(jìn)行了模擬，得到其阻力并與文獻(xiàn)[17]中試驗結(jié)果對比如表7.

從表7可以看出，文中CFD計算阻力系數(shù)和試驗結(jié)果相差小于2.89%，精度很高.

表7 無限水域準(zhǔn)確性驗證比

(2) 近自由液面

利用劃分的網(wǎng)格，對SUBOFF AFF-1模型軸線距自由液面距離h為1.1D、2.2D、3.3D情況下0攻角時傅汝德數(shù)Fn為0.512時的流場進(jìn)行模擬.

圖6為潛器縱剖面底部y+值沿x方向的變化.文中計算所得阻力與文獻(xiàn)[2]中試驗結(jié)果對比，如表8.

圖6 潛器縱對稱面底部曲線y+Fig.6 y+ of the bottom line of SUBOFF

表8 近自由液面準(zhǔn)確性驗證

從圖6可以看出，y+在30～110，適合對潛器流場進(jìn)行求解.

從表8可以看出，文中近自由液面CFD計算阻力與試驗相比誤差在2.53%以內(nèi)，精度很高，滿足計算精度要求.文中CFD計算結(jié)果精度很高的主要原因是網(wǎng)格劃分方式采用結(jié)構(gòu)化網(wǎng)格，并對潛器邊界層網(wǎng)格及自由液面網(wǎng)格進(jìn)行了加密，網(wǎng)格更加細(xì)致.

4 結(jié)果與分析

文中數(shù)值模擬了Fn分別為0.205、0.359、0.512時，潛深h*為0.50、0.55、0.60、0.65、0.70、0.75、0.80、0.90、1.00、1.10、2.20、3.30時共36種工況下的流場.文中對近水面潛體與自由液面耦合作用的阻力、升力、俯仰力矩及自由液面情況進(jìn)行了總結(jié)分析.

4.1 阻力

Fn分別為0.205、0.359、0.512時潛器總阻力系數(shù)CT隨無量綱潛深h*的變化曲線如圖7.

圖7 總阻力系數(shù)隨潛深的變化曲線Fig.7 Total drag coefficient over depth

潛器所受總阻力可以分為粘性摩擦阻力及剩余阻力，為了便于分析自由液面對阻力的影響，將總阻力系數(shù)CT分為摩擦阻力系數(shù)Cf及剩余阻力系數(shù)Cp，其關(guān)系為：

CT=Cf+Cp

(4)

將近自由液面時的總阻力系數(shù)CT減去無限水深時的總阻力系數(shù)CT∞可以得到由于自由液面影響而產(chǎn)生的阻力系數(shù)，即興波阻力系數(shù)Cw，其關(guān)系為：

Cw=CT-CT∞

(5)

無限水深時的總阻力系數(shù)CT∞可由CFD計算得到，其阻力系數(shù)如表9.Fn=0.205、0.359、0.512時不同水深時潛器所受的總阻力系數(shù)CT、摩擦阻力系數(shù)Cf、剩余阻力系數(shù)Cp及興波阻力系數(shù)Cw結(jié)果如表10.

表9 無限水深阻力系數(shù)

表10 近水面阻力系數(shù)結(jié)果

從圖7及表9、10可以看出，當(dāng)潛深h*從0.50增大到3.3時：

(1) 當(dāng)潛器高速航行時(Fn=0.512)，總阻力系數(shù)先增大后減小，在潛深h*為0.65時有最大值，其值為無限水域阻力系數(shù)的4.79倍；

(2) 當(dāng)潛器航速中等時(Fn=0.359)，總阻力系數(shù)呈先減小后增大再減小的趨勢，在潛深h*為0.80時有極大值，其值為無限水域阻力系數(shù)的1.85倍；最大值發(fā)生在h*為0.5時，其值為無限水域阻力系數(shù)的2.32倍；

(3) 當(dāng)潛器低速航行時(Fn=0.205)，總阻力系數(shù)變化趨勢與航速中等時類似，也先減小后增大最后減小，在h*為1.00時有最大值，其值為無限水域阻力系數(shù)的2.34倍；當(dāng)h*為0.5時，其值為無限水域阻力系數(shù)的2.21倍；

(4) 3種航速下摩擦阻力系數(shù)變化均較小，說明潛深對摩擦阻力影響較??；

(5) 興波阻力系數(shù)：① 興波阻力受自由液面影響很大，其變化趨勢與總阻力變化趨勢幾乎一致，并非離自由液面越近，興波阻力越大；總阻力增加主要是由于興波引起的；② 當(dāng)潛深h*為3.30時，興波阻力占總阻力的比例非常小，在Fn=0.205時為3.08%，在Fn=0.359時為7.62%，在Fn=0.512時為16.14%；③ 當(dāng)潛深h*距自由液面很近時，興波阻力占總阻力的比例變的很大，在Fn=0.205時，興波阻力占總阻力的比例最大值出現(xiàn)在h*為1.00時，其比例為52.74%；在Fn=0.359時，最大值出現(xiàn)在h*為0.50時，其比例為56.87%；在Fn=0.512時，最大值出現(xiàn)在h*為0.65時，其比例為79.14%.

4.2 升力

文中升力為潛器垂直方向上的所有力的合力，其由水動升力與靜水浮力的和減去重力所得.為了保證深水時潛器處于零浮力狀態(tài)，重力與潛器完全浸沒在水中時的靜水浮力數(shù)值相等.

傅汝德數(shù)Fn為0.205、0.359、0.512時潛器升力隨無量綱潛深h*的變化曲線如圖8.

圖8 升力隨潛深的變化曲線Fig.8 Lift over depth

從圖8可以看出，當(dāng)潛深h*從0.50增大到3.3時：

(1) 潛器升力總體呈先增大后趨于0的趨勢.在h*較大時，潛器距自由液面較遠(yuǎn)，受自由液面影響較小，由于SUBOFF潛器為回轉(zhuǎn)體，故其升力在0附近；而當(dāng)h*較小時，其離水面較近，出現(xiàn)明顯的露背現(xiàn)象，導(dǎo)致升力減小.

(2) 當(dāng)潛深h*較小時，速度越高，潛器升力受自由液面影響越大.當(dāng)潛深h*為0.5、Fn=0.512時潛器所受升力與Fn=0.205時相差444.7 N.

4.3 俯仰力矩

傅汝德數(shù)Fn為0.205、0.359、0.512時潛器俯仰力矩系數(shù)隨無量綱潛深h*的變化曲線如圖9.

圖9 俯仰力矩系數(shù)隨潛深的變化曲線Fig.9 Pitching moment coefficient over depth

從圖9可以看出，當(dāng)潛深h*從0.50增大到3.3時，俯仰力矩系數(shù)絕對值先增大后趨于0.在h*較大時，潛器距自由液面較遠(yuǎn)，受自由液面影響較小，由于SUBOFF潛器為回轉(zhuǎn)體，故其俯仰力矩系數(shù)在0附近；而當(dāng)h*較小時，其離水面較近，受自由液面影響較大.

潛深h*為0.8時，潛器升力系數(shù)、俯仰力矩系數(shù)沿潛器軸向分布曲線如圖10.Fn=0.512時，潛器升力系數(shù)、俯仰力矩系數(shù)沿潛器軸向分布曲線如圖11.潛器浮心縱向位置距頭部2.008 6 m，從圖10、11可以看出，潛器產(chǎn)生俯仰力矩主要是由于升力分布不均引起的.

圖10 潛深h*為0.8時升力系數(shù)、俯仰力矩系數(shù)沿潛器軸向分布Fig.10 Distribution of the lift coefficient and pitching moment coefficient along the length of the SUBOFF at the submergence depth h*=0.8

圖11 Fn=0.512時升力系數(shù)、俯仰力矩系數(shù)沿潛器軸向分布Fig.11 Distribution of the lift coefficient and pitching moment coefficient along the length of the SUBOFF at Fn=0.512

4.4 自由液面

潛器近水面航行時，自由液面最大波面升高A與艇體最大直徑D的比值定義為A*，即A*=A/D.

傅汝德數(shù)Fn為0.205、0.359、0.512時自由液面最大波面升高A*隨無量綱潛深h*的變化曲線如圖12.

傅汝德數(shù)Fn為0.512時潛器縱對稱面波形圖如圖13，潛深h*為0.8時潛器縱對稱面波形圖如圖14.

Fn為0.512時，潛深h*為0.5、0.8、1.1及3.3時潛器縱對稱面動壓分布云圖分別如圖15.

圖12 無量綱最大波面升高隨潛深的變化曲線Fig.12 Normalized maximum wave height over depth

圖13 Fn=0.512時潛器縱對稱面波形Fig.13 Wave pattern on longitudinal plane at Fn=0.512

圖14 潛深h*為0.8時潛器縱對稱面波形Fig.14 Wave pattern on longitudinal plane at h*=0.8

圖15 Fn=0.512時潛器縱對稱面動壓分布云圖Fig.15 Distribution of the dynamic pressure on longitudinal plane at Fn=0.512

對自由液面興波進(jìn)行分析，可以得到以下結(jié)論：

(1) 從圖12可以看出，無量綱最大波面升高A*總體呈先增大后減小的趨勢，在h*在0.80左右時有最大值.速度為Fn=0.512時最大波面升高達(dá)到了潛器直徑D的0.4倍.當(dāng)潛深h*從0.80增大時最大波面升高減小的主要原因是潛器離自由液面越遠(yuǎn)，受自由液面影響越小，產(chǎn)生的興波越小；當(dāng)潛深h*從0.50增大到0.80時最大波面升高增大的主要原因是潛器近水面興波的最大波面升高并不是發(fā)生在頭部附近，而是發(fā)生在尾流場，當(dāng)h/D為0.50時，潛體上表面與靜水自由液面相切，上面流體較少，導(dǎo)致對尾部流場影響較小，最大波面升高相對較小，此從圖13、14、15可以看出.

(2) 從圖12、14可以看出，相同潛深時，速度越大，自由液面升高越大，說明速度越大，潛器對自由液面作用越明顯；

(3) 從圖13、14、15可以看出，潛器近水面興波的最大波面升高并不是發(fā)生在頭部附近，而是發(fā)生在尾流場；

(4) 從圖13、14可以看出，興波波長受速度影響較大，受潛深影響較小.

興波理論中波長λ與速度U的關(guān)系為[3]：

(6)

從圖13、14可以看出，在Fn=0.512時，不同潛深時產(chǎn)生的波長均在6 m左右；在Fn=0.359時，產(chǎn)生的波長約為3 m；在Fn=0.205時，產(chǎn)生的波長約為1 m；此波長計算結(jié)果與興波理論中波長結(jié)果一致.

5 結(jié)論

當(dāng)潛深h*從0.50增大到3.30時：

(1) 阻力：當(dāng)潛器高速航行時，總阻力系數(shù)先增大后減小，當(dāng)潛器中、低速航行時，總阻力系數(shù)呈先減小后增大再減小的趨勢；自由液面對潛器總阻力影響很大，各速度下最大阻力系數(shù)是其無限水深時的2.32～4.79倍；隨著潛深的變化，摩擦阻力變化很小，而興波阻力變化很大；興波阻力變化趨勢與總阻力的幾乎一致；總阻力增加主要是由于興波產(chǎn)生的.

(2) 升力：潛器升力總體呈先增大后趨于0的趨勢；而當(dāng)潛深h*較小時，其離水面較近，出現(xiàn)明顯的露背現(xiàn)象，導(dǎo)致升力減小，此時，速度越高，潛器升力受自由液面影響越大.

(3) 俯仰力矩系數(shù)絕對值先增大后趨于0.

(4) 無量綱最大波面升高A*先增大后減小，在h*在0.80左右時有最大值；相同潛深時，速度越大，自由液面升高越大；潛器近水面興波的最大波面升高并不是發(fā)生在頭部附近，而是發(fā)生在尾流場；自由表面興波波長受速度影響較大，受潛深影響較小.