工程材料这一特点，更多的考虑到冰的力学性质，因而对理论分析法而言有了较大改进。
近年来，离散元方法得到广泛的关注，尤其在对冰的研究领域内发挥着重要作用。该方法在细观上能够合理地描述冰材料的结构性质，在宏观上可以准确地模拟出冰与结构作用中的破碎过程[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.