河冰单轴压缩破坏过程细观数值仿真

邓 宇，王 娟，李志军

（1.黄河水利委员会 黄河水利科学研究院，河南 郑州 450003；2.郑州大学 水利与环境学院，河南 郑州 450001；3.大连理工大学 海岸和近海工程国家重点实验室，辽宁 大连 116024）

1 研究背景

河冰内部细观微裂纹的变化是导致其复杂宏观表象的重要因素，对河冰内部细观微裂纹开裂与扩展机理的研究将有助于从本质上解释河冰开河期的宏观力学行为［1］。随着细观力学理论的不断完善和数值模拟方法的成熟，利用数值仿真技术，研究冰的细观结构与其宏观力学特性的关系已成为冰研究的一个热点。相关学者围绕冰的晶粒位错堆积［2］、晶界滑移［3］和晶体弹性各向异性的机制［4］等细观结构理论已开展了大量的研究工作。近年来，关于高原冰［5］、冰川冰［6］、极地冰［7］、海冰［8］、大气冰［9］等有关冰材料细观结构观测研究成果也不断涌现，与上述冰类型相比较，河冰由于自身材料的特殊性、形成及演化过程的多变性及有效数据获取的季节性和时效性，其细观结构、物理力学性能及开裂演变规律的研究仍处于起步阶段［10-11］。相关学者在冰的数值计算方面也取得了长足的进步，王军等［12］对冰塞形成及其演变机理问题、杨开林等［13］对渠道冰问题、季顺迎等［14］对海冰与海洋平台结构物相互作用问题、李志军等［15］对海冰断裂问题、韩雷等［16］对冰板与锥面结构相互作用中的弯曲破坏问题、茅泽育等［17］对冰盖下水流流动的掺混特性问题均开展过数值模拟研究，但基于细观冰力学理论，采用数值模拟方法，从细观尺度出发对河冰结构特性、材料特性及微裂纹的演变规律进行研究，进而构建河冰宏细观多尺度的响应机制，现有研究成果仍有待于进一步丰富。

本文采用数值模拟与试验实测资料分析相结合的方法，以黄河河道内封冻期的天然河冰为研究对象，构建了黄河冰强度的细观计算模型，从细观尺度对河冰的单轴压缩强度进行了分析，对不同初始缺陷含量变化、不同晶粒尺寸变化下河冰压缩破坏过程进行了模拟。该方法不仅能够直观地描述冰晶体破坏的全过程，还可反映出冰晶体的形状、分布形式和界面过渡区等细观组分对河冰宏观力学性能的影响，有效突破河冰宏观力学性能试验因取样、时间、地域、环境等各种客观因素及试验过程中人工操作主观因素的影响，从本质上去认知河冰材料的力学性能特点。

2 河冰的采样观测及物理试验概况

黄河的冰凌问题主要集中在内蒙古河段。2014年黄河封冻期间，在黄河内蒙古典型河段头道拐现场钻取了冰样、并依照垂直冰面方向和平行冰面方向分层切取制作了1 mm厚的冰薄片，将加工好的冰薄片放置于费氏台正交偏光镜下，观测黄河冰的晶体结构特征和冰内气泡的分布和尺寸。现场观测结果表明，沿冰面的垂直方向，黄河冰晶体由颗粒冰层过渡到柱状冰层，再由柱状冰层过渡到颗粒冰层，而后又过渡到柱状冰层（图1）。将黄河河道内取出的部分冰坯加工成7.0 cm×7.0 cm×17.5 cm的标准单轴压缩试样，运回低温实验室进行试验。试验设备采用CSS—44100型电子万能试验机，冰试样在压缩前放入低温箱中恒温12 h以上以保证达到热平衡。试验温度分别为-3℃、-5℃、-7℃、-10℃和-15℃，加载方向分为垂直冰面方向和平行冰面方向。试验得到的数据及具体分析结果详见文献［18］。

图1 垂直冰面方向冰晶粒分布结构

3 河冰细观数值模型

本文研究的河冰晶体结构形式为柱状冰。从细观层面看，河冰的组成包括冰晶粒、泥沙颗粒、气泡、杂质等，本文将河冰看作由冰晶粒和气泡、微裂纹等初始缺陷组成的两相材料，基于ANSYS软件的二次开发平台，编程建立了河冰细观数值计算模型。主要内容包括：河冰细观结构生成、本构关系和破坏准则建立、细观参数确定。

3.1 冰晶粒的生成 基于Voronoi多边形随机性及图形多样性的特征，并结合黄河冰细观结构试验的观测结果，采用Voronoi多边形对河冰晶粒尺寸及分布进行了模拟，生成了随机分布的冰晶粒棱柱体试件（图2（b），其中晶粒颜色随机生成以便区别晶粒轮廓），同黄河冰典型横断面晶体结构相比较（图2（a）），具有较相似的形状。

图2 黄河冰晶粒分布结构

3.2 晶界和初始缺陷的生成 由于晶界厚度非常小，采用无厚度的弹簧单元模拟：以切向弹簧模拟晶粒间的滑移，以法向弹簧模拟晶粒间的挤压（图3（b））。对于多晶冰而言，其内部初始缺陷包括：气泡、杂质、微裂缝等。其中气泡最为常见，影响也相对较大。本文假设河冰内部的初始缺陷分布在冰晶的交界面上，且为随机分布（图3（c）），采用弱化的弹簧单元模拟初始缺陷，弱化的弹簧单元等效作用区域记为初始裂纹面积，则初始缺陷含量为初始裂纹总面积与晶界总面积之比（二维情况下为边长），含初始缺陷的河冰有限元模型如图3（a）所示。

3.3 本构关系及破坏准则 河冰在破坏时表现出较强的脆性［1］，其宏观性能的影响因素主要包括：细观结构和细观组分性能。为了分析细观结构对河冰宏观性能的影响，将细观各相看作是均质的线弹性材料，假设晶粒、晶界均为弹脆性材料，且在同一横截面内各向同性。晶粒的破坏准则采用最大拉应力准则，晶界的破坏考虑为拉断和剪切两种破坏形式，采用法向和切向弹簧模拟晶界界面处开裂和滑移，晶界的力学性能参数，如弹模、强度等，根据晶粒的性能选取。以晶界处弹模和强度以及弹簧作用面积估算法向弹簧的刚度及最大允许拉伸值，切向最大滑移值取为法向拉伸值的10倍，破坏准则如式（1）所示：

式中：FN、FS分别为法向和切向弹簧单元的恢复力（下标N为晶界表面法线方向，S为晶界切线方向）；为弹簧受拉极限承载力；为弹簧剪切极限承载力，取值如式（2）所示：

式中：σt为冰晶界抗拉强度；τb为冰晶界剪切强度；Ae为弹簧等效作用面积。

假设单元破坏后失去承载能力，对其材料属性进行修改。其中，破坏后的晶粒单元弹性模量取为极小值（原始值的10-5～10-4），破坏后的晶界采用弱化的弹簧单元模拟，弹簧单元的相关参数取值与初始缺陷单元相同。

图3 黄河冰有限元模型

3.4 细观参数的确定 河冰细观计算模型中需确定两类参数：细观结构参数和细观材料参数。其中，细观结构参数主要包括晶粒尺寸、晶粒分布情况以及晶界的构成方式、初始缺陷的确定等。细观材料参数主要包括冰晶粒的弹性模量、强度、晶界的力学性能、晶粒各向异性等参数，河冰细观计算过程中主要的参数取值见表1［19］。

3.5 解算方法 河冰细观模型在加载方向和试件底面中心区域水平方向施加约束。单轴压缩强度模拟时，荷载采用位移加载方式逐级施加，每步荷载求解后，根据破坏准则判断模型中单元的破坏情况，有单元破坏后需重新计算当前荷载步。解算的具体过程如下：施加边界条件，根据最大荷载位移和求解步数分级加载；每步荷载施加后进行应力分析，得到整体的位移场和应力场；若有单元破坏，改变单元的材料属性，使其在后面的计算中对刚度矩阵不再贡献，并返回同一荷载步重新计算，直至无新的破坏单元出现；通过约束端面的节点反力计算外荷载大小，除以加载端面积得到构件的名义应力；若名义应力下降为最大荷载的30%以下判断为结构破坏，结束计算；否则施加下一步荷载。

表1 河冰细观计算主要参数取值

4 河冰单轴压缩强度及开裂过程模拟

4.1 计算结果分析 为了验证模型的适用性，进行了黄河冰单轴压缩破坏过程的数值模拟，计算分析了河冰单轴压缩荷载作用下的力学响应。在加载初始阶段，在晶界处初始缺陷位置出现应力集中现象（图4（a）），裂纹沿着晶界向外扩展，并逐渐聚集；随着荷载的进一步增加，裂纹积聚成核，并向晶粒处扩展，河冰晶粒开始起裂（图4（b））；随后在荷载的持续作用下，出现穿晶破坏（图4（c）），宏观裂缝发展（图4（d））至最终破坏（图4（e））。在整个模拟过程中，出现的裂缝多平行于荷载方向，裂缝扩展比较凌乱，多条小裂缝贯穿并导致了最终承载力的丧失，在加载最后阶段出现了大面积的破坏单元。图4（f）为河冰物理试验的最终破坏试样图，试样上部较亮部分为主要破坏区域，在可见光下直视可观察到内部裂缝多平行于荷载方向，试样下部分布有比较散乱的裂缝和微裂纹。河冰单轴压缩物理试验的这种破坏模式与模拟分析观察的情况近似。

河冰单轴压缩强度计算结果如图5所示，图5中列举了河冰模拟计算结果、黄河冰部分物理试验及典型文献试验［20］的河冰单轴压缩强度。通过对比分析可知，河冰的单轴压缩强度模拟结果与黄河冰物理试验结果的误差在5%左右，两者结果相差较小。

图5 数值模拟计算结果与试验结果对比

4.2 单轴压缩强度影响因素分析（1）初始缺陷的影响。初始缺陷的存在是材料破坏的关键因素。河冰初始缺陷的含量在10%以内，但其分布尚未有明确的结论。针对黄河冰，其内部初始缺陷包含气泡、杂质、泥沙颗粒、微裂缝等。通过对黄河冰内气泡含量观测结果可知，黄河封冻期冰内气泡含量集中在1%～3%之间。根据黄河冰的结构特征并结合试验观测结果，本文中数值分析的初始缺陷含量取为2%～8%，假定其随机分布在晶界上，文中主要分析了初始缺陷随机分布的位置、初始缺陷的含量等对河冰破坏过程及强度的影响。图6列举了初始缺陷含量为6%时，4种不同初始缺陷分布的河冰破坏形态，从4幅图中可看出，试件最终破坏均为多条纵向平行裂缝的压碎破坏，初始缺陷的随机位置对破坏模式没有明显影响。

不同初始缺陷位置下河冰单轴压缩强度的模拟值见表2，初始缺陷位置对河冰单轴压缩强度有一定的影响，但是波动范围不大，在可接受的误差范围内。图7描述了不同初始缺陷含量下河冰单轴压缩强度的变化。随着初始缺陷含量的增加，河冰轴压强度呈下降趋势，当初始缺陷含量由2%增加到5%时，河冰轴压强度降低了约50%，当初始缺陷含量继续增大，河冰轴压强度下降幅度渐缓。

图6 初始缺陷对河冰破坏形态的影响

表2 不同初始缺陷位置下河冰单轴压缩强度模拟值

图7 初始缺陷含量对河冰单轴压缩强度的影响

（2）晶粒尺寸的影响。河冰的晶粒尺寸对河冰的力学性能影响很大。通过对黄河冰细观结构的观测，在柱状冰层，黄河冰晶粒尺寸在4～8 mm之间，且随着深度的增加，冰晶粒尺寸逐步增大。为分析晶粒尺寸对河冰宏观性能的影响，假定河冰细观结构中所有晶粒尺寸为均一值，对不同晶粒尺寸的河冰轴压开裂过程进行了模拟。图8分别针对晶粒尺寸d=4 mm、d=6 mm、d=8 mm、d=10 mm的试件典型破坏模式进行了计算分析，由图8可看出，晶粒尺寸对河冰最终破坏形态无明显影响。图9描述了河冰强度与晶粒尺寸之间的变化关系，分析可知，轴压强度与晶粒尺寸d-1/2呈线性相关关系。

图8 不同晶粒尺寸的河冰细观结构及典型破坏

图9 河冰抗压强度模拟值随d-1/2变化的趋势线

5 结论

本文在现有研究成果及河冰物理试验研究的基础上，利用细观冰力学理论，建立了考虑冰晶粒、冰晶界和初始缺陷的河冰细观数值计算模型，模拟了河冰在单轴受压荷载作用下的开裂过程，分析了晶粒尺寸、初始缺陷含量等各相组分对河冰单轴压缩强度的影响，成果可用于河冰材料细观尺度上断裂过程和断裂机理的了解和认知，也为分析宏观冰盖体的断裂过程提供了依据和参考，主要结论如下：

（1）模拟了河冰在单轴压缩荷载作用下的开裂破坏过程，其最终破坏形态与物理试验结果基本吻合，裂缝发展过程与物理试验观察情况较为一致，河冰的单轴压缩强度模拟值与试验结果相差较小。

（2）河冰数值模拟结果表明，随着初始缺陷含量的增加，河冰单轴压缩强度呈下降趋势，在初始缺陷含量为5%时，河冰单轴压缩强度约为初始缺陷含量2%时的一半，初始缺陷含量继续增大后，河冰单轴压缩强度下降幅度渐缓。

（3）河冰数值模拟结果表明，河冰单轴压缩强度与晶粒尺寸的d-1/2呈线性相关关系。

本文结果仅针对河冰准静态单轴压缩强度加载情况进行了数值模拟计算，对于复杂加载条件下的力学行为仍需进一步研究。