密度分層對(duì)近水面航行潛體興波阻力影響分析

何廣華， 劉雙， 王威， 潘雁甲

(1.哈爾濱工業(yè)大學(xué) 機(jī)電工程學(xué)院，黑龍江 哈爾濱 150001；2.山東船舶技術(shù)研究院，山東 威海 264209；3.哈爾濱工業(yè)大學(xué)(威海) 船舶與海洋工程學(xué)院，山東 威海 264209)

潛體在密度分層流中航行時(shí)，其水動(dòng)力特性較為復(fù)雜。眾所周知，潛體在單一均勻流體中航行時(shí)，會(huì)在自由液面處興起波浪，形成凱爾文波系；實(shí)際海水中由于溫、鹽差的存在會(huì)出現(xiàn)分層現(xiàn)象，潛體航行于多層流之間時(shí)，由于體積效應(yīng)及湍流尾跡效應(yīng)會(huì)在密度分層處產(chǎn)生內(nèi)波[1]，其周期較長，振幅較大，會(huì)對(duì)潛體的航行性能產(chǎn)生較大影響[2]。故研究潛體在密度分層流這一特殊海況下航行時(shí)的受力特性具有十分重要的意義。

文獻(xiàn)[3-6]對(duì)于運(yùn)動(dòng)物體在分層流體中激發(fā)內(nèi)波的現(xiàn)象進(jìn)行了研究。Gou等[7]采用時(shí)域高階邊界元法研究了兩層流體中的波衍射問題。Song等[8]對(duì)內(nèi)孤立波與海洋結(jié)構(gòu)物的相互作用進(jìn)行了數(shù)值模擬。Lin等[9]、Bonneton等[10]對(duì)密度分層流進(jìn)行了試驗(yàn)研究。針對(duì)物體在密度分層流體中的受力情況，Motygin等[11]基于勢流理論推導(dǎo)了二維物體分別在上、下層流體中運(yùn)動(dòng)所產(chǎn)生的興波阻力公式。Miloh等[12]采用格林函數(shù)法得到了半潛扁球體在2層流體中阻力及自由表面擾動(dòng)的解。Grue[13]對(duì)“弗雷姆”號(hào)的航行過程進(jìn)行研究，發(fā)現(xiàn)的受力變化規(guī)律與Ekman[14]一致。勾瑩等[15]采用模型實(shí)驗(yàn)的方法研究了箱型結(jié)構(gòu)在2層流中拖動(dòng)時(shí)的阻力特性，并與單層流中拖航阻力實(shí)驗(yàn)結(jié)果進(jìn)行了對(duì)比研究。

綜上可知，針對(duì)潛體在分層流體中的研究主要以運(yùn)動(dòng)激發(fā)內(nèi)波的研究居多，而對(duì)由于內(nèi)波存在引起的潛體受力變化研究并不多見。與受力相關(guān)的研究也大多關(guān)注水面艦船，本文著重研究潛體航行于密度分層流中的受力特性。本研究基于粘流理論、壓力-速度耦合項(xiàng)求解采用Simple算法，通過定義各流體項(xiàng)體積分?jǐn)?shù)以及計(jì)算域的壓力分布，建立了密度分層流數(shù)值水池模型，針對(duì)近水面潛體航行特性開展了模擬研究。

1 密度分層流數(shù)值模型基本理論

1.1 控制方程

本模型基于雷諾平均的N-S方程，不可壓縮流體的控制方程為：

(1)

(2)

式(1)為連續(xù)性方程，式(2)為動(dòng)量方程，兩者均為時(shí)均計(jì)算處理后的形式。

1.2 湍流模型

需要引入湍流模型使方程(2)封閉。本研究采用Realizablek-ε湍流模型，該模型穩(wěn)定性良好，壓力梯度求解精度高，工程應(yīng)用比較廣泛，其湍動(dòng)能及湍動(dòng)能耗散率輸運(yùn)方程為：

(3)

(4)

式中：μ為分子擴(kuò)散所造成的動(dòng)力粘性；μt為湍流粘性系數(shù)；C1=max[0.43,η/(η+5)]；η=Sk/ε；S為平均應(yīng)變率；C2和C1z為常數(shù)；C2z為浮力對(duì)耗散率影響的函數(shù)；Gk為平均速度梯度產(chǎn)生的湍動(dòng)能；Gb為由浮力產(chǎn)生的湍動(dòng)能；YM表示湍流脈動(dòng)膨脹對(duì)總耗散率的影響；σk和σz分別是湍動(dòng)能k和耗散率ε的湍流普朗特?cái)?shù)，其他參數(shù)取值為：C1z=1.44，C2=1.9，σk=1.0，σz=1.2。

1.3 平面進(jìn)行波理論

航行體興波存在首波系和尾波系，且兩波系中的橫波在艉部之后相遇而疊加，稱為興波干擾。艏橫波的第1個(gè)波峰和艉橫波第1個(gè)波峰之間的距離稱為興波長度，用mL表示，m為與船型和Fr有關(guān)的系數(shù)。由平面進(jìn)行波理論可得，興波阻力系數(shù)Cw為：

(5)

式中：ρ為水密度；g為重力加速度；λ為波長；v為航速；S為濕表面面積；L為潛體長；C、D為常數(shù)。

由式(5)可知，對(duì)于給定船型，其興波阻力系數(shù)Cw與長度Fr的4次方成比例。

2 密度分層流數(shù)值模型設(shè)置與驗(yàn)證

2.1 數(shù)值模擬相關(guān)設(shè)置

鑒于所研究的問題及結(jié)構(gòu)均對(duì)稱，故在潛體中縱剖面處采用對(duì)稱面條件來提高計(jì)算效率。計(jì)算域總長為15L，寬為3.5L?；跉W拉多相流模型，通過UDF的方式分別給定空氣、淡水、鹽水三相的密度及初始體積分?jǐn)?shù)。

計(jì)算域前端邊界條件為速度入口，后端為壓力出口，使用切割體網(wǎng)格，對(duì)自由表面、密度分界面以及潛體表面曲率變化較大等處進(jìn)行網(wǎng)格加密以準(zhǔn)確捕捉流場特性，在潛體重心處進(jìn)行網(wǎng)格定位。調(diào)節(jié)網(wǎng)格增長率、棱柱層厚度等參數(shù)控制網(wǎng)格質(zhì)量，經(jīng)多次數(shù)值實(shí)驗(yàn)，最終總網(wǎng)格數(shù)為400萬左右，其中棱柱層厚度為0.018 m，棱柱層數(shù)為6，增長率為1.2。計(jì)算域網(wǎng)格如圖1所示。

圖1 計(jì)算域網(wǎng)格分布Fig.1 Grid distribution of the domain

模型采用全附體標(biāo)準(zhǔn)SUBOFF標(biāo)模，其總長L為3 m，最大直徑0.174 m。計(jì)算域中流體密度沿垂向變化，分層情況如圖2所示。研究中SUBOFF潛體以不同速度U航行于不同位置處。潛體潛深為d1，指潛體中心到自由液面的垂向距離；d2為潛體位于淡水層中時(shí)其中心距內(nèi)波面的距離；h1、ρ1分別為淡水層流體的深度及密度；h2、ρ2分別為鹽水層流體的深度及密度。

圖2 計(jì)算流域分層Fig.2 Diagram of submarine in the density stratified fluid

2.2 網(wǎng)格及時(shí)間步長收斂性驗(yàn)證

采用航速Fr=0.5的全附體潛體進(jìn)行收斂性研究，此時(shí)各參數(shù)為淡水層h1=0.24L，密度取為997.561 kg/m3；鹽水層高h(yuǎn)2=0.67L，密度為1 020 kg/m3，d1=d2=0.12L。通過改變網(wǎng)格的基礎(chǔ)尺寸來控制網(wǎng)格數(shù)量，分別取網(wǎng)格基礎(chǔ)尺寸Mesh A為0.07 m，Mesh B為0.065 m，Mesh C為0.056 m，Mesh D為0.052 mm進(jìn)行網(wǎng)格收斂性研究，取時(shí)間步長0.01 s。

表1是網(wǎng)格數(shù)量改變對(duì)阻力計(jì)算結(jié)果的影響。其中誤差是指采用某一網(wǎng)格量計(jì)算結(jié)果與采用最精細(xì)網(wǎng)格Mesh D所得計(jì)算結(jié)果之間的差別，以百分比形式表示；表中Rt代表總阻力，Rf代表摩擦阻力。由表1可以看出，當(dāng)網(wǎng)格基礎(chǔ)尺寸由0.056 m降低到0.052 m級(jí)別時(shí)，對(duì)總阻力及摩擦阻力計(jì)算結(jié)果帶來的影響并不大，說明此時(shí)網(wǎng)格收斂?？紤]計(jì)算效率，在之后的計(jì)算中均采用Mesh C網(wǎng)格。

表1 網(wǎng)格收斂性分析Table 1 Mesh convergence analysis

采用4種時(shí)間步長方案對(duì)上述計(jì)算模型進(jìn)行了模擬，此時(shí)網(wǎng)格劃分采用Mesh C方式進(jìn)行，結(jié)果見表2，表中誤差為采用某一時(shí)間步長計(jì)算結(jié)果與采用最小步長Time D所得計(jì)算結(jié)果之間的差別?？梢姰?dāng)時(shí)間步長從0.010 s降低到0.005 s，總阻力的誤差為0.352%，摩擦阻力誤差為0.022%，說明此時(shí)時(shí)間步長的改變對(duì)計(jì)算結(jié)果的影響不大，故后續(xù)的研究中將會(huì)采用Time C方案。

表2 時(shí)間步長收斂性分析Table 2 Time step convergence analysis

2.3 計(jì)算模型有效性驗(yàn)證

關(guān)于分層流中潛體受力特性的CFD研究較少，本文基于STARCCM+進(jìn)行仿真，對(duì)數(shù)值模型進(jìn)行驗(yàn)證。根據(jù)文獻(xiàn)[15]中的實(shí)驗(yàn)裝置參數(shù)進(jìn)行建模。其中水深為0.6 m，上層淡水密度ρ1=997 kg/m3；下層是密度ρ2=1 024 kg/m3的鹽水；上層淡水深h1為0.3 m，下層鹽水深h2為0.3 m。實(shí)驗(yàn)?zāi)Ｐ蜑?0 cm×45 cm×35 cm的方型浮箱，吃水為0.2 m，速度分別為0.14、0.16、0.18、0.2、0.24 m/s。

數(shù)值研究結(jié)果與試驗(yàn)結(jié)果[15]的對(duì)比見圖3，其中數(shù)值模擬計(jì)算結(jié)果為箱體在雙層流中拖航時(shí)所受的總阻力。由圖3可以看出，數(shù)值模擬結(jié)果與文獻(xiàn)中的試驗(yàn)數(shù)據(jù)吻合較好，說明本文所采用的CFD計(jì)算模型對(duì)潛體受力的模擬具有較好的精度。

圖3 密度分層流數(shù)據(jù)對(duì)比Fig.3 Comparison of the data in density stratified fluid

3 模擬結(jié)果展示與參數(shù)影響分析

3.1 航速變化

3.1.1 航速變化對(duì)潛體所受阻力的影響

潛體位于淡水層中，潛體與交界面的距離會(huì)對(duì)其水動(dòng)力性能帶來較大影響[16]。為系統(tǒng)地研究潛體所處位置對(duì)其阻力曲線的影響，本文模擬采用3個(gè)工況來展示與各交界面距離不同時(shí)潛體的受力特征，見表3，各算例的Fr均取為0.2、0.3、0.4、0.5、0.6、0.7、0.8、0.9、1.0。

表3 各工況參數(shù)設(shè)置Table 3 Parameters of each conditions

設(shè)計(jì)工況A、B對(duì)比來觀察淡水層厚度對(duì)阻力曲線的變化影響；工況B、C的對(duì)比可以觀察潛體在與內(nèi)波面距離不變情況下，與自由液面距離(即潛深)變化對(duì)阻力曲線變化的影響如圖4所示。

圖4 各工況下阻力隨航速的變化Fig.4 Resistances of submerged bodies with different forward speeds under various conditions

由圖4(a)可知，在3種工況下，隨著航速增加，潛體的總阻力均呈現(xiàn)不斷增加的趨勢。對(duì)于工況A，即淡水層的厚度較小、潛體距離2個(gè)交界面均較近時(shí)，其總阻力較其他工況最大，且其增幅會(huì)在高航速區(qū)(0.7≤Fr≤0.9)變得平緩，當(dāng)速度繼續(xù)加大，阻力會(huì)繼續(xù)加速上升。結(jié)合圖4(b)、(c)分析原因?yàn)樵诿芏确謱恿髦?，?dāng)潛體在與交界面距離較近處航行時(shí)會(huì)興起較大波浪，進(jìn)而引發(fā)較大的興波阻力，興波阻力為剩余阻力的主要成分，當(dāng)速度較低時(shí)，興波阻力與摩擦力均隨速度增加而加大，故總阻力加速上升，而當(dāng)航速超過Fr=0.5時(shí)，興波阻力系數(shù)不斷下降，這會(huì)在一定程度上導(dǎo)致總阻力的增速降緩。而隨著航速的繼續(xù)加大，摩擦阻力的持續(xù)增加又導(dǎo)致了總阻力的持續(xù)上升。

由圖4(b)可知，隨著航速增加，潛體所受摩擦阻力會(huì)不斷上升；由于摩擦阻力主要受航速及流體密度的影響，故3個(gè)工況下摩擦阻力隨速度的變化幾乎一致，阻力曲線重合度很高。由圖4(c)可見，3個(gè)工況下的興波阻力系數(shù)表現(xiàn)出相同的趨勢，即隨著潛體航速的增加，阻力曲線先上升后下降，在Fr=0.5附近達(dá)到峰值，這是由于在式(6)中有一余弦項(xiàng)cos(2πmL/λ)的值是在-1.0～+1.0變動(dòng)，這是由潛體艏、艉的興波干擾情況所決定的，故興波阻力系數(shù)Cw與Fr曲線存在波動(dòng)，當(dāng)Fr=0.5時(shí)潛體艏、艉興波發(fā)生了明顯的不利干擾。

為直觀的展現(xiàn)圖4(c)中現(xiàn)象，圖5為工況B中的潛體在以Fr=0.3，0.5，0.7航行時(shí)在自由液面處的興波?？梢?，當(dāng)Fr=0.5時(shí)，自由液面處的興波具有較高的峰、谷值，說明此時(shí)在潛體艏、艉興波發(fā)生了不利干擾，這與圖4(c)中的分析一致。

圖5 工況B中不同航速下自由液面興波分布Fig.5 Wave distribution on free surface with different forward speeds under condition B

3.1.2 航速變化對(duì)潛體表面壓力的影響

針對(duì)工況A、B、C，分別選擇Fr=0.3，0.5，0.7速度展示潛體表面壓力分布情況。圖6中潛體沿X軸負(fù)方向前進(jìn)；圖6為3種工況下潛體整體壓力隨航速的變化情況，3個(gè)工況中潛體表面的一些拐點(diǎn)處壓力值變化相近且集中，不易區(qū)分，故通過圖7展示潛體以工況B航行時(shí)，其艏、艉及附體周圍存在的一些拐點(diǎn)處(圖6(b)中方框)的壓力隨航速的變化情況以直觀分析。

綜合圖6、圖7可見，不同工況下，潛體表面壓力分布展現(xiàn)出相同的形狀。潛體表面壓力主要可分為靜壓、動(dòng)壓2部分，動(dòng)壓主要由速度變化引起，而靜壓則是指潛體表面單位面積上方流體所具有的靜壓能，主要是由潛深決定。同一航速下，工況A的潛體的表面壓力會(huì)低于工況B，而工況C時(shí)潛體表面壓力最大；分析其原因主要在于，工況A相比B、C，其潛體距離自由液面比較近，是潛深的加大增加了潛體表面的靜壓能從而導(dǎo)致了壓力的增加。

圖6 不同航速下各工況潛體表面壓力分布Fig.6 Surface pressure distribution of submerged bodies with different forward speeds under various conditions

圖7 不同航速下工況B中潛體分段表面壓力分布Fig.7 Subsection pressure distribution of submerged bodies with different forward speeds under condition B

此外，隨著航速的增加，潛體表面壓力逐漸加大，其中潛體艏、艉部以及駕駛艙、尾翼等表面曲率變化較大處壓力隨航速的變化較為劇烈，而潛體主體部分的壓力基本維持穩(wěn)定；分析原因?yàn)樗俣鹊母淖冎饕绊憹擉w表面一些拐點(diǎn)處的流體動(dòng)壓，進(jìn)而使得這些部位的壓力變化較大。

3.2 潛體與上下交界面距離的變化

較單一流體，航行于密度分層流中潛體，涉及到自由液面(上交界面)與內(nèi)波面(下交界面)；本節(jié)研究了潛體與上交界面距離d1(即潛深)及與下交界面的距離d2變化對(duì)其受力特征的影響。

3.2.1 潛體與上下交界面距離對(duì)阻力的影響

工況D、E中相同的設(shè)置參數(shù)有ρ1=997.561 kg/m3，ρ2=1 020 kg/m3，h2=0.67L，F(xiàn)r=0.5。不同的參數(shù)設(shè)置見表4，工況D、E分別用于研究潛體與自由液面距離d1、與內(nèi)波面距離d2的變化對(duì)其水動(dòng)力性能的影響。

表4 與交界面距離工況設(shè)置Table 4 Conditions of the distance from the interfaces

圖8為潛體與上、下界面的距離(d1、d2)改變時(shí)其阻力變化情況，其中各圖實(shí)線為d1改變時(shí)的阻力曲線，虛線為d2改變時(shí)的阻力曲線。

由圖8中d1變化可見，隨著潛深的增加，潛體所受總阻力逐漸下降；而在位置變化時(shí)并沒有改變流體的密度以及潛體的航速，故潛體所受摩擦阻力變化極??；此外，隨著潛深的增加，潛體的興波阻力系數(shù)逐漸下降，結(jié)合圖8(a)中總阻力的下降，說明當(dāng)潛深改變時(shí)，興波阻力為總阻力的主要成分。

圖8 阻力隨d1、d2改變的變化Fig.8 Resistances of submerged bodies with different d1 and d2

由于興波阻力是潛體對(duì)其興起的波浪提供能量而產(chǎn)生的，故隨著潛深的增加，潛體的興波也將逐漸減弱。這也說明對(duì)于額定功率下航行的潛體，下潛越深，其航速越高?？梢灶A(yù)計(jì)，當(dāng)潛深足夠大時(shí)，其阻力將變得特別小。

由圖8中d2變化可見，隨著d2的加大，總阻力出現(xiàn)下降的趨勢，但相比d1改變帶來的阻力變化，d2改變導(dǎo)致的阻力下降幅度很低。與d1改變時(shí)的情況相同，潛體所受的摩擦阻力并未隨潛體與內(nèi)波面距離d2的改變發(fā)生較大的變化。由圖8(c)可見，潛體的興波阻力系數(shù)隨著d2的增加也出現(xiàn)降低的趨勢，興波阻力系數(shù)曲線與總阻力曲線一致的變化說明d2變化時(shí)興波阻力同樣為總阻力的主要成分。與圖8(a)中一致，d2改變帶來的阻力系數(shù)變化幅度低于d1。

同時(shí)發(fā)現(xiàn)，在距離d2達(dá)到0.15L之前，潛體所受的總阻力及興波阻力系數(shù)下降速率遠(yuǎn)大于d2超過0.15L之后，當(dāng)d2超過0.15L后，再繼續(xù)增加d2，總阻力及興波阻力系數(shù)幾乎沒有變化。這說明，當(dāng)潛體與內(nèi)波面之間的距離較小時(shí)，距離的變化對(duì)其阻力的影響比較明顯，而當(dāng)潛體航行位置離內(nèi)波交界面較遠(yuǎn)時(shí)，海洋環(huán)境中的密度分層將幾乎不會(huì)對(duì)潛體阻力產(chǎn)生影響。

綜上可知，潛體與2個(gè)交界面的距離改變對(duì)其阻力特性影響較大，尤其當(dāng)潛體靠近交界面時(shí)，其興波阻力會(huì)變大，進(jìn)而導(dǎo)致總阻力的增加。比較來看，潛體與自由液面的距離d1(即潛深)變化對(duì)阻力所帶來的影響幅度要大于潛體與內(nèi)波面距離d2變化所帶來的影響。

3.2.2 潛體與上下交界面距離潛體表面壓力的影響

圖9為針對(duì)工況D、E研究潛體表面壓力分布情況。由圖9可見，不同位置工況下潛體表面壓力分布形狀基本相同，觀察每一個(gè)工況下潛體的表面壓力分布情況，潛體艏、艉及潛體附體處壓力數(shù)值會(huì)比較高，而其余部分比較低。

由圖9(a)可知，當(dāng)d2固定時(shí)，隨著潛深d1增加，潛體表面壓力整體逐漸增加，這主要是由于潛深增加，作用在潛體表面的靜水壓強(qiáng)增加導(dǎo)致的。與圖6、7中航速增加所帶來的影響不同，潛深的增加并沒有使?jié)擉w艏、艉等拐點(diǎn)處的壓力數(shù)值產(chǎn)生較大的改變，而是使?jié)擉w壓力整體上增加。說明航速主要影響動(dòng)壓，進(jìn)而對(duì)潛體表面拐點(diǎn)處的壓力影響較大；而潛深主要影響靜壓，進(jìn)而影響潛體整體壓力數(shù)值；相比來說，動(dòng)壓帶來的壓力數(shù)值變化比較明顯。

圖9 不同位置下潛體表面壓力分布Fig.9 Pressure distribution on submerged bodies surface at different positions

由圖9(b)可見，當(dāng)d1不變時(shí)，d2的改變對(duì)潛體表面壓力的影響很小。這是由于d1不變使得潛體所受靜壓不變，且航速固定，動(dòng)壓也不會(huì)產(chǎn)生較大的變化。

4 結(jié)論

1)潛體的受力特性隨航速的影響較大。隨著航速的增加，潛體總阻力及摩擦阻力均不斷上升，由于潛體艏、艉興波產(chǎn)生了不利干擾使得在Fr=0.5附近時(shí)潛體的興波阻力系數(shù)最大，此時(shí)潛體會(huì)興起較大的波浪。

2)潛體與上下交界面之間距離的改變同樣會(huì)對(duì)潛體水動(dòng)力性能產(chǎn)生較大的影響。當(dāng)潛體靠近交界面時(shí)，其遭受的阻力逐漸增加，興波阻力為總阻力的主要成分。潛體與自由液面的距離即潛深的變化對(duì)潛體受力的影響要大于其與內(nèi)波面距離變化時(shí)帶來的影響。

3)潛體航速及位置的不同會(huì)對(duì)潛體表面壓力帶來影響。其中速度變化主要通過改變動(dòng)壓進(jìn)而對(duì)潛體艏、艉及附體等表面曲率變化較大處的壓強(qiáng)產(chǎn)生影響；而位置的變化主要通過改變靜壓，進(jìn)而影響潛體整體壓強(qiáng)。

4)本數(shù)值模型具有較好的精度，可為密度分層流中潛體受力研究提供一種有效的分析手段。