航行體近水面航行附加質(zhì)量與阻尼的數(shù)值研究

邱緒明，柯貴喜，路 偉，江興隆

（1.中國船舶集團有限公司第七一〇研究所，湖北 宜昌 443003；2.清江創(chuàng)新中心，湖北 武漢 430076）

0 引言

近水面無人有動力自航式航行體是一種先進的海洋工程裝備，研究其運動特性具有很高的工程實用價值。航行體在工作中會頻繁地產(chǎn)生橫搖、縱搖運動，其搖蕩的周期與幅值會影響其工作的穩(wěn)定性，一旦其搖蕩幅值超過標(biāo)準(zhǔn)范圍，航行體很容易出現(xiàn)難以控制的偏航與失穩(wěn)現(xiàn)象[1-3]。其運動的狀態(tài)主要取決于工作環(huán)境條件和自身運動特性，而在描述其在水中運動特性的運動方程中，附加質(zhì)量及阻尼是2個關(guān)鍵參數(shù)[4-7]。

本文針對橫搖、縱搖運動受粘性作用較大的特點對這2種運動狀態(tài)進行研究，運用STAR-CCM+軟件，建立了帶興波自由面的數(shù)值水池。在粘性流環(huán)境下模擬航行體在阻尼作用下的自由橫搖和縱搖的衰減過程，得到航行體在水中橫搖與縱搖的固有周期，根據(jù)公式計算出航行體在空氣中的自由橫搖和縱搖周期。航行體在空氣與水中運動的周期之差是由水的附加質(zhì)量引起的，可以根據(jù)周期之差計算出橫搖和縱搖的附加質(zhì)量，航行體在水中運動橫搖和縱搖的幅值衰減是由水阻尼產(chǎn)生的能量耗散。本文采用峰值處理法計算出航行體橫搖和縱搖的對數(shù)衰減率以及阻尼系數(shù)[8-9]。

1 仿真基礎(chǔ)條件

1.1 模型參數(shù)

選取航行體簡化模型作為仿真對象，與原模比例為1∶1，航行體主尺度參數(shù)如表1所示。

表1 航行體主尺度參數(shù)Table 1 Main parameters of underwater vehicle

航行體的網(wǎng)格模型及其坐標(biāo)軸對應(yīng)方位如圖1所示。

圖1 航行體網(wǎng)格模型Fig.1 Mesh model of underwater vehicle

1.2 計算域處理

將整個數(shù)值水池分為內(nèi)域和外域，內(nèi)域剛好包含整個航行體，外域范圍大致為-3Lx＜x＜3Lx，-9Ly＜y＜9Ly，-2Lz＜z＜2Lz，其中Lx、Ly、Lz分別為航行體自身的長、寬、高尺度。如圖2所示，中間的內(nèi)域為動網(wǎng)格區(qū)域，外部則為靜態(tài)網(wǎng)格區(qū)域，當(dāng)物體運動時，動網(wǎng)格區(qū)域隨其同步運動。

圖2 三維數(shù)值水池空間劃分Fig.2 Spatial division of three-dimensional numerical wave tank

動網(wǎng)格區(qū)域需要采用重疊網(wǎng)格方法。重疊網(wǎng)格由背景網(wǎng)格和部件網(wǎng)格相互重疊組成。各網(wǎng)格區(qū)域在空間上重疊，但不存在連通關(guān)系，各區(qū)域分別獨立采用離散控制方程進行隱式不定常求解，并通過邊界點的插值和數(shù)據(jù)傳遞來實現(xiàn)區(qū)域連通，最終得到整個計算域的流場信息[10-11]。

本文采用流體體積方法（Volume of Fluid，VOF）建立了興波自由面，使用六自由度運動求解器DFBI和重疊網(wǎng)格技術(shù)求解航行體運動，并對波面大致區(qū)域和重疊區(qū)域內(nèi)的網(wǎng)格進行了加密處理，處理后的網(wǎng)格截面如圖3所示。

圖3 數(shù)值水池網(wǎng)格截面Fig.3 Grid section of numerical wave tank

2 計算原理

STAR-CCM+針對不可壓縮粘性流體，在整個計算域內(nèi)運動規(guī)律采用連續(xù)性方程和 N-S方程作為控制方程：

式中：ρ為流體密度；μ為流體黏度；p為靜壓；fi為單位質(zhì)量力；ui、uj為時均化后的速度分量。

本文采用隱式非定常計算，在每個時間步長對計算域內(nèi)流場進行控制方程的求解，然后對航行體表面進行積分以獲取當(dāng)前時刻的航行體受到的橫搖或縱搖力矩[9]。以橫搖為例，取三維右手直角坐標(biāo)系，航行體前進方向為x軸，水平運動方向為y軸，垂直水平面方向為z軸，其表達式為

式中：Fz、Fy分別為對應(yīng)網(wǎng)格位置處流場作用于航行體表面的水動力作用力在z軸與y軸方向的分量；Gy、Gz分別為航行體重心在y、z軸的坐標(biāo)，顯然，當(dāng)坐標(biāo)原點取為航行體重心時，Gy、Gz均等于0。

由航行體所受的橫搖力矩Mt及其橫搖慣性矩與橫搖附加慣性矩Jxx、ΔJxx求解出航行體的橫搖角加速度，并進一步迭代求解出橫搖角速度θ˙t以及橫搖角度θt等運動信息，計算表達式如下：

根據(jù)迭代求解出的下一步運動參數(shù)信息，返回給求解器結(jié)合重疊網(wǎng)格技術(shù)對航行體的運動狀態(tài)進行更新，選取合適的時間步長，在每個時間步長內(nèi)重新以上迭代計算過程，即可準(zhǔn)確地模擬航行體的橫搖自由衰減運動，縱搖自由衰減運動求解過程與之同理。

3 仿真與數(shù)據(jù)處理

本文分別給定15°、10°、5°的初始橫搖、縱搖角度，釋放航行體自由運動，計算至適當(dāng)時長，計算結(jié)果如圖4-9所示。

圖4 初始橫搖角15°自由衰減曲線Fig.4 Free attenuation curve with initial roll angle of 15°

圖5 初始橫搖角10°自由衰減曲線Fig.5 Free attenuation curve with initial roll angle of 10 °

圖6 初始橫搖角5°自由衰減曲線Fig.6 Free attenuation curve with initial roll angle of 5 °

圖7 初始縱搖角15°自由衰減曲線Fig.7 Free attenuation curve with initial pitch angle of 15 °

圖8 初始縱搖角10°自由衰減曲線Fig.8 Free attenuation curve with initial pitch angle of 10 °

圖9 初始縱搖角5°自由衰減曲線Fig.9 Free attenuation curve with initial pitch angle of 5 °

本次CFD仿真目的是得到航行體的橫搖與縱搖的附加質(zhì)量與阻尼系數(shù)，但仿真直接得到的結(jié)果只有橫搖與縱搖的自由衰減曲線，因此本文對附加質(zhì)量與阻尼系數(shù)的計算過程進行說明。

3.1 附加質(zhì)量計算

附加質(zhì)量與能量輻射無關(guān)，它是由物體運動引起的周圍流場中流體的慣性運動所致，當(dāng)航行體作簡諧搖蕩時，在半個加速周期中因物體運動而積聚在流體中的能量（指與輻射能量無關(guān)的那一部分能量）在接下來半個減速周期中又全部還給了物體。

以橫搖為例，航行體自由橫搖運動時，作用于航行體上有慣性力矩和恢復(fù)力矩[12]，其表達式如下：

式中：Mxx、Cxx分別為航行體受到的慣性力矩與恢復(fù)力矩；Jxx、ΔJxx分別為航行體的橫搖慣性矩與橫搖附加慣性矩；D為航行體自身重量；h為航行體的穩(wěn)心高，即航行體的重心與浮心的垂向距離。

根據(jù)動平衡原理，航行體在自由橫搖運動狀態(tài)下任意瞬時慣性力矩與恢復(fù)力矩大小相等，方向相反，則

由該運動方程求解得出

本文中，航行體自身重量D、橫搖慣性矩Jxx、橫搖慣性矩ΔJxx、穩(wěn)心高h(yuǎn)為已知量，固有周期可以通過仿真得到的自由橫搖衰減曲線得出，因而可依據(jù)式9得到航行體的橫搖附加慣性矩。

觀察橫搖仿真結(jié)果圖4-6，可以發(fā)現(xiàn)橫搖周期基本保持不變，隨著橫搖初始角度的減小，橫搖周期僅有不明顯的縮短。

參考橫搖周期受阻尼影響的公式

式中：Tθ1為阻尼影響下的橫搖周期；v為液體的運動粘度。由式10可以看出，橫搖阻尼使橫搖周期稍有增大，但是作用甚小，以至于可以忽略不計。理論上當(dāng)橫搖初始角度減小時，橫搖阻尼大幅減小，因而橫搖周期會有不明顯的縮短。仿真結(jié)果的該項特征符合理論預(yù)期。

航行體縱搖附加慣性矩的計算過程與橫搖同理，而觀察縱搖仿真結(jié)果圖7-9，卻發(fā)現(xiàn)當(dāng)縱搖角度較大時，縱搖周期大幅縮短，且在初始縱搖角度為10°和15°時均出現(xiàn)該現(xiàn)象。

本文結(jié)合模型特征分析，縱搖角度超過一定程度時，航行體有較大部分體積位于水面以上，形成“出水”現(xiàn)象，與其正常工作形態(tài)有較大差異。在“出水”形態(tài)下，航行體縱搖運動所能帶動的流域內(nèi)隨其作慣性運動的流體質(zhì)量大幅減少，即縱搖附加慣性矩ΔJyy大幅減小。根據(jù)式9，縱搖的固有周期Tθ會隨ΔJyy大幅減小而大幅縮短。因此，本文統(tǒng)計縱搖周期時不考慮“出水”狀態(tài)下的數(shù)據(jù)，僅統(tǒng)計縱搖幅值大幅衰減后的縱搖周期。

統(tǒng)計圖4-9的橫搖與縱搖運動周期并取均值，得到航行體的橫搖固有周期為 6.8 s，縱搖固有周期為11.4 s。代入其他參數(shù)，由式9可以得到航行體的橫搖、縱搖附加慣性矩分別為15 873 kg·m2和36 120 kg·m2。

3.2 阻尼系數(shù)計算

由于實際上水是有粘性的，航行體會受到水的非線性阻尼作用，因此以橫搖為例，航行體在靜水中的非線性橫搖運動方程為[13]

式中，B1、B2分別為航行體非線性橫搖運動中的線性、非線性阻尼系數(shù)。該方程適用于較大角度的橫搖運動，而本文中考慮較小的航行體的橫搖角度，因此可以簡化為線性橫搖運動，其運動方程為

式中，B為航行體線性橫搖運動中的線性阻尼系數(shù)，理論上小角度橫搖運動時線性阻尼系數(shù)B保持不變。

對橫搖衰減曲線圖4-6采用峰值處理法，將自由橫搖衰減曲線上的橫搖角相對于平衡位置的幅值按照正負(fù)分為2組橫搖峰值角度θ1、θ2、…、θn，根據(jù)2個相鄰的橫搖峰值角度θi、θi+1可以算出橫搖對數(shù)衰減率

將一組橫搖峰值角度所有相鄰數(shù)值分別代入求出橫搖對數(shù)衰減率，再取平均值，即

由式（14）可以得到航行體橫搖對數(shù)衰減率。

航行體橫搖無因次衰減系數(shù)為

航行體橫搖線性阻尼系數(shù)為

依據(jù)式（16）代入已知量便可得到航行體的橫搖線性阻尼系數(shù)。

分別統(tǒng)計并處理圖 4-6的橫搖運動仿真數(shù)據(jù)，得到初始橫搖角分別為15°、10°、5°時，航行體的橫搖線性阻尼系數(shù)分別為 221.909（N·m）/（（°）/s），201.766（N·m）/（（°）/s），191.374（N·m）/（（°）/s）。可以看出，隨著航行體初始橫搖角的增大，橫搖線性阻尼系數(shù)有一定程度的增大，這與小角度橫搖狀態(tài)下的線性阻尼系數(shù)保持不變的理論預(yù)期不相符。

本文結(jié)合航行體水下運動的特點進行分析，產(chǎn)生這種現(xiàn)象的原因在于航行體橫搖運動的阻尼主要成分是漩渦阻尼，漩渦阻尼的阻尼系數(shù)大小與橫搖幅值一般近似呈平方關(guān)系，所以航行體在小角度橫搖時也不會保持恒定的線性阻尼系數(shù)。

同理可對圖 7-9縱搖運動仿真數(shù)據(jù)采用峰值處理法，得到初始縱搖角分別為15°、10°、5°時，航行體的縱搖線性阻尼系數(shù)分別為 385.666 （N·m）/（（°）/s），372.151 （N·m）/（（°）/s），257.381 （N·m）/（（°）/s）。可以看出，當(dāng)初始縱搖角超過一定程度時，縱搖線性阻尼系數(shù)存在突變，縱搖衰減速度大幅增加。

這種現(xiàn)象可以結(jié)合航行體模型結(jié)構(gòu)特征進行分析。當(dāng)初始縱搖角超過一定程度時，航行體如前述有“出水”現(xiàn)象，此時其縱搖運動會造成水線面的劇烈變化，而這個水線面的變化過程由于要克服水的表面張力，會產(chǎn)生大量的能量耗散，因此會出現(xiàn)縱搖線性阻尼系數(shù)的突變。當(dāng)縱搖幅值沒有超過這個程度，例如圖 9初始縱搖角為 5°時，縱搖幅值衰減特征與橫搖相類似，說明此時縱搖運動阻尼的主要成分也是漩渦阻尼，阻尼系數(shù)大小與縱搖阻尼近似呈平方關(guān)系。

4 結(jié)束語

本文運用STAR-CCM+軟件，在粘性環(huán)境下對航行體受粘性影響較大的橫搖與縱搖 2種運動進行了仿真，并采用重疊網(wǎng)格與VOF波方法保證了仿真結(jié)果的嚴(yán)謹(jǐn)性。根據(jù)粘流仿真結(jié)果計算得到相應(yīng)自由度運動時的附加質(zhì)量與線性阻尼系數(shù)，對比初始角度分別為15°、10°、5°時自由衰減運動的特征，得出了阻尼系數(shù)隨初始角度的變化特性，并分析了縱搖“出水”運動狀態(tài)對附加質(zhì)量與阻尼系數(shù)的影響。