加筋土挡墙地震稳定性破裂面随机搜索法

李 倩 ，凌天清 ，韩林峰 ，张瑞刚

（1.重庆交通大学土木工程学院，重庆 400074；2.重庆交通大学重庆市山区道路建设与维护技术重点实验室，重庆 400074；3.水利部长江水利委员会河湖保护与建设运行安全中心，湖北 武汉 430010）

加筋土挡墙因节省成本、易于施工、设计灵活、承受大变形且不易破坏的优点而广泛适用于土木工程中.加筋土挡墙必须具有足够的稳定性，以抵抗各种可能的外部荷载作用.目前，地震荷载作用下加筋土挡墙的稳定性分析多采用拟静力法，该方法假设墙后填土的破裂面为平面，将挡土墙所受到地震加速度等效为水平和竖直方向的单向加速度[1].在应用拟静力法进行稳定性计算时，由于加筋土挡墙的水平成层特点，一般结合水平条分法（HSM）进行分析，该方法由Lo 等[2]基于极限平衡理论首次提出.Shahgholi 等[3]提出了一种简化的HSM，将整个加筋区分成若干个水平土条，只考虑每个土条的竖向平衡和整体水平向平衡（不考虑力矩平衡），将N层水平土条的平衡方程个数减少为 2N+1 个.Nouri 等[4]提出了一种不同的HSM，称为 5N-1 公式，采用该方法计算了加筋土挡墙维持内部稳定所需的筋材临界长度.Shekarian 等[5]提出了一种用于加固回填土刚性挡土墙抗震分析的HSM 方法.

虽然拟静力法应用广泛，但该方法没有考虑地震波特性，为此，提出了拟动力计算方法.拟动力计算方法多采用正弦波模拟地震加速度，但只考虑垂直向上传播的正弦波，没有考虑向上传播的波与自由地面反射波之间的相互作用，由于该方法不考虑地面对加速度的放大影响，需在分析中假设一个合适的放大系数.Choudhury 等[6]采用拟动力法确定了地震作用下加筋土挡墙抵抗直接滑移破坏和倾覆破坏所需的加筋长度.Bellezza[7]将均质无黏性土体假设为Kelvin-Voigt 固体，通过求解一维波动方程得到填土体内的加速度，并计算了墙体所受土压力.徐鹏等[8]分析了黏土地基对加筋土挡墙稳定性的影响，提出了黏土地基加筋土挡墙的设计建议.

采用HSM 进行稳定性分析的关键是确定加筋土挡墙的破裂面位置和形态，Sil 等[9]采用圆弧状破裂面分析了均匀黏性土质边坡在地震作用下的动力性能.Prater[10]假定滑动面为圆形和对数螺旋形，研究了均质土坡的屈服加速度.在实际情况下，加筋土挡墙破裂面的位置往往是未知的，破裂面也不是光滑的，学者们研究了加筋土结构破裂面的确定方法.蒋薇等[11]基于极限平衡理论采用优化算法求得加筋土边坡临界滑动面.林永亮等[12]以极限分析理论为基础，假定破裂面为对数螺旋形，研究了加筋土边坡的临界高度.邓东平等[13]以简化Janbu 法为基础采用滑动面搜索法计算了边坡的整体稳定性.Hu 等[14]利用变尺度混沌优化算法确定了非圆弧滑移面.Malkawi 等[15]采用蒙特卡罗优化方法计算了边坡临界滑动面.Sarma 等[16]采用一种附加应力容许准则的极限平衡方法确定了边坡的临界滑动面.万文等[17]将加速混合遗传算法应用于边坡最危险滑动面求解之中.Ausilio 等[18]采用极限分析原理结合对数螺旋破裂面研究了加筋土挡墙的地震稳定性.Nimbalkar等[19]以水平条分法为基础结合回归分析对加筋土挡墙进行了拟动力稳定性分析，Basha 等[20]采用对数螺旋破裂面对加筋土挡墙进行了拟动力计算.

本文以一个随机变量为基础，采用随机搜索法结合水平条分法，对加筋土挡墙破裂面进行搜索计算，生成了一个多线段形式的破裂面，得到了加筋土挡墙临界破裂面的位置.考虑墙后填土为黏弹性介质，采用黏弹性地震加速度计算地震惯性力，对加筋土挡墙进行了地震稳定性拟动力分析.

1 加筋土挡墙破裂面的确定方法

1.1 水平土条的划分方法

水平条分法计算模型如图1 所示，图中：点A、C分别为坡面顶点和坡脚点；点B为墙后填土表面的破裂点；α为挡墙倾角；H为挡墙高度；D为水平土条厚度；n为筋材总层数；Vi、Vi+1分别为作用于第i个水平土条上、下两侧竖直条间力；Hi、Hi+1分别为作用于第i个水平土条上、下两侧水平条间力；qhi、qvi分别为第i个水平土条水平、竖直向地震惯性力；zi和Wi分别为第i个水平土条的位置深度和质量；Si、Ni分别为第i个水平土条破裂面上的切向与法向力；βi为第i个水平土条破裂面与水平方向的夹角；Ti为第i个水平土条中筋材的拉力；li=D/sin βi为第i个水平土条的破裂面长度.

图1 水平条分法计算模型Fig.1 Calculation model of horizontal slice method

根据极限平衡原理，第i个水平土条的水平向与竖直向平衡方程为

考虑平衡条件，水平土条在破裂面处存在式（3）表示关系.

式中： τf为失效剪应力；τr为维持平衡所需剪应力；Fs为安全系数，假定对所有水平土条安全系数取相同值.

根据Mohr-Coulomb 破坏准则，第i个水平土条破裂面上Si和Ni满足式（4）所示关系.

式中： φ 为土的内摩擦角；c为土的黏聚力.

滑动楔体水平方向的整体平衡方程为

将式（1）～（4）代入式（5），令Fs=1，筋材总拉力为

将墙后填土视为黏弹性介质，采用Kelvin-Voigt模型，根据横波与纵波在平面内传播的运动方程可得Kelvin-Voigt 黏弹性介质运动方程为[21-23]

式中：uh为土体水平位移；uv为土体竖向位移；t为时间；z为土体深度；ρ 为土体密度；G为土体剪切模量；λ 为拉梅常数；ηs为横波传播黏性系数；ηp为纵波传播黏性系数.

Bellezza[21]利用阻尼比公式和归一化的横波频率公式，并假设挡墙模型下方的刚性层受到谐波激励，地表的剪切应力和法向应力为0，挡墙底部与刚性基底的位移一致，得出地震波在Kelvin-Voigt 黏弹性土体中的水平加速度为

式中：ωs为横波频率；kh为水平地震加速度；

同理，地震波在Kelvin-Voigt 黏弹性土体中的竖向加速度为

式中： γ 为土的重度.

作用在第i个水平土条上的水平向地震惯性力和竖直向地震惯性力分别如式（11）和式（12）所示.

将式（11）、（12）代入式（6），可求出筋材总拉力，筋材总拉力可用一类似于土压力的无量纲参数K表示，即

1.2 多线段破裂面形状及生成

假定加筋土挡墙滑动破裂面由多个倾斜线段组成，这些线段在一个平面内以不同的长度和角度相互连接，如图2 所示.多线段破裂面的倾斜线段数与筋材层数相同，若加筋层数有n层，即有n个水平土条，则滑动破裂面中的倾斜线段有n条，随着层数的增加，多线段滑动破裂面逐渐逼近曲线形态.

图2 多线段破裂面形式Fig.2 Form of polyline fracture surface

多线段破裂面采用产生随机角的方式生成，如图3 所示.Lc为填土表面破裂点与坡面顶部之间的距离，∠CBE为生成随机角所需的最大角度，其值可由已知几何条件求得.图3 中挡墙加筋层数有5 层，则程序采用产生的随机数生成5 条虚线，每条虚线与相应的水平土条界面的交点为1、2、3、4 以及坡脚点C，将所有的交点相连可得到多线段破裂面形状.利用产生的随机角计算得到 θ1～ θ4的值，得到每条倾斜线段的倾斜角度，可计算每层水平土条中的筋材拉力.

图3 多线段破裂面生成方法Fig.3 Generation method of polylinefracture surface

1.3 多线段破裂面随机搜索

加筋土挡墙多线段破裂面搜索程序采用Fortran语言编制.程序首先采用random_seed 函数产生随机数所需的种子，然后在随机角生成循环步的初始时刻，根据种子产生随机数，采用random_number函数将各随机数与最大角度相乘，便可得到各随机角的角度值，进而可得到每条倾斜线段的倾斜角度.根据倾斜线段利用水平条分法便可计算各水平土条中的筋材拉力.

如图4 所示，多线段破裂面搜索程序有两层循环计算，外层循环为墙后水平填土表面破裂点位置循环，内层循环为随机角度循环.为计算准确，在初始计算时需合理选择初始破裂点，在选定初始破裂点之后程序便开始产生随机角并计算每层筋材的拉力，进一步计算所有筋材的总拉力，为保证计算的精确性，需适当设置随机角循环次数.

图4 破裂面搜索计算流程Fig.4 Flow chart of fracture surface search calculation

在一次内循环计算结束时记录筋材总拉力最大值，最大值所对应的破裂面形状即为此次循环所得破裂面形状.一次内循环计算结束后，程序便选择新的破裂点进行下一次外循环计算，并记录此次外循环中筋材总拉力的最大值，如此往复，每选择一次破裂点便会计算得到一次筋材总拉力的最大值.程序最后对计算得到的筋材总拉力进行比较，最大值所对应的破裂面形状即为所求的加筋土挡墙临界破裂面.

1.4 多线段破裂面计算特点

在以往的研究中，加筋土挡墙地震稳定性计算多是预先假定滑动面形状或采用数学优化算法对破裂面位置进行求解.相比于预先确定破裂面形状的计算方法，本文提出的滑动破裂面随机搜索计算方法更具灵活性，且无需进行优化计算.

采用随机角度结合水平条分法的计算方法便于计算机编程，且计算方法中只有一个随机变量，计算过程简单，各计算因素易于掌握，相比于采用数学优化方法求解，具有计算简便的特点.

水平条分法不是本文所提破裂面计算方法应用的必需条件，该破裂面计算方法也可用于竖向条分分析方法，结合Janbu 法、Sarma 法等经典方法可进行边坡稳定性分析.对数螺旋破裂面的滑动轨迹与土的内摩擦角有关，该破裂面的应用仅限于均质土层问题，而本文所提出的计算方法可用于非均质土层的稳定性分析.

2 多线段破裂面的合理性验证与分析

2.1 与筋材总拉力和Lc/H的比较

为便于将本文方法与现有方法的计算结果进行对比分析，计算参数取值相同，设置 H=5 m，α=0°，φ=30°，Vs=100m/s，γ=18kN/m3，ξ=10%，kv=0，计算结果如表1 所示.

表1 筋材总拉力对比Tab.1 Comparison of total reinforcement forces kN/m

由表1 可知：地震作用下本文算法得到的筋材总拉力最大，产生这种差异的原因在于本文考虑了地震加速度的非线性分布.文献[3]采用拟静力计算方法进行计算，没有考虑加速度放大的因素；文献[19]的拟动力方法虽考虑了加速度相位的变化，但假设加速度变化幅值为常数；文献[20]认为加速度沿地面线性增加，但方法不满足边界处的应力条件，即加速度的放大值不受约束.

为分析计算方法的可靠性，在 kh=0.2时计算Lc/H随 φ的变化规律，并与文献[4,18]进行了对比，如图5 所示.由图可知：本文计算方法所得参数变化规律与文献[4,18]相同，即在相同条件下，Lc/H随 φ的增大而减小，挡墙倾角越小，Lc/H越小.本文算法计算结果偏小，所得墙后填土表面破裂点到坡面顶部的距离最小.文献[4]采用对数螺旋破裂面结合水平条分法，文献[18]采用对数螺旋破裂面结合能量原理对 Lc/H进行了分析.文献[4,18]虽然采用了相同的破裂面形式，但由于采用了不同的计算方法造成计算结果略有差距.本文计算结果小于文献[4,18]，主要是破裂面形态不同所致.

图5 Lc/H 与φ关系的比较Fig.5 Comparison of relationship between Lc/H and φ

2.2 与对数螺旋破裂面位置的比较

对数螺旋破裂面是加筋土挡土结构在稳定性计算时常用到的一种破裂面形状，如图6 所示.对数螺旋破裂面方程可用式（14）表示.

图6 对数螺旋破裂面形式Fig.6 Form of logarithmic spiral fracture surface

式中： θ 为极坐标系中的旋转角度；θ0为初始角度；θh为总角度；r为对数螺旋线的半径；r0为初始半径.

对数螺旋破裂面采用 θ0与 θh两个最重要的参数进行定义.在计算时，通过确定维持平衡所需筋材拉力最大值，来搜寻对数螺旋破裂面的位置.利用MATLAB 工具箱中fminsearch 函数，采用优化方法求解筋材拉力计算方程，使得筋材拉力达到最大值，求解过程中的约束条件为 0°＜θ0＜90°，0°＜θh＜90°.当 α=15°时，对数螺旋破裂面与多线段破裂面位置计算结果如图7 所示.由图可知：多线段破裂面计算方法所得加筋土挡墙的破裂面位置与对数螺旋破裂面较接近，但更靠近临空面，随着填土内摩擦角的减小两者计算结果差距略有增大，在此算例中最大相差9.6%.

图7 多线段破裂面与对数螺旋破裂面位置的比较Fig.7 Position comparison of polyline fracture surfaces and logarithmic spiral fracture surfaces

本文所提多线段破裂面计算方法与对数螺旋破裂面的不同之处在于，多线段破裂面采用随机角的方式生成，随机角的产生不受破裂面方程的制约.对数螺旋破裂面在优化过程中，主要优化θ0与θh两个角度参数，优化过程中破裂面受到对数螺旋方程的限制.而采用随机角生成的多线段破裂面，每一段破裂形状不受破裂面方程的限制，求解过程中形式变化灵活，对筋材总拉力的求解更为精确.由此导致多线段破裂面位置的求解结果相对于对数螺旋破裂面更接近挡墙临空面.

3 稳定性影响因素分析

假定H=5m，ξ=10%，γ=18 kN/m3，α 为0°～30°，φ为25°～45°，kh为0.1～0.2，计算不同地震加速度和挡墙倾角下，φ 与K、Lc/H的关系分别如图8、9 所示.

图8 α不同时K 随φ的变化规律Fig.8 Trend of K changing with φ under different α

3.1 α、kh 与K 的关系

由图8 可知：地震作用下，K随着 φ 的增大而减小，随着加筋土挡墙坡度的减小而减小，随着地震加速度的增大而增大.当 α=0°时，φ由2 5°增大到45°，kh由0.1增加到0.2 时，K增加了约19.7%和23.2%；在φ为25°且kh为0.2 情况下，挡墙倾角从0°增加到15°和从15°增加到 30°时，K分别减小了约21.3% 和37.7%；kv对K的影响不大，尤 其当φ较大时；当φ=45°，α=0°时，kh= 0.2，相比于kh=0.1时，K仅增加了8.3%.

3.2 α、kh与Lc/H 的关系

由图9 可知：地震作用下，Lc/H随着加筋土挡墙坡度的减小而减小，随着 φ 的增大而减小，该结论也可从文中对比分析算例图5 中得到验证.当不考虑α与φ的变化时，当kh=0.1，kv对Lc/H的影响可忽略；当kh=0.2 时，kv的变化对Lc/H的影响相比于K更明显.

图9 α不同时 Lc/H 随φ的变化规律Fig.9 Trend of L c/H changing with φ under different α

4 结 论

1） 将加筋土挡墙破裂面视为多线段形式，基于水平条分法提出了破裂面生成新方法，以筋材总拉力作为目标函数，采用角度随机的方式对目标函数进行搜索计算，确定加筋土挡墙破裂面的位置.

2） 所提出的方法由于考虑了地震加速度的非线性分布，计算得地震作用下筋材总拉力较大.与对数螺旋破裂面相比，多线段破裂面计算得加筋土挡墙破裂面位置更接近临空面，当填土内摩擦角为25°时，两者位置最大相差 9.6%.

3） 地震作用下，筋材总拉力随着加筋土挡墙坡度和地震加速度的增大而增大，筋材长度随着填土内摩擦角的增大而减小.当填土内摩擦角较大时，竖向加速度的变化对筋材总拉力和筋材长度的影响较小.