浮體上浮姿態的穩定性分析與試驗驗證

趙橋生，方 田，肖冬林，潘廣善，徐萌萌

（1.天津大學，天津 300072；2.中國船舶科學研究中心，江蘇 無錫 214082）

0 引 言

由于有效容積及浮力需求，工程上浮體往往為鈍體。為降低技術復雜程度及成本，也無姿態控制系統。僅通過自身浮力實現上浮，顯然其水動力特性以及主體自身的靜力特性（重心、浮心、重量等參數）成為影響其運動穩定性的主要因素。一旦兩者匹配不當將會使航行體在上浮運動中存在“枯葉”現象，即劇烈搖晃甚至翻轉，將大大影響浮體應急上浮以保障乘員生命力或儀器設備安全性的功能。所以工程上迫切需要為保證浮體上浮姿態穩定提供設計指導原則。盡管大攻角上浮鈍體會產生嚴重分離，非線性影響較為復雜，其穩定性可能因浮體而異，可能沒有一般性原則。但小攻角是大攻角必經之路，故小攻角下的穩定性準則仍然具有非常重要的工程意義，至少應該是大攻角浮體遵循的必要條件。

評價穩定性離不開浮體的水動力特性。在流線型的水下航行體數值計算研究方面的研究較多，如周廣禮[1]基于粘性流場直接數值模擬方法開展了潛艇上浮運動數值模擬研究。Atkins[2]利用數值模擬手段研究了俯仰角在-16°～16°之間三個不同狀態下的垂直面操縱性水動力。Roddy[3]開展了水動力CFD 技術處理其設計問題研究，其目標是為CFD 代碼的驗證提供一個試驗數據庫。Wu（2005）[4]考慮了海底對航行體的影響，研究了水下航行體在接近海底過程中的水動力數值模擬方法。林兆偉等（2016）[5]對潛器斜航拖曳、純升沉等典型操縱運動進行了數值模擬。黃苗苗（2019）[6]基于DFBI 方法結合重疊網格技術，研究了內孤立波作用下水下航行體自由運動特性的數值模擬方法。

當然也可以對特定的浮體，將流場及運動用數值方法耦合求解來揭示它的穩定特性，但這耗時、費力，且無一般性。本文在小攻角假設下，建立了浮體上浮的姿態運動方程，給出解析解，分析姿態收斂的條件，形成了穩定性判據。通過數值計算給出了特定浮體的流體動力，并進行了上浮穩定性水池試驗。試驗結果反映了數值計算流體動力及穩定性判據的工程價值。

1 軸對稱浮體上浮運動特性及穩定性分析

1.1 軸對稱浮體上浮穩定性判據

實際使用浮體幾乎均為軸對稱體，研究其上浮運動及穩定性具有重要工程意義。在工程應用中，研究人員關心浮體姿態甚于浮體上浮的軌跡。事實上，當姿態非常穩定時，軌跡必然比較鉛垂。在此給出了浮體上浮姿態的穩定性判據，如下所示。

浮體上浮運動坐標如圖1 所示，為隨體坐標系xyz-G（右手系），oo′為鉛垂線，原點設在浮體的質心G上，記B為浮心，G為質心，研究xGz垂直平面中浮體上浮穩定性。浮心位于重心之上，有xB＞0，浮力FB=mB·g，mB為浮體排開水的質量。設浮體受擾動，x軸與垂線oo′之間產生夾角記為傾角θ，規定繞y軸順時針旋轉θ為正（圖1中方向），對應旋轉角速度、對應旋轉角加速度、特征長度L，攻角α遵從右手法則，圖1中自上浮速度u到Gx順時針轉為正。

圖1 浮體上浮運動坐標示意圖Fig.1 Schematic diagram of buoyant body’s coordinate

則在xyz-G系中，攻角α=θ。記浮體質量為m，繞y軸慣量為Jy，附加質量為A55，上浮速度為u。M'w、M'q分別為對浮體質心的位置導數和旋轉導數。對于浮體運動，動力學方程為

浮體上浮過程中，小擾動情況下有：sinα≈α，sinθ≈θ，則由方程（1）可得

其中，方程（3）中的A和B分別為

當A＞0，B＞0時，K＞A指數發散；當A＞0，B＜0時，K＜A指數收斂；當A＜0，指數發散。

綜上，按狀態1和狀態2可得穩定分布域，如圖2所示。

圖2 穩定分布域Fig.2 Stable distribution domain

結論：A＞0，B＜0為運動穩定條件，即

1.2 浮體上浮運動特性

圖3 浮體上浮運動坐標示意圖Fig.3 Schematic diagram of buoyant body’s coordinate

亦即

微分方程（9）的解為

微分方程（8）的解為

限于純浮力，往往浮體上浮速度達到平衡狀態也不會很大，在這一區間，阻力系數Cd往往不是常數，所以公式（10）和（11）是近似解，若將實際阻力系數Cd=Cd(u)代入式（8），則可以得到精確解。

2 浮體流體動力數值計算

2.1 浮體實例的流體外形

本文研究的浮體外形為鈍體，該浮體的外形及坐標系如圖4所示，具體參數如表1所示。

圖4 浮體示意圖Fig.4 Schematic diagram of the buoyant body

表1 浮體參數Tab.1 Parameters of buoyant body

2.2 控制方程

不可壓縮粘性流體連續性方程和時均化的N-S方程，如式(12)～(13)所示:

根據Boussinesq假設，雷諾應力可表示為

2.3 湍流模型

按作者及其他研究者經驗，本文對于浮體上浮阻力CFD 計算采用RNGk-ε模型。針對位置水動力計算采用標準k-ε模型,旋轉水動力計算采用SSTk-ε模型。

（1）RNGk-ε模型

而Rε則可表示為

式中，η=Sk/ε，η0=4.38，β=0.012，常數ck=cε=1.393，C1=1.42，C2=1.68。

（2）SSTk-ω模型

在k-ω模型中將湍流渦粘度μt表示成湍流動能k和特殊湍流動能耗散率ω的函數：

a*為低湍流雷諾數修正系數，

Gk、Gω為湍流產生項，Yk、Yω為湍流耗散項，均考慮了低湍流雷諾數的修正。

k和ω分別滿足方程（20）和方程（21）：

擴散系數Γk，Γω以下式定義：

式中，σk，σw分別為湍流動能k和湍流動能耗散率ω的普朗特數。速度及湍流耗散率ω以經驗的形式給出：

式中，u+=u/uτ，y+=ρuτy/μ，uτ為剪切速度，κ=0.418 7，E=9.793。

（3）標準k-ω模型

湍流動能k方程為

湍流耗散率ε方程為

式中，G1ε=1.44，C2ε=1.92，Gμ=0.09，σk=1.0，σε=1.3，Gk和Gb為湍流動能生成項，σk和σε為湍流普朗特數，Sk和Sε為源項。

2.4 離散格式及求解

對流項采用二階迎風格式，擴散項采用中心差分格式，離散的流體運動方程求解采用單獨求解器的求解方法，壓力速度耦合問題以SIMPLE 法解決，離散方程以Gauss-Seidel 迭代方法求解，以代數多重網格技術加速迭代的收斂速度。

2.5 邊界條件

（1）無窮遠邊界條件：在數值計算中以入口條件和出口條件描述。

入口條件：定義入口處的速度為無窮遠處的速度，u=u∞；

出口條件：定義壓力出口處的壓力大小為無窮遠處的壓力值，p=p∞；

（2）物面條件：物面的速度為零，u→=0；

（3）對稱面條件：在對稱面處的法向速度為零，un=0，所有物理量在對稱面法向上的梯度為零。

2.6 浮體上浮阻力預報

本文采用實尺度計算，以避免尺度效應。計算對象外形表面網格劃分如圖5所示，采用H-O型結構化網格。近壁面第一層網格內y＋≈40，雷諾應力及湍流耗散標量的分布是以經驗公式近似處理的。

圖5 計算模型網格劃分示意圖Fig.5 Schematic diagram of computational model meshing

圖6 給出了浮體表面壓力系數分布圖，圖7 給出了阻力隨速度的變化曲線。

圖6 浮體表面壓力系數分布Fig.6 Distribution of surface pressure coefficients of thebuoyant body

圖7 浮體阻力計算結果Fig.7 Calculation results of resistance of the buoyant body

2.7 浮體位置導數及旋轉導數預報

基于本文的數值計算方法，在圖1 的xyz-G坐標系中xGz平面內，分別在攻角α=0°～9°，無量綱角速度q'=-0.15～0.15 范圍內，計算浮體上浮受到的力矩M，對應的無量綱力矩系數M'= 2M/ρu2L3，其中ρ為流體介質質量密度，u為前方來流速度，L為浮體長度。它相對于α及q'的導數即力矩的位置導數和旋轉導數M'w和M'q，其具體數值列于表2。圖8 為α=1°條件下的浮體繞流流線圖，圖9為q'=0.1時浮體的瞬時流線圖。

表2 水動力導數計算結果Tab.2 Calculation results of hydrodynamic derivatives

圖8 浮體流線圖（攻角1°）Fig.8 Streamline of the buoyant body at attack angle of 1°

圖9 浮體旋轉運動瞬時流線圖Fig.9 Streamline of the buoyant body with rotating rate at certain instant

3 浮體上浮運動及穩定性判定

3.1 上浮極限速度計算

按1.2 節的上浮速度計算方法和圖7 給出的阻力R(u)，得到該浮體上浮極限速度為1.38 m/s，圖10 給出了浮體的上浮速度與時間的關系曲線。

圖10 上浮速度隨時間變化曲線Fig.10 Buoying velocity versus time

3.2 浮體上浮姿態穩定性預估

前文1.1 節中判據式（7b）的物理本質是，若M'w＞0，即浮體在有攻角時會產生發散力矩，則隨著速度增加此力矩越來越大，而浮力產生的力矩為恢復力矩，大小與上浮速度無關，所以恢復力矩應在浮體任何速度下都大于發散力矩。因此上浮姿態穩定性判定僅需判斷當浮體達到最大速度亦即極限速度時能否滿足式（7b）判據。

4 浮體上浮水池試驗驗證

4.1 浮體上浮試驗概況

為了對該全尺度浮體上浮運動的姿態穩定性結果進行試驗驗證，在中國船舶科學研究中心的水池內進行該浮體上浮試驗。水池內裝滿淡水，浮體布置在位于水下一定深度的安裝座上，安裝座可調整初始姿態角。

浮體上浮試驗時，正浮狀態（零度傾角）試驗的模型安裝照片如圖11所示，帶傾角狀態上浮試驗的模型照片如圖12 所示。試驗過程中，通過水下相機拍攝浮體的上浮運動過程及姿態角變化，同時采用壓力傳感器記錄浮體特定位置的壓力變化，通過換算可獲得浮體上浮過程的速度和位移隨時間變化的曲線。

圖11 試驗模型（正?。〧ig.11 Test buoyant body model at vertical postion

圖12 模型試驗初始狀態（傾斜上?。〧ig.12 Test buoyant body model at inclination postion

4.2 上浮極限速度驗證

浮體勻速段上浮速度數值預報結果與試驗結果進行了對比，如表3所示。由于水池內為淡水，對該淡水中的上浮試驗結果進行了修正，修正為對應海水的上浮結果。從表3 中可以得到上浮勻速段的速度為1.33 m/s，計算值與試驗值相比較為一致，相對偏差為3.8%。

表3 上浮極限速度計算與試驗結果比較Tab.3 Comparison of floating velocity between numerical calculation and test results

4.3 上浮姿態穩定性驗證

在水池中對該浮體進行了0°傾角和30°傾角全尺度上浮試驗，通過在深度方向的水下高速攝像機拍攝浮體的運動姿態。上浮試驗過程中的模型照片如圖13所示，浮體上浮出水時的照片如圖14所示。

圖13 上浮模型試驗Fig.13 Buoying of the body in the tank

圖14 模型試驗上浮出水Fig.14 Body floating to the surface

圖15給出了浮體在30°初始傾角工況下，自23 m水深釋放，浮體姿態角隨不同深度變化的試驗數據。從水池試驗上浮運動的姿態角試驗數據和上浮運動錄像可以看出，該浮體在整個上浮運動過程姿態是穩定的。

圖15 上浮過程中不同深度下姿態角的試驗結果Fig.15 Test results of attitude angle at different depths during buoying

5 結 論

對于實際應用浮體，工程上特別注重其上浮姿態的穩定性。本文的研究可得到以下結論：

（1）本文在小攻角狀態下，對任意外形浮體進行了運動穩定性分析，并給出了姿態的穩定性判據。對大攻角情況，該判據至少應是穩定的必要條件。

（2）針對文中特定浮體，數值計算給出了阻力及力矩的位置導數與旋轉導數，從而預報了該浮體的極限上浮速度，給出了姿態穩定的判斷。

（3）理論預報結果得到了水池試驗驗證，表明本文提出的穩定性判據具有工程應用價值。