自然界的岩土材料一般均具有很强的各向异性性质,而横观各向同性则为较常见的一种体现形式.该类岩土介质的弹性行为通常采用5个弹性参数加以描述;其塑性行为和破坏机制与其他岩土材料也存在显著差异,即随着加载方向和层面倾角的变化,材料强度具有很强的各向异性^[1-2].为此,近年来建立了大量的各向异性强度准则^[3-4],很多学者认为岩土材料强度的各向异性与其微结构有密切关系,可以考虑将组构张量和加载方向及应力状态相联系,作为各向异性参数引入到各向同性强度准则中以模拟岩土材料的各向异性性质.Pietruszczak等提出了一个各向异性屈服准则^[4],给出了各向异性程度及加载方向对材料强度的影响的理论分析,其各向异性参数是应力和微结构张量联合不变量的显式函数,使得实验验证更易实现,该准则也被应用到了多种岩土材料及结构的分析中^[5].

岩土结构的破坏一般最先由微裂纹、微孔洞或软弱夹层等内部缺陷诱发局部变形,并逐渐发展到一个具有有限宽度的窄带,即剪切带.人们长期关注剪切带产生的条件^[6-7],以及影响剪切带发展演化的各种因素,并基于有限元软件作了大量的数值模拟.然而由于经典连续体理论不包含任何内部长度尺度,导致数值模拟中常常存在网格依赖性.为此人们发展了多种正则化方法,微极理论便为其中的一种，该理论又称为Cosserat连续体理论^[8],引入了独立的转动自由度和偶应力及相应的微曲率,使得本构方程中包含了一个内部长度参数,可以有效地解决网格依赖性问题.

本文的目的是发展适用于横观各向同性岩土材料的无网格依赖性的弹塑性有限元模型,从而对横观各向同性岩土结构的力学行为做出比较全面合理的预测与解释.对横观各向同性岩土材料多强调强度的各向异性^{[3-5, 9]},本文在文献^[9]的基础上,通过考虑弹性阶段的各向异性,建立了一个完整的横观各向同性的弹塑性模型，模拟了横观各向同性岩土结构的应变局部化现象,并对各向异性参数及材料主方向对剪切带及结构承载力的影响进行了分析.

1 基于Cosserat连续体的横观各向同性弹塑性模型 1.1 Cosserat连续体的横观各向同性弹性本构描述

对于平面应变问题,Cosserat连续体的弹性应力应变关系可以表示为

(1)

其中，

(2)

(3)

其中：κ_xz和κ_yz为微曲率；m_xz和m_yz为与微曲率共轭的偶应力，如图 1所示.

D_e为弹性本构矩阵,当采用如图 2所示局部坐标时,可表示为

(4)

式中:

其中:E₁为各向同性面方向的弹性模量；E₂为各向同性面法向的弹性模量；ν₁为各向同性面方向的泊松比；ν₂为各向同性面法向的泊松比,G₂为各向同性面法向的剪切模量；G_c为Cosserat剪切模量.

在图 2中,1o2为局部坐标系,xoy为整体坐标系.3轴垂直于1o2面,1o3面为各向同性面,2为其法线方向,则整体坐标系下的弹性本构矩阵为

(5)

其中R为转换矩阵:

(6)

1.2 Cosserat连续体的横观各向同性屈服准则

除了弹性阶段变形参数的各向异性外,材料进入塑性阶段,需考虑其强度准则的各向异性.Pietruszczak^[4]基于应力不变量与组构张量给出了一个各向异性材料强度参数的描述途径:

(7)

其中:η为与应力状态有关的强度参数；η₀为η的平均值;Ω_ij为组构张量的偏张量,描述η偏离η₀的程度;l_i为加载向量的各个分量.式(7)用Ω_ij的主值表示并考虑横观各向同性材料满足Ω₁=Ω₃,又由于Ω₁+Ω₂+Ω₃=0及l²₁+l²₂+l²₃=1,式(7)可进一步表示为

(8)

加载分量可按下式计算：

(9)

式中:表示应力张量,M=e⁽²⁾⊗e⁽²⁾,e⁽²⁾_i(i=1,2,3)表示2方向与整体坐标系下3个坐标轴夹角的余弦,如图 2所示:

(10)

将上述横观各向同性强度参数的描述方式与Drucker-Prager准则相结合，则可得到一个适用于描述横观各向同性岩土材料的扩展的Drucker-Prager准则：

(11)

其中：分别为偏应力第二不变量和应力第一不变量,在Cosserat连续体中分别表示为

(12)

(13)

A_φ,B_φ为内摩擦角φ和黏聚力c的函数,这里考虑内摩擦角的各向异性:

(14)

若Drucker-Prager准则为Mohr-Coulomb准则的内接圆,并考虑黏聚力为线性硬化或软化,则可将其分别表示为

(15)

其中:c₀为初始黏聚力;h_c^p为硬化(软化)参数;ε^p为等效塑性应变.

式(11),式(14),式(15)构成了适用于横观各向同性岩土材料的屈服准则.

对岩土材料一般采用非关联流动法则,塑性势函数及其参数取与屈服函数类似的形式,即

(16)

(17)

(18)

其中ψ为膨胀角.

记塑性流动方向为

(19)

式(19)最后两项体现了各向异性的影响.

2 数值算例

以ABAQUS软件为平台,开发UEL用户单元子程序,采用平面应变8节点等参元(CPE8R)4积分点减缩积分进行数值计算.计算模型如图 3所示.上下边界固定x方向位移,左右边界自由.上下边界施加对称的位移荷载.表 1为计算参数,其中弹性变形参数取自文献^[10],强度参数均取自黏土材料合理范围内,Cosserat参数G_c,l_c根据现有经验取值,β和Ω₁的值^[5]会在具体的算例中给出.

2.1 网格依赖性及数值收敛性检验

由于经典连续体中未引入任何内部长度参数,导致材料在伴随应变软化出现应变局部化现象时剪切带宽度病态地依赖于网格尺寸.数学角度的解释是,随着应变软化的出现,控制方程丧失了椭圆性.图 4给出了三种网格密度等效塑性应变分布图,等效塑性应变由UVARM15表示.可以看出随着网格加密，等效塑性应变峰值有所提高,但是剪切带的宽度基本保持不变.正是由于引入转动自由度和内部长度参数的原因,使得Cosserat连续体较好地解决了网格依赖性问题.数值结果显示本文在应力更新过程中所采用的向后欧拉算法具有较好的收敛性,能够达到预期的效果.

2.2 材料主方向及各向异性程度对剪切带及结构承载力的影响

图 5为β=30°,Ω₁=-0.05时,剪切带的形成过程演化图,其中u_y为上边界施加的位移.这里Ω₁取负值意味着材料各向同性面法向具有较高的强度^[5].可以看出剪切带率先在较弱的方向形成,之后逐渐发展并形成另外一支较弱的与之共轭的剪切带.

图 6为Ω₁=-0.05时等效塑性应变随材料主方向的变化云图,可见,β=0°和β=90°时两条共轭的剪切带基本呈对称分布,β=30°和β=60°时两条剪切带出现强弱之分,且随着β的增大等效塑性应变峰值逐渐增大,结构趋向于沿较弱方向最先发生破坏.剪切带倾角约为α=56°,如图 6d所示.图 7为β=30°时等效塑性应变随各向异性程度变化的云图,可以看出等效塑性应变随|Ω₁|的增大而减小,且弱剪切带对结构破坏的影响逐渐降低,使得强剪切带成为结构破坏的主要原因,说明各向异性程度的增加使结构破坏模式对方向更为敏感.

图 8为D_e和f均为横观各向同性时不同材料主方向下的承载力-位移曲线，若D_e为横观各向同性，f为各向同性，则无法体现结构强度的方向相关性；若D_e为各向同性，f为横观各向同性，则无法体现结构刚度的方向相关性；只有按照本文的D_e和f均为横观各向同性建立的弹塑性模型才能全面合理地模拟横观各向同性岩土结构强度和刚度随材料主方向的演化.

由图 9可以看出，β=0°时结构承载力随|Ω₁|的增大而增大,而β=90°时则随|Ω₁|的增大而降低.为了更好地描述结构承载力随材料主方向和各向异性程度的演化规律,图 10给出了归一化的曲线,可以看出整体上关于β=90°对称,当β∈[0,90°]时,承载力随β的增大而降低,随Ω₁的演化则存在一个反转点,约为β=45°,该点结构承载力与各向异性程度无关,该点之前承载力随|Ω₁|的增大而增大,该点之后则随|Ω₁|的增大而降低;当β∈[90°,180°]时,反转点为β=135°,该规律与文献^[5]中的一致,验证了程序的正确性.

3 结论

1) 材料主方向直接影响横观各向同性岩土结构的刚度和承载力,本文考虑弹性阶段本构矩阵的坐标旋转及塑性阶段屈服准则的各向异性而建立的弹塑性模型可以对横观各向同性岩土结构的刚度和承载力随层面倾角的演化做出较为合理的模拟.

2) 材料主方向和各向异性程度对应变局部化 模式有较大影响.对于β≠0°和β≠90°的情况,结 构趋向于沿较弱方向形成单一剪切带而破坏；对于β=0°或β=90°的情况,由于对称性,结构则趋向于形成两条共轭的成X形分布的剪切带破坏.

3) 由于Cosserat连续体的控制方程中包含了内部长度参数,能够对剪切带的发展演化及其宽度进行合理的预测.