毕业论文边坡数值模拟

通过对滑坡变形破坏现场的地质分析和相似材料模型模拟试验研究，对滑坡变形破坏及滑坡对管道的作用机理有了一定的认识。为了更深入地研究滑坡与管道相互作用的机理，验证前述认识，采用有限元计算方法对管道沿线典型滑坡进行了数值模拟研究。

研究方法

对于典型滑坡的滑体、滑动面、滑床、滑坡周界等可以用有限元方法进行离散，着重考虑滑坡潜在滑动面的摩擦接触情况，采用摩擦接触力学方法，对接触单元进行改进，建立接触力学模型，对滑坡进行研究分析。

在考虑可能滑坡面的法向抗张强度和切向抗剪强度下，通过多步迭代计算的方法，获得非线性接触问题的收敛解。通过接触力学计算分析，模拟滑坡的孕育；从滑动面的接触状态的改变（粘着、滑动、未接触状态），获得滑坡的滑动情况，滑坡的动态过程（启动、滑动、停止），同时还进行坡体和管道的应力、变形，以及边坡滑面接触关系分析。在这一过程中重点研究管道在滑坡发生过程中受力破坏的力学机理。

数值模拟

改进接触单元

数值模拟研究中，接触单元的法向与切向模量不一样的，法向模量EN可由两接触物间不可能相互挤入这一事实来确定。在三维情况下，如取N表示接触单元的法向、T1表示切向方向1、T2表示切向方向2，则可推导出：

（1）法向模量EN可由接触物间不可相互挤入的条件确定；

（2）而模量 、剪切模量 ；

（3）在切向平面内的剪切模量 可由接触状态（粘着、滑动等）确定。于是可确定接触单元的弹性矩阵[D]，进而可获得三维接触单元的刚度矩阵为：

山区油气管道地质灾害防治研究

式中：[B]是接触单元的几何矩阵。

对包含潜在滑面的整个边坡结构进行有限元离散，可得到边坡结构的刚度矩阵：

山区油气管道地质灾害防治研究

潜在滑面的强度

对于未出现任何破坏的边坡，其潜在滑面是具有一定强度的（实际上滑面还未形成。一般情况下是具有法向强度 和切向强度τ，的，一般情况下 。及 。法向强度的数值可取为边坡潜在滑面处的材料抗张强度，而对于切向强度通常在岩土力学中是取Mohr-Coulomb模型，认为切向强度τs与法向压应力之间σN是一个线性关系：

山区油气管道地质灾害防治研究

式中：c是边坡潜在滑面处的内聚力；ψ是内摩擦角。由上式可确定边坡潜在滑面处的切向强度。而在滑坡发生后，潜在滑面已经变为真实滑面，此时真实滑面位置退化为纯摩擦接触问题。

斜坡潜在滑面的接触状态

在滑坡发生前，潜在滑面处于完好的结合状态，在斜坡自重及外部条件的作用下，接触状态会改变。

非接触状态；当接触单元的法向应力是张应力，且该张应力大于潜在滑面处材料的抗张强度，即 ，此时，滑面将被拉开。接触单元的法向间隙会大于零，即Zg>0。接触单元的模量需取为一个很小的值EN=10-6～10-4Ee,Ee是接触单元的初始模量

山区油气管道地质灾害防治研究

当接触单元的法向应力是张应力，且当接触单元的切向应力τ大于潜在滑面材料的剪切强度τs时，潜在滑面将会破坏，此时Zg>0且EN=10-6～10-4Ee。

山区油气管道地质灾害防治研究

粘着状态：当接触单元的法向应力为压应力，且接触单元的剪应力小于潜在滑面处的抗剪强度时，此时Zg≤0。这意味着接触单元的法向隙Zg是一个很小的近似于零的值

山区油气管道地质灾害防治研究

当接触单元的法向应力为张应力，且接触单元的剪应力小于潜在滑面处的抗剪强度τs，此时

山区油气管道地质灾害防治研究

滑动状态：当接触单元的法向应力为压应力，且接触单元的剪应力大于潜在滑面处的抗剪强度时，此时潜在滑面将破坏，变为真实滑面，而Zg≤0。

山区油气管道地质灾害防治研究

由于剪应力的作用，潜在滑面变为真实滑面，滑坡发生了。接触单元的剪切模量 由下式决定：

山区油气管道地质灾害防治研究

非线性迭代过程的收敛计算

接触状态的最终确定，需满足以下条件：

非接触状态应满足Zg>0且EN=10-6～10-4Ee，而接触状态（粘着、滑动）应满足Zg=0（在指定精度ε下）。

每两次迭代循环计算时，接触单元模量的改变应满足

山区油气管道地质灾害防治研究

其中ΔE、ΔG是相对的模量改变量，ΔEp、ΔGp是许可的相对模量改变量，

山区油气管道地质灾害防治研究

对于滑动状态的接触单元，当下式满足时，滑动才处于稳定状态：

山区油气管道地质灾害防治研究

在以上3个条件同时满足时，该接触问题的计算才是收敛的。

滑坡过程的逐步计算

在边坡滑动过程的数值模拟中，使用了逐步计算的方法。用一系列时间步t1,t2，…，tK，…，tn来分解边坡的滑动过程。而对其中的某一个时间步tK，根据边坡下滑的加速度，对边坡的滑动部分施加相应的惯性力，所以，边坡将保持平衡状态。

在某一个时间步tK时，滑坡系统的势能改变ΔΠ可写为：

山区油气管道地质灾害防治研究

其中：

山区油气管道地质灾害防治研究

在tk时间步，上式中ΔWK是体积力bi和边界力Fs（不包含接触边界—滑面）所做的功，重力势能的改变量为ΔwE,ui为边坡系统的位移量，ΔUK是应变能改变量，σij和εij是应力和应变，Ω为边坡系统的总体积，ΔAK为滑动面上接触应力所做的功，p。和pt为接触面上的法向接触压力和切向力，r为接触面 和 分别为接触面的上、下边界， 和 分别为接触面上的上、下边界的切向位移，ΔKK为动能改变量。由于问题的复杂性，研究中未考虑材料破坏所消耗的能量。

在时间步从t1=0到tn计算完成后，可获得整个滑坡动态过程的应力、应变、接触状态改变情况、滑动距离等重要参数。