Loading [MathJax]/jax/output/SVG/jax.js

    基于几何细分的三维地质模型自适应精细化构建方法

    孙黎明, 刘禹杉, 张睿卓, 杨玉坚, 严俊

    孙黎明, 刘禹杉, 张睿卓, 杨玉坚, 严俊. 基于几何细分的三维地质模型自适应精细化构建方法[J]. 岩土工程学报, 2023, 45(S1): 244-248. DOI: 10.11779/CJGE2023S10049
    引用本文: 孙黎明, 刘禹杉, 张睿卓, 杨玉坚, 严俊. 基于几何细分的三维地质模型自适应精细化构建方法[J]. 岩土工程学报, 2023, 45(S1): 244-248. DOI: 10.11779/CJGE2023S10049
    SUN Liming, LIU Yushan, ZHANG Ruizhuo, YANG Yujian, YAN Jun. Adaptive smooth building method for 3D geological model based on geometric subdivision[J]. Chinese Journal of Geotechnical Engineering, 2023, 45(S1): 244-248. DOI: 10.11779/CJGE2023S10049
    Citation: SUN Liming, LIU Yushan, ZHANG Ruizhuo, YANG Yujian, YAN Jun. Adaptive smooth building method for 3D geological model based on geometric subdivision[J]. Chinese Journal of Geotechnical Engineering, 2023, 45(S1): 244-248. DOI: 10.11779/CJGE2023S10049

    基于几何细分的三维地质模型自适应精细化构建方法  English Version

    基金项目: 

    国家自然科学基金项目 U19A2049

    详细信息
      作者简介:

      孙黎明(1987—),男,博士,高级工程师,主要从事智慧水利方面的研究工作。E-mail: sunlm@iwhr.com

    • 中图分类号: TP332

    Adaptive smooth building method for 3D geological model based on geometric subdivision

    • 摘要: 为钻孔采样数据设计的广义三棱柱体元是地质建模中的常用方法,它利用相邻的三棱柱体元表达地层空间分布情况,而通常相邻钻孔采样间隔较大,导致这些区域模型的精度较低,地层面不光滑。基于此提出了三棱柱三维地质模型的插值平滑的理论和方法,在不该变模型数据结构的情况可以自适应地细分三棱柱模型实现高精度插值。①提出了三棱柱三维地质模型几何平滑度的定义与计算方法,用于判断一个三棱柱体元是否需要插值。②详细讨论了基于三棱柱几何平滑度的三维地质模型插值虚拟钻孔设计与计算方法。③基于虚拟钻孔设计计算和基于现有三棱柱体元的三维地质模型自适应平滑插值方法研究。
      Abstract: The generalized triangular prism (GTP) voxel designed for the borehole data is the commonly method in the geological modelling, and it can be used to build the 3D solid geological model to express the stratigraphic space through the adjacent GTP voxel based on the borehole sampling data. Its accuracy is lower, and the stratigraphic surface is non-smooth when the interval of the adjacent borehole is long, so the theory and method of interpolation and smoothing of the GTP 3D geological model is proposed, which can interpolate the GTP model adaptively but not change the model structure. The main content is as follows: (1) The definition of GTP model smoothness and the relevant method are proposed to measure the smooth degree of the GTP model, which is used to deicide whether the GTP voxel needs interpolation. (2) The design of virtual boreholes and the corresponding method for the interpolation of the GTP geological model according to the GTP geometric smoothness are discussed. (3) The interpolation method for the adaptive GTP geological model based on the design and calculation of the virtual boreholes, can be studied based on the existing GTP building method.
    • 自然界中存在着大量的饱和边坡,例如强降雨引起堰塞坝、冰湖冰碛坝暂态性饱和以及水底边坡等,此类边坡不同于一般的无水边坡,当地震发生时,饱和边坡在地震和地下水共同作用下会产生超静孔隙水压力从而导致边坡发生滑塌。例如,2003年7月日本宫城县地震中,大量岸坡出现严重损坏,震前强降雨导致边坡内部接近饱和,使得岸顶部出现了约3 m左右的高低差,坡脚出现了3~5 m的堆积体。由此可知,饱和边坡内由于含水率高,加之地震诱发的超静孔隙水压力作用,在地震作用下更容易发生滑塌破坏,确定饱和边坡地震稳定性分析的简易评价方法,从而采取有效的防治措施,对于阻止灾害发生具有重要意义。

      目前,有限元方法是边坡地震滑移变形分析中常用的数值分析方法[1-2]。然而,由于在有限元法中需要使用网格划分方法,当目标变形较大时,存在网格扭曲且精度大幅降低,或难以继续进行分析的问题。因此,有限元方法在预测滑移大变形方面还存在不足。例如,黄帅等[3]基于砂质边坡的弹塑性有限元模型研究了地下水对边坡地震动力响应的影响规律,但是没有给出边坡的滑移形式;董士杰等[4]在不考虑地下水影响时,给出了地震作用下加筋土边坡变塑性分布情况,但是其不能反映边坡的大变形过程;同样,王飞等[5]和Yin等[6]研究了地震作用下边坡的变形和破坏特征,但是均只能给出边坡的最大剪应变和变形的云图,没能实现边坡滑移大变形模拟。而光滑粒子流体动力学方法(SPH)采用拉格朗日坐标系的运动方程进行计算,不需要创建网格,因此它具有适用于大变形问题的优点。近年来,SPH在弹塑性体大变形分析方面取得了一定进展,已被应用于岩土工程领域[9],例如唐宇峰等[10]研究发现SPH计算的边坡安全系数与极限平衡法及有限元法十分接近,均可作为失稳判据。Nonoyama等[11]进一步推动了SPH方法在边坡稳定性分析中的应用。Zhang等[12]提出了一种基于SPH的黏性土边坡滑坡模拟方法,并与有限元方法进行对比,其模拟效果吻合较好。张卫杰等[13]基于Biot固结原理提出了一种基于正则化修正的水土耦合SPH数值模型,并将该模型用于模拟三维水体溃坝模拟,但是该改进模型采用弹性模型描述土骨架的应力-应变关系,并不涉及外部动力荷载。另外,还有研究者采用SPH方法,提出了适用于土体大变形流滑分析的三维SPH仿真模型,并应用于土体流滑特性的模拟[14-16]。因此,由于SPH方法具有大变形分析的优点,已被国内外学者应用在了边坡的大变形分析中,并验证了其可行性;然而对于地震作用下边坡的临界滑移面形状及其大变形特性分析,尤其是含有地下水的边坡在地震作用下会产生超静孔隙水压力,如何在SPH方法中进行考虑,需要进一步深入研究。

      综上,本研究通过在SPH方法的地震响应分析中引入了有效应力的本构模型和瑞利衰减,建立了基于改进SPH方法的边坡地震滑移大变形的数值模型。之后,进行了饱和边坡地震滑移行为的再现分析,并研究了改进的SPH方法中粒子密度对饱和边坡地震滑移面和破坏模式的影响规律,最后将本文提出的数值模型与Newmark方法进行了对比分析。

      SPH方法中核近似是通过以下公式进行计算[17]

      f(x)=f(x)W(|xx|,h)dx (1)

      式中,W为权函数,xx为粒子之间距离,h为平滑距离。核函数要求在整个区域中积分的值为1以及h→0的极限为狄拉克δ函数。而且在式(1)中,权函数W被狄拉克δ函数代替的那个数值是关于δ函数之和相关的恒等式,核近似可以解释为近似核等式。

      将式(1)的积分式通过粒子近似的方法进行离散化:

      f(xi)=Nj=1mjρjf(xj)W(|xjxi|,h) (2)

      式中,xi为SPH粒子坐标,N为距离粒子i的距离内存在的粒子数。虽然SPH粒子在它们的坐标值xi处有质量mi,密度ρi,体积mj/ρj,函数值f(xi),但没有明确的形状。即,SPH方法中的粒子是与有限元法中的高斯点和节点相当的概念,而且与离散单元法中的粒子完全不同。通常,在对连续介质进行分析时,除了函数f(x),还要对其微分形式f(x)进行评价[18]。根据式(2)可以得到

      f(x)=f(x)W(xx,h)dx (3)

      此时,

      f(x)W(xx,h)={f(x)W(xx,h)}f(x)W(xx,h) (4)

      所以可以将式(3)转化为

      f(x)={f(x)W(xx,h)}dxf(x)W(xx,h)dx (5)

      根据散度定理,可得

      {f(x)W(xx,h)}dx=sf(x)W(xx,h)ˆndS=0 (6)

      因此,函数f(x)微分形式f(x)内核近似为

      <f(x)>=f(x)W(xx,h)dx (7)

      利用粒子将式(7)离散化,可得到

      <f(x)>=Nk=1f(xk)ρ(xk)W(xxk,h)mk (8)

      针对W函数目前已经提出了很多形式,本文采用比较成熟的三次样条函数[12]。内核函数W(xxk,h)满足

      limh0W(xx,h)=δ(xx) (9)
      W(x,h)dx=1 (10)

      在计算平滑距离h内的点的物理量f(xi)时,该粒子平滑距离h内的粒子的值的加权平均作为核函数的值,相当于求出一个权重;超出h的粒子,核函数值为零,不计入计算[19]。对于平滑距离h计算精度会随着h的值而变化问题,岩本哲也等[20]应用Chen等[21]提出的Corrective Smoothed Particle Method(CSPM)方法,可以大幅提高分析精度对h的依赖性。因此,在本研究中也将采用CSPM方法代替原始SPH方法的公式化进行计算。在CSPM方法中,f(xi)及其微分形式为

      f(xi)=Nj=1mjρjf(xj)WijNj=1mjρjWij (11)
      f(xi)xαNj=1mjρj(xjαxiα)Wijxβ=Nj=1mjρj{f(xj)f(xi)}Wijxβ (12)

      岩土体等弹性体或弹塑性体的运动方程为

      ¨u(z)=1ρ(z)dτ(z)dz (13)

      式中,¨u(z)为深度z上水平方向加速度,τ(z)为简剪切应力,ρ(z)为密度。如果是弹性体,可以用剪切刚性G与剪切应变γ(z)以来表示剪切应力τ(z)

      τ(z)=Gγ(z) (14)

      或者用增量的形式来表示:

      dτ(z)=G0dγ(z) (15)

      另外,剪切应变γ(z)可以表示为

      γ(z)=du(z)dz (16)

      本文中,笔者采用SHP方法针对式(13)近似求解。首先,通过内核近似[22]来求du(z)/dt。此时,考虑到通常存在的数值稳定性问题,不把u(z)直接代入式(13),而是根据

      f(x)=1ρ[{ρf(x)}f(x)ρ] (17)

      进行改写的部分运用内核近似。

      关于粒子a的剪切应变利用式(18)求解:

      γa=duadz=bmbρb(uaub)dWabdz (18)

      式中,下标a,b分别为求物理量的粒子与包含于其影响半径内的粒子。另外,其核函数满足

      Wab=W(xaxb,h) (19)

      因此,由式(18)求得的剪切应变γa,通过式(14)来求剪切应力τa

      接下来,针对式(13)运用恒等式(20)求边坡不同深度的加速度:

      f(x)=ρ(f(x)ρ+f(x)2ρ) (20)
      ¨u(z)=bmb(τaρ2a+τbρ2b)dWabdz (21)

      在求得粒子a的加速度后,根据得到的¨ua来求粒子a的移动量,更新密度ρa

      ρa=bmbWab (22)

      如果对于非饱和边坡在地震作用下需要考虑超静孔隙水压力的作用,因此需要采用有效应力分析,土的应力-应变关系属于非线性模型,则可以通过将式(14)或式(15)转换为有效应力模型,然后运用SPH方法对边坡模型进行分析。

      为了能模拟地震作用下饱和边坡的动力行为,采用有效应力法进行加速度的求解。对于松散砂土循环剪切试验时,超静孔隙水压力随着剪切力的增加而单调增加;而密砂情况下,砂土在达到一定的剪应力比时变得松散,此时,超静孔隙水压力像松散砂土一样随循环剪切力而单调增加;但当剪应力比超过某一应力比时,因剪切而产生的正的剪胀性会导致超静孔隙水压力减小。该应力比一般被称为相变角,在循环剪切过程中,数值一定。如果是紧密砂土,在达到相变角后,随着剪切应力比的增加,由于剪胀效应,孔隙水压力会减小,但如果不超过一定数值,当剪切应力比大到一定程度后,水压力不变。在此项研究中,有效应力的本构模型,采用了社本模型[23]

      有效应力模型的应力-应变关系采用Ramberg- Osgood模型[24],其可用于非线性分析的三参数模型。在大变形过程中,有效应力会随着循环剪切的进行而发生变化。因此,应力-应变关系可以用增量形式表示为

      dγ=dτG0(1+αβ|ΔτG0|β1) (23)

      式中:γ为剪切应变;τ为剪切应力;Δτ为剪切应力卸载时的平均剪切应力增量(参见图 1)。该关系由增量显示表示,如式(24),(25)所示。αβ分别表示由方程(24)给出,G0是由方程(25)计算的初始剪切弹性系数。另外,G0i在骨架曲线上的值为G0,在滞回曲线上的值为2G0

      α=(2γrf)β1β=2+π hmax2π hmax (24)
      G0=G0i(σmσmi)mγrf=γrfi(σmσmi)mr (25)
      图  1  有效应力路径
      Figure  1.  Effective stress path

      式中:σm为有效围压,σmi为初始有效围压;γrfσm中的标准应变;hmax为最大衰减常数,mmr为初期剪切刚度及标准剪切应变围压依赖性的常数,可取0.5。

      对于社本模型,在发生循环剪切现象的阶段,超静孔隙水压力上升程度的变化以及循环剪切的进行,要小于以往的试验结果,所以本文采用初始有效围压σmi替有效围压σm

      在大变形分析中,目标物体中发生的旋转有时是不能忽略的。由于简单地将柯西应力对时间微分得到的应力速度不具有客观性,因此存在应力值因坐标轴的旋转而明显变化的问题。Gray等[25]引入Jaumann应力速度来消除这种明显的应力变化:

      ˆσαβ=˙σαβωαγσγβ+σαγωγβ (26)

      式中,ˆσαβ为Cauchy应力增量,ωαβ为一个自旋张量,

      ωαβ=12(ναxβνβxα) (27)

      由于岩土体结构存在阻尼,因此地震波传播过程中会存在衰减行为,在地震工程学领域广泛使用的有限元法的地震响应分析中,大部分情况下使用瑞利衰减。在有限元法中,瑞利衰减矩阵[C]使用αRβR表示为

      [C]=αR[M]+βR[K] (28)

      式中:[M]为质量矩阵;[K]为刚度矩阵。在有限元法中,质量矩阵中有匹配质量矩阵和集中质量矩阵,当分析对象的物体具有材料非线性时,式(28)中的刚度矩阵[K]对应于切线刚度矩阵。在集中质量矩阵中,矩阵的对角分量为正值,所有非对角分量均为零。因为对角分量的值是分配给相应节点的质量,所以在SPH方法中可以认为是对应粒子的质量。

      在有限元方法中,式(28)的衰减矩阵作用于节点速度矢量。即,当时间t处的节点速度矢量为{vt}时,由瑞利衰减引起的衰减力矢量为

      [C]{vt}=αR[M]{vt}+βR[K]{vt} (29)

      式(39)的右边第二项满足

      [K]{vt}[K]{ut}{utΔt}Δt={ft}{ftΔt}Δt (30)

      式中:{ut},{ft}分别表示时刻t处的节点位移矢量、恢复力矢量;Δt为时间增量。由此,式(29)的右边第二项可以解释为在从要素的应力状态得到的节点力的时间增量乘以系数β

      综上,为了在SPH方法中引入瑞利衰减,在式(13)运动方程的右边增加衰减项:

      ¨u(z)=1ρ(z)dτ(z)dz + αRνα+βRΔgαΔt (31)

      在SPH方法中,每个时间步长通过运动方程(31)求出粒子的加速度,并更新速度和粒子坐标。在时间t中,当粒子的加速度为¨u(z,t),速度用˙u(z,t)表示,位移用u(z,t)表示时,下一个时间步长t+Δt中的粒子速度˙u(z,t+Δt)、粒子位置u(z,t+Δt)可以通过下式得到[26]

      ˙u(z,t+t)=˙u(z,t)+¨u(z,t+t)t (32)
      u(z,t+Δt)=u(z,t)+˙u(z,t+Δt)Δt (33)

      综上,在光滑粒子流体动力学方法(SPH方法)中加入瑞利衰减阻尼,基于有效应力本构模型改进SPH方法的运动方程,然后基于改进的开源代码程序即可实现地震作用下饱和边坡的大变形模拟。SPH方法的计算流程图如图 2所示。

      图  2  SPH方法的计算流程图
      Figure  2.  Flow chart of SPH method

      图 2可知,基于SPH方法的计算在初始运行阶段会从输入文件中读入包含有粒子类型、坐标、速度、半径、密度、体积和压强等物理属性信息。随后系统进入时间步长积分的主循环阶段,时间步数istep从0开始每循环一次增加1,每当系统完成一次循环后,会判断是否满足可视化条件或输出条件,如果满足则系统将绘制出计算结果的粒子分布图和等值线图,如果满足输出条件则系统将输出并储存该时刻的相关数据信息。最后系统将判断是否满足循环条件,从而决定是否继续循环或者结束程序。

      以笔者参与的某构筑物地震安全性评价项目中的饱和基坑边坡为研究对象,由于需要进行抗震加固设计,因此需要确定地震作用下的滑移面位置及主要的滑移破坏形式。土体以砂质土为主。坡长为17 m,高度约为12 m,坡度45°,概化分析模型如图 3所示。

      图  3  边坡概化模型
      Figure  3.  Conceptual model for slope

      通过现场取样,开展了动态空心圆柱扭剪试验(图 4),并结合地质勘察报告确定了该边坡的砂质土物理力学参数取值,如表 1所示。

      图  4  室内试验测试
      Figure  4.  Laboratory tests
      表  1  土体物理力学参数
      Table  1.  Physical and mechanical parameters of soils
      泊松比 弹性量/
      MPa
      重度/
      (kN·m-3)
      黏聚力/
      kPa
      内摩擦角/(°) 残余黏聚力/
      kPa
      残余内摩擦角/(°)
      0.3 100.4 16.8 50.42 15 5.95 15
      下载: 导出CSV 
      | 显示表格

      在SPH分析中,为了确定数值模型中粒子密度对分析结果的影响,构建了3种粒子密度不同的模型,SPH粒子的间隔dp和每个模型所用的SPH粒子的总数np图 5所示。边界条件的设定:底部为固定约束,侧面为水平向约束,垂直方向自由。模型的平滑距离h = 2.0dp,衰减阻尼通过有限元方法首先获得该模型的第一阶和第二阶的自振频率,进而可以求出瑞利阻尼。在SPH数值积分过程中所用的时间步是通过反复试验设定,当粒子在间隔最小的情况下也能得到稳定的分析结果时候的值为1.0×10-5 s。

      图  5  SPH计算模型
      Figure  5.  Computational model for SPH

      本文选用地震波为El Centro地震波,该地震波震级为7.1级,震中距为12 km,记录的最大峰值加速度为3.42 m/s2。其是用1940年位于加州南部的埃尔森特罗的一台强震仪由人类记录的第一条地震波,记录的地震动时程如图 6所示。

      图  6  地震动记录
      Figure  6.  Ground motion recording

      为了验证基于SPH方法的有效应力分析的合理性,针对该边坡砂质土取样(直径为50 mm,高度为100 mm),并进行了空心圆柱扭剪试验,与本文改进的SPH方法计算结果进行了对比。其中模型的粒子数为50个,粒子间距0.01 m,施加频率2 Hz的正弦波作为剪切力,时间为10 s。该循环剪切试验中的有效应力路径和剪切应力-应变关系如图 78所示。

      图  7  有效应力路径对比
      Figure  7.  Comparison of effective stress paths
      图  8  剪切应力-应变关系对比
      Figure  8.  Comparison of shear stress-strain relationship

      图 78可知,基于改进的SPH方法计算的有效应力路径和剪应力-应变关系,从初始状态到超静孔隙水压力上升、有效应力减小的过程得到了很好的再现。且有效应力归零后剪切应变急剧增加的循环滞回现象也基本上得到再现,也就是说此时土体发生了振动液化现象,有效应力降低为0,而应变会急剧增加产生大变形。从图 8还可以发现,平均主应力在25~50 kPa,SPH与试验结果差距较大,主要原因是由于动载荷作用过程中方向的改变,以及动三轴试验过程测试以及操作误差等原因,总体而言,测试值和数值结果误差在允许范围之内。另外,文中平均主应力是恒定的,排水过程中应力路径的变化对砂土的抗剪强度几乎没有影响,但会影响小变形范围内的应力-应变特性,表现为不同的弹性轨迹,可参见参考文献[1]。因此,对于模拟中的大变形影响不大。另外,剪应力-应变的循环滞回现象有一些差异,图 8(a)试验中的滞回曲线呈现V形,而图 8(b)数值方法中的滞回曲线呈现倒S形,但整体是吻合较好的。

      本节将通过SPH方法再现边坡的地震滑移破坏行为,确定出边坡的临界滑移面,并于与Newmark方法的对数螺旋临界滑移面形状进行对比。采用Newmark方法计算该边坡的临界滑移面。基于变分原理对边坡的对数螺旋状临界滑面的计算流程进行推导[28],如图 9所示。

      图  9  边坡临界滑移面计算示意图
      Figure  9.  Calculation diagram of critical slip surface of slope

      图 9所示,在进行边坡螺旋曲线滑移面的计算时,把滑移体简化为不同的土柱体,其相对转动中心的阻力矩和下滑力矩满足转动平衡方程:

      M=MRMD=MR(MDV+MDH)=0 (34)

      式中:MR为滑移面的摩擦阻力矩;MD为滑动转矩,其为滑移面土柱自重而形成的滑动转矩MDV与滑移面土柱惯性力而形成的滑动转矩MDH之和;H为边坡的高度。

      图 9中的边坡高度H为基础,建立其标准化坐标系,以坡脚为坐标原点,可以求得不同的转动力矩:

      MR=1FSni=1{τisinθidX(XiXc)}τicosθidY(YcYi)=ni=1(Nm+S(Yi)φm)[(YcYi)(XiXc)Yi]dXi = ni=1cm[(YcYi)(XiXc)Yi]dXi (35)
      MDV=ni=1[(¯YiYi)(XiXc)]dXi (36)
      MDH=a2ni=1[(¯YiYi)(Yi + ¯Yi2Yc)]dXi (37)
      Xi=xiHYi=yiH¯Yi=¯yiHYi=dYidXi (38)
      cm=cFs1γHψm=tanφFs (39)

      式中:XiYi¯Yi分别为以边坡高H来进行标准化的滑移面位置xy的坐标;τi为滑移面切向剪应力;Fs为边坡稳定性安全系数;Y为临界滑移面坐标Xi的斜率;c为滑移面的黏聚力;S(Yi)为滑移面的垂直应力;φ为内摩擦角;cmφm分别为标准化的黏聚力与内摩擦角;γ为滑移面上土体重度;a为作用于dXi区间土柱的地震加速度。

      地震与地下水的耦合作用会产生超静孔隙水压力,因此在转动平衡公式(34)中需要考虑超静孔隙水压力诱发的滑移体的下滑力的作用。本文超静孔隙水压力usis的计算式主要采用Huang等[29]中提出的简化计算方法。临界滑移面位置处的孔隙水压力计算式为

      u=u0+usis (40)

      式中,u为孔隙水压力,u0为静孔隙水压力,usis为超静孔隙水压力。

      进而转动平衡方程(40)变为

      M=MRMD=MR(MDV+MDH+MDu)=0 (41)

      式中,MDu为孔隙水压力产生的转动矩,其中MDu的表达式为

      MDu=ni=1{ucosθidX(XiXc)+usinθidY(YcYi)} (42)

      对数螺旋临界滑面的形状,计算表达式为[30]

      X=Xc+Aexp(βψm)sinβ (43)
      Y=YcAexp(βψm)sinβ (44)

      式中:XCYC分别为对数螺旋极坐标;A为对数螺旋常数;β为与斜面位置相关的角。此处β为以对数螺旋极坐标为中心并沿铅直方向至滑移面上的点(XY)的逆时针方向旋转的角。屈服加速度asy可由式(35)中边坡稳定性安全系数Fs为1.0时,代入平衡方程式(34)求得。

      通过Newmark方法确定了边坡的临界滑移面如图 10所示。

      图  10  由Newmark方法计算的滑移面
      Figure  10.  Slip surface calculated by Newmark method

      图 10可知,该边坡滑移面为典型的对数螺旋形滑移面,且剪出口位于坡脚位置。

      接下来采用SPH方法计算边坡滑移面形状。其中每个粒子的本构方程是有效应力模型。内摩擦角和黏聚力在首次达到断裂状态之前使用表 1中的峰值强度,之后使用残余强度值。在分析过程中,不同峰值加速度时边坡的状态是:在施加小于1 m/s2峰值加速度地震波时,边坡不会发生大的变形,在施加到2,3 m/s2峰值加速度时各个点逐渐出现位移,在施加到4 m/s2峰值加速度时,边坡内部出现滑移线,并导致完全滑坡。考虑到水平地震是边坡滑塌的主要影响因素,因此本文只考虑水平地震。

      根据不同粒子密度,分别计算3个模型在地震作用下的滑移面以及滑塌模式,如图 11所示。

      图  11  地震持续10 s时滑移面
      Figure  11.  Slip surface under earthquake at time 10 s

      图 11所示,从地震作用开始10 s后不同模型出现的剪切应变最大的区域可以看出,在dp为0.5 m时的滑移面与Newmark方法计算的滑移面出现位置非常相似;而dp为2 m时剪切应变最大分布区域在坡顶面比Newmark方法位置更深的位置处。dp为5 m时未能出现明显的滑移面。基于此,进一步计算了地震作用结束后的边坡滑移面分布情况,如图 12所示。

      图  12  地震结束后滑移面
      Figure  12.  Slip surface after earthquake

      图 12所示,地震作用结束后,3个模型的滑移面形状大致相同。在dp为5,2 m的情况下剪应变最大值分布超过80%的区域几乎出现在相同的地方,与Newmark方法相比,其滑移面在接近坡顶位置更偏向内部,接近坡脚位置滑移面形状吻合较好,说明粒子密度较小的情况下,容易高估边坡的滑移规模。在dp为0.5 m的情况下,滑移面与Newmark方法计算的滑移面整体吻合较好,且能较好地模拟边坡地震作用下滑移大变形过程。因此,在使用SPH方法进行边坡地震破坏分析中,选择合适的粒子密度是准确模拟饱和边坡地震滑移大变形的关键。分析其原因主要是由于对于SPH方法,粒子运动是由压力梯度的作用产生,粒子压力是通过状态方程由粒子自身的密度和内能来计算。SPH方法考虑了相邻粒子之间的影响作用,因此粒子的运动速度和相邻粒子的平均速度相近。在计算过程中对连续密度场的估计值会受到粒子分布的显著影响,即粒子质量分布的聚集/稀疏区域。

      综上,本文改进的SPH方法可以较好地模拟饱和边坡地震临界滑移面的形状,且能再现饱和边坡地震作用下的滑移大变形过程。

      (1)本文通过在光滑粒子流体动力学方法(SPH方法)中加入瑞利衰减阻尼,基于有效应力本构模型,对SPH方法进行了改进,建立了用于饱和边坡地震滑移大变形的数值计算方法;并将该方法与饱和土样的空心圆柱扭剪试验结果进行了对比,得出改进后的SPH方法能很好地再现饱和土体的有效应力路径和剪应力-应变滞回性能,说明其适用于饱和土体的地震大变形计算。

      (2)将超静孔隙水压力简便计算方法代入Newmark方法中,确定了考虑超静孔隙水压力影响临界滑移面计算方法,并与改进的SPH方法确定的饱和边坡地震滑移临界面形状与进行了对比,得出本文改进的SPH方法能较好地再现饱和边坡地震作用下的滑移过程。

      (3)分析了SPH数值计算方法中粒子密度对边坡地震滑移大变形的影响。结果表明,采用SPH方法进行边坡地震破坏分析时,结果受设置的粒子密度的影响。在dp = 0.5 m的情况下获得了非常相似的结果,而在dp = 5 m时,在滑移面的发展过程中出现了一定差异,但形状整体吻合较好。

    • 图  1   广义三棱柱地质模型

      Figure  1.   GTP model built by boreholes

      图  2   广义三棱柱平滑度

      Figure  2.   Smoothness of GTP

      图  3   精细化三维地质模型构建方法

      Figure  3.   Building method for smooth GTP model

      图  4   基于52个钻孔构建的地质曲面和模型,增加138个虚拟钻孔后的地质曲面和本文方法构建的精细化模型

      Figure  4.   The GTP surface and model based on 52 borehole data, the surface based on the 138 virtual borehole added by the adaptive GTP interpolation method and the model

    • [1]

      WANG Z, QU H, WU Z, et al. Formal representation of 3D structural geological models[J]. Computers & Geosciences, 2016, 90: 10-23.

      [2]

      JACQUEMYN C, JACKSON M D, HAMPSON G J. Surface-based geological reservoir modelling using grid-free NURBS curves and surfaces[J]. Mathematical Geosciences, 2019, 51(1): 1-28. doi: 10.1007/s11004-018-9764-8

      [3]

      Abdul-Rahman, Alias, and Morakot Pilouk. Spatial data modelling for 3D GIS[M]. Springer Science & Business Media, 2007.

      [4]

      WANG Y, FERNÀNDEZ-GARCIA D, SAALTINK M W. Modeling reactive multi-component multi-phase flow for Geological Carbon Sequestration (GCS) with Matlab[J]. Computers & Geosciences, 2023, 172: 105300.

      [5]

      Lixin W, Wu L, Lixin W, et al. Topological relations embodied in a generalized tri-prism (GTP) model for a 3D geoscience modeling system[J]. Computers and Geosciences, 2004, 30(4): 405-418. doi: 10.1016/j.cageo.2003.06.005

      [6] 齐安文, 吴立新. 基于类三棱柱的三维地质模拟与拓扑研究[J]. 矿山测量, 2003(3): 65-66, 64. https://www.cnki.com.cn/Article/CJFDTOTAL-KSCL200303020.htm

      QI Anwen, WU Lixin. Three-dimensional geological simulation and topological research based on quasi-triangular prism[J]. Mine Surveying, 2003(3): 65-66, 64. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-KSCL200303020.htm

      [7] 陈学习, 吴立新, 车德福, 等. 基于钻孔数据的含断层地质体三维建模方法[J]. 煤田地质与勘探, 2005, 33(5): 5-8. https://www.cnki.com.cn/Article/CJFDTOTAL-MDKT200505001.htm

      CHEN Xuexi, WU Lixin, CHE Defu, et al. 3D modeling method of geological bodies including faults based on borehole data[J]. Coal Geology & Exploration, 2005, 33(5): 5-8. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-MDKT200505001.htm

      [8]

      LI Q, LIAO G, CHEN X, et al. Approach on area coordinate, volume coordinate and their usage in true 3dgis[J]. International Conference on Geo-spatial Solutions for Emergency Management and the 50th Anniversary of the Chinese Academy of Surveying and Mapping, 2009, C(5): 158-164.

      [9] 李青元, 张洛宜, 曹代勇, 等. 三维地质建模的用途、现状、问题、趋势与建议[J]. 地质与勘探, 2016, 52(4): 759-767. https://www.cnki.com.cn/Article/CJFDTOTAL-DZKT201604018.htm

      LI Qingyuan, ZHANG Luoyi, CAO Daiyong, et al. Usage, status, problems, trends and suggestions of 3D geological modeling[J]. Geology and Exploration, 2016, 52(4): 759-767. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-DZKT201604018.htm

      [10] 车德福, 吴立新, 殷作如. 三维地质建模中地层错断处自动构建的新方法[J]. 煤炭学报, 2009, 34(10): 1305-1309. https://www.cnki.com.cn/Article/CJFDTOTAL-MTXB200910002.htm

      CHE Defu, WU Lixin, YIN Zuoru. The new method for creating 3D model automatically of bad break within geological body[J]. Journal of China Coal Society, 2009, 34(10): 1305-1309. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-MTXB200910002.htm

      [11]

      WANG W, SUN L, LI Q, et al. Representing the geological body heterogeneous property field using the quadratic generalized tri-prism volume function model[J]. Arabian Journal of Geosciences, 2017, 10(5): 115.

      [12]

      YAN B, LIU J, QI Q, et al. Study on chain relationship and risk assessment model of coal mine geological disasters[J]. Arabian Journal of Geosciences, 2022, 15(9): 900.

    • 期刊类型引用(3)

      1. 黄帅,王中根,李昱,冯宇翔. 基于改进SPH方法的滑坡涌浪对大坝结构冲击响应规律. 工程科学与技术. 2025(01): 120-131 . 百度学术
      2. 魏星,程世涛,谢相焱,陈睿. 考虑强度速率衰减效应的地震滑坡SPH-FEM模拟. 岩土工程学报. 2024(08): 1753-1761 . 本站查看
      3. 郑厚国,周永强,刘烨. 高水头作用下饱和边坡破坏物质点法模拟研究. 能源与环保. 2023(05): 276-283+288 . 百度学术

      其他类型引用(4)

    图(4)
    计量
    • 文章访问数:  158
    • HTML全文浏览量:  26
    • PDF下载量:  46
    • 被引次数: 7
    出版历程
    • 收稿日期:  2023-07-04
    • 网络出版日期:  2023-11-23
    • 刊出日期:  2023-10-31

    目录

    /

    返回文章
    返回