基于離散元方法的豎直群樁冰荷載分析

李輝輝+張寶森+翟必垚+季順迎

摘要：在高緯度海域或高寒地區的江河中，碎冰對海洋平臺樁腿和橋墩等直立群樁結構產生很大的沖擊作用，并在一定條件下造成不同程度的破壞。依據碎冰在自然條件下的離散分布特性，采用Voronoi切割算法構造隨機分布和具有非規則幾何形態的碎冰初始狀態，在此基礎上采用若干個具有黏結-破碎功能的球形顆粒構造具有一定厚度、尺寸和形狀等物理參數的河冰離散單元模型，將豎直群樁簡化為剛體結構，從而對碎冰的運移及其對群樁的沖擊作用進行數值模擬。在此基礎上可確定碎冰撞擊下的群樁冰荷載，并分析冰塊尺寸、冰速和冰厚對樁結構所受冰荷載的影響。

關鍵詞：冰荷載；離散元方法；豎直群樁

中圖分類號：TU473文獻標志碼：A文章編號：1672-1683（2017）01-0126-06

Abstract：In the seas and rivers in high-latitude regions，the broken ice floes can generate obvious ice load on offshore platform piles and bridge piers，and even cause damage of the vertical piles.Based on the discrete distribution characteristics of ice floes under natural conditions，the Voronoi cutting algorithm was used to construct the initial state of randomly-distributed broken ice floes with irregular geometric shapes.The discrete element method was adopted to simulate the drifting process of ice floes and the interaction between ice floes and vertical piles.The ice floe model with certain thickness，size，shape and other physical parameters was constructed with bonded spheres with consideration to their failure effect.The vertical piles were simplified as rigid structures.A numerical simulation was conducted to determine the ice load on the pile group under the impact of ice floes，and analyze the effects of ice size，ice velocity and ice thickness on the ice load of pile structures.

Key words： ice load；discrete element method；vertical pile group

對冰的認識和研究最初始于北方內河航道交通的需要。黃河以北的地區，包括新疆、內蒙和東北三省，秋末冬初，氣溫逐漸降低，太陽輻射減弱，江河出現水凍結成冰的現象，浮冰隨流水順流而下形成秋季流冰[1]。春季由于氣候變暖、太陽輻射加強，封凍期形成的大面積的冰蓋融化破碎形成尺寸大小不等的冰塊，順流而下，形成春季流冰[2]。流冰在經過水工建筑物時，會造成對水工建筑物有嚴重作用的流冰荷載。

冰荷載是寒區內河水工建筑物的主要荷載之一。水工建筑物因冰凌的靜力動力作用遭到破壞的實例時有發生[3]。20世紀50年代官廳水庫中的橋梁墩臺被冰力剪斷；20世紀60年代渤海灣中海洋平臺的烽火臺被流冰撞倒。國外，20世紀60年代美國阿拉斯加庫克灣兩座采油平臺在劇烈的冰振下倒塌；日本雅內港的聲向崎燈標被冰推倒；20世紀70年代芬蘭波茲尼亞灣的鋼制燈塔在流冰撞擊下倒塌；20世紀90年代美國緬因州的一座鋼筋混凝土大橋被流冰撞毀[4]。隨著社會的發展，流冰引發的災害現象也越來越得到人們的重視[5-7]。

目前，國內外對冰荷載的研究包括試驗和理論分析兩種方法。試驗方法主要是在原型結構上進行現場測量和室內實驗室試驗，如加拿大Beaufort海域Molikpaq平臺[8]、加拿大聯邦大橋橋墩以及Baltic海Kemi-I燈塔上的冰荷載測量[9-10]，國內主要集中在渤海的海洋平臺的測量工作[11-13]。理論分析主要有理論模型分析和數值模型分析兩種方法，理論模型有Matlock首次提出的Matlock模型，Karna引入失效區域對其進行了改進[14]，Withalm等其他學者也對其進行不斷的發展與應用[15-16]。數值模型主要包括有限元方法和離散元方法，其更加關注冰作為復雜的工程材料這一特點，更多的考慮到冰的力學性質，因而對理論分析法而言有了較大改進。

近年來，離散元方法得到廣泛的關注，尤其在對冰的研究領域內發揮著重要作用。該方法在細觀上能夠合理地描述冰材料的結構性質，在宏觀上可以準確地模擬出冰與結構作用中的破碎過程[17]。如Hopkins采用擴展圓盤[18-19]來模擬近似圓形的浮冰，并將其應用于對Alaska北岸浮冰的運動、斷裂情況的分析，Hansen采用二維圓盤離散單元模型對錨系船舶在矩形浮冰中的航行進行了分析[20]，另外還有很多學者利用離散元方法對河冰或者海冰進行了不斷深入的研究[21-22]。

為此，本文采用離散元方法，將冰離散為若干個具有黏結-破碎功能的球體顆粒單元，將豎直群樁簡化為剛體結構，對冰運移以及冰對群樁的作用過程進行數值模擬。分析不同工況（冰厚、尺寸、冰速）對樁結構所受冰荷載的影響。

1 冰的離散單元模型

在離散元模擬中，顆粒單元之間視為在法向上由彈簧與阻尼器連接，在切向上由彈簧、阻尼器和滑動摩擦器連接，見圖1（a）。法向彈性力為Fn=Knxn-Cnn，切向彈性力為Fs=min（Ksxs-Csx，μFn）。式中Kn和Ks分別是法向和切向剛度系數，Cn和Cs是法向和切向阻尼系數，x和分別為顆粒的變形和變形率。法向阻尼Cn=ζ2MKn，其中ζn為無量綱黏滯系數，ζn=-lneπn+ln2e，e為回彈系數。切向與法向的剛度、阻尼系數有如下關系：Ks=αKn，Cs=βCn。這里取α=0.5，β=0.5。

考慮單元間的凍結作用，建立了顆粒單元的黏結-破碎功能。兩個黏結單元之間采用平行黏結模型，如圖1（b）所示，即在兩個顆粒單元之間設定彈性黏結圓盤，傳遞軸向力、剪力、彎矩和扭矩，力和力矩都用法向和切向分量來表示，i=ni+si，i=ni+si。式中，ni、ni、si、si分別為力與力矩的法向分量和切向分量。

在顆粒相互作用時，相對位置和相對轉角的變化引起黏彈性力的變化，并通過判斷最大軸向力和剪切力分別與給定的法向黏結強度和切向黏結強度做比較，判斷顆粒間的黏結是否發生破壞，黏結作用是否消失，宏觀的表現是冰塊是否出現裂紋，是否發生破碎。

2 豎直樁冰荷載的離散元分析

本文計算模擬的矩形水道見圖2，長500 m寬100 m，左側冰區長300 m寬100 m。群樁位于冰區前方100 m處，共有5個直立樁，樁的直徑為2 m，樁間距20 m。浮冰是采用Voronoi切割算法構造的，隨機分布且形狀非規則，由單層球體單元黏結，單元直徑為0.3 m，豎直樁采用圓柱體單元。整個冰區的浮冰密集度為80%，厚度為0.3 m，平均尺寸是25 m2，水面高度是20 m，水流速度為1 m/s。整個冰區有1 600個浮冰塊，共由172 695個球體顆粒組成。具體離散元參數見表1。

碎冰在水流速度的拖曳力作用下基本保持勻速運動，到達群樁處時撞擊樁腿，有冰堆積，部分冰塊出現破碎，在直立樁前堆積的冰塊能迅速繞過樁腿繼續向前運動，整個過程如圖3所示。初始時刻碎冰在直水道的分布情況見圖3（a），不同顏色代表不同冰塊。浮冰會發生三個相互垂直方向上的平動和轉動，其中規定浮冰前進的方向為x方向，垂直于水面向上為z方向，x、y、z三個方向符合右手定則。經過計算可以提取得到5個豎直樁的三個方向的冰荷載加以研究。

統計剩下4個豎直樁的x方向的冰力時程曲線如圖5所示，5個樁腿冰力峰值分別是860.77 kN、919.98 kN、1 208.59 kN、1 201.78 kN、1 203.47 kN，相應的均值為85.40 kN、69.22 kN、91.18 kN、94.94 kN、86.14 kN。因為水道邊界形狀規則，碎冰尺寸相近，碎冰運移時冰速保持恒定，所以5個豎直樁所受冰力具有很高的相似性，冰力峰值和均值也都很相近。因此本文在分析冰荷載時，只取水道中間的1號樁腿所受的冰荷載作為研究對象，這會和研究其他樁腿的冰力得到相同的結論。

3 豎直樁上冰荷載的影響因素分析

3.1 冰速對冰荷載的影響

冰速對冰荷載的影響顯而易見，不論是動冰能量理論和動冰強度極限理論，都將冰速作為影響冰力的關鍵因素加以考慮。如工程中流動冰塊對豎直樁沖擊荷裁P=VhARa，V為冰塊流速，h為冰厚，A為冰蓋面積，Ra為壓縮強度，冰速V作為單獨的一項加以研究[4]。而且冰的強度值與加載速率緊密相關，高應變率下表現為脆性材料，低應變率下表現為彈性材料，因而冰速又間接影響冰的壓縮強度與彎曲強度。為研究冰速與冰荷載的關系，本文將保持冰厚為0.4 m，尺寸為25 m2，冰速分別設置成1 m/s、2 m/s、3 m/s、4 m/s、5 m/s、6 m/s對比，得到冰荷載時程。

計算提取不同速度下冰荷載最大值和相應的均值，由此做出曲線見圖6，對比可以看出冰速對于冰荷載的影響很大，冰力隨著冰速增加大體呈線性增長的趨勢。根據動冰能量理論推導的流冰對河道上的橋墩作用冰載荷計算公式：P=0.4VhmsΩR，式中，V是流冰速度；h是冰厚，ms是冰塊質量；Ω冰塊面積；R是局部擠壓強度。在工程應用中，將冰速與冰荷載大小視為成線性相關，數值模擬結果與工程規范很好的對應。這主要是因為冰速的增加會導致浮冰與樁腿碰撞時的動量增大，產生的冰力也隨之增大。模擬結果很好的驗證了根據動冰能量理論所建立規范的合理性。

3.2 冰塊尺寸對冰荷載的影響

浮冰尺寸作為流冰動能的重要影響要素，在動冰與樁撞擊過程中也會影響冰荷載的大小。本文為研究尺寸與冰荷載的關系，將冰速設置為1 m/s，冰厚為0.4 m，生成冰塊平均面積分別為15 m2、20 m2、25 m2、30 m2、35 m2、40 m2，尺寸大小見圖7，其他計算參數同前。

計算模擬了浮冰250 s運移過程，記錄了其x方向的冰力，計算提取得到不同尺寸下冰荷載最大值和相應的均值，并由此做出曲線見圖8。結果表明，冰塊尺寸對冰荷載有一定影響，隨著冰塊尺寸的增加，樁所受冰荷載隨之增加。很多工程應用中都只將冰與建筑物接觸的水平投影寬度作為影響因素考慮，冰塊尺寸間接影響碎冰碰撞樁腿的接觸面積，從而影響冰力大小；另一方面浮冰尺寸的增加會導致整塊浮冰質量增加，動量也隨之增大，產生的冰荷載也隨之增大。由冰荷載均值可以看出，在一定尺寸范圍內，冰荷載遵循動冰能量理論，即冰荷載受碎冰動能決定，隨著尺寸的增大，冰塊的動能增大，冰荷載均值也隨之增大；尺寸超過某個臨界值之后，冰荷載遵循動冰強度極限理論，即統計的冰荷載是在碎冰碰撞樁腿發生破壞的時刻產生，由碎冰的破壞強度決定，尺寸的變化并不影響其破壞強度，因此在圖8（b）后半段的冰荷載均值受尺寸的影響很小。

3.3 冰厚對冰荷載的影響

冰厚與冰速一樣，是影響冰力的重要因素。為了研究冰厚與冰荷載的關系，將保持冰塊尺寸為25 m2，冰速為4 m/s，冰厚分別設置成0.3 m、0.4 m、0.5 m、0.6 m做了四個算例對比，得到的冰荷載時程，計算提取得到每個尺寸下冰荷載最大峰值和相應的均值，由此做出曲線見圖9。可以看出，冰厚對冰荷載影響較大，隨著冰塊尺寸的增大，樁腿所受垂向冰荷載也隨之增大，之間近似成線性關系。主要原因是冰厚增加，冰塊質量增大，從而引起碎冰對樁腿的沖擊動量的增加；另一方面，隨著冰厚增加，模擬冰塊的顆粒粒徑相應增加，因而顆粒之間的黏結力也增大，宏觀上即表現為冰塊破壞強度的增大，冰塊破壞前給樁腿造成的冰荷載的增大。

4 結論

本文采用了離散元方法對冰的運移以及冰對群樁的作用過程進行數值模擬，討論了不同影響因素下，冰荷載的變化規律。以上研究表明，本文建立的離散元法能夠合理計算流凌在水域中的運移過程，可對冰與樁的動力作用過程進行模擬，能得到很好的冰力時程曲線。經過分析可以發現冰厚、冰塊尺寸和冰速都是影響冰力大小的重要因素。

由于影響冰荷載的因素諸多，加之工程規范中的計算公式各系數項的意義及取值原則，在實際工程中不可照搬，應該根據具體情況設置。用離散元法對冰的模擬更加接近實際情況，根據實際工況設置參數進行數值模擬的計算結果可以為規范計算提供很好的參考，也可為實際工程問題提供很好的參照。

參考文獻（References）：

[1] 陸欽年，段忠東，歐進萍，等.黑龍江省的凌汛災害及其區域劃分[J].自然災害學報，2000（3）：54-58.（LU Qin-nian，DUAN Zhong-dong，OU Jin-ping，et al.Ice flood and its zonation in Heilongjiang Province [J].Journal of Natural Disasters，2000（3）：54-58.（in Chinese））DOI：10.3969/j.issn.1004-4574.2000.03.008

[2] 于天來，張礫元，雷俊卿，等.春季流冰對橋墩產生動冰壓力時冰抗壓強度取值探討[J].中外公路，2010，30（3）：168-172.（YU Tian-lai，ZHANG Li-yuan，LEI Jun-qing，et al.Value of ice compressive strength when spring ice produces dynamic ice pressure on the pier [J].Journal of China & Foreign Highway，2010，30（3）：168-172.（in Chinese）） DOI：10.14048/j.issn.1671-2579.2010.03.028

[3] 蘇盛奎.冰荷載對水工建筑物的作用[J].水運工程，1980（5）.（SU Sheng-kui.The effect of ice load on hydraulic structures [J].Port & Waterway Engineering，1980（5）.（in Chinese）） DOI：10.16233/j.cnki.issn1002-4972.1980.05.008

[4] 蔡之瑞，孫柏濤，郭世榮，等.冰荷載的實驗研究與計算方法[J].地震工程與工程振動，1997（4）：49-56.（CAI Zhi-rui，SUN Bai-tao，GUO Shi-rong，et al.Experimental research on ice load and its calculation method [J].Earthquake Engineering and Engineering Dynamics，1997（4）：49-56.（in Chinese）） DOI：10.13197/j.eeev.1997.04.007

[5] 宋波，牛立超，黃帥，等.強震作用下冰水域橋墩動力響應研究分析[J].巖土力學，2014（S1）：285-291.（SONG Bo，NIU Li-chao，HUANG Shuai，et al.Research on dynamic response of pier column structures in ice-plagued area under highly seismic action [J].Rock and Soil Mechanics，2014（S1）：285-291.（in Chinese）） DOI：10.16285/j.rsm.2014.s1.041

[6] 賈玲玲.波浪力對深水橋墩結構在地震和冰荷載作用下的響應影響分析[J].震災防御技術，2010，05（2）：263-269.（JIA Ling-ling.Seismic and ice response analysis of bridge pier in deep water with the wave effect [J].Technology for Earthquake Disaster Prevention，2010，05（2）：263-269.（in Chinese）） DOI：10.3969/j.issn.1673-5722.2010.02.015

[7] 吳云鳳.黃河海勃灣水利樞紐施工期防凌設計[J].人民黃河，2013（11）：90-92.（WU Yun-feng.Ice prevention design during construction of Haibowan water control project on the Yellow River [J].Yellow River，2013（11）：90-92.（in Chinese））DOI：10.3969/j.issn.1000-1379.2013.11.032

[8] Timco G W，Johnston M.Ice Loads on the Molikpaq in the Canadian Beaufort Sea [J].Cold Regions Science and Technology，2003，37：51-68.DOI：10.1016/S0165-232X（03）00035-1

[9] Brown T G，Mttnen M.Comparison of Kemi-I and Confederation Bridge Cone Ice Load Measurement Results [J].Cold Regions Science and Technology，2009，55：3-13.DOI：10.1016/j.coldregions.2008.04.005

[10] Brown T G，Tibbo J S，Tripathi D，et al.Extreme ice load events on the Confederation Bridge [J].Cold Regions Science and Technology，2010，60：1-14.DOI：10.1016/j.coldregions.2009.08.004

[11] Qu Y，Yue Q，Bi X，et al.A random ice force model for narrow conical structures[J].Cold Regions Science and Technology，2006，45：148-157.DOI：10.1016/j.coldregions.2006.05.008

[12] Yue Q，Bi X.Ice-induced jacket structure vibrations in Bohai Sea [J].Journal of Cold Regions Engineering，2000，14：81-92.DOI：10.1061/（ASCE）0887-381X（2000）14：2（81）

[13] 岳前進，畢祥軍，于曉，等.錐體結構的冰激振動與冰力函數[J].土木工程學報，2003，36：16-19.（YUE Qian-jin，BI Xiang-jun，YU Xiao，et al.Ice-induced vibration and ice force function of conical structure [J].China Civil Engineering Journal，2003，36：16-19.（in Chinese）） DOI：10.3321/j.issn：1000-131X.2003.02.004

[14] Krn T，Kamesaki K，Tsukuda H.A numerical model for dynamic ice-structure interaction [J].Computers & Structures，1999，72：645-658.DOI：10.1016/S0045-7949（98）00337-X

[15] Withalm M，Hoffmann N P.Simulation of full-scale ice-structure-interaction by an extended matlock-model [J].Cold Regions Science and Technology，2010，60：130-136.DOI：10.1016/j.coldregions.2009.09.006

[16] Yue Q，Guo F，Krn T.Dynamic ice forces of slender vertical structures due to ice crushing [J].Cold Regions Science & Technology，2009，56（2-3）：77-83.DOI：10.1016/j.coldregions.2008.11.008

[17] 季順迎，李紫麟，李春花，等.碎冰區海冰與船舶結構相互作用的離散元分析[J].應用力學學報，2013（4）：520-526.（JI Shun-ying，LI Zi-lin，LI Chun-hua，et al.Discrete element modeling of ice loads on ship hulls in broken ice fields [J].Chinese Journal of Applied Mechanics，2013（4）：520-526.（in Chinese）） DOI：10.11776/cjam.30.04.D032

[18] Hopkins M A.Discrete element modeling with dilated particles [J].Engineering Computations，2004，21：422-430.DOI：10.1108/02644400410519866

[19] Hopkins M A，Shen H H.Simulation of pancake-ice dynamics in a wave field [J].Annals of glaciology，2001，33：355-360.DOI：10.1016/j.coldregions.2012.02.003

[20] Hansen E H，Lset S.Modelling floating offshore units moored in broken ice：model description [J].Cold Regions Science and Technology，1999，29：97-106.DOI：10.1016/S0165-232X（99）00023，3

[21] Lau M，Lawrence KP，Rothenburg L.Discrete element analysis of ice loads on ships and structures [J].Ships and Offshore Structures，2011，6（3）：211-221.DOI：10.1080/17445302.2010.544086

[22] Lu W，Lubbad R，Lset S.Simulating ice-sloping structure interactions with the cohesive element method [J].Journal of Offshore Mechanics and Arctic Engineering，2014，136：031501.