0.   引言

单剪试验被广泛应用于评价土体的剪切强度,和模拟现场加载条件下的力学响应（如坡度,打桩,滑坡,地震,顶管隧道）。然而,由于单剪过程中某些力学特性无法直接测量,比如确定试样剪切带位置、侧向边界的径向应力、主应力/应变旋转程度等^[1-2],均是单剪试验研究所亟需解决的问题。此外,与理想纯剪试验相比,单剪试验侧向边界不存在互补的剪切应力,这将导致圆柱形试样的应力不均匀性^[3],并使试验研究变的复杂。因此,建议提出一种有效的方法来评估单剪过程中试样的不均匀性。

有限元法（FEM）是研究土样应力/应变分布的理想分析工具。已经有研究者利用二维模拟再现了单剪试验中的应力不均匀特性^[4-5]。然而,由于二维平面应变假设只适用于立方体试样,因此只能提供有限的结果和关于边界效应的讨论。最近,Doherty等^[6]使用修正剑桥模型提出了一种更先进的三维有限元建模方法。但他的是采用高岭土力学参数,建模中忽略了土体的各向异性。

本文利用考虑初始各向异性的改进的SIMSAND模型和ABAQUS大变形计算模块,提出了一种模拟单剪试验的数值方法。首先将改进的SIMSAND模拟写入ABAQUS/explicit自定义用户子程序vumat来实现有限元的计算,然后基于GDS型单剪试验,建立相应等尺寸的三维计算模型模拟一系列枫丹白露砂试验,最后通过模拟结果研究试样不均匀性的边界效应。

1.   考虑各向异性的本构模型改进说明

本文基于一种简单的临界状态砂土模型^[7-11]（SIMSAND模型）,其基本本构方程如表1所示,综合考虑非线性弹性、非线性塑性硬化法则、修正剪胀剪缩特征、三维强度和临界状态等特征,来有效地描述砂土特征。

表  1  SIMSAND模型的基本本构方程

Table  1.  Basic constitutive equations of SIMSAND

组成部分	本构方程
弹性准则	${\dot{\epsilon }}_{ij}^{\text{e}}=\frac{1+\upsilon }{3K(1-2\upsilon )}{{\sigma }^{\prime }}_{ij}-\frac{\upsilon }{3K(1-2\upsilon )}{{\sigma }^{\prime }}_{kk}{\delta }_{ij}$ $K=K_0\cdot p_{\text{at}}\frac{{(2.97-e)}^2}{(1+e)}{\left(\frac{p^{\prime }}{p_{\mathrm{at}}}\right)}^n$
屈服面	$f=\frac{q}{p^{\prime }}-H$
塑性势面	$\frac{\partial g}{\partial p^{\prime }}=A_{\text{d}}\left(M_{\text{pt}}-\frac{q}{p^{\prime }}\right),\frac{\partial g}{\partial q}=1$
硬化参数	$H=\frac{M_{\text{p}}{\epsilon }_{\text{d}}^{\text{p}}}{k_{\text{p}}+{\epsilon }_{\text{d}}^{\text{p}}}$
临界状态参数	$e_{\text{c}}=e_{\text{ref}}-\lambda {\left(\frac{p^{\prime }}{p_{\text{at}}}\right)}^{\xi }$ $\tan {\varphi }_{\text{p}}={\left(\frac{e_{\text{c}}}{e}\right)}^{n_{\text{p}}}\tan {\varphi }_{\mu }\tan {\varphi }_{\text{pt}}={\left(\frac{e_{\text{c}}}{e}\right)}^{-n_{\text{d}}}\tan {\varphi }_{\mu }$
三维强度标准	$\begin{array}{l}{M}_{\text{p}}=\frac{6\sin {\varphi }_p}{3-\sin {\varphi }_p}{\left[\frac{2c_1^4}{1+c_1^4+(1-c_1^4)\sin 3\theta }\right]}^{\frac{1}{4}},其中c_1=\frac{3-\sin {\varphi }_{\text{p}}}{3+\sin {\varphi }_{\text{p}}}\\ M_{\text{pt}}=\frac{6\sin {\varphi }_{\text{pt}}}{3-\sin {\varphi }_{\text{pt}}}{\left[\frac{2c_2^4}{1+c_2{}^4+(1-c_2^4)\sin 3\theta }\right]}^{\frac{1}{4}},其中c_2=\frac{3-\sin {\varphi }_{\text{pt}}}{3+\sin {\varphi }_{\text{pt}}}\end{array}$

下载: 导出CSV 

| 显示表格

组构各向异性已被广泛证实会影响颗粒材料的强度^[12-14]。在单剪条件下,引入组构各向异性联合不变量,可将原始各向同性强度扩展到砂的各向异性强度。在这种方法中,利用由Oda等^[15]提出的组构张量F_ij来描述组构各向异性的程度：

式中,Δ为描述颗粒材料组构各向异性分布的参数,0<Δ<1为典型砂土的典型范围。

根据张量理论^[16]采用联合不变量A,将组构张量的偏差部分与偏应力张量联系起来：

式中,s_ij为偏应力张量(s_ij = σ_ij-p′δ_ij),d_ij为组构张量的偏差部分(d_ij= F_ij-F_kkδ_ij/3),δ_ij为Kronecker delta。

结合各向异性强度的方式类似于Li等^[17]、Gao等^[14]的方法,是将各向异性的影响函数引入到强度准则中。各向异性插值函数由联合不变量A和组构各向异性参数Δ构成。利用各向异性修正函数修正原始各向同性峰值强度M_p和相变强度M_pt,

式中,$g\text{(}A,\Delta \text{)}$为各向异性修正函数,该函数采用与Pietruszczak等^[18-19]的研究相同的形式,

其中,c₁为材料参数。对于各向同性的砂样,Δ= 0,g(A,Δ) ≡ 0。对于0<Δ<1的组构各向同性砂样,g(A, Δ)随A和Δ的改变而改变。在大多数方法中,组构各向异性的程度可能不便于测量。为了解决这个困难,组合项c₁被认为是一个单一的材料参数,因此只需要校准而不是确切的值。

另外,Vaid等^[20]的三轴试验和单剪试验结果的比较表明,剪切刚度退化的主要原因是由于单剪阶段,主应力旋转过程中产生的应力各向异性。Yang等^[21-23]通过引入非共轴塑性应变模拟了这种退化。参考各向异性影响函数$g\text{(}A,\Delta \text{)}$修正强度准则的方法^[14],将各向异性影响函数$g\text{(}A,\Delta \text{)}$式（4））引入到塑性硬化方程中,即对控制$q/p-{\epsilon }_{\text{d}}^{\text{p}}$曲线斜率的参数k_p进行修正,来实现单剪过程中由主应力轴偏转引起的剪切刚度弱化现象,如各向异性校正函数如下所示：

式中 对于各向同性的砂样,$\Delta$=0,$g\text{(}A,\Delta \text{)}$≡ ${\tilde{c}}_2\cdot (1+$$A)$≡0,${k^{\prime }}_{\text{p}}=k_{\text{p}}\cdot \exp \left[{\tilde{c}}_2\cdot (1+A)\right]$=0；对于组构各向同性砂样,0<Δ<1,${\tilde{c}}_2\cdot (1+A)$≠0,其中${\tilde{c}}_2$是材料参数,其数值大小可控制剪切刚度${k^{\prime }}_{\text{p}}$的变化程度。

根据ABAQUS/Explicit所提供的用户材料子程序接口VUMAT^[24],将SIMSAND模型写入有限元计算程序,其基本流程如图1所示。

图  1  ABAQUS/Explicit分析的基本流程

Figure  1.  Flow chart of explicit finite element analysis based on ABAQUS/Explicit

下载: 全尺寸图片 幻灯片

ABAQUS/Explicit计算程序中,首先进行节点的位移计算,通过内力矢量I和外力矢量P的动力平衡方程,求解节点的在t时刻下的加速度${\ddot{u}}_{\left(t\right)}$：

节点的位移u可以通过时间中心积分法求得

式中,Δt为时间增量,其中单元的应变增量$\Delta \epsilon$可通过节点的位移增量$\Delta u=\Delta t_{(t+\Delta t)}{\dot{u}}_{\left(t+\frac{\Delta t}{2}\right)}^{}$算得。

然后通过用户材料子程序VUMAT（SIMSAND本构关系）,采用切面算法^{[8, 25-26]}计算单元的应力增量Δσ,如图2所示,其中塑性算子为

图  2  切面法塑性修正迭代计算

Figure  2.  Schematic diagram of a general cutting plane algorithm

下载: 全尺寸图片 幻灯片

式中,f为屈服函数,g为塑性势函数,D为弹性矩阵,κ^*为硬化参数,ε^p为属性应变。

假设应变增量Δε完全为弹性变形,弹性试算后的应力点为${\sigma }_n$。当试算点${\overline{\sigma }}_{n+1}$在屈服面外时f >0,此时进入塑性修正,此过程需要跟新应力点σ和硬化参数κ^*,

式中,$\partial {\epsilon }_{}^{\text{p}}$为屈服面向外发展过程中得到的塑性应变增量。若此时应力点第一次塑性修正后依然没与屈服面重合,需要进行再次塑性修正,跟新应力点和硬化参数,屈服面向外扩展。这个过程不断重复直至修正后的应力点和屈服面重合。

2.   单剪试验的有限元分析模型

2.1   室内单剪试验

为了使用改进模型验证有限元方法,使用商用GDS单剪设备进行了一系列单剪试验,其位移传感器,量程为50 mm,分辨率可达0.001 mm。可选的荷载范围为5 kN或10 kN,其荷重传感器精度可达到0.1%的传感器量程。试验材料是枫丹白露标准砂（NE34）,为法国岩土工程应用的参考材料,其颗粒的平均粒径为0.21 mm,级配如图3所示,不均匀性系数为1.53,最小孔隙比为0.510,最大孔隙比为为0.882,密度为2.65 g/cm³。伺服控制系统能够在水平和垂直两个方向上传导应力或应变以控制加载路径。图3总结了该装置的基本特征：①圆柱形试样的直径D₀ = 70 mm,高度H₀ = 25 mm（高宽比H₀/D₀ = 0.36）；②将试样密封在厚度为0.2 mm的软丁基薄膜中,其自身置于具有70 mm内径的刚性圆形聚四氟乙烯涂覆的薄环（每个1 mm）上,保持恒定的横截面积但允许单剪变形；③在试样的顶部和底部,砂与粗烧结不锈钢板接触,以防止基座和试样之间的任何滑动。

图  3  GDS单剪试验

Figure  3.  Simple shear tests

下载: 全尺寸图片 幻灯片

2.2   有限元单剪模型

如表2所示,进行了一系列单调加载恒定体积和恒定正应力单剪试验。在K₀固结后,通过将试样的高度限制在恒定值来确保体积恒定条件（相当于“不排水条件”）。值得注意的是,所有试验均在干燥条件下完成。同时,恒定法向应力条件（相当于“排水条件”）通过固定K₀固结后试样顶部的法向应力来实现。所有的试验结果都被用于有限元分析。

表  2  枫丹白露砂单调加载的单剪试验列表

Table  2.  Summary of monotonic simple shear tests on fontainebleau sand

编号	加载方式	e0(初始)	e′(K0固结后)	${{\sigma }^{\prime }}_{\text{n0}}$ /kPa
S-1	法向应力恒定	0.70	0.691	104
S-2	法向应力恒定	0.70	0.688	208
S-3	法向应力恒定	0.70	0.678	416
S-4	法向高度恒定	0.68	0.666	208
S-5	法向高度恒定	0.68	0.654	416

下载: 导出CSV 

| 显示表格

图4为单剪试验的三维模型。即该模型的尺寸与图4所示GDS单剪设备的尺寸相同。由于对称性,只建立一半的圆柱体试样模型。如图4所示,在ABAQUS中通过C3D8R单元对圆柱体试样进行划分,共有11000个单元。金属环和顶/底板是圆柱柱试样的边界,用刚体模拟。金属环的内径为70 mm,厚度为1 mm,而顶板和底板的直径均为70 mm,厚度为2.5 mm。总共有29个刚性环堆叠用以限制试样的横向位移,并且使用两个刚性（顶部和底部）平板来限制其竖向位移。在圆柱体试样和刚性边界（环和摩擦基座）之间施加无摩擦的接触。试样的顶部和底部表面分别连接到顶部和底部刚性板,以防止它们和试样之间发生任何滑动。

图  4  单剪试验的三维有限元建模

Figure  4.  Three-dimensional FEM model for simple shear tests

下载: 全尺寸图片 幻灯片

单剪试验的模拟过程分为K₀固结和剪切两个过程。在第一步（K₀固结）中,分别在顶板上施加等于0.2,0.4,0.8 kN的力以产生分别对应于104,208,416 kPa的初始法向应力${{\sigma }^{\prime }}_{\text{n0}}$。然后,通过对底板施加0.05 mm/s的位移来模拟剪切过程,同时将顶板在x,y方向上的位移固定为零。整个模拟过程花了100 s的时间产生了20%的剪应变。模拟的边界条件被认为与基于物理实验室的GDS简单剪切试验类似。恒定法向应力简单剪切试验通过保持顶板上的恒定法向力来模拟；通过在剪切过程中将试样保持在恒定高度来模拟恒定体积简单剪切试验。

位移和力的测量方法与实际室内试验相似。有效法向应力${{\sigma }^{\prime }}_{\text{n}}$和剪应力$\tau$分别通过校准底板上测得的垂直和水平力来获得。体积应变${\epsilon }_{\text{v}}$=ΔH/H₀,通过测量圆柱试样高度的变化获得（其中ΔH为试样高度的改变量,H₀为剪切开始时的试样高度）。剪应变γ= Δd/H₀（其中Δd是水平位移）。分别监测x轴（Path-1）,y轴（Path-2）和25 mm半径圆（Path-3）的应力状态,来评估剪切过程中的应力不均匀性,如图4（c）,（d）所示。

2.3   模型参数确定

为验证SINSAND模型的可应用性,选用了法国常用的枫丹白露标准砂的系列三轴试验进行模拟^[27-28]。砂土的模型参数：K₀=100,$\nu$=0.25,n=0.51,${\varphi }_{\text{c}}$=33.2,e_ref=0.811,$\lambda$=0.055,$\xi$=0.46,k_p =0.0022,A_d =0.39,n_p =1.9,n_d =4,${\tilde{c}}_1$=0.18,${\tilde{c}}_1$=2.0。其中弹性相关参数通过拟合各向同性压缩试验得到,临界状态和塑性颗粒咬合效应参数均采用Jin等^[29-31],提出的参数优化识别方法获得。具体标定如图5所示。

图  5  用于参数标定的枫丹白露砂三轴试验

Figure  5.  Simulated results of triaxial tests using determined parameters for Fontainebleau NE34 sand

下载: 全尺寸图片 幻灯片

通过有限元模拟初始状态（${{\sigma }^{\prime }}_{\text{n0}}$= 104 kPa,e₀ = 0.7,S-1）为恒定正应力条件下的单剪试验来确定枫丹白露砂的各向异性参数（${\tilde{c}}_1$,${\tilde{c}}_2$）。此外,还进行了类似的没有考虑各向异性的模拟（即${\tilde{c}}_1$=0和${\tilde{c}}_2$= 0）,以便强调在考虑参数${\tilde{c}}_1$,${\tilde{c}}_2$时模型的性能。图6为模拟结果与试验结果的比较,结果表明引入参数${\tilde{c}}_1$= 0.18和${\tilde{c}}_2$= 2.0的模拟结果与试验结果一致。因此,参数（${\tilde{c}}_1$= 0.18和${\tilde{c}}_2$= 2.0）可用于后续的单剪试验模拟。

图  6  用于参数标定的枫丹白露砂单剪试验

Figure  6.  Calibration of anisotropic parameters

下载: 全尺寸图片 幻灯片

3.   试验验证

为了验证改进模型的有限元分析的性能,基于单剪试验的足尺三维模型（图4）,分别模拟了表2中所示的4个附加的单剪试验。图7为模拟结果和试验结果的比较,其模拟曲线中的数值波动是由有限元显式计算的应力波传递所造成。恒定体积试验是通过固定试样高度来实现的,即ΔH = 0（试样高度变化为零）,Δr = 0（试样径向变化为零）。在恒定体积剪切的初始阶段,颗粒重新排列出现压缩密实（即剪缩）过程,此时由于体积应变恒定,砂样首先出现法向应力变小的趋势,随后当应力比达到相变状态时,土体开始发生剪胀,法向应力逐渐增加并达到最终的极限状态。因此可以基于足尺三维单剪模型来获得恒定正应力和恒定体积条件下的单调加载特性（剪缩或剪胀）。由于室内直接试验研究困难^[3],可以结合有限元模拟结果来分析试样的应力不均匀性。

图  7  三维单剪试验模拟结果

Figure  7.  Simulated results by three-dimensional simple shear tests

下载: 全尺寸图片 幻灯片

图8（a）为在恒定正应力单剪试验（S-3）下试样顶部表面沿着三条路径的垂直应力分布的演变。这3条路径先前在图4中进行了描述,这些路径用于评估沿x轴（路径-1）,y轴（路径-2）和半径为25 mm的圆（路径-3）的应力不均匀性。在剪切过程中,应力不均匀性随剪切应变水平的增加而增加,最低和最高垂直应力位于侧向（右侧和左侧）边界。结果还表明,垂直应力在试样的中心区域几乎可以保持恒定。

图  8  垂直有效应力不均匀分布

Figure  8.  Nonuniform distribution of vertical effective stress

下载: 全尺寸图片 幻灯片

恒定体积试验的垂直应力分布（S-5）的结果也在图8（b）中给出。与恒定法向应力条件不同,由于恒定体积的边界条件,垂直应力首先会由于低剪切应变水平相对减小（剪缩）,然后增加（剪胀）。与恒定法向应力试验相比,可以发现垂直分布的类似分布特征,最低和最高值也位于横向边界中并且相当均匀地位于中心区域中。为了研究侧向边界互补剪应力缺失引起的应力不均匀特性,根据网格化之后的有限元模型汇总了11000个高斯点。图9为对于恒定正应力和恒定体积单剪试验在剪切过程中基于有限元中的所有高斯点的竖向应力${{\sigma }^{\prime }}_z$的正态分布：

图  9  恒法向应力试验的垂直应力的正态分布（S-3：${{\sigma }^{\prime }}_{\text{n}}$=416 kPa）

Figure  9.  Normal distribution of vertical stress under constant normal stress tests

下载: 全尺寸图片 幻灯片

式中,f(x)为可能的密度函数,μ为试样x分布的平均至；$\sigma$为标准差,其中试样x对应竖向应力${{\sigma }^{\prime }}_z$。

对于恒定正应力单剪试验,平均值μ对应于平均垂直应力${{\sigma }^{\prime }}_z$,其保持稳定并接近初始法向应力（即${{\sigma }^{\prime }}_n$ _{=416 kPa）。偏差$\sigma$对应于应力不均匀程度,随着剪切应变的增加,应力逐渐增大,如图9所示。另外,如图10所示,也可以在恒定体积单剪试验中找到类似的应力不均匀性演变（标准偏差随剪切应变$\gamma$增加）。}

图  10  恒体积试验的垂直应力的正态分布（S-5：${{\sigma }^{\prime }}_{\text{n}}$=416 kPa）

Figure  10.  Normal distribution of vertical stress under constant volum tests

下载: 全尺寸图片 幻灯片

4.   应力应变不均匀性分析

根据对恒定法向应力单剪试验（S-3）的模拟结果,在图11中绘制了剪切应力$\tau$的演变（相应于有限元模拟中的${{\sigma }^{\prime }}_{xz}$）,塑性偏应变${\epsilon }_{\text{d}}^{\text{p}}$和各向异性变量A。对于早期的剪切阶段（应变$\gamma$小于或等于1%）,试样中的剪切应力、塑性偏应变和各向异性变量A相当均匀,特别是在中央或中部区域。对于应变大于5%的剪切阶段,剪应力峰值逐渐从左上角和右下角扩展到中心区。在试样中可以找到塑性偏应变的对角线区域。此外,由于侧边界和试样之间没有摩擦力,在剪切过程中,侧向边界处的剪切应力和偏应力几乎为零,这可能是导致应力/应变不均匀性根本的原因,从而最终导致逐渐破坏。这些模拟当中的剪应力和塑性应变的对角线累积区域的趋势与Budhu^[32]的射线照片所记录的实际物理破坏区域相似。

图  11  恒定法向应力为416 kPa试验的连续单剪过程

Figure  11.  Profiles of successive simple shearing process under constant normal stress at ${{\sigma }^{\prime }}_{\text{n}}$= 416 kPa

下载: 全尺寸图片 幻灯片

图12绘制了恒定体积单剪（S-5）试验后的剪切应力$\tau$,塑性偏应变${\epsilon }_{\text{d}}^{\text{p}}$和各向异性变量A。类似于通过恒定法向应力试验获得的图形（S-3）,从左上角到右下角形成了累积塑性偏应变的对角区域。由于试样的高度恒定的边界条件,对角区塑性偏应变的积累比恒定正应力条件下更明显。如图11（c）和图12（c）所示,各向异性变量A也沿试样的对角线区域积聚。

图  12  法向应力为416 kPa的恒定体积试验的连续单剪过程

Figure  12.  Profiles of successive simple shearing process under constant volume at ${{\sigma }^{\prime }}_{\text{n}}$ = 416 kPa

下载: 全尺寸图片 幻灯片

图1给出了有限元建模中所有高斯点的3个模拟状态变量（$\tau$,${\epsilon }_{\text{d}}^{\text{p}}$和A）的不均匀性特性的概率分析。可以看到,正态分布的演变对于剪应力和塑性偏应变是相似的,因为偏差随应变水平而增大,对应于剪切应力或塑性应变的不均匀性在大大增加时。此外,对于如图13,14所示的各向异性变量A,分布的偏差在低剪切应变水平下变化很大,如$\gamma$=0.5%,这意味着对剪应变的增加不敏感（超过1%）。

图  13  基于恒定法向应力试验过程中的概率分布（S-3：${{\sigma }^{\prime }}_{\text{n}}$=416 kPa）

Figure  13.  Probability distribution with increase of shear strain based on constant normal stress

下载: 全尺寸图片 幻灯片

图  14  恒定体积试验过程中的概率分布（S-5：${{\sigma }^{\prime }}_{\text{n}}$=416 kPa）

Figure  14.  Probability distribution based on constant volume

下载: 全尺寸图片 幻灯片

5.   结语

本文考虑砂土初始各向异性的影响,提出了一个模拟全尺寸单剪试验的数值模拟方法。

该方法能引入横观各向异性联合不变量的参数,将模型从原始各向同性强度扩展到各向异性强度,由此获得了单剪条件下固有的各向异性特性。

考虑组构各向异性的SIMSAND模型子程序写入有限元模型,通过建立等尺寸的三维有限元模型,通过建立等尺寸的三维有限元模型,可以很好的模拟枫丹白露砂的GDS单剪试验。

本研究可以提高对单剪条件下试样剪切特性的理解,为分析试样剪切过程中应力应变不均匀性提供有效的计算工具。