（岩土工程专业论文）重力坝深层抗滑稳定研究.pdf

摘要 重力坝抗滑稳定是重力坝设计的一个经典领域。本研究将具有坚实物理理论 背景的、基于非相关联流动法则的三维边坡稳定机动分析方法和可靠度理论应用 到重力坝深层抗滑稳定分析中，使确定性分析方法和可靠度分析的成果互为补 充，从而科学地解决规范中遇到的问题。 本研究报告继续完善基于非相关联流动法则的三维机动分析方法。通过与已 有算例进行比较，对该方法在三维边坡稳定分析和重力坝深层抗滑稳定领域应用 的可行性进行验证。算例表明，基于非相关联流动法则的三维机动分析方法避免 了上限法中由于摩擦角过大而导致速度场不收敛的局限性。 然后对岩土工程中的可靠度理论和数值方法进行研究分析，提出适用于三维 边坡稳定分析的数值方法。通过在三维边坡稳定分析中的应用实例，验证了该方 法在边坡稳定分析中的实用性和可靠性。 最后，将三维机动分析法和可靠度分析应用到向家坝深层抗滑稳定分析研究 中，通过最后求解得到的可靠度指标作为设计可行性依据，为三维边坡稳定方法 在重力坝深层抗滑稳定中的应用提供了一种新思路。采用这一新的方法，可以较 好地解决三维安全系数在现有的重力坝设计规范中无法参照的问题。 稳定 关键词：三维边坡稳定分析，非相关联流动法则，可靠度理论，重力坝抗滑 第1 章前言 第1 章前言 1 1重力坝深层抗滑稳定分析有关问题的提出 重力坝抗滑稳定是重力坝设计的一个经典领域。其相应的分析方法是建立在 极限平衡与极限分析方法基础上的。极限平衡与极限分析法已经广泛应用到二维 边坡稳定分析中，而在三维边坡稳定方面的工程应用还比较少，尤其是在重力坝 深层抗滑稳定领域的应用。 使用极限平衡方法，相应“抗剪断”的强度指标，要求得到一个超过3 的安 全系数。半个世纪以来，在应用这一设计准则进行深层抗滑稳定分析时，频繁地 遇到复核的成果无法过关的困境。然而，这些工程完工后却绝大多数在正常工作。 因而，也迫使人们重新思考这一方法包含的内在缺陷。 同时由于岩土工程问题存在着天然的不确定性，处理不确定性也是岩土工程 师一项重要的工作内容。在进行工程设计和安全评价时，不仅要很好地了解已经 掌握的各种分析、判断手段，而且更需要把握这些分析过程中所包含的各项不确 定因素及其影响。工程建设中的重大决策实际上是对各项不确定因素造成的风险 的评价。在结构工程领域，风险分析和可靠度理论已经形成了一个成熟的手段， 将这一理论引入水工结构的分析，有利于科学地解决规范修正中遇到各个问题。 确定性方法和可靠度分析的成果能够互为补充，从而提高分析成果的精度。 本研究希望通过以下途径来解决这些问题。 ( 1 ) 引入理论体系严格和使用范围广泛的极限分析方法 ( 2 ) 在安全控制标准和强度指标方面引入风险分析的原理，但是并不抛 弃和否定建立在安全系数基础上传统方法 以下将对边坡稳定的三维分析方法和可靠度方面的研究现状进行综述。 1 2边坡稳定三维分析方法的研究现状 自从d u n c 锄( 1 9 9 6 ) 归纳总结了三维边坡稳定分析的主要研究成果后，基于条 柱法的三维分析方法取得了很大进展。s e e d 等( 1 9 9 0 ) 和m o r g 吼s t e m ( 1 9 9 2 ) 讨论了 发展三维边坡稳定分析在实际中的必要性，认为应该把它作为主要研究方向。三 l 中国水利水电科学研究院博士学位论文 维边坡稳定分析可以更加真实地反映边坡的实际形态，比如由于岩体间断的存在 会对安全系数产生影响。特别是当滑裂面已经确定时，使用三维分析可以恰当地 考虑滑体内由于滑裂面的空间变异特征对边坡稳定安全系数的影响，使结果更具 有实际价值。 1 2 1 三维边坡稳定的极限平衡法 我们知道，为了使问题静定可解，任一个二维极限平衡方法都引入了一些关 于界面内力的假定。为了使数值计算可行，通常进一步采取一些简化方法使静力 和力矩平衡条件不完全满足，并且简化滑面为圆形。在将这些极限平衡方法扩展 到三维情况时，所有这些局限将遇到更加巨大的挑战。这些方法区别在于：( 鸪关 于条块内力的假定；( b ) 平衡方程；( c ) 滑面形状的简化。 三维极限平衡法首先将滑动土体剖分为具有垂直界面的条柱，并建立如图 1 1 所示的坐标系，其中x 和y 的正方向分别与滑坡方向和重力方向相反，z 轴的 正方向按右手法则确定。中间沿滑动方向尺寸最大的断面称之为中性面，如果滑 裂面对称，中性面就是对称轴所在的平面。如果不作任何静力假定，每个条柱上 的作用力应当如图1 2 所示，其中平行于y o z 平面的界面称之为行界面，即图1 2 中的a b f e 和d c g h ；平行于x o y 平面的界面称之为列界面，即图1 2 中的a d h e 和b c g f 。 图1 1具有垂直界面条柱的滑动土体 图1 2未作任何假定的条柱上的作 用力 第】章前言 1 简化b i s h o p 法 i h u n g r ( 1 9 8 9 ) 本方法主要包括以下假定( 参见图1 1 和图1 2 ) ： 1 ) 忽略条柱间沿y 轴方向的所有剪力g 和k 2 ) 滑裂面上的剪切力平行于x o y 平面，即忽略平行于y o z 平面的剪切力矗：。 本方法满足每个条柱沿y 轴方向力的平衡，并根据绕z 轴的整体力矩平衡求 解安全系数。但由于未满足沿x 轴和z 轴方向的整体力的平衡，所以根据整体力 矩平衡条件求解安全系数时与坐标轴的位置有关，故本方法比较适合于滑裂面为 旋转面的情况。 i i h u a n g & t s a i ( 2 0 0 0 ) h u a n g & t s a i 首先根据 m o h r - - c o u l o m b 抗剪强度准则，定义了每 个条柱的安全系数咒；，如图1 3 所示。 抗滑力乃由两部分组成，即平行于 x o y 平面的，和平行于y o z 平面的瓢， 分别定义沿x 轴的安全系数最。和沿z 轴 的安全系数尼对于每个条柱都是相等 的，但每个条柱的足f 随岛和啦的变化而 不同。 留1 3 条柱底部抗滑力示意圈 本方法忽略了条柱间沿y 轴方向的所有剪力，只考虑条柱沿y 轴方向的静力 平衡，故条柱间的其他作用力e 只eq 均未出现在建立的方程中。 本方法特点是每个条柱的滑动方向也作为求解的一部分。但本方法由于未满 足沿工轴和z 轴方向的整体力的平衡，所以在建立力矩平衡方程时就与所绕坐标 轴的位置有关，这也就要求滑裂面至少有一部分为球形。 2 简化j a n b u 法 i h u n g r ( 1 9 8 9 ) 本方法的假定同b i s h o p 法，满足每个条柱沿y 轴方向力的平衡，并根据沿z 轴方向的整体力的平衡求解安全系数。但由于未考虑滑裂面上平行于j 脚平面的 剪切力瓦：和沿z 轴方向整体力的平衡，所以比较适合于滑裂面对称的情况；同 时h u n g r 还指出此方法由于忽略了所有条间剪力，所以不适合计算条间的抗剪强 ! 望查型变皇型兰旦窒堕塑主兰竺堡壅 度比较大而底滑面的抗剪强度相对较小的情况。另外，h u n g r 发现，此方法通常 给出比较小的安全系数。 i i 冯树仁等( 1 9 9 9 ) 其假定主要包括( 参见图1 2 ) ： 1 ) 忽略条柱问沿y 轴方向的所有剪力g 和k 2 ) 忽略平行于y o z 平面的剪切力瓦；。 本方法满足每个条柱沿y 轴方向力的平衡，并根据沿x 轴方向的整体力的平 衡求解安全系数，比较适合滑裂面对称的情况。冯树仁等也指出若滑面为球面， 旋转椭圆面等规则光滑曲面，上述假定比较真实地反映了实际情况。若剪切面形 状不规则，上述假定可能会导致一定的误差。 3 三维s p e n c e r 法 i z h a n gx i n g ( 1 9 8 8 ) 本方法可看作二维s p e n c e r 法在三维条件下的扩展。其假定主要包括( 参见 图1 2 ) ： ( 1 ) 滑裂面为对称的椭球面，滑裂面在x o y 平面内为圆弧，在z 轴方向由椭 圆面生成。 ( 2 ) 滑裂面上的剪切力平行于x o y 平面，即忽略平行于y o z 平面的剪切力瓦：。 ( 3 ) 将所有的条间力简化为一个平行于x o y 平面的作用力震，其与工轴的倾 角为常量，这假定相当于二维领域中的s p e n c e r 法。 ( 4 ) 假定条柱底部作用一个端部力只并假定其大小和作用方向。 根据上述假定，安全系数可由求解各条柱在x 和y 方向力的平衡方程以及滑 块绕z 轴的弯矩平衡方程而获得。 本方法通过假定的条柱底部端部力p 反映三维效应，显得比较简单，同时也 使得本方法的理论基础不够严格。另外，本方法由于假定滑裂面对称，故其很难 在实际工程中应用。 i i c h e n & c h a m e a u ( 1 9 8 2 ) 其假定主要包括( 参见图1 2 ) ： 1 ) 滑裂面对称； 2 1 滑动方向只沿着x o y 平面，则滑裂面上的剪切力平行于x o y 平面，即忽略 4 第1 章前占 平行于y o z 平面的剪切力瓦：； 3 ) 忽略条柱问行界面上沿z 轴方向的剪切力日； 4 ) 假定条柱间行界面上的作用力和g 的合力平行于x o y 平面，其与x 轴 的倾角为常量，这一假定相当于二维领域中的s p e n c e r 法： 5 ) 假定条柱间列界面上的剪切力p 和y 的合力平行于条柱底部，并假定其 分布。 本方法通过在工轴和y 轴方向的整体力的平衡和绕z 轴的整体力矩平衡求解 安全系数。由于假定滑裂面对称，故其很难在实际工程中应用。另外，在计算非 粘性土坡时，某些情况下，本方法得出的三维安全系数比二维的安全系数还低， 因此本方法的可靠性有一些问题( h u t c h i s o n & s a r m a ( 1 9 8 5 ) ) 。 i i i 陈祖煜等( 2 0 0 1 ) 本方法首先建立如图1 1 所示的坐标系，图1 4 所示为垂直条柱上的受力情 况。主要引入了以下几点假设： ( 1 ) 作用在行界面( 平行于y o z 平面的界面，图1 4 中的a b f e 和d c g h ) 的条间力g 平行于x o y 平面，其与x 轴的倾角肋常量，这一假定相当于二维领 域中的s p e n c e r 法。 ( 2 ) 作用在列界面( 平行于 x o y 平面的界面，图1 4 中的a d h e 和b c g f ) 的作用力q 为水平方 向，与z 轴平行。 ( 3 ) 作用在底滑面的剪切力r 与x o y 平面的夹角为p 。规定剪切 力的z 轴分量为正时砌正值。 假定同一列条柱( z = 常量) 的 p 值相同，对不同z 坐标的条柱， 假定p 的一个分布形状。 1 ) p = r = 常量； 图l 4 作用在具有垂直界面的条柱上的力 2 ) 在x o y 平面的左、右两侧假定p 的方向相反，并线性分布，假定此分布形 状为缃，则有 5 中国水利水电科学研究院博士学位论文 , 胪0 r = k z - r 。：，。z 三? 0 ( 1 ) 反2 z o 正方向i = l ，广，负方向f _ 一1 肌。从现在开始，所有的推导都针对o z 正方向的条块进行。这些推导同样适用于o z 负方向的条块。中性面上的条块沿 o x 方向从1 到肛编排。某一o z 方向第i 列及删方向第行的条块编号为u 。 1 三维速度场的求解 对编号为u 的条块，它的速度为。该条块相对i - l d 条块的速度，即列 与列界面上的相对速度记为k l ，。与此类似，k ，代表条块i j 相对条块驴l 的 速度，如图2 7 所示。 c 篇， 条块 ( 。j ) 条块 ( i j - 1 ) i ”【w 图2 7 依据圪l ，和l 确定， 根据流动法则和位移协调的要求，某一条块的速度可由其相邻条块的速度 来确定。具体的求解方法可根据条柱所处的不同位置分为三类： ( 1 ) 中性面条块的速度场 中性面各条块( 声0 ) 的速度场的求解可视为一个二维问题，根据式( 2 9 ) 与式 ( 2 1 0 ) 可以很容易地求出。，与p o ，。 ( 2 ) 破坏体内条柱的速度场 对于编号为( i ，力的条块，其速度可由与该条块相邻的左侧与下侧条块的 速度l ，和p 抽求出，如图2 7 所示。 根据m o h r - c o u l o m b 屈服准则与相关联流动法则，有： m ( j ，f j ) = s i n c e , ， ( 2 1 3 ) ( ，一l ) = s i n # , 3 j ( 2 1 4 ) m ( w ，n ，。，) = s i n 谚。 ( 2 - 1 5 ) 第2 章三维边坡稳定机动分析方法 速度协调条件要求 k + ，= k ，一，川 ( 2 1 6 ) = k , j - v ，- l ( 2 1 7 ) 当v , j j 、1 已知时，联立式( 2 1 3 ) - - ( 2 1 5 ) ，即可采用迭代法求出v , j 在墨弘 z 轴上的三个分量。 ( 3 ) 滑坡体边界条块的速度 对于每列的第一个条块( 产= 1 ) ，如图2 7 i f ) 所示，由于不存在与其相邻的行界 面条块，此时式( 2 1 5 ) 不存在，因此方程的个数比未知数的个数少1 个，此时需 要引入一个假定才能求解该条块的速度。 假定第一个条块速度的绝对值与相邻的产1 条块的速率存在如下关系。 i 巧，。i - - e , i v 吐一 ( 2 1 8 ) 有了这个附加假设，即可联立式( 2 1 3 ) ，( 2 1 4 ) ，( 2 1 8 ) 确定每一列边界条块的 速度，作为一种简化计算方法，通常取卣= 手，= 六一一1 0 。 2 安全系数的求解 整个滑坡体的功能方程可以表示为： 。+ 磁，+ 广胛+ r 旷 ( 2 19 ) 式中矿为以上求得的应变速率。方程左边三项分别表示作用在条柱行界面、 歹0 界面和底滑面的内能耗散。这里所有未知的应力砖全部被消去，这样就可以 获得包含于式( 2 1 9 ) 中的唯一的未知变量一安全系数f 或极限承载力兀 如上所述，在边坡稳定的上限分析方法中，包括了一个重要的假定，即岩 土材料在破坏时，其变形速率与滑面的夹角为，即服从相关联流动法则。然而， 大量的试验表明，如果摩擦角过大，式( 2 3 ) 有可能过高地估计了土体发生塑性流 动时的剪胀，违背了物理实际情况，从而有可能过高地估计了边坡破坏时的极限 荷载或安全系数。 因此希望可以建立一个与底滑面成任意指定剪胀角的速度场玎满足真实 的物理背景。但是这样就引起塑性速度v 与“组合摩擦力”，不再保持垂直，组 2 l 中国水利水电科学研究院博士后研究报告 合摩擦力p 在塑性速度v 上做功不再为零，使得( 2 ，1 9 ) 式的求解变得困难。下 面我们希望通过引入第二个速度场( 虚拟速度场) 来解决该问题。 2 2虚拟速度场的概念 陈祖煜等( 2 0 0 5 ) 利用图2 8 说明了极限分析理论中假定分析单元的主应力 轴与主应变率轴重合的假定。“两轴重合”假定意味着在破坏面上，法向应力、 剪应力与剪应变在同一平面内。由于极限分析是从建立一个协调的速度场开始 的。一旦假定主应力轴和主应变率轴重合，就可以由应变率轴的方向来确定应力 的方向。我们知道，在滑动面上的剪应力本来是包含了其大小与方向两个未知分 量的。一旦应变速率场确定，其方向分量变为已知，则问题可由静不定变为静定 可解。 剪 在此方向上 力与剪应变 a 图2 8 说明“重合”假定示意图( 陈祖煜等，2 0 0 5 ) 在上述推导内力功的过程中，我们假定法向力剪切力c ，组合摩擦力p 和塑性速度y 都在同一个平面a 内，如图2 9 所示。当塑性速度y 与“组合摩擦 力”p 不垂直时，我们希望引入一个与底滑面成咖的虚拟速度场q ，然后对滑坡 体建立如下的虚功方程： p o ，+ 璐j + d ：j = w q + r q ( 2 2 0 ) 第2 章三维边坡稳定机动分析方法 p l a n ea va n d n 图2 9n c ，p , v 和q 构成的平面a t h es l i ps u r f a c e y 和q 这两个速度场使得基于非相关联流动法则的边坡稳定分析方法得以 实现： ( 1 ) 破坏面上的塑性速度v 可以与滑面成任意剪胀角奶即速度法向与切向 满足 善：譬丝：一t a n f ， ( 2 2 1 ) 圪汐f d 7、7 ( 2 ) 一旦塑性速度v 确定，就可以得到在平面a 内，与破坏面成磊角的虚拟 速度场q ，使得其与组合摩擦力p 垂直。 ( 3 ) 建立如式( 2 2 0 ) 的虚功方程得到安全系数。 这样，该方法通过建立满足非相关联流动法则的速度场矿和满足 m o h r - c o u l o m b 屈服准则的虚拟速度场q 得以实现。下面我们希望通过楔体稳定 这一特殊的三维问题，继续对虚拟速度场的概念进行说明论证。 2 3楔体稳定问题的广义三维极限分析方法 c h e r t ( 2 0 0 4 ) 提出了计算剪胀角为任意值的楔体稳定广义解，并且使用严格的 数学力学方法证明了当楔体以摩擦角发生剪胀时，安全系数取最大值，从而也实 现了潘家铮最大原理在楔体稳定领域的理论证明。实现这一突破的关键技术就是 引入了两个速度场的概念。这里简要回顾一下广义解的推导。 中国水利水电科学研究院博士后研究报告 1 对如图2 1 0 所示的楔体建立一个分别按与两滑面的夹角为p ，和n 进行剪 胀的速度场，此时楔体两滑面的“组合摩擦力”的方向矢量p f 与肼的计算公式 分别为： p l ：一些肼+ c o s ( p , - # , i ) 砷 ( 2 2 2 ) c o s p tc o s p l p ，：一兰堕珊+ 竺幽珥 c o s p rc o s p ， 式中：m 为楔体的塑性速度，玛，刀，分别表示两滑面的法向量。 图2 1 0 楔体破坏时的速度场 这里，速度场m 是楔体发生塑性位移的真实速度场，也就是应用非相关联 流动法则时材料发生的真实位移。在“两轴重合”的假定下，我们通过p ，和办 获得滑面上的剪力方向，使问题变为静定可解。 2 将作用在楔体上的力投影到与左、右两个“组合摩擦力”相垂直的轴上， 该轴的表达式为： q ：a 。n ：一皇坐红里堕垒型m x t t r 一! ! 亟旦二盘韭堕监啊。埘 竺8 易：譬，。、 5 向8 b ( 2 2 4 ) + 竺垫 垃竺鲤立打，。打， 。 此时楔体的静力平衡方程的表达式为： qr 巳+ q c e r + q w = 0( 2 2 5 ) 式中，矽为楔体所受的外力矢量。 对于式( 2 2 s ) ，可以将q 理解为第二个速度场，它与左右两渭面的组合摩擦 2 4 第2 章三维边坡稳定机动分析方法 力p ，和肼垂直，它与速度场m 组成楔体的两个速度场，如果令作用于楔体上的 所有的力在这一虚速度场q 上进行如式( 2 2 5 ) 所示的虚功原理操作，那么，作用 于楔体上的两个未知内力，即左、右两个“组合摩擦力”所做的功为0 ，这样就 获得仅包含安全系数一个未知数的控制方程式( 2 2 6 ) 。 t a ( p , ，岛，d = 一c o s 九c o s ( p , 一丸) 白4 一c o s 屯c o s ( p t 一九) c 盯4 + 【s i n 九c o s ( , , 一九) 丸+ e o s ( p ，九) s i n 丸 + 口，c o s ( 岛一九) c o s ( 屏一丸) 】矿2 0 ( 2 2 6 ) 式中：q ，虬，q 为楔体所受的外力分别在m 、吩与一上的投影，即： 形= a w m + 6 1 ，啊+ c w n ， ( 2 2 7 ) 这样通过建立两个“速度场”的方法成功地获得了采用非相关联流动法则的 楔体稳定分析的广义解。下面将把该方法推广到建立一般非相关联流动法则的三 维边坡稳定机动分析方法上。 2 4非相关联流动法则的三维边坡稳定机动分析方法 离散模式与三维边坡稳定分析上限方法相同，参考图2 5 ，2 6 ，这里同样引 入了中性面的假定。离散条柱速度场的求解主要包括以下两个步骤： 第一步，建立一个满足非相关联流动法则的速度场k 该速度应满足与条柱 底滑面的夹角为，。，在行界面与列界面的相对速度与侧面的夹角应分别为 ，与”l ，。这是一个真实的速度场，一旦建立了这一真实的速度场，根据主 应力轴与主应变率轴重合的假定，可以确定条柱底滑面与界面上“组合摩擦力” 的方向。 r 第二步，建立一个虚拟的速度场q 。由于这一速度场与底滑面与界面上的“组 合摩擦力”的方向垂直，条柱底滑面与界面上的内力中的“组合摩擦力”在这一 速度场上做功为0 ，因而可通过功能平衡方程式非常简单地求出安全系数。 2 4 1 条柱真实速度场y 的求解 1 中性面速度场的确定 这里中性面各条块的速度场的求解不再视为一个二维问题，而将中性面视为 2 5 中国水利水电科学研究院博士后研究报告 一系列具有一定厚度的柱体，作为一个三维问题进行求解，如图2 1 1 所示。 对于中性面的第一个条块，其速度场与左右滑面的关系如图2 1 1 ( ”所示，由 此可得计算公式为： ( a ) 中性面条柱与左右底滑面示意图 ( b ) 中性面第一条块速度场的求解 ( c ) 中性面其它条块速度场的求解 图2 1 1 中性面速度场求解示意图 i o ( p - ) = s i n 虻u 。 由( j ，l ，1 ) = s i n l 。 ( 2 2 8 ) 【i = l 式中，n i l ，l ，1 分别为编号为( l ，1 ) 与( 1 ，1 ) 的条柱底滑面的法向量，“j v o ，l f ”表示 矢量v o ，1 的模，p “。表示条柱( f ，d 底滑面经缩减后的剪胀角。 对于中性面的中间条块，其速度场与左右滑面以及前一条柱速度场的关系如 图2 1 l ( c ) 所示，计算公式为： 第2 章三维边坡稳定机动分析方法 jo ( - ) = s i n 屯。 m ( ，u ) 2s i nu ， ( 2 2 9 ) 【西( p 以，“+ j ) 2 s i n j ，。 式中，m 1 1 o 。，分别表示编号为( 1 ，力与( 1 ，d 条柱的底滑面以及条柱界面 的法向量，+ 加表示条柱( o ，力行界面经缩减后的剪胀角。根据中性面上各点的 速度均在该面内的假定，此时中性面不再一定是一个竖直面，其法向量可由中性 面的底滑面法向向量与其塑性速度确定。 2 条柱速度场的确定 对于编号为以力的条块，其塑性速度y 可由与相邻编号为( i q ，力与( i , j - 1 ) 条柱 的速度确定，如图2 7 所示。即： o ( 巧，i v , j ) = s i n u , 。 ( 2 3 0 ) 似， l l ) = s i n g , , t ” ( 2 3 1 ) m ( ，1 w ) = s i n ”叫， 式中，k 分别为条柱( f 力与条柱( f - 1 ，力及( f ，产1 ) 的相对速度，在已知吃。及 。的条件下，联立方程式( 2 3 0 ) - - ( 2 3 2 ) ，即可求出匕在与弘z 轴方向的分量。 对于每列的第一个条块，由于无相邻行界面条柱，根据喜假定，式( 2 3 2 ) n - f f t c a 下式代替。 叫= 考瞻一 ( 2 3 3 ) 2 4 2 虚拟速度场q 的求解 为了得到虚拟速度场，首先要确定条柱底滑面、行界面与列界面“组合摩擦 力”p 的方向。 1 滑面上的组合摩擦力方向 若条柱滑面的真实塑性速度方向用单位矢量肼表示，现将滑面上的内力分 为两部分，第一部分为凝聚力，用c 表示，其大小为e a ，第二部分为法向力 中国水利水电科学研究院博士后研究报告 与摩擦力t a i l 晚的合力，称为“组合摩擦力”，用户表示。 则滑面上的组合摩擦力单位矢量可以表示为m 与法向量露的线性组合，如 图2 1 2 所示，有： p = l a b i - i b q m =! 型垒二盟栉一些! 堕m c o s p ec o s _ p 。 ( 2 3 4 ) 图2 1 2 滑面上的“组合摩擦力”、法向力与塑性速度示意图 2 虚拟速度场q 的求解 ( 1 ) 中性面虚拟速度场的求解 这里中性面虚拟速度场的求解也作为一个三维问题进行求解，如图2 1 1 所示。 其虚拟速度场为 幺j = 正ux p l ( 2 3 5 ) 这里正u ，p l 分别表示条柱( 一1 ，力和( 1 ，力的底滑面组合摩擦力的单位方向矢量。 ( 2 ) 条柱虚拟速度场的求解 对编号为( f ，力构造一个机动许可的虚拟速度场q j 其协调方程满足： q ，j p j j = 0 q “只l = 0 ( 2 3 6 ) q m p = 0 特别地，对于每列的第一个条柱，其速度协调方程要求 第2 章三维边坡稳定机动分析方法 q 。p = 0 q 2 只。，= 0 ( 2 3 7 ) 陋= 毒jq ，一 式中6 为大于0 的系数，通常取6 = 1 。见。，p 由，a 。分别表示条柱 d 的底滑面、 行界面与列界面组合摩擦力的单位方向矢量，( ，j 表示条柱编号p 1 ，力的塑性速 度。各条柱的二次虚拟速度场可用显式求出。 2 4 3 安全系数的求解 根据功能平衡方程式( 2 2 0 ) 计算边坡的安全系数。条柱底滑面的内能耗散的 表达式为： 口= 巳2 + q 加q ，= q 加4 | q ，i c o s 啦妇 ( 2 3 8 ) 式中，j 0 为条柱( i ，j ) 底滑面上的组合摩擦力，c 加为底滑面上的凝聚阻力，仑 为条柱的虚拟速度场，4 。为底滑面的面积，q 加为底滑面上的抗剪力与其 塑性速度q 。的夹角。 同理，条柱列界面与行界面的内能耗散为： 式中：只“，霉。，分别表示条柱( f ，力行界面与列界面上的组合摩擦力，q 。g 卅， 分别表示条柱以力行界面与列界面上缩减后的凝聚阻力，其大小为c a ，珐+ 分别表示这两个面上的虚拟速度场，如j ，4 。分别表示这两个面的面积， n h 卅，分别为这两个面上的剪切力r 与各晃面虚拟速度q 的夹角。 2 5确定临界滑面的最优化方法 与二维领域相比，三维问题的目标函数具有多自由度、多极值、高度非线性 的特点。对于此类复杂问题，寻求安全系数的整体极值将变得更加困难。 最优化问题的一般提法是：对于一个具有以个自变量z r 吖z j ，幼，纠的目标 2 9 功q ， 魄 舾 宝致钆a =，驴以氐酚既 中国水利水电科学研究院博士后研究报告 函数f ，确定使f 获得最小值的自变量z 。在边坡稳定分析中，临界滑裂面 的搜索问题实质上为求解目标函数( 安全系数或加载系数) 的极值问题。因此，首 先需要将安全系数f 表述为一系列待优化的自变量的函数。 2 5 1 任意形状三维滑裂面的构筑 c h e ne ta 1 ( 2 0 0 0 ) 提出下面具有一般意义的任意形状三维滑裂面的构成方法。 首先，在中性面，即0 断面上建立一个任意形状的二维滑裂面，中性面的滑 裂面由一系列的控制点由直线或光滑曲线连结，各控制点之间的节点的坐标通过 内插确定，以进一步细分滑裂面。然后，根据中性面的数据在z 轴方向用椭球定 义滑裂面的空间形状，如图2 1 3 所示。滑体由该断面上一系列倾斜界面来剖分， 倾斜界面在该断面的高度用h 表示，其倾斜角碇义为与y 轴的夹角。界面在边 坡和底滑面上的点的坐标分别表示为( ，圯) 和( 矸，玎) 。每个界面在z 方向伸展形 成垂直x - y 平面的行界面，如图2 1 3c o ) 中断面即4 卅惭面所示。相应某一z 坐标 值的界面，其h 值由下式给定： 矿h 2 + 丽z 2 = l ( 2 4 0 ) 日2 k ( 日+ b ) 1 2 ”7 式中两个参数f 和口可以描述一般性的滑裂面( 将随后讨论) 。该断面滑裂面上的 点的坐标( 薪，j 彳) 是根据中性轴上的坡顶坐标( ，躬) 和通过上式得到的h 确定的 ( 图2 1 4 ) ： ：一一h s i n 5 l j 一h c o s g j ( o ) ( 2 4 1 ) 第2 章三维边坡稳定机动分析方法 ( b ) 图2 1 3任意形状三维滑裂面的定义方法 ( a ) 中性面上的滑裂面；c o ) 断面a a ，用椭球定义滑裂面的空间形状 图2 1 4某一“f 常数”的列界面上底滑面点的确定 注意中性面上滑裂面的两个端部h = 0 ，因而有h = 0 。对于一个地形起伏较大 的坡面，根据中性面坡顶的坐标( ，瑶) 来确定其它列界面的滑面坐标耐，硝) ，有 可能导致该坐标值处于其它剖面的边坡面上的情况，此时，既可通过调整值来 避免出现这类情况，也可由程序予以判断，删去处于边坡面以上的条块。 式( 2 4 0 ) 中的强行断面上椭圆的长短轴比例系数。如果刍等于1 ，该断面上的 滑动面为圆形。对于0 断面上控制滑裂面的每一个控制点，我们都指定一个螽值。 与孙弘和西的处理相同，两个控制点中间的各界面的当由插值确定。每一行断 面上椭球与边坡面相交最大的z 值，由下式给出： = f ( 日+ b ) ( 2 4 2 ) 参数占使得滑裂面在坡顶和坡趾处可以具有一定长度，以满足地基承载力问 题的需要。因为坡顶和坡趾处h 为0 ，故有z m x = t ；b ，如果f 和b 都不为零，则 其值非零，这就可以使滑体的冠部具有一定长度。特别地，当b - - - o 时，问题可 回归到一般的三维椭球形滑裂面的问题。 3 l 中国水利水电科学研究院博士后研究报告 2 5 2 目标函数的确定 当滑裂面与破坏模式按上述方法建立后，则安全系数f 可表述为与中性面相 关的一组变量的函数： f = f x l ，y l ，心，儿，x 。，4 ，磊，以，- ，白，白，厶)( 2 4 3 ) 边坡稳定极限分析要求寻求一临界滑裂面与临界破坏模式，使安全系数f 达到最 小。这个破坏模式包括： ( 1 ) 中性面的各控制点的坐标值工” ( 2 ) 各控制点的倾角舀 ( 3 ) 用于描述椭球面形状的各控制点的箍。 在实际中，我们可以针对不同的问题，将上述某些变量固定不动，以减少优 化过程中的自由度的个数。 2 5 3 常见的最优化分析方法 1 单形法 单纯形法是一种无约束的最优化直接方法。所谓单纯形是指n 维空间中具有 时1 个顶点的凸多边体。其计算方法为： 对某一初始向量z 口，按下面模式构筑 个向量z ( 卢1 ，2 ，功，组成单形， z 1 = 【z ? 十p ，2 ；+ g ，一，o + g 】r z 2 = 咖纠+ ”，o 7 ( 2 卿 z ”= 【z ? + g ，z + g ，一，z ，o + p r 其中 口：逊! 鎏! 二! 口 1 n 4 2 痒：鲤4 4 n 4 2 式中：口为选定的步长按照一定的方式通过反射， 新逼近极值点。收敛准则为 ( 2 4 5 ) ( 2 4 6 ) 扩充和收缩，使单形不断更 第2 章三维边坡稳定机动分析方法 岛黔约一，2 卜 其中 y z z 4 伽k = o + 1 ) ( 2 4 8 ) + l j 2 牛顿法 牛顿法通过解析手段寻找使目标函数f 对自变量z i 的偏导数为零的极值点 ( 0 2 w o 老i = o , i = 1 ，2 ，功。同时，从理论上讲，还需要满足由二阶导数形成的h e s s i a n 矩阵正定( 汐2 汐两2 d ，i = 1 ，2 ，哟。这是达到极小值的充分条件。此类方法中 以导数为研究的主要对象，因此，也称为以导数为基础的方法( g r a d i e n t - - b a s e d m e t h o d ) 。一般认为，当自由度较多时，直接搜索法效率很低。此时需要考虑牛 顿法体系的分析方法。这里对负梯度法作简单的介绍。 负梯度法的基本思想是对一个初始滑裂面，寻找一个使安全系数减少速率最 大的方向。在数学上，就是匿，善，善) 7 这个向量。在这个方向上，进行 一次搜索，找到这一方向安全系数的低谷点。完成了这第一次迭代后，再在这个 新的起点( 即上述低谷区) 重复这样的运算，直到收敛至极值点。 3 随机搜索法 随机搜索法是一种最原始的非数值优化算法。这种算法的基本思想是，首先 确定待优化变量的变化范围，然后让计算机产生一系列的在( o ，1 ) 范围内均匀分 布的随机数，对各变量的搜索范围进行扫描，然后找出其中的最优解。例如，若 目标函数中一待优化变量蔚，其搜索区域的轴线用r 表示，宽度为岛，半带宽 d ，= d f l 2 。在搜索区域内，该变量可按下式取值： 墨= r + ( o 5 一口) z( 2 4 9 ) 式中口为一系列在( o ，1 ) 范围内均匀分布的伪随机数。 根据m o n t ec a r l o 原理，当搜索次数趋于无穷大时，随机搜索法获得的最优解 就是整体的极值点。尽管直接采用随机搜索作为一种优化算法也有成功的例子 ( s i e g e l 1 9 8 1 ) ，然而，对一些多自由度的问题，通常需要进行数万次的扫描，才 中国水利水电科学研究院博士后研究报告 能达到足够的精度。在实际应用中，我们通常将随机搜索与确定性方法结合起来， 以弥补这两种算法的缺点。通常的做法是，首先利用随机搜索法在最优解的邻域 内寻求一个点，然后再以该点作为最优化计算的初值。c h e r t ( 2 0 0 1 ) 采用随机搜索 与单纯形相结合在边坡稳定三维极限分析上限法中的成功应用，表明了这种混合 优化算法在求解一些多自由度的复杂问题是完全可行的。 4 模拟退火法 模拟退火算法( s a ) 是一种利用m o n t ec a r l o 迭代求解策略的随机寻优算法，它 基于物理中固体物质的退火过程的“熵最大”原理。固体退火过程一般是这样的， 先将固体加热至融化，然后使温度“徐徐”降低。在每一个“徐徐”降温的过程， 都会使系统在这一温度下达到平衡状态。从统计学的观点来看，随着温度的降低， 物质的能量将逐渐趋近于一个较低的状态，并最终达到某种平衡状态。1 9 8 3 年 k i r k p a t r i c k 等意识到物理退火与优化问题的相似性，提出了模拟退火算法。与传 统的确定性方法相比，模拟退火算法的优点在于能在局部极小点处根据 m e t r o p o l i s 准则来判断是否接受新的状态，因而是一种全局优化算法。 5 遗传算法 遗传算法是一种基于生物界自然选择和进化机制发展起来的高度并行、随机 与自适应性的搜索算法。它主要借鉴了达尔文关于生物进化的基本规律，即为保 证物种的稳定性，子代都将继承父代的一些特征，但是子代与父代之间及子代的 不同个体间总有些差异。在自然选择的过程中，具有适应度变异的个体将保留下 来，而不具有适应度变异的个体将被淘汰。经过一代代的生存环境的选择作用， 物种变异将最终朝着向适应度大的方向进化。 遗传算法模拟了自然界中生物进化的过程，是一种通用的全局优化算法。它 是从代表问题可能的潜在解集的一个初始种群开始的，一个种群是由若干个个体 组成，每个个体实质上是待优化问题的一个可行解。遗传算法首先要对优化问题 的解进行编码，我们称一个解的编码称为个染色体，组成编码的元素称为基因。 目前最常用的编码方式是采用二进制编码。当初始种群产生以后，就可以计算出 每个个体的目标函数值，进而可求出其对应的适应度。当初始种群的适应度确定 后，自然选择规律按照适应度的大小以概率方式决定其中生存下来的染色体，以 形成一个可以繁衍下一代的种群。借助于自然遗传学的遗传算子，父代的两个染 第2 章三维边坡稳定机动分析方法 色体相互交叉，形成了下一代的种群，即形成了优化问题的一个新解。新解在产 生的过程中也可能发生基因变异，使解有更大的遍历性。以新解作为父代，重复 上述步骤，直至寻找到最优解或达到某一个迭代终止的标准为止。 2 6算例验证 2 6 。l 楔体算例 为对本文的方法进行验证，考察了具有闭合解的两个楔体算例f 具体几何参 数见表2 1 。显然，对于本节所讨论的四面体楔体，中性面应通过两个滑面的交 棱线。为了与传统解进行对比，本节还引入以下两个假定： 1 楔体作为一个整体呈刚体运动，即条柱间无相对位移，故可认为界面的i 脚， c - - 0 ； 。2 楔体滑动时，与左右两滑面分别按剪胀角y 。，沙。进行剪胀，其中， 满足0 九，0 s 丸。 表2 i楔形体算例的倾向倾角 1 对称楔形体算例 本算例为一几何形状对称楔体，其剖面划分如图2 1 5 所示，其物理力学参数 为：左右两滑面采用相同的抗剪强度指标：c = 5 0 k p a ，矿- - 2 0 。，岩石容重为2 6 k n m 3 ( 数据文件s 6 d a 0 。楔体不受除重力外的其它外荷载作用。若假定楔体左右两滑 面的剪胀角相等，且均在( 0 ，力的范围内变化，表2 2 列出了滑面不同剪胀角对应 的楔体的安全系数。考虑计算误差影响，旷0 。的结果与传统方法的静力解1 5 5 6 一致。 从表中不难看出，楔体的安全系数随剪胀角渺的增加而逐渐增大，这与楔体 广义解的结果是完全一致的。当 u = 2 0 。时楔体的安全系数达到最大，与基于相关 中国水利水电科学研究院博士后研究报告 联流动法则的上限解的计算结果完全一致。 图2 1 6 - - 2 1 8 为计算得到的不同剪胀角对应的楔体各条柱速度场情况。由得到 的速度场可知整个楔体作为刚体进行滑动，速度场与楔体广义解是一致的。 图2 1 5 几何形状对称的楔体计算剖面 表2 2楔体的不同剪胀角洌应的安全系数f 剪胀角妒( o ) 051 01 5 2 0 极限平衡解 安全系数， 1 5 5 81 5 6 71 5 7 31 5 7 71 5 7 91 5 5 6 ( a ) 一次速度场v = ( - 0 们z , - 0 3 9 ) t( ”二次速度场e - - ( - o 9 - 0 1 ， 图2 1 6 对称楔体第一列条柱速度场的示意图( 剪胀角p 卸。) ( a ) 一次速度场v - 0 9 7 , 0 2 6 ) 7( ”二次速度场q = ( - 0 9 8 ，- 0 1 ) t 图2 1 7 对称楔体第二列条柱速度场的示意图( 剪胀角p t l o 。) 3 6 第2 章三维边坡稳定机动分析方法 一次速度场v 鼍- 0 9 9 , - - 0 1 ) t ( ”二次速度场q 鼍0 9 9 | - 0 1 ) t 图2 1 8 对称楔体第二列条柱速度场的示意图( 剪胀角￥- z 0 。) 2 不对称楔体算例 本算例考察一个几何形状不对称的楔形体的情况，计算剖面如图2 1 9 所示。 楔体左右滑面的力学性质为：c = 5 0 k p a ，萨3 0 。楔体的平均容重为2 6 k n m 3 ( 数 据文件s 3 d a 0 。 表2 3 列举了该楔体算倒机动解的计算结果，与楔体广义解完全一致。表2 4 为计算得到的不对称楔体的速度场情况，与广义解中的解析速度场是完全相同 的。 h = i0 0 m 图2 1 9 不对称楔体的计算削面 表2 3 不对称楔体的不同剪胀角对应的安全系数， 剪胀角妒( o ) 051 01 52 02 53 0 传统方法 安全系期 1 7 2 21 7 6 0i 7 8 8 i | 8 1 3 i 8 3 2i 8 4 5 l | 8 5 1 i 6 4 0 3 7 中国水利水电科学研究院博士后研究报告 表2 4 不对称楔体的不同剪胀角v 对应的速度场 真实速度场y 虚拟速度场q 剪胀角 = 燎方法 广义解法 = 维方溃 广义解法 栌：0 0 ( - 0 8 3 ，一0 5 6 ，0 )( - 0 ，8 3 2 , - 0 5 5 4 ，o )( 一0 9 8 5 ，一0 0 9 2 ，0 o g )( - 0 9 8 5 ，- 0 0 9 2 , - 0 0 8 ) 配户1 0 。 ( - 0 9 1 ，一04 2 ，0 0 3 )( o 9 0 9 ，- 04 i6 , 0 0 3 ) ( - 0 9 7 2 ，- 0 0 9 ，0 0 8 )( 一0 9 7 2 , - 0 0 9 ，0 o s ) g = 2 0 。 ( 一o 9 6 2 ，- 0 2 6 5 ，0 0 5 ) ( - 0 ，9 6 2 ，一0 2 6 5 ，0 0 5 )( - 0 9 5 2 ，- 0 0 8 ，0 o s )( - 0 9 5 2 ，一0 0 8 ，0 0 8 ) 戽3 0 0 ( - 0 9 9 3 ，- 0 0 9 ，0 0 8 ) ( - o 9 9 3 ，- 0 0 8 ，0 0 8 )( - 0 9 2 2 ，- 0 0 8 ，0 0 7 )( - 0 9 2 2 ，- 0 0 8 ，0 0 7 ) 2 6 2 椭球体算例 图2 2 0 为z h a n gx i n g 提供的椭球体滑面算例，其中包括两个算例情况，算例 1 为一滑裂面为椭球型的均质边坡，算例2 的滑面由部分椭球体与部分软弱夹层 组成。国内外很多学者都选择本例题来检验各自三维程序的合理性。 、 算倒l ；圜强屑裂：毡奶算倒l ；圜强屑裂面、j 嘞 算例2 ： ，耆 2 0 m 1 5 1 0 5 1 00 图2 。2 0z h a n gx i n g 椭球体滑面算例 图2 2 1 算例1 中椭球形滑裂面算例 ( a ) 滑块三维示意图；泐滑动后的三维地形图 第2 章三维边坡稳定机动分析方法 算例l 的滑裂面为一个对称的椭球体。在三维边坡稳定分析中，以滑体的对 称轴为其中性面，然后在中性面的两侧建立一系列与x o y 面平行的列界面，如图 2 2 1 ( a ) 所示，所构造的椭球形滑面的形状见图2 2 1 ( b ) 。令滑面上的抗剪强度指标 为：4 = 2 0 0 ，c = 2 9 o k p a ，