孫 哲, 眭旭鵬, KOROBKIN Alexander 鄧彥增, 張桂勇, 宗 智, 姜宜辰

(1. 大連理工大學(xué) 船舶工程學(xué)院,遼寧 大連 116024; 2. 東安格利亞大學(xué) 數(shù)學(xué)學(xué)院,英國 諾維奇 NR47TJ; 3. 中國空氣動力研究與發(fā)展中心 計算空氣動力研究所,四川 綿陽 621000)

在現(xiàn)如今的船舶與海洋工程領(lǐng)域,對于多體船或者剖面曲率很大的船舶,較大的航行速度會使船體受到很大的砰擊力,這會導(dǎo)致十分嚴(yán)重的結(jié)構(gòu)損傷甚至發(fā)生海上事故[1].因此,砰擊現(xiàn)象的研究對于船舶和海洋結(jié)構(gòu)物的性能以及強度具有重要的意義.

砰擊入水的開創(chuàng)性工作可以追溯到von Karman[2]和Wagner[3]對于楔形體入水過程的研究.基于經(jīng)典Wagner模型,眾多學(xué)者對其進行了不同程度的發(fā)展和完善.Dobrovol’skaya[4]研究了具有自由液面流體的相似流動,并提出了入水砰擊問題的數(shù)學(xué)積分表達式.在Zhao等[5]的模型中,自由液面是線性的而物體邊界和伯努利方程保持非線性.Cointe等[6]、Howison等[7]和Oliver[8]通過匹配漸近展開法將整個流場分成3部分并進行了深入的研究.段文洋等[9]提出了通過輔助線來構(gòu)建在流動分離發(fā)生后的虛擬物體表面,以便快速地預(yù)報砰擊力.除了對于楔形體的一系列研究,Tassin等[10]對拋物線物體的入水模型開展研究,而Korobkin[11-12]則研究了具有任意形狀物體的入水過程.

試驗方法是研究砰擊過程的特征以及驗證解析理論模型的重要手段[13],其中楔形體是被研究最多的模型.典型的研究工作包括Panciroli等[14-15]對剛性和彈性對稱楔形體的砰擊過程進行試驗和計算.近年來,眾多學(xué)者對楔形體的傾斜砰擊入水以及非對稱楔形體的砰擊入水進行了試驗[16-18].對于更復(fù)雜的船體結(jié)構(gòu),比如船艏部分,Aarsnes[19]所進行的自由落體試驗是其后數(shù)值方法驗證[20-21]中引用最多的研究成果之一.

如今,數(shù)值模擬方法以數(shù)據(jù)可靠、計算速度快等特點成為砰擊基礎(chǔ)理論研究和實際應(yīng)用中最常用的工具.邊界元法(BEM)是最早應(yīng)用于砰擊模擬的數(shù)值方法之一,比如Sun[20]、Bao等[22]和Wu等[23]的研究工作.盡管在數(shù)值模擬中可以考慮自由液面的非線性效應(yīng)和流體的黏度,但計算成本過高使其并不適用于海洋結(jié)構(gòu)物設(shè)計的初期.

目前的解析理論研究大多數(shù)集中在勻速入水或控制入水.然而,在這種情況下物體所受的砰擊力可能會與其在自由落體入水情況下所受的砰擊力大不相同.同時,自由落體入水情況更適用于實際工程領(lǐng)域中的結(jié)構(gòu)強度評估.此外,雖然數(shù)值模擬方法同樣能夠應(yīng)用于結(jié)構(gòu)物自由落體砰擊入水的模擬,但解析理論模型更易于從物理學(xué)的角度分析砰擊力的組成和特性,在時間成本方面解析理論模型相較于數(shù)值模擬方法更具有明顯的優(yōu)勢,同樣也能避免在數(shù)值模擬的過程中進行網(wǎng)格劃分操作或者物理模型選取不當(dāng)帶來的主觀誤差.盡管大多只適用于流動分離之前的砰擊入水階段,但是目前的解析理論模型依舊是工程領(lǐng)域預(yù)報結(jié)構(gòu)物砰擊入水初期階段的運動響應(yīng)以及砰擊力峰值最有效率的工具.

鑒于上述背景,將時域上的精細積分法與砰擊力計算的解析理論模型相結(jié)合,將物體的運動控制方程轉(zhuǎn)換為一組常微分方程組,在提升計算精度和減少計算時間的同時,對二維結(jié)構(gòu)物在砰擊過程中的砰擊載荷進行研究.

1 砰擊入水的解析理論模型

1.1 理論基礎(chǔ)

本文使用的解析理論模型均基于勢流理論,即流體被視為不可壓縮且無黏流體.由于砰擊是一個在極短時間內(nèi)發(fā)生的高速動態(tài)過程,所以可以合理忽略流體的重力對于物體砰擊過程的影響[24-25].同時,流體的表面張力也被忽略.

圖1 二維對稱物體砰擊示意圖

t時刻物體表面任意一點的垂向坐標(biāo)可表示為

z=f(y)-ξ(t)

物體關(guān)于z軸對稱,因此其濕表面在y軸上的投影長度是2c(t),其中c(t)可以由Wagner條件[11, 25]計算得到:

(1)

流體動力學(xué)采用基于平板撞擊假設(shè)的Wagner理論[3]進行求解,計算區(qū)域由物體濕表面在y軸上的投影和未受擾動的自由液面構(gòu)成.流體的速度勢函數(shù)φW可通過如下控制方程和邊界條件求解:

(2)

式(2)的解φW(y,0,t)可表示為

(3)

-c≤y≤c

1.2 砰擊壓力計算

物體沿著濕表面的壓力分布是通過伯努利方程計算的:

式中:φ(y,z,t)為z≤0區(qū)域的速度勢函數(shù);ρ為液體的密度.在基于Wagner理論的不同砰擊解析理論模型中,通過不同的近似方法處理非線性效應(yīng)和濕表面形狀來計算物體所受的壓力.

(1) 經(jīng)典Wagner模型(OWM)[26-27].

在OWM中,僅考慮伯努利方程的線性項,將速度勢函數(shù)式(3)代入,得到壓力表達式:

(2) 基于非線性伯努利方程的Wagner模型(WN)[26].

在WN模型中,考慮了伯努利方程的非線性項,同樣將速度勢函數(shù)式(3)代入,得到壓力表達式:

(3) Logvinovich模型(OLM)[25,28].

在OLM中,在考慮了非線性項的同時,將壓力的計算放置在物體浸沒深度所處的平面上(即z=-ξ(t)), 而不是在未受擾動的自由液面上.同時,將速度勢函數(shù)φ表示為Wagner解φW(y,0,t)在z=0處的泰勒展開式,并將式(2)中的第3個方程作為邊界條件.于是,速度勢函數(shù)φ可以改寫為

物體沿濕表面的OLM壓力計算公式可最終表示為

p(y,t)=

(4) 改進的Logvinovich模型(MLM)[25].

在MLM中,Korobkin[25]提出通過考慮物體的形狀來進一步提高OLM的計算精度,即拋開之前的平板碰撞假設(shè),使壓力的計算在真實的物體濕表面上進行.與OLM的處理方法類似,速度勢函數(shù)φ也表示為φW(y,0,t)在z=0處的泰勒展開式,并將式(2)中的第3個方程作為邊界條件.于是,速度勢函數(shù)φ可以改寫為

物體沿濕表面的MLM壓力計算公式可最終表示為

式中:fy為f函數(shù)對y的導(dǎo)數(shù).

(5) 廣義Wagner模型(GWM)[25].

在GWM中,速度勢和壓力的計算放置在物體和流體交界高度附近的物體濕表面上,其余對于速度勢φ函數(shù)的處理與MLM一致.于是,速度勢函數(shù)φ可以改寫為

φ(y,t)≈φW(y,0,t)+(f(y)-f(c(t)))×

物體沿濕表面的GWM壓力計算公式可最終表示為

1.3 砰擊力計算

在物體砰擊過程中,作用在物體上的砰擊力是通過沿物體濕表面對其所受的壓力進行積分獲得的.不過,壓力表達式中的速度相關(guān)項包含了一個奇異項,1/(c2-y2), 這使得在y=c附近壓力變?yōu)樨撝登也豢煞e分[27].根據(jù)Korobkin[25]的建議,這一速度相關(guān)項的砰擊壓力應(yīng)在區(qū)域-c*≤y≤c*進行積分,其中c*是距y=c處最近且滿足pv(c*,t)=0的點.

對于具有任意形狀的物體,其形狀函數(shù)可以表示為物面上均勻分布的離散點的線性差值,如下式所示:

(4)

之后,根據(jù)1.2節(jié)所述各個解析理論模型的壓力計算公式和式(4),逐個對分段進行積分,便可得到物體所受的砰擊力.

為了描述物體在砰擊入水過程中的運動狀態(tài),借助牛頓第二定律并對其進行改寫:

(5)

式中:m為物體的質(zhì)量;g為重力加速度,本文取9.81 m/s2;Fv和Fa分別為組成物體所受砰擊力FH的兩部分,Fv為速度v的相關(guān)項,Fa為加速度a的相關(guān)項,

(6)

(7)

pv為物體所受砰擊壓力中的速度相關(guān)項部分;pa為物體所受砰擊壓力中的加速度相關(guān)項部分;av1、av2、aa1和aa2為砰擊力系數(shù).對于不同的解析理論模型,av1、av2、aa1和aa2的表達式各不相同,如表1所示.表中:

表1 不同理論模型下的砰擊力系數(shù)

2 時域上的精細積分法

在計算獲得物體所受的砰擊力Fv和Fa之后,物體的動力學(xué)控制方程,即式(5),可以通過精細積分法[29]進行求解.在精細積分法中,動力學(xué)控制方程可重新表示為

(8)

M=m+ρaa1

G=ρav1

K=0

接著,引入一個新的相空間變量:

v=[ξζ]T

其中ζ表示為

于是,式(8)可進一步表示為

(9)

f=[0r]T

對于式(9)所描述的物體動力學(xué)系統(tǒng),若已知第k個時間步的變量vk,則下一個時間步中的變量vk+1可通過如下所示的Duhamel積分得到:

vk+1=exp(Hη)vk+

(10)

式中:η為時間間隔.

在時域上使用精細積分法求解物體運動特征的關(guān)鍵在于計算指數(shù)矩陣exp(Hη).該計算包含了兩個步驟.首先,使用指數(shù)矩陣的附加定理將矩陣擴展為如下形式:

exp(Hη)=exp(Hδ)m,δ=η/m,m=2N

(11)

通常選擇20作為N的取值,于是應(yīng)用附加定理之后的時間間隔δ可被視為一個非常小的值,因此前4階泰勒展開式可以用來作為exp(Hδ)的表達形式:

exp(Hδ)≈I+R

R=Hδ+(Hδ)2/2!+(Hδ)3/3+(Hδ)4/4!

(12)

將式(12)代入式(11), exp(Hδ)可近似表示為

exp(Hη)=(I+R)2N=(I+2R+R2)2N-1

(13)

為了避免由于I的值遠大于Rn的增量而導(dǎo)致計算精度的下降,并不是直接對式(13)連續(xù)地進行平方,而是使用下式來表征增量的部分:

在重復(fù)N次上述操作后,指數(shù)矩陣exp(Hδ)最終可以表示為

exp(Hη)=I+RN

如式(10)所示,為了得到第k+1個時間步中的變量vk+1,還需要確定f的值,而這可以通過“預(yù)估-校正”的方法獲得.方法的具體細節(jié)詳見文獻[30].

3 結(jié)果與討論

3.1 模型驗證

首先驗證精細積分法在砰擊問題上的時間收斂性.選取5個不同的時間步長Δt并結(jié)合MLM模擬了一個二維楔形體的自由落體砰擊過程.楔形體的半寬為0.11 m、底升角為25°.加速度隨時間變化的預(yù)報結(jié)果如圖2所示,當(dāng)時間步長小于 0.000 1 s 后,所得結(jié)果幾乎沒有差異.因此,對于之后的所有模擬,均使用 0.000 1 s 作為時間步長.

圖2 基于MLM模型的時間步長收斂性驗證

接著,將1.2節(jié)所述的各類解析理論模型與精細積分法相結(jié)合,研究楔形體在砰擊過程中的運動狀態(tài).具體而言,將上文所述形狀的楔形體置于3個不同的高度處,使其進行自由落體運動.在此,定義物體初始入水速度為v0,則這3個不同初始高度對應(yīng)的物體初始入水速度分別為v0= 2.2,3.1,3.6 m/s.在本文中,所研究的物體均為二維物體,因此取楔形體的單位長度質(zhì)量為2.125 kg/m.預(yù)報楔形體在砰擊過程中加速度隨時間的變化,并將所得結(jié)果與文獻[14]中的實驗結(jié)果進行對比,對比結(jié)果如圖3所示.

觀察圖3可以發(fā)現(xiàn),OWM給出了與其他4種模型相比過高的預(yù)測值,這與文獻[14, 25]中所指出的現(xiàn)象相一致.而其他考慮了非線性因素或物體表面影響的模型得到了更加準(zhǔn)確的預(yù)報結(jié)果,尤其是MLM和GWM.當(dāng)v0= 2.2,3.1 m/s時,MLM預(yù)測的加速度峰值與實驗結(jié)果基本一致,而當(dāng)v0= 3.6 m/s時,非線性模型預(yù)測的加速度峰值較試驗結(jié)果高出10%～20%.

為了檢驗該方法是否適用于形狀更加復(fù)雜的二維剖面,選取如圖4所示的船體剖面模型進行驗證.同樣地,選取3個不同的高度使其進行自由落體運動,對應(yīng)的入水速度分別為v0= 0.61,1.48,2.43 m/s.船體剖面的單位長度質(zhì)量為261 kg/m.將不同解析理論模型得到的垂向砰擊力F隨時間變化的預(yù)報結(jié)果與Aarsnes[19]的試驗結(jié)果以及Sun[20]的BEM數(shù)值模擬結(jié)果進行對比,對比結(jié)果如圖5所示.

圖4 船體剖面示意圖

觀察圖5可以發(fā)現(xiàn),雖然在入水初期,解析理論模型與試驗以及BEM得到的結(jié)果吻合良好,但如果濕表面到達船體剖面上部的末端, 即當(dāng)c等于剖面半寬之后,解析理論模型便不再能夠預(yù)報流體和船體剖面的運動變化.然而,砰擊理論指出,砰擊力會在流體發(fā)生流動分離之后立刻達到峰值,這意味著解析理論模型能夠相對可靠地預(yù)報砰擊過程中物體所受的砰擊力峰值的大小.對于不同模型的計算結(jié)果,OWM依然明顯地過度預(yù)測了垂向砰擊力,而其他模型的結(jié)果更接近Aarsnes[19]的試驗結(jié)果和Sun[20]的BEM結(jié)果.總體來說,如同上文根據(jù)楔形體模型得到的推論一致,對于所有的工況,MLM提供了最準(zhǔn)確的預(yù)測.

3.2 自由落體狀態(tài)下的砰擊力特性研究

解析理論模型的優(yōu)點之一是易于研究砰擊過程中主要的影響因素和物理機理.Korobkin[31]通過解析方法研究了橢圓拋物線物體在自由落體狀態(tài)下的砰擊載荷和運動狀態(tài)變化,Scolan等[32]通過研究發(fā)現(xiàn),無論是在二維模型還是三維模型,當(dāng)物體的入水速度恒定,流體能量在砰擊過程中會均勻地傳遞到抬升與飛濺的液體中.在本節(jié)中,研究物體在自由落體狀態(tài)下所受砰擊力的特性,即砰擊力和物體浸沒深度之間的關(guān)系.

其中:系數(shù)Av和Aa是浸沒深度ξ的函數(shù), 即Av=Av(ξ),Aa=Aa(ξ).

對于砰擊水動力在物體所受的流體載荷中占主導(dǎo)地位的情況,即對于輕質(zhì)物體或具有較大入水速度的物體,物體自身的重力可以合理地忽略不計.因此,通過系數(shù)Av和Aa,便可以將自由落體狀態(tài)下物體的動力學(xué)方程式(5)重新表述為

(14)

之后,式(14)可以進一步表示為

(15)

(16)

Av=ξCv

(17)

Aa=ξ2Ca

(18)

式中:Cv和Ca為兩個僅與楔形體底升角有關(guān)的常數(shù).在OWM中,Cv和Ca可通過下式計算得到:

而在MLM中,Cv和Ca的表達式非常復(fù)雜,具體形式可參照Korobkin等[33]的研究成果,這里不再贅述.通過對式(16)的兩側(cè)在時域上進行積分,物體在砰擊入水過程中的速度[31]可表示為

(19)

之后,將式(19)代入式(14)或式(15),便可得到物體的加速度表達式:

進一步地,加速度絕對值的最大值|ξ|max可表示為

|ξ|max=

(20)

觀察式(20)可以發(fā)現(xiàn),在不考慮物體重力的情況下,對于指定底升角的楔形體,其在砰擊入水過程中最大加速度或砰擊力峰值發(fā)生的浸沒深度由質(zhì)量決定,而不受初始入水速度的影響(假設(shè)楔形體在此浸沒深度之前不會發(fā)生流動分離).

選取半寬為0.11 m,底升角分別為β= 15°,25°,35° 的二維楔形體模型作為研究對象,使用MLM計算了其從不同的高度處在自由落體狀態(tài)下的砰擊過程,并預(yù)報了物體加速度隨浸沒深度的變化.在計算過程中,物體自身的重力忽略不計.不同下落高度所對應(yīng)的物體初始入水速度分別為v0= 2.2,3.1,3.6 m/s.同時,針對不同底升角的楔形體模型,選取了兩個不同的單位長度質(zhì)量進行研究,分別為2.125和8.5 kg/m.

計算所得到的結(jié)果如圖6所示,其展示的結(jié)果與式(20)一致,即雖然物體在砰擊過程中最大加速度或砰擊力峰值的數(shù)值會隨初始入水速度改變,但是其所對應(yīng)的浸沒深度卻與初始入水速度無關(guān).此外,在砰擊入水過程中,物體的加速度或受到的砰擊力會隨著底升角的減小或者自身質(zhì)量的增大而減小.

圖6 自由落體狀態(tài)下楔形體的加速度變化

不過隨著物體質(zhì)量的增加或者初始入水速度的減小,物體自身的重力相較于所受到的砰擊水動力不再是小量.此時,忽略物體重力會導(dǎo)致計算結(jié)果存在誤差.若以5%作為誤差的閾值來判斷預(yù)報計算結(jié)果的準(zhǔn)確性,則針對每一個不同底升角的楔形體模型,依據(jù)其質(zhì)量和初始入水速度繪制能夠忽略物體自身質(zhì)量的臨界曲線.在此,物體單位長度質(zhì)量的研究范圍選取為0.05～50 kg/m,各楔形體模型的臨界曲線如圖7所示.從展示的結(jié)果來看,隨著物體質(zhì)量或者底升角的增大,滿足誤差要求的臨界初始入水速度也隨之增大.

圖7 各楔形體模型的臨界曲線

不過,對于任意復(fù)雜形狀的物體,不同于上述研究的楔形體模型,很難再從式(15)得到類似式(19)這樣的簡單表達式,因為式(17)和式(18)不再是必要的.然而,通過對式(16)進行多次分部積分,可以得到:

(21)

(22)

觀察式(22)可以發(fā)現(xiàn),上述對于楔形體的結(jié)論對于任意復(fù)雜形狀的物體仍然適用,即在忽略物體重力的情況下,最大加速度或砰擊力峰值對應(yīng)的浸沒深度僅與物體的質(zhì)量和形狀相關(guān).

選擇半寬均為0.11 m的兩個不同形狀的曲線楔形體作為驗證此結(jié)論的模型,其中一個為凸面形狀,另一個為凹面形狀,模型的具體構(gòu)型如圖8所示.與上述方法類似,在整個砰擊入水的過程中,物體同樣處于自由落體狀態(tài),物體的單位長度質(zhì)量分別取2.125和8.5 kg/m,物體初始入水速度分別為v0= 2.2,3.1,3.6 m/s,選擇MLM進行預(yù)報物體加速度隨浸沒深度的變化,同時不考慮物體自身的重力.

圖8 曲線型楔形體形狀示意圖

計算所得到的結(jié)果如圖9所示,所得到的結(jié)論與前文對于楔形體模型的分析一致,即雖然最大加速度或砰擊力峰值的數(shù)值取決于物體的初始入水速度,但是其發(fā)生的浸沒深度卻與初始入水速度無關(guān).此外,對比圖8中3類不同形狀楔形體所得的結(jié)果可以發(fā)現(xiàn),剖面的肥瘦程度在很大程度上會影響物體加速度(或砰擊力)和加速度峰值(或砰擊力峰值)所對應(yīng)的浸沒深度.剖面越“肥胖”,物體加速度或砰擊力越大,加速度峰值或砰擊力峰值越早來臨.

圖9 自由落體狀態(tài)下曲線型楔形體的加速度變化

對于此類曲線形狀的物體,不同于楔形體模型,增大物體質(zhì)量或者減小初始入水速度并非會使得計算結(jié)果的誤差單調(diào)增加,因此繪制出能夠忽略物體自身質(zhì)量的臨界曲線十分困難.于是,在此僅給出“物體的初始入水速度需要大于1 m/s”作為能忽略物體自身質(zhì)量的近似前提條件.

4 結(jié)論

采用基于勢流理論的砰擊解析理論模型,研究了自由落體條件下任意二維對稱剖面在砰擊過程中的動力學(xué)特性.使用了幾種典型的解析理論模型,從經(jīng)典Wagner模型到更復(fù)雜的MLM和GWM等.

為了求解物體在自由落體狀態(tài)下的運動狀態(tài),將精細積分法與解析模型相結(jié)合,并在時域上進行積分.若二維剖面的形狀為復(fù)雜的曲線形狀,那么該物體的表面則通過在物面上生成均勻離散的點并在點與點之間進行線性插值進行處理.

之后,選取兩種不同的二維模型,即楔形體和船體剖面,將按照上述方法得到的結(jié)果與現(xiàn)有文獻中的試驗和其他數(shù)值模擬結(jié)果進行了對比.對比結(jié)果表明,考慮非線性因素和濕表面形狀的解析理論模型通常能夠得到更為精確的結(jié)果,尤其是MLM.不過解析理論模型固有的局限性使其無法分析流體在發(fā)生流動分離之后的動態(tài)變化,但這些模型仍然可以正確地預(yù)報流動分離發(fā)生前物體所受的砰擊力峰值.

此外,對于在砰擊過程中可以忽略物體重力的情況,即對于輕質(zhì)物體而言,砰擊力占其所受水動力載荷的主要成分,如果保證物體在自由落體入水過程中產(chǎn)生的最大加速度或砰擊力峰值發(fā)生在流體發(fā)生流動分離之前,那么物體在最大加速度或砰擊力峰值發(fā)生時的浸沒深度僅僅由物體的形狀和質(zhì)量有關(guān),而與其初始的入水速度無關(guān).