0 引言

自Zienkiewicz^[1]于1965年第一次将强度折减系数引入有限元计算中后，数值模拟法在边坡工程中发展迅猛，如今已形成了多种成熟的算法^[2]，如离散元法^[3]、统一强度折减法、双强度折减法^[4]、非饱和土降雨入渗稳定系数求解^[5-9]等。边坡模型在采用数值模拟方法进行稳定性求解时学者们的关注点更多是在初始应力平衡后的强度折减问题上，根据最终计算结果确定边坡的稳定状态及位移、塑性区监测依据等，但初始应力场搭建对边坡稳定性计算的影响是非常大的，数值模拟法相比极限平衡法的优势就是可以查看位移、应力应变分布等多种判据，但采用错误的初始应力场后，边坡位移、应力应变分布将会有很大的不同，并且边坡的稳定系数分布会出现一定的偏差，但由于边坡稳定系数主要受抗剪强度参数的影响，变化相对较小，最不利稳定系数几乎没有变化，但对塑性区、边坡整体稳定系数分布及节点位移分布影响较大，甚至可能出现量级的差距，目前大部分学者均没有注意此问题。

边界尺寸的确定是建立数值模型的重要准备工作，边界尺寸的大小影响模型的应力应变分布，至今学者们没有统一的边界条件选取依据，一部分学者弱化模型计算的时长缺陷，将模型建的足够大，将尺寸影响降到最小^[10]；一部分学者将尺寸距离与坡高关系进行了对比，最新的研究成果以潘永坚等^[11]为代表对二维模型和三维模型的尺寸选取均进行了讨论，但文章中的三维模型实质为假三维模型，均属于平面应变问题，却得到了二维模型与三维模型不同的边界尺寸选取建议，网格划分尺寸过大，普遍适用性较差，且既有成果均是以稳定系数为分析依据，未将应力应变等判据考虑在内，因此笔者认为，应以典型案例模型为例，除了对上、下、左、右各边界距离进行单变量尺寸变化条件下边坡稳定系数的影响分析，还应分析土体的应力应变分布规律，进一步优化边界尺寸确定原则。

网格划分是边坡模型进行有限元法及有限差分法的关键步骤，网格节点是插值计算的基础和载体，位移与应力应变关系是由节点传递得到的，如施加应力、体力等实质均是施加在节点之上，网格尺寸大小影响数值模型节点数，这关系着模型的准确性。划分过密增加计算时间，且模型不容易收敛，划分过疏，影响计算结果的准确性，但到底影响多大关注相对较少。

基于此，引申了3个需进一步研究的内容，即边界尺寸确定、初始应力场搭建的方法选择以及边坡模型的网格尺寸确定，基于FLAC3D^[12]，以典型边坡案例为研究对象进行分析研究，旨在为边坡数值模型的建立与求解提供一些有意义的参考。

1 模型边界尺寸的确定 1.1 目前模型边界尺寸确定存在的问题

极限平衡法计算模型相比于数值计算模型边界尺寸的要求并不高，只要滑动面在边界之内就可以得到准确的结果，但由于数值计算模型需要考虑其整体的应力应变分布情况，因此必须要考虑边界效应的影响。

目前模型边界尺寸确定存在的问题主要为既有成果均是以稳定系数为分析依据，未将应力应变等判据考虑在内，应力应变分布受限会导致计算结果失真，产生位移突变、剪应变增量突变等情况，且分析边界尺寸时未制定标准模型进行对比，并未采用极限平衡法进行校核，采用从小边界尺寸往大边界尺寸推算的方法。此种推算方法在单变量调整过程中，模型的应力应变分布会因为部分边界尺寸过小而使结果失真。为了证明此结论，采用Dawson边坡模型^[13](图 1)以碎石土为例(表 1)进行计算分析，仅更改模型的下边界距离U(坡脚至模型底面距离)与坡高H的比值，其他边界距离与坡高比值(单变量调整)不变，均取0.5。采用极限平衡法中的Corps of Engieers#1法和Bishop法，结果如表 2所示，在U/H为0.25，0.5及0.75时，边坡稳定系数均为1.082，但当U/H为1.0和1.25时，边坡稳定系数产生了突变，变化为1.090，稳定系数突变点对应的比值为0.75。但当维持U/H为0.75，将右边界距离与坡高比值R/H从0.25扩大至0.5时，稳定系数由1.082变化为1.090，这说明产生此种情况主要是由于右边界尺寸偏小而无法真实体现应力分布，出现应力分布受限的情况从而使计算结果失真。

基于此，对Dawson边坡模型进行边界尺寸修改，模型推算方法采用从大边界尺寸往小边界尺寸进行推算，模型的下边界U/坡高H分别取1.75，1.5，1.25，1，0.75和0.5进行对比；模型的右边界R/坡高H分别取2.5，2，1.75，1.5，1.25，1，0.75和0.5进行对比；模型的左边界L/坡高H分别取2.0，1.75，1.5，1.25，1，0.75和0.5进行对比。采用张鲁渝，郑颖人等^[14]推荐的尺寸建立标准边坡对比模型，模型的下边界U/坡高H取1，右边界R/坡高H取2.5，左边界L/坡高H取1.5。将稳定系数、复合位移作为边界尺寸的判别指标，建立边坡模型进行计算分析，模型网格尺寸为0.5 m×0.5 m。

根据边坡相关监测规范及实际工程中对边坡位移监测的精度要求，可取稳定系数突变后的最小尺寸且复合位移变化相差在1 mm范围内的尺寸作为初拟最小单变量边界尺寸，据此建立新的模型边界尺寸标准，通过计算多个边坡模型的稳定性对比验证其正确性。

1.2 边坡模型边界尺寸确定

根据表 1土体物理力学参数建立边坡模型分析不同边界尺寸对边坡稳定系数的影响，对比计算结果如表 3所示。

(1) 下边界U/坡高H的边界尺寸确定

当单变量为下边界U/坡高H时，边界尺寸比值从1.75缩小至0.75的稳定系数均为1.090，当比值缩小至0.5时边坡模型的稳定系数变为1.082，模型计算精度受到影响。从复合位移判据分析，复合位移值随着模型下边界U/坡高H比值的减小而减小，由于通常边坡模型的应力集中均发生在坡脚处，因此在满足“安全距离”的前提下，下边界对边坡模型的影响较小，边坡模型位移没有发生突变异常。

根据稳定系数和复合位移结果，在其他尺寸不变的条件下，将模型下边界的尺寸U不小于0.75H作为下边界尺寸的初拟安全距离。

(2) 右边界R/坡高H的边界尺寸确定

当单变量为右边界R/坡高H时，边界尺寸比值从2.5缩小至1.75的稳定系数为1.090，1.5和1.25的稳定系数分别为1.086和1.082，这表明当右边界R/坡高H比值大于1.75时稳定系数趋于稳定。但结合复合位移判据来看，比值为1.75时的位移产生了突变，与实际不符，这是因为右边界尺寸R对边坡模型的精度影响较大，模型右边界R/坡高H比值不满足“安全距离”时，由于尺寸效应导致模型应力应变分布受限，因此在比值为1.75时产生了位移突变，当比值大于2.0后模型位移变化正常。

根据稳定系数和复合位移结果，在其他尺寸不变条件下，将模型右边界至坡顶边界的尺寸R不小于2H作为右边界尺寸的初拟安全距离。

(3) 左边界L/坡高H的边界尺寸确定

当单变量为左边界L/坡高H时，边界尺寸比值从2.0缩小至0.5的稳定系数均为1.090，这是因为边坡土质为非黏性土体，潜在滑动面为浅层滑动，破坏区域均在坡顶至坡脚范围内，因此稳定系数几乎无变化，但观察复合位移计算结果，比值为0.75和0.5时的位移均相较比值为1.0时的位移大，与右边界R/坡高H的边界尺寸对比结果相同，这是由于尺寸效应产生了应力应变分布受限区，导致位移发生了突变异常，这说明边坡模型计算精度不仅受土层性质影响，受尺寸效应影响同样明显，确定边界尺寸需要同时符合稳定系数与模型位移的实际变化规律。

因此，考虑土质边坡模型尺寸通用在其他尺寸不变条件下，将模型左边界至坡脚边界的尺寸L不小于1H作为左边界尺寸的初拟安全距离。

综上所述，根据初拟尺寸确定了土质边坡模型的推荐通用边界尺寸，即坡脚至下边界尺寸U不小于0.75H，右边界至坡顶边界尺寸R不小于2H，左边界至坡脚边界尺寸L应小于1H。据此建立新标准尺寸模型，计算结果如图 2所示，稳定系数为1.090，最大复合位移为0.56 mm，计算结果正确。

为进一步论证该尺寸的正确性，除本研究模型外，对文献[11, 13-14]中的3种模型均按照本研究推荐边界尺寸进行重新建模后分别进行计算对比，计算结果和极限平衡法计算结果对比如表 4所示，在该边界尺寸确定原则下建立的模型均与极限平衡法计算结果相差微小。

2 初始应力场的搭建 2.1 目前初始应力场搭建存在的问题

初始应力场最常用的搭建方式为弹性求解方式^[15-17]，该种求解方式通常是对其弹性模量增大一个量级进行弹性求解，收敛快计算时长短，是学者们最普遍的采用方法，但由于采用弹性本构模型，初始应力场求解完成后不会出现塑性屈服区域；第2种方法是更改黏聚力参数的弹塑性求解，通过大幅度提高黏聚力和抗拉强度进行计算平衡，这种方式求解处的初始应力场存在塑性屈服区域，比弹性求解更接近实际，但计算时长要比弹性求解法高出好几倍；第3种方法是直接弹塑性求解法，采用摩尔库伦模型参数先进行弹性求解获得初始应力场后再进行塑性求解稳定系数，该方法从实际应用上与水平应力和竖向应力之比的K0法(模型设置竖向及水平初始应力进行求解)相同，并且和第2种方法从本质上来说是等效的，只不过是分阶段求解，不用再人为增加一行提高强度参数的命令流。

目前初始应力场搭建存在的问题主要集中于没有结合研究问题而搭建相应的初始应力场，即初始应力场的适用性问题。初始应力场求解目的是为了真实地描述边坡目前阶段岩土体已存在的应力状态，同时为了更好地分析破坏阶段与现阶段的相对位移变化通常对初始应力场求解后进行位移清零操作。

2.2 常用的初始应力场搭建方法适用性分析

为分析常用的初始应力场搭建方法的适用性，以上文刚确定的边界尺寸模型为研究对象进一步进行分析, 岩土体参数同文献[1]。由于第3种方法和第2种方法等效，主要比较弹性求解与分阶段弹塑性求解，模型网格划分、参数幅值、尺寸大小(1 m×1 m)等均相同，为更好地对比方法所得结果差异，基于Flac3D自带的FISH语言编写边坡稳定系数分区命令流^[18]，由于采用强度折减法得到的破坏时位移并不具备参考性，采用应力平衡法进行二次计算，从塑性区应变增量、坡顶的监测位移曲线及稳定系数分布云图综合分析边坡模型的现状稳定性及不同初始应力求解方法的差别。

首先采用弹性求解方式搭建初始应力场，搭建初始应力场时步很短，共5 469步，耗时12 s，各向节点位移速度清零后赋予摩尔库伦本构模型及土体参数后进行强度折减法求解，共计算17 085步，总耗时10 min。稳定系数分布云图(图 3)可知最小稳定系数为1.020，从塑性区应变增量云图(图 4)可知边坡目前接近整体贯通，稳定系数较低，最大剪应变增量为5.146×10^-2，计算平衡时坡顶监测最大位移为5 cm，处于欠稳定-不稳定状态。采用极限平衡法中的Bishop法进行计算验证，稳定系数为1.028，结果相差不大。

采用分阶段弹塑性求解方式搭建初始应力场，首先将岩土体的黏聚力提高到100 kPa进行弹性求解，搭建初始应力场时步共7 180步，耗时28 s，用时为弹性求解法的2倍多。强度折减法求解共计算18 673步，总耗时17 min，虽计算时步相差不多，但比弹性求解法多耗时7 min。得到初始应力场后求得最小稳定系数为1.012，稳定系数分布云图(图 5)与弹性求解法分布结果基本相同，但稳定系数相比弹性求解法小0.008。最大剪应变增量为0.592×10^-2，是弹性求解法的11.5%，计算平衡时坡顶监测最大位移为0.58 cm，是弹性求解法的11.6%(图 6)。处于欠稳定-不稳定状态，由于最不利稳定系数主要是由岩土体抗剪强度决定，因此最不利稳定系数相差并不大，但由于初始应力屈服区域的存在与否决定了整个模型的稳定系数区域分布有一定的改变，弹性求解初始应力场没有塑性屈服区域，计算平衡时的位移与最大剪应变增量均是分阶段弹塑性求解法求得结果的8.5倍多，且潜在最不利稳定系数区域与剪应变增量区域也比分阶段弹塑性求解法要大一些。

基于上述分析，当仅需要确定边坡稳定系数判定边坡稳定状态及位移、剪应变等判据时，可采用弹性求解法求解初始应力场平衡，计算时步较少且用时节省接近1/2但精度相对较差；若除了以上判据还需要获取边坡位移监测依据，或抗滑桩、挡墙等支挡结构设置有位移要求时，建议获取边坡模型初始应力场时采用分阶段弹塑性求解法，如表 5所示。

3 模型网格尺寸的划分 3.1 目前模型网格尺寸划分存在的问题

网格尺寸的划分对边坡模型稳定系数的影响非常巨大，数值模拟的计算结果通常受限于网格数量与节点数量，网格越密精度越高，但当网格密集到一定程度后，再加密网格时计算结果几乎没有变化，但计算时长却在成倍增加。

网格的形式也影响模型的计算时长，二维模型主要采用楔形和四边形网格，三维模型对应采用四面体和六面体网格，一般来说，四边形和六面体网格更容易收敛且划分效果较好。

目前模型网格尺寸划分主要存在的问题集中于网格尺寸大小的确定，一个普通的边坡模型，网格采用1 m×1 m网格与0.5 m×0.5 m网格计算结果可以相差10%以上，但目前通常的做法是根据模型的尺寸选取一个中等网格尺寸进行划分求解，逐渐加密直至计算结果相差不大时，将此状态下的网格大小作为模型的最终网格尺寸，但此方法耗费时间较长，且要经若干次的网格划分才能得到。

笔者认为，可以找寻边坡的模型尺寸与网格及节点数量的关系，边坡模型坡体部分应加密网格，而其他部分如边界尺寸部位可通过尺寸梯度变化适当变疏，通过大量的算例验证，得到一个可以普遍适用的网格尺寸的确定方法。

3.2 基于边坡模型尺寸与网格数量的关系确定

为确定边坡模型尺寸与网格及节点数量的关系，建立不同高度的边坡模型，相关物理力学参数仍选取表 1，一级坡坡率采用1∶1，二级及以上坡率采用1∶1.25，8 m和12 m边坡一级坡到顶，其他坡高分级高度为8 m，分级台阶为2 m，基于前文的成果确定其边界尺寸，制定5种网格划分形式(表 6)，考虑稀疏尺寸过渡，一般边界区域的网格个数取边坡坡体区域网格个数的0.5倍。以网格及节点数量为模型内部变量，稳定系数作为主要判别标准进行关系分析。

根据表 6的计算结果可知，坡体区域每平米网格数为4个时，每个边坡模型的稳定系数与极限平衡法计算结果是最为接近的，平均误差仅2.31%，计算精度较为理想，平均耗时约10 min。当坡体区域每平米网格数为2个时，每个边坡模型的稳定系数与极限平衡法计算结果的平均误差为5.03%，平均耗时4 min。其他划分形式下的边坡模型由于网格较为稀疏，虽然耗时均在1 min以内，但因计算结果误差较大导致结果精度较低。因此，建议在边坡模型中，划分网格时应逐线划分，坡体区域每平米划分应不小于4个网格，相应不小于5个节点，当精度要求不太高且模型较大时，可将距离坡体区域每平米划2~4个网格，远离坡体的区域应采用尺寸变化梯度过渡增大，该区域网格对于边坡模型的稳定系数影响不大，但由于各区域网格节点需要较完美的衔接，因此变化梯度不应过大，建议一般边界区域的网格个数取边坡坡体区域网格个数的0.5倍。

4 结论

基于土质边坡数值模型稳定性精度的研究现状，引申出了边界尺寸的确定、初始应力场搭建的方法选择以及边坡模型的网格尺寸确定这3个问题，以典型边坡案例为研究对象对以上问题进行了分析研究，在以下方面取得了一些进展：

(1) 将不同边界尺寸下数值模型中的应力应变分布情况引入到模型边界尺寸的确定中，采用边坡稳定系数、复合位移作为判据，分析不同边界尺寸对边坡稳定系数的影响，取稳定系数突变点且复合位移变化相差在1 mm范围内的尺寸作为初拟单变量边界尺寸，并采用既有文献中的5组模型进行验证，局部优化了模型边界尺寸确定原则，在张鲁渝、郑颖人研究成果的基础上将模型坡顶至下边界尺寸U由2H建议为1.75H，右边界至坡顶边界尺寸R由2.5H建议为2.0H，左边界至坡顶边界尺寸L由1.5H建议为1.0H。

(2) 论述了弹性求解、分阶段弹塑性求解初始应力场的原理及方法，从塑性区应变增量、坡顶的监测位移曲线及稳定系数分布云图综合分析了不同初始应力求解方法的适用范围。弹性求解法计算收敛快速、稳定系数精度与分阶段弹塑性求解法相差不大，但由于边坡应力应变变化不能体现边坡的真实情况，建议仅在需要稳定系数判定边坡稳定性的情况下使用；分阶段弹塑性求解法稳定系数精度最高，初始应力场阶段已产生部分屈服区域，与K0法一致，原理明晰，模型应力应变分布更接近实际情况，但求解时间较长，建议在除需要稳定系数判定外，需要获取边坡位移监测依据、或抗滑桩、挡墙等支挡结构设置有位移要求时的情况下使用。

(3) 通过找寻边坡的模型尺寸与网格及节点数量的关系，确定了一个普遍适用的网格划分方法，建议在边坡模型中，划分网格时应逐线划分，坡体区域每平米划分应不小于4个网格，相应不小于5个节点，离坡体区域较远的区域可采用尺寸梯度变化渐变。

(4) 土质边坡的稳定性与边坡土体的物理力学参数取值关系密切，实际公路工程中根据不同土体的物理力学性质结合规范推荐坡率法及工程类比法综合确定坡率，边坡土体根据抗剪强度可分为黏性土和非黏性土两大类，而边坡模型边界尺寸与边坡土体特征关系较大，为确保所得结论的准确性及适用性，不仅以非黏性碎石土边坡为研究对象讨论了土质边坡的模型尺寸确定，而且采用多位学者的黏性土边坡论文模型进行稳定性计算案例比较，计算结果表明，该边界尺寸优化原则对土质边坡普遍适用。