横观各向同性岩土结构极限承载力的数值研究

删除或更新信息，请邮件至freekaoyan#163.com(#换成@)

本站小编 Free考研考试/2020-03-23

常江芳, 楚锡华, 徐远杰
武汉大学土木建筑工程学院, 湖北武汉 430072
收稿日期: 2015-03-16
基金项目: 国家自然科学基金资助项目(11172216,11472196).
作者简介: 常江芳(1988-),女,河北石家庄人,武汉大学博士研究生;
徐远杰(1956-),男,湖北武汉人,武汉大学教授。

摘要: 基于Cosserat连续体,建立了描述横观各向同性岩土材料的弹塑性模型,在该模型中,弹性本构关系由5个变形参数描述,塑性阶段通过引入组构张量及加载方向建立了修正的Drucker-Prager准则.针对上述模型推导了一致性映射返回算法的迭代格式及模量矩阵,结合适用于Cosserat连续体的平面应变8节点减缩单元应用ABAQUS用户自定义单元接口UEL进行了数值实现.数值算例分析了各向异性程度和材料主方向对横观各向同性岩土结构的极限承载力和变形局部化模式的影响.结果表明,所建立模型能较好地模拟横观各向同性岩土结构的应变局部化行为,且较好地解决了伴随应变局部化现象出现的网格依赖性问题.
关键词：岩土结构横观各向同性Cosserat连续体应变局部化极限承载力数值模拟
Numerical Investigation of Ultimate Bearing Capacity for Transversely Isotropic Geostructures
CHANG Jiang-fang, CHU Xi-hua, XU Yuan-jie
School of Civil and Architectural Engineering, Wuhan University, Wuhan 430072, China
Corresponding author: CHANG Jiang-fang, E-mail: cjf881024@163.com
Abstract: An elastic-plastic model for transversely isotropic geostructures is developed based on the Cosserat continuum. In the model, the elastic constitutive relationship is described by 5 deformation parameters and the Drucker-Prager yield criterion is extended by introducing fabric tensor and loading direction in the plastic stage. The iterative format for return mapping algorithms and the tangent modulus matrix are formulated for the proposed model. The numerical simulation is implemented based on the user subroutine (UEL) of ABAQUS, and a plane strain 8-noded reduced integrated element is used. The influence of the material principal direction and the anisotropic degree on the strain localization and the bearing capacity of the structure are analyzed. Numerical results show the validity and performance of the proposed model in simulating the strain localization behavior of transversely isotropic geostructures. Furthermore, the mesh dependency accompanied with strain localization is effectively solved.
Key Words: geostructurestransversely isotropyCosserat continuumstrain localizationultimate bearing capacitynumerical simulation
自然界的岩土材料一般均具有很强的各向异性性质,而横观各向同性则为较常见的一种体现形式.该类岩土介质的弹性行为通常采用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所示.
图 1(Fig. 1)

图 1 二维Cosserat连续体应力和偶应力示意图Fig.1 Stress and couple-stress in a two-dimensional Cosserat continuum

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

(4)

式中:

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

图 2 横观各向同性材料主方向示意图Fig.2 Schematic of material principal orientation

在图 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]会在具体的算例中给出.
图 3(Fig. 3)

图 3 计算模型示意图Fig.3 Sketch of computing model

表 1(Table 1)

表 1 材料参数Table 1 Material parameters

E₁	ν1	E₂	ν₂	G₂	φ₀	ψ₀	c₀	h_p	β	Ω₁	G_c	l_c
MPa	ν1	MPa	ν₂	MPa	(°)	(°)	kPa	kPa	(°)	Ω₁	MPa	m
30.5	0.139	16.8	0.146	8.92	20	18	11.4	-30	-	-	4.46	0.01

表 1 材料参数 Table 1 Material parameters

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

图 4 不同网格密度下等效塑性应变分布图Fig.4 Distribution of equivalent plastic strain with different mesh densities (a)-8×16； (b)-10×20; (c)-15×30.

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

图 5 剪切带随加载的演化图Fig.5 Variation of shear band with the increasing displacement (a)-u_y=6.6mm； (b)-u_y=8.4mm; (c)-u_y=12mm.

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

图 6 不同材料主方向下等效塑性应变分布图Fig.6 Distribution of equivalent plastic strain with different material principal directions (a)-β=0°； (b)-β=30°;(c)-β=60°； (d)-β=90°.

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

图 7 不同各向异性程度下等效塑性应变分布图(β=30°)Fig.7 Distribution of equivalent plastic strain with different anisotropic degrees(β=30°) (a)-Ω₁=0； (b)-Ω₁=-0.05;(c)-Ω₁=-0.1; (d)-Ω₁=-0.15.

图 8(Fig. 8)

图 8 承载力-位移曲线随材料主方向的演化Fig.8 Variation of bearing capacity-displacement with material principal direction

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

图 9 承载力-位移曲线随各向异性程度的演化Fig.9 Variation of bearing capacity-displacement with different anisotropic degree (a)-β=0°； (b)-β=90°.

图 10(Fig. 10)

图 10 承载力随主方向β及各向异性程度Ω₁Fig.10 Variation of bearing capacity with β and Ω₁

3 结论1) 材料主方向直接影响横观各向同性岩土结构的刚度和承载力,本文考虑弹性阶段本构矩阵的坐标旋转及塑性阶段屈服准则的各向异性而建立的弹塑性模型可以对横观各向同性岩土结构的刚度和承载力随层面倾角的演化做出较为合理的模拟.
2) 材料主方向和各向异性程度对应变局部化模式有较大影响.对于β≠0°和β≠90°的情况,结构趋向于沿较弱方向形成单一剪切带而破坏；对于β=0°或β=90°的情况,由于对称性,结构则趋向于形成两条共轭的成X形分布的剪切带破坏.
3) 由于Cosserat连续体的控制方程中包含了内部长度参数,能够对剪切带的发展演化及其宽度进行合理的预测.
参考文献

[1]	Duncan J M, Seed H B. Strength variation along failure surfaces in clay[J].Journal of Soil Mechanics&Foundations Division, 1966, 92(6): 81–104.(0)
[2]	Nova R. The failure of transversely isotropic rocks in triaxial compression[J].International Journal of Rock Mechanics and Mining Sciences&Geomechanics Abstracts, 1980, 17(6): 325–332.(0)
[3]	Gao Z W, Zhao J D, Yao Y P. A generalized anisotropic failure criterion for geomaterials[J].International Journal of Solids and Structures, 2010, 47(22): 3166–3185.(0)
[4]	Pietruszczak S, Mroz Z. Formulation of anisotropic failure criteria incorporating a microstructure tensor[J].Computers and Geotechnics, 2000, 26(2): 105–112.(0)
[5]	Azami A, Pietruszczak S, Guo P. Bearing capacity of shallow foundations in transversely isotropic granular media[J].International Journal for Numerical and Analytical Methods in Geomechanics, 2010, 34(8): 771–793.(0)
[6]	Rice J R, Rudnicki J W. A note on some features of the theory of localization of deformation[J].International Journal of Solids and Structures, 1980, 16(7): 597–605.(0)
[7]	Dietsche A, Steinmann P, Willam K. Micropolar elastoplasticity and its role in localization[J].International Journal of Plasticity, 1993, 9(7): 813–831.(0)
[8]	Cosserat E, Cosserat F. Théorie des corps déformables[M].Paris: Hermann, 1909.(0)
[9]	Chang J F, Chu X H, Xu Y J. Finite-element analysis of failure in transversely isotropic geomaterials[J].International Journal of Geomechanics, 2014, 15(6): 04014096.(0)
[10]	金耀华, 孙炎, 钱玉林. 水平条形均布荷载下横观各向同性弹性地基的位移分析[J].岩土力学, 2008, 26(sup 1): 163–167. ( Jin Yao-hua, Sun Yan, Qian Yu-lin. Analysis of displacements of transversely isotropic foundation under level strip uniformly distributed load[J].Rock and Soil Mechanics, 2008, 26(sup 1): 163–167.)(0)