海洋內波影響水下航行體水動力特性數值模擬

汪超，杜偉，李廣華，杜鵬*，趙森，李卓越，陳效鵬，胡海豹

1 西北工業大學 航海學院，陜西 西安 710072

2 武漢第二船舶設計研究所，湖北 武漢 430205

3 鄭州機電工程研究所，河南 鄭州 450015

0 引 言

地球表面四分之三的海洋長年受到太陽熱輻射的作用，海水的溫度、鹽度和密度在鉛直方向上形成了穩定的層化結構，當外力擾動破壞了這種分層結構時，就可能產生內波[1]。海洋表面的大氣壓力、風場、地震以及水下航行體的運動等均可成為引發內波的擾動源。內孤立波是一種特殊的非線性內波，其非線性和頻散相平衡，通常以單個波包或波包列的形式出現，最顯著的特征是可以維持一定的波形、波速長距離傳播，距離可達上百公里。此外，內孤立波還攜帶有巨大能量，在海洋內部產生大幅垂向波動的同時，會導致密度躍層上下的海水流動呈現剪切狀態，引起海水強烈的幅聚幅散和突發性的強流[2]。

內孤立波在我國附近海域十分活躍，南海東南部呂宋海峽至大陸近海均為內孤立波高發海域[3]。水下航行體遭遇內孤立波時，其運動穩定性、懸停姿態和水動力特性都有可能會發生突變或失控，導致水下航行體運動失穩、操縱失控等[4]。海洋內孤立波已成為水下航行體安全航行中必須考慮的一個災害性環境因素，迫切需要針對內孤立波與水下航行體的強非線性作用問題開展研究工作，分析內孤立波作用下航行體的水動力特性，并給出內孤立波載荷的作用規律。

實驗研究是客觀認識內孤立波與水下結構物相互作用的重要途徑，但由于受到內波水槽尺度、分層流體制備以及測量技術的限制，數值模擬逐漸成為國內外學者研究內孤立波與水下結構物耦合作用問題的重要方法，并且開展了大量的相應研究。陳杰等[5]以兩層流體中內孤立波的KdV-mKdV 理論解為基礎，模擬了內孤立波與運動潛體之間的相互作用，得到了帶速潛體所受的荷載特性，為以后研究內波與運動潛體的相互作用開辟了一條新路徑；關暉等[6]應用雙推板造波法模擬內孤立波數值水槽，采用有限體積自適應半結構多重網格法研究了上凸、下凹兩種內孤立波的載荷效應，并分析了水下航行體在海洋中遭遇內孤立波時的流場形態和受力特性。Gou 等[7]通過時域數值模型計算了水平淹沒圓柱遇到界面孤立波時的垂直運動，結果表明下凹內孤立波接近潛體時，自身的垂直力和位移變化較大。黃苗苗等[8]基于RANS 方法和KdV 理論建立了分層流中的內孤立波與水下航行體相互作用的數值模擬方法，發現內孤立波會使得航行體周圍流場變得異常復雜，從而引起突變的流體動力和運動，需對此時的航行安全性應給予極大重視。Zou 等[9]對海洋內孤立波與水下浮洞（SFT）之間的相互作用進行了數值研究，分析了內孤立波振幅、SFT與密躍層的相對距離、SFT 截面幾何形狀以及兩流體層密度比對水平力、垂向力系數的影響。Liu等[10]研究了分層流中潛艇的阻力和尾跡，結果表明潛艇的下潛深度顯著影響其水動力性能，當其接近內界面時，會產生更大的阻力。

當前國內外相關研究大多針對某一工況下內孤立波與海洋結構物相互作用的載荷特性開展，對于不同內波環境對海洋結構物水動力特性的影響仍缺乏進一步研究。鑒于此，本文擬通過CFD仿真軟件Fluent 模擬分層流中內孤立波與固定、懸浮航行體的相互作用，對不同內波環境下航行體的水動力特性變化進行深入研究，以得到波體耦合的作用規律和影響機理。

1 數學模型

1.1 內孤立波理論

實際海洋密度呈現多層不均勻分布，為簡化求解，通常將海水假設為兩層、均勻分布的流體，忽略密躍層厚度，將分界面看作密度躍層。設上層流體的深度和密度分別為h1和 ρ1，下層流體的深度和密度分別為h2和 ρ2，假設流體無旋且不可壓縮。考慮到上、下層流體密度差異不大，內孤立波在傳播過程中自由面誘導的垂向速度基本可以忽略，故流體域頂部采用了剛蓋假設，底部設為剛性不可滲透的固壁。建立直角坐標系xoz，振幅為a的內孤立波沿ox軸向右傳播，以豎直向上為z軸正方向，水平向右為x軸正方向，如圖1 所示。

圖1 兩層流體之間的內孤立波Fig. 1 Internal solitary wave between two layers of fluid

設內孤立波為永形波，以波速c沿ox軸正方向傳播，Ui(i=1,2)表示第i層流體質點的平均速度，對波面位移ζ(x,t)做如下變換：

則兩層Green-Naghdi 模型可化為如式(2)所示的KdV 方程。

式中：c0，c1和c2分別為線性項系數、非線性項系數和色散項系數。KdV 理論要求非線性和色散之間是平衡的，是描述兩層流體下內孤立波最早的理論模型。其在方程中保留了非線性項，對于大振幅內孤立波來說，非線性與色散性之間不一定平衡。現有的實驗研究表明，內孤立波振幅較小，即a/h<0.1時（其中h表示總水深），采用KdV 方程進行數值模擬與實驗吻合良好，而當a/h>0.1時，數值結果誤差較大。這主要是由于內孤立波的KdV 理論模型不存在理論極限振幅，因而其不適用于大振幅的內孤立波。

為了彌補KdV 理論只能適用于小波幅內孤立波的缺陷，在式(2)中引入一個立方非線性項，得到eKdV 方程：

eKdV 方程適用于波幅不超過極限波幅的內孤立波。而mKdV 方程[11-12]適用于波幅從0 到hˉ變化的情況（hˉ是界面與臨界水平hc之間的距離）。其界面位移ζ(x,t)的解為上式中：cm為波的相速度； λm為特征波長；g為重力加速度。

mKdV 方程模擬內孤立波垂向分布的水平速度最接近實驗結果，適用于各種上、下層厚度比，具有較強的非線性[13]。因本文模擬的尺度接近實驗室尺度，且數值造波方法是利用水平速度入口造波，故采用mKdV 方程作為造波理論的基礎。

1.2 數值方法

計算模型采用美國國防高級研究計劃局（DARPA）的SUBOFF 潛艇標準模型，縮尺比為1:10。縮小后艇體總長為0.435 m，最大回轉直徑為0.051 m，附體包括指揮臺和垂直艉翼，艉翼剖面為NACA 0020 翼型，對稱布置于艇艉部，潛體重心水平位置距離艏部0.203 m。為了提高計算效率，懸浮、自由運動航行體與內孤立波耦合模擬中選擇光體SUBOFF 模型，詳細參數見表1。

表1 SUBOFF 模型幾何參數（縮尺比1:10）Table 1 Geometric parameters of SUBOFF model (Scale 1:10)

數值水槽長H=15 m，總水深h=1 m，設上層流體的深度和密度分別為h1=0.2 m和ρ1=995 kg/m3，下層流體的深度和密度分別為h2=0.8 m 和ρ2=1 023 kg/m3，流體動力黏度為μ=0.001 N·s/m2，重力加速度g=9.81 m/s2。SUBOFF 模型放置在工作區，其艏部距離入口邊界7 m，計算域的后5 m為消波區，計算域示意圖如圖2 所示。

圖2 計算域示意圖Fig. 2 Schematic diagram of the computational domain

重疊網格方法不涉及網格變形及重構，在模擬瞬時大幅運動方面比傳統的變形網格方法更有優勢，而且更容易調整航行體的初始位置，因此，本文采用重疊網格方法模擬內孤立波與潛艇的相互作用。對數值水槽和潛艇模型分別進行結構化網格劃分，在潛艇表面處使用漸變網格，使嵌套網格和背景網格能夠均勻過渡，邊界層網格劃分保證y+<1（如圖3 所示）。消波區網格逐漸粗化，以增強數值耗散。邊界條件設置如下：水槽左側為速度入口，右側為壓力出口，上表面為Symmetry邊界，底部邊界為無滑移壁面。

圖3 重疊網格示意圖Fig. 3 Schematic diagram of the overset grid

通過mKdV 理論方程計算得到內孤立波在傳播過程中誘導的上、下層流體平均水平速度U1和，將速度控制方程通過UDF 加載到入口邊界上進行數值造波。

流場通過求解不可壓黏性流動納維?斯托克斯（N-S）方程組得到，其控制方程為：

式中：U為流場速度； μ為動力黏性系數； ρ為流體密度；p為流場壓強；fs為添加到消波海綿層中的源項。

航行體多自由度運動控制方程為：

式中：m為 航行體的質量；vG，F，IG， ωG以及MG分別為航行體的速度矢量、所受合外力、轉動慣量、轉動角速度和所受慣性矩，下標G 表示重心。

通過求解不可壓縮、有黏N-S 方程得到流場中航行體的受力和力矩，再將其代入多自由度運動方程中，求解獲得航行體的加速度與角加速度，積分得到航行體各角速度與速度分量，代入控制方程后進行下一時間步網格更新與數值計算。如此反復迭代，直至達到設定的收斂精度或給定時間。在計算過程中，流體作用在航行體上的力和力矩使航行體運動，同時航行體壁面也對流場產生影響。

初始時水槽中流場速度均取為零。在入口處施加速度入口邊界條件，采用壓力出口方法處理定常來流出流問題，同時采用海綿層方法對內孤立波進行數值消波。內孤立波形成于兩層流體之間的分界面處，采用流體體積 (volume of fluid,VOF)方法來捕捉兩相界面的變化。數值處理求解器為非定常隱式求解器，控制方程采用有限體積法（finite volume method, FVM）進行離散，對流項以及湍流輸運項采用二階迎風格式離散，時間離散方式采用一階隱式，壓力速度耦合設置為PISO，節點壓力梯度計算采用最小二乘法，動量方程為二階迎風格式，時間步長取0.005 s。

1.3 算例驗證

為驗證后續數值結果的準確性，對照黃文昊等[15]的實驗工作設置圓柱型結構內孤立波載荷模擬算例。數值水池大小為30 m×1 m×1 m (長×寬×高)，水池中豎直放置一直徑D=0.15 m的圓柱，吃水深度l=0.535 m，距離入口LS=7 m，如圖4 所示。上下層流體深度、密度、黏度系數以及數值設置同1.2 節。

圖4 驗證算例示意圖Fig. 4 Schematic diagram of the verification example

處理后的結果與文獻[15]試驗結果的對比如圖5 所示。由圖5 可知，分層水槽中豎直圓柱上內孤立波載荷特性與試驗結果基本吻合，表明利用該數值方法來模擬分層流中內孤立波與水下航行體相互作用是合理可行的。

圖5 模擬結果與試驗結果對比Fig. 5 Comparison of simulations with experimental results

1.4 算例設置

下文將探究不同內波環境對水下固定、懸浮航行體水動力特性的影響，在模擬中改變波幅、潛深、航速和俯仰角。研究波幅對航行體水動力特性的影響時，分別設置波幅η為0.1、0.15 和0.2 m，航行體在流體分界面下方0.1 m 處保持水平狀態，即潛深d=0.1 m、俯仰角α=0?和航速u=0 m/s；研究潛深對航行體水動力特性的影響時，設置航行體的潛深d分別為分界面以下0.02、0.06 和0.12 m，保持波幅η=0.1 m、俯仰角α=0?和航速u=0 m/s；對于固定航行體，增加了俯仰角對其水動力特性影響的算例，分別設置航行體初始俯仰角 α為0°、6°和10°，保持波幅η=0.1 m、潛深d=0.1 m和航速u=0 m/s；而對于懸浮航行體，研究航速對其水動力特性的影響時，在不同時間段內分別給定航行體的運動狀態。0～20 s 是造波階段，航行體處于固定狀態，保持波幅η=0.1 m、潛深d=0.1 m和初始俯仰角α=0?，工況設置示意如圖2 所示；20～30 s航行體加速至0.5 m/s；30 s 后，在航行體將要遇到內孤立波時開啟航行體的縱向、垂向和俯仰運動3 個自由度。

為簡化下文數據分析，這里對各物理量進行無因次變換。

式中：Fx，Fz和My分別為內孤立波作用下航行體所受到的水平力、垂向力和力矩，力矩轉動中心為航行體重心，水平力向右為正方向，垂向力向上為正方向，力矩逆時針為正方向； ρˉ為上下層流體平均密度；A為 航行體最大橫截面積；t為計算時間。

2 內孤立波與固定航行體

2.1 波幅的影響

由圖6(a)可知，航行體所受水平力的趨勢基本一致，均為先向右增大，后向左增大，但隨著內孤立波波幅的增加，航行體所受水平力的最大值不斷增加，航行體受力曲線開始變化的時間也隨之提前，受力周期也不斷變短。其原因在于，當增大內孤立波波幅時，航行體具有的能量增加，波速也會隨之增加，但特征波長會隨之減小；那么，大波幅內孤立波會更快接觸到航行體，在更短時間內通過航行體。由圖6(b)可以看出：航行體所受垂向力的趨勢也基本一致，均為先向下增大后向上增大，隨著內孤立波波幅的增加，航行體所受垂向力的最大值不斷增加，航行體受力曲線開始變化的時間也隨之提前，且所受垂向力的受力周期也不斷變短，其原因與水平力曲線變化的原因一致；波幅越大，垂向受力波動越明顯。由圖6(c)可以看出，航行體所受力矩的整體趨勢基本一致，均為先順時針增大，后逆時針增大，但隨著內孤立波波幅的增加，航行體所受力矩的變化曲線振蕩更加明顯，受力開始時間也隨之提前，且力矩的受力周期也不斷變短。

圖6 不同波幅內孤立波作用下固定航行體水動力載荷數值結果Fig. 6 Numerical results of hydrodynamic loads on fixed vehicle under the action of internal solitary waves with different wave amplitudes

2.2 潛深的影響

由圖7(a)航行體所受水平力變化曲線可以看出，航行體位于流體分界面下0.02 和0.06 m，即d′=0.02和0.06時，其所受水平力的變化趨勢類似，但位于流體分界面下0.12 m（d′=0.12）時，航行體受力在幅值和相位上存在差異，且初始先受到向左的水平力。這主要是因為潛深0.02 和0.06 m的固定航行體穿過內孤立波，內孤立波向右運動使航行體先受到向右的水平作用力，而潛深為0.12 m的航行體沒有直接與內孤立波作用，受到的是內孤立波誘導的流場，下層流體質點運動方向均為向左，所以航行體先受到向左的作用力；潛深為0.02 m 的航行體受到的水平力幅值最大。由圖7(b)垂向力的變化曲線可以看出，航行體在不同潛深處所受到的垂向力的變化趨勢相同，航行體穿越內孤立波時所受到的垂向力大于未穿過內孤立波情況下的垂向力，d′=0.02時航行體所受垂向力的幅值和作用時間均大于d′=0.06時的航行體。由圖7(c)力矩的變化曲線可以看出，航行體在穿越內孤立波和位于波谷下方時，航行體所受到的力矩的變化方向先順時針增大，后逆時針增大。在航行體穿越內孤立波波面的工況下，內孤立波靠近航行體導致流場中出現渦結構，在波面右側流場中存在向下的垂向速度，對航行體艏部產生向下的作用力，從而出現一個小幅度逆時針的力矩，之后航行體穿越內孤立波右波面。然而，由于波面是向下凹的形狀，航行體穿越波面時，其艇體艏部上方應當先接觸到上層流體，因而其艏部上方受到的壓力變小，航行體艏部下方的壓力大于上方壓力，航行體受到順時針的力矩，并且隨著內波向右運動，其所受到的力矩不斷增大；當航行體重心通過波面后，航行體艉部下表面受到的壓力大于其上表面的壓力，航行體艉部受到向上的作用力，因而順時針力矩會出現小幅度的減少；當航行體穿出內孤立波左波面時，其艏部下方先接觸到下層流體，且受到的壓力變大，而大于艏部上方的壓力，仍然受到順時針的力矩作用。之后，內孤立波逐漸遠離航行體，但在波面左側，由于渦結構的存在，流場中存在向上的速度，航行體艉部受到向上的作用力，出現逆時針的力矩。且隨著航行體深度的增加，力矩幅值波動逐漸減小，力矩曲線也更加平滑，說明航行體潛深越接近兩層流體分界面，其力矩變化情況越復雜，存在更多的不確定性，對航行體的航行安全會產生嚴重影響。

圖7 不同潛深處固定航行體的水動力載荷數值結果Fig. 7 Numerical results of hydrodynamic loads on fixed vehicle at different depths

2.3 俯仰角的影響

由圖8(a)可以看出，不同俯仰角下航行體受到內孤立波作用產生的水平力的變化趨勢與周期基本一致，但隨著俯仰角的增加，航行體所受水平力的最大值不斷增加。其原因在于，航行體俯仰角增加后，其水平迎流面積也隨之變大，受到流速及流速方向不同而導致的摩擦力也會隨之增大；而在水平方向上航行體的受力僅由內孤立波的流動引起，所以水平力的最大值隨著俯仰角增加而不斷增大。同時，從水平力的曲線變化還可以看出，當俯仰角不斷增大時，水平力曲線振蕩程度也隨之增強，俯仰角越大，曲線越不平滑，一方面說明二維情況下航行體對流體的阻塞對于流場影響較大，使內孤立波波形不穩定，并在一定程度上也體現出帶有俯仰角的航行體受力情況更加復雜，對航行體航行安全影響更大。由圖8(b)可以看出，航行體受到內孤立波作用產生的垂向力的變化趨勢與周期也幾乎相同，但隨著俯仰角的增加，航行體所受到的負方向垂向力的最大值不斷減小，而正方向的最大值基本保持不變。由圖8(c) 可以看出，隨著航行體俯仰角的增加，航行體所受力矩的趨勢更加明顯，并且負方向力矩的最大值也隨之增大，而正方向力矩則隨之有一定的減小。

圖8 不同俯仰角姿態下固定航行體的水動力載荷數值結果Fig. 8 Numerical results of hydrodynamic load on fixed vehicle under different pitch angles

3 內孤立波與懸浮航行體

3.1 波幅的影響

圖9 所示為不同波幅內孤立波作用下懸浮航行體水動力載荷數值結果。由圖9(a)可以看出，不同波幅內孤立波作用下航行體所受水平作用力的變化趨勢基本相同。由于航行體都是位于分界面下0.1 m 處，在內孤立波誘導的順時針流場作用下，懸浮航行體先受到向左的水平力，這個水平力較小。在航行體逐漸接近內孤立波波谷位置，向左的作用力逐漸減小，到達波谷所在的中垂線附近時，此時水平力為0，此后，航行體背離內孤立波，受到的水平力向右并逐漸減小；且波幅越大，航行體受到的水平力幅值越大。由圖9(b)可知：航行體在內孤立波流場的作用下，其初始狀態就受到一個向下的力，壓著航行體朝下方加速運動；當內孤立波越過航行體后，垂向力變成向上，航行體繼續朝下方做減速運動；波幅越大，航行體受到的垂向力也越大。由圖9(c)可知：航行體的俯仰力矩與垂向力的關系比較大，η′=0.15和η′=0.20結果曲線變化相近，且力矩的幅值要略大于η′=0.10的算例，力矩曲線振蕩較大；航行體先低頭再快速抬頭，內孤立波過去之后又開始低頭。

圖9 不同波幅內孤立波作用下懸浮航行體水動力載荷數值結果Fig. 9 Numerical results of hydrodynamic loads on suspended vehicle under the action of internal solitary waves with different wave amplitudes

3.2 潛深的影響

從圖10 整體水動力載荷數值結果來看，懸浮航行體所受的水平力、垂向力和俯仰力矩變化趨勢一致，僅載荷幅值略有差異。由圖10(a)可知：在下層流體中的航行體受到的水平力先向左，后向右，這是由于渦的速度方向是順時針的；當航行體接近內孤立波波谷時，此時水平力向左逐漸增大，過了波谷之后水平力朝右增大。從圖10(b)可以看出：航行體在豎直方向先受到向下不斷波動的垂向力，當航行體到達內孤立波波谷位置時，垂向力變為0；隨著內孤立波繼續向右行進，越過航行體后短暫時間內垂向力朝上，但是在慣性力作用下，航行體還是低頭朝下運動。從圖10(c)可以發現：航行體初始位于下層流體，其艏部受到流場的作用力低頭向下，因此初始力矩為低頭力矩；隨著波的向前推進，航行體尾部流體向下的速度不斷增加，航行體受到抬頭力矩；內孤立波波谷經過航行體，流體想要把航行體往上推，由于二維的阻塞作用，只能使航行體逆時針低頭。

圖10 不同潛深處懸浮航行體的水動力載荷數值結果Fig. 10 Numerical results of hydrodynamic loads on suspended vehicle at different depths

3.3 航速的影響

由之前的結果可以發現，位于分界面下層的航行體遭遇內孤立波會產生劇烈的晃動并導致“掉深”現象的發生，本案例就是具有初始航速的航行體，本節分析其在內孤立波作用下的受力和運動情況，并討論潛艇發生劇烈晃動和“掉深”的原因。

圖11 所示為運動航行體與內孤立波相互作用過程的變化情況。由圖中可以看出，航行體經過內孤立波波谷中垂線之前，其俯仰角度基本沒有變化，而在經過波谷中垂線后開始低頭向下運動。

圖11 運動航行體與內孤立波的作用Fig. 11 Interaction between moving vehicle and internal solitary wave

動壓表示流體顆粒每單位體積的動能，也可以理解為驅動流體流動所需要的壓力。圖12 所示為36.5 s 時刻動壓場Dp中航行體附近的流線圖。從動壓場來看，航行體艏部上下動壓相差不大，但航行體艉部上下動壓相差較大；將航行體附近的流場放大并繪制流線，可以清楚地看到，內孤立波經過航行體之后，其逆時針轉動的流場會作用到航行體的艉部，航行體受到一個逆時針旋轉的力矩，有低頭的趨勢；此外，內孤立波行進過程中，常伴隨著尾波的產生，尾波會消耗主波一部分能量，且與主波流場旋向相反，此時航行體剛好處于主波和尾波之間，若繼續向左運動，其艏部會受到內孤立波尾渦的“拍擊”作用，航行體不能維持自身平衡，進而會發生“掉深”現象。

圖12 36.5 s 時刻動壓場中航行體附近流線圖Fig. 12 Streamline near the vehicle in the dynamic pressure field at the moment of 36.5 s

4 結 論

本文主要研究了不同內波環境對航行體水動力特性的影響規律。主要結論如下：

1） 對于固定航行體而言，內孤立波波幅越大，航行體所受水平力、垂向力和力矩的峰值越大，航行體受載曲線開始變化的時間也隨之提前，并且受力周期也不斷變短；內孤立波波幅越大，航行體受載的變化曲線振蕩越明顯。航行體潛深越接近兩層流體分界面時，其所受水平力和垂向力峰值越大、力矩變化情況變得更加復雜，存在更多的不確定性，對航行體的安全會產生嚴重威脅。隨著航行體俯仰角的增大，其所受水平力、力矩峰值也隨之增加，相反地，垂向力隨著俯仰角的增加而減小。

2） 對于懸浮航行體而言，內孤立波波幅越大，航行體上的水平力、垂向力峰值也越大，力矩曲線變化振蕩更大。潛深越接近于內孤立波波谷縱深位置，懸浮航行體的水動力載荷峰值越大，受影響的時間越長。

3） 具有初始航速的航行體穿過內孤立波后，會發生“掉深”現象，其原因可能是由于航行體艉部受到內孤立波主波流場向上的作用力，艏部受到尾渦的“拍擊”，導致航行體失去平衡進而發生“掉深”。