地震中重力式挡土墙抗滑分析的离散颗粒‑SPH耦合模拟方法与流程

本发明属于岩土工程研究领域,尤其涉及一种地震中重力式挡土墙抗滑分析的离散颗粒-sph耦合模拟方法。

背景技术：

重力式挡土墙作为岩土工程中的支护结构，在我国边坡支护中有着广泛应用。但是对其抗震设计的理论并不成熟，目前的分析方法主要是拟静力法，地震作用下重力式挡土墙有水平向加速度，即将水平向的惯性力简化为静力来分析挡土墙的抗滑稳定性。为了深入获得挡土墙在地震全过程中的抗滑稳定性，可以通过建立有限元或有限差分等连续体模拟挡土墙的地震响应，但是这类连续体模型难以分析离散颗粒与挡土墙的非连续动力接触作用。因此为了精细化分析重力式挡土墙在地震全过程中的抗滑稳定性，可以基于离散颗粒模拟土体，但由于受离散颗粒数量和计算机性能的限制，很难在全域上采用离散颗粒模拟土体的力学行为，为此人们发展了固体离散-连续土体耦合方法以减少离散元模拟规模，并用这种方法模拟挡土墙的地震响应(见文献“周健，金炜枫.基于耦合方法的挡土墙地震响应的数值模拟[j].岩土力学，2010，31(12)：3949-3957.”)。但是这些耦合模型中连续固体基于有限元或有限差分网格离散，在地震大变形模拟中存在网格畸变等数值失稳问题，因此不如无网格法及离散颗粒适应地震大变形分析，另外连续土体和离散颗粒在耦合边界上容易分离。

光滑粒子动力学(smoothparticlehydrodynamicsmethod,sph)是无网格法的一种，可以适应土体的地震大变形模拟，因此可以考虑用sph方法与离散颗粒耦合，靠近挡土墙的土体用离散颗粒模拟，远离挡土墙的土体用sph方法模拟，这样可以考虑挡土墙与离散颗粒的非连续接触作用，又可以有效降低离散元模拟规模，同时连续土体也可以适应大变形分析，但是尚需构建跨尺度衔接方法来耦合挡土墙周围离散颗粒和sph连续土体，从而保证连续土体和离散颗粒在耦合过渡域上不会分离。

技术实现要素：

本发明为了克服现有离散颗粒-连续土体耦合方法在精细化分析地震中重力式挡土墙抗滑稳定性时，难以顾及连续土体的大变形问题、以及连续土体和离散颗粒在耦合边界上容易分离的问题，本发明提出了一种地震中重力式挡土墙抗滑分析的离散颗粒-sph耦合模拟方法，从而不仅可以考虑挡土墙与离散颗粒的非连续接触作用并有效降低离散元模拟规模，同时克服已有方法中连续土体难以适应大变形分析的困难、以及连续土体和离散颗粒在耦合边界上会分离的问题。

本发明涉及一些简写和符号，以下为注解：

sph：指无网格法中的光滑粒子动力学方法(smoothparticlehydrodynamicsmethod,sph)，利用这种方法建立的连续土体模型可以适应地震大变形分析。

{u^c}int：耦合过渡域上的sph质点位移组合的列向量，这里用上标c表示连续土体上的变量，下标int表示耦合。

{f^c}int：耦合过渡域上的sph质点受到离散颗粒的耦合力组合的列向量。

{u^c}uin：不在耦合过渡域上的sph质点位移组合的列向量，这里下标uin表示不耦合。

不在耦合过渡域上的sph质点速度组合的列向量。

不在耦合过渡域上的sph质点加速度组合的列向量。

{f^c}uin：不在耦合过渡域上的sph质点受到的外力列向量。

[m]：sph离散连续固体微分方程得到的质量矩阵。

[c]：sph离散连续固体微分方程得到的阻尼矩阵。

[k]：sph离散连续固体微分方程得到的刚度矩阵。

耦合过渡域上第i个离散颗粒的平动位移列向量，这里上标m表示离散颗粒上的变量。

耦合过渡域上第i个离散颗粒的平动加速度列向量。

耦合过渡域上n个离散颗粒的位移列向量联立起来组成的总体位移列向量。

{fi^m}contact：耦合过渡域上第i个离散颗粒受到其他离散颗粒的作用合力列向量，这里下标contact表示颗粒接触。

{f^m}contact：耦合过渡域上n个离散颗粒的{fi^m}contact联立起来组成的列向量。

{fi^m}int：耦合过渡域上第i个离散颗粒受到第i个sph质点耦合力列向量。

{f^m}int：耦合过渡域上n个离散颗粒的{fi^m}int联立起来组成的列向量。

[m]：耦合过渡域上n个离散颗粒的平衡方程联立起来形成的总体质量矩阵。

n^slip：挡土墙与颗粒接触的侧边上有限元边界段数，这里上标slip表示滑动；

fi^slip：挡土墙与颗粒接触的侧边上，其中第i段有限元边界受到的颗粒水平推力。

f^slip：挡土墙受土体的水平推力合力，即挡土墙滑动力合力。

n^anti-slip：挡土墙底部与颗粒接触的有限元边界段数，这里上标anti-slip表示抗滑动。

fi^anti-slip：挡土墙底部第i段有限元边界受离散颗粒水平摩擦力。

f^anti-slip：挡土墙底部受土体的抗滑力合力。

k：挡土墙的抗滑稳定安全系数。

为解决上述技术问题采用以下技术方案：

本发明提供一种地震中重力式挡土墙抗滑分析的离散颗粒-sph耦合模拟方法，包括散颗粒集合体宏观参数和连续土体匹配数值试验，分区域建立挡土墙有限元模型、挡土墙附近离散颗粒模型、远离挡土墙的sph连续土体模型，同时离散颗粒和连续土体之间有部分重叠的耦合过渡域，从而实现地震中重力式挡土墙的抗滑稳定安全系数计算，具体步骤如下：

步骤1：颗粒集合体宏观参数和连续土体匹配数值试验；

步骤2：建立挡土墙有限元模型、建立挡土墙附近的土体离散颗粒模型、建立远离挡土墙的sph连续土体模型；

步骤3：建立离散颗粒和连续土体之间跨尺度衔接；

步骤4：挡土墙有限元模型与离散颗粒模型耦合；

步骤5：计算地震中挡土墙的抗滑稳定安全系数。

其中，

(1)所述步骤1中，颗粒集合体宏观参数和连续土体匹配数值试验包括：

对重力式挡土墙周围取土样，进行室内三轴试验确定土体的内摩擦角和黏聚力，进行室内动三轴试验确定土体的动模量，这里得到的内摩擦角、黏聚力和动模量赋值给连续土体模型使用；不断调整离散颗粒的细观参数进行离散颗粒的双轴压缩试验和循环双轴数值试验，让双轴压缩数值试验得到的内摩擦角和黏聚力和实际试验一致，让循环双轴数值试验得到的动模量和实测动模量一致，从而确定颗粒的细观参数如法相刚度kn、切向刚度ks和摩擦系数fc，将这些细观参数赋值给离散颗粒模型使用。

(2)所述步骤2中，挡土墙有限元模型、土体离散颗粒模型、sph连续土体模型分述如下：

(2.1)挡土墙有限元模型用有限元方法中的四节点等参单元。

(2.2)土体离散颗粒模型中，颗粒的运动方程包含平动方程和转动方程。

(2.3)sph连续土体模型中，用光滑粒子动力学(smoothparticlehydrodynamicsmethod,sph)建立连续土体模型，所用sph核函数采用已有的高斯型核函数。

(3)所述步骤3中，建立离散颗粒和连续土体之间跨尺度衔接，具体方法如下：

离散颗粒和连续土体之间有部分重叠的耦合过渡域，耦合过渡域上每个离散颗粒和连续土体sph质点一一对应耦合且位移、速度和加速度相等，并建立离散颗粒和连续土体sph质点耦合过渡域之间的耦合力算法如下：

设耦合过渡域上共有n个sph质点和对应的n个离散颗粒，然后将连续土体sph质点分为两部分，在耦合过渡域上的n个sph质点位移组合为列向量{u^c}int，sph质点受到离散颗粒的耦合力组合为列向量{f^c}int，这里用上标c表示连续土体上的变量，下标int表示耦合；对不在耦合过渡域上的sph质点，其位移、速度、加速度和外力分别组合为列向量{u^c}uin、和{f^c}uin，这里下标uin表示不耦合；同时用[m]、[c]和[k]分别表示sph离散连续固体微分方程得到的质量矩阵、阻尼矩阵和刚度矩阵。然后基于sph离散可以得到动力方程：

设耦合过渡域上第i个离散颗粒的质量为mi，平动位移列向量为相应加速度列向量为受到其他离散颗粒的作用合力为列向量{fi^m}contact，受到第i个sph质点耦合力列向量为{fi^m}int，这里上标m表示离散颗粒上的变量，这里下标contact表示颗粒接触；不考虑颗粒转动，单个离散颗粒的平衡方程为：

将耦合过渡域上n个离散颗粒的平衡方程联立起来，设总体质量矩阵为[m]，n个颗粒的{fi^m}contact和{fi^m}int分别联立组合为{f^m}contact和{f^m}int，这样耦合过渡域上离散颗粒的总体平衡方程可写为：

按如下步骤求解耦合力：

(a)通过消去未知耦合力求耦合质点上的位移

耦合过渡域上sph质点和离散颗粒加速度相同且耦合力为一对反力，有{f^c}int＝-{f^m}int，将式(3)代入式(1)可消去未知的耦合力{f^c}int和{f^m}int，如下所示：

上式中[mint]为[m]代入[m]后的耦合质量矩阵，上式中已消去未知的耦合力{fi^c}int和{fi^m}int，可求解耦合过渡域中sph质点的加速度速度和位移{u^c}int。

(b)求耦合过渡域上的耦合力

将求解得到的速度和位移{u^c}int代入式(1)可以计算得到作用于sph质点上的耦合力列向量{f^c}int，然后可得离散颗粒质点上作用的耦合力列向量{f^m}int＝-{f^c}int。

(4)所述步骤4中，挡土墙有限元模型与离散颗粒模型耦合的具体方法为：

通过保证耦合边界上力和速度的连续性来实现有限元-离散颗粒耦合，在挡土墙有限元模型与离散颗粒接触面上，接触面为有限元模型的线边界，接触面在离散元模型中为离散颗粒的墙边界，提取接触面上有限元节点的速度作为离散元模型中墙边界节点的速度；有限元模型统计接触面上受颗粒的合力fx、fy和合力矩mxy，将fx、fy和mxy转化为节点力，作为力边界施加至挡土墙有限元模型边界节点上：

(5)所述步骤5中，计算地震中挡土墙的抗滑稳定安全系数的具体方法为：

在挡土墙-离散颗粒-sph连续土体耦合体系的底部施加载荷，即在连续土体底部加载地震加速度，设挡土墙与颗粒接触的侧边上共有n^slip段有限元边界，设其中第i段有限元边界受到的颗粒水平推力为fi^slip，设挡土墙受土体的水平推力合力f^slip，此合力即为挡土墙滑动力合力，有这里上标slip表示滑动；在挡土墙底部有n^anti-slip段有限元边界与颗粒接触，设其中第i段墙受离散颗粒水平摩擦力为fi^anti-slip，这样挡土墙受土体的抗滑力合力为这里上标anti-slip表示抗滑动；这样可计算挡土墙的抗滑稳定安全系数为抗滑力合力和滑动力的比值，设抗滑稳定安全系数为k，有

本发明采用以上技术方案与现有技术相比，具有以下技术效果：同已有挡土墙地震响应的离散颗粒-连续土体耦合方法相比，有两个优势(a):本方法采用离散颗粒和连续土体上用光滑粒子动力学(smoothparticlehydrodynamicsmethod,sph)方法模拟，更适应大变形分析，克服已有耦合方法中有限元或有限差分网格在大变形分析中网格畸变难以适应地震大变形分析的问题；(b)本方法在离散颗粒-sph连续土体中设置耦合过渡域，耦合的离散颗粒和sph质点位移相同，克服已有方法中离散颗粒和连续土体在耦合边界上会脱离的问题。

附图说明

图1为挡土墙与土体原型图；

图2为挡土墙-离散颗粒-sph连续土体耦合模型图；

图3为离散颗粒-sph连续土体耦合过渡域图；

图4为离散颗粒与挡土墙有限单元边界耦合图。

图中1.土体，2.挡土墙，3.地震载荷，4.土体离散颗粒模型，5.离散颗粒-sph连续土体耦合过渡域，6.连续土体sph质点，7.挡土墙有限元模型。

具体实施方式

为了使本发明实现的技术手段、创新特征、达成目的与功效易于明白了解，下面结合具体图示，进一步阐述本发明。

本发明涉及一些简写和符号，以下为注解：

sph：指无网格法中的光滑粒子动力学方法(smoothparticlehydrodynamicsmethod,sph)，利用这种方法建立的连续土体模型可以适应地震大变形分析。

{u^c}int：耦合过渡域上的sph质点位移组合的列向量，这里用上标c表示连续土体上的变量，下标int表示耦合。

{f^c}int：耦合过渡域上的sph质点受到离散颗粒的耦合力组合的列向量。

{u^c}uin：不在耦合过渡域上的sph质点位移组合的列向量，这里下标uin表示不耦合。

不在耦合过渡域上的sph质点速度组合的列向量。

不在耦合过渡域上的sph质点加速度组合的列向量。

{f^c}uin：不在耦合过渡域上的sph质点受到的外力列向量。

[m]：sph离散连续固体微分方程得到的质量矩阵。

[c]：sph离散连续固体微分方程得到的阻尼矩阵。

[k]：sph离散连续固体微分方程得到的刚度矩阵。

耦合过渡域上第i个离散颗粒的平动位移列向量，这里上标m表示离散颗粒上的变量。

耦合过渡域上第i个离散颗粒的平动加速度列向量。

耦合过渡域上n个离散颗粒的位移列向量联立起来组成的总体位移列向量。

{fi^m}contact：耦合过渡域上第i个离散颗粒受到其他离散颗粒的作用合力列向量，这里下标contact表示颗粒接触。

{f^m}contact：耦合过渡域上n个离散颗粒的{fi^m}contact联立起来组成的列向量。

{fi^m}int：耦合过渡域上第i个离散颗粒受到第i个sph质点耦合力列向量。

{f^m}int：耦合过渡域上n个离散颗粒的{fi^m}int联立起来组成的列向量。

[m]：耦合过渡域上n个离散颗粒的平衡方程联立起来形成的总体质量矩阵。

n^slip：挡土墙与颗粒接触的侧边上有限元边界段数，这里上标slip表示滑动；

fi^slip：挡土墙与颗粒接触的侧边上，其中第i段有限元边界受到的颗粒水平推力。

f^slip：挡土墙受土体的水平推力合力，即挡土墙滑动力合力。

n^anti-slip：挡土墙底部与颗粒接触的有限元边界段数，这里上标anti-slip表示抗滑动。

fi^anti-slip：挡土墙底部第i段有限元边界受离散颗粒水平摩擦力。

f^anti-slip：挡土墙底部受土体的抗滑力合力。

k：挡土墙的抗滑稳定安全系数。

如图1-4所示，本发明设计一种地震中重力式挡土墙抗滑分析的离散颗粒-sph耦合模拟方法，通过颗粒集合体宏观参数和连续土体匹配数值试验，分区域建立挡土墙有限元模型、挡土墙附近离散颗粒模型、远离挡土墙的sph连续土体模型，同时离散颗粒和连续土体之间有部分重叠的耦合过渡域，基于以上内容模拟得到地震中挡土墙的抗滑稳定安全系数，具体步骤如下：

步骤1：离散颗粒集合体宏观参数和连续土体匹配数值试验：对重力式挡土墙周围取土样，进行室内三轴试验确定土体的内摩擦角和黏聚力，进行室内动三轴试验确定土体的动模量，这里得到的内摩擦角、黏聚力和动模量赋值给sph连续土体模型使用；不断调整离散颗粒的细观参数进行离散颗粒的双轴压缩试验和循环双轴数值试验，让双轴压缩数值试验得到的内摩擦角和黏聚力和实际试验一致，让循环双轴数值试验得到的动模量和实测动模量一致，从而确定颗粒的细观参数如法相刚度kn、切向刚度ks和摩擦系数fc，将这些细观参数赋值给离散颗粒模型使用；

步骤2：如图2所示，建立挡土墙有限元模型7、建立挡土墙附近的土体离散颗粒模型4、建立远离挡土墙的sph连续土体模型6：

(2.1)挡土墙有限元模型7用有限元方法中的四节点等参单元，设四个节点上的形函数分别为n1、n2、n3和n4，设形函数中的变量为r和s，等参单元四个节点上形函数具体形式为：

(2.2)土体离散颗粒模型4中，颗粒的运动方程包含平动方程和转动方程：

平动方程：

转动方程：

上式中m为颗粒质量，和为颗粒x和y方向平动加速度，fx和fx为作用在颗粒上x和y方向合力，mxy为颗粒所受xy方向转矩，i为转动惯量，为角加速度。

(2.3)sph连续土体模型6中，用光滑粒子动力学(smoothparticlehydrodynamicsmethod,sph)建立连续土体模型，所用sph核函数采用已有的高斯型核函数，设sph核函数为w(r,h)，r为两个sph质点之间的距离，h为光滑长度，核函数为w(r,h)形式如下：

步骤3：建立离散颗粒4和连续土体6之间跨尺度衔接，具体方法如下：建立离散颗粒4和连续土体6之间有部分重叠的耦合过渡域5，耦合过渡域5的特点及耦合力计算方法如下：

如图3所示，离散颗粒与连续土体部分重叠的耦合过渡域5中，连续土体用光滑粒子动力学(smoothparticlehydrodynamicsmethod,sph)模拟，要求在这个耦合过渡域5上每个连续土体sph质点与离散颗粒重合且位移相等；这时耦合问题的关键是在耦合过渡域上强化力的相容性来实现离散颗粒-连续土体的跨尺度耦合，由于在耦合过渡域上离散颗粒和连续土体sph质点对应，因此这个问题可以表述为求解离散颗粒和连续土体在耦合过渡域上的耦合力，且这个耦合力使sph质点与离散颗粒位移相同；

设耦合过渡域上共有n个sph质点和对应的n个离散颗粒，然后将连续土体sph质点分为两部分，在耦合过渡域上的n个sph质点位移组合为列向量{u^c}int，sph质点受到离散颗粒的耦合力组合为列向量{f^c}int，这里用上标c表示连续土体上的变量，下标int表示耦合；对不在耦合过渡域上的sph质点，其位移、速度、加速度和外力分别组合为列向量{u^c}uin、和{f^c}uin，这里下标uin表示不耦合；同时用[m]、[c]和[k]分别表示sph离散连续固体微分方程得到的质量矩阵、阻尼矩阵和刚度矩阵。然后基于sph离散可以得到动力方程：

设耦合过渡域上第i个离散颗粒的质量为mi，平动位移列向量为相应加速度列向量为受到其他离散颗粒的作用合力为列向量{fi^m}contact，受到第i个sph质点耦合力列向量为{fi^m}int，这里上标m表示离散颗粒上的变量，这里下标contact表示颗粒接触；不考虑颗粒转动，单个离散颗粒的平衡方程为：

将耦合过渡域上n个离散颗粒的平衡方程联立起来，设总体质量矩阵为[m]，n个颗粒的{fi^m}contact和{fi^m}int分别联立组合为{f^m}contact和{f^m}int，这样耦合过渡域上离散颗粒的总体平衡方程可写为：

按如下步骤求解耦合力：

(a)通过消去未知耦合力求耦合质点上的位移

耦合过渡域上sph质点和离散颗粒加速度相同且耦合力为一对反力，有{f^c}int＝-{f^m}int，将式(7)代入式(5)可消去未知的耦合力{f^c}int和{f^m}int，如下所示：

上式中[mint]为[m]代入[m]后的耦合质量矩阵，上式中已消去未知的耦合力{fi^c}int和{fi^m}int，可求解耦合过渡域中sph质点的加速度速度和位移{u^c}int。

(b)求耦合过渡域上的耦合力

将求解得到的速度和位移{u^c}int代入式(5)可以计算得到作用于sph质点上的耦合力列向量{f^c}int，然后可得离散颗粒质点上作用的耦合力列向量{f^m}int＝-{f^c}int。

步骤4：挡土墙有限元模型7与离散颗粒模型4耦合，具体耦合方式为：如图4所示，通过保证耦合边界上力和速度的连续性来实现有限元-离散颗粒耦合，在挡土墙有限元模型7与离散颗粒4的接触面上，接触面为有限元模型的线边界，接触面在离散元模型中为离散颗粒的墙边界，提取接触面上有限元节点的速度作为离散元模型中墙边界节点的速度；有限元模型统计接触面上受颗粒的合力fx、fy和合力矩mxy，将fx、fy和mxy转化为节点力，作为力边界施加至挡土墙有限元模型边界节点上：

步骤5：计算地震中挡土墙的抗滑稳定安全系数，具体方法为：在如图2所示的挡土墙-离散颗粒-sph连续土体耦合体系的底部施加载荷3，即在连续土体底部加载地震加速度，设挡土墙与颗粒接触的侧边上共有n^slip段有限元边界，设其中第i段有限元边界受到的颗粒水平推力为fi^slip，设挡土墙受土体的水平推力合力f^slip，此合力即为挡土墙滑动力合力，有这里上标slip表示滑动；在挡土墙底部有n^anti-slip段有限元边界与颗粒接触，设其中第i段墙受离散颗粒水平摩擦力为fi^anti-slip，这样挡土墙受土体的抗滑力合力为这里上标anti-slip表示抗滑动；这样可计算挡土墙的抗滑稳定安全系数为抗滑力合力和滑动力的比值，设抗滑稳定安全系数为k，有