海冰與柱狀結構作用的近場動力學數值模擬

賈賓， 王慶， 李文頡， 王鍵偉

(1.哈爾濱工程大學 船舶工程學院，黑龍江 哈爾濱 150001; 2.中國船舶重工集團公司 第七〇三研究所, 黑龍江 哈爾濱 150078)

在渤海等覆冰海域，海冰會受風力、潮流力等作用發生漂移，在與柱狀結構作用的過程中，海冰可能出現劈裂、擠壓以及屈曲等破壞模式，其中擠壓破壞最為常見。同時，此作用過程中極易形成交變載荷，引發結構振動。當前對于這種振動的機理解釋大致可分為2類觀點：1)認為冰激振動是一種強迫振動，實質是一種共振現象，與冰的破壞長度有關[1]；2)認為冰激振動是一種自激振動，實質在于能量供給與耗散之間的平衡關系。典型的自激振動模型是Maattanen[2]提出的負阻尼模型。

海冰與海洋結構相互作用的研究主要有理論方法、試驗方法和數值模擬3種方法。數值方法憑借其應用靈活、成本低廉且效率高等優勢，應用逐漸趨于廣泛。其中，在有限元(FEM)方法應用領域，Martonen等[3]應用ANSYS對層冰與剛性錐體結構的作用進行了非線性有限元模擬，其中冰的破壞采用多面失效準則判斷；Hilding等[4]應用內聚力單元法對冰與直立結構的相互作用進行了研究；武文華等[5]應用LS-DYNA計算了海冰與JZ20-2NW北高平臺的相互作用過程中的動冰力和平臺振動響應；龔榆峰等[6]應用ABAQUS的UEL平臺開發出模擬失效的海冰單元，通過冰排與斜坡、圓錐結構相互作用算例進行了驗證。而在離散元(DEM)方法應用領域，Paavilainen等[7]運用二維DEM方法模擬了冰與斜板作用下破碎、堆積現象，研究了作用過程中冰載荷和碎塊形狀特點；Lau[8]運用三維離散元方法模擬了冰排與錐體橋墩的相互作用，并同模型試驗結果進行了對比驗證；狄少丞等[9]基于GPU并行算法，根據梁單元理論構建海冰粘結模型，模擬了海冰與直立、錐體等平臺結構的相互作用；邵帥等[10]采用DEM-FEM耦合算法模擬冰-結構相互作用，采用DEM模擬海冰的破碎，同時運用FEM計算平臺結構的振動響應。

近場動力學(peridynamics，PD)理論由美國的Silling教授提出，并給出了相應的數值算法[11-12]。與經典連續介質力學不同，近場動力學基于非局部作用思想，采用積分型控制方程進行求解。其避免了傳統數值方法中連續性假設和求解偏微分方程的困難[13]，在求解材料的裂紋擴展等非連續性問題上優勢較為明顯，在混凝土、復合材料等破壞問題的分析上取得了良好的效果[13]。

本文采用鍵型近場動力學方法(bond-based perdynamics)，模擬了海冰與柱狀結構的相互作用過程。通過建立PD模型，同時考慮結構振動的影響，將單自由度振動方程加入到結構的位移計算中，求解相互作用過程中冰載荷和柱狀結構的振動響應，將結果同文獻[14]中離散元方法的計算結果進行對比以驗證其合理性。討論分析了樁徑和冰速等參數對冰載荷和結構振動響應的影響。

1 近場動力學理論

如圖1所示，對于某一時刻t下占據空間區域R的物質體，其動力學控制方程為：

(x′-x)dVx′+b(x,t)

(1)

式中：ρ為材料密度；f為物質點x與物質點x′之間的“力密度函數”；u代表物質點的位移矩陣；Hx表示以物質點x為中心且鄰域半徑δ劃定的近場范圍；b為施加在物質點x上的體積力。

彈脆性材料的作用力密度f(η,ξ)可表示為[15]：

(2)

式中：ξ、η分別代表物質點對的初始相對位置和相對位移。

圖1 近場動力學理論相互作用示意Fig.1 Schematic diagram of interaction in peridynamic frame

ξ=x-x′

(3)

η=u(x′,t)-u(x,t)

(4)

c為微彈性模量。對于三維問題的數值模擬[12]，有：

(5)

式中：K為體積模量；s為“鍵”伸長率[12]，為：

(6)

對于各向同性材料，當物質點對的“鍵”處于拉伸狀態時，s>0；反之若為壓縮狀態，s<0。

而μ(t,ξ)為一用來表示“鍵”狀態的標量。

(7)

式中：s0為臨界伸長率。當“鍵”伸長率達到s0時，認為“鍵”發生斷裂。反之則沒有發生破壞。

此外，Silling還提出了描述物質點局部損傷程度的物理量[12]：

(8)

φ取值范圍為0～1。當φ=0時，表明物質點保持完好狀態；當φ=1時，說明該物質點與其近場范圍內的所有物質點之間均不再存在相互作用。

在數值計算時，需將連續體在笛卡爾坐標系下離散為空間內的一系列物質點，同時將求解空間積分方程轉化為求解有限和的形式。其基本方程為：

(9)

2 數值模型建立

2.1 海冰PD材料模型

海冰的力學性質易受溫度、孔隙度、鹽度以及加載速率等影響。如圖2所示，應變率對海冰的壓縮強度影響顯著[16]。在海冰與柱狀結構的相互作用中，當加載速率較低時，冰力呈準靜態特點，海冰體現為韌性材料；當加載速率較高時，冰力呈隨機性特點，海冰材料呈脆性破壞；而在中等加載速率下，冰力則呈穩態特性，海冰材料處于韌-脆轉換狀態[17]?？紤]到本文算例中所選取的冰速，重點考慮海冰的脆性破壞，在計算模型中將海冰視為彈脆性材料，滿足線彈性本構關系。

圖2 海冰材料特性隨應變速率的變化特征Fig.2 Characteristic variety of sea ice with strain rates

此外，海冰在較高應變率下的脆性破壞呈現出較明顯的拉壓異性。相關實驗數據顯示，冰的壓縮強度一般為拉伸強度的3～5倍[18]，此處取sc=-4·st，其中st和sc分別為冰材料臨界拉伸、壓縮率。海冰的PD材料本構關系如圖3所示。

圖3 海冰材料本構關系Fig.3 Constitutive relation of sea ice

海冰發生韌-脆轉換時對應的應變率約在1.0×10～3.0×10-3s-1[19]。根據文獻[20-23]中渤海冰的相關實驗數據統計發現，海冰在韌-脆轉換范圍內壓縮強度極值近似為脆性區壓縮強度極值的1.3～2倍。本文取韌-脆轉換狀態下海冰的臨界壓縮率為脆性狀態下的1.8倍。

2.2 臨界伸長率s0

近場動力學方法通過s0判定材料是否發生斷裂，可由其與斷裂面的能量釋放率G0關系進行推導。

(10)

式中：E為彈性模量；G0可通過線彈性斷裂力學進行推導[24]：

(11)

式中：KIC為斷裂韌度，可通過實驗測量等方法確定，以此建立s0與KIC之間的關系：

(12)

季順迎等基于渤海海冰斷裂韌度實驗結果，歸納出斷裂韌度KIC與溫度T之間的擬合關系[25]：

KIC=-9.293T+76.366

(13)

2.3 接觸模型

根據Parks等所述，計算物體間的接觸力時需引入短程排斥力，以避免不同物體的物質點占據相同空間位置[15]。當兩粒子間距小于短程排斥力臨界范圍dp時，粒子間相互排斥力密度函數fs為：

(14)

式中：p和i分別為來自不同物體上的粒子；cs為短程力常數。依據經驗，

(15)

dp=min{0.9|xp-xi|,1.35(rp+ri)}

(16)

式中rp和ri分別代表粒子p和i的物質點半徑。

2.4 樁腿振動模型

如圖4所示，樁腿被簡化為具有一定質量M、結構剛度K和阻尼系數C的振動結構，受海冰作用產生水平x方向的振動，滿足單自由度振動的動力學方程：

(17)

式中Fc表示柱狀結構與海冰之間的接觸力。

2.5 計算方法

海冰粒子的運動采用向后差分法進行求解。在第n步時，

(18)

(19)

(20)

根據文獻[12]，時間步長需滿足:

(21)

圖4 結構振動模型Fig.4 Vibration model of the structure

圖5 數值模擬計算流程Fig.5 Flow chart of numerical simulation

3 結果與分析

3.1 參數選取

數值模擬主要參數取自文獻[14]，如表1所示。海冰厚度方向劃分3層粒子，粒子尺寸Δx約為0.087 m，取近場范圍δ為3·Δx，時間步長Δt為1.0×10-4s。為體現海冰半無限寬特性并忽略邊界的影響，除與樁腿發生接觸的一面處于自由狀態之外，其他三面均作剛固處理。取樁腿移動的方向為坐標系x軸正向，計算模型如圖7所示。

圖6 “鍵”力計算方法Fig.6 Computation method of the bond force

表1 計算參數Table 1 Calculation parameters

圖7 海冰與柱狀結構相互作用的數值模型Fig.7 Numerical model of the interaction between the sea ice and cylindrical structure

3.2 計算結果

圖8、9分別給出了黃焱等的海冰與柱狀結構相互作用的實驗現象以及本文數值計算得到的效果圖。通過對比發現，近場動力學數值模型能夠較為合理地反映海冰在此過程中的擠壓破碎現象。

圖8 海冰與直立樁腿相互作用的實驗現象[26]Fig.8 Experimental phenomenon of interaction between sea ice and cylindrical structure[26]

圖9 數值模擬可視化效果圖(第21 400步)Fig.9 Virtual reports of the numerical simulation(step 21 400)

0.5 m/s冰速下柱狀結構所受水平冰載荷時程變化如圖10(a)所示。此外，增加采樣頻率，在圖10(b)中給出了4.8～6 s內的冰載荷細節。根據圖像，水平冰力呈無規律的加、卸載隨機性特點。計算得到的柱狀結構振動響應時程如圖11所示。

圖10 水平冰力的時程變化Fig.10 Horizontal ice force during the simulation

圖11 振動響應的時程變化Fig.11 Vibration response of structure during the simulation

表2給出了在相互作用過程中計算得到的水平冰力、柱狀結構振動位移及振動加速度的數據分析結果，同時與文獻[14]中計算得到的結果進行了對比。經分析，通過PD方法計算得到的冰載荷與柱狀結構的振動響應與文獻[14]中的計算結果吻合較好，驗證了計算方法的合理性。

表2 PD計算結果與文獻[14]計算結果的對比Table 2 Comparision of results calculated by Peridynamics and from literature [14]

4 參數影響分析

4.1 樁徑對水平冰力與振動位移的影響

對1.2 m、1.6 m和2.4 m樁徑下海冰與柱狀結構相互作用進行了數值模擬，冰速仍取0.5 m/s，其他計算參數不變，同3.1節所述。水平冰力與結構振動位移計算結果分別如表3～4所示。

表3 不同樁徑作用下水平冰力的計算結果Table 3 Horizontal ice forces with various diameters of cylindrical structure

表4 不同樁徑作用下結構振動位移的計算結果Table 4 Structure vibration displacements with various diameters of cylindrical structure

如圖12(a)所示，水平冰力整體隨樁腿直徑增大呈增加趨勢?？紤]當樁徑較大時，兩者間接觸區域較大，導致了隨機冰力的增大。如圖12(b)所示，對比不同樁腿直徑下結構的振動響應發現，結構振動位移極值和均值均隨樁徑增加而增加?？紤]樁徑增大提升了冰力值總體水平，增大了振動系統外力輸入，進而導致了振動響應的增大。

圖12 樁徑的影響規律Fig.12 The effect of diameter of cylindrical structure

此外，如圖13所示，計算發現Fmax/Fmean隨徑厚比D/h增大逐漸減小，與Sodhi的實驗結果一致，符合Kry提出的冰載荷極值與均值隨著結構寬度增加逐漸接近的觀察結果[26]。

圖13 Fmax/Fmean與D/h的關系Fig.13 Ratio of maximum force to mean force Fmax/Fmean versus the ratio of diameter to ice thickness D/h

4.2 冰速對水平冰力與振動位移的影響

分別取冰速0.2、0.8和1.1 m/s進行計算，樁徑仍取2 m，其他參數保持不變。水平冰力、結構振動位移計算結果分別如表5～6所示。

表5 不同冰速下水平冰力的計算結果Table 5 Horizontal ice forces with various moving velocities of sea ice

表6 不同冰速下結構振動位移的計算結果Table 6 Structure vibration displacements with various moving velocities of sea ice

如圖14(a)所示，水平冰力總體隨冰速增大而增大，但增速逐漸放緩，0.5、0.8和1.1 m/s下冰力均值十分接近。而冰速0.2 m/s時水平冰力極值較0.5 m/s時較高，與Sodhi和Morris的平板冰與剛性立柱結構接觸實驗結果具有一定的相似性[27]。如圖14(b)所示，振動位移均值隨冰速增大而增大，但極值卻隨冰速呈負增長趨勢。通過對比各冰速下結構的振動位移時程曲線發現，0.8和1.1 m/s下結構振動位移曲線同0.5 m/s時類似，均呈隨機性。但0.2 m/s下結構振動位移幅值持續維持在0.2～0.4 mm，且具有較明顯的周期性特點，如圖15所示。因此，認為在0.2 m/s冰速下柱狀結構的振動為穩態振動。

圖14 冰速的影響規律Fig.14 The effect of moving velocities of sea ice

圖15 結構振動位移的時程變化(0.2 m/s)Fig.15 Vibration displacement of cylindrical structure during the simulation (0.2 m/s)

黃焱等[27]的實驗結果顯示，當柱狀結構與海冰相互作用發生穩態振動時，結構振動位移主頻與結構的自振頻率近似。對0.2 m/s冰速下的柱狀結構振動位移時程數據進行快速傅里葉變換(FFT)，計算結果如圖16所示。由圖可知，結構振動位移在頻域上具有6.625 Hz的主頻，與系統自振頻率6.497 Hz相近，驗證了結構發生穩態振動。

圖16 結構振動位移頻域分析結果Fig.16 Analysis results of structure vibration displacements in frequency domain

此外，計算了冰載荷極值與均值的比值Fmax/Fmean隨冰速-厚度比v/h的變化關系，其結果如圖17所示。

圖17 Fmax/Fmean與v/h的關系Fig.17 Ratio of maximum force to mean force Fmax/Fmean versus the ratio of velocity to ice thickness v/h

由圖17可知，Fmax/Fmean隨著v/h增大呈微幅減小的趨勢，這與Sodhi的實驗結果相一致[27]。

5 結論

1)總體來說，海冰與柱狀結構相互作用引發的水平冰載荷與樁徑、冰速呈正相關關系。結構的振動位移隨著樁徑增大而增大，隨著冰速增大而減小。

2)柱狀結構發生穩態振動時，其結構振動位移具有較為明顯的主頻且與系統自振頻率近似。

3)擠壓過程中，冰載荷極值與均值的比值Fmax/Fmean隨著D/h和v/h的增大逐漸減小，其中隨v/h的變化幅度較小。