边坡稳定有限元分析

54边坡稳定有限元分析本例将演示如何使用有限元方法分析边坡稳定性并计算其安全系数。任务首先，分析无超载作用下的边坡稳定性，然后分析在大小为q=35.0kN/m2的条形超载作用下的边坡稳定性，最后为边坡施加预应力锚杆，并分析其稳定性。边坡的几何尺寸（包括各点的坐标）如下图所示。图25.1边坡几何尺寸（多段线上各点的坐标）土层剖面包含两种类型的土，其参数如下：表25.1岩土材料参数列表岩土材料参数/分类岩土材料1岩土材料2天然重度：γ[kN/m3]1820弹性模量：E[MPa]21300泊松比：υ[-]0.30.2粘聚力：ceff[kPa]9120内摩擦角：φeff[°]2338剪胀角：ψ[°]00饱和重度：γsat[kN/m3]202255计算我们使用“GEO5岩土工程有限元分析计算模块”（以下简称“有限元模块”）（v18版）来分析该问题。下面为建模和分析步骤：-建模阶段：分析设置和几何建模-工况阶段[1]：分析边坡无超载作用时的稳定性-工况阶段[2]：分析加入超载后边坡的稳定性-工况阶段[3]：分析加入锚杆后边坡的稳定性-结论建模阶段：分析设置和几何建模在分析设置界面中设置“分析类型”为“边坡稳定分析”，保持其他选项不变。图25.2【分析设置】界面注：选择“边坡稳定分析”时和选择“应力应变分析”时的设置以及建模过程几乎完全一样。在【分析】界面点击“开始分析”按钮即可以分析并计算边坡的安全系数。在“有限元-边坡稳定分析”模块中，各个工况阶段之间是相互独立的，即当前工况阶段的分析结果不受上一工况阶段分析结果的影响。下一步，设置全局坐标范围。设置的坐标范围要足够大，这样才能使得所要分析的区域不受边界条件的影响。对于该算例，设置全局坐标范围0m,40m，设置底边界距离多段线最低点距离为10m。设置各个多段线和土层剖面，其参数如下表所示。图25.3全局坐标对话框56表25.2各多段线及其节点的坐标列表多段线点编号多段线1多段线2x(m)z(m)x(m)z(m)10.00-9.000.00-11.50211.00-9.0014.00-11.00321.00-2.5021.00-9.25429.50-2.5040.00-9.00532.25-4.00640.00-4.00设置各个岩土材料的参数并将其指定到相应的分区。在本算例中，我们选择Drucker-Prager（DP）模型（见注）。设置两种岩土材料的剪胀角ψ均为0°，即当材料受到剪力作用时，其体积不发生改变。注：分析边坡稳定性时，必须选择非线性弹塑性模型作为岩土材料的本构模型，因为在边坡稳定分析过程中岩土材料会产生塑性应变，且塑性应变的产生是和岩土材料的强度参数c和φ相关的。在本算例中，我们选择Drucker-Prager作为本构模型，该模型和经典的莫尔-库伦模型相比，允许产生更多的塑性应变。在本章的最后，将给出不同本构模型下计算得到的安全系数的对比。图25.4添加岩土材料对话框57下图为指定岩土材料到相应分区之后的剖面。图25.5【指定材料】界面建模阶段的最后一步是生成网格。网格密度对边坡的稳定性分析影响很大，所以必须设置一个足够大的网格密度。对本算例，设置网格边长为1.5m，点击“启动网格生成”按钮生成网格。在本章最后，将给出采用网格边长为1.0、1.5、2.0m时计算得到的结果。图25.6【生成网格】界面–网格边长为1.5m工况阶段[1]：分析边坡无超载作用时的稳定性切换到工况阶段[1]，点击“开始分析”按钮进行分析，并保持分析设置类型为“标准”，即默认设置。58图25.7“分析设置”对话框注：边坡安全系数是通过对岩土材料强度参数c和φ进行折减得到的（强度折减法）。分析时，不断折减c和φ直到边坡达到临界状态（破坏），此时c和φ的折减系数即为安全系数（详细信息见帮助文档-F1）。因此，程序通过下式得到边坡的安全系数：FS=tanφs/tanφp其中：φs–真实内摩擦角φp–边坡失稳时的内摩擦角在边坡稳定性分析中，非常重要的两个计算结果为位移矢量图和等效塑性应变εeq.,pl图。等效塑性应变区域即为边坡的潜在破坏区域（见下图）。图25.8【分析】界面–工况阶段[1]（等效塑性应变εeq.,pl）59注：“有限元-边坡稳定分析”模块的分析结果中包括位移等值图（Z向和X向）和应变等值图（总应变、等效塑性应变）。软件得到的计算结果是折减参数之后的情况，因此该结果和实际的情况是不同，该结果表示的是边坡失稳时的情况。工况阶段[2]：分析加入超载后边坡的稳定性添加工况阶段[2]，在【超载】界面中定义如下图所示类型的超载。图25.9“添加超载”对话框–工况阶段[2]点击“开始分析”按钮得到分析结果，观察等效塑性应变等值图。图25.10【分析】界面–工况阶段[2]（等效塑性应变εeq.,pl）工况阶段[3]：分析加入锚杆后边坡的稳定性添加工况阶段[3]，在【锚杆】界面中点击“添加”按钮，弹出添加锚杆对话框，输入如下所示的锚杆参数，锚杆预应力为F=50kN。-锚杆长度：l=16m60-倾角：α=17°-锚杆水平间距：b=1m-锚杆直径：d=20mm图25.11添加锚杆对话框–工况阶段[3]注：在进行边坡稳定性分析时，预应力锚杆中的预应力采用作用在锚头处的等效集中压力代替。集中压力作用点处的岩土材料可能会发生塑性变形。由于等效塑性应变的位置即为潜在滑动面的位置，因此，用户应当对得到的塑性应变分布进行认真分析。如果得到的锚头附近塑性应变对边坡稳定性起着决定性作用，那么这时用户必须对初始输入参数进行一些修改。其他参数保持不变，点击“开始分析”按钮分析得到结果，并观察得到的等效塑性应变等值图（和前一工况阶段类似）。图25.12【分析】界面–工况阶段[3]（等效塑性应变εeq.,pl）61现在，将分析得到的安全系数记录到一个表格中，然后采用其他的本构模型（莫尔-库伦模型（MC），修正莫尔-库伦模型（MCM））再次分析该问题。注：检验计算得到的滑面形状（等效塑性应变带）是非常重要的，因为，在某些情况下，局部的失稳可能发生在我们认为不会发生失稳的区域（详细信息见帮助文档-F1）。当采用DP模型并设置网格边长为1m时，将会得到下图中的结果。从图中可以看出，锚头附近的土体出现了局部破坏，但是边坡本身却没有发生破坏，这时得到的安全系数并不是边坡的安全系数。如果这种情况发生，可以对于模型的参数进行一定的修改，例如：-增大网格边长；-将锚头附近岩土材料的强度参数提高（具有更高的c和φ），或采用弹性模型；-在锚杆锚头下设置一个梁单元（这样可以让预应力能均匀的传递到岩土材料中）。图25.13【分析】界面–工况阶段[3]（锚头附近产生局部塑性应变，DP模型，网格边长1m）图25.14【分析】界面–工况阶段[3]（使用极限平衡法中的Bishop法分析得到的圆弧滑面）62结论下表列出的是各个工况阶段采用不同的非线性本构模型和网格密度得到的边坡安全系数。为了对比，我们还列出了“GEO5-土质边坡稳定分析”模块中采用Bishop法和Spencer法得到的安全系数。表25.3使用不同本构模型和方法得到的边坡安全系数本构模型网格边长[m]工况阶段[1]FS工况阶段[2]FS工况阶段[3]FS注DP1.01.671.441.03锚头附近区域发生局部塑性应变DP1.51.691.421.67DP2.01.741.481.69MC1.01.561.350.90锚头附近区域发生局部塑性应变MC1.51.581.351.56MC2.01.601.411.56MCM1.01.781.561.14锚头附近区域发生局部塑性应变MCM1.51,811.541.78MCM2.01.851.601.81BISHOP极限平衡法－1.511.331.47见下文SPENCER极限平衡法－1.511.321.52见下文注：采用极限平衡法时，在“土质边坡稳定分析”模块中选择“标准-安全系数法”作为分析设置，并分别使用Bishop法和Spencer法搜索得到最危险圆弧滑面（不限制搜索区域）。从上面的分析结果可以得出如下结论：-对某些重点区域的网格进行局部加密会得到更为精确的结果，但是另一方面，会使得各个工况阶段的分析时间增加。-边坡稳定分析中必须使用允许发生塑性变形的非线性弹塑性模型。-最大等效塑性应变εeq.,pl发生的区域为潜在破坏区域。-DP本构模型得到的安全系数比MC模型得到的稍大，即DP模型允许更大的塑性应变。