【仿真技术】简单的极限平衡计算的编程实现

笔者从教时间不长，只有三四年，每一次教学都能从中获取很多新的发现。正如徐芝纶老先生所言“学无止境、教亦无止境”。

极限平衡是土力学中计算边坡安全系数最常用的方法，笔者在教学过程中给学生编写了一个简单的极限平衡计算程序，程序非常简单，只有几十行。后面有时间的话会加入图形界面进行进一步完善。

计算采用最经典的瑞典条方法，具体原理公式如下：

图1 瑞典条方法计算基本原理

如下图所示有一个均值边坡，尺寸如下图。假定土体的重度为23.5kN，土体的黏聚力25kPa，内摩擦角30°，假定圆弧的圆心坐标为（10，40），半径为30m，试计算该滑弧对应的安全系数。

本程序在几何图形计算中使用了shapely库，shapely是一个python包，用于设置平面特征的理论分析和操作（通过python的 ctypes 模块）来自著名和广泛部署的地理类库的功能。本程序可以为三个部分：

1. 计算滑坡体的几何信息：采用圆心和半径绘制圆，然后将这个圆与边坡多边形相切，交集部分即为滑坡体，代码如下：

polygon = Polygon([[0, 0], [50, 0], [50, 40], [40, 40], [10,10],[0,10], [0, 0]])circle = Point(center[0], center[1]).buffer(radius)pintersect = polygon.intersection(circle)

2. 将滑坡体分条：我们提供的思路是将滑坡体x方向的最大和最小值找出来，然后绘制一系列不同x坐标的竖线，通过计算竖线与滑坡体的交点，得到每一个滑块的信息。具体代码如下：

''' get the x coordinates limit the pintersect '''xlimit = pintersect.bounds[0], pintersect.bounds[2]ylimit = pintersect.bounds[1], pintersect.bounds[3]'''slice the pinterset into 10 slices '''slice_num = 20xstep = (xlimit[1] - xlimit[0])/slice_num''' get the slices using intersction'''slice_points = []for i in range(slice_num+1):    line_tmp = LineString([(xlimit[0]+xstep*i, ylimit[0]), (xlimit[0]+xstep*i, ylimit[1])])    pts = pintersect.intersection(line_tmp)    slice_points.append(pts.xy)

3. 计算该滑面对应的安全系数：这一块稍微麻烦，需要把每一个滑块上的物理量算出来，然后用公式（1）计算边坡的安全系数，这一块代码略长一点，有兴趣的同学可以自己看代码。

最终实现的效果如下所示，边坡的安全系数Fos = 1.13。其实这个程序的计算代码并不长，但是把计算过程绘制出来还是毕竟麻烦的，采用了Matplotlib图形库，有心的朋友可以仔细看看代码呀。

图2 程序初步运行结果