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

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

0 引 言

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

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

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

1 軸對稱浮體上浮運動特性及穩(wěn)定性分析

1.1 軸對稱浮體上浮穩(wěn)定性判據(jù)

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

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

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

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

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

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

當A＞0，B＞0時，K＞A指數(shù)發(fā)散；當A＞0，B＜0時，K＜A指數(shù)收斂；當A＜0，指數(shù)發(fā)散。

綜上，按狀態(tài)1和狀態(tài)2可得穩(wěn)定分布域，如圖2所示。

圖2 穩(wěn)定分布域Fig.2 Stable distribution domain

結(jié)論：A＞0，B＜0為運動穩(wěn)定條件，即

1.2 浮體上浮運動特性

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

亦即

微分方程（9）的解為

微分方程（8）的解為

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

2 浮體流體動力數(shù)值計算

2.1 浮體實例的流體外形

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

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

表1 浮體參數(shù)Tab.1 Parameters of buoyant body

2.2 控制方程

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

根據(jù)Boussinesq假設，雷諾應力可表示為

2.3 湍流模型

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

（1）RNGk-ε模型

而Rε則可表示為

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

（2）SSTk-ω模型

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

a*為低湍流雷諾數(shù)修正系數(shù)，

Gk、Gω為湍流產(chǎn)生項，Yk、Yω為湍流耗散項，均考慮了低湍流雷諾數(shù)的修正。

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

擴散系數(shù)Γk，Γω以下式定義：

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

式中，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和σε為湍流普朗特數(shù)，Sk和Sε為源項。

2.4 離散格式及求解

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

2.5 邊界條件

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

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

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

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

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

2.6 浮體上浮阻力預報

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

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

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

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

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

2.7 浮體位置導數(shù)及旋轉(zhuǎn)導數(shù)預報

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

表2 水動力導數(shù)計算結(jié)果Tab.2 Calculation results of hydrodynamic derivatives

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

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

3 浮體上浮運動及穩(wěn)定性判定

3.1 上浮極限速度計算

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

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

3.2 浮體上浮姿態(tài)穩(wěn)定性預估

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

4 浮體上浮水池試驗驗證

4.1 浮體上浮試驗概況

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

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

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

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

4.2 上浮極限速度驗證

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

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

4.3 上浮姿態(tài)穩(wěn)定性驗證

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

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

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

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

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

5 結(jié) 論

對于實際應用浮體，工程上特別注重其上浮姿態(tài)的穩(wěn)定性。本文的研究可得到以下結(jié)論：

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

（2）針對文中特定浮體，數(shù)值計算給出了阻力及力矩的位置導數(shù)與旋轉(zhuǎn)導數(shù)，從而預報了該浮體的極限上浮速度，給出了姿態(tài)穩(wěn)定的判斷。

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