多晶冰的二元介质本构模型

Binary Medium Constitutive Model for Polycrystalline Ice

摘要

自然界中冰川广布，在外界环境的扰动下，冰川冰会发生变形，诱导冰川灾害的发生。因此，研究多晶冰的力学特性和本构模型对于理解冰川灾害的发生具有重要的理论价值。本文基于多晶冰试样在不同温度和围压试验条件下的三轴试验数据及连续介质力学等理论框架建立了多晶冰的二元介质本构模型。在所建立的本构模型中，多晶冰试样由胶结元和摩擦元共同构成，其中，胶结元表现为线弹性力学特性，摩擦元表现为弹塑性力学特性。在外荷载作用下，胶结元逐渐破损转变为摩擦元，两者共同承担外部荷载。为了描述胶结元和摩擦元在代表性单元体中的关系，借助均匀化理论并引入破损率和应变集中系数，以描述试样在受荷载过程中内部结构的演化和应变分布的不均匀性。试验结果表明，多晶冰试样的应力‒应变曲线呈应变软化，体积应变曲线随着围压的增加从先体缩后体胀转变为完全的体缩趋势。通过试错法等方法给出在-2 ℃条件下各参数与围压的关系及在-6 ℃部分试验条件下的参数取值，并对比了试验和模拟曲线，验证了本文模型的合理性。对部分结构参数进行敏感性分析，表明了模型的通用性。

Abstract

Objective Glaciers are distributed widely in nature. Under the disturbance of external environmental factors, glacier ice becomes unstable, which triggers the occurrence of glacier disasters. Therefore, studying the mechanical properties and constitutive models for polycrystalline ice has important theoretical value for understanding the triggering mechanism of glacier disasters. Methods For preparing the polycrystalline ice samples, ice particles with a diameter of less than 2 mm were selected and layered in the mold for compaction. The samples were then degassed and saturated using the bottom-in water saturation method. On the mesoscopic level, based on the preparation method, the samples consist of both skeleton ice and pore ice. The confining pressures at -2 and -6 ℃ in the triaxial tests were 0.5, 1.0, 3.0, and 4.0 MPa, respectively. The instrument used for the tests was an MTS triaxial testing machine. Based on the framework of breakage mechanics and homogenization, a Binary Medium Constitutive Model was established to describe the mechanical properties of polycrystalline ice samples. This model adopted a mesoscopic approach to describe the failure process of polycrystalline ice. The samples in the model consist of two parts: bonded elements and frictional elements. At the beginning of the tests, the samples were intact and composed entirely of bonded elements, which exhibited elastic-brittle mechanical behavior. With the increase in external load, the local pore ice in the samples is destroyed, and the skeleton ice near the pore ice slips and crushes, showing the characteristics of loose ice particles. The frictional elements are composed of the broken parts of the samples, exhibiting elastic-plastic mechanical behavior. In the RVE (representative volume element), during loading, the bonded elements gradually fracture and transform into frictional elements, both of which jointly bear external loads. Finally, the polycrystalline ice becomes entirely composed of frictional elements. The breakage ratio and the strain concentration factor were introduced using homogenization theory to describe the evolution of the internal structure and the non-uniform strain distribution under external load in the RVE. The linear elastic constitutive relation in the model was applied to the bonded elements. For the frictional elements, the ideal elastoplastic constitutive relation was adopted, in which the material was linear elastic before failure and ideally plastic at failure. The yield criterion of the frictional elements follows the Mohr-Coulomb yield criterion, assuming a non-associated flow rule. The breakage ratio is expressed as a function of strain, describing the transformation from bonded to frictional elements in the RVE. The breakage ratio function adopted an empirical formula. The strain concentration factor is also expressed as a function of strain, representing the relationship between the strain of the bonded elements and that of the RVE. The parameters of the model include two types: material parameters and structural parameters. The material parameters were obtained from test results and previous research, whereas the structural parameters were derived from macro analysis of test data or trial-and-error methods. Results and Discussions The deviatoric stress-axial strain curves of the polycrystalline ice samples showed a strain-softening trend, which can be divided into the elastic-brittle stage, plastic yield stage, softening stage, and plastic flow stage. During the plastic flow stage, the curves stabilized after reaching the residual strength. Under different test conditions, the samples reached peak stress when the axial strain approached 1%, indicating that polycrystalline ice exhibited brittle material. As the confining pressure increases, the peak stress of polycrystalline ice rises, but when the confining pressure exceeds a specific value, the peak stress decreases. The results confirmed that a threshold exists when the confining pressure is less than 4 MPa. Before and after this threshold, the response of peak stress to confining pressure was opposite. At low confining pressures, the volumetric strain-axial strain curves of polycrystalline ice samples showed slight compression at first, followed by significant dilatancy. As the confining pressure increases, the curves display complete compression under high confining pressures. The tested curves at -2 °C under different confining pressures were compared to the simulated curves, and both exhibited similar variation trends. In addition, the evolution law of parameter values with confining pressure was established. The tested curves at -6 °C under low and high confining pressures were compared to the simulated curves. Conclusions The effect of confining pressure on the samples is manifested in two ways. On one hand, the confining pressure exerts a compacting effect and lateral restraint on the samples, which reduces the pores, strengthens the connection between the ice particles, and constrains the radial displacement of the samples, resulting in an increase in peak stress and volume compression of the polycrystalline ice. On the other hand, the confining pressure produces melting and crushing effects on the ice particles, weakening the bonding between them and causing fracturing and sliding of the particles, which decreases the peak stress of the samples. These two effects coexist and together determine the peak stress and the volumetric strain pattern of the samples. A decrease in temperature reduces the amount of unfrozen water in the samples and intensifies the bonding between ice particles, which increases peak stress and decreases the tendency for volumetric contraction. The rationality of the constitutive model is verified by comparing the curves at different temperatures. In addition, a sensitivity analysis of structural parameters is conducted. The universality of the model is verified by comparing the curves under different parameters.

Graphical abstract

关键词

多晶冰 / 二元介质 / 破损率 / 应变集中系数

Key words

polycrystalline ice / binary medium / breakage ratio / strain concentration factor

本刊网刊

冰川在全球范围内广泛分布于高纬度和高海拔地区。中国作为世界上中、低纬度冰川最发育的国家，约有0.54%的面积被冰川覆盖，主要分布在新疆、甘肃、青海等地区^[1]。随着近年来气候变暖及高海拔地区地壳运动强烈，冰川内部薄弱面在受到扰动后会进一步发育，形成完整横断面，引发冰川滑动、冰崩等冰川灾害^[2‒3]。例如：2016年，西藏地区的阿鲁错湖附近发生冰川崩塌，碎屑行进了7 km，覆盖了沿途8～9 km²的区域；2021年，印度北阿坎德邦喜马拉雅山楠达德维冰川发生了冰崩灾害，摧毁了下游两座水电站，造成80人死亡^[4‒5]。冰川灾害会淹没田地，冲毁建筑，诱发冰崩堵江、堰塞湖等次生灾害，对经济建设和生态环境造成毁灭性打击。针对冰川灾害的研究，一方面，是基于实地监测、遥感技术、数值模拟等方法，评估灾害的冰崩体积和行进距离等^[6]；另一方面，是基于理论对一定环境条件下冰的力学特性进行研究，探究冰川灾害的启动原因。冰川灾害的产生与冰川冰在外界荷载扰动下的失稳密不可分，因此，研究多晶冰的力学特性对理解冰川灾害的产生机理具有重要的理论价值。

目前，研究者对冰试样进行了大量试验，对冰的力学特性和变形特性进行了研究。针对冰的峰值应力和断裂强度的研究适用于破冰船工作等短期工程；针对冰的蠕变行为研究适用于时间、空间维度较大的问题，如极地冰川的退化^[7‒9]。三轴试验是研究材料力学行为的经典方法，有学者通过改变三轴试验中的围压、温度等变量进行了相关研究，提出：多晶冰试样在低围压和高应变率的试验条件下表现出脆性特征^[10]；多晶冰试样的峰值应力在较低围压水平下与围压成正相关，但是在高围压下峰值应力反而降低^[11]；温度对多晶冰试样的强度具有增强作用^[12‒13]；应变率的降低导致冰样从剪切破坏过渡到蠕变破坏^[14‒16]；此外，试样的孔隙率与试样的破坏形式密切相关，随着孔隙率的增加，试样的峰值强度出现削弱^[17]。

本构模型是描述材料力学性质的基础，通过试验数据建立冰的本构模型，可以预测自然环境下冰的破坏行为。针对多晶冰试样的本构模型，Wong等^[18]采用摩尔库伦屈服准则，构建了正常固结破碎冰与时间无关的塑性模型。李洪升^[19]等认为海冰是一种非线性黏性材料，并建立了适用于单轴试验条件的海冰非线性黏弹性本构方程。Gol'dshtein等^[20]采用非相关联的流动法则，给出描述海冰损伤的经验方程。王刚^[21]和Wang^[22]等提出包含最大主剪应力和中间主应力的广义双剪黏弹塑性本构模型，又将Drucker‒Prager屈服准则引入到黏弹塑性本构模型，并对渤海海冰进行了相关模拟。徐洪宇等^[23]对人造多晶冰试样进行了恒定加载速率三轴试验，发现应力‒应变曲线表现为应变软化。胡志强等^[24]利用理想弹塑性材料模拟冰的载荷特性，可较合理地模拟冰受挤压失效的行为特性。Xu等^[25]基于冰与结构之间的相互作用提出非线性黏弹性部分和塑性部分串联的本构模型，并验证了冰与钢板相撞的工况。Yu等^[26]基于冰与结构之间的作用，提出一个包含3部分的现象学模型，研究冰的应力‒应变关系，并考虑了冰破裂长度。

以往研究都将试样整体假定为均匀介质材料，并在热力学、弹塑性力学等知识框架下给出本构关系。冰作为自然界中一种常见的材料，内部结构并非完全胶结在一起，也会存在应力、应变的不均匀分布。二元介质框架通过均匀化理论，将微观角度由破损和未破损两部分组成的状态转换为宏观上的均匀状态，从而构建了完整的二元介质本构模型^[27‒28]。在二元介质模型中，代表性单元体（representative volume element，RVE）由胶结元和摩擦元两部分组成，在受到外荷载作用时，胶结元逐渐转化为摩擦元，宏观力学性质也会发生改变。目前，二元介质模型已经被应用于描述岩石^[29]、融土^[30]、冻土^[31‒33]等材料。本文基于均匀化理论，建立冰的二元介质本构模型，描述多晶冰的力学特性，同时给出确定结构参数的方法，并将试验曲线与模拟曲线进行对比，验证模型合理性。

1 多晶冰试样的力学特性试验

1.1 三轴试样制样方法

三轴试验已被广泛应用于融土、冻土的研究，但传统的抽真空饱和对于冰并不适用，有学者对三轴试验制样方法进行了改进，进而对多晶冰试样进行三轴试验^[21]。试验采用MTS三轴材料试验仪，搭配冷浴循环，压力室温度范围恒定在-30 ℃至室温之间。多晶冰试样尺寸为61.8 mm×125.0 mm，试验材料为人造多晶冰。试验中，选取的温度为（-2 ± 0.2） °C和（-6.0 ± 0.2） ℃，围压分别为0.5、1.0、3.0和4.0 MPa，应变率为1.25 mm/min。

多晶冰试样的制备流程为：

1）将纯净水置于扁平容器中，并在-6 ℃环境下保持24 h，直至完全冻结，将冻结后的单晶冰捣碎成块状。随后，经粉碎机粉碎后的块状冰过2 mm圆孔筛，取粒径小于2 mm的冰晶颗粒进行试验，以保证多晶冰三轴试样的均匀性。

2）将冰颗粒放入模具中并分层压实。压实过程中，将多晶冰样品分为3层，每层冰颗粒的重量接近90 g。压实后，对模具顶部进行光滑处理。最后，将蒸馏水从底部向上缓慢注入模具，得到多晶冰的饱和样品。蒸馏水必须提前放入冰块，确保其温度为0，防止在注入过程中引起模具内冰的融化。同时，为了减少蒸馏水注入的冲刷作用，在模具底部依次铺设滤纸和多孔塑料片。

3）将整套模具置于负温环境下冻结24 h，直至试样完全冻结。冻结完成后，拆模并打磨至试验所需尺寸且表面光滑，用橡胶膜套住试样并将其置于预设为试验温度的恒温箱内备用。

多晶冰试样三轴试验示意图如图1所示。

试验过程保持剪切速率恒定，多晶冰试样在轴向应变较小时就达到峰值应力，因此在试样达到4%轴向应变时停止试验。试验停止后，多晶冰试样出现鼓胀破坏，破坏形态如图2所示。

1.2 三轴试验结果分析

试验过程中应变率保持1.25 mm/min。依据试验数据，不同围压下多晶冰的应力‒应变曲线及体积应变曲线如图3所示。

由图3（a）、（c）可以看出，在-2 ℃和-6 ℃条件下，不同围压下多晶冰试样的应力‒应变关系均呈现应变软化。在初始阶段，应力‒应变曲线呈线性增加的趋势，可以近似为线弹性阶段；随着轴向应变的进一步增大，偏差应力增长趋于缓慢，这是由于试验过程中试样局部发生破坏，试样内部出现冰颗粒的碎裂、滑移、错动，降低了试样整体抵抗变形的能力，进入塑性屈服阶段；试样达到峰值应力后，偏应力迅速下降，表现出明显的软化；随后，随着轴向应变的增加，进入塑性流动阶段，并在达到残余强度后趋于稳定。在不同围压试验条件下，试样均呈现出脆性材料的特点，在轴向应变达到1%左右时就已经达到对应的峰值应力。

由图3（b）、（d）可知：在围压较小的试验条件下，多晶冰试样先发生体缩现象，随后发生较大的体胀，最后体积应变为负值，表现为试样的体胀；随着围压的上升，围压对侧向变形的限制增加，试样的体缩趋势不变，体胀趋势降低，最终体积应变表现为试样的体缩；随着围压的进一步增大，试样的体积应变没有体胀趋势，表现为完全的体缩。对比不同温度下的体积应变曲线可以发现，温度的降低对试样的体缩趋势具有抑制作用。

图4为不同试验条件下多晶冰试样的峰值应力。由图4可知：试验开始时，峰值应力随围压的增加而增大，两者成正相关；当围压增加到一定程度后，试样的峰值应力反而减小，两者成负相关。这是因为当围压较低时，围压对试样内部空隙起到一定的压实作用，强化颗粒间的连接，围压对试样有侧限约束作用，因此随着围压的增加，多晶冰试样的峰值应力同样增加。但在围压较大的试验条件下，冰样会发生压融、压碎，引起试样内部冰颗粒的破裂、滑移，削弱多晶冰颗粒之间的连接，导致多晶冰试样峰值应力下降。基于本文试验结果，可以确定在围压σ₃达到4 MPa之前，总是存在一个围压阈值，在达到该值之前，试样的峰值应力与围压成正相关；当达到或超过这个阈值时，峰值应力与围压成负相关。在本文的后续表述中，定义超过该阈值的围压试验条件为高围压试验条件。

综上，在外荷载作用下，多晶冰试样首先发生局部冰颗粒破坏，出现冰颗粒的压碎、压融，进而引起内部冰颗粒的滑移、错动，并逐渐发展为整个试样的破坏。这种力学行为可以采用二元介质模型进行描述。

2 多晶冰试样的二元介质本构模型

2.1 二元介质本构模型

基于断裂力学，二元介质通过均匀化理论将材料微观的不均匀状态转化为宏观的均匀状态。二元介质模型认为RVE在加载前完全由胶结元构成，宏观上表现为胶结元的力学性质；在受到外荷载作用后，局部受力产生破坏，胶结元向摩擦元转化。在未完全破坏时，RVE微观结构由胶结元和摩擦元共同构成，其微观组分处于动态变化过程；在完全破坏时，RVE在宏观上表现为摩擦元的力学性质。

为描述多晶冰的力学性质，依据试验条件和制样方法，对冰的三轴试样进行分析，依据二元介质模型概念，给出胶结元和摩擦元的定义，如图5所示。在宏观上，多晶冰三轴试样是均匀的、各向同性的材料；在微观上，多晶冰试样由孔隙冰和骨架冰共同构成，骨架冰构成试样的骨架，孔隙冰填充在试样孔隙中，将试样完全胶结成一个整体。在没有承受外荷载时，多晶冰试样完全由胶结元构成，表现出弹性材料的力学特性；随着施加外荷载，试样内局部孔隙冰出现压融、压碎，失去胶结能力，引起范围内骨架冰的破碎和错动。对于出现破坏的局部试样，其表现出松散的冰颗粒力学特性，本文将出现破损的局部整体认为是摩擦元。在加载过程中，孔隙冰逐渐破损，引起骨架冰的错动和嵌挤，使冰颗粒间相互接触，补偿丧失的胶结能力，此时胶结元向摩擦元转化。整个试验过程中，胶结元和摩擦元共同承担外荷载。

在对多晶冰的二元介质本构模型推导的过程中，将未破损材料定义为胶结元（以下标b表示），将破损材料定义为摩擦元（以下标f表示）。

在三轴试验条件下，须遵守：

{\sigma }_{\mathrm{11}}={\sigma }_{\mathrm{1}}，{\sigma }_{\mathrm{22}}={\sigma }_{\mathrm{33}}={\sigma }_{\mathrm{3}}（1）

{\epsilon }_{\mathrm{11}}={\epsilon }_{\mathrm{1}}，{\epsilon }_{\mathrm{22}}={\epsilon }_{\mathrm{33}}={\epsilon }_{\mathrm{3}}（2）

{\epsilon }_{\mathrm{v}}={\epsilon }_{kk}={\epsilon }_{\mathrm{1}}+{\epsilon }_{\mathrm{2}}+{\epsilon }_{\mathrm{3}}={\epsilon }_{\mathrm{1}}+\mathrm{2}{\epsilon }_{\mathrm{3}}（3）

{\epsilon }_{\mathrm{s}}=\sqrt[]{\frac{\mathrm{2}}{\mathrm{3}}\left({\mathit{\epsilon }}_{ij}-\frac{\mathrm{1}}{\mathrm{3}}{\delta }_{ij}{\epsilon }_{kk}\right)\left({\mathit{\epsilon }}_{ij}-\frac{\mathrm{1}}{\mathrm{3}}{\delta }_{ij}{\epsilon }_{kk}\right)}=\frac{\sqrt[]{\mathrm{2}}}{\mathrm{3}}{\left[({\epsilon }_{\mathrm{1}}-{\epsilon }_{\mathrm{2}}{)}^{\mathrm{2}}+({\epsilon }_{\mathrm{2}}-{\epsilon }_{\mathrm{3}}{)}^{\mathrm{2}}+({\epsilon }_{\mathrm{3}}-{\epsilon }_{\mathrm{1}}{)}^{\mathrm{2}}\right]}^{\frac{\mathrm{1}}{\mathrm{2}}}（4）

式（1）～（4）中，{\sigma }_{\mathrm{1}}、{\sigma }_{\mathrm{2}}、{\sigma }_{\mathrm{3}}分别为大、中、小主应力，{\epsilon }_{\mathrm{1}}为轴向应变，{\epsilon }_{\mathrm{2}}和{\epsilon }_{\mathrm{3}}为侧向应变，{\sigma }_{\mathrm{11}}、{\sigma }_{\mathrm{22}}、{\sigma }_{\mathrm{33}}分别为大、中、小主应力方向的应力大小，{\epsilon }_{\mathrm{11}}、{\epsilon }_{\mathrm{22}}、{\epsilon }_{\mathrm{33}}分别为大、中、小主应变方向的应变大小，{\mathit{\epsilon }}_{ij}为应变张量，{\epsilon }_{kk}为对应变张量{\mathit{\epsilon }}_{ij}中i=j的项进行求和，{\delta }_{ij}为克罗内克符号，{\epsilon }_{\mathrm{v}}为体积应变，{\epsilon }_{\mathrm{s}}为广义剪应变，i、j、k的取值范围为1～3。

取代表性单元RVE，其局部应力‒应变关系为：

{\sigma }_{\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}}=D_{\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}}{\epsilon }_{\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}}（5）

式中，{\mathit{\sigma }}_{\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}}、{\mathit{D}}_{\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}}、{\mathit{\epsilon }}_{\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}}分别代表RVE的局部应力张量、局部刚度矩阵和局部应变张量。

依据均匀化理论，平均应力和平均应变可以表示为：

{\left[\begin{array}{c}{\sigma }_{\mathrm{1}}\\ {\sigma }_{\mathrm{3}}\end{array}\right]}_{\mathrm{a}\mathrm{v}\mathrm{e}}=\frac{\mathrm{1}}{V}\int {\left[\begin{array}{c}{\sigma }_{\mathrm{1}}\\ {\sigma }_{\mathrm{3}}\end{array}\right]}_{\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}}\mathrm{d}V（6）

{\left[\begin{array}{c}{\epsilon }_{\mathrm{1}}\\ {\epsilon }_{\mathrm{3}}\end{array}\right]}_{\mathrm{a}\mathrm{v}\mathrm{e}}=\frac{\mathrm{1}}{V}\int {\left[\begin{array}{c}{\epsilon }_{\mathrm{1}}\\ {\epsilon }_{\mathrm{3}}\end{array}\right]}_{\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}}\mathrm{d}V（7）

为了表征荷载作用下试样内部胶结元的破损程度，引入破损率{\lambda }_{\mathrm{v}}对摩擦元在RVE中的体积占比进行描述，其定义为：

{\lambda }_{\mathrm{v}}=\frac{V_{\mathrm{f}}}{V}（8）

式（6）～（8）中，下标\mathrm{a}\mathrm{v}\mathrm{e}代表平均，下标local代表局部，V代表RVE总体积，V_{\mathrm{f}}代表摩擦元所占体积。

依据{\lambda }_{\mathrm{v}}的定义，将胶结元、摩擦元的应力、应变与RVE的应力、应变关系表示为：

\left[\begin{array}{c}{\sigma }_{\mathrm{1}}\\ {\sigma }_{\mathrm{3}}\end{array}\right]=(\mathrm{1}-{\lambda }_{\mathrm{v}}){\left[\begin{array}{c}{\sigma }_{\mathrm{1}}\\ {\sigma }_{\mathrm{3}}\end{array}\right]}_{\mathrm{b}}+{\lambda }_{\mathrm{v}}{\left[\begin{array}{c}{\sigma }_{\mathrm{1}}\\ {\sigma }_{\mathrm{3}}\end{array}\right]}_{\mathrm{f}}（9）

\left[\begin{array}{c}{\epsilon }_{\mathrm{1}}\\ {\epsilon }_{\mathrm{3}}\end{array}\right]=(\mathrm{1}-{\lambda }_{\mathrm{v}}){\left[\begin{array}{c}{\epsilon }_{\mathrm{1}}\\ {\epsilon }_{\mathrm{3}}\end{array}\right]}_{\mathrm{b}}+{\lambda }_{\mathrm{v}}{\left[\begin{array}{c}{\epsilon }_{\mathrm{1}}\\ {\epsilon }_{\mathrm{3}}\end{array}\right]}_{\mathrm{f}}（10）

依据式（9）、（10），可得相应的增量表达式为：

\left[\begin{array}{c}\mathrm{d}{\sigma }_{\mathrm{1}}\\ \mathrm{d}{\sigma }_{\mathrm{3}}\end{array}\right]=(\mathrm{1}-{\lambda }_{\mathrm{v}}^{\mathrm{i}}){\left[\begin{array}{c}\mathrm{d}{\sigma }_{\mathrm{1}}\\ \mathrm{d}{\sigma }_{\mathrm{3}}\end{array}\right]}_{\mathrm{b}}+{\lambda }_{\mathrm{v}}^{\mathrm{i}}{\left[\begin{array}{c}\mathrm{d}{\sigma }_{\mathrm{1}}\\ \mathrm{d}{\sigma }_{\mathrm{3}}\end{array}\right]}_{\mathrm{f}}+\mathrm{d}{\lambda }_{\mathrm{v}}\left({\left[\begin{array}{c}{\sigma }_{\mathrm{1}}\\ {\sigma }_{\mathrm{3}}\end{array}\right]}_{\mathrm{f}}^{\mathrm{i}}-{\left[\begin{array}{c}{\sigma }_{\mathrm{1}}\\ {\sigma }_{\mathrm{3}}\end{array}\right]}_{\mathrm{b}}^{\mathrm{i}}\right)（11）

\left[\begin{array}{c}\mathrm{d}{\epsilon }_{\mathrm{1}}\\ \mathrm{d}{\epsilon }_{\mathrm{3}}\end{array}\right]=(\mathrm{1}-{\lambda }_{\mathrm{v}}^{\mathrm{i}}){\left[\begin{array}{c}\mathrm{d}{\epsilon }_{\mathrm{1}}\\ \mathrm{d}{\epsilon }_{\mathrm{3}}\end{array}\right]}_{\mathrm{b}}+{\lambda }_{\mathrm{v}}^{\mathrm{i}}{\left[\begin{array}{c}\mathrm{d}{\epsilon }_{\mathrm{1}}\\ \mathrm{d}{\epsilon }_{\mathrm{3}}\end{array}\right]}_{\mathrm{f}}+\mathrm{d}{\lambda }_{\mathrm{v}}\left({\left[\begin{array}{c}{\epsilon }_{\mathrm{1}}\\ {\epsilon }_{\mathrm{3}}\end{array}\right]}_{\mathrm{f}}^{\mathrm{i}}-{\left[\begin{array}{c}{\epsilon }_{\mathrm{1}}\\ {\epsilon }_{\mathrm{3}}\end{array}\right]}_{\mathrm{b}}^{\mathrm{i}}\right)（12）

式中，上标i代表初始状态下的参数取值。

定义胶结元和摩擦元模量刚度矩阵分别为{\mathit{D}}_{\mathrm{b}}、{\mathit{D}}_{\mathrm{f}}，则有：

{\left[\begin{array}{c}\mathrm{d}{\sigma }_{\mathrm{1}}\\ \mathrm{d}{\sigma }_{\mathrm{3}}\end{array}\right]}_{\mathrm{b}}=D_{\mathrm{b}}{\left[\begin{array}{c}\mathrm{d}{\epsilon }_{\mathrm{1}}\\ \mathrm{d}{\epsilon }_{\mathrm{3}}\end{array}\right]}_{\mathrm{b}}（13）

{\left[\begin{array}{c}\mathrm{d}{\sigma }_{\mathrm{1}}\\ \mathrm{d}{\sigma }_{\mathrm{3}}\end{array}\right]}_{\mathrm{f}}=D_{\mathrm{f}}{\left[\begin{array}{c}\mathrm{d}{\epsilon }_{\mathrm{1}}\\ \mathrm{d}{\epsilon }_{\mathrm{3}}\end{array}\right]}_{\mathrm{f}}（14）

分析式（12），给出摩擦元应变增量的表达式：

{\left[\begin{array}{c}\mathrm{d}{\epsilon }_{\mathrm{1}}\\ \mathrm{d}{\epsilon }_{\mathrm{3}}\end{array}\right]}_{\mathrm{f}}=\frac{\mathrm{1}}{{\lambda }_{\mathrm{v}}^{\mathrm{i}}}\left\{\left[\begin{array}{c}\mathrm{d}{\epsilon }_{\mathrm{1}}\\ \mathrm{d}{\epsilon }_{\mathrm{3}}\end{array}\right]-(\mathrm{1}-{\lambda }_{\mathrm{v}}^{\mathrm{i}}){\left[\begin{array}{c}\mathrm{d}{\epsilon }_{\mathrm{1}}\\ \mathrm{d}{\epsilon }_{\mathrm{3}}\end{array}\right]}_{\mathrm{b}}-\mathrm{d}{\lambda }_{\mathrm{v}}\left({\left[\begin{array}{c}{\epsilon }_{\mathrm{1}}\\ {\epsilon }_{\mathrm{3}}\end{array}\right]}_{\mathrm{f}}^{\mathrm{i}}-{\left[\begin{array}{c}{\epsilon }_{\mathrm{1}}\\ {\epsilon }_{\mathrm{3}}\end{array}\right]}_{\mathrm{b}}^{\mathrm{i}}\right)\right\}（15）

将式（15）代入式（14）可得摩擦元的应力增量表达式，再将其代入式（11）可得到总应力的表达式为：

\left[\begin{array}{c}\mathrm{d}{\sigma }_{\mathrm{1}}\\ \mathrm{d}{\sigma }_{\mathrm{3}}\end{array}\right]=(\mathrm{1}-{\lambda }_{\mathrm{v}}^{\mathrm{i}})(D_{\mathrm{b}}-D_{\mathrm{f}}){\left[\begin{array}{c}\mathrm{d}{\epsilon }_{\mathrm{1}}\\ \mathrm{d}{\epsilon }_{\mathrm{3}}\end{array}\right]}_{\mathrm{b}}+D_{\mathrm{f}}\left[\begin{array}{c}\mathrm{d}{\epsilon }_{\mathrm{1}}\\ \mathrm{d}{\epsilon }_{\mathrm{3}}\end{array}\right]+D_{\mathrm{f}}\mathrm{d}{\lambda }_{\mathrm{v}}\left({\left[\begin{array}{c}{\epsilon }_{\mathrm{1}}\\ {\epsilon }_{\mathrm{3}}\end{array}\right]}_{\mathrm{f}}^{\mathrm{i}}-{\left[\begin{array}{c}{\epsilon }_{\mathrm{1}}\\ {\epsilon }_{\mathrm{3}}\end{array}\right]}_{\mathrm{b}}^{\mathrm{i}}\right)-\mathrm{d}{\lambda }_{\mathrm{v}}\left({\left[\begin{array}{c}{\sigma }_{\mathrm{1}}\\ {\sigma }_{\mathrm{3}}\end{array}\right]}_{\mathrm{b}}^{\mathrm{i}}-\left[\begin{array}{c}{\sigma }_{\mathrm{1}}\\ {\sigma }_{\mathrm{3}}\end{array}\right]\right)（16）

为了描述胶结元应变与总应变的关系，引入应变集中系数C：

{\left[\begin{array}{c}{\epsilon }_{\mathrm{1}}\\ {\epsilon }_{\mathrm{3}}\end{array}\right]}_{\mathrm{b}}=C\left[\begin{array}{c}{\epsilon }_{\mathrm{1}}\\ {\epsilon }_{\mathrm{3}}\end{array}\right]（17）

根据式（17）给出胶结元应变张量的增量形式：

{\left[\begin{array}{c}\mathrm{d}{\epsilon }_{\mathrm{1}}\\ \mathrm{d}{\epsilon }_{\mathrm{3}}\end{array}\right]}_{\mathrm{b}}=C^{\mathrm{i}}\left[\begin{array}{c}\mathrm{d}{\epsilon }_{\mathrm{1}}\\ \mathrm{d}{\epsilon }_{\mathrm{3}}\end{array}\right]+{\left[\begin{array}{c}{\epsilon }_{\mathrm{1}}\\ {\epsilon }_{\mathrm{3}}\end{array}\right]}^{\mathrm{i}}\mathrm{d}C（18）

将式（17）代入式（10），将C引入摩擦元应变表达式得到：

{\left[\begin{array}{c}{\epsilon }_{\mathrm{1}}\\ {\epsilon }_{\mathrm{3}}\end{array}\right]}_{\mathrm{f}}=A\left[\begin{array}{c}{\epsilon }_{\mathrm{1}}\\ {\epsilon }_{\mathrm{3}}\end{array}\right]（19）

A=\frac{\mathrm{1}}{{\lambda }_{\mathrm{v}}}\left[I-(\mathrm{1}-{\lambda }_{\mathrm{v}})C\right]（20）

式中，\mathit{I}为单位矩阵，A为中间变量。

将式（18）代入式（16）得到RVE的应力增量为：

\left[\begin{array}{c}\mathrm{d}{\sigma }_{\mathrm{1}}\\ \mathrm{d}{\sigma }_{\mathrm{3}}\end{array}\right]=\left[(\mathrm{1}-{\lambda }_{\mathrm{v}}^{\mathrm{i}})(D_{\mathrm{b}}-D_{\mathrm{f}})C^{\mathrm{i}}+D_{\mathrm{f}}\right]\left[\begin{array}{c}\mathrm{d}{\epsilon }_{\mathrm{1}}\\ \mathrm{d}{\epsilon }_{\mathrm{3}}\end{array}\right]+(\mathrm{1}-{\lambda }_{\mathrm{v}}^{\mathrm{i}})·

\begin{array}{c}(D_{\mathrm{b}}-D_{\mathrm{f}}){\left[\begin{array}{c}{\epsilon }_{\mathrm{1}}\\ {\epsilon }_{\mathrm{3}}\end{array}\right]}^{\mathrm{i}}\mathrm{d}C-D_{\mathrm{f}}\mathrm{d}{\lambda }_{\mathrm{v}}\left({\left[\begin{array}{c}{\epsilon }_{\mathrm{1}}\\ {\epsilon }_{\mathrm{3}}\end{array}\right]}_{\mathrm{f}}^{\mathrm{i}}-{\left[\begin{array}{c}{\epsilon }_{\mathrm{1}}\\ {\epsilon }_{\mathrm{3}}\end{array}\right]}_{\mathrm{b}}^{\mathrm{i}}\right)-\\ \mathrm{d}{\lambda }_{\mathrm{v}}\left({\left[\begin{array}{c}{\sigma }_{\mathrm{1}}\\ {\sigma }_{\mathrm{3}}\end{array}\right]}_{\mathrm{b}}^{\mathrm{i}}-{\left[\begin{array}{c}{\sigma }_{\mathrm{1}}\\ {\sigma }_{\mathrm{3}}\end{array}\right]}_{\mathrm{f}}^{\mathrm{i}}\right)\end{array}（21）

化简式（21），以总应变表示胶结元和摩擦元的应变：

\begin{array}{c}\left[\begin{array}{c}\mathrm{d}{\sigma }_{\mathrm{1}}\\ \mathrm{d}{\sigma }_{\mathrm{3}}\end{array}\right]=\left[(\mathrm{1}-{\lambda }_{\mathrm{v}}^{\mathrm{i}})(D_{\mathrm{b}}-D_{\mathrm{f}})C^{\mathrm{i}}+D_{\mathrm{f}}\right]\left[\begin{array}{c}\mathrm{d}{\epsilon }_{\mathrm{1}}\\ \mathrm{d}{\epsilon }_{\mathrm{3}}\end{array}\right]+(\mathrm{1}-{\lambda }_{\mathrm{v}}^{\mathrm{i}})·\\ (D_{\mathrm{b}}-D_{\mathrm{f}}){\left[\begin{array}{c}{\epsilon }_{\mathrm{1}}\\ {\epsilon }_{\mathrm{3}}\end{array}\right]}^{\mathrm{i}}\mathrm{d}C+\left[(D_{\mathrm{f}}^{\mathrm{i}}-D_{\mathrm{f}})A^{\mathrm{i}}+(D_{\mathrm{f}}-D_{\mathrm{b}}^{\mathrm{i}})C^{\mathrm{i}}\right]{\left[\begin{array}{c}{\epsilon }_{\mathrm{1}}\\ {\epsilon }_{\mathrm{3}}\end{array}\right]}^{\mathrm{i}}\mathrm{d}{\lambda }_{\mathrm{v}}\end{array}（22）

2.2 胶结元和摩擦元的本构关系

2.2.1 胶结元

在试验前及加载初期，试样处于未破损状态，整体由胶结元构成，其力学特性表现为线弹性，因此拟用线弹性本构模型来描述胶结元的力学行为。

式（13）中{\mathit{D}}_{\mathrm{b}}的具体形式为：

{\mathit{D}}_{\mathrm{b}}=\frac{E_{\mathrm{b}}(\mathrm{1}-v_{\mathrm{b}})}{(\mathrm{1}+v_{\mathrm{b}})(\mathrm{1}-\mathrm{2}{v}_{\mathrm{b}})}\left[\begin{array}{cccccc}\mathrm{1}& \frac{v_{\mathrm{b}}}{(\mathrm{1}-v_{\mathrm{b}})}& \frac{v_{\mathrm{b}}}{(\mathrm{1}-v_{\mathrm{b}})}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \frac{v_{\mathrm{b}}}{(\mathrm{1}-v_{\mathrm{b}})}& \mathrm{1}& \frac{v_{\mathrm{b}}}{(\mathrm{1}-v_{\mathrm{b}})}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \frac{v_{\mathrm{b}}}{(\mathrm{1}-v_{\mathrm{b}})}& \frac{v_{\mathrm{b}}}{(\mathrm{1}-v_{\mathrm{b}})}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \frac{\mathrm{1}-\mathrm{2}{v}_{\mathrm{b}}}{\mathrm{2}(\mathrm{1}-v_{\mathrm{b}})}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \frac{\mathrm{1}-\mathrm{2}{v}_{\mathrm{b}}}{\mathrm{2}(\mathrm{1}-v_{\mathrm{b}})}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \frac{\mathrm{1}-\mathrm{2}{v}_{\mathrm{b}}}{\mathrm{2}(\mathrm{1}-v_{\mathrm{b}})}\end{array}\right]（23）

式中，v_{\mathrm{b}}、E_{\mathrm{b}}分别为胶结元的泊松比和弹性模量，v_{\mathrm{b}}由前人经验确定，E_{\mathrm{b}}可以通过试验初期的应力‒应变曲线确定。

对式（23）进行三轴试验应力状态下的化简，可以得到新的应力‒应变关系表达式为：

{\left[\begin{array}{c}\mathrm{d}{\sigma }_{\mathrm{1}}\\ \mathrm{d}{\sigma }_{\mathrm{3}}\end{array}\right]}_{\mathrm{b}}=\frac{E_{\mathrm{b}}(\mathrm{1}-v_{\mathrm{b}})}{(\mathrm{1}+v_{\mathrm{b}})(\mathrm{1}-\mathrm{2}{v}_{\mathrm{b}})}\left[\begin{array}{cc}\mathrm{1}-v_{\mathrm{b}}& \mathrm{2}{v}_{\mathrm{b}}\\ v_{\mathrm{b}}& \mathrm{1}\end{array}\right]{\left[\begin{array}{c}\mathrm{d}{\epsilon }_{\mathrm{1}}\\ \mathrm{d}{\epsilon }_{\mathrm{3}}\end{array}\right]}_{\mathrm{b}}（24）

2.2.2 摩擦元

摩擦元采用理想弹塑性本构模型，屈服准则采用摩尔库伦屈服准则。根据理想弹塑性本构模型的定义，当摩擦元受外荷载作用而未达到屈服时，应力状态表现为线弹性，即三轴压缩试验简化后的应力‒应变关系为：

{\left[\begin{array}{c}\mathrm{d}{\sigma }_{\mathrm{1}}\\ \mathrm{d}{\sigma }_{\mathrm{3}}\end{array}\right]}_{\mathrm{f}}=\frac{E_{\mathrm{f}}}{(\mathrm{1}+v_{\mathrm{f}})(\mathrm{1}-\mathrm{2}{v}_{\mathrm{f}})}\left[\begin{array}{cc}\mathrm{1}-v_{\mathrm{f}}& \mathrm{2}{v}_{\mathrm{f}}\\ v_{\mathrm{f}}& \mathrm{1}\end{array}\right]{\left[\begin{array}{c}\mathrm{d}{\epsilon }_{\mathrm{1}}\\ \mathrm{d}{\epsilon }_{\mathrm{3}}\end{array}\right]}_{\mathrm{f}}=D_{\mathrm{f}}^{\mathrm{e}}{\left[\begin{array}{c}\mathrm{d}{\epsilon }_{\mathrm{1}}\\ \mathrm{d}{\epsilon }_{\mathrm{3}}\end{array}\right]}_{\mathrm{f}}（25）

式中：v_{\mathrm{f}}、E_{\mathrm{f}}分别为摩擦元的泊松比和弹性模量，v_{\mathrm{f}}由前人经验确定，E_{\mathrm{f}}可以通过试验曲线确定；{\mathit{D}}_{\mathrm{f}}^{\mathrm{e}}为摩擦元的弹性刚度矩阵，上标\mathrm{e}代表为与弹性相关的参数。

摩尔库伦屈服准则的屈服函数形式为：

f=\frac{\mathrm{1}}{\mathrm{2}}({\sigma }_{\mathrm{1}}-{\sigma }_{\mathrm{3}})-\frac{\mathrm{1}}{\mathrm{2}}({\sigma }_{\mathrm{1}}+{\sigma }_{\mathrm{3}})\mathrm{s}\mathrm{i}\mathrm{n}\phi -c\mathrm{c}\mathrm{o}\mathrm{s}\phi（26）

式中，f为屈服函数，c为摩擦元的黏聚力，\phi为摩擦元的内摩擦角。

采用非相关联的流动法则，给出塑性势函数形式：

g=\frac{\mathrm{1}}{\mathrm{2}}({\sigma }_{\mathrm{1}}-{\sigma }_{\mathrm{3}})-\frac{\mathrm{1}}{\mathrm{2}}({\sigma }_{\mathrm{1}}+{\sigma }_{\mathrm{3}})K\mathrm{s}\mathrm{i}\mathrm{n}\phi -c\mathrm{c}\mathrm{o}\mathrm{s}\phi（27）

式中，g为塑性势函数，K为与围压相关的模型参数。

对于理想弹塑性模型，现有应变关系符合：

{\left[\begin{array}{c}\mathrm{d}{\epsilon }_{\mathrm{1}}\\ \mathrm{d}{\epsilon }_{\mathrm{3}}\end{array}\right]}_{\mathrm{f}}={\left[\begin{array}{c}\mathrm{d}{\epsilon }_{\mathrm{1}}\\ \mathrm{d}{\epsilon }_{\mathrm{3}}\end{array}\right]}_{\mathrm{f}}^{\mathrm{e}}+{\left[\begin{array}{c}\mathrm{d}{\epsilon }_{\mathrm{1}}\\ \mathrm{d}{\epsilon }_{\mathrm{3}}\end{array}\right]}_{\mathrm{f}}^{\mathrm{p}}（28）

式中，上标\mathrm{p}代表与塑性相关的参数。摩擦元的塑性应变张量满足：

{\left[\begin{array}{c}\mathrm{d}{\epsilon }_{\mathrm{1}}\\ \mathrm{d}{\epsilon }_{\mathrm{3}}\end{array}\right]}_{\mathrm{f}}^{\mathrm{p}}=\mathrm{d}\lambda \partial g/{\left[\begin{array}{c}\partial {\sigma }_{\mathrm{1}}\\ \partial {\sigma }_{\mathrm{3}}\end{array}\right]}_{\mathrm{f}}（29）

式中，\mathrm{d}\lambda为非负标量函数。

将式（28）、（29）代入式（25），有

{\left[\begin{array}{c}\mathrm{d}{\sigma }_{\mathrm{1}}\\ \mathrm{d}{\sigma }_{\mathrm{3}}\end{array}\right]}_{\mathrm{f}}=D_{\mathrm{f}}^{\mathrm{e}}\left({\left[\begin{array}{c}\mathrm{d}{\epsilon }_{\mathrm{1}}\\ \mathrm{d}{\epsilon }_{\mathrm{3}}\end{array}\right]}_{\mathrm{f}}^{\mathrm{e}}-\mathrm{d}\lambda \partial g/{\left[\begin{array}{c}\partial {\sigma }_{\mathrm{1}}\\ \partial {\sigma }_{\mathrm{3}}\end{array}\right]}_f\right)（30）

在塑性变形时，应力点停留在屈服面上，基于一致性条件，可以给出\mathrm{d}\lambda的表达式为：

\mathrm{d}\lambda =\frac{\mathrm{1}}{H}\left(\partial f/{\left[\begin{array}{c}\partial {\sigma }_{\mathrm{1}}\\ \partial {\sigma }_{\mathrm{3}}\end{array}\right]}_f\right)D_{\mathrm{f}}^{\mathrm{e}}{\left[\begin{array}{c}\mathrm{d}{\epsilon }_{\mathrm{1}}\\ \mathrm{d}{\epsilon }_{\mathrm{3}}\end{array}\right]}_{\mathrm{f}}（31）

式中，中间变量H的表达式为：

H=\left(\partial f/{\left[\begin{array}{c}\partial {\sigma }_{\mathrm{1}}\\ \partial {\sigma }_{\mathrm{3}}\end{array}\right]}_{\mathrm{f}}\right)D_{\mathrm{f}}^{\mathrm{e}}\left(\partial g/{\left[\begin{array}{c}\partial {\sigma }_{\mathrm{1}}\\ \partial {\sigma }_{\mathrm{3}}\end{array}\right]}_f\right)（32）

将式（31）代入式（30）有

{\left[\begin{array}{c}\mathrm{d}{\sigma }_{\mathrm{1}}\\ \mathrm{d}{\sigma }_{\mathrm{3}}\end{array}\right]}_{\mathrm{f}}=\left[D_{\mathrm{f}}^{\mathrm{e}}-\frac{\mathrm{1}}{H}{H}_{*}{H}_{\#}\right]{\left[\begin{array}{c}\mathrm{d}{\epsilon }_{\mathrm{1}}\\ \mathrm{d}{\epsilon }_{\mathrm{3}}\end{array}\right]}_{\mathrm{f}}（33）

式中，H_{*}、H_{\#}的表达式分别为：

H_{*}=D_{\mathrm{f}}^{\mathrm{e}}\left(\partial g/{\left[\begin{array}{c}\partial {\sigma }_{\mathrm{1}}\\ \partial {\sigma }_{\mathrm{3}}\end{array}\right]}_f\right)（34）

H_{\#}=\left(\partial f/{\left[\begin{array}{c}\partial {\sigma }_{\mathrm{1}}\\ \partial {\sigma }_{\mathrm{3}}\end{array}\right]}_f\right)D_{\mathrm{f}}^{\mathrm{e}}（35）

2.3 破损率与应变集中系数

破损率与材料的应变水平有关，在试验过程中，破损率描述了胶结元向摩擦元转化的过程。本文中，破损率函数形式采用经验公式：

{\lambda }_{\mathrm{v}}=\mathrm{1}-\mathrm{e}\mathrm{x}\mathrm{p}\left[-(a_{\mathrm{v}}{\epsilon }_{\mathrm{v}}{)}^{r_{\mathrm{v}}}-(a_{\mathrm{s}}{\epsilon }_{\mathrm{s}}{)}^{r_{\mathrm{s}}}\right]（36）

式中，a_{\mathrm{v}}、a_{\mathrm{s}}、r_{\mathrm{v}}、r_{\mathrm{s}}为需要确定的破损率参数。

应变集中系数与材料的应变水平有关，它可以表征RVE中的胶结元应变与RVE总应变的关系，其函数形式采用经验公式：

C=\mathrm{e}\mathrm{x}\mathrm{p}\left[-(t_{\mathrm{v}}{\epsilon }_{\mathrm{v}}{)}^{u_{\mathrm{v}}}-(t_{\mathrm{s}}{\epsilon }_{\mathrm{s}}{)}^{u_{\mathrm{s}}}\right]（37）

式中，t_{\mathrm{v}}、t_{\mathrm{s}}、u_{\mathrm{v}}、u_{\mathrm{s}}为需要确定的应变集中系数参数。

3 模型验证与参数敏感性分析

3.1 模型验证

在冰的二元介质本构模型中，胶结元的参数为弹性模量E_{\mathrm{b}}和泊松比v_{\mathrm{b}}；摩擦元的参数为弹性模量E_{\mathrm{f}}、泊松比v_{\mathrm{f}}、黏聚力c和内摩擦角\phi；破损率参数为a_{\mathrm{v}}、a_{\mathrm{s}}、r_{\mathrm{v}}、r_{\mathrm{s}}；应变集中系数的参数为t_{\mathrm{v}}、t_{\mathrm{s}}、u_{\mathrm{v}}、u_{\mathrm{s}}；塑性势函数中引入参数K。其中，胶结元、摩擦元的参数可由试验结果给出，破损率、应变集中系数的参数通过模拟试验数据得到。

在-2 ℃的试验条件下，当应变率为1.25 mm/min时，不同围压下的试验与数值模拟结果对比如图6所示；表1、2分别为本构模型材料参数和结构参数的取值，其中，P_a为标准大气压，为0.101 MPa。由表1、2可知，大部分参数的取值与围压相关。

图6中，试验曲线与模拟曲线具有相似的变化趋势。在初始阶段，应力‒应变曲线近似线弹性关系，该模型可以很好地模拟试样的应变软化特性；对于体积应变曲线，模型可以模拟低围压下体积应变从剪缩到剪胀的现象及高围压下体积应变呈现体缩的现象，这与试验规律一致。试验曲线与模拟曲线的对比表明，本文建立的二元介质本构模型是合理的。

为验证-6 ℃试验条件下模型的可行性，分别选取低围压0.5 MPa和高围压4.0 MPa的两组试验对比试验和数值模拟结果，如图7所示。由图7可知，模拟值与试验值较为接近，多晶冰的二元介质本构模型可以描述-6 ℃试验条件下试样的应力‒应变曲线和体积应变曲线规律，进一步证明了模型的合理性。表3为-6 ℃条件下，两种围压条件下的部分参数取值。

3.2 参数敏感性分析

在围压为0.5 MPa、温度为-2 ℃的试验条件下，对二元介质本构模型中的参数a_{\mathrm{v}}、a_{\mathrm{s}}、r_{\mathrm{s}}、t_{\mathrm{v}}、t_{\mathrm{s}}、u_{\mathrm{v}}、K进行敏感性分析，结果如图8所示。

由图8（a）可知，随着a_{\mathrm{v}}的增加，应力‒应变曲线的峰值应力和残余强度均表现为降低趋势，体积应变曲线虽然均表现为先体缩后体胀的趋势，但是体胀量逐渐降低。由图8（b）可知：随着a_{\mathrm{s}}的增加，应力‒应变曲线的峰值应力降低，但残余强度增加，出现软化程度降低的趋势；对于体积应变，参数变化对前期的体缩趋势影响不大，但体胀趋势出现明显的降低。由图8（c）可知，峰值应力、残余强度和体胀量均随r_{\mathrm{s}}的增加而增大。由图8（d）可知：t_{\mathrm{v}}对应力‒应变曲线影响较小；随着t_{\mathrm{v}}的增加，体积应变曲线的体胀趋势逐渐增大，依旧表现为先体缩后体胀的趋势。由图8（e）可知：t_{\mathrm{s}}的增加引起应力‒应变曲线从应变硬化向应变软化过渡，而峰值应力逐渐降低；对于体积应变曲线则表现出促进体胀产生的作用。由图8（f）可知：u_{\mathrm{v}}的增加同样可以促进应力‒应变曲线从硬化趋势向软化趋势转变，且随着u_{\mathrm{v}}的增加，残余强度出现减小的规律；u_{\mathrm{v}}对体积应变曲线影响较小。由图8（g）可知：K的降低引起应力‒应变曲线峰值应力和残余强度的降低，而且削弱了试样的软化趋势；随着K的减小，体积应变曲线由先体缩后体胀的趋势变为完全体缩的趋势。综上，通过改变参数取值，模拟曲线可以描述应力‒应变曲线的硬化规律或软化规律；而体积应变曲线可以表现出明显的体缩或体胀趋势，表明多晶冰的二元介质本构模型具有通用性。

4 结 论

1）在给定的试验条件下，多晶冰试样的应力‒应变曲线均表现出应变软化的曲线特征，而且在较小的轴向应变范围内就已经达到峰值应力。峰值应力曲线表明，多晶冰试样的峰值应力随着围压的增加呈现出先增加后降低的趋势，随着温度的降低呈现出增加的趋势。对于体积应变曲线，围压的增加会减弱多晶冰试样的体胀趋势：在低围压下，试样表现出先体缩后体胀的趋势；在高围压试验条件下，多晶冰试样表现出完全体缩的规律，这是由于围压对试样存在侧向约束的作用。温度的降低对试样的体缩趋势具有抑制作用。

2）基于二元介质框架建立了多晶冰的本构模型以描述多晶冰的力学行为。在模型中，多晶冰试样由线弹性的胶结元和弹塑性的摩擦元共同构成，在受到荷载作用时，胶结元和摩擦元相互转化，依靠破损率建立联系，两者共同承担外界荷载。二元介质本构模型中涉及到胶结元参数、摩擦元参数及结构参数，文中选取-2 ℃试验环境下的试验数据进行试验曲线与数值模拟曲线的对比，并给出了在-2 ℃条件下各参数的取值依据；选取部分-6 ℃试验环境下的试验数据进行对比，表明了二元介质本构模型的合理性；最后，通过模型参数的敏感性分析，证明了模型的通用性。

参考文献

[1] Liu Shiyin, Yao Xiaojun, Guo Wanqin,et al.The contemporary glaciers in China based on the second Chinese glacier inventory[J].Acta Geographica Sinica,2015,70(1):3‒16.

[2] 刘时银,姚晓军,郭万钦,等.基于第二次冰川编目的中国冰川现状[J].地理学报,2015,70(1):3‒16.

[3] Zhang Yingsong, Liu Shiyin, Shangguan Donghui,et al.Thinning and shrinkage of laohugou No.12 glacier in the western Qilian Mountains,China,from 1957 to 2007[J].Journal of Mountain Science,2012,9(3):343‒350. doi:10.1007/s11629-009-2296-4

[4] Li Kunzhong, Zhang Mingzhe, Xing Aiguo.Numerical runout modeling and dynamic analysis of the ice avalanche-debris flow in Sedongpu Basin along Yarlung Zangbo River in Tibet[J].The Chinese Journal of Geological Hazard and Control,2021,32(1):18‒27.

[5] 李昆仲,张明哲,邢爱国.雅鲁藏布江色东普沟崩滑‒碎屑流过程模拟及运动特征分析[J].中国地质灾害与防治学报,2021,32(1):18‒27.

[6] Kääb A, Leinss S, Gilbert A,et al.Massive collapse of two glaciers in western Tibet in 2016 after surge-like instability[J].Nature Geoscience,2018,11(2):114‒120. doi:10.1038/s41561-017-0039-7

[7] Singh R, Aryan V, Joshi M.Understanding the flash flood event of 7th February 2021 in Rishi Ganga basin,Central Himalaya using remote sensing technique[J].Remote Sensing Applications:Society and Environment,2022,26:100744. doi:10.1016/j.rsase.2022.100744

[8] Yin Yueping, Li Bin, Gao Yang,et al.Geostructures,dynamics and risk mitigation of high-altitude and long-runout rockslides[J].Journal of Rock Mechanics and Geotechnical Engineering,2023,15(1):66‒101. doi:10.1016/j.jrmge.2022.11.001

[9] Sayles F H.Triaxial and creep tests on frozen Ottawa sand[C]//Proceedings of 2nd International Permafrost Conference.Washington D C:National Academy of Sciences,1973:384‒391.

[10] Cole D M.The microstructure of ice and its influence on mechanical properties[J].Engineering Fracture Mechanics,2001,68(17/18):1797‒1822. doi:10.1016/s0013-7944(01)00031-5

[11] Breton D J, Baker I, Cole D M.Microstructural evolution of polycrystalline ice during confined creep testing[J].Cold Regions Science and Technology,2016,127:25‒36. doi:10.1016/j.coldregions.2016.03.009

[12] Durham W B, Heard H C, Kirby S H.Experimental deformation of polycrystalline H₂O ice at high pressure and low temperature:Preliminary results[J].Journal of Geophysical Research:Solid Earth,1983,88(Supp1):B377‒B392. doi:10.1029/jb088is01p0b377

[13] Singh S K, Jordaan I J.Triaxial tests on crushed ice[J].Cold Regions Science and Technology,1996,24(2):153‒165. doi:10.1016/0165-232x(95)00017-6

[14] Iliescu D, Schulson E M.Brittle compressive failure of ice:Monotonic versus cyclic loading[J].Acta Materialia,2002,50(8):2163‒2172. doi:10.1016/s1359-6454(02)00060-5

[15] Farid H, Saeidi A, Farzaneh M.Prediction of failure in atmospheric ice under triaxial compressive stress[J].Cold Regions Science and Technology,2017,138:46‒56. doi:10.1016/j.coldregions.2017.03.005

[16] Deng Kun, Feng Xiaowei, Tan Xiaojun,et al.Experimental research on compressive mechanical properties of ice under low strain rates[J].Materials Today Communications,2020,24:101029. doi:10.1016/j.mtcomm.2020.101029

[17] Peternell M, Wilson C J L, Hammes D M.Strain rate dependence for evolution of steady state grain sizes:Insights from high-strain experiments on ice[J].Earth and Planetary Science Letters,2019,506:168‒174. doi:10.1016/j.epsl.2018.10.037

[18] Deng Kun, Feng Xiaowei, Tan Xiaojun,et al.Experimental research on compressive mechanical properties of ice under low strain rates[J].Materials Today Communications,2020,24:101029. doi:10.1016/j.mtcomm.2020.101029

[19] Qiu Wenliang, Peng Rongxin.Research on the numerical simulation for plastic model of ice as building materials under triaxial compression[J].Construction and Building Materials,2021,268:121183. doi:10.1016/j.conbuildmat.2020.121183

[20] Wong T T, Morgenstern N R, Sego D C.A constitutive model for broken ice[J].Cold Regions Science and Technology,1990,17(3):241‒252. doi:10.1016/s0165-232x(05)80004-7

[21] Li Hongsheng, Yue Qianjin, Shen Wu.Viscous effects and model theory of sea ice[J].Journal of Dalian University of Technology,1993,33(2):145‒149.

[22] 李洪升,岳前进,沈梧.海冰的粘性效应及模型理论[J].大连理工大学学报,1993,33(2):145‒149.

[23] Gol'dshtein R V, Marchenko A V.The choice of constitutive relations for an ice cover[J].Journal of Applied Mathematics and Mechanics,1999,63(1):73‒78. doi:10.1016/s0021-8928(99)00011-8

[24] Wang Gang, Ji Shunying, Lv Hexiang,et al.Generalized twin shear viscoelastic-plastic constitutive model with application to Bohai Sea ice dynamics[J].Chinese Journal of Applied Mechanics,2006,23(3):357‒361.

[25] 王刚,季顺迎,吕和祥,等.广义双剪粘弹塑性本构模型在渤海海冰动力学中的应用[J].应用力学学报,2006,23(3):357‒361.

[26] Wang Gang, Ji Shunying, Lv Hexiang,et al.Drucker‒Prager yield criteria in viscoelastic-plastic constitutive model for the study of sea ice dynamics[J].Journal of Hydrodynamics,2006,18(6):714‒722. doi:10.1016/s1001-6058(07)60011-0

[27] Xu Hongyu, Lai Yuanming, Yu Wenbing,et al.Experimental research on triaxial strength of polycrystalline ice[J].Journal of Glaciology and Geocryology,2011,33(5):1120‒1126.

[28] 徐洪宇,赖远明,喻文兵,等.人造多晶冰三轴压缩强度特性试验研究[J].冰川冻土,2011,33(5):1120‒1126.

[29] Hu Zhiqiang, Gao Yan, Yao Qi.A new constitutive model of ice material for ship-ice interaction based on ideal elasto-plastic property[J].Naval Architecture and Ocean Engineering,2016,32(1):65‒73.

[30] 胡志强,高岩,姚琪.一种理想弹塑性模拟的冰材料本构模型[J].船舶与海洋工程,2016,32(1):65‒73.

[31] Xu Ying, Hu Zhiqiang, Ringsberg J W,et al.Nonlinear viscoelastic-plastic material modelling for the behaviour of ice in ice-structure interactions[J].Ocean Engineering,2019,173:284‒297. doi:10.1016/j.oceaneng.2018.12.050

[32] Yu Bowei, Low Y M.A phenomenological model for simulating ice loads on vertical structures incorporating strain rate-dependent stress‒strain characteristics[J].Applied Ma-thematical Modelling,2022,101:132‒156. doi:10.1016/j.apm.2021.07.021

[33] 沈珠江,陈铁林.岩土破损力学:基本概念、目标和任务[C]//岩石力学新进展与西部开发中的岩土工程问题——中国岩石力学与工程学会第七次学术大会论文集.西安:中国岩石力学与工程学会,2002:22‒25.

[34] Shen Z J.Progress in binary medium modeling of geological materials[M]//Modern Trends in Geomechanics.Berlin:Springer Berlin Heidelberg,2006:77‒99. doi:10.1007/978-3-540-35724-7_5

[35] Liu Enlong, Zhang Jianhai.Binary medium model for rock sample[M]//Constitutive Modeling of Geomaterials.Berlin:Springer Berlin Heidelberg,2013:341‒347. doi:10.1007/978-3-642-32814-5_47

[36] Liu Enlong, Nie Qing, Zhang Jianhai.A new strength criterion for structured soils[J].Journal of Rock Mechanics and Geotechnical Engineering,2013,5(2):156‒161. doi:10.1016/j.jrmge.2013.03.002

[37] Wang Pan, Liu Enlong, Zhang De,et al.An elastoplastic binary medium constitutive model for saturated frozen soils[J].Cold Regions Science and Technology,2020,174:103055. doi:10.1016/j.coldregions.2020.103055

[38] Wang Dan, Liu Enlong, Zhang De,et al.An elasto-plastic constitutive model for frozen soil subjected to cyclic loading[J].Cold Regions Science and Technology,2021,189:103341. doi:10.1016/j.coldregions.2021.103341

[39] Wang Pan, Liu Enlong, Zhi Bin,et al.A rate-dependent constitutive model for saturated frozen soil considering local breakage mechanism[J].Journal of Rock Mechanics and Geotechnical Engineering,2023,15(9):2458‒2474. doi:10.1016/j.jrmge.2022.11.017

基金资助

国家自然科学基金项目(U22A20596)