颗粒材料在宏观上具有非常复杂的力学行为，如非线性、应变软化等.为描述颗粒材料的力学行为，基于连续介质模型已发展了大量弹塑性模型.这些模型需至少包含屈服面、流动法则和硬化(软化)法则3个基本要素.然而上述3个基本要素的确定通常具有一定难度.除弹塑性模型外，近年来，超弹性模型亦被用来描述岩土颗粒材料的力学行为^{[1, 2, 3]}.如Jiang和Liu基于热力学失稳和弹性模量的压力相关性推导了一个能够计算屈服和剪胀性的GE模型^{[2, 4, 5]}.Einav和Puzrin则将压力相关性的超弹性模型和相应塑性模型联合模拟土的力学行为^[6].Houlsby等基于弹性模量为平均应力的幂函数推导了满足热力学要求的超弹性模型^[3].考虑到颗粒材料力学特性，超弹性模型仅适用于应变非常小的颗粒材料体系，当需要描述颗粒材料的破坏行为时，通常将超弹性模型与塑性模型或损伤模型相结合.Volokh^[7]提出的基于超弹性软化模型的概念则给出了另外一条可能的途径.超弹性软化模型的主要思想为^[7]:①真实材料的应变能存在一个极限值；②真实材料的应变能为理想材料应变能的函数.随着应变增大(假定趋于无穷大)，真实材料的应变能将趋于极限值，其应力-应变关系表现出软化现象.必须指出Volokh超弹性模型是针对脆性材料或生物体材料提出的，且其模型假定真实材料和理想材料的应变相同，但应力不同.因此在超弹性软化模型中通常以应变描述应变能.

超弹性软化模型概念清晰，且形式简单，易于数值实现.本文试图将针对脆性材料和生物体材料的Volokh超弹性模型扩展至颗粒材料.考虑颗粒材料主要以剪切变形和破坏为主，文献[8]发展一个颗粒材料的超弹性软化模型，本文在此基础上，进一步建议体积应变能和剪切应变能区别处理，得到了两种修正的超弹性软化模型.并基于ABAQUS平台开发了相应的程序代码.数值算例考察了修正模型模拟应变局部化及应变软化现象的能力，验证了所建立模型的可行性与正确性.

为保持论文完整性，简单介绍Volokh超弹性软化模型的概念^[7].以ψ(C)表示真实材料的应变能，其应力可表示为

假定存在一个与真实材料具有相同应变的理想材料，该理想材料的应变能以W表示，若真实材料性质与理想材料相同，则

此时，真实材料的应变能具有下列特征：

Volokh给出的一个表达式为

将式(7)代入式(1)，相应地有

对比式(1)可以看到，式(9)存在指数因子exp(-W/Φ)，随应变的增加，理想材料的应变能W不断增加，指数因子exp(-W/Φ)从1向0趋减，应力增加的趋势不断减小，并且应力在达到某一极限值后开始减小.这使得式(9)能够反应材料的应变软化现象.基于式(9)的模型称为超弹性软化模型.

对于线性各向同性超弹性材料，在小变形的情况下其应变能函数可表示为

考虑到岩土颗粒材料在破坏过程中体变部分与形变部分所起作用不同，因此在软化模型中采用不同的处理方式.

考虑颗粒材料的变形和破坏特性，通常认为颗粒材料是以剪切变形和破坏为主，为此仅在应变能的剪切部分引入了软化函数，则有

进一步可认为颗粒材料的应变能由部分修正后的应变能(表征软化部分的材料)和部分理想材料的应变能组成，即

可见因子α在一定程度上可视为扰动状态因子^{[9, 10]}.式中α₀为初始软化参数，β为模型参数，表征了软化对等效应变的依赖程度，ε为等效应变，采用计算.

对应变能函数ψ^M1求导可得到相应的应力，进一步对应变求导可得到相应的模量张量.

对于颗粒材料而言，考虑虽然其更容易发生剪切破坏，但也可发生体变部分的破坏.为此本文假定极限体积应变能大于极限剪切应变能，分别考虑其软化，将第二种软化函数表示为

算例1    考虑图 1所示平面应变条件下的颗粒材料集合体，尺寸为0.6 m×0.8 m，上边界单元节点在水平方向被约束，下边界x,y 方向的自由度被约束，左右边界自由，颗粒集合体夹在两刚性板之间，承受由位移控制竖向加载.数值模拟所采用的参数为：体积模量K=60 MPa,剪切模量G=30 MPa，参数Φ=7.5 kJ/m³，参数Φ₁=15.0 kJ/m³，参数Φ₂=7.5 kJ/m³，α₀=0.5，β=10.为考察式(12)与式(15)所建立模型对颗粒材料力学行为的模拟能力及参数α₀和β的影响，图 2给出了Volokh软化模型与第一种和第二种修正的软化模型所预测的位移承载曲线.由图 2可以看出，随着位移加载，Volokh软化方式与两种软化方式所预测的颗粒集合体的承载力均先逐渐升高，当加载位移达到9 mm时，Volokh软化方式的承载力开始降低，材料进入软化阶段，当加载位移达到11 mm时，两种软化方式的承载力开始降低，材料进入软化阶段.图 3给出了加载位移为15 mm时两种软化方式所预测的等效应变分布图.可以看出，两种软化方式均预测出了明显的剪切带，第一种软化方式的高等效应变区域较第二种软化方式大.考察α₀的变化对于材料软化的影响，图 4所示为第一种软化方式在加载位移为15 mm，α₀分别取0.25，0.50和0.75且β取10时的位移承载曲线图.可以看到，α₀在取0.50时软化明显.由于结果的相似性，未讨论第二种软化方式在不同α₀取值时的位移承载曲线图.考察β的变化对于材料软化的影响，图 5给出第一种软化方式在加载位移为15 mm，β分别取5，10和15且α₀取0.50时的位移承载曲线，当β取10时，材料软化明显.第二种软化方式结果类似，故未讨论第二种软化方式在不同β取值时的位移承载曲线图.从图 4可以看出α₀取0.50时等效应变区域最大，剪切带最为明显.从图 5可以看到β取10时等效应变区域最大，且剪切带最为明显.图 6所示为不同加载时刻时第一种软化方式的等效应变分布图，可以看到，随着加载位移的增大，应变局部化逐渐增加，剪切带逐渐产生.由于第二种软化方式在不同α₀与β取值时的等效应变分布与第一种软化方式类似，故未给出第二种软化方式的等效应变分布图.

图1(Figure 1)

图1 平板压缩问题网格图 Fig. 1 Diagram of mesh for the rectangle plane subjectedto a prescribed vertical displacement

图2(Figure 2)

图2 Volokh软化方式与两种软化方式所预测的位移承载曲线图(α0=0.5，β=10) Fig. 2 Load-displacement curves predicted by the Volokh softening model and the two softeningmodels (α0=0.5,β=10)

图3(Figure 3)

(a)—第一种软化方式； (b)—第二种软化方式. 图3 等效应变分布图(α0=0.5，β=10，加载位移=15 mm) Fig. 3 Equivalent strain distribution (α0=0.5,β=10load displacement=15 mm)

图4(Figure 4)

(a)—位移-承载； (b)—α0=0.25； (c)—α0=0.50; (d)—α0=0.75. 图4 第一种软化方式所预测的等效应变分布图 Fig. 4 Equivalent strain distribution predicted by the first softening model

图5(Figure 5)

(a)—位移-承载; (b)—β=5； (c)—β=10; (d)—β=15. 图5 第一种软化方式所预测的等效应变分布图 Fig. 5 Equivalent strain distribution predicted by the first softening model with β=5, 10, 15, α0=0.50under the load displacement of 15 mm

图6(Figure 6)

(a)—加载位移=1.5 mm； (b)—加载位移=6.0 mm; (c)—加载位移=10.5 mm； (d)—加载位移=14.5 mm. 图6 第一种软化方式所预测的不同时刻的等效应变分布图(α0=0.5，β=10) Fig. 6 Equivalent strain distribution predicted by the first softening model at different step time (α0=0.5,β=10)

算例2    考虑图 7所示地基变形，尺寸为34 m×17 m.忽略自重，地基通过位于其上的刚性基础而受到荷载作用，地基和基础之间的接触面假定为理想粘接.所施荷载通过位于基础中心的有限元网格节点垂直作用于地基.地基的底部边界固定，左右两边界水平方向固定.模拟所用参数为体积模量K=60 MPa,剪切模量G=30 MPa，参数Φ=4 kJ/m³，参数Φ₁=10 kJ/m³，参数Φ₂=4 kJ/m³，α₀=0.5，β=10.

图7(Figure 7)

图7 地基平面网格图 Fig. 7 Diagram of mesh for the soil foundation

图 8为加载位移为30 cm时理想模型与两种软化模型所预测的等效应变分布图.从图中看到，第一种软化方式的等效应变局部化区域最大，其应变局部化最明显.理想模型的局部化区域最小，应变局部化不明显.

图8(Figure 8)

(a)—第一种软化方式； (b)—第二种软化方式; (c)—理想模型. 图8 等效应变分布图(α0=0.5，β=10，加载位移=30 cm) Fig. 8 Equivalent strain distribution (α0=0.5,β=10, load displacement=30 cm)

1) 本文将Volokh超弹性软化模型推广到脆性颗粒材料，所建议的模型具有模拟颗粒材料应变软化及应变局部化的能力.

2) 表征软化部分比例的控制因子α本质上类似于扰动状态概念中的扰动因子，其演化应伴随着材料从理想状态向软化状态调整.软化因子α采用随等效应变线性变化的演化方程，能够较好地模拟颗粒材料应变软化及应变局部化现象.

3) 第二种软化方式实质上包含了第一种软化方式.当第二种软化方式的参数Φ₁取到无穷大时，意味着未对体变部分进行软化，这时第二种软化方式退化为第一种软化方式.

4) 相对于颗粒材料的其他本构的模型，超弹性软化模型主要以真实材料的应变能存在极限为基础，概念清晰，表达式简单，易于数值实现，希望能给颗粒材料本构行为的研究工作者一些借鉴.