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

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

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

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

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

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

1 近場(chǎng)動(dòng)力學(xué)理論

如圖1所示，對(duì)于某一時(shí)刻t下占據(jù)空間區(qū)域R的物質(zhì)體，其動(dòng)力學(xué)控制方程為：

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

(1)

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

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

(2)

式中：ξ、η分別代表物質(zhì)點(diǎn)對(duì)的初始相對(duì)位置和相對(duì)位移。

圖1 近場(chǎng)動(dòng)力學(xué)理論相互作用示意Fig.1 Schematic diagram of interaction in peridynamic frame

ξ=x-x′

(3)

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

(4)

c為微彈性模量。對(duì)于三維問(wèn)題的數(shù)值模擬[12]，有：

(5)

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

(6)

對(duì)于各向同性材料，當(dāng)物質(zhì)點(diǎn)對(duì)的“鍵”處于拉伸狀態(tài)時(shí)，s>0；反之若為壓縮狀態(tài)，s<0。

而μ(t,ξ)為一用來(lái)表示“鍵”狀態(tài)的標(biāo)量。

(7)

式中：s0為臨界伸長(zhǎng)率。當(dāng)“鍵”伸長(zhǎng)率達(dá)到s0時(shí)，認(rèn)為“鍵”發(fā)生斷裂。反之則沒(méi)有發(fā)生破壞。

此外，Silling還提出了描述物質(zhì)點(diǎn)局部損傷程度的物理量[12]：

(8)

φ取值范圍為0～1。當(dāng)φ=0時(shí)，表明物質(zhì)點(diǎn)保持完好狀態(tài)；當(dāng)φ=1時(shí)，說(shuō)明該物質(zhì)點(diǎn)與其近場(chǎng)范圍內(nèi)的所有物質(zhì)點(diǎn)之間均不再存在相互作用。

在數(shù)值計(jì)算時(shí)，需將連續(xù)體在笛卡爾坐標(biāo)系下離散為空間內(nèi)的一系列物質(zhì)點(diǎn)，同時(shí)將求解空間積分方程轉(zhuǎn)化為求解有限和的形式。其基本方程為：

(9)

2 數(shù)值模型建立

2.1 海冰PD材料模型

海冰的力學(xué)性質(zhì)易受溫度、孔隙度、鹽度以及加載速率等影響。如圖2所示，應(yīng)變率對(duì)海冰的壓縮強(qiáng)度影響顯著[16]。在海冰與柱狀結(jié)構(gòu)的相互作用中，當(dāng)加載速率較低時(shí)，冰力呈準(zhǔn)靜態(tài)特點(diǎn)，海冰體現(xiàn)為韌性材料；當(dāng)加載速率較高時(shí)，冰力呈隨機(jī)性特點(diǎn)，海冰材料呈脆性破壞；而在中等加載速率下，冰力則呈穩(wěn)態(tài)特性，海冰材料處于韌-脆轉(zhuǎn)換狀態(tài)[17]?？紤]到本文算例中所選取的冰速，重點(diǎn)考慮海冰的脆性破壞，在計(jì)算模型中將海冰視為彈脆性材料，滿足線彈性本構(gòu)關(guān)系。

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

此外，海冰在較高應(yīng)變率下的脆性破壞呈現(xiàn)出較明顯的拉壓異性。相關(guān)實(shí)驗(yàn)數(shù)據(jù)顯示，冰的壓縮強(qiáng)度一般為拉伸強(qiáng)度的3～5倍[18]，此處取sc=-4·st，其中st和sc分別為冰材料臨界拉伸、壓縮率。海冰的PD材料本構(gòu)關(guān)系如圖3所示。

圖3 海冰材料本構(gòu)關(guān)系Fig.3 Constitutive relation of sea ice

海冰發(fā)生韌-脆轉(zhuǎn)換時(shí)對(duì)應(yīng)的應(yīng)變率約在1.0×10～3.0×10-3s-1[19]。根據(jù)文獻(xiàn)[20-23]中渤海冰的相關(guān)實(shí)驗(yàn)數(shù)據(jù)統(tǒng)計(jì)發(fā)現(xiàn)，海冰在韌-脆轉(zhuǎn)換范圍內(nèi)壓縮強(qiáng)度極值近似為脆性區(qū)壓縮強(qiáng)度極值的1.3～2倍。本文取韌-脆轉(zhuǎn)換狀態(tài)下海冰的臨界壓縮率為脆性狀態(tài)下的1.8倍。

2.2 臨界伸長(zhǎng)率s0

近場(chǎng)動(dòng)力學(xué)方法通過(guò)s0判定材料是否發(fā)生斷裂，可由其與斷裂面的能量釋放率G0關(guān)系進(jìn)行推導(dǎo)。

(10)

式中：E為彈性模量；G0可通過(guò)線彈性斷裂力學(xué)進(jìn)行推導(dǎo)[24]：

(11)

式中：KIC為斷裂韌度，可通過(guò)實(shí)驗(yàn)測(cè)量等方法確定，以此建立s0與KIC之間的關(guān)系：

(12)

季順迎等基于渤海海冰斷裂韌度實(shí)驗(yàn)結(jié)果，歸納出斷裂韌度KIC與溫度T之間的擬合關(guān)系[25]：

KIC=-9.293T+76.366

(13)

2.3 接觸模型

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

(14)

式中：p和i分別為來(lái)自不同物體上的粒子；cs為短程力常數(shù)。依據(jù)經(jīng)驗(yàn)，

(15)

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

(16)

式中rp和ri分別代表粒子p和i的物質(zhì)點(diǎn)半徑。

2.4 樁腿振動(dòng)模型

如圖4所示，樁腿被簡(jiǎn)化為具有一定質(zhì)量M、結(jié)構(gòu)剛度K和阻尼系數(shù)C的振動(dòng)結(jié)構(gòu)，受海冰作用產(chǎn)生水平x方向的振動(dòng)，滿足單自由度振動(dòng)的動(dòng)力學(xué)方程：

(17)

式中Fc表示柱狀結(jié)構(gòu)與海冰之間的接觸力。

2.5 計(jì)算方法

海冰粒子的運(yùn)動(dòng)采用向后差分法進(jìn)行求解。在第n步時(shí)，

(18)

(19)

(20)

根據(jù)文獻(xiàn)[12]，時(shí)間步長(zhǎng)需滿足:

(21)

圖4 結(jié)構(gòu)振動(dòng)模型Fig.4 Vibration model of the structure

圖5 數(shù)值模擬計(jì)算流程Fig.5 Flow chart of numerical simulation

3 結(jié)果與分析

3.1 參數(shù)選取

數(shù)值模擬主要參數(shù)取自文獻(xiàn)[14]，如表1所示。海冰厚度方向劃分3層粒子，粒子尺寸Δx約為0.087 m，取近場(chǎng)范圍δ為3·Δx，時(shí)間步長(zhǎng)Δt為1.0×10-4s。為體現(xiàn)海冰半無(wú)限寬特性并忽略邊界的影響，除與樁腿發(fā)生接觸的一面處于自由狀態(tài)之外，其他三面均作剛固處理。取樁腿移動(dòng)的方向?yàn)樽鴺?biāo)系x軸正向，計(jì)算模型如圖7所示。

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

表1 計(jì)算參數(shù)Table 1 Calculation parameters

圖7 海冰與柱狀結(jié)構(gòu)相互作用的數(shù)值模型Fig.7 Numerical model of the interaction between the sea ice and cylindrical structure

3.2 計(jì)算結(jié)果

圖8、9分別給出了黃焱等的海冰與柱狀結(jié)構(gòu)相互作用的實(shí)驗(yàn)現(xiàn)象以及本文數(shù)值計(jì)算得到的效果圖。通過(guò)對(duì)比發(fā)現(xiàn)，近場(chǎng)動(dòng)力學(xué)數(shù)值模型能夠較為合理地反映海冰在此過(guò)程中的擠壓破碎現(xiàn)象。

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

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

0.5 m/s冰速下柱狀結(jié)構(gòu)所受水平冰載荷時(shí)程變化如圖10(a)所示。此外，增加采樣頻率，在圖10(b)中給出了4.8～6 s內(nèi)的冰載荷細(xì)節(jié)。根據(jù)圖像，水平冰力呈無(wú)規(guī)律的加、卸載隨機(jī)性特點(diǎn)。計(jì)算得到的柱狀結(jié)構(gòu)振動(dòng)響應(yīng)時(shí)程如圖11所示。

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

圖11 振動(dòng)響應(yīng)的時(shí)程變化Fig.11 Vibration response of structure during the simulation

表2給出了在相互作用過(guò)程中計(jì)算得到的水平冰力、柱狀結(jié)構(gòu)振動(dòng)位移及振動(dòng)加速度的數(shù)據(jù)分析結(jié)果，同時(shí)與文獻(xiàn)[14]中計(jì)算得到的結(jié)果進(jìn)行了對(duì)比。經(jīng)分析，通過(guò)PD方法計(jì)算得到的冰載荷與柱狀結(jié)構(gòu)的振動(dòng)響應(yīng)與文獻(xiàn)[14]中的計(jì)算結(jié)果吻合較好，驗(yàn)證了計(jì)算方法的合理性。

表2 PD計(jì)算結(jié)果與文獻(xiàn)[14]計(jì)算結(jié)果的對(duì)比Table 2 Comparision of results calculated by Peridynamics and from literature [14]

4 參數(shù)影響分析

4.1 樁徑對(duì)水平冰力與振動(dòng)位移的影響

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

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

表4 不同樁徑作用下結(jié)構(gòu)振動(dòng)位移的計(jì)算結(jié)果Table 4 Structure vibration displacements with various diameters of cylindrical structure

如圖12(a)所示，水平冰力整體隨樁腿直徑增大呈增加趨勢(shì)。考慮當(dāng)樁徑較大時(shí)，兩者間接觸區(qū)域較大，導(dǎo)致了隨機(jī)冰力的增大。如圖12(b)所示，對(duì)比不同樁腿直徑下結(jié)構(gòu)的振動(dòng)響應(yīng)發(fā)現(xiàn)，結(jié)構(gòu)振動(dòng)位移極值和均值均隨樁徑增加而增加?？紤]樁徑增大提升了冰力值總體水平，增大了振動(dòng)系統(tǒng)外力輸入，進(jìn)而導(dǎo)致了振動(dòng)響應(yīng)的增大。

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

此外，如圖13所示，計(jì)算發(fā)現(xiàn)Fmax/Fmean隨徑厚比D/h增大逐漸減小，與Sodhi的實(shí)驗(yàn)結(jié)果一致，符合Kry提出的冰載荷極值與均值隨著結(jié)構(gòu)寬度增加逐漸接近的觀察結(jié)果[26]。

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

4.2 冰速對(duì)水平冰力與振動(dòng)位移的影響

分別取冰速0.2、0.8和1.1 m/s進(jìn)行計(jì)算，樁徑仍取2 m，其他參數(shù)保持不變。水平冰力、結(jié)構(gòu)振動(dòng)位移計(jì)算結(jié)果分別如表5～6所示。

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

表6 不同冰速下結(jié)構(gòu)振動(dòng)位移的計(jì)算結(jié)果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時(shí)水平冰力極值較0.5 m/s時(shí)較高，與Sodhi和Morris的平板冰與剛性立柱結(jié)構(gòu)接觸實(shí)驗(yàn)結(jié)果具有一定的相似性[27]。如圖14(b)所示，振動(dòng)位移均值隨冰速增大而增大，但極值卻隨冰速呈負(fù)增長(zhǎng)趨勢(shì)。通過(guò)對(duì)比各冰速下結(jié)構(gòu)的振動(dòng)位移時(shí)程曲線發(fā)現(xiàn)，0.8和1.1 m/s下結(jié)構(gòu)振動(dòng)位移曲線同0.5 m/s時(shí)類似，均呈隨機(jī)性。但0.2 m/s下結(jié)構(gòu)振動(dòng)位移幅值持續(xù)維持在0.2～0.4 mm，且具有較明顯的周期性特點(diǎn)，如圖15所示。因此，認(rèn)為在0.2 m/s冰速下柱狀結(jié)構(gòu)的振動(dòng)為穩(wěn)態(tài)振動(dòng)。

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

圖15 結(jié)構(gòu)振動(dòng)位移的時(shí)程變化(0.2 m/s)Fig.15 Vibration displacement of cylindrical structure during the simulation (0.2 m/s)

黃焱等[27]的實(shí)驗(yàn)結(jié)果顯示，當(dāng)柱狀結(jié)構(gòu)與海冰相互作用發(fā)生穩(wěn)態(tài)振動(dòng)時(shí)，結(jié)構(gòu)振動(dòng)位移主頻與結(jié)構(gòu)的自振頻率近似。對(duì)0.2 m/s冰速下的柱狀結(jié)構(gòu)振動(dòng)位移時(shí)程數(shù)據(jù)進(jìn)行快速傅里葉變換(FFT)，計(jì)算結(jié)果如圖16所示。由圖可知，結(jié)構(gòu)振動(dòng)位移在頻域上具有6.625 Hz的主頻，與系統(tǒng)自振頻率6.497 Hz相近，驗(yàn)證了結(jié)構(gòu)發(fā)生穩(wěn)態(tài)振動(dòng)。

圖16 結(jié)構(gòu)振動(dòng)位移頻域分析結(jié)果Fig.16 Analysis results of structure vibration displacements in frequency domain

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

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

由圖17可知，F(xiàn)max/Fmean隨著v/h增大呈微幅減小的趨勢(shì)，這與Sodhi的實(shí)驗(yàn)結(jié)果相一致[27]。

5 結(jié)論

1)總體來(lái)說(shuō)，海冰與柱狀結(jié)構(gòu)相互作用引發(fā)的水平冰載荷與樁徑、冰速呈正相關(guān)關(guān)系。結(jié)構(gòu)的振動(dòng)位移隨著樁徑增大而增大，隨著冰速增大而減小。

2)柱狀結(jié)構(gòu)發(fā)生穩(wěn)態(tài)振動(dòng)時(shí)，其結(jié)構(gòu)振動(dòng)位移具有較為明顯的主頻且與系統(tǒng)自振頻率近似。

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