■成小飛

(1.福建省交通科學(xué)技術(shù)研究所，福州 350004；2.福建省公路、水運工程重點實驗室，福州 350004)

本文對航標結(jié)構(gòu)開展水動力特性分析，采用光滑粒子流體動力學(xué)(SPH)方法進行分析計算。

1 數(shù)值計算方法

1.1 光滑粒子流體動力學(xué)方法

光滑粒子流體動力學(xué)(SPH)方法是一種先進的數(shù)值離散方法，通過積分近似和粒子近似兩個過程將任意的函數(shù)及其空間導(dǎo)數(shù)離散成一系列具有獨立質(zhì)量和體積的粒子間的相互作用。該方法因其完全拉格朗日特性，在處理流體自由液面大變形問題時無需借助自由液面追蹤技術(shù)，且特別適用于具有復(fù)雜結(jié)構(gòu)邊界和移動流固界面的物理問題，在計算流體力學(xué)領(lǐng)域中具有廣闊的應(yīng)用前景。

在SPH方法中，空間內(nèi)任意點處的場函數(shù)可以轉(zhuǎn)換為如下形式：

式中：i為目標點處的粒子；j為粒子i支持域內(nèi)與粒子i相互作用的其它粒子，數(shù)量為N；x為粒子的位置矢量；m 為粒子的質(zhì)量；Wij=W(xi-xj，h)是核函數(shù)，其數(shù)學(xué)表達式為：

其中，q=|xi-xj|/h，h為光滑長度，取粒子間距的1.5倍。

參照式(1)的形式，得到空間內(nèi)任意點處的場函數(shù)空間導(dǎo)數(shù)的表達式：

式中，rij為粒子i和j的距離。

眾所周知，粘性流體運動可以由連續(xù)性方程和粘性Navier-Stokes(N-S)方程來描述，流體控制方程寫作：

式中：p為壓強；v為速度矢量；g為重力加速度。

利用式(3)分別對式(5)和式(6)進行離散，得到SPH形式的流體控制方程：

式中：∏ij為人工粘性。

另外，為閉合以上公式(8)，還需要引入狀態(tài)方程：

式中：B為與數(shù)值流體壓縮率有關(guān)的參數(shù)；ρ0是參照密度。

粒子的運動位置通過式(10)計算：

1.2 邊界條件

在浮體結(jié)構(gòu)水動力學(xué)問題中涉及到自由液面條件和物面條件兩種。其中，由于SPH方法的拉格朗日無網(wǎng)格特性，自由液面條件是自動滿足的。物面條件是通過某種方式對固壁邊界附近的流體施加排斥力，具體實施方法如圖1所示，在固壁邊界上布置兩排動力學(xué)邊界粒子，動力學(xué)邊界粒子參與粒子間連續(xù)性方程和動量方程的計算，但相對固壁邊界的位置不改變。

圖1 動力學(xué)邊界粒子法示意圖

1.3 數(shù)值積分模式

SPH方法將流體控制方程組轉(zhuǎn)化為一個常微分方程組，簡寫作如下形式：

式中：Θ為計算流體密度的標量算式，F(xiàn)為粒子受到的外力矢量。

采用具有二階精度的預(yù)報矯正法對式(11)進行積分。首先在預(yù)報步計算中間時刻粒子的密度(ρn+1/2)、速度(vn+1/2)和位置(xn+1/2)：

由式(12)計算粒子的Θn+1/2和Fn+1/2，再計算校正步粒子的密度、速度和位置：

并由下式得到下一時刻的粒子信息：

數(shù)學(xué)模型中，計算時間步長需分別滿足包含粘性條件的CFL條件和考慮粒子加速度限制的時間步長約束條件。

式中：F為粒子受到的外力矢量；ξ為小于0.5的安全系數(shù)。

得到時間步長的計算公式：

2 數(shù)值波浪水槽

2.1 主動吸收式造波

數(shù)值波浪水槽是在時域內(nèi)模擬非線性波浪與浮體結(jié)構(gòu)相互作用的重要工具。在基于SPH方法所建立的數(shù)值波浪水槽中，采用推板式造波機生產(chǎn)波浪，并且為了減小波浪二次反射帶來的影響，引入了主動吸收式造波技術(shù)。造波板的運動速度按下式計算：

式中：vwm為造波板的運動速度；ω=2π/T為波浪角頻率，T為波周期；η=Hsin(ωt)/2為目標波浪函數(shù)，H為波高；ηwm為造波板運動平衡位置處的波面高度；Q為水力傳遞函數(shù)；D為波浪頻率的函數(shù)。

2.2 阻尼層消波

由于數(shù)值波浪水槽長度的限制，波浪將在水槽末端發(fā)生反射，反射波浪重新進入試驗區(qū)域，對試驗結(jié)果產(chǎn)生干擾，因此在數(shù)值波浪水槽中采用了阻尼層消浪技術(shù)，其原理是通過在流體動量方程中添加阻尼項來減小水質(zhì)點的運動速度。阻尼層內(nèi)的流體動量方程寫作：

其中：μs為衰減函數(shù)，采用線性形式的計算公式：

式中：αs為阻尼層系數(shù)，取為4；xs為阻尼層的起始位置坐標；ls為阻尼層長度，取入射波長的1.5倍。

3 浮體結(jié)構(gòu)數(shù)學(xué)模型

3.1 浮體運動方程

根據(jù)剛體動力學(xué)理論，剛體結(jié)構(gòu)的運動可以分解成質(zhì)心的平移和結(jié)構(gòu)繞質(zhì)心的轉(zhuǎn)動。本次研究的浮標結(jié)構(gòu)同時受到流體力、重力和錨鏈力的作用，其運動方程寫作：

式中：M、I、V和ω分別為浮體的質(zhì)量、轉(zhuǎn)動慣量、線速度和角速度；Fw和Tw分別為流體力及其力矩；Ft和Tt分別為錨鏈力及其力矩。

3.2 流體力計算方法

浮體結(jié)構(gòu)占據(jù)流場空間，通過對流體施加排斥力來阻止流體穿透結(jié)構(gòu)邊壁并向浮體結(jié)構(gòu)內(nèi)部空間聚集的趨勢。浮體結(jié)構(gòu)的流體力計算示意如圖2。根據(jù)粒子間動量方程，得到浮體結(jié)構(gòu)固壁粒子i對流體粒子j的作用力為：

根據(jù)作用力與反作用力原理，流體粒子j對固壁粒子i的反作用力為：

在粒子i的支持域內(nèi)有N個流體粒子對固壁粒子起作用，因此固壁粒子i受到的流體力為：

累加浮體結(jié)構(gòu)所有固壁粒子受到的流體力，得到結(jié)構(gòu)總的流體力：

式中：Xc為浮體結(jié)構(gòu)質(zhì)心位置矢量。

圖2 錨泊浮式結(jié)構(gòu)SPH模型示意圖

3.3 系泊纜繩模擬方法

采用集中質(zhì)量法模擬浮體結(jié)構(gòu)的系泊系統(tǒng)，如圖3所示。集中質(zhì)量法是將連續(xù)的纜繩離散為有限個(Nt+1)質(zhì)點和有限個(Nt)長度相等、質(zhì)量不計的彈簧，纜繩的質(zhì)量集中在每個質(zhì)點上，纜繩長度的變化通過彈簧的伸長和縮短來實現(xiàn)。首先對纜繩的各質(zhì)點建立靜力平衡方程：

式中：Fti和Ft(i-1)分別是第i個質(zhì)點上端彈簧和下端彈簧的拉力；ti和t(i-1)分別是第i個質(zhì)點上端彈簧和下端彈簧與水平方向的夾角；Wti是第i個質(zhì)點在水中的重量。

第1個質(zhì)點的靜力平衡方程為：

式中：FMx和FMz分別為水平和垂向錨定力；Wt1為第1個質(zhì)點在水中的重量。

圖3 集中質(zhì)量法模型示意圖

由式(27)和式(28)得到第i個質(zhì)點上端的水平拉力和垂向拉力：

則第i個質(zhì)點上端的總拉力為：

第i個質(zhì)點的坐標位置為：

式中：ltj為第j段彈簧的長度。

得到第i個質(zhì)點坐標的計算公式：

上述集中質(zhì)量法需要通過迭代計算來實施，迭代流程如圖4所示。首先假設(shè)纜繩的初始水平錨定力FMx和垂向錨定力FMz，由式(32)計算每個質(zhì)點的垂向位置，可知除纜繩錨固點處的第1個質(zhì)點外，纜繩上的所有質(zhì)點均被拉離開水底。對纜繩上第Nt+1個質(zhì)點的垂向位置z(Nt+1)和浮體系纜點的垂向位置zf進行比較：當z(Nt+1)-zf＞Δ(Δ為設(shè)定的誤差范圍)時，纜繩在垂向方向上過度抬起，實際上應(yīng)有部分質(zhì)點處于拖地狀態(tài)，因此將第1個質(zhì)點的垂向錨定力設(shè)為零，并從第2個質(zhì)點開始逐個將其垂向位置設(shè)為零，每將一個質(zhì)點的垂向位置賦為零值后，重新由式(32)計算 z(Nt+1)，直至|z(Nt+1)-zf|≤Δ；當 zf-z(Nt+1)＞Δ時，說明纜繩在垂向方向上的抬起程度不充分，因此通過增大垂向錨定力來增加各質(zhì)點的高度，直至|zf-z(Nt+1)|≤Δ。在迭代得到各質(zhì)點的垂向高度后，由式(32)計算每個質(zhì)點的水平位置，對纜繩上第Nt+1個質(zhì)點的水平位置x(Nt+1)和浮體系纜點的水平位置xf進行比較：當x(Nt+1)-xf＞Δ時，纜繩在水平方向被過度拉伸，因此通過減小水平錨定力來減小各段彈簧的水平拉伸程度，直至|x(Nt+1)-xf|≤Δ；當xf-x(Nt+1)＞Δ時，纜繩在水平方向的拉伸程度不足，因此通過增大水平錨定力來增大各段彈簧的水平拉伸程度，直至|x(Nt+1)-xf|≤Δ。在迭代得到纜繩的水平和垂向錨定力后，通過式(30)計算第(Nt+1)個質(zhì)點的彈簧拉力，即所求的錨鏈力。

圖4 集中質(zhì)量法實施流程

4 數(shù)學(xué)模型驗證

采用2節(jié)所介紹的數(shù)值波浪水槽設(shè)置方法，建立了如圖5所示的數(shù)值波浪水槽。直角坐標系o-xz的定義如圖中所示。水槽高25m，水深20m。在x=0.5m處設(shè)置推板式造波機，造波機采用主動吸收式造波方式運動。在水槽末端設(shè)置阻尼層，阻尼層長度為入射波長的1.5倍。在水槽內(nèi)，沿程依次布置四個浪高儀：3#浪高儀位于浮體結(jié)構(gòu)的中心位置，用于標定入射波浪要素；1#和2#浪高儀分別設(shè)置在距離浮體結(jié)構(gòu)0.75，通過兩點法分離結(jié)構(gòu)的入射、反射波浪；4#浪高儀設(shè)置在距離浮體結(jié)構(gòu)0.75L處，用來測量結(jié)構(gòu)的透射波浪。

圖5 數(shù)值波浪水槽

在數(shù)值波浪水槽，分別測試了波高H=2m、波周期T=4s、6s和7s三種規(guī)則波浪情況。圖6為不同波況下3#浪高儀處測得的波面高度歷時曲線，可以看到，數(shù)值波浪形態(tài)基本平穩(wěn)，波浪高度和波浪周期的計算結(jié)果與預(yù)期基本一致，尤其當波周期較長時，計算效果更好。因此，所建立的數(shù)值波浪水槽具有生成目標波浪的能力，可以用于后續(xù)關(guān)于航標動力響應(yīng)的研究。

圖6 標定點處波面歷時曲線

5 水動力特性分析結(jié)果

在圖5所示的數(shù)值波浪水槽中，基于3節(jié)所介紹的浮體結(jié)構(gòu)建模方法，建立了波浪與浮標相互作用的數(shù)學(xué)模型。應(yīng)用所建立的數(shù)學(xué)模型，計算了不同波況(H=1.0m、1.5m 和 2.0m、T=4s) 和不同錨泊角度 (θ=30°、45°、60°和90°)下浮標的運動姿態(tài)。

圖7給出不同錨泊角度下浮標運動響應(yīng)的比較，分別以浮標質(zhì)心的橫蕩、垂蕩運動以及浮標繞質(zhì)心的橫搖運動表示，圖中各方向上的運動幅值均除以入射波浪振幅A來作無因次化?？梢钥吹剑^泊角度為90°、即采用垂直錨泊時，浮標的運動幅值最小。對于橫蕩和橫搖運動，當錨泊角度為30°時運動幅值也相對較小，當錨泊角度為45°和60°時運動幅值較大。而對于垂蕩運動，錨泊角度為30°和60°時運動幅值基本一致，當錨泊角度為45°時運動幅值最大。此外，隨著入射波高的增大，浮標橫蕩和橫搖的無因次化幅值大幅度減小，無因次化垂蕩幅值微弱減小，說明小波高波浪作用下波能向浮標動能轉(zhuǎn)換的比例更高，浮標運動產(chǎn)生的輻射波能就較大。

圖7 不同錨泊角度下浮標運動幅值的比較

圖8給出不同錨泊角度下浮標錨鏈力的比較，圖中迎浪側(cè)和背浪側(cè)的錨鏈力幅值均除以gA3來作無因次化?？梢钥吹?，錨泊角度為90°即采用垂直錨泊時，浮標的錨鏈力最大，尤其是背浪側(cè)方向。當錨泊角度為30°時兩側(cè)錨鏈力也相對較大，當錨泊角度為45°和60°時，錨鏈力較小。入射波高對錨鏈力的影響表現(xiàn)為錨鏈力隨入射波高的增大而增大，這與大波高作用下浮標運動響應(yīng)較大的情況相對應(yīng)。

綜上分析，采用90°錨泊角度時浮標的消浪性能最好，航標的穩(wěn)定性也最佳，但由于錨泊系統(tǒng)荷載過大，浮標的安全性受到極大威脅。而采用30°錨泊角度時，雖然浮標的消浪性能略有降低，穩(wěn)定性也稍有下降，但有效緩解了錨泊系統(tǒng)荷載過大的問題，浮標的安全性得到提升。因此，采用30°錨泊角度，航標的穩(wěn)定性和安全性才更有保證。

圖8 不同錨泊角度下浮標錨鏈力的比較

6 結(jié)論

本文基于光滑粒子流體動力學(xué)(SPH)方法對航標結(jié)構(gòu)進行水動力分析，研究不同錨泊角度對航標穩(wěn)定性及安全性的影響。首先建立了航標結(jié)構(gòu)的數(shù)值波浪水槽模型，并進行了驗證，表明其數(shù)值模型可以用來分析計算波浪對航標結(jié)構(gòu)的水動力情況。經(jīng)計算，結(jié)果表明采用30°錨泊角度，航標的穩(wěn)定性和安全性才更有保證。