Influences of constitutive model for rockfill materials on calculated stress and deformation of concrete-faced dams
-
摘要: 堆石料本构模型选择是影响混凝土面板坝应力变形计算结果的主要因素。以某典型面板坝为例，采用邓肯E-B非线性弹性模型和“南水”双屈服面弹塑性模型，对坝体填筑过程和蓄水过程进行了三维有限元模拟；研究了两种模型计算的坝体和面板的位移和应力分布差异；分析了两种模型计算结果呈现差异的原因。两种模型计算结果最显著的差异体现在两个方面：①“南水”模型计算的坝体沉降、面板挠度等位移指标小于E-B模型，其原因是“南水”模型计算的坝体小主应力和变形模量更高。②E-B模型计算结果显示蓄水后面板底部顺坡向受拉；而“南水”模型计算的面板顺坡向应力为全断面受压。前者是由E-B模型的各向同性弹性本质决定的，垫层料因水压力作用顺坡向膨胀，产生沿坡面向上的位移；后者是其F1屈服面持续扩张产生塑性体积收缩导致的，塑性体缩抵消了顺坡向弹性膨胀，使垫层料产生沿坡面向下的位移。Abstract: Choice of constitutive model for rockfill materials is the most influencing factor that affects the computed stress and deformation results of concrete-faced dams (CFD). In this study, the Duncan's E-B (EB) nonlinear elastic model and the Nanshui (NS) double-yield surface elastoplastic model are used for rockfill materials to simulate the construction and impounding processes of a typical CFD by using three-dimensional finite element method. The differences in the displacement and stress of rockfill zones and concrete slabs are compared, and the reasons were analyzed. Two most important differences are found. First, the displacement quantities such as the dam settlement and slab deflection predicted by the NS model are smaller than those obtained by the EB model, and this difference can be attributed to the higher minor principal stresses and therefore higher deformation moduli within rockfill zones by the former. Second, the EB model predicts a tensile zone along the upstream slope near the bottom of face slabs, while the results by the NS model shows that the slope stresses of face slabs are completely compressive. The tensile stress predicted by the EB model is due to its isotropic elasticity nature, i.e., the cushion materials expand along the slope under the huge hydrostatic pressure, leading to a displacement upward towards the dam crest. The NS model, on the other hand, predicts an additional plastic volume contraction within the upstream rockfill zones due to the expansion of F1 yield surface, which completely counteracts the elastic expansion along the slope and results in a downward displacement of the cushion layer.
-
Keywords:
- concrete-faced dam /
- rockfill material /
- constitutive model /
- stress and deformation
-
0. 引言
边坡及坝基深层的抗滑稳定计算是岩土工程重要研究内容。目前，三维滑动面极限平衡抗滑稳定计算方法[1]主要是二维方法的拓展，由于极限平衡方法存在不同假设条件使得三维扩展失去了严格理论基础[2-3]。即复杂滑动面（非单一滑动平面）极限平衡方法本质属于力学静不定问题，通过静力平衡条件确定每个分块底面的法向力和切向力必须对分块界面上的作用力要素进行经验简化假设。
如果采用有限元方法计算得到的真实应力场，通过应力积分就可以确定每个分块底面的实际法向力Nm和切向滑动力Tm，由法向力可以确定分块底面的切向理论极限抗滑力Sm=fNm+cAm。把每个分块底面切向滑动力Tm矢量求和可得滑动力合成矢量T，同样把每个底面的切向抗滑力Sm矢量求和得出抗滑力合成矢量S。若能够得到潜在滑动面以矢量力表达的抗滑稳定极限平衡安全度系数，则所有抗滑稳定计算中存在的缺陷、简化假设问题可以彻底消除。然而，滑动力合矢量T与抗滑力合矢量S作用方向存在夹角δ为非共线不平衡力，非共线不平衡矢量力条件下的极限平衡抗滑稳定安全度系数计算研究成为世界级理论空白难题。
潜在滑动面的抗滑稳定安全度系数是表达实际受力情况的切向滑动力与滑面在实际受力条件下能够提供的理论极限概念状况抗滑力的接近程度，即两种状况的极限平衡接近程度。为了依据有限元应力场实际受力情况的滑动力矢量Tm与理论概念状况的极限抗滑力矢量Sm，研究矢量抗滑稳定极限平衡安全度系数，作者发现并提出了投影方向极值原理[4-5]，揭示了潜在滑动面应力矢量非共线不平衡特征条件的抗滑稳定极限平衡计算内在极值规律。
对于二维滑动面问题已在文献[5]中详细论证，本文根据三维实际应力场特性重点讨论投影方向极值原理的三维扩展模型（以下简称理论模型）及其具体使用程序步骤。三维潜在滑动面研究对象任意点i上的真实应力、极限平衡计算受力如图 1。图 1中τ与s/K为共线应力。
![]()
图 1 滑动面任意点的真实应力状况与极限平衡受力计算简图Figure 1. Real stress states of any point on slip surface and calculation of limit equilibrium forces单一滑动平面的极限平衡法为静定问题（安全系数K时的极限平衡概念表达式$ \frac{{(fN + cA)/K}}{T} = 1 $，分子为理论可用极限概念状况，分母为真实受力状况）。对于复杂滑动面极限平衡法本质为力学超静定问题的经验简化假设条件下的按静定平衡近似求解方法。其三维扩展时所增加的简化假设条件大于平衡方程数目，使得所有三维极限平衡计算法[6-17]都难于摆脱简化假设的窘境。但是，理论模型的三维扩展却十分简便，根本原因是理论模型已经囊括了抗滑稳定计算问题的普遍力学本质——复杂滑动面的滑动应力合矢量与理论极限抗滑应力合矢量普遍是非共线不平衡力的客观规律，只有特例情况才是矢量共线极限平衡问题（单一滑动平面或复杂滑动面在完全理想塑性的临界共线极限状态时，合矢量夹角δ＝0）。非共线不平衡力系的平衡计算投影方向及力矩中心点位置选取对安全系数的影响敏感，且其影响随着合矢量夹角的大小而变化，必须采用投影方向极值原理才能理论求解安全程度。相关问题的阐述及各种计算方法存在的缺陷参阅文献[5]，需要重点说明4点。
（1）投影方向极值原理揭示了复杂滑动面的滑动力矢量与极限抗滑力矢量为非共线不平衡力，其势能比具有正交极值特性。即确定安全系数K的极限平衡研究计算方向滑动势能比最大、确定极限状态接近程度系数λ的耗散方向势能最小。
（2）理论模型是根据滑动面应力矢量不平衡力的固有势能极值正交特性，把矢量非共线不平衡问题的极限平衡抗滑稳定计算分成两个独立正交极值方向分别研究矢量力的抗滑性能贡献。潜在滑动面研究对象的应力矢量极限平衡抗滑稳定安全度计算成为极值正交双向计算理论。
（3）有限元应力场计算成果原本是平衡的，为了根据真实应力场计算极限平衡抗滑稳定安全度，使滑动面两侧的切向滑动应力合矢量与理论可用极限抗滑应力合矢量成为非共线不平衡力系统。那么，在抗滑稳定完整计算过程中必须使该系统达到极限平衡，理论模型通过极限状态程度系数$ \lambda $归一化——得到塑性耗散力Fts为0的共线极限平衡安全度系数$ {K_\lambda } = K \div $ $ \lambda $，客观理论表达了真实受力应力状况和极限抗滑力概念状况的矢量力共线极限平衡抗滑稳定安全度。
（4）考虑三维滑动面应力矢量极限平衡时，由于滑动面与单元棱线每个交点的应力状态不同，极限抗滑应力方向按交点滑动剪应力方向确定时，由交点及其形心点组成的单元内小三角形滑平面（简称微面）上的切向滑动力矢量与抗滑力矢量作用方向在微面内存在微小偏角。因此，三维滑动面应力矢量极限平衡必须考虑两者作用方向偏角影响（二维情况单元的滑动力与抗滑力方向均沿滑面方向不存在偏角）。
1. 三维滑动面投影方向极值原理模型
1.1 空间滑动面投影方向极值原理
理论模型研究对象是潜在滑动面，在真实受力条件的整体极限平衡安全度（为了便于理解和表述可把三维潜在滑动面视为具有微小厚度的空间壳体）。滑面上侧受力为考虑了滑动体应力应变本构关系和满足变形协调要求的有限元计算成果应力，即切向应力τ和法向应力$ \sigma $，滑面分布应力与滑动体真实受力特性具有平衡等价功效；滑面下侧受力为理论极限概念的可用切向应力$ s/K = (f\sigma ' + c)/K $和法向应力$ \sigma ' $，其中$ \sigma ' = \sigma $（大小相等分布相同、作用方向相反）；滑面两侧的切向应力τ与s/K为共线应力。滑面通过实体单元的三角形微面矢量极限受力简图如图 2。
![]()
图 2 滑面通过单元三角形微面矢量极限平衡计算简图Figure 2. Calculation of vector limit equilibrium on triangle micro-plane of sliding surface over element微面上侧的切向合力为Tm、法向合力Pm；微面下侧的切向合力为Sm/K、法向合力Nm；微面两侧的法向合力为平衡力Nm=Pm，而切向合力Tm和Sm/K为非共线不平衡力。法向平衡合力中不包括考虑开裂时需要释放的拉力。
潜在滑动面两侧的切向力、法向力矢量求和，得切向滑动力合矢量$ T = \sum\limits_m {\overrightarrow {{T_m}} } $、可用极限抗滑力合矢量$ S/K = \sum\limits_m {\overrightarrow {{S_m}} } /K $、法向平衡合矢量$ - N = P = \sum\limits_m {\overrightarrow {{P_m}} } $。绘制潜在滑动面各微面三维矢量图，沿合矢量S和T所确定平面（记为ST平面）的法线方向观察，即把滑动面的所有矢量力投影到ST平面上，如图 3所示（因法向平衡矢量Nm=-Pm为反向重叠闭环，对矢量平衡没有影响，图中未示出）。
![]()
图 3 三角形微面矢量力在合矢量平面的投影示意图Figure 3. Projection sketch of vector forces of triangular micro-plane in resultant vector plane图 3中$ {T'_m} $，$ {S'_m} $分别是微面m的滑动力矢量Tm、抗滑力矢量Sm在ST平面的投影；$ {\alpha_m} $，$ {\varepsilon_m} $分别为$ {T'_m} $，$ {S'_m} $与合矢量平面局部坐标轴$ x $方向的倾角（二维情况$ {\alpha_m} = {\varepsilon_m} $）；微面m的i点的滑动应力$ {\tau_i} $和抗滑应力$ {s_i} $在ST平面的投影为$ {\tau '_i} $和$ {s'_i} $，则矢量积分$ {T'_m} = \int_m {\overrightarrow {{{\tau '}_i}} \text{d}{A_m}} $、$ {S'_m} = \int_m {\overrightarrow {{{s'}_i}} \text{d}{A_m}} $；$ {F'_n} $和$ {\omega_n} $分别为微面n拉应力大于允许值考虑开裂时需要释放的拉力矢量在合矢量平面的投影和倾角，矢量$ {{F}^{\prime }}_{n}={\displaystyle {\int }_{n}\overrightarrow{{{\sigma }^{\prime }}_{i拉}}}\text{d}{A}_{n} $；$ \theta $为合矢量平面内投影极值方向与局部坐标轴x的夹角。
矢量图 3表明：
（1）真实应力场时给定三维滑动面的合矢量夹角δ为固定值与安全系数K无关，但不平衡耗散力的方向和大小随安全系数K的变化而变化。
（2）计算安全系数K的平衡投影方向θ应该垂直于不平衡耗散力（K2时垂直于AB，K1时垂直于AC），才能消除耗散力的影响。但安全系数是未知量，平衡投影方向θ与安全系数有关必然也是未知量。
（3）投影方向只能得到1个平衡方程，不可能理论定解K和θ两个未知量，若要得到理论解必须采用投影方向极值原理的各点不平衡力势能极值条件方程。
（4）投影方向极值原理在图 3中的具体表达：在垂直于AB的投影方向各矢量力达到极限平衡，同时该方向的滑动势能达到最大；沿AB方向的不平衡耗散力最小，同时AB方向的耗散势能也最小。
（5）如果不平衡耗散力Fts为0，则A点与B点重合、合矢量夹角δ=0，滑动面在安全系数K达到整体塑性临界共线极限平衡状态，滑动面的矢量极限平衡安全度系数等于K。
1.2 空间滑动面理论模型建立
依三维滑动面的矢量图 3，在耗散方向的耗散力
$$ {F}_\text{ts}={\displaystyle \sum\limits_{m}\left[{{T}^{\prime }}_{m}\mathrm{sin}(\theta +{\alpha }_{m})-\frac{{{S}^{\prime }}_{m}}{K}\mathrm{sin}(\theta +{\varepsilon }_{m})\right]}+{\displaystyle \sum\limits_{n}{{F}^{\prime }}_{n}\mathrm{sin}(\theta +{\omega }_{n})}。 $$ (1) 将矢量力$ {\boldsymbol{T}'_m} = \int_m {\overrightarrow {{{\tau '}_i}} \text{d}{A_m}} $，$ {S'_m} = \int_m {\overrightarrow {{{s'}_i}} \text{d}{A_m}} $，$ {{F}^{\prime }}_{n}={\displaystyle {\int }_{n}\overrightarrow{{{\sigma }^{\prime }}_{i拉}}}\text{d}{A}_{n} $代入，并把滑动面应力矢量积分写成标量积分可得
$$ {F_\text{ts}} = \sum\limits_m {\int_m {} \left[ {{{\tau '}_i}\sin (\theta + {\alpha_m} + {\beta_{\tau i}}) - \frac{{{{s'}_i}}}{K}\sin (\theta + {\varepsilon_m} + {\beta_{\text{s}i}})} \right]} \text{d}{A_m} + $$ $$ {\displaystyle \sum\limits_{n}{\displaystyle {\int }_{n}{{\sigma }^{\prime }}_{i拉}}\mathrm{sin}(\theta +{\omega }_{n}+{\beta }_{\sigma i拉})}\text{d}{A}_{n} 。 $$ (2) 式中：βτi，βsi，βσi拉分别为微面的应力矢量与其合矢量方向的夹角。
令$ {\alpha_i} = {\alpha_m} + {\beta_{\tau i}} $，$ {\varepsilon_i} = {\varepsilon_m} + {\beta_{si}} $，$ {\varpi }_{i}={\varpi }_{m}+{\beta }_{\sigma i拉} $，并把三维滑动面各微面求和式$ \sum\limits_m {\int_m {f(x)\text{d}{A_m}} } $简写成$ \int_A {f(x)\text{d}A} $，得不平衡耗散力为
$$ {F_\text{ts}} = \int_A {\left[ {{{\tau_i}'}\sin (\theta + {\alpha_i}) - \frac{{{{s'}_i}}}{K}\sin (\theta + {\varepsilon_i})} \right]} \text{d}A + $$ $$ {\displaystyle {\int }_{A}{{\sigma }^{\prime }}_{i拉}\mathrm{sin}(\theta +{\omega }_{i})}\text{d}A 。 $$ (3) 在投影方向的平衡条件为
$$ {\displaystyle \sum\limits_{m}\left[{{T}^{\prime }}_{m}\mathrm{cos}(\theta +{\alpha }_{m})-\frac{{{S}^{\prime }}_{m}}{K}\mathrm{cos}(\theta +{\varepsilon }_{m})\right]}+{\displaystyle \sum\limits_{n}{{F}^{\prime }}_{n}\mathrm{cos}(\theta +{\omega }_{n})}=0。 $$ (4) 同样代入矢量力的应力积分和式（3）相同的简写变化过程得到：
$$ \int_A {\left[ {{{\tau '}_i}\cos (\theta + {\alpha_i}) - \frac{{{{s'}_i}}}{K}\cos (\theta + {\varepsilon_i})} \right]} \text{d}A + $$ $$ {\displaystyle {\int }_{A}{{\sigma }^{\prime }}_{i拉}\mathrm{cos}(\theta +{\omega }_{i})}\text{d}A=0 。 $$ (5) 由耗散力最小极值条件$ \frac{{\partial {F_\text{ts}}}}{{\partial \theta }} = 0 $或投影方向平衡条件均可得
$$ K=\frac{{\displaystyle {\int }_{A}{{s}^{\prime }}_{i}\mathrm{cos}(\theta +{\varepsilon }_{i})\text{d}A}}{{\displaystyle {\int }_{A}{{\tau }^{\prime }}_{i}\mathrm{cos}(\theta +{\alpha }_{i})\text{d}A}+{\displaystyle {\int }_{A}{{\sigma }^{\prime }}_{i拉}\mathrm{cos}(\theta +{\omega }_{i})\text{d}A}} 。 $$ (6) 滑动面的不平衡耗散力Fts是由于三角形微面的共线不平衡应力矢量$ ({\tau_i} - {s_i}/K) $引起的，那么，该剪切应力矢量在ST平面内沿投影方向的分量为$ \left[ {{{\tau '}_i}\cos (\theta + {\alpha_i}) - \frac{{{{s'}_i}}}{K}\cos (\theta + {\varepsilon_i})} \right] $，此分量对应的剪切应变为$ \left[ {{{\tau '}_i}\cos (\theta + {\alpha_i}) - \frac{{{{s'}_i}}}{K}\cos (\theta + {\varepsilon_i})} \right]/G $。三角形微面的滑动势能为$ {\int_m {\left[ {{{\tau '}_i}\cos (\theta + {\alpha_i}) - \frac{{{{s'}_i}}}{K}\cos (\theta + {\varepsilon_i})} \right]} ^2}\frac{{\text{d}{A_m}}}{G} $，则滑动面在投影方向的滑动势能（并考虑释放拉力）为
$$ {W_\theta } = {\int_A {\left[ {{{\tau '}_i}\cos (\theta + {\alpha_i}) - \frac{{{{s'}_i}}}{K}\cos (\theta + {\varepsilon_i})} \right]} ^2}\frac{{\text{d}A}}{G} + $$ $$ {\displaystyle {\int }_{A}\frac{{{\sigma }^{\prime }}_{i拉}^{\text{2}}{\mathrm{cos}}^{2}(\theta +{\omega }_{i})}{E}}\text{d}A 。 $$ (7) 同理可得耗散方向耗散势能为
$$ {W_{ \bot \theta }} = \int_A {{{\left[ {{{\tau '}_i}\sin (\theta + {\alpha_i}) - \frac{{{{s'}_i}}}{K}\sin (\theta + {\varepsilon_i})} \right]}^2}} \frac{{\text{d}A}}{G} + $$ $$ {\displaystyle {\int }_{A}\frac{{{\sigma }^{\prime }}_{i拉}^{\text{2}}{\mathrm{sin}}^{2}(\theta +{\omega }_{i})}{E}}\text{d}A 。 $$ (8) 总势能为
$$ W = {W_\theta } + {W_{ \bot \theta }} = \int_A {\left( {{{\tau '}_i}^2 + \frac{{{{s'}_i}^2}}{{{K^2}}}} \right)} \frac{{\text{d}A}}{G} - $$ $$ 2{\displaystyle {\int }_{A}\frac{{{\tau }^{\prime }}_{i}{{s}^{\prime }}_{i}}{KG}\mathrm{cos}({\alpha }_{i}-{\varepsilon }_{i})\text{d}A}+{\displaystyle {\int }_{A}{{\sigma }^{\prime }}_{i拉}^{\text{2}}}\frac{\text{d}A}{E} 。 $$ (9) 由滑动势能最大极值条件$ \frac{{\partial {W_\theta }}}{{\partial \theta }} = 0 $或耗散势能最小极值条件$ \frac{{\partial {W_{ \bot \theta }}}}{{\partial \theta }} = 0 $均可得
$$ a\frac{1}{{{K^2}}} - 2b\frac{1}{K} + c = 0 。 $$ (10) 式中：$ a = \int_A {\frac{{s_i{'}^2}}{G}} \sin (2\theta + 2{\varepsilon_i})\text{d}A $；
$$ b = \int_A {\frac{{{{\tau '}_i}{{s'}_i}}}{G}\sin (2\theta + {\varepsilon_i} + {\alpha_i})} \text{d}A \text{；} $$ $$ c={\displaystyle {\int }_{A}\frac{{{\tau }^{\prime }}_{i}^{2}}{G}}\mathrm{sin}(2\theta +2{\alpha }_{i})\text{d}A+{\displaystyle {\int }_{A}\frac{{{\sigma }^{\prime }}_{i拉}^{\text{2}}}{E}\mathrm{sin}(2\theta +2{\varpi }_{i})\text{d}A} 。 $$ 由极值条件方程式（10）得
$$ K = \frac{a}{{b \pm \sqrt {{b^2} - ac} }} 。 $$ (11) 联立式（6），（11）两个极值条件超越方程，可以唯一定量求得安全系数K和相对于合矢量平面局部坐标轴的投影极值方向角θ。
G，E分别为第i计算点处的剪切模量和弹性模量。当考虑拉应力大于允许拉应力开裂时，开裂点处仅有需要释放的应力没有任何抗力；当不考虑拉应力大于允许拉应力开裂时，有关$ {{\sigma }^{\prime }}_{i拉} $项从模型中删除即可。本文及文献[5]为开裂后理论模型，文献[4]为开裂前考虑抗拉理论模型。
把联立方程求解得到的投影极值方向角θ及安全系数K代入式（3）可得耗散力Fts，则耗散方向极限状态接近程度系数：
$$ \lambda = 1 - \frac{{{F_\text{ts}}}}{{T\cos (\theta - \mu )}} 。 $$ (12) 滑动面应力矢量极限平衡抗滑稳定安全度系数：
$$ {K_\lambda } = K \div \lambda 。 $$ (13) 式中：μ为滑动力合矢量T与合矢量平面局部坐标轴$ x $的夹角。
安全系数K为实际受力条件下在投影极值方向的极限平衡理论解，但此时耗散方向仍然存在不平衡耗散力Fts≠0，由其表达的极限状态接近程度系数λ＜1，当极限状态接近程度系数为1时不平衡耗散力Fts=0所对应的考虑滑面塑性发展安全储备矢量力共线极限平衡安全系数就是Kλ。
由式（9）可知总势能W与投影方向角θ无关，因此，滑动势能$ {W_\theta } $或耗散势能$ {W_{ \bot \theta }} $对投影方向θ求偏导的极值条件$ \frac{{\partial {W_\theta }}}{{\partial \theta }} = 0 $或$ \frac{{\partial {W_{ \bot \theta }}}}{{\partial \theta }} = 0 $，等价于势能比极值条件$ \frac{{\partial {W_\theta }/W}}{{\partial \theta }} = 0 $或$ \frac{{\partial {W_{ \bot \theta }}/W}}{{\partial \theta }} = 0 $，即势能极值条件方程式（10）等价于势能比极值条件方程。
共线平衡力系计算问题的投影方向和力矩中心位置可以任意选取计算结果不变，但非共线不平衡矢量力情况两者都不能随意选取。当使用力矩平衡条件时，力矩平衡中心点位置必须在不平衡耗散力作用线上，同时还必须远离滑动面处才行。由力矩平衡条件方程得到的安全系数K与投影极值方向解相同。力矩中心点位置极值原理是对投影方向极值原理的进一步补充和限定，是理论模型封闭完备的必要条件—阻断理论模型的非议漏洞。
极值条件不是笔者经验假设，是客观自身存在的规律表现，笔者仅仅是发现和揭示了事物内在本质联系特性，并利用这个普世规律研究滑动面应力矢量极限平衡抗滑稳定理论模型。
2. 理论模型计算步骤
因有限元应力场为单元节点的数值解，因此理论模型计算过程只能先依据三角形微面各角点的应力坐标分量计算微面合矢量力坐标分量，具体有10个步骤。
（1）采用现行可靠有限元软件（如ANSYS、ABAQUS等），确定考虑塑性影响的真实单元应力场。
（2）根据应力场计算结果，确定滑动面与其通过单元棱线的交点坐标值及单元应力场在该点的应力数值$ {\left\{\sigma \right\}}_\text{e}={\left\{{\sigma }_{x}，{\sigma }_{y}，{\sigma }_{z}，{\tau }_{xy}，{\tau }_{yz}，{\tau }_{zx}\right\}}_\text{e} $；由单元棱线各交点的坐标及其应力确定单元所有交点的形心点坐标及其应力数值，并计算形心点对应的滑面上落点坐标。根据有限元网格相邻单元关系确定各交点处相邻单元在交点的应力平均数值$ \left\{\sigma \right\}=\{{\sigma }_{x}，{\sigma }_{y}，{\sigma }_{z}，{\tau }_{xy}， $ $ {\tau }_{yz}，{\tau }_{zx}\} $。
（3）确定滑动面通过单元三角形微面的局部坐标系与整体坐标系的转换矩阵[B]：
$$ [\boldsymbol{B}] = \left[ {\begin{array}{*{20}{c}} {{l_1}}&{{l_2}}&{{l_3}} \\ {{m_1}}&{{m_2}}&{{m_3}} \\ {{n_1}}&{{n_2}}&{{n_3}} \end{array}} \right] 。 $$ (14) 式中：li，mi，ni分别为微面局部坐标轴$ x' $，$ y' $，$ z' $在整体坐标系xyz中的方向余弦，并且ni为微面的法线方向余弦。
（4）根据局部坐标应力与整体坐标应力转换关系$ {[\boldsymbol{\sigma} ]_{i'}} = {[\boldsymbol{B}]^\text{T}}[\boldsymbol{\sigma}][\boldsymbol{B}] $确定微面各角点在局部坐标系的法向应力$ {\sigma_{iz'}} $、切向应力分量$ {\tau_{ix'}} $和$ {\tau_{iy'}} $。其中
$$ [\boldsymbol{\sigma}] = \left[ {\begin{array}{*{20}{c}} {{\sigma_x}}&{{\tau_{xy}}}&{{\tau_{zx}}} \\ {{\tau_{xy}}}&{{\sigma_y}}&{{\tau_{yz}}} \\ {{\tau_{zx}}}&{{\tau_{yz}}}&{{\sigma_z}} \end{array}} \right] 。 $$ (15) （5）由单元棱线交点的形心点落点和单元的2个相邻棱线交点构成一个三角形微面，由切向应力$ {\tau_{ix'}} $和$ {\tau_{iy'}} $计算三角形微面的滑动力分量$ {T_{ix'}} $，$ {T_{iy'}} $；按每个点的滑动应力方向和法向应力$ {\sigma_{iz'}} $计算极限抗滑应力$ {s_i} = {f_i}{\sigma_{iz'}} + {c_i} $及抗滑应力坐标分量$ {s_{ix'}} $和$ {s_{iy'}} $，计算三角形微面局部坐标方向抗滑力分量$ {S_{ix'}} $和$ {S_{iy'}} $。分别把微面的滑动力、抗滑力局部坐标分量转换成整体坐标分量：$ {F_x} = {F_{ix'}}{l_1} + {F_{iy'}}{m_1} $，$ {F_y} = {F_{ix'}}{l_2} + {F_{iy'}}{m_2} $，$ {F_z} = {F_{ix'}}{l_3} + {F_{iy'}}{m_3} $。
（6）把各微面的滑动力、抗滑力整体坐标分量分别累加，得到滑动力总合矢量T的坐标分量（Tx，Ty，Tz）及抗滑力S总合矢量的坐标分量（Sx，Sy，Sz）。根据2个总合矢量的坐标分量计算合矢量所在平面的局部坐标-整体坐标转换矩阵，形式如同[B]。
（7）把各微面的滑动力、抗滑力整体坐标分量分别投影到合矢量平面内，得到各微面滑动力$ {T'_i} = {\tau '_i}\text{d}A $、抗滑力$ {S'_i} = {s'_i}\text{d}A $、相对于合矢量平面局部坐标轴x的倾角αi和εi。
（8）当滑动面上有需要释放的拉应力时，计算方法过程同上。
（9）检查各微面的滑动力是否有与滑动力合矢量反向的，若有把该微面上的抗滑力方向进行反向（微面上滑动力与滑动力合矢量反向时，该微面的抗滑力与其同向）。全部检查操作完成后重新计算抗滑力合矢量坐标分量及合矢量所在平面的局部坐标-整体坐标转换矩阵，返回到第（7）→（10）步。若无则进行第（10）步。
（10）在合矢量平面内，根据计算理论模型确定安全系数K、投影极值方向角θ、极限状态接近程度系数λ、滑动面应力矢量安全度系数Kλ。
计算程序主框图见图 4。
3. 三维典型算例理论模型解
3.1 椭球滑动面边坡Zhang-X算例[6]
算例坡高12.2 m，坡比1︰2，包括无软夹层椭球滑动面和坡脚以下1.45 m处0.1 m厚水平软夹层组合椭球面2种情况，土层物理力学参数如表 1。
表 1 椭球滑动面土层力学参数Table 1. Mechanical parameters of soil on ellipsoidal surface参数 内摩擦角φ/(°) 黏聚力c/(kN·m-2) 重度/(kN·m-3) 弹性模量E/kPa 泊松比μ 剪胀角/(°) 土层 20 28.74 18.84 5.0×104 0.30 4 软层 10 0. 18.84 2.0×103 0.25 2 有限元应力场计算采用ANSYS软件，塑性屈服准则按主应力空间莫尔库仑强度条件的六棱锥等底面积圆锥Druker-Prager塑性屈服条件计算[18]。抗滑稳定影响区域Solid45六面体单元边长最大不超过2 m，自重分5个荷载增量子步施加。有限元模型坡面倾向的水平方向为x坐标轴方向范围0~50 m，坡面走向为z坐标轴方向0~150 m，高度方向为y坐标0~22.2 m。对称面处坡肩点坐标（14.6，22.2，75）、坡脚点（39，10，75）。有软夹层情况共计剖分单元34650个，单元剖分网格及水平位移云图见图 5。
![]()
图 5 有限元网格及水平位移云图Figure 5. Finite element meshes and nephogram of horizontal displacement椭球滑动面中心点坐标(32.9, 31.3, 75), 短轴24.4 m，长轴66.9 m，椭球面方程式为
$$ \frac{{{{(x - 32.9)}^2}}}{{{{24.4}^2}}} + \frac{{{{(y - 31.3)}^2}}}{{{{24.4}^2}}} + \frac{{{{(z - 75)}^2}}}{{{{66.9}^2}}} = 1 。 $$ 两种计算情况对称面处等效塑性应变云图及滑动面如图 6所示，理论模型解及相关文献计算结果见表 2，3。
![]()
图 6 椭球滑动面对称面处无软夹层情况、有软夹层情况等效塑性应变云图Figure 6. Equivalent plastic strain nephogram of ellipsoidal sliding surface without and with soft strata on symmetric plane表 2 椭球滑动面两种情况抗滑稳定计算结果Table 2. Calculated results of stability of ellipsoidal sliding surface in two cases计算方法 情况1无软夹层 情况2有软夹层 安全系数 与本文Kλ相比/% 安全系数 与本文Kλ相比/% Zhang-X极限平衡法[6] 2.122 -2.08 1.553 -14.13 Hungr简化Bishop法[6] 2.167 +0.04 1.620 -9.41 陈祖煜极限平衡法[2-3] 2.187 +0.96 1.640 -8.08 郑宏滑面极限平衡法[7-8] 2.140 -1.22 1.706 -3.89 Lam通用极限平衡法[9] — 1.603 -10.57 Huang极限平衡法[10] — 1.665 -6.46 Jiang主滑向极限平衡法[11] 2.127 -1.84 1.766 -0.37 极限分析上限法[3, 12] 2.262 +4.24 1.717 -3.23 应力数值解代数和法[13] — 1.607 -10.30 有限元应力主滑向极限平衡法[14] 2.226 +2.69 1.730 -2.46 表 3 矢量理论模型计算结果Table 3. Calculation results of vector theoretical model参数 本文矢量理论解 滑动力合矢量倾角和合矢量夹角δ 投影极值方向倾角和投影极值方向角θ 应力代数和法安全系数 滑面通过单元数和三角形微面数 情况1无软夹层 Kλ=2.1661，K=2.0546，λ=0.9485 倾角24.5607°，夹角δ=2.8630° 倾角32.8839°，θ=-8.3232° 代数和法K=2.1006，与本文Kλ相比为-3.12% 单元1606个，微面6420个 情况2有软夹层 Kλ=1.7725，K=1.7363，λ=0.9796 倾角25.8634°，夹角δ=1.1573° 倾角20.8109°，θ=-5.0528° 代数和法K=1.7404，与本文Kλ相比为-1.84% 单元1684个，微面6776个 由表 2，3计算结果可以看出，三维极限平衡法与理论模型解偏差为-14.13%~+4.24%。无软夹层情况1的均质边坡极限平衡法计算结果与理论模型解偏差较小-2.08%~+4.24%，但有软夹层的情况2时都是小于理论模型解偏于安全的估计值，如Zhang-X极限平衡法偏差-14.13%，并且同一方法两种情况偏差变幅最大。
矢量理论模型安全度系数表达两种滑面情况的绝对安全度Kλ为2.1661，1.7725级差为0.3936；极限状态接近程度系数表达滑面距离完全塑性临界共线极限平衡状态的相对安全度1/$ \lambda $为1/0.9458，1/0.9796级差为0.0345。
由图 6可知，无软夹层情况的椭球滑动面通过区域的等效塑性应变小于38×10-6；有软夹层情况的椭球面通过区域等效塑性应变最大值150×10-6，软夹层内为35×10-3~75×10-3。有软夹层滑动面塑性程度远大于无软夹层，实际应力状况更接近临界共线极限平衡状态，滑动可能性大。
有软夹层情况，按滑动面20个单元面的应力合成矢量为一组数据，采用CAD精准绘制的合矢量平面内矢量平衡关系如图 7。
随投影方向角变化的安全系数、势能、势能比曲线见图 8。
![]()
图 8 随投影方向角变化的安全系数、势能、势能比曲线Figure 8. Curves of factor of safety, potential energy and potential energy ratio with respect to projection direction矢量平衡关系图 7中合矢量T与S的夹角与表 3计算结果一致，滑动面的滑动力矢量与抗滑力矢量为非共线不平衡矢量；垂直于耗散力的投影极值方向与滑动力合矢量方向（合矢量平面局部坐标轴x）夹角，与表 3计算值θ=-5.0528°一致。
图 8中随投影方向角θ增大变化，安全系数式（6）按直线规律增大（深绿色），安全系数式（11）按曲线快速减小（浅绿色），两式曲线交点为理论解。耗散势能曲线（粉色）、耗散势能比曲线（橘红色）、变形势能比曲线（紫色）的极值点与理论解投影方向角θ=-5.0528°一致。而变形能曲线（深蓝色）受安全系数变化影响敏感，沿投影方向轴的极值点位置与理论解不对应，但变形势能比曲线减小了安全系数的变化影响，精准地反应了极值条件方程物理意义。如果变形能曲线以总势能线（青色）为方向轴时的切点位置与投影极值方向角理论解对应。
3.2 楔体抗滑稳定算例理论解
楔体抗滑稳定算例1选自文献[2]、算例2选自文献[14]。几何形体及物理力学参数指标见表 4，5所示。
表 4 楔体几何形体特征Table 4. Geometric features of wedges部位 楔体算例1 楔体算例2 倾向/(°) 倾角/(°) 倾向/(°) 倾角/(°) 左结构面 115 45 120 40 右结构面 245 45 240 60 顶面 180 10 180 0 坡面 180 60 180 60 楔体坡高/m 64.89 98.4 表 5 模型物理力学参数Table 5. Physical and mechanical parameters of models算例工况 内摩擦角φ/(°) 黏聚力c/(kN·m-2) 重度/(kN·m-3) 弹性模量E/kPa 泊松比$ \nu $ 剪胀角/(°) 楔体1对称工况左右结构面 20 50 25.48 1.0×105 0.25 4 弹性岩体 25.48 5.0×106 0.20 楔体1非对称工况左结构面 20 50 25.48 1.0×105 0.25 4 右结构面 10 30 25.48 1.0×105 0.25 3 楔体2左右结构面相同 30 50 26.0 1.0×105 0.25 4 弹性岩体 26.0 5.0×106 0.20 有限元应力场计算采用ANSYS软件，结构面按厚度0.4 m的薄层六面体单元Solid45，单元网格剖分最大尺寸3 m，结构面塑性屈服按主应力空间莫尔库仑强度条件六棱锥等底面积圆锥理想弹塑性D-P准则。岩体为Solid45四面体单元，滑动体单元边长最大不超过5 m，滑动面下部岩体最大边长10 m。自重分5个荷载增量子步施加。有限元模型坡面倾向水平为x坐标方向，坡面走向为z坐标轴方向，高度方向为y坐标。算例1有限元模型坐标范围x=0~250 m、y=5~110 m、z=0~145 m共计剖分419858个单元，其中结构面单元11616个，如图 9；算例2坐标范围x=0~185 m、y=30~150 m、z=17~135 m剖分342623个单元，其中结构面单元6936个。
![]()
图 9 楔体算例1与算例2有限元剖分网格图Figure 9. Finite element meshes of wedges (example 1 and example 2)结构面单元的等效塑性应变云图见图 10。理论模型解及相关文献计算结果汇总见表 6，7。
![]()
图 10 楔体1对称材料、楔体1非对称材料、楔体2结构面单元的等效塑性应变云俯视图Figure 10. Top views of equivalent plastic strain nephogram for wedge 1 (symmetric material), wedge 1 (asymmetric material) and wedge 2 in soft rock elements表 6 楔体抗滑稳定算例计算结果Table 6. Calculated results of anti-sliding stability of wedges计算方法 楔体1材料对称 楔体1材料非对称 楔体2交线倾向11.3380° 安全系数 投影倾角/(°) 与本文比较/% 安全系数 投影倾角/(°) 与本文比较/% 安全系数 投影倾角/(°) 与本文比较/% 楔体极限平衡[15, 16] 1.5578 22.9098 +10.17 1.1815 22.9098 -7.91 1.6400 28.9964 +4.30 陈祖煜简化极限平衡法[2] 1.556 22.9098 +10.06 1.167 22.9098 -9.25 — 郑宏滑面严格平衡法[8] — 1.636 +4.07 极限平衡法[11] — 1.629 +3.66 主滑向应力代数和[14] — 1.497 -4.84 表 7 矢量理论模型计算结果Table 7. Calculation results of vector theoretical model参数 楔体1材料对称 楔体1材料非对称 楔体2交线倾向11.3380° 安全系数 投影倾角 安全系数 投影倾角 安全系数 投影倾角 本文矢量理论模型解 Kλ=1.3994 40.2520° Kλ=1.2750 28.7662° Kλ=1.5694 42.2410° 安全系数K和程度系数λ K=1.3365 λ=0.9551 K=1.2384 λ=0.9713 K=1.5034 λ=0.9580 合矢量夹角δ和投影倾向 δ=2.5649° 倾向0° δ=1.6438° 倾向1.6267° δ=2.3986° 倾向16.7218° 滑动力矢量倾角和倾向 倾角39.7084° 倾向0° 倾角29.1022° 倾向3.1532° 倾角44.0963° 倾向10.4716° 代数和法安全系数K K=1.3392 -4.50% K=1.2375 -3.03% K=1.4978 -4.78% 滑面通过单元和微面数 单元11616个 微面46493个 单元6936个 微面28014个 由有限元等效塑性应变云图 10可知楔体结构面最大塑性应变都位于楔体的坡面处最低点，楔体1对称情况塑性应变值为0.015909最小、非对称情况为0.025739最大、楔体2为0.023369居中，塑性应变值的大小表现了结构面的塑性发展程度。塑性程度排序与表 7中极限状态程度系数的排序顺次一样，楔体1对称情况λ=0.9551相对安全度1/λ最大、非对称情况λ=0.9713相对安全度1/λ最小、楔体2的λ=0.9580居中。理论模型的极限状态接近程度系数与滑动面的塑性程度具有关联性。
由楔体抗滑稳定计算结果表 6，7可以看出，三维极限平衡法与理论模型解的偏差-9.25%~+10.17%。楔体1结构面材料力学参数对称时大于理论模型解10%，而非对称时又小于理论模型解-7.91%~-9.29%，并且同一方法两种情况偏差最大变幅达到19.3%。楔体2结构面几何特性非对称力学参数对称情况的极限平衡法与理论模型解偏差-4.84%~+4.30%。
楔体算例2的合矢量平面内，按50个单元面合矢量为一组数据精准绘制的矢量平衡关系及投影极值方向见图 11。合矢量平面内随投影方向角变化的安全系数、势能、势能比曲线见图 12。两图的解读分析与图 7，8相同不做赘述。
![]()
图 11 楔体算例2合矢量平面内的矢量平衡图Figure 11. Vector equilibrium of wedge 2 in resultant vector plane![]()
图 12 楔体算例2随投影方向变化的安全系数、势能曲线Figure 12. Variation of curves of factor of safety and potential energy of wedge 2 with projection direction岩质楔体极限平衡稳定计算方法[15-16]，是以结构面交线方向为极限平衡计算方向，即假设结构面的滑动力和抗滑力方向是结构面交线方向，这种简化假设当楔体处于临界极限平衡状态是合理的（即安全度系数为1时具有合理性）；当安全度系数不为1时的非临界极限状态时，滑动面上的剪切应力方向与结构面交线方向不同，必然导致了楔体极限平衡方法与理论模型解存在较大偏差，而且安全系数越远离1时假设计算方向影响的偏差就越大。楔体极限平衡抗滑稳定计算方法在安全系数非1时的计算结果是用临界极限平衡时的滑动方向研究一般极限平衡问题，具有极限平衡基本概念认知混淆固有偏差。特别是很多文献方法都以该方法的计算结果为评价标准，相关问题请岩土工程业内学者积极关注。
滑面抗滑稳定极限平衡的安全度系数可以不为1，而临界极限平衡的安全度系数为1（即将滑动时刻），前者是一般极限平衡力学概念且包含后者特定情况，但绝不仅是临界极限平衡特定情况。两者的应力状态、力学基本概念具有本质区别，绝不能混淆基本概念用特例取代一般。土力学研究存在的诸多错误问题[19-21]根源都与该力学基本概念的混淆有关。
理论模型不受滑面交线方向影响，采用每个计算点处的实际滑动剪切应力方向进行计算，没有人为简化假设限定和相应的适用条件制约。
4. 应力场参数变化对理论解的敏感性
（1）结构面单元剖分尺寸
采用非对称岩质楔体边坡算例2进行分析，对结构面剖分时控制最大尺寸分别采用3 m四面体单元和3，4，5 m的六面体单元进行有限元应力场计算，理论模型解如表 8所示。
表 8 单元剖分尺寸对理论模型解的影响Table 8. Influences of element size on solution of theoretical model单元类型及尺寸 3 m四面体 3 m六面体 4 m六面体 5 m六面体 安全度系数Kλ 1.5711（+0.108%） 1.5694（0.0%） 1.5703（+0.057%） 1.5712（+0.115%） 安全系数K 1.5055 1.5034 1.5045 1.5054 程度系数λ 0.9582 0.9580 0.9581 0.9581 三角形微面数 28014 16417 10742 结构面单元数 13673 6936 4056 2646 单元总数 135732 342623 190545 127207 注：括号内数字为相对于3 m六面体单元的安全度系数Kλ变化率。 结构面有限元剖分尺寸对理论模型安全度系数的最大影响仅0.115%。
（2）弹性模量变化影响
采用非对称岩质楔体边坡的5 m六面体剖分网格，分别把结构面弹性模量1×105 kPa、岩体弹性模量5×106 kPa乘以0.5，0.75，1.25，1.5后计算有限元应力场，得到矢量理论模型抗滑稳定安全度系数计算结果如表 9。
表 9 弹性模量变化对理论模型解的影响Table 9. Influences of change of elastic modulus on solution of theoretical model变化部位 弹性模量乘数 0.50 0.75 1.00 1.25 1.50 仅结构面 安全度系数Kλ 1.6042（+2.06%） 1.5858（+0.92%） 1.5712 1.5592（-0.77%） 1.5497（-1.39%） 程度系数λ 0.9623 0.9599 0.9581 0.9566 0.9553 仅岩体 安全度系数Kλ 1.5330（-2.49%） 1.5555（-1.01%） 1.5712 1.5825（+0.71%） 1.5920（+1.31%） 程度系数λ 0.9535 0.9562 0.9581 0.9595 0.9608 结构面
及岩体安全度系数Kλ 1.5711 1.5711 1.5712 1.5711 1.5712 程度系数λ 0.9581 0.9581 0.9581 0.9581 0.9584 注：括号内数字为相对于弹性模量乘数1.00的安全度系数Kλ变化率。 结构面弹模增大、或岩体弹模减小时，滑动面的抗滑力矢量减小幅度大于滑动力，矢量理论安全度系数减小。而结构面与岩体的弹模按相同比例变化时，理论安全度系数几乎不变。
（3）泊松比变化影响
分别把结构面泊松比μ=0.25、岩体泊松比μ=0.20按0.05级差进行变化确定有限元应力场，理论模型安全度系数见表 10。
表 10 泊松比变化对理论模型的影响Table 10. Influences of Poisson's change on solution of theoretical model部位 泊松比μ μ=0.15 0.20 0.25 0.30 0.350 结构面 安全度系数Kλ 1.5849（+0.86%） 1.5787（+0.75%） 1.5712（0.0%） 1.5619（-0.60%） 1.4827（-1.29%） 程度系数λ 0.9601 0.9593 0.9581 0.9571 0.9559 岩体 安全度系数Kλ 1.5577（-0.87%） 1.5712（0.0%） 1.5856（+0.91%） 1.6016（+1.90%） 1.6171（+2.84%） 程度系数λ 0.9570 0.9581 0.9593 0.9609 0.9622 两者同比例 泊松比乘数 0.75μ 0.8μ 1.0μ 1.2μ 1.4μ 安全度系数Kλ 1.5685（-0.17%） 1.5689（-0.15%） 1.5712（0.0%） 1.5750（+0.24%） 1.5805（+0.59%） 程度系数λ 0.9584 0.9584 0.9581 0.9581 0.9583 注：括号内数字为相对于结构面泊松比μ=0.25、岩体泊松比μ=0.20的安全度系数Kλ变化率。 结构面泊松比增大、或岩体泊松比减小时滑动面的理论安全度系数减小。而结构面与岩体的泊松比按相同比例变化时，理论安全度系数随泊松比变化比例数值的增大而增大，变化幅度不大于0.6%。
（4）塑性流动法则参数剪胀角的变化影响
在结构面D-P模型的塑性本构流动法则中，塑性屈服应变参数剪胀角的变化对有限元应力场的影响，按照剪胀角为0°至8°级差2°变化，理论安全度系数如表 11。
表 11 剪胀角变化对理论模型解影响Table 11. Influences of change of dilatancy angle on theoretical model剪胀角 0° 2° 4° 6° 8° 安全度系数Kλ 1.5800（+0.56%） 1.5753（+0.26%） 1.5712（0.0%） 1.5681（-0.20%） 1.5646（-0.42%） 安全系数K 1.5207 1.5126 1.5054 1.4997 1.4935 程度系数λ 0.9625 0.9602 0.9581 0.9564 0.9546 注：括号内数字为相对于剪胀角4°的安全度系数Kλ变化率。 随结构面的塑性剪胀角增大，理论模型安全度系数减小，且安全系数和极限状态接近程度系数也减小。
（5）结构面的塑性准则影响
分别采用主应力空间六棱锥等底面积圆锥D-P塑性准则、六棱锥外接圆锥D-P塑性准则、弹性模型计算应力场，相应的理论安全度系数计算结果见表 12。
表 12 不同塑性准则对理论模型解的影响Table 12. Influences of different plastic criteria on solution of theoretical model塑性准则 等底面积D-P 外接圆D-P 弹性 安全度系数Kλ 1.5712 1.4200（-10.64%） 1.2960（-21.23%） 安全系数K 1.5054 1.3315 1.1801 程度系数λ 0.9581 0.9397 0.9106 滑动力合矢量/(104 kN) 292.525 306.053（+4.42%） 321.222（+8.93%） 滑动力合矢量倾角/(°) 44.1091 48.5751 51.6391 抗滑力合矢量/(104 kN) 439.202 405.907（-8.20%） 379.014（-10.63%） 合矢量夹角/(°) 2.3900 3.5533 5.1188 极值方向倾角/(°)
极值方向倾向/(°)42.3137
16.560146.8366
20.854551.4313
19.4447注：括号内数字为相对于等底面积D-P准则的变化率。 结构面的塑性准则条件由等底面积D-P准则→外接圆D-P准则→弹性的变化过程，实质是对结构面塑性屈服条件的放宽过程，即滑动面上实际剪应力虚假增大、而法向压应力虚假减小过程。剪应力的增大导致了滑动力合矢量的增大，而法向压应力的减小产生抗滑力的减小，因此使安全度系数显著减小。
笔者除探讨以上应力场影响因素规律外，还研究了结构面的抗剪断力学参数内摩擦角φ和黏聚力c逐步减小情况，得到规律：合矢量夹角δ、安全度系数Kλ相应减小（逐步接近临界极限状态）；极限状态接近程度系数λ逐步变大接近1.0；合矢量夹角δ逐步接近0，滑动力合矢量的倾角和投影极值方向的倾角逐渐趋近于结构面交线的倾角（接近临界极限平衡状态）。但是，随抗剪断力学参数的减小，等底面积D-P塑性准则有限元应力场的塑性迭代计算收敛逐步变得困难（迭代次数增大机时增大），矢量理论模型安全度系数Kλ和安全系数K也逐步接近；滑动力合矢量倾角和投影极值方向倾角由陡于结构面交线倾角最终逐步变成了缓于结构面交线倾角。因此，本文认为对于三维岩质楔体边坡的抗滑稳定计算当采用等底面积D-P塑性准则时，矢量理论模型安全度系数的临界值是投影极值方向倾角不缓于结构面交线倾角的对应数值。对于等底面积D-P塑性准则应力计算无法收敛的楔体边坡（或计算收敛困难时），可采用外接圆D-P塑性准则应力场计算成果，此时矢量理论模型安全度系数的临界值建议以1.05为限。
5. 结论
三维滑动面应力矢量极限平衡抗滑稳定理论模型的建立和计算结果是岩土工程创新性基础力学理论研究成果。滑动面应力矢量势能极值特征曲线的存在性和极值点与计算结果的吻合关系是理论模型投影方向极值原理合理可靠性的具体证明。
（1）理论模型的三维拓展考虑每个计算点的应力状态差异性，使三角形微面切向滑动力与抗滑力的作用方向在微面内存在夹角是客观存在的合理表达。
（2）理论模型没有人为简化假设限定和适用条件制约，既适用于土质边坡抗滑稳定计算，也适用于岩质边坡抗滑稳定计算。
（3）理论模型的滑动面矢量抗滑稳定安全度系数Kλ=K÷λ，由投影极值方向的极限平衡安全系数K和耗散方向的极限平衡状态接近程度系数λ确定。三维理论模型解是只依赖于空间有限元单元应力场每个计算点的应力状态数值，通过2个极值条件方程理论定解。
（4）目前已有的三维极限平衡计算方法，在均质土边坡的椭球滑动面问题与理论模型解的偏差较小-2.0%~+4.2%，而有软夹层情况的土质边坡或岩质楔体边坡稳定计算结果与理论模型偏差较大-14%~+10%。相关问题应该引起岩土工程抗滑稳定研究学者的积极关注。
（5）有限元应力场的影响因素（单元剖分尺寸、弹性模量、泊松比、剪胀角、塑性准则）敏感性计算表明，三维岩质楔体稳定问题塑性准则的影响最大。
-
![]()
图 3 邓肯E-B模型计算的竣工期堆石坝体位移与应力分布等值线
Figure 3. Contours of displacements and stresses of rockfill zones predicted by E-B model
![]()
图 4 “南水”模型计算的竣工期堆石坝体位移与应力分布等值线
Figure 4. Contours of displacements and stresses of rockfill zones predicted by Nanshui model
![]()
图 5 竣工时坝体内主应力作用面（E-B模型）
Figure 5. Principal stress planes after construction (E-B model)
![]()
图 6 堆石料单向压缩试验结果与两种模型模拟结果
Figure 6. Oedometric compression results of rockfill materials and predictions by two models
![]()
图 8 邓肯E-B模型计算的蓄水期面板位移与应力分布等值线
Figure 8. Contours of displacements and stresses of concrete slabs predicted by E-B model
![]()
图 9 “南水”模型计算的蓄水期面板位移与应力分布等值线
Figure 9. Contours of displacements and stresses of concrete slabs predicted by Nanshui model
![]()
图 10 邓肯E-B模型计算的蓄水引起的垫层表面位移增量等值线
Figure 10. Contours of impounding-induced displacements of cushion layer predicted by E-B model
![]()
图 11 “南水”模型计算的蓄水引起的垫层表面位移增量等值线
Figure 11. Contours of impounding-induced displacements of cushion layer predicted by Nanshui model
![]()
图 12 邓肯E-B模型计算的蓄水引起的堆石坝体位移增量等值线
Figure 12. Contours of impounding-induced displacements of rockfill zones predicted by E-B model
![]()
图 13 “南水”模型计算的蓄水引起的堆石坝体位移增量等值线
Figure 13. Contours of impounding-induced displacements of rockfill zones predicted by Nanshui model
![]()
图 20 “南水”模型计算的典型单元应变与应力发展过程
Figure 20. Evolution of strain and stress of element predicted by.Nanshui model
![]()
图 21 “南水”模型塑性柔度矩阵非对角元变化规律
Figure 21. Off-diagonal elements of plastic flexibility matrix of Nanshui model
表 1 筑坝堆石料本构模型参数
Table 1 Constitutive model parameters of damming rockfill materials
材料 ρ/ (g·cm-3) φ0/
(°)Δφ/
(°)Rf k n kb m cd/
%nd Rd 堆石料A 2.20 54.3 9.8 0.61 1230 0.26 611 0.12 0.28 0.77 0.55 堆石料B 2.20 52.3 8.0 0.61 892 0.33 351 0.27 0.52 0.59 0.57 
下载: 导出CSV
-
[1] 混凝土面板堆石坝设计规范: NB/T 10871—2021E[S]. 北京: 中国水利水电出版社, 2022. Code for Design of Concrete Face Rockfill Dams: NB/T 10871—2021E[S]. Beijing: China Water & Power Press, 2022. (in Chinese)
[2] 程展林, 姜景山, 丁红顺, 等. 粗粒土非线性剪胀模型研究[J]. 岩土工程学报, 2010, 32(3): 460-467. CHENG Zhanlin, JIANG Jingshan, DING Hongshun, et al. Nonlinear dilatancy model for coarse-grained soils[J]. Chinese Journal of Geotechnical Engineering, 2010, 32(3): 460-467. (in Chinese)
[3] XIAO Y, LIU H L, YANG G. A constitutive model for the state-dependent behaviors of rockfill material considering particle breakage[J]. Science China (Technological Sciences), 2014, 57(8): 1636-1646. doi: 10.1007/s11431-014-5601-6
[4] FU Z Z, CHEN S S, LIU S H. Hypoplastic constitutive modelling of the wetting induced creep of rockfill materials[J]. Science China (Technological Sciences), 2012, 55(7): 2066-2082. doi: 10.1007/s11431-012-4835-4
[5] 孔宪京, 邹德高. 紫坪铺面板堆石坝震害分析与数值模拟[M]. 北京: 科学出版社, 2014. KONG Xianjing, ZOU Degao. Seismic Damage Analysis and Numerical Simulation of Zipingpu Concrete Face Rockfill Dam[M]. Beijing: Science Press, 2014. (in Chinese)
[6] DUNCAN J M, CHANG C Y. Nonlinear analysis of stress and strain in soils[J]. Journal of Soil Mechanics and Foundations Division, ASCE, 1970, 96(SM5): 1629-1653.
[7] 沈珠江. 理论土力学[M]. 北京: 中国水利水电出版社, 2000. SHEN Zhujiang. Theoretical Soil Mechanics[M]. Beijing: China Water & Power Press, 2000. (in Chinese)
[8] 土工试验方法标准: GB/T 50123—2019[S]. 北京: 中国计划出版社, 2019. Standard for geotechnical testing method: GB/T 50123—2019[S]. Beijing: China Planning Press, 2019. (in Chinese)
[9] 殷宗泽. 土工原理[M]. 北京: 中国水利水电出版社, 2007. YIN Zongze. Geotechnical Principle[M]. Beijing: China Water & Power Press, 2007. (in Chinese)
[10] FU Z Z, CHEN S S, PENG C. Modeling cyclic behavior of rockfill materials in a framework of generalized plasticity[J]. International Journal of Geomechanics, ASCE, 2014, 14(2): 191-204. doi: 10.1061/(ASCE)GM.1943-5622.0000302
[11] 方维凤. 混凝土面板堆石坝流变研究[D]. 南京: 河海大学, 2003. FANG Weifeng. Study on Rheology of Concrete Face Rockfill Dam[D]. Nanjing: Hohai University, 2003. (in Chinese)
[12] 王永明, 汝璇卿. 高面板堆石坝工作性态的弹性和弹塑性对比分析[J]. 三峡大学学报(自然科学版), 2009, 31(6): 24-28. doi: 10.3969/j.issn.1672-948X.2009.06.006 WANG Yongming, RU Xuanqing. Comparative analysis of high concrete faced rockfill dam with elastic analysis and elastoplastic analysis methods[J]. Journal of China Three Gorges University (Natural Sciences), 2009, 31(6): 24-28. (in Chinese) doi: 10.3969/j.issn.1672-948X.2009.06.006
[13] 顾淦臣, 束一鸣, 沈长松. 土石坝工程经验与创新[M]. 北京: 中国电力出版社, 2004. GU Ganchen, SHU Yiming, SHEN Changsong. Experience and Innovation of Earth-Rock Dam Engineering[M]. Beijing: China Electric Power Press, 2004. (in Chinese)
[14] 殷宗泽, 张坤勇, 朱俊高. 面板堆石坝应力变形计算中考虑土的各向异性[J]. 水利学报, 2004, 35(11): 22-26, 32. doi: 10.3321/j.issn:0559-9350.2004.11.004 YIN Zongze, ZHANG Kunyong, ZHU Jungao. Computation for stress and deformation of concrete slab inrock-fill dam in consideration of soil anisotropy[J]. Journal of Hydraulic Engineering, 2004, 35(11): 22-26, 32. (in Chinese) doi: 10.3321/j.issn:0559-9350.2004.11.004
[15] 殷宗泽, 徐志伟. 土体的各向异性及近似模拟[J]. 岩土工程学报, 2002, 24(5): 547-551. doi: 10.3321/j.issn:1000-4548.2002.05.001 YIN Zongze, XU Zhiwei. Anisotropy of soils and its approximate simulation[J]. Chinese Journal of Geotechnical Engineering, 2002, 24(5): 547-551. (in Chinese) doi: 10.3321/j.issn:1000-4548.2002.05.001
-
其他相关附件
计量
- 文章访问数: 457
- HTML全文浏览量: 51
- PDF下载量: 104
