自由落體狀態下二維結構物砰擊載荷計算

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

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

在現如今的船舶與海洋工程領域,對于多體船或者剖面曲率很大的船舶,較大的航行速度會使船體受到很大的砰擊力,這會導致十分嚴重的結構損傷甚至發生海上事故[1].因此,砰擊現象的研究對于船舶和海洋結構物的性能以及強度具有重要的意義.

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

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

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

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

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

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

1.1 理論基礎

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

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

t時刻物體表面任意一點的垂向坐標可表示為

z=f(y)-ξ(t)

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

(1)

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

(2)

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

(3)

-c≤y≤c

1.2 砰擊壓力計算

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

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

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

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

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

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

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

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

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

p(y,t)=

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

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

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

式中:fy為f函數對y的導數.

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

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

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

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

1.3 砰擊力計算

在物體砰擊過程中,作用在物體上的砰擊力是通過沿物體濕表面對其所受的壓力進行積分獲得的.不過,壓力表達式中的速度相關項包含了一個奇異項,1/(c2-y2), 這使得在y=c附近壓力變為負值且不可積分[27].根據Korobkin[25]的建議,這一速度相關項的砰擊壓力應在區域-c*≤y≤c*進行積分,其中c*是距y=c處最近且滿足pv(c*,t)=0的點.

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

(4)

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

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

(5)

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

(6)

(7)

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

表1 不同理論模型下的砰擊力系數

2 時域上的精細積分法

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

(8)

M=m+ρaa1

G=ρav1

K=0

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

v=[ξζ]T

其中ζ表示為

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

(9)

f=[0r]T

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

vk+1=exp(Hη)vk+

(10)

式中:η為時間間隔.

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

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

(11)

通常選擇20作為N的取值,于是應用附加定理之后的時間間隔δ可被視為一個非常小的值,因此前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的增量而導致計算精度的下降,并不是直接對式(13)連續地進行平方,而是使用下式來表征增量的部分:

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

exp(Hη)=I+RN

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

3 結果與討論

3.1 模型驗證

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

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

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

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

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

圖4 船體剖面示意圖

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

3.2 自由落體狀態下的砰擊力特性研究

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

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

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

(14)

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

(15)

(16)

Av=ξCv

(17)

Aa=ξ2Ca

(18)

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

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

(19)

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

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

|ξ|max=

(20)

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

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

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

圖6 自由落體狀態下楔形體的加速度變化

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

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

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

(21)

(22)

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

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

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

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

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

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

4 結論

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

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

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

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