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

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

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

0 引言

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

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

1 仿真基礎(chǔ)條件

1.1 模型參數(shù)

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

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

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

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

1.2 計(jì)算域處理

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

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

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

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

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

2 計(jì)算原理

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

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

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

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

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

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

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

本文分別給定15°、10°、5°的初始橫搖、縱搖角度，釋放航行體自由運(yùn)動(dòng)，計(jì)算至適當(dāng)時(shí)長(zhǎng)，計(jì)算結(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é)果只有橫搖與縱搖的自由衰減曲線，因此本文對(duì)附加質(zhì)量與阻尼系數(shù)的計(jì)算過(guò)程進(jìn)行說(shuō)明。

3.1 附加質(zhì)量計(jì)算

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

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

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

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

由該運(yùn)動(dòng)方程求解得出

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

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

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

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

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

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

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

3.2 阻尼系數(shù)計(jì)算

由于實(shí)際上水是有粘性的，航行體會(huì)受到水的非線性阻尼作用，因此以橫搖為例，航行體在靜水中的非線性橫搖運(yùn)動(dòng)方程為[13]

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)束語(yǔ)

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