• 全国中文核心期刊
  • 中国科技核心期刊
  • 美国工程索引(EI)收录期刊
  • Scopus数据库收录期刊

对数动态骨架本构模型的多维实现及其在ABAQUS中的应用

董青, 陈苏, 李小军, 董云, 陈亚东, 周正华, 朱俊

董青, 陈苏, 李小军, 董云, 陈亚东, 周正华, 朱俊. 对数动态骨架本构模型的多维实现及其在ABAQUS中的应用[J]. 岩土工程学报, 2025, 47(1): 192-199. DOI: 10.11779/CJGE20231021
引用本文: 董青, 陈苏, 李小军, 董云, 陈亚东, 周正华, 朱俊. 对数动态骨架本构模型的多维实现及其在ABAQUS中的应用[J]. 岩土工程学报, 2025, 47(1): 192-199. DOI: 10.11779/CJGE20231021
DONG Qing, CHEN Su, LI Xiaojun, DONG Yun, CHEN Yadong, ZHOU Zhenghua, ZHU Jun. Multi-dimensional implementation of logarithmic dynamic skeleton constitutive model and its application in ABAQUS[J]. Chinese Journal of Geotechnical Engineering, 2025, 47(1): 192-199. DOI: 10.11779/CJGE20231021
Citation: DONG Qing, CHEN Su, LI Xiaojun, DONG Yun, CHEN Yadong, ZHOU Zhenghua, ZHU Jun. Multi-dimensional implementation of logarithmic dynamic skeleton constitutive model and its application in ABAQUS[J]. Chinese Journal of Geotechnical Engineering, 2025, 47(1): 192-199. DOI: 10.11779/CJGE20231021

对数动态骨架本构模型的多维实现及其在ABAQUS中的应用  English Version

基金项目: 

国家自然科学基金青年基金项目 52408366

国家自然科学基金重大项目 52192675

详细信息
    作者简介:

    董青(1992—),女,助理研究员,主要从事岩土力学方面研究工作。E-mail:2458810997@qq.com

    通讯作者:

    陈亚东, E-mail: 354227312@qq.com

  • 中图分类号: TU43

Multi-dimensional implementation of logarithmic dynamic skeleton constitutive model and its application in ABAQUS

  • 摘要: 对数动态骨架本构以对数函数作为骨架曲线,引入“修正动骨架曲线”和“阻尼比退化系数”的概念,考虑试验阻尼对滞回曲线的影响,实现试验阻尼和滞回阻尼在数学上自洽的非线性动力本构。该本构只适用于一维土层场地时域非线性地震反应分析。基于对数动态骨架本构函数关系式,给出对数动态骨架参数的求解方法,并推导出加卸载曲线的时变切线剪切模量。通过ABAQUS软件的操作平台和等效剪应变算法,开发了基于对数动态骨架本构模型的显式子程序模块,该模型适用于二维和三维土层场地非线性地震反应分析。选取含软弱淤泥质土层场地作为计算模型,将从基岩输入不同加速度峰值的EI centro地震动的数值模拟结果进行对比,验证可考虑土体阻尼效应的对数动态骨架本构以及开发子程序的合理性和可用性。
    Abstract: The logarithmic dynamic skeleton constitutive model uses the logarithmic function as the skeleton curve and introduces the concepts of "modified dynamic skeleton curve" and "damping ratio degradation coefficient". Considering the influences of the test damping on the hysteresis curve, the nonlinear dynamic constitutive model for the test damping and hysteresis damping is realized. This constitutive method is only suitable for the nonlinear seismic response analysis of one-dimensional soil layers. In this study, based on the logarithmic dynamic skeleton constitutive function relation, the method of solving the logarithmic dynamic skeleton parameters is given, and the time-varying tangent shear modulus of the loading and unloading curve is derived. Based on the operating platform of ABAQUS software and the equivalent shear strain algorithm, a display subroutine module based on the logarithmic dynamic skeleton constitutive model is developed, which is suitable for nonlinear seismic response analysis of two-dimensional and three-dimensional soil layers. The numerical simulation results of the input EI centro ground motion with different acceleration peaks from bedrock are compared to verify the logarithmic dynamic skeleton constitutive model considering damping effects of soils and the rationality and availability of the developed subprogram.
  • 自1960年Henri Vidal[1]首次提出了“加筋土”的概念和设计理论以来,土工合成材料广泛应用于岩土工程领域[1]。在边坡中铺设土工合成材料能够优化土体应力分布,提高边坡的稳定性和承载力[2-4]。有别于二维土工合成材料如土工格栅、土工布。张孟喜[5]首先提出H-V立体筋的概念。作为一种新型加筋方式,H-V立体筋在筋–土相互作用上展现出更好的优势,除了具有水平筋作用外,竖筋能约束土体的变形并在竖筋间形成“加强土楔体”或“加强土柱”,增强土体的抗拉强度和抗剪强度[6],具有较好的应用前景。

    常用的加筋边坡稳定性分析方法主要有极限平衡法、有限元法等。改进最小势能法作为边坡稳定性分析方法之一,研究成果表明,相对极限平衡法,该方法计算方便快捷,具有不需要条块划分、不需要反复迭代求解等优点[7-9]。目前,改进最小势能法主要应用于有加固措施如锚杆、锚索、抗滑桩的边坡稳定性分析中。邹治来[10]考虑土工格栅与土体接触面之间摩擦势能和土工格栅弹性拉伸势能,建立土工格栅加筋边坡的最小势能稳定性分析方法。但是分析过程中,由于土工格栅与土体之间摩擦阻力计算过程复杂,参数众多,需要进行拉拔试验才能确定参数值,实现有很大难度。采用最小势能法分析了普通土工格栅加筋边坡稳定性,李铀等[11]、孙加平等[12]提出了考虑滑床剪切势能影响的近似计算模型,改进最小势能对锚杆(索)加固边坡稳定分析方法。温树杰等[13]分析桩体变形特征,考虑抗滑桩与滑体相互接触处的弹性压缩势能、挠曲势能提出基于最小势能法加固边坡稳定性分析方法。未检索到最小势能法应用于H-V立体筋加筋边坡稳定性分析。

    H-V立体筋边坡稳定性分析主要采用条分法,张孟喜等[14]采用瑞典条分法对立体加筋边坡的稳定性进行探讨,在计算过程中未考虑条间力的作用,容易导致计算结果偏小。马学宁等[15]建立水平条间力的假设条件,推导出水平条分瑞典法和简布法分析H-V立体加筋边坡的稳定性,该方法需要在水平筋材处于土条的中线上的条件下进行,且水平条分简布法需要进行反复迭代,计算量较大。其它对于H-V立体筋边坡稳定性分析文献尚不多见。

    因此,本文基于改进最小势能法,对H-V立体筋作用的基本原理和受力特性进行分析,构建出水平筋和竖筋的弹性势能方程,并对加筋边坡单元体虚位移方向的剪切势能进行推导,通过计算边坡最危险滑裂面系统势能最小时的虚位移,求解出边坡安全系数,并通过算例对该方法的适用性和可行性进行验证。

    本文利用改进最小势能法从能量的角度对H-V立体筋边坡的稳定性展开分析,在系统中涉及的能量种类主要有滑体储存的弹性势能,H-V筋因水平筋线弹性变形产生弹性势能、筋土之间因摩擦储存摩擦势能、竖筋产生侧向变形储存弹性势能和滑面剪切势能,概念清晰,可直接计算加筋边坡稳定性分析结果,相对于计算复杂的有限元法和需反复迭代、需对超静定问题进行多次假定的极限平衡法更有实用价值,为H-V立体筋边坡稳定性分析提供新的思路。

    图1所示的H-V立体筋均质边坡,设潜在滑裂面为AC,其曲线方程表示为y=f(x),土体黏聚力为c,内摩擦角为φ,重度为γ,土体的重力为W,第i排筋条与水平方向的夹角为αi,水平方向布置则αi=0;第i排筋条与竖直方向的夹角为βi,每排筋条的单位方向向量记为T=(cosαi,cosβi)。为合理建立力学模型,做如下说明:①假设滑动体整体为刚性体,在滑裂面上存在弹性压缩变形和剪切变形,储存的弹性压缩势能为Ve,剪切势能为Vτ;②在合外力R的作用下加筋边坡与滑体产生相同的虚位移d=(d1,d2);③微分段的长度为dl,弹性压缩变形可用刚度系数为Ki的弹簧模拟,其大小与微段dl呈正比,即Ki=midl,mi为地基系数,可通过试验测定或查阅相关资料获得。

    图  1  H-V加筋均质边坡计算图示
    Figure  1.  Calculation of H-V ribs-reinforced homogeneous slope

    滑裂面方程为y=f(x),则求得曲线上任意切线上法线的斜率为k=1/f(x)。滑裂面外法线单位向量和滑面切向单位向量可表示为

    n(x)=[n1(x),n2(x)]={f(x)1+[f(x)]2,11+[f(x)]2}, (1)
    t(x)=[t1(x),t2(x)]={11+[f(x)]2,f(x)1+[f(x)]2} (2)

    假设在各种力的作用下,滑体发生虚位移,使系统势能最小,滑体此时的势能即是系统达到极限平衡状态时的最小势能,则任意微段储存的弹性势能为

    Ve1=i=1N12Ki(dn)2=12lm(dn)2dl (3)

    H-V立体筋结构在传统水平筋的基础上设置竖筋,竖筋宽度可与水平筋相同,其截面形状可根据工程需要采用不同厚度变化的矩形、多边形等,形成一种空间形式的三维组合加筋结构,如图2所示。

    图  2  H-V立体筋三维结构示意图
    Figure  2.  Three-dimensional structural diagram of H-V ribs

    H-V立体筋结构的基本原理:将一定量具有抗拉强度的H-V立体筋铺设在土体中,水平筋的拉力作用会逐渐延伸到路堤边坡自由面的筋条中,水平筋土之间摩擦力提供锚固力;竖筋的作用则表现为增加筋材拔出过程中竖筋的抗阻力,使土体结合更紧密,将整个H-V筋锚固在了土体中,提高边坡承载力。相对普通均质边坡,H-V筋的这种拉力状态也会使上部土体传递来的荷载较为集中,相对更均匀地分布到筋材下部的土体,使下部土体受力更合理。

    构建H-V立体筋势能函数时需要确定水平筋与滑面交点坐标和竖筋有效数目,因此首先需要确定筋条的特征坐标。加筋边坡在发生失稳破坏时,根据坡体的滑动状态可分为滑动区和锚固区,如图3所示,设坡面倾角为θ,筋材各层的垂直间距为hi,第一层水平筋距离坡脚的竖直距离为h1。则第i层水平筋与滑动面的交点纵坐标yi可表示为

    图  3  筋条的特征坐标计算简图
    Figure  3.  Simplified calculation of characteristic coordinates of geogrid
    yi=y1+hi(i1) (4)

    将式(4)与滑动面方程y=f(x)联立即可确定各排筋条与滑面的交点坐标,若方程组有解,则筋条产生变形储存了势能,确认筋条储存势能且取x>0时,第i层水平筋与滑动面的交点横坐标xi可表示为

    xi=f(yi)          (x>0) (5)

    土工格栅长度为l,则第i层水平筋右末端横坐标xi表示为

    xi=x1+l+hi(i1)/tanθ (6)

    由式(5),(6)可得第i层土工格栅的锚固区长度为

    li=xixi (7)

    设竖筋间距为D,则锚固区内第i层内竖筋的个数ni可表示为

    ni=li/D+1 (8)

    则可计算出锚固区内竖筋的总个数n=i=1mni

    在土体滑动的过程中,水平筋与土体之间摩擦力作功会产生摩擦势能,水平筋自身也会由于拉伸变形产生线弹性应变能。

    对水平筋进行作用力分析时,由于筋土之间接触复杂,水平筋和填土分开考虑计算筋材的摩擦力参数较多,实现起来有很大的难度。加筋路堤稳定性计算时,土工合成材料应用手册将水平筋与路堤填土分开考虑,所以在计算稳定性系数时水平筋的作用力一般直接取其抗拉强度值。本文将水平筋与其附近土体考虑为一种复合体,所以筋条的抗拉力应该分散到其作用范围以内的土体中[14]。假设筋条的抗拉力平均分配到每层筋条的土体中,若路堤中每层筋条水平间距为Sx,竖向间距为Sy图4)。

    图  4  水平筋复合抗拉力计算示意图
    Figure  4.  Schematic diagram of calculation of horizontal ribs

    则第i个土条中单位厚度复合土体中抗拉力为

    Tst=ThSxSylisinθi, (9)

    式中,Th为单根水平筋条抗拉力,Sx为筋条水平间距,Sy为筋条竖向间距,li为第i条底弧长度,θi为第i条底弧仰角。

    图5所示,对单位长度为li的土体单元进行分析,水平筋的摩擦阻力Tki与水平筋的作用力Tst平衡,所以满足Tki=Tst,由于土工格栅所受摩擦阻力与位移呈钝角关系,所以摩擦所储存的势能为力所做功的负值。

    图  5  土体单元水平筋受力分析图
    Figure  5.  Analysis diagram of horizontal reinforcement of soil element
    Ve2=i=1nTki(dT), (10)

    式中,Tki为第i排水平筋摩擦阻力,d为滑体虚位移,T为水平筋水平轴向单位向量。

    根据以往文献资料结果显示,基于筋-土界面加筋摩擦机理的圆弧滑动分析方法可以计算加筋后的安全系数,在抗滑项中加入界面摩擦阻力,但结果可能偏小[16]。本文将水平筋与其附近土体考虑为一种复合体,因为筋土间的界面被忽略,在计算水平筋所受摩擦阻力时缺少界面摩擦阻力,抗滑力中缺少界面摩擦项,且土工格栅摩擦阻力是在抗拉强度条件下计算的,因此可能导致安全系数计算结果偏大,具体影响还需要进一步研究。

    分析水平筋的线应变势能时,滑体在向最小势能方向移动时会产生虚位移,因筋材均必须要伸入到较深坡体内部,因此可假设筋材与右端与滑裂面的交点g处是固定的[11]。设布置了m排筋材,筋条的刚度为K,如图6所示,在合外力作用下滑体发生虚位移d,水平筋的坡面端头由e移至e,伸长量为Δk,假设e坐标为(xe, ye),g的坐标为(xg, yg)。可得

    图  6  立体筋材结构受力示意图
    Figure  6.  Schematic stress diagram of three-dimensional rib structure
    Δk=|gege|=(d1+vkx)2+(d2+vky)2vkx2+vky2, (11)

    式中,vkx=xgxe,  vky=ygye

    水平筋的线应变能为

    Ve3=i=1m[12KΔk2]=K2i=1m[(d1+vkx)2+(d2+vky)2+(vkx2+vky2)]Ki=1m[(d1+vkx)2+(d2+vky)2vkx2+vky2] (12)

    当H-V立体筋在荷载作用下竖筋首先发挥作用,在其自身刚度下对土体产生侧阻力,阻挡土体向两侧运动,这相当于给土体增加多个侧向约束,竖筋对土体的势能主要是侧阻力对竖筋作用产生的变形势能。

    对竖筋的作用力进行分析,竖筋的侧向抵抗作用在竖筋间将形成“土体加强区”,以有效改变加筋土的受力状况[6]。如图7所示,在竖筋两侧分别作用有主动土压力σa和被动土压力σp,沿竖筋表面作用有剪切应力τv,τh,在竖筋顶部作用有竖向压力σb。由于竖筋表面积与厚度及其顶部面积相对路堤尺寸可忽略,因此竖筋的主要作用力为竖向筋条两侧的主动土压力及被动土压力,在稳定性分析时,忽略竖筋表面的剪应力作用,竖筋的侧阻力可由被动土压力与主动土压力之差计算得到[14]

    图  7  竖筋单元受力图
    Figure  7.  Force diagram of three-dimensional tendons

    则主动、被动土压力计算公式为

    Ta=(γHiKa2cKa)A, (13)
    TP=(γHiKP+2cKP)A (14)

    竖筋作用力为Tsh,

    Tsh=TpTa=[TpTa]=[γHi(KpKa)+2c(Kp+Ka)]A, (15)

    式中,c为土体黏聚力,γ为土体重度,Hi为竖筋上覆土高度,φ为内摩擦角,Ka为主动土压力系数,Kp为被动土压力系数,且Ka=tan2(45°φ/2),Ka=tan2(45°+φ/2),A为竖筋侧面积,A=Bh,其中B为竖筋宽度,h为竖筋高度。

    由于竖筋的侧向变形对H-V筋中竖筋作用发挥起重要作用,所以尽管竖筋高度相对边坡尺寸较小,也需要进行势能分析。在进行H-V加筋边坡稳定性分析时,考虑到竖筋的排列会对微分圆弧作用面的受力产生影响,所以在计算抗滑力时只计算锚固区(即破裂面外部)的有效竖筋阻力,锚固区内竖筋总个数为n,则竖筋在侧阻力作用下产生的变形势能为

    Ve4=i=1nTshΔx=i=1nTshd1, (16)

    式中,Tsh为竖筋作用力,Δx为竖筋水平方向的变形量,大小为虚位移d在水平方向的分量。

    通过式(15)可看出竖筋作用力Tsh的作用点位于竖筋高度的中部,但由于竖筋高度h较小,在计算立体筋合力时,竖筋作用力的方向可近似认为作用在水平筋上。则立体筋的合力为竖筋作用力Tsh和水平筋作用力Tst的和,即

    Ti=Tst+nTsh (17)

    边坡的失稳往往是沿着土体抵抗力最小的方向滑动,即为使坡体系统势能最小的虚位移方向。由文献[8]可知采用弹性势能以及合外力做功的方式表示边坡的剪切势能,可大大简化计算量,同时便于程序的实现。为便于势能公式的表达,令Δ=R12+R22,Δ=d12+d22,则合外力R和虚位移d的方向向量可表示为

    Rr=1Δ(R1,R2) ,dd=1Δ(d1,d2} (18)

    图1中任取土条作为研究对象,土条的受力分析如图8所示。

    图  8  边坡土条受力图
    Figure  8.  Force diagram of slope strips

    作用在土条上的力有法向力Ni,剪应力τi,合外力为Ri,H-V立体筋合力为Ti,土条左右面上的条间力为Ei,Ei+1。由于边坡失稳沿着岩土体抵抗力最小的方向,即虚位移d,因此各个土条左右面上的条间力作用方向与虚位移方向一致,在计算过程中,将第i个土条条间力合力记为Ei,Ei+1,将第i+1个土条条间力合力记为Ei+1,Ei+2,ΔE=Ei+1Ei,则各个土条沿虚位移方向满足平衡方程:

    Nicosα+τilicosβ+Ticosω+ΔEi=Ricosϕ (19)

    式中 α为任意微面上法向力与虚位移方向的夹角,且cosα=n(dd)li为滑裂面微段长度;β为剪应力与虚位移的夹角,且cosβ=t(dd)ω为水平筋与虚位移的夹角,cosω=T(dd)ϕ为合外力与虚位移的夹角,cosϕ=Rrdd

    由式(19)可知滑面上剪应力的解析解为

    τi=RiΔRidd+TiTddNiniddΔEiddtli (20)

    为得到线弹性范围内剪切应变势能产生的函数,需要先构建剪应变和剪应力的力学求解模型,土体由于剪切变形发生的剪应变为δi,在剪切力τi的作用下滑体沿着滑裂面切向会产生相对位移,即由滑裂面上的点f移动至点f,受剪切变形影响的土体深度为hi,剪切位移记为d,则

    δi=dhi (21)

    滑面上的剪切势能为

    Vτ=12i=1nτiδihi=12i=1nRiRidd+TiTddNiηiddΔEiddtlidthihili=12i=1n(RiRd+TiTdNiηidΔEi)=Ve112Rd12i=1nTiTd+12i=1nΔEi, (22)

    式中,i=1nΔEi的物理意义为各个条块条间力的矢量和,由于各个条块之间的作用力属于作用力与反作用力,因此i=1nΔEi=0,则上式可化简为

    Vτ=Ve112Rd12i=1nTiTd (23)

    综合式(3)、(10)、(12)、(17)、(23)可得H-V立体筋边坡系统的总势能为

    V=Ve1+Ve2+Ve3+Ve4+VτRd=2Ve1+Ve2+Ve3+Ve432Rd12i=1nTiTd (24)

    安全系数Fs为沿着虚位移d方向投影得到的抗滑力Fm与下滑力Fn的比值[8],即

    Fs=FmFn , (25)

    式中,抗滑力主要由微分圆弧dl的极限抗滑力和H-V立体筋复合抗拉力产生的抗滑力提供,下滑力由合外力(包括滑体自重)提供。本文的抗滑力Fm及下滑力Fn均与虚位移相关,因此需先求解虚位移,再进行安全系数求解。

    根据最小势能原理可知,H-V立体筋边坡系统势能最小时,总势能函数在虚位移的一阶偏导处需满足

    Vd1=0 ,Vd2=0 } (26)

    将式(24)代入式(26),即可求出系统最小势能时坡体的虚位移d。其中令

    a1=i=1nkiηx2=mx1x2[f(x)]21+[f(x)]2dx, (27)
    a2=a3=i=1nkiηx2=mx1x2f(x)1+[f(x)]2dx, (28)
    a4=i=1nkiηy2=mx1x211+[f(x)]2dx, (29)
    A=K[d1+vkx2(vkx2+vky2)(d1+vkx)[(d1+vkx)2+(d2+vky)2](vkx2+vky2)], (30)
    B=K[d2+vky2(vkx2+vky2)(d2+vkx)[(d1+vkx)2+(d2+vky)2](vkx2+vky2)], (31)
    C=(12i=1nTi+i=1nTkii=1nTsh)cosαi+32Rx, (32)
    D=(12i=1nTi+i=1nTki)cosβi+32Ry (33)

    则式(26)可表示为

    2(a1d1+a2d2)+A=C ,2(a4d2+a3d1)+B=D } (34)

    根据方程组计算得出d1,d2

    抗滑力Fm由滑裂面上各微段的法向力Ni、极限剪切力τili以及水平加筋等效应力或H-V立体加筋应力沿着虚位移方向提供的抗滑力投影之和组成。本文研究的土工格栅对抗滑力的贡献主要体现在增加了滑体的剪切力。整个滑裂面上任意微段的法向力Ni提供的下滑力为

    Fm1=i=1nNicosα=i=1nm(dn)licosα=i=1nm[d1f(x)d2]ddnli, (35)

    式中,Tsh为任意微面上法向力与虚位移方向的夹角,且cosα=n(dd)

    由莫尔-库仑破坏准则,微分圆弧dl滑动面上极限抗剪力在虚位移方向上所能提供的抗滑力Fm2

    Fm2=x1x2dFm2=x1x2[c+m(dn)tanφ]dlcosβ, (36)

    式中,β为剪应力与虚位移的夹角,cosβ=t(dd)

    任意微段上的H-V立体筋复合抗拉力提供的抗滑力Fm3

    Fm3=i=1nTicosω=i=1n(Tst+Tsh)cosω, (37)

    式中,ω为水平筋与虚位移的夹角,cosω=T(dd)

    在虚位移方向上总的抗滑力Fm

    Fm=Fm1+Fm2+Fm3 (38)

    作用在滑体上的所有力中,只有合外力R沿着虚位移c=10方向上提供下滑力Fn,则总下滑力为

    Fn=i=1nRicosϕ=|R|R1d1+R2d2ΔΔ (39)

    图9所示的H-V立体筋均质边坡[15],边坡坡度为1∶1,边坡高度6 m,土体黏聚力c=10kPa,重度γ=18kN/m3,内摩擦角φ=20。该边坡采用的H-V立体筋材料为单向拉伸高密度聚乙烯拉筋,筋材的的抗拉刚度为10000 kN/m,土体的地基系数为100 kN/m3,共铺设5层,每层筋材长为6 m,第一层距坡顶0.6 m,最后一层距坡底0.6 m,竖筋间距D=20 cm,筋材之间的间距为1.2 m,竖筋高度H=10 mm。

    图  9  H-V立体筋边坡示意图
    Figure  9.  Schematic diagram of reinforced slope

    建立如图9所示直角坐标系,通过4.5H法搜索滑动面圆心线,通过改变边坡坡度,得到不同的滑动面的方程,利用MATLAB编程试算求解出安全系数的数值最小的滑裂面即为最危险滑动面。则试算得出的边坡最危险滑动面圆心o的坐标为(-1.02,10.72),半径R=10.77 m,滑动面方程AC

    (x+1.02)2+(y10.72)2=10.772

    本文采用改进最小势能法计算未加筋、加水平筋及加H-V立体筋边坡的安全系数,并与文献[14,15]的水平条分瑞典法的计算结果进行比较,计算结果如表1所示。

    表  1  与文献计算结果对比
    Table  1.  Comparison between calculated results and those in published papers
    项目改进最小势能法文献[14]的方法文献[15]的水平条分瑞典法
    未加筋1.1781.2041.215
    加水平筋1.5631.5761.749
    加立体筋1.7281.7192.121
    下载: 导出CSV 
    | 显示表格

    表1可见,对未加筋边坡,本文采用的改进最小势能法计算的安全系数与文献[14,15]的计算结果相差最大值为0.037,相差最小值0.026,相对误差小于3.05%,证明本文的计算方法是合理可行的。对于加水平筋,加立体筋边坡稳定性提高了14%,说明立体筋加筋效果更明显。对比文献[15]的水平条分瑞典法计算结果,本文的计算结果明显偏小,原因在于对加筋边坡稳定性分析时考虑了单个条块水平条间合力为0和整体力矩平衡的条件,考虑了竖向条间力作用,使计算结果偏大。

    目前,文献[15]采用的水平条分法还是一种相对新的计算方法,为了更好地验证本文计算方法的可靠性,针对算例在表2中设置不同的参数与极限平衡法(Bishop法和Janbu法)在无筋(不加筋)、加筋(加水平筋和H-V立体筋)状态下的计算结果进行对比。

    表  2  算例参数
    Table  2.  Parameters of examples
    编号黏聚力c/kPa重度γ/(kN·m-3)内摩擦角φ/(°)地基系数m/(kN·m-3)
    1101820100
    2101822100
    3101824100
    4101826100
    552030200
    6102030200
    7152030200
    8202030200
    下载: 导出CSV 
    | 显示表格

    根据以往文献资料表明,当边坡为单一地层且不考虑加固措施时,通过最小势能法计算边坡稳定性系数过程中可将地基系数m约去[8],即单一均质边坡稳定性不受地基系数的影响。因为极限平衡法分析边坡稳定性参数中没有地基系数m,为避免地基系数测试结果取值误差产生干扰,以及为后文验证H-V立体筋对提高边坡稳定性的作用,设置表3进行改进最小势能法和极限平衡法无筋边坡计算结果对比。

    表  3  无筋边坡算例计算结果对比
    Table  3.  Comparison of calculated results of examples of unreinforced slopes
    对比分组表2编号改进最小势能法极限平衡法备注
    BishopJanbu
    A11.1781.1171.098φ渐增
    21.1801.1441.137
    31.1821.1791.160
    41.1841.2071.192
    B51.1171.2421.231c渐增
    61.2381.3371.329
    71.3441.4151.411
    81.5941.5881.584
    下载: 导出CSV 
    | 显示表格

    表3计算结果可知,其他参数保持不变,设置两组对比方式,A组内摩擦角φ逐渐增大,B组黏聚力c逐渐增大,从计算结果可以看出无筋边坡利用改进最小势能法与极限平衡法计算结果相对误差保持在10%范围以内,证明本文分析计算边坡稳定性的方法是合理可行的,且改进最小势能法计算结果偏大于Janbu法,偏小于Bishop法。

    表4可看出:①改进最小势能法计算结果跟极限平衡法对比,施加水平筋相对误差小于4.7%,施加立体筋相对误差小于2.2%,证明改进最小势能法与极限平衡方法的计算结果具有很好的一致性。②加筋边坡的安全系数均随土体强度指标的增加呈现递增的趋势,其他参数保持不变,对照组C内摩擦角φ逐渐增大,对照组D黏聚力c逐渐增大,分析数值变化规律可知时D组加筋边坡安全系数值变化更大,可知黏聚力c对边坡稳定性影响较大。③综合表3,4计算结果还可分析得出,对比无筋边坡,施加水平筋边坡的安全系数提高33.4%,对比施加水平筋,施加立体筋边坡的安全系数提高16.5%,说明施加立体筋在提高边坡稳定性方面的效果更显着。在不同的条件下,采用改进最小势能法的计算结果对比极限平衡法(Bishop法和Janbu法)还存在微小差异,主要因为本文的计算方法与极限平衡法采用的计算模型不一样,且对安全系数的定义和计算仍存在差异性。

    表  4  加筋边坡算例计算结果对比
    Table  4.  Comparison of calculated results of examples of reinforced slope
    对比分组表2编号改进最小势能法极限平衡法备注
    BishopJanbu
    水平筋H-V筋水平筋H-V筋水平筋H-V筋
    C11.5631.7491.5851.7351.5591.723φ渐增
    21.5651.7611.6131.7621.5581.758
    31.5671.7841.6451.7911.5611.782
    41.5701.7921.6781.8301. 5761.804
    D51.6301.9191.6461.9271.6271.930c渐增
    61.8242.0781.8442.1051.8152.111
    71.8792.1251.9022.1731.8802.123
    81.9792.2212.0202.2501.9722.218
    下载: 导出CSV 
    | 显示表格

    本文在改进最小势能边坡稳定性分析基础上,提出H-V立体筋边坡稳定性分析新方法,并通过算例验证该方法的有效性和合理性,得出以下两点结论。

    (1)本文基于改进的最小势能算法,考虑水平筋的摩擦势能、线弹性应变能以及竖筋变形势能的影响,进一步分析了加筋边坡的剪切势能,并将其表示为弹性势能的函数,简化了计算过程。通过计算边坡最危险滑裂面沿虚位移方向静力平衡方程,得到适用于施加水平筋和H-V立体筋的边坡稳定性分析方法。

    (2)利用改进最小势能法求解文献算例,将其计算结果与极限平衡法(Bishop法和Janbu法)的计算的结果相比较,相对误差在10%以内,表明本文分析方法是可行且合理的。相比施加水平筋,施加立体筋时边坡安全系数提高了14%,说明H-V立体筋在提高边坡稳定性方面效果显着。

  • 图  1   对数动骨架曲线构造规则

    Figure  1.   Construction rules of logarithmic dynamic skeleton curve

    图  2   阻尼比修正骨架曲线示意图

    Figure  2.   Diagram of damping ratio-modified skeleton curve

    图  3   动骨架参数b求解示意图

    Figure  3.   Schematic diagram of solving dynamic skeleton parameters

    图  4   VUMAT子程序三维算法实现流程图

    Figure  4.   Flow chart of VUMAT subroutine 3D algorithm implementation

    图  5   EI centro地震加速度时程

    Figure  5.   Time histories of EI centro earthquake acceleration

    图  6   正弦波形输入下两本构模拟地表加速度时程和软弱层应力应变曲线

    Figure  6.   Nonlinear seismic responses of a site with weak soil layer under input ground motion of sinusoidal acceleration

    图  7   输入0.5 m/s2 EI centro地震动

    Figure  7.   Input ground motions of 0.5m /s2 EI centro

    图  8   输入1 m/s2 EI centro地震动

    Figure  8.   Input ground motions of 1 m /s2 EI centro

    图  9   输入1.5 m/s2 EI centro地震动

    Figure  9.   Input ground motions of 1.5 m /s2 EI centro

    图  10   地表加速度归一化反应谱

    Figure  10.   Normalized response spectra of surface acceleration

    表  1   土体动力剪切非线性特性参数

    Table  1   Dynamic shear nonlinear characteristic parameters of soils

    土类 a0/10-9 b0/10-5 a1/10-3 b1
    1 52.61 9.82 1.05 5.80
    2 51.92 13.45 0.90 5.78
    3 42.68 7.82 1.20 6.06
    2 35.43 9.34 0.90 5.78
    4 26.57 4.63 0.80 5.49
    5 20.14 3.31 1.00 5.68
    6 15.64 1.80 1.10 5.52
    7 9.36 0.93 1.70 6.42
    8 9.36 1.08 0.80 5.93
    9 8.28 0.95 1.00 5.44
    10 1.00 1.72 7.18 3.70
    下载: 导出CSV

    表  2   含软弱层场地条件

    Table  2   Site conditions of soft layers

    土类 土层名称 埋深/m 平均vs/
    (m·s-1)
    密度/
    (t·m-3)
    1 淤泥 2.95 109 1.60
    2 淤泥质粉质黏土 4.45 110 1.62
    3 淤泥质黏土 8.95 119 1.68
    2 淤泥质粉质黏土 11.95 134 1.62
    4 粉质黏土 13.45 140 1.92
    5 淤泥质粉质黏土 16.50 144 1.78
    6 粉质黏土 20.50 171 1.93
    7 粉质黏土混碎石 22.60 230 2.02
    8 黏土 25.00 234 1.95
    9 粉质黏土 27.20 247 1.98
    10 计算基底 666 2.65
    下载: 导出CSV
  • [1]

    DUNCAN J M, CHANG C Y. Nonlinear analysis of stress and strain in soils[J]. Journal of the Soil Mechanics and Foundations Division, 1970, 96(5): 1629-1653. doi: 10.1061/JSFEAQ.0001458

    [2]

    CHEN G X, RUAN B, ZHAO K, et al. Nonlinear response characteristics of undersea shield tunnel subjected to strong earthquake motions[J]. Journal of Earthquake Engineering, 2020, 24(3): 351-380. doi: 10.1080/13632469.2018.1453416

    [3]

    DAFALIAS Y, HERMANN L. Bounding surface formulation of soil plasticity[C]// John Wiley & Sons, Chichester, 1982.

    [4]

    BOULANGER R W, ZIOTOPOULOU K. A constitutive model for clays and plastic silts in plane-strain earthquake engineering applications[J]. Soil Dynamics and Earthquake Engineering, 2019, 127: 105832. doi: 10.1016/j.soildyn.2019.105832

    [5]

    YAO E L, WANG S Y, RUAN B, et al. Numerical study on site response considering ground motion spatial variation[J]. Soil Dynamics and Earthquake Engineering, 2019, 127: 105836. doi: 10.1016/j.soildyn.2019.105836

    [6]

    SHAHIR H, MOHAMMADI-HAJI B, GHASSEMI A. Employing a variable permeability model in numerical simulation of saturated sand behavior under earthquake loading[J]. Computers and Geotechnics, 2014, 55: 211-223. doi: 10.1016/j.compgeo.2013.09.007

    [7]

    PREVOST J H, KEANE C M. Shear stress-strain curve generation from simple material parameters[J]. Journal of Geotechnical Engineering, 1990, 116(8): 1255-1263. doi: 10.1061/(ASCE)0733-9410(1990)116:8(1255)

    [8] 荣棉水, 卢滔, 李小军, 等. 一种基于动态骨架曲线的土层时域积分方法及其验证[J]. 应用基础与工程科学学报, 2013, 21(1): 78-89.

    RONG Mianshui, LU Tao, LI Xiaojun, et al. A direct time-domain integral method based on dynamic skeleton curve constitutive model and its verification[J]. Journal of Basic Science and Engineering, 2013, 21(1): 78-89. (in Chinese)

    [9]

    YNIESTA S, BRANDENBERG S J, SHAFIEE A. ARCS: a one dimensional nonlinear soil model for ground response analysis[J]. Soil Dynamics and Earthquake Engineering, 2017, 102: 75-85. doi: 10.1016/j.soildyn.2017.08.015

    [10]

    ZHUANG H, CHEN G. A viscous-plastic model for soft soil under cyclic loadings[J]. Geotechnical Special Publication, 2006: 343-350.

    [11] 庄海洋, 陈国兴, 朱定华. 土体动力黏塑性记忆型嵌套面本构模型及其验证[J]. 岩土工程学报, 2006, 28(10): 1267-1272.

    ZHUANG Haiyang, CHEN Guoxing, ZHU Dinghua. Dynamic visco-plastic memorial nested yield surface model of soil and its verification[J]. Chinese Journal of Geotechnical Engineering, 2006, 28(10): 1267-1272. (in Chinese)

    [12] 赵丁凤, 阮滨, 陈国兴, 等. 基于Davidenkov骨架曲线模型的修正不规则加卸载准则与等效剪应变算法及其验证[J]. 岩土工程学报, 2017, 39(5): 888-895. doi: 10.11779/CJGE201705013

    ZHAO Dingfeng, RUAN Bin, CHEN Guoxing, et al. Validation of modified irregular loading-unloading rules based on Davidenkov skeleton curve and its equivalent shear strain algorithm implemented in ABAQUS[J]. Chinese Journal of Geotechnical Engineering, 2017, 39(5): 888-895. (in Chinese) doi: 10.11779/CJGE201705013

    [13]

    RUAN B, ZHAO K, WANG S Y, et al. Numerical modeling of seismic site effects in a shallow estuarine bay (Suai Bay, Shantou, China)[J]. Engineering Geology, 2019, 260: 105233. doi: 10.1016/j.enggeo.2019.105233

    [14] 董青, 周正华, 苏杰, 等. 基于对数动骨架考虑可逆孔压的有效应力本构研究[J]. 岩土工程学报, 2020, 42(12): 2322-2329. doi: 10.11779/CJGE202012020

    DONG Qing, ZHOU Zhenghua, SU Jie, et al. Constitutive model for effective stress based on logarithmic skeleton curve considering reversible pore pressure[J]. Chinese Journal of Geotechnical Engineering, 2020, 42(12): 2322-2329. (in Chinese) doi: 10.11779/CJGE202012020

    [15]

    DONG Q, CHEN S, JIN L G, et al. Time-domain dynamic constitutive model suitable for mucky soil site seismic response[J]. Earthquake Engineering and Engineering Vibration, 2024, 23(1): 1-13. http://www.cnki.com.cn/Article/CJFDTotal-EEEV202401001.htm

    [16] 李小军, 廖振鹏, 张克绪. 考虑阻尼拟合的动态骨架曲线函数式[J]. 地震工程与工程振动, 1994, 14(1): 30-35.

    LI Xiaojun, LIAO Zhenpeng, ZHANG Kexu. A functional formula of dynamic skeleton curve taking account of damping effect[J]. Earthquake Engineering and EngineeringVibration, 1994, 14(1): 30-35. (in Chinese)

    [17]

    LIAO Z P, YANG G. Multi-directional transmitting boundaries for steady-state SH waves[J]. Earthquake Engineering & Structural Dynamics, 1995, 24(3): 361-371.

    [18]

    CHEN G X, ZHAO D F, CHEN W Y, et al. Excess pore-water pressure generation in cyclic undrained testing[J]. Journal of Geotechnical and Geoenvironmental Engineering, 2019, 145(7): 4019022. http://www.xueshufan.com/publication/2941705810

  • 期刊类型引用(4)

    1. 邢亮. 基于判别分析法和机器学习对边坡稳定性的预测研究. 广东土木与建筑. 2024(12): 25-29+50 . 百度学术
    2. 郝广杰. 考虑土层变形特征的边坡稳定性分析方法研究. 粉煤灰综合利用. 2023(02): 62-67 . 百度学术
    3. 温树杰,赖光甜,孙自立. 改进的边坡最小势能法及潜在滑动方向研究. 铁道科学与工程学报. 2022(08): 2249-2258 . 百度学术
    4. 何军. 深层岩土开采中边坡水土流失控制方法研究. 环境科学与管理. 2022(09): 148-153 . 百度学术

    其他类型引用(3)

图(10)  /  表(2)
计量
  • 文章访问数:  283
  • HTML全文浏览量:  43
  • PDF下载量:  62
  • 被引次数: 7
出版历程
  • 收稿日期:  2023-10-16
  • 网络出版日期:  2024-05-29
  • 刊出日期:  2024-12-31

目录

/

返回文章
返回