劉 源,王玉紅,王赤忠*

(1.陸軍軍事交通學(xué)院 鎮(zhèn)江校區(qū),鎮(zhèn)江 212003) (2.浙江大學(xué) 海洋學(xué)院,舟山 316021)

波-流-體相互作用是船舶與海洋工程中最具挑戰(zhàn)性的課題之一．流的存在會顯著改變波浪參數(shù)，如波高、頻率等．船舶和鉆井平臺等海洋結(jié)構(gòu)物在波-流聯(lián)合作用下,其受力和運(yùn)動發(fā)生顯著變化,甚至?xí)?dǎo)致結(jié)構(gòu)物的破壞從而影響其安全性．因此，對波-流-體相互作用研究在海洋工程領(lǐng)域備受關(guān)注．

文獻(xiàn)[1]中 基于時域線性理論,研究了位于均勻流中半圓柱的二維輻射和繞射問題,分析了均勻流對漂浮物體處的波浪爬高及作用在物體上的水動力的影響．隨后,將此方法拓展到了三維垂直圓柱繞射問題[2]．文獻(xiàn)[3]中也作了類似的研究，研究表明:當(dāng)入射波傳播方向和水流方向相同時,上游處的波浪爬高會顯著增大,同時水流對作用在結(jié)構(gòu)物上的水動力載荷有明顯的影響,特別是對二階平均力的影響尤為顯著．文獻(xiàn)[4-5]中采用時域二階理論,計算了三維垂直圓柱處的波浪爬高及其所受水動力．國內(nèi)學(xué)者也作了類似的工作[6-7]．目前完全非線性理論的波-流和結(jié)構(gòu)物相互作用方法已經(jīng)建立并得到應(yīng)用[8-13]．相比較線性或二階理論而言,基于非線性理論的數(shù)值計算需要更多的計算資源．

以上波-流-體相互作用研究主要是關(guān)于繞射問題,而對于波-流-體輻射問題的研究則很少,主要工作基于線性理論，如文獻(xiàn)[1]．文中基于時域二階理論,采用高階有限元法研究均勻流中半圓柱的輻射問題,將文獻(xiàn)[1]中的位于均勻流中的線性輻射問題擴(kuò)展到二階,計算自由表面波和半圓柱所受到的波浪力,分析水流對一、二階波及作用在半圓柱上的一、二階水動力的影響,同時討論波和水動力在不同流速下的非線性特性．

1 控制方程

定義坐標(biāo)系oxy(圖1),x軸和靜水面重合,y軸垂直向上,在無波動時水深為h．Sf為瞬時自由表面,Sb為半圓柱表面(物面),其單位外方向矢量為N=(Nx,Ny)．水底表面Sbot位置為y=-h,Sc為流體域?的人工截斷邊界．假設(shè)流體為均勻、不可壓縮和無粘性的理想流體,且流體的運(yùn)動為無旋．半圓柱運(yùn)動時流場速度勢φ滿足Laplace方程:

2φ=0 in ?

(1)

若流場內(nèi)存在大小為U0并沿x軸正方向的穩(wěn)態(tài)流,則流場總的速度勢為:

Φ=U0x+φ

(2)

圖1 坐標(biāo)系Fig.1 Coordinate system

Φ在流場內(nèi)也滿足拉普拉斯方程以及邊界條件:

ηt+Φxηx-Φy=0 onSf

(3)

(4)

ΦN=φN+U0NxonSb

(5)

ΦN=0 ony=-h

(6)

以及在Sc上滿足輻射條件．方程中t、g和η分別為時間、重力加速度和波高,下標(biāo)表示對其的導(dǎo)數(shù)．將方程式(2)代入式(3)～(6)并將φ和η進(jìn)行攝動，展開到二階:

φ=φ(0)1+εφ(1)+ε2φ(2)+…

(7)

η=εη(1)+ε2η(2)+…

(8)

式中：ε為小量(一般可取為波陡)；φ(0)為由水流作用引起的擾動勢,而由其造成的波高忽略不計．φ(k)、η(k)(k=1,2)分別表示第k階速度勢和波高．同時,將水平位移X和角位移Θ也展開到二階,在求解給定簡諧運(yùn)動的輻射問題時X和Θ只取到一階,即X=εX(1),Θ=εΘ(1)．將方程式(7～8)代入方程式(1～6)可得各階(k=0, 1, 2)定解問題:

2φ(k)=0 in ?(0)

(9)

(10)

(11)

(12)

(13)

g(0)=-U0nx

(14)

式中，p為壓力,由Bernoulli方程有：

(15)

將式(2)代入式(15)再代入式(14)中并進(jìn)行攝動，展開得:

2 數(shù)值方法與計算分析

混合邊值問題式(9～13)通過有限元法求解:

Kφ(k)=F(k)(k=0,1,2)

(16)

式中：K為有限元系數(shù)矩陣；F(k)為右端向量,其元素表達(dá)式分別為:

式中，Ni為單元插值函數(shù)，文中采用12節(jié)點(diǎn)等參單元．由于該類單元形函數(shù)為三次插值,對其二次求導(dǎo)后仍保持連續(xù),故可以方便而精確地計算出二階導(dǎo)數(shù)，如方程式(10、11、12)中的φxx、φyy和φxy．

由于在時域內(nèi)計算,每一時刻自由面上的波高和速度勢必須知道,文中采用四階Adams-Bashforth格式計算(t+Δt)時步自由表面波高和速度勢．同時,采用人工阻尼區(qū)(圖1)來吸收反射波[14],保證外傳波條件．

對均勻流中半徑為a的半圓柱繞射問題進(jìn)行計算, 水深h=10a．半圓柱做水平或垂直運(yùn)動：

Y=Asinωt

(17)

式中：A為運(yùn)動幅值；ω為運(yùn)動頻率；t為時間．

圖2 一階波幅值隨Ka變化關(guān)系比較Fig.2 Comparisons of amplitudes of firstorder waves versus Ka

圖3給出了Ka=1.0時半圓柱左右端(x=±a)線性(一階)波和二階波時間歷程．3種不同付氏數(shù)Fn=0、0.064和0.128的計算結(jié)果在圖中同時給出．從圖中可以看出,半圓柱左端(x=-a處)線性波和二階波均隨付氏數(shù)增大而增大,而右端(x=a處)的情況正好相反．左端波變化更明顯．這是因?yàn)榘雸A柱左端波向負(fù)x軸方向傳播,而流的速度朝正x軸方向,兩者相反;而在右端,波向正x軸方向傳播,和流的方向相同,即多普勒效應(yīng)．圖4給出了Ka=1.0時兩種付氏數(shù)Fn=0.064和0.128下線性與線性和二階的合成波的對比．兩者在波峰和波谷處均差異較大,這也表明在不同付氏數(shù)下波的非線性效應(yīng)均較明顯．

圖3 Ka=1.0時一階波和二階波時間歷程Fig.3 First-and second-order wave histories at Ka=1.0

圖4 波浪時間歷程Fig.4 Histories of waves

圖5、6給出半圓柱左右兩端線性波和二階波幅值隨無因次波數(shù)Ka的變化關(guān)系．

圖5 一階波幅值隨Ka變化關(guān)系Fig.5 Amplitudes of first order waves versus Ka

圖6 二階波幅值隨Ka變化關(guān)系Fig.6 Amplitudes of second order waves versus Ka

從圖5中可以看出,在同一付氏數(shù)下,左、右兩端線性波幅值隨Ka增大而穩(wěn)定增大;而在同一無因次波數(shù)Ka下,左端線性波幅值隨付氏數(shù)增大而增大,而右端情況相反．而二階波(圖6),其在同一Ka下隨付氏數(shù)變化情況和一階波類似,但在同一付氏數(shù)下隨Ka增大而緩慢增加,隨后近似趨于常數(shù)．

圖7為時間t從4T到8T(T為線性波周期)時間內(nèi)的線性波以及線性和二階合成波的自由表面輪廓．從圖中可以看出,半圓柱左端自由表面波的幅值較大而波長較短,而右端波幅值較小而波長較長,即多普勒效應(yīng)．從圖中可以看出,由于二階波的影響,半圓柱左、右兩端區(qū)域特別是半圓柱附近的波形存在一定的差異．

圖7 A/a=0.1、Ka=2.0、Fn=0.128時從t=4T到8T時間內(nèi)的波形輪廓Fig.7 Wave profiles at A=0.1a, Ka=2.0 and Fn=0.128for t=4T to 8T with time interval Δt=0.1T

圖8給出了線性力(一階力)和二階力幅值在Fn=0、0.064和0.128時隨Ka變化關(guān)系．

圖8 線性和二階力幅值隨Ka變化關(guān)系Fig.8 Amplitudes of first and second orderforce components versus Ka

在Fn=0時,水平方向的線性力和二階力因?yàn)閷ΨQ性為零．在Fn=0.064和0.128時,水平方向的線性力和二階力均隨Ka增大而緩慢增大．這主要是因?yàn)樗鞯拇嬖趯?dǎo)致流場內(nèi)的速度勢不再關(guān)于y軸對稱,進(jìn)而在半圓柱表面上的速度分布也不再左右對稱;而在垂直方向,在同一付氏數(shù)下線性力和二階力隨Ka增加較快．同時,3種不同付氏數(shù)下線性力和二階力差別很?。畧D9為二階平均力在3種不同付氏數(shù)下隨Ka變化關(guān)系．在付氏數(shù)大于零時,在無因次頻率0.4～2.0范圍內(nèi),水平和垂直方向的平均力均為負(fù)值,且隨著付氏數(shù)的增大而變?。?/p>

現(xiàn)在考慮半圓柱做水平運(yùn)動時的情況．運(yùn)動幅值仍取A/a=0.1．圖10給出x=±5a處線性波幅值和文獻(xiàn)[1]中的對比,在兩種不同付氏數(shù)Fn=0.064和0.128下其對比結(jié)果均符合較好．圖11給出了線性力以及線性和二階的合成力的對比．在兩種不同付氏數(shù)Fn=0.064和0.128下水平力幾乎沒有差別,而垂直力,在兩種付氏數(shù)下其形態(tài)雖有所不同(圖12(b,d)),但兩者的非線性特性均非常明顯．圖13為一階力和二階力隨Ka變化關(guān)系．3種不同付氏數(shù)下的水平一階力差異較小并最終趨近于一常數(shù)．而垂直力有明顯的差異．與一階力相反,二階力在水平方向差異明顯,而在垂直方向則幾乎沒有區(qū)別．

圖9 二階平均力隨Ka變化關(guān)系Fig.9 Second order mean force components versus Ka

圖10 比較一階波隨Ka變化關(guān)系Fig.10 Comparisons of amplitudes offirst order waves versus Ka

圖11 Ka=1.0時水動力時間歷程Fig.11 Histories of forces at Ka=1.0

圖12 一階力和二階力幅值隨Ka變化關(guān)系Fig.12 Amplitudes of first-and second-orderforce components versus Ka

平均二階力(圖13)在水平方向變化較復(fù)雜．在Ka<1.3時小于零,其后隨著Ka的增加其值增加為正值;而垂直方向的二階平均力則均為負(fù),而隨著Fn的增大略微減?。?/p>

圖13 二階平均力幅值隨Ka變化關(guān)系Fig.13 Second order mean force components versus Ka

3 結(jié)論

(1) 半圓柱上游區(qū)域(左端)和下游區(qū)域(右端)的一階波和二階波隨付氏數(shù)的增大而增大,而下游波隨付氏數(shù)變化情況則正好相反,而上游波在不同付氏數(shù)下非線性特性均較明顯．

(2) 在做垂直運(yùn)動時,水平力隨付氏數(shù)不同而變化明顯,而垂直力差異很小;在水平運(yùn)動時,垂直線性力和水平二階力隨付氏數(shù)變化而明顯不同,但水平線性力和垂直二階力的差異很?。?/p>