非共轴特性对岩土材料应变局部化行为影响的数值研究

常江芳，王伟，牛庆合，闻磊

(1.石家庄铁道大学力学系，石家庄 050043；2.石家庄铁道大学土木工程学院，石家庄 050043)

岩土体的应变局部化现象通常表现为宏观剪切带，是导致边坡、堤坝、基坑和挡土墙等岩土工程结构渐进破坏的根本原因[1-2]。岩土颗粒材料一般具有很强的各向异性特征[3]，在主应力方向旋转的加载条件，将出现主应力与主塑性应变率方向不共轴的现象[4]。这导致传统共轴理论无法准确预测应变局部化分叉的开始，非共轴理论可以较好地解决这一问题。早期，Rudnicki等[5]建立了具有角点结构屈服面的非共轴弹塑性本构模型，考虑了切向应变率的屈服效应。之后，Papamichos等[6]在此基础上给出了预测剪切带萌生和倾角的理论解。王兴等[7-8]提出了一种改进的角点理论并将其应用到了砂土状态相关剪胀模型中，该模型只在主应力方向改变的条件下产生非共轴塑性变形,克服了传统角点模型的不足。钱建固等[9]推导了有限变形条件下土体变形分叉的非共轴理论。Qian等[10-11]将Rudnicki等[5]的非共轴模型推广到了三维应力空间，考虑应力第三不变量的影响，在平面应变条件下对剪切带进行了分析，更好地预测了剪切带的萌生、倾角以及材料的剪胀性，随后对平面应变状态下土体的软化特性进行了模拟。吕玺琳等[12]采用非共轴理论对真三轴状态下砂土的分叉行为进行了分析，发现中主应力对分叉特性有显著影响。Du等[13]建立了黏土的三维非共轴本构模型，很好的描述了单调剪切加载下土体非共轴塑性流动特点。Yang等[14-16]模拟了简单剪切条件下非共轴特性对岩土颗粒材料的应力应变相应的影响，发现非共轴作用的发挥受硬化定律、流动法则及初始侧向压力等多种因素的影响，并将其非共轴模型应用到了一些结构问题的分析，如浅基础等，证明如果不考虑非共轴特性的影响，将导致不安全的设计这一共识[17-20]。

研究表明非共轴特性不单影响应变局部化分叉的开始，其对剪切带的后续发展演化也有影响。这是因为，当剪切带形成后，其内部发生了剧烈的剪切变形，是一种类似于简单剪切的变形，存在主应力方向旋转的情况，进而也存在非共轴现象[21-22]。已经证实采用传统共轴的弹塑性本构模型不能准确预测土体分叉的开始，建立非共轴本构模型是合适的。但是，基于经典连续体理论的非共轴本构模型，不包含任何内部长度参数，在模拟边值问题中的应变局部化现象时往往存在网格依赖性而不能很好地描述剪切带的宽度，更无法体现非共轴特性对剪切带发展演化的影响。为此，现建立基于微极理论的非共轴弹塑性本构模型，并自主开发相应的用户自定义单元(user-defined element,UEL)子程序，研究土体非共轴特性对剪切带的萌生、发展和演化这一渐进失效全过程的影响，克服传统本构模型的局限性，对进一步完善本土本构模型具有重要意义。

1 基于微极理论的非共轴本构模型

(1)

图1 非共轴塑性应变率的方向Fig.1 Direction of the non-coaxial plastic strain rate

(2)

共轴塑性应变率与传统塑性流动理论相同，用张量指标表示为

(3)

(4)

(5)

非共轴塑性应变率[6-7]可表示为

(6)

(7)

式中：hnc为非共轴塑性模量；sij为偏应力张量。

最终，非共轴塑性本构关系可以表示为

(8)

(9)

式中：hnc为非共轴塑性模量；Nijkl/hnc为非共轴项应变率的柔度矩阵；(De)ijkl为弹性本构张量。在经典连续体理论中，弹性本构关系表示为

(10)

(11)

微极理论中引入了物质点的微转动自由度ω，其几何方程除了位移和应变的关系外，还存在微转动位移和微曲率的关系[22]为

εij=uj,i-eijkωk，κij=ωj,i

(12)

式(12)中：κij为由微转动位移产生的微曲率张量。

同样，微极理论中的物理方程除了应力与应变关系之外，还存在与偶应力与微曲率的关系，其弹性本构关系表示为

(13)

(14)

式中：μij为偶应力张量；a、b和c为材料参数；G为传统剪切模量，G与a的乘积为Gc，称为微极剪切模量；l为内部长度参数,l与b的乘积、l与c的乘积亦可分别表示为lt和lb,分别代表与扭转和弯曲相关的内部长度尺度，在平面应变条件下参数c不起作用。

在应微极理论中，与非共轴相关的式(5)需要变为对应柯西应力和偶应力的两部分，即

(15)

最后，在微极理论框架下，基于D-P(Drucker-Prager)屈服准则和非关联流动法则，引入非共轴塑性应变率，形成了一个基于微极理论的非共轴弹塑性本构模型，对岩土材料应变局部化分叉和剪切带演化行为进行模拟。

2 应变局部化分叉理论

应变局部化条件建立在弱不连续性分叉理论的基础之上，认为声学张量行列式为零意味着分叉的开始[5]，即

det[Q]=det[nDn]=0

(16)

式(16)中：Q为声学张量；D为本构张量；n为剪切带的外法线方向矢量，n=(cosθ,sinθ,0)，θ为剪切带方向角。

在经典连续体理论中，对于平面应变问题，局部化条件可表示为

(17)

式(17)中：Qjk=niDijklnl，Qjk=n1D1jk1n1+n1D1jk2n2+n2D2jk1n1+n2D2jk2n2。

在微极连续体中，三维条件下，经推导Q扩展为一个6×6的矩阵，即

(18)

3 数值模拟

考虑到非共轴流动理论中，存在一个与屈服面相切的塑性应变率，采用完全隐式欧拉算法存在收敛困难的问题，因此采用自带误差控制的修正欧拉算法[23]。基于ABAQUS的UEL子程序二次开发，实现了以上本构模型的程序代码。

3.1 程序验证

为了验证程序的有效性，对Leighton Buzzard砂[24]进行了简单剪切试验的数值模拟，并与Roscoe的实验结果[4]进行了对比。简单剪切试验模拟如图2所示，采用了一个20节点减缩单元(C3D20R)，下边界固定约束，上边界节点施加水平的位移荷载，竖直向施加垂直压力134.45 kPa，侧向压力系数0.5，前后表面约束法向位移模拟平面应变加载条件。Leighton Buzzard砂的摩擦角38°，剪切模量100 MPa，泊松比0.3，采用基于D-P准则的理想弹塑性模型。

图2 简单剪切试验Fig.2 Simple shear test

图3为数值结果与Roscoe的实验结果的对比，从应力应变响应曲线和大主应力方向、大主塑性应变率方向的对比中可以看出，两者吻合较好，能够较好地反映Leighton Buzzard砂在主应力旋转的加载条件下的非共轴特性。

图3 数值结果与试验对比Fig.3 Comparison of the numerical results and the experiment results

3.2 平面应变压缩试验

为了研究非共轴特性对应变局部化的影响，对平面应变压缩试验进行了数值模拟。为了还原真实试验条件，建立10 cm×20 cm×2 cm的三维实体模型，采用C3D20R单元剖分网格，上下边界约束水平向自由度，施加大小相等方向相反的竖向位移荷载uy，前后边界约束法向位移，实现平面应变加载环境，左右边界施加围压100 kPa，如图4所示。计算中调用用户单元子程序(UEL子程序)。

图4 几何模型与边界条件Fig.4 Geometric model and boundary condition

模型采用的材料参数如表1所示[25]。

表1 材料参数Table 1 Material parameters

为了进一步验证程序的正确性，将Gc、hnc设为零，使程序退化为基于经典连续体的共轴本构模型，并将其与ABAQUS自带的程序进行了对比验证，结果如图5所示，其中图5(a)左侧图中的PEEQ代表ABAQUS自带程序计算得到的等效塑性应变，右侧图中的UVARM42是用户自定义程序得到的等效塑性应变。可见由等效塑性应变表征的剪切带模式和承载力位移曲线均体现了较好的一致性。

图5 退化为经典共轴本构模型的UEL子程序与ABAQUS自带程序的对比Fig.5 Comparison of the results obtained by using an UEL and ABAQUS built-in program

3.2.1 应变局部化分叉预测

如第3节所述，声学张量行列式为零意味着分叉的开始，图6和图7给出了分别基于微极连续体理论建立的共轴和非共轴本构模型所得到的一些典型区域的分叉情况。其中纵轴代表声学张量行列式与弹性声学张量行列式的比值，是一个归一化的结果，横轴为加载增量步。单元a、c和e位于剪切带内部，单元b、d位于剪切带两侧附近，单元f位于远离剪切带的位置。可以看出，在加载初期所有位置处归一化的声学张量行列式均为1，因为加载初期材料处于弹性阶段；随后纵坐标值出现了骤降，意味着加载进入塑性阶段，且数值最终趋于一个定值。共轴模型中曲线纵坐标趋近于一个非零常数，而非共轴模型则趋近于零。原因是非共轴模型中考虑了切向应变率，其弹塑性本构张量中增加了非共轴项，导致声学张量行列式最终趋于零而出现奇异性，这也是数值计算收敛困难的原因所在。同时观察到，不同位置处，分叉的时刻不同，最早出现分叉是在剪切带中心位置单元a，随后是位于剪切带内部的单元c和e，单元b和d对应曲线在加载进入塑性阶段之后出现了回弹，意味着剪切带出现后，其内部呈现继续加载状态，而外部则出现了卸载，单元f对应曲线一直保持为水平线，说明远离剪切带的位置材料一直处于弹性阶段。

图6 共轴模型模拟分叉Fig.6 Bifurcation in the coaxial model

图7 非共轴模型模拟分叉Fig.7 Bifurcation in the non-coaxial model

图8给出了由共轴模型和非共轴模型得到的结构承载力位移曲线，结果表明共轴模型得到的承载力峰值要高于非共轴模型，同时，分叉点均出现在峰值点之前，且共轴模型的分叉时刻早于非共轴模型，分叉时所对应的四种情况下剪切带模式如图9所示，随着非共轴塑性模量hnc的降低，非共轴程度加强，承载力峰值进一步降低，分叉时刻更滞后，这与目前的普遍认识一致[12,19-20]。这是因为共轴模型无法考虑由于屈服面切向的塑性应变率所引起的屈服效应，导致共轴模型预测的结果要比非共轴模型预测结果偏危险。

图8 共轴模型和非共轴模型模拟的分叉时刻对比Fig.8 Bifurcation in the coaxial and the non-coaxial models

图9 分叉时对应的等效塑性应变云图Fig.9 Distribution of the equivalent plastic strain when bifurcation occurs

3.2.2 剪切带宽度预测

图10～图12分别为根据共轴塑性应变、非共轴塑性应变和总体塑性应变计算得到的等效塑性应变。可以观察到，图10中共轴等效塑性应变随着非共轴程度的增强，其峰值逐渐降低，图11中非共轴等效塑性应变则随着非共轴程度的增强，其峰值逐渐增大，且在共轴模型中为零，但是注意到非共轴等塑性应变要比共轴等效塑性应变的值小一个数量级，最终导致图12中根据总体塑性应变计算得到的等效塑性应变随着非共轴程度的增强是逐渐减低的趋势。

图11 非共轴等效塑性应变云图Fig.11Distribution of the non-coaxial part of the equivalent plastic strain

图13、图14给出了微极理论中一些物理量的分布云图，其中图13为微转动位移，可以看出，随着非共轴程度的增强，剪切带内微转动位移峰值逐渐降低，转动变形趋于均匀。图14为微曲率的分布图，共轴模型中剪切带两侧微曲率呈现大小相等方向相反的规律，随着非共轴程度加强，其峰值也逐渐降低，且剪切带两侧的微曲率出现不对称性。这说明非共轴特性从某种程度上削弱了微极转动效应。为了验证基于微极理论非共轴本构模型在克服网格依赖性方面的有效性，对10×20、14×28、15×30以及16×32 4种网格密度下的剪切带宽度进行了比较(图15、图16，由等效塑性应变表征的剪切带)，可以看出随着网格加密，经典连续体理论模拟得到的剪切带出现了分支，宽度也随之改变，而基于微极理论得到的剪切带宽度几乎不变。

图13 微转动位移分布图Fig.13 The distribution of the micro rotation displacement in z direction

图14 微曲率κzx分布图Fig.14 The distribution of the micro curvature κzx

图15 网格敏感性调查(微极理论)Fig.15 Investigation of the mesh sensitivity(micropolar theory)

图16 网格敏感性调查(经典连续体理论)Fig.16 Investigation of the mesh sensitivity(classical continuum theory)

4 结论

本文建立了一个基于微极理论的非共轴弹塑性本构模型，并采用自带误差控制的修正的欧拉算法进行了数值实现。通过简单剪切试验数值结果与已有室内试验结果对比，验证了程序的正确性。重点研究了非共轴特性对应变局部化分叉预测和剪切带宽度的影响。结果表明非共轴模型预测的局部化分叉点滞后于共轴模型的预测结果，与已有文献结论一致，且随非共轴程度的增强，局部化变形趋于均匀，塑性区的范围扩大，剪切带宽度有变宽的趋势。等效塑性应变峰值呈现出随非共轴程度的增强而减小的趋势。可以看出本文建立的本构模型既可预测局部化分叉，又可体现非共轴特性对剪切带宽度的影响。综合比较下，基于微极理论的非共轴本构模型比传统的基于经典连续体理论的非共轴模型在模拟岩土材料边值问题中的应变局部化现象时更具优越性，可为深入探究岩土结构由应变局部化导致的渐进破坏的内在机制提供科学指导。