Spherical cavity expansion-based method for cone factor considering nonlinear characteristics of clay
-
摘要: 当前规范中根据静力触探结果解释黏土不排水剪切强度采用的是推荐的锥尖系数,但该系数存在经验性取值的问题,更重要的是当前锥尖系数均是基于理想弹塑性土体模型导得,未考虑土体应变应变的非线性特性。为此,首先在ABAQUS有限元软件中模拟了圆锥贯入过程,评估锥尖阻力qc,饱和黏土采用双曲线硬化弹塑性模型,并使用ALE重划分技术防止网格发生畸变。有限元结果表明:土体刚度、锥面粗糙度和破坏准则等因素均影响着锥尖系数。同时在有限元中模拟了不排水圆孔扩张过程,得到球孔扩张极限承载系数的表达式。随后建立了锥尖阻力与球孔极限扩孔压力之间的等效换算关系,最后对工程实例进行了计算分析。Abstract: According to the current specifications, the undrained shear strength of clay is interpreted based on the recommended cone factor derived from empirical values. However, this factor has limitations in terms of its empirical nature. Moreover, the existing cone factors are derived from the elastic-perfectly plastic model, without considering the nonlinear characteristics of soil deformation. To tackle these issues, a simulation of the cone penetration process is performed using the finite element software ABAQUS. The cone resistance qc is assessed using the hyperbolic hardening model for clay, and the arbitrary Lagrangian Eulerian (ALE) remeshing technique is employed to prevent grid distortion. The finite element results show that the factors such as soil stiffness, cone roughness and failure criteria affect the cone factor. Furthermore, the undrained expansion for a circular cavity is simulated using the finite element method to obtain the expression for the ultimate bearing capacity factor. Subsequently, an equivalent conversion relationship between the smooth cone tip resistance and the limit expansion pressure of the cavity is established. Finally, a verification against an engineering case is conducted.
-
0. 引言
静力触探(CPT)广泛运用于陆上及海上岩土工程勘察,黏土不排水剪切强度可根据锥尖阻力进行确定,即
su=qc−σv0Nk。 (1) 式中:σv0为当前深度上覆总压力,Nk为锥尖承载系数。正常固结黏土中锥尖系数通常为11~19,平均约为15。锥尖系数直接影响不排水强度评估的可靠性。工程中有时根据黏土塑性指数估算Nk,但地域经验性太强,且没有统一的结论。
Yu等[1]总结了承载力理论、圆孔扩张理论、应变路径法、大变形有限元等不同方法推导的锥尖系数表达式,其中最重要的特征是Nk与土体刚度系数G/su相关,其次受锥面粗糙度、初始应力各向异性等因素影响。其中,所有Nk与G/su的关系式均是基于理想弹塑性本构假设,尽管文献[2,3]也曾考虑应力应变非线性,但并未直接给出计算公式,而是侧重分析超固结比、应变软化特性和孔隙水压力的消散等特性。
在CPT贯入问题分析中,土体应力应变的非线性可能对锥尖贯入阻力qc有不可忽略的影响。为此本文采用双曲线硬化模型描述饱和黏土非线性应力应变特性。首先开展大变形有限元数值计算,分析了土体刚度系数、锥面粗糙度等因素对锥尖系数的影响关系。进一步地,提出椭圆形锥尖承载机构将锥尖贯入阻力与球孔扩张极限扩孔压力建立联系,给出锥尖系数的计算式。与有限元和一工程实例进行对比验证,显示了本文方法的适用性。
1. 锥尖贯入ALE大变形有限元计算
1.1 ABAQUS中双曲线硬化模型定义
在商业化有限元软件ABAQUS中模拟圆锥贯入过程。选用内置的广义塑性模型,其认为材料服从Mises屈服准则、各向同性硬化及相关联流动法则。通过定义偏应力q-塑性应变εpq关系曲线从而确定土体硬化关系,双曲线型硬化关系表示为
q=εq1E+εqquRf。 (2) 式中:E为土体初始弹性模量,对不排水分析,E=3G;qu为极限偏应力值,qu=2su;εq为偏应变,εpq=εq−q/E;Rf为土体破坏比。
本文还计算了Tresca屈服的情况,以讨论不同土体屈服准则对锥尖阻力的影响,Tresca屈服准则采用Mohr-Coulomb塑性模型实现。
1.2 有限元数值建模
CPT贯入属于空间轴对称问题,以探杆中心线为对称轴建模,如图 1。以直径36 mm,锥顶角60°的标准锥为计算对象,土体计算区域宽0.9 m,深1.4 m,经计算该计算区域足够大,边界条件对贯入阻力结果没有影响。初始时刻在土体表面预留一锥角,使锥尖与土体预先接触,可避免初始贯入时单元过度扭曲。
土单元采用CAX4R,轴对称模型共12840个单元。探杆和探头视为离散刚体,在贯入过程中与土体产生相对滑动,为此在有限元中设置接触性质,采用Coulomb摩擦接触描述界面的接触性质。土域上表面为应力边界条件,底面和外侧为位移边界。探杆与土体接触始终设为光滑,锥面与土体界面粗糙度可不同。CPT探头分为3部分(sleeve,cone tip以及retaining wall),retaining wall为锥尖以下与对称轴重合的部分,将其与锥头、探杆设置为一体,有助于减少贯入过程中锥尖处土体单元发生过度扭曲。图 1中红色刚体代表CPT探头。
采用Abaqus/Explcit求解器,取部分土域做ALE自适应网格重划分。ALE网格重划分强度设置对这类问题有显著影响,每个增量步网格重划分一次;每次重划分扫掠3次。材料密度设置取一较大值1000可缩短计算时间,同时一定程度上减小数值振荡。另外,土体泊松比对贯入阻力有一定的影响,泊松比越大数值振荡现象越明显,且随着土体刚度增大更加显著,文中泊松比统一取v=0.49。
初始时刻设置均匀分布的初始应力场,围压50 kPa。贯入速率2 cm/s,累计贯入0.72 m。前处理中建立历史变量,后处理中输出其接触力CFT2,表示锥面上总接触力的竖向分量,即得到qc。
1.3 锥尖阻力分析
提取达到稳定时的锥尖阻力值,根据式(1)得到锥尖系数,如图 2所示,半对数坐标系中,锥尖系数随刚度系数具有良好的线性关系。同时列出了本文理想弹塑性土体的结果以作对比,结果显示,相同初始刚度系数E/su情况下,理想弹塑性模型的锥尖系数高于双曲线硬化模型的锥尖系数,说明CPT锥尖极限阻力不仅与土体刚度有关,还与土体应力应变曲线有关联,这是以往的数值研究中鲜有提及的。观图中差异可知,在CPT贯入问题中考虑土体应力应变的非线性是有必要的。
图 3对比了Mises屈服准则和Tresca屈服准则饱和黏土中的锥尖系数,尽管在三轴压缩路径下土体的应力应变响应一致,但确实导致了锥尖系数Nk的差异。基于Mises准则的计算结果高于基于Tresca准则的计算结果,相差约0.3~0.6。整体上两者的差异与土体刚度系数关系不大,平均相差0.45。此外,基于理想弹塑性模型计算也得到了一致的结论。
图 4绘出了锥面摩擦系数αc对锥尖系数Nk的影响。锥面和土体间的接触通过摩擦系数αc和极限剪应力τmax设置,锥面极限剪应力设为{\alpha_{\mathrm{c}}}{\tau _{\max }},且{\tau _{\max }} = 2 s_{\mathrm{u}} / \sqrt{3} 。根据图 4中结果,摩擦系数{\alpha_{\mathrm{c}}}对锥尖系数Nk的影响约为\sqrt 3 {\alpha_{\mathrm{c}}},该值与以往的基于理想弹塑性土体的计算结果基本一致,表明锥面摩擦系数{\alpha_{\mathrm{c}}}与土体应力应变没有直接联系。
天然土层应力状态呈各向异性,采用Te等[4]定义的原位应力比\varDelta ,\varDelta = \left(\sigma_{\mathrm{v} 0}-\sigma_{\mathrm{h} 0}\right) / 2 s_{\mathrm{u}}表征天然土层中的土体应力的初始各向异性对锥尖系数的影响,通常\varDelta \in (-0.5,0.5)。图 5中随着原位应力比\varDelta 的增大,由竖向总应力\sigma_{\mathrm{v} 0}表示的锥尖系数Nk明显减小,与原位应力比呈线性关系。Yu等[2]、Lu等[5]和Liyanapathirana[6]针对理想弹塑性黏土的计算结果分别给出了1.83\varDelta ,1.9\varDelta 和1.7\varDelta 的比例关系,较为接近。本文同样针对理想弹塑性土体得到了相近的比例关系,约为1.9\varDelta 。而采用双曲线硬化模型,计算的比例系数约为1.0\varDelta ,显著小于理想弹塑性模型的1.9\varDelta 。
ALE贯入计算结果表明:土体破坏比Rf对极限锥尖阻力也是有影响的。图 6显示,随着破坏比Rf减小,锥尖系数逐渐增大,不同土体刚度时Rf的影响规律大致相同,即刚度系数E/su与破坏比Rf不存在相互影响。
2. 基于球孔扩张的锥尖阻力分析方法
观察贯入后锥尖附近土体的主应力分布情况发现,主应力方向由锥肩处水平向逐渐过渡为锥尖处的竖向。本文假设锥尖处土体在贯入作用下呈椭圆形扩张,椭圆周上径向应力{\sigma _r}等于球孔扩张极限扩孔压力plim,见图 7。
由于锥尖系数与应力水平无关,对于完全光滑的标准圆锥,假设存在椭圆锥面,其上径向应力{\sigma _r}指向中心O点,边界A点和B点处的径向应力分别为{\sigma _r}和\sigma_r^B,在各向同性土体中\sigma_r^A=\sigma_r^B,且等于球孔扩张极限扩孔压力p_{\lim }。计算半径r满足椭圆方程:
\frac{(r \cos \theta)^2}{\left(r_0\right)^2}+\frac{(r \sin \theta)^2}{\left(\sqrt{3} r_0\right)^2}=1 \text { 。 } (3) 式中:{r_0}=D/2,D为圆锥直径;夹角\theta \in (0,{\mathtt{π}}/2)。
扩张面A—B上任意一点的竖向应力分量为
\sigma_y=\sigma_r \cdot \sin \theta \ 。 (4) 沿路径A—B作二重积分得到锥头净贯入阻力q_{\mathrm{t}}(q_{\mathrm{t}}=q_{\mathrm{c}}-\sigma_{\mathrm{v} 0})为
q_{\mathrm{t}}=\frac{\int_0^{\frac{\mathtt{π}}{2}} \sigma_y}{l_{O A}}=\int_0^{\frac{\mathtt{π}}{2}} \sigma_r \cdot \sin \theta \sqrt{\frac{3}{2 \cos ^2 \theta+1} \mathrm{~d} \theta}。 (5) 根据式(5),得到q_{\mathrm{t}}=1.4038{p_{\lim }},且对于各向同性土体,有\sigma_{\mathrm{v} 0}={p_0},故锥尖系数进一步表示为
N_{\mathrm{k}}^*=1.4038 N_{\mathrm{s}}^* \text { 。 } (6) 式中:N_{\mathrm{k}}^*,N_{\mathrm{s}}^*分别为理想弹塑性土体中的锥尖系数和球孔扩张承载系数。根据文献[7]有如下表达式:
N_{\mathrm{s}}^*=\frac{4}{3}\left[1+\ln \left(G / s_{\mathrm{u}}\right)\right]。 (7) 图 8通过对比验证式(6)的准确性。可见该方法计算的锥尖系数与基于Tresca屈服的ALE计算值比较吻合,略低于基于Mises屈服的ALE计算值。同时列出了现有的几种基于小孔扩张理论的锥尖阻力分析方法[4, 8-9],这些方法研究均基于理想弹塑性土体,与ALE结果对照,可以看出,本文方法具有更高的精确度。
注意,上述推导的球孔与锥尖系数换算关系仅针对理想弹塑性土体,因为在该土体中,锥尖破坏区内土体为完全塑性,破坏区外土体为纯弹性,此时正好可以用极限平衡方法分析。而对于非理想弹塑性土体,破坏区外土体并非纯弹性,理论上极限平衡方法是不太适合的。所以,对于本文双曲线硬化模型,采用如下修正的关系式:
N_{\mathrm{k}}+\eta_{\mathrm{k}}\left(N_{\mathrm{s}}^*-N_{\mathrm{s}}\right)=1.4038 N_{\mathrm{s}}^*。 (8) 式中:\eta_{\mathrm{k}}为计算得到的锥尖系数之差N_{\mathrm{k}}^*-{N_{\mathrm{k}}}与球孔扩张承载系数之差N_{\text{s}}^*-{N_{\mathrm{s}}}的比值,根据计算结果,系数{\eta_{\mathrm{k}}}随刚度系数E/su的波动较小,可取{\eta_{\mathrm{k}}}平均值为0.67。双曲线硬化模型时球孔扩张承载系数{N_{\mathrm{s}}}的计算涉及到刚度系数E/su和破坏比Rf,由于该问题参数少且规律性强,可以直接根据有限元结果将其承载系数归纳为
N_{\mathrm{s}}=\frac{4}{3}\left[\ln \left(E / s_{\mathrm{u}}\right)-0.65 R_{\mathrm{r}}^4-0.42\right] \text { 。 } (9) 图 9验证了式(9)的合理性,图中本文方法计算的锥尖系数与基于Tresca屈服的ALE计算值吻合良好,同样地,略低于基于Mises屈服的ALE计算值。式(8)的优势之处在于,对于不同的土体模型,其与理想弹塑性模型的球孔扩张承载系数之差很容易计算,从而快速计算出锥尖系数Nk。最后,在这基础上再考虑锥面粗糙度、原位应力比等实际影响因素即可分析现场工程问题。
3. 现场试验验证
采用Nash等[10]在苏格兰Bothkennar黏土中开展的原位试验进行验证。现场采用多种测试手段对不排水剪切强度进行了测定,如图 10(a),其中原位十字板剪切所测强度沿深度具有很好的一致性,图中直接绘出了其均值线,并取原位十字板剪切试验结果作为不排水剪切强度的代表值。根据实测锥尖阻力{q_{\mathrm{c}}}、竖向总应力{\sigma_{\mathrm{v} 0}}和不排水剪切强度{s_{\mathrm{u}}}可得到实际的锥尖系数Nk,如图 10(d)。实测锥尖系数为11~12.5。
采用本文方法确定锥尖系数Nk需确定土体刚度系数,此处不排水强度已知,而剪切模量G采用由SBPM试验测得的水平模量,如图 10(c)。破坏比Rf取0.9即可。根据Lu等[5],锥尖摩擦系数一般为0.2~0.6,本文取αc=0.5。本场址初始应力均匀性较好,原位应力比的影响可忽略。最终本文计算的锥尖系数Nk与实测结果对比如图 10(d),可见本文方法较为准确。
4. 结论
当前静力触探工程中并未考虑土体非线性应力应变关系的影响,将黏土一贯视作理想弹塑性以简化参数和推广应用。针对该问题,在ABAQUS中进行了CPT模拟,发现土体刚度系数、锥面粗糙度、屈服准则、土体原位应力比和破坏比等因素均对锥尖阻力的评估有影响,进而影响黏土原位不排水剪切强度的确定。接着提出椭圆形锥尖承载机构,基于球孔扩张极限扩孔压力推导了锥尖系数的表达式。该方法用于计算锥尖系数,能够考虑土体应力应变非线性的影响。相同初始刚度情况下,理想弹塑性模型假设确实高估了锥尖系数,进而偏保守的测定出黏土的原位不排水剪切强度。建议在静力触探与土体参数的关系研究中以及工程应用中更多地考虑土体非线性特性的影响。
-
-
[1] YU H S, MITCHELL J K. Analysis of cone resistance: review of methods[J]. Journal of Geotechnical and Geoenvironmental Engineering, 1998, 124(2): 140-149. doi: 10.1061/(ASCE)1090-0241(1998)124:2(140)
[2] YU H S, HERRMANN L R, BOULANGER R W. Analysis of steady cone penetration in clay[J]. Journal of Geotechnical and Geoenvironmental Engineering, 2000, 126(7): 594-605. doi: 10.1061/(ASCE)1090-0241(2000)126:7(594)
[3] MAHMOODZADEH H, RANDOLPH M F. Penetrometer testing: effect of partial consolidation on subsequent dissipation response[J]. Journal of Geotechnical and Geoenvironmental Engineering, 2014, 140(6): 04014022. doi: 10.1061/(ASCE)GT.1943-5606.0001114
[4] TEH C I, HOULSBY G T. An analytical study of the cone penetration test in clay[J]. Géotechnique, 1991, 41(1): 17-34. doi: 10.1680/geot.1991.41.1.17
[5] LU Q, RANDOLPH M F, HU Y, et al. A numerical study of cone penetration in clay[J]. Géotechnique, 2004, 54(4): 257-267. doi: 10.1680/geot.2004.54.4.257
[6] LIYANAPATHIRANA D S. Arbitrary Lagrangian Eulerian based finite element analysis of cone penetration in soft clay[J]. Computers and Geotechnics, 2009, 36(5): 851-860. doi: 10.1016/j.compgeo.2009.01.006
[7] VESIĆ A S. Expansion of cavities in infinite soil mass[J]. Journal of the Soil Mechanics and Foundations Division, 1972, 98(3): 265-290. doi: 10.1061/JSFEAQ.0001740
[8] VESIC A S. Design of pile foundations[R]. Washington D C: National Cooperative Highway Research Program, Transportation Research Board, 1977.
[9] LADANYI B, JOHNSTON G H. Behavior of circular footings and plate anchors embedded in permafrost[J]. Canadian Geotechnical Journal, 1974, 11(4): 531-553. doi: 10.1139/t74-057
[10] NASH D F T, POWELL J J M, LLOYD I M. Initial investigations of the soft clay test site at Bothkennar[J]. Géotechnique, 1992, 42(2): 163-181. doi: 10.1680/geot.1992.42.2.163