Local large deformation between cut-off wall and core wall on deep overburden by meshless method
-
摘要: 土质心墙坝作为覆盖层地基上的一种重要坝型，其上部土质心墙与下部混凝土防渗墙直接相连，由于材料刚度差异悬殊，两者在接头区域往往存在防渗墙“贯入”心墙的局部大变形现象。在任意拉格朗日-欧拉(ALE)框架下，开发了无网格大变形弹塑性分析方法，协同发挥了ALE框架在大变形分析中稳定和高精度的优势以及无网格法(Meshless)无需单元拓扑的特点，解决了有限元法因网格扭曲引发的精度降低或复杂重剖分问题。基于面向对象的C++语言，将该方法集成到自主开发的GEODYNA计算系统中，实现无网格大变形方法与有限元(FEM)、比例边界有限元(SBFEM)的无缝耦合。将耦合Meshless-FEM-SBFEM方法应用于某深厚覆盖层上土质心墙坝，联合土体广义弹塑性模型，再现了坝体在填筑过程中防渗墙“贯入”心墙的局部大变形发展过程，发现在竣工期小变形分析低估防渗墙竖向应力约4.1 MPa(占13%)，并指出了防渗墙顶部两侧土体间存在剪切带。Abstract: As a significant type of dam on deep overburden, the earth core rockfill dam connects the concrete cut-off wall and soil core wall directly, and the "penetration" phenomenon with local large deformation may exist between cut-off wall and core wall due to the difference in material stiffness. Using an arbitrary Lagrangian-Euler (ALE) framework, an elastic-plastic meshless large deformation method is developed in this research. The meshless method has the benefit of flexible nodal distribution, and the ALE framework shows the advantage of accuracy and stability in large deformation analysis, thus avoiding the precision reduction or re-meshing procedure in the mesh-based large deformation method. The introduced approach is incorporated into the self-developed calculation platform GEODYNA and combined with the finite element method (FEM) and the scale boundary finite element method (SBFEM). Finally, the coupled meshless-FEM-SBFEM applied to an earth core rockfill dam on deep overburden, and combined with the generalized elasto-plastic model, the "penetration" phenomenon with local large deformation is simulated. The results indicate that the large deformation analysis can capture the stress distribution of the cut-off wall and the soil deformation near the joint zone more reasonably, the vertical stress of the cut-off wall calculated by the small deformation analysis is underestimated by about 4.1 MPa (13%), and a shear zone exists between the soil area at the top of the cut-off wall.
-
0. 引言
随着中国西部能源的逐步开发，在以深厚覆盖层为代表的不良地质条件下建设高土石坝已是“难以避让”的问题。而挖除覆盖层不但对地质环境有着较大的影响，也将大幅度的延长建坝工期，相比之下直接将高土石坝建在深厚覆盖层上成为更符合实际施工情况的选择[1]。土质心墙坝对深厚覆盖层有着较好的适应性，近几十年来在国内取得了长足的发展并积累了大量的经验[2-4]。由于上部土质心墙与下部混凝土防渗墙直接相连，两者材料刚度相差悬殊，在填筑过程中连接处往往会发生混凝土防渗墙“贯入”心墙的局部大变形现象。这无疑给深厚覆盖层上土质心墙坝的数值模拟[4-5]和安全评价带来新的挑战。
大变形分析是岩土工程中常见但最具挑战性的课题之一[6-8]。早期基于网格的大变形分析方法主要是在单一的拉格朗日框架下实现的，包括完全拉格朗日(TL)[9]和更新拉格朗日(UL)[10]，上述两种方法在单元发生严重扭曲后，均出现精度急剧降低甚至出现求解失败现象。随着数值模拟方法的不断完善，TL和UL逐步被任意拉格朗日-欧拉(ALE)框架取代，其中最具代表性且在岩土工程中应用广泛的是RITSS[11]和CEL[12]。RITSS法需要不断网格重剖分及复杂的人机交互过程，无疑大幅度增加计算成本；CEL法内置于商用Abaqus中，无自主知识产权，难以实现深度二次开发。在基于网格的大变形模拟方法发展的同时，另一类基于节点的大变形方法也取得了长足的进步。其中最具代表性的是由Lucy在1977年提出的无网格法，随后广泛地应用于岩土工程问题中[13-17]。该方法仅需节点离散几何模型，无需网格的拓扑结构，在本质上回避了网格扭曲的问题，其在大变形模拟中有着天然的优势。然而相比于有限元法，无网格法计算效率较低，一直无法应用到大规模岩土工程模拟中。
发生大变形区域往往伴随着应力梯度的突变，因此对潜在大变形区建立精细化模型是取得精准、合理模拟结果的前提。研究表明[6-8]，岩土工程中的大变形现象往往集中在局部区域。以深厚覆盖层上土质心墙坝为例，其潜在大变形区域仅为防渗墙-心墙接头局部[2]，其余区域坝体单元应变较小。若整个坝体均采用无网格大变形分析无疑放大了该方法效率较低的劣势，增加了不必要的计算成本。因此可仅对潜在大变形局部区域采用无网格大变形分析，其余坝体部分采用FEM-比例边界有限元(SBFEM)跨尺度分析，协同发挥Meshless在大变形、FEM在计算效率及SBFEM[18]在跨尺度建模方面的优势，进而高效地开展深厚覆盖层上土质心墙坝防渗墙-心墙接头局部大变形分析。
本文协同发挥了ALE框架在大变形分析中高精度、稳定的优势及无网格方法无需单元拓扑的特点，开发了弹塑性无网格大变形分析方法。采用面向对象的C++语言，设计构造了超单元类数据结构，将该方法集成到自主开发的大型FEM计算系统GEODYNA中，建立了Meshless-FEM-SBFEM多数值耦合分析方法，提高了大变形分析效率。将该方法应用于某深厚覆盖层上土质心墙坝填筑期的模拟中。首先建立FEM-SBFEM土质心墙坝跨尺度模型，通过小变形分析确定潜在大变形区域。进一步在整体模型网格保持不变的基础上，采用无网格节点替换该区域网格开展大变形分析。最终再现了坝体填筑过程中防渗墙“贯入”心墙的局部大变形发展过程，指出了防渗墙顶部两侧土体间存在剪切带，发现了小变形分析低估防渗墙应力。
1. 无网格大变形分析方法理论推导
无网格法的核心思想是通过一系列的节点来离散和求解控制方程，无需网格的拓扑结构，有效地回避了网格扭曲的问题。根据插值函数、积分方式不同，无网格分为多种[13-17]，本文选取的基于径向基函数(RBF)插值及全局背景网格积分的RPIM法，在边界处可实现与SBFEM, FEM无缝耦合。每一个计算步由拉格朗日步和欧拉步组成。通过在拉格朗日步求解的节点位移，计算高斯点在变形过程中位置的改变，随后在每个欧拉步生成新高斯点，采用RPIM实现应力、应变及本构内变量在新旧高斯点间的传递，进而建立弹塑性无网格大变形分析方法。
在每个拉格朗日计算步中，流程与基于全局弱形式的无网格法相似。首先采用节点离散模型，引入覆盖全域的背景网格，在每个背景网格内生成相应高斯点用于数值积分。为方便计算，文本采用矩形背景网格，其位置计算与有限元高斯点积分点相同，如图 1所示。相比于有限元通过单元内节点计算高斯点形函数，无网格法则根据高斯点支持域内节点计算形函数。对边界内高斯点进行数值积分得到相应的刚度阵和内、外力向量(边界外高斯点对数值积分无贡献)，进而建立离散化的平衡方程，求解节点及高斯点相关信息。
本文采用径向基点插值函数(RPIM)构造形函数，插值点位移可表示为
$$ u(x) = \sum\limits_{i = 1}^n {{R_i}(x){a_i}} + \sum\limits_{j = 1}^m {{p_j}(x){b_j}} 。 $$ (1) 通过上式可以看出RPIM形函数由径向基Ri(x)和附加基pj(x)组成。n表示插值点支持域内节点个数，m表示添加附加基的个数，ai和bj是形函数待求常量，本文采用复合二次径向基(MQ)，可表示为
$$ {R}_{i}(r)={({r}_{i}{}_{}^{2}+{(c{d}_{av})}^{2})}^{q} \text{，} $$ (2) $$ {r}_{i}=\sqrt{{(x-{x}_{i})}^{2}+{(y-{y}_{i})}^{2}} 。 $$ (3) Liu等[19]通过多组数值试验得出结论：当ac=3，q=0.98或1.03时，在模拟固体及流体问题中有较好的效果，本文研究中ac=3，q=1.03。对于附加基pj(x)，现有研究成果表明：在RPIM形函数的基础上引入多项式附加线性基可以提高形函数模拟精度、降低对形状参数的敏感性及保证形函数的再生性。本文采用完备的线性基作为附加基。
$$ {p}_{j}(x)=[\begin{array}{ccc}1& x& y\end{array}] \text{，} $$ (4) 式(1)可表示为矩阵形式：
$$U = {R^T}a{\bf{ + }}{p^T}b$$ (5) 其中径向基矩阵：
$$ {R}^{T}={\left[\begin{array}{cccc}{R}_{1}({r}_{1})& {R}_{2}({r}_{1})& \mathrm{...}& {R}_{n}({r}_{1})\\ {R}_{1}({r}_{2})& {R}_{2}({r}_{2})& \mathrm{...}& {R}_{n}({r}_{2})\\ \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}\\ {R}_{1}({r}_{n})& {R}_{2}({r}_{n})& \mathrm{...}& {R}_{n}({r}_{n})\end{array}\right]}_{n\times n} \text{，} $$ (6) 附加基矩阵：
$$ {P}^{T}={\left[\begin{array}{cccc}1& 1& \cdots & 1\\ {x}_{1}& {x}_{2}& \cdots & {x}_{n}\\ {y}_{1}& {y}_{2}& \cdots & {y}_{n}\\ ⋮& ⋮& \ddots & ⋮\\ {p}_{m}({x}_{1})& {p}_{m}({x}_{2})& \cdots & {p}_{m}({x}_{n})\end{array}\right]}_{m\times n} 。 $$ (7) 方程组共有m+n个变量，而仅有n个方程，因此需要引入m个约束条件，添加m个方程：
$$ {\displaystyle \sum _{i=1}^{n}{p}_{j}({x}_{i}){a}_{i}}={P}^{T}a=0,j=1,2,\cdots ,m \text{，} $$ (8) 联立式(1)，(8)，可得到以下矩阵方程
$$ {U}_{t}=\left[\begin{array}{c}U\\ 0\end{array}\right]=\left[\begin{array}{cc}R& P\\ {P}^{T}& 0\end{array}\right]\left\{\begin{array}{c}a\\ b\end{array}\right\}=G{a}_{0} \text{，} $$ (9) 将待求常数代入，式(1)可改写为
$$ u(x) = \left\{ {{R^T}(x){P^T}(x)} \right\}{{\bf{G}}^{ - 1}}{U_t} = {\phi ^T}(x){U_t} 。 $$ (10) 最终插值点位移及位移偏导可表示为
$$ u(x)={\varphi }^{T}(x){U}_{t} \text{，} $$ (11) $$ {u}_{,l}(x)={\varphi }_{,l}{}^{T}(x){U}_{t} 。 $$ (12) 式中：l为x方向坐标或y方向坐标，逗号为对相应坐标求偏导。进一步该点应变、应力可表示为
$$ \varepsilon = \left[ {\begin{array}{*{20}{c}} {\frac{{\partial {\phi _1}}}{{\partial x}}}&0& \cdots &{\frac{{\partial {\phi _n}}}{{\partial x}}}&0 \\ 0&{\frac{{\partial {\phi _1}}}{{\partial y}}}& \cdots &0&{\frac{{\partial {\phi _n}}}{{\partial y}}} \\ {\frac{{\partial {\phi _1}}}{{\partial y}}}&{\frac{{\partial {\phi _1}}}{{\partial x}}}& \cdots &{\frac{{\partial {\phi _n}}}{{\partial y}}}&{\frac{{\partial {\phi _n}}}{{\partial x}}} \end{array}} \right]\left\{ {\begin{array}{*{20}{c}} {{u_1}} \\ {{v_1}} \\ \vdots \\ {{u_n}} \\ {{v_n}} \end{array}} \right\} = Bu 。 $$ (13) $$ \sigma =D\epsilon =DBu 。 $$ (14) 引入本构矩阵，该高斯点处刚度阵可表示为
$$ {K}_{e}=BDB 。 $$ (15) 同样的，可根据节点自由度组装整体体力向量和外力向量，问题域平衡方程离散形式可表示为
$$ Ku={F}^{b}+{F}^{t} 。 $$ (16) 随后在每个欧拉计算步中，根据节点位移增量，更新节点坐标，在此基础上根据式(11)中插值函数计算高斯点位移增量并更新其位置。图 2表示欧拉步中节点、高斯点在固定的背景网格上流动的过程。
![]()
图 2 节点、高斯点在背景网格上流动Figure 2. Update location of Gauss points and nodes on background mesh根据节点位移计算高斯点变形：
$$ {\overline{(\begin{array}{cc}x& y\end{array})}}_{node}={(\begin{array}{cc}x& y\end{array})}_{node}+\Delta {(\begin{array}{cc}x& y\end{array})}_{node} \text{，} $$ (17) $$ {\overline{(\begin{array}{cc}x& y\end{array})}}_{gauss}={(\begin{array}{cc}x& y\end{array})}_{gauss}+\varphi \Delta {(\begin{array}{cc}x& y\end{array})}_{node} 。 $$ (18) 式中下标表示更新后的节点或高斯点坐标。其中高斯点携带相关信息(应力、应变、本构内变量)变形后，由于位置改变，不再适用于下一个计算步的数值积分，需在固定的背景网格上引入一套新的高斯点，并将变形后高斯点相关信息映射回新的高斯点，如图 3所示。该过程同样采用RPIM插值函数，计算旧高斯点到新高斯点的映射关系。
$$u(x) = \sum\limits_{i = 1}^n {{R_i}(x){a_i}} + \sum\limits_{j = 1}^m {{p_j}(x){b_j}} $$ (19) $$ \kappa ' = {\phi _{gauss}}\kappa 。 $$ (20) 式中：$ {\phi _{gauss}} $为高斯点间插值函数；下标为新高斯点；κ为本构内变量，其维度由不同的本构模型决定。
2. 无网格大变形分析方法算例验证
地基承载力研究是典型的大变形算例，基于多种大变形分析方法(包括RITSS、CEL等)已开展对该问题的研究。本文采用无网格大变形分析方法模拟地基承载力问题，并通过与已有结果对比验证该方法的精度。图 4(a)为地基几何模型(考虑对称性仅给出一半)，地基土区域尺寸为4B×4B，其中B=2 m，为顶部基础宽度。土体采用基于Mohr-Coulomb准则的理性弹塑性模型模拟，弹性模量为100 kPa，泊松比为0.495，摩擦角为0°，黏聚力为1 kPa，并忽略自重。基础部分则视为刚体，在顶部分2000步施加0.8 m(0.4B)的位移边界条件，模拟其贯入过程。图 4(b)为网格节点模型，在基础附近1B×1B的区域节点平均间距控制为0.05B，并逐步过渡到0.2B。
图 5给出归一化地基反力f=q/cu与归一化贯入深度w=z/B的关系曲线。为更好验证方法的精度，RITSS[20]、CEL[21]计算结果及基于小变形假设的地基承载力理论解均汇总于图 5。结果表明：基于无网格大变形分析的计算结果与其他大变形分析方法结果曲线趋势一致，在归一化贯入深度达0.1前，归一化地基反力迅速增长，达0.1后反力增长缓慢，到0.3时反力基本稳定，最终求解的归一化承载力为6.24，验证了本文大变形分析方法具备足够的精度。
3. Meshless-FEM-SBFEM耦合分析
基于Meshless-FEM-SBFEM耦合大变形分析方法，可以协同发挥Meshless在大变形分析、SBFEM在跨尺度及FEM在计算效率中的优势。由于3种方法的形函数在节点处具备Kronecker δ性质，在各区域边界处位移连续，可通过共用节点的形式，实现在刚度阵层次的无缝耦合，如图 6所示。各区域首先计算各自刚度矩阵，然后根据节点自由度，组装到耦合刚度矩阵中，同理可得到耦合外力向量。
Meshless-FEM-SBFEM耦合大变形分析方法同样是在任意拉格朗日-欧拉框架下进行的，其每个计算步可分为拉格朗日步和欧拉步，在拉格朗日步中各方法建立平衡方程及求解过程中，与小变形分析一致均采用Cauchy应力作为应力度量。在随后的欧拉步中更新节点坐标、插值应力及应变信息，而由插值带来的内外力偏差将作为不平衡力带入下一计算步。
目前国际上还没有将SBFEM、FEM、Meshless耦合到一起的商用软件。本文采用C++语言及面向对象设计方法，将本文开发的无网格大变形分析方法集成到笔者自主开发的大型岩土工程计算系统GEODYNA。在笔者已有研究的基础上，在统一框架下实现了Meshless、FEM和SBFEM耦合分析，并联合土体广义弹塑性模型，为深厚覆盖层上土质心墙坝防渗墙-心墙接头局部大变形模拟提供有力的技术支撑。图 7以无网格大变形分析法与FEM耦合为例，给出数值实现的流程图。
![]()
图 7 Meshless-FEM耦合数值实现流程Figure 7. Flow chart of numerical application of Meshless-FEM4. 深厚覆盖层上土质心墙坝防渗墙-心墙接头局部无网格大变形模拟
将开发的无网格大变形分析方法应用于某深厚覆盖层上土质心墙坝中，再现其防渗墙“贯入”心墙的局部大变形发展过程。文献研究表明[2, 5]：在土质心墙坝填筑阶段，由于刚度差异，混凝土防渗墙与周围土体不均匀沉降现象最为明显。进而防渗墙“贯入”心墙现象也最为强烈，因此本文仅对坝体填筑阶段进行模拟。首先开展基于FEM-SBFEM小变形分析，确定潜在大变形区域，在此基础上，采用局部无网格节点替换潜在大变形区的有限元网格，开展大变形分析。通过对比模拟结果，揭示小变形分析的局限性。最后讨论无网格节点密度对模拟结果的影响。
4.1 算例模型
如图 8所示，坝高240 m，坝顶宽度16 m，覆盖层深度58 m。防渗墙长65.5 m，底部宽度1.2 m，顶部宽度0.4 m。坝体分55步进行填筑，为保证计算效率及精度，坝体平均网格尺寸控制约为5 m。同时为实现对心墙-防渗墙局部精细化分析，采用SBFEM实现从整体到局部的跨尺度建模，控制该区域结构及周围土体平均网格尺寸约为0.25 m，最终模型共计节点数18255，共计单元数17368。
![]()
图 8 坝体网格模型及心墙-防渗墙接头局部网格Figure 8. Mesh of dam and local mesh of cut-off wall-core wall zone4.2 材料参数
为更精确地捕捉填筑过程中坝体材料的力学行为，本节中堆石料(A区)、心墙料(B区)、过渡料(C区)及覆盖层(D区)区域均采用堆石体广义塑性模型[22-24]。具体参数汇总于表 1。表中G0，K0，ms，mv为与弹性模量相关系数，Mg，Mf，af，ag为与塑性流动方向相关系数，其余为与塑性模量相关系数。防渗墙及基岩部分采用线弹性模型，详细参数见表 2。这里需要指出在土体与防渗墙之间设置薄层实体单元材料与周围土体本构一致。
表 1 堆石料、心墙料、过渡料及覆盖层区域广义塑性模型参数Table 1. Parameters of generalized plastic model in rock-fill zone, core-wall, transition zone and overburden zone分区 G0 K0 Mg Mf af ag H0 HU0 A 1000 1400 1.80 1.38 0.45 0.40 1800 3000 B 720 662 1.45 1.20 0.20 0.20 400 400 C 1321 1585 1.64 1.15 0.20 0.45 1800 1300 D 382 351 1.35 1.00 0.25 0.60 500 500 ms mv m1 mu rd 𝛾DM 𝛾u 𝛽0 𝛽1 0.50 0.50 0.20 0.20 105 50 4 35 0.022 0.47 0.47 0.40 0.40 7.5 20 3.5 20 0.013 0.50 0.50 0.41 0.41 120 20 7 22 0.018 0.612 0.612 0.44 0.44 10 12 3.5 20 0.01 表 2 线弹性模型参数Table 2. Parameters of linear elastic model线弹性材料 E/MPa 𝜇 基岩 13000 0.25 防渗墙 30000 0.167 4.3 小变形结果分析——确定潜在大变形区
整个填筑过采用非线性迭代法，最大迭代次数为50，收敛条件为不平衡力小于5%。最终给出竣工期心墙-防渗墙接头区域网格变形图及相应的剪应变分布情况，如图 9所示。
![]()
图 9 竣工期局部网格变形图及剪应变分布Figure 9. Deformation mesh and shear strains after dam construction随着填筑步的增加，心墙-防渗墙接头局部区域土体网格出现翻转畸形，同时防渗墙顶部两侧土体的局部剪应变逐步累积，在竣工期已超400%。在如此强烈的剪切变形下，小变形分析难以满足精度要求。
4.4 无网格大变形结果分析
基于小变形分析结果，对剪应变超过30%的矩形区域(20 m×20 m)，用无网格节点代替有限元网格开展无网格大变形分析，其余位置保持网格不变，在每个计算步完成后更新节点位置，如图 10所示。
在保证与小变形分析材料参数、分析步数及收敛精度一致的条件下，开展基于Meshless-FEM-SBFEM模型的大变形分析。首先给出两者在竣工期坝体沉降及水平位移的对比，如图 11所示。通过对比发现两者在坝体变形模拟中计算结果一致。
![]()
图 11 竣工期坝体变形对比(小变形分析vs大变形分析)Figure 11. Comparison of dam deformations after dam construction (small deformation analysis vs. large deformation analysis)图 12给出大变形分析在典型时刻心墙-防渗墙接头土体区域节点变形图。结果表明：采用无网格大变形分析，在本质上避免了网格畸形问题。图中蓝色节点为初始位置在同一水平线的节点对(仅左侧，右侧同理)，图中随着填筑步增加，两者发生明显的剪切滑移，最终在防渗墙顶部土体间产生剪切带。通过与图 9(b)的对比，可以发现剪切带位置与小变形分析中高剪应变区基本吻合。
![]()
图 12 大变形分析典型时刻心墙-防渗墙接头区域节点变形Figure 12. Nodal deformations of cut-off wall-core wall zone at typical time by large deformation analysis图 13给出竣工期心墙-防渗墙接头附近土体水平位移对比(为方便对比，均采用变形前区域)。考虑大变形分析后，防渗墙顶部土体水平位移明显增加，反映了伴随混凝土防渗墙“贯入”土质心墙过程，顶部土体发生向两侧排土的现象。
小变形与大变形分析局部土体应力场也存在明显的差异，图 14给出竣工期心墙-防渗墙接头附近土体应力水平对比。相比于小变形分析，大变形分析中防渗墙顶部局部土体出现较高应力水平区(大于0.9)，且基本与图 12中剪切带位置吻合，反映了防渗墙顶部土体发挥较大的强度阻碍防渗墙“贯入”心墙。
防渗墙应力是整个防渗系统的重中之重，由于顶部土体变形模式及应力水平不同，基于小变形分析及大变形分析的防渗墙应力也存在相应的差异，图 15给出竣工期防渗墙竖向应力对比。
![]()
图 15 竣工期防渗墙竖向应力对比Figure 15. Comparison of vertical stresses of cut-off wall after dam construction两种方法的应力偏差主要来自防渗墙顶部，考虑大变形分析后，顶部土体向两侧位移明显，土体发挥较大的强度阻碍防渗墙贯入，增大了防渗墙顶部与周围土体的摩擦作用，进而增加了防渗墙竖向应力。考虑大变形分析后，防渗墙在竣工期竖向应力增大约4.1 MPa，相对偏差约为13%。
最后开展心墙-防渗墙接头附近土体节点密度对防渗墙应力影响的研究。在保证其它区域网格不变的前提下，控制心墙-防渗墙接头附近区域节点为0.125，0.5，1 m，开展大变形分析。表 3汇总了各模型防渗墙顶部及最大竖向应力。结果表明：节点平均密度对防渗墙应力模拟有一定的影响，且随着节点间距减小，模拟结果趋于收敛，0.125 m和0.25 m的节点平均间距模型计算结果接近(相对偏差 < 2%)，因此综合考虑计算效率及精度，本文建议节点平均间距为0.25 m(约为防渗墙厚度的0.2倍)。
表 3 节点平均间距对防渗墙应力影响Table 3. Influences of average nodal displacement on stress of cut-off wall计算模型 顶部应力/MPa 最大应力/MPa 平均间距0.125 m -13.6 -35.9 平均间距0.25 m -13.8 -36.1 平均间距0.5m -14.7 -37.5 平均间距1 m -16.1 -39.2 5. 讨论
随着岩土工程对大变形模拟需求的不断提高，近几十年来众学者开发了RITSS、CEL、光滑粒子流(SPH)、物质点法(MPM)、粒子有限元(PFEM)以及一系列在此基础上发展并衍生的大变形分析方法。然而相比于大变形方法高速发展，其在大型岩土工程中应用的研究相对较少(尤其在三维问题模拟中)。其原因在于商用程序调试困难且专业性较差，难以满足岩土工程计算需求；而自主研发软件工作量大且计算能力有限，无法开展大规模计算。因此本文建议在现有自主研发的有限元程序FEM框架下开发大变形分析方法或实现与FEM的耦合分析，既能够集成FEM中丰富的材料本构及成熟的静、动力分析方法，又可以通过局部大变形分析的方法提高计算效率，进而满足大型岩土工程的计算需求。同时地震引发土体液化现象也是岩土工程中常见的大变形问题，因此考虑孔隙水压力的累积-消散过程与土骨架变形、应力耦合的大变形分析有着较高理论价值及工程意义。
6. 结论
本文在ALE框架下，开发了弹塑性无网格大变形分析方法，摆脱了FEM在大变形过程中因网格畸形带来的精度降低或复杂的重剖分问题。将该方法集成到自主开发的GEODYNA计算系统中，联合土体广义弹塑性模型，再现了深厚覆盖层上土质心墙坝心墙-防渗墙接头局部大变形现象。
(1)本文建立了基于Meshless-FEM-SBFEM的多数值耦合分析方法，协同发挥了Meshless在大变形、SBFEM在跨尺度、FEM在计算效率方面的优势。
(2)基于小变形分析与本文大变形分析在大坝整体变形模拟中计算结果一致，而在防渗墙-心墙接头局部大变形区域，两者计算结果出现一定差异。
(3)在防渗墙贯入高塑性土过程中，墙顶土体间出现了明显的剪切带，位置与小变形高剪应变区基本吻合。但相比小变形分析结果，大变形顶部土体侧向位移更大。
(4)考虑大变形分析后，墙顶土体出现较高的应力水平区域，反映了土体发挥了较大的强度，阻碍防渗墙贯入。因此加剧了防渗墙顶部与周围土体的摩擦作用，增大了防渗墙竖向应力。对比结果表明，小变形分析低估了4.1 MPa防渗墙竖向应力(约为13%)。
(5)心墙-防渗墙接头附近土体节点密度对防渗墙应力模拟存在一定的影响，且随着节点间距的减小，模拟结果趋于收敛。综合考虑计算效率及精度，本文建议节点平均间距为0.25 m(约为防渗墙厚度的0.2倍)。
-
![]()
图 2 节点、高斯点在背景网格上流动
Figure 2. Update location of Gauss points and nodes on background mesh
![]()
图 7 Meshless-FEM耦合数值实现流程
Figure 7. Flow chart of numerical application of Meshless-FEM
![]()
图 8 坝体网格模型及心墙-防渗墙接头局部网格
Figure 8. Mesh of dam and local mesh of cut-off wall-core wall zone
![]()
图 9 竣工期局部网格变形图及剪应变分布
Figure 9. Deformation mesh and shear strains after dam construction
![]()
图 11 竣工期坝体变形对比(小变形分析vs大变形分析)
Figure 11. Comparison of dam deformations after dam construction (small deformation analysis vs. large deformation analysis)
![]()
图 12 大变形分析典型时刻心墙-防渗墙接头区域节点变形
Figure 12. Nodal deformations of cut-off wall-core wall zone at typical time by large deformation analysis
![]()
图 15 竣工期防渗墙竖向应力对比
Figure 15. Comparison of vertical stresses of cut-off wall after dam construction
表 1 堆石料、心墙料、过渡料及覆盖层区域广义塑性模型参数
Table 1 Parameters of generalized plastic model in rock-fill zone, core-wall, transition zone and overburden zone
分区 G0 K0 Mg Mf af ag H0 HU0 A 1000 1400 1.80 1.38 0.45 0.40 1800 3000 B 720 662 1.45 1.20 0.20 0.20 400 400 C 1321 1585 1.64 1.15 0.20 0.45 1800 1300 D 382 351 1.35 1.00 0.25 0.60 500 500 ms mv m1 mu rd 𝛾DM 𝛾u 𝛽0 𝛽1 0.50 0.50 0.20 0.20 105 50 4 35 0.022 0.47 0.47 0.40 0.40 7.5 20 3.5 20 0.013 0.50 0.50 0.41 0.41 120 20 7 22 0.018 0.612 0.612 0.44 0.44 10 12 3.5 20 0.01 
下载: 导出CSV
表 2 线弹性模型参数
Table 2 Parameters of linear elastic model
线弹性材料 E/MPa 𝜇 基岩 13000 0.25 防渗墙 30000 0.167
下载: 导出CSV
表 3 节点平均间距对防渗墙应力影响
Table 3 Influences of average nodal displacement on stress of cut-off wall
计算模型 顶部应力/MPa 最大应力/MPa 平均间距0.125 m -13.6 -35.9 平均间距0.25 m -13.8 -36.1 平均间距0.5m -14.7 -37.5 平均间距1 m -16.1 -39.2
下载: 导出CSV
-
[1] 金峰, 周虎, 李玲玉, 等. 堆石混凝土系列技术在西部水电工程中的应用[J]. 水电与抽水蓄能, 2021, 7(1): 16-22. https://www.cnki.com.cn/Article/CJFDTOTAL-DBGC202101005.htm JIN Feng, ZHOU Hu, LI Lingyu, et al. Application of innovative technologies based on rock-filled concrete in hydropower projects in western China[J]. Hydropower and Pumped Storage, 2021, 7(1): 16-22. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-DBGC202101005.htm
[2] YU X, ZOU D G, KONG X J, YU L. Large-deformation finite element analysis of the interaction between concrete cut-off walls and high-plasticity clay in an earth core dam[J]. Engineering Computations, 2017, 34(4): 1126-1148. doi: 10.1108/EC-04-2016-0118
[3] 丁艳辉, 张其光, 张丙印. 高心墙堆石坝防渗墙应力变形特性有限元分析[J]. 水力发电学报, 2013, 32(3): 162-167. https://www.cnki.com.cn/Article/CJFDTOTAL-SFXB201303029.htm DING Yanhui, ZHANG Qiguang, ZHANG Bingyin. FEM analysis of stress-deformation characteristics of cut-off walls in high core rockfill dam[J]. Journal of Hydroelectric Engineering, 2013, 32(3): 162-167. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-SFXB201303029.htm
[4] 沈振中, 田振宇, 徐力群, 等. 深覆盖层上土石坝心墙与防渗墙连接型式研究[J]. 岩土工程学报, 2017, 39(5): 939-945. doi: 10.11779/CJGE201705019 SHEN Zhenzhong, TIAN Zhenyu, XU Liqun, et al. Reasonable connection type for cutoff wall and core wall of earth-rock dams on deep overburden layers[J]. Chinese Journal of Geotechnical Engineering, 2017, 39(5): 939-945. (in Chinese) doi: 10.11779/CJGE201705019
[5] 温立峰, 李炎隆, 柴军瑞. 坝基混凝土防渗墙力学性状的统计分析[J]. 水利学报, 2021, 52(2): 241-254. https://www.cnki.com.cn/Article/CJFDTOTAL-SLXB202102012.htm WEN Lifeng, LI Yanlong, CHAI Junrui. Statistical analysis of mechanical properties of dam foundation concrete cutoff wall[J]. Journal of Hydraulic Engineering, 2021, 52(2): 241-254. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-SLXB202102012.htm
[6] WANG D, BIENEN B, NAZEM M, et al. Large deformation finite element analyses in geotechnical engineering[J]. Computers & Geotechnics, 2015, 65(S): 104-114.
[7] AUGARDE C E, LEE S J, LOUKIDIS D. Numerical modelling of large deformation problems in geotechnical engineering: a state-of-the-art review[J]. Soils and Foundations, 2021, 61(6): 1718-1735.
[8] 王栋, 年廷凯, 陈煜淼. 边坡稳定有限元分析中的三个问题[J]. 岩土力学, 2007, 28(11): 2309-2313, 2318. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201007045.htm WANG Dong, NIAN Tingkai, CHEN Yumiao. Three problems in slope stability analyses with finite element method[J]. Rock and Soil Mechanics, 2007, 28(11): 2309-2313, 2318. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201007045.htm
[9] MILLER K, JOLDES G, LANCE D E, et al. Total Lagrangian explicit dynamics finite element algorithm for computing soft tissue deformation[J]. Communications in Numerical Methods in Engineering, 2007, 23(2): 121-134.
[10] BATHE K J, RAMM E, WILSON E L. Finite element formulations for large deformation dynamic analysis[J]. International Journal for Numerical Methods in Engineering, 1975, 9(2): 353-386.
[11] HU Y, RANDOLPH M F. A practical numerical approach for large deformation problems in soil[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 1998, 22(5): 327-350.
[12] BENSON D J. Computational methods in Lagrangian and eulerian hydrocodes[J]. Computer Methods in Applied Mechanics and Engineering, 1992, 99(2/3): 235-394.
[13] LUCY L B. A numerical approach to the testing of the fission hypothesis[J]. The Astronomical Journal, 1977, 82: 1013.
[14] 黄雨, 郝亮, 野々山人. SPH方法在岩土工程中的研究应用进展[J]. 岩土工程学报, 2008, 30(2): 256-262. http://www.cgejournal.com/cn/article/id/12763 HUANG Yu, HAO Liang, NONOYAMA Hideto. The state of the art of SPH method applied in geotechnical engineering[J]. Chinese Journal of Geotechnical Engineering, 2008, 30(2): 256-262. (in Chinese) http://www.cgejournal.com/cn/article/id/12763
[15] BELYTSCHKO T, LU Y Y, GU L. Element-free Galerkin methods[J]. International Journal for Numerical Methods in Engineering, 1994, 37(2): 229-256.
[16] 邹德高, 龚瑾, 孔宪京, 等. 基于无网格界面模拟方法的面板坝防渗体跨尺度分析[J]. 水利学报, 2019, 50(12): 1446-1453, 1466. https://www.cnki.com.cn/Article/CJFDTOTAL-SLXB201912003.htm ZOU Degao, GONG Jin, KONG Xianjing, et al. The cross-scale analysis of concrete-face rock-fill dam anti-seepage structure based on the simulation of meshfree interface[J]. Journal of Hydraulic Engineering, 2019, 50(12): 1446-1453, 1466. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-SLXB201912003.htm
[17] JIE Y, TANG X W, LUAN M T, et al. Adaptive element free Galerkin method applied to analysis of earthquake induced liquefaction[J]. Earthquake Engineering and Engineering Vibration, 2008, 7(2): 217-224.
[18] CHEN K, ZOU D G, KONG X J, et al. An efficient nonlinear octree SBFEM and its application to complicated geotechnical structures[J]. Computers and Geotechnics, 2018, 96: 226-245.
[19] LIU G R, GU Y T. A local point interpolation method for stress analysis of two-dimensional solids[J]. Structural Engineering and Mechanics, 2001, 11(2): 221-236
[20] TIAN Y H, CASSIDY M J, RANDOLPH M F, et al. A simple implementation of RITSS and its application in large deformation analysis[J]. Computers and Geotechnics, 2014, 56: 160-167.
[21] 赵燕兵, 刘海笑. 条形基础嵌入问题的大变形有限元对比分析[J]. 天津大学学报(自然科学与工程技术版), 2016, 49(10): 1055-1061. https://www.cnki.com.cn/Article/CJFDTOTAL-TJDX201610008.htm ZHAO Yanbing, LIU Haixiao. Comparison of large deformation finite element methods for penetration of strip foundations[J]. Journal of Tianjin University (Science and Technology), 2016, 49(10): 1055-1061. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-TJDX201610008.htm
[22] 于玉贞, 卞锋. 高土石坝地震动力响应特征弹塑性有限元分析[J]. 世界地震工程, 2010, 26(增刊1): 341-345. https://www.cnki.com.cn/Article/CJFDTOTAL-SJDC2010S1069.htm YU Yuzhen, BIAN Feng. Elasto-plastic FEM analysis of dynamic response of high earth-rockfill dams during earthquake[J]. World Earthquake Engineering, 2010, 26(S1): 341-345. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-SJDC2010S1069.htm
[23] 周墨臻, 张丙印, 张宗亮, 等. 超高面板堆石坝面板挤压破坏机理及数值模拟方法研究[J]. 岩土工程学报, 2015, 37(8): 1426-1432. doi: 10.11779/CJGE201508010 ZHOU Mozhen, ZHANG Bingyin, ZHANG Zongliang, et al. Mechanisms and simulation methods for extrusion damage of concrete faces of high concrete-faced rockfill dams[J]. Chinese Journal of Geotechnical Engineering, 2015, 37(8): 1426-1432. (in Chinese) doi: 10.11779/CJGE201508010
[24] ZOU D G, XU B, KONG X J, et al. Numerical simulation of the seismic response of the Zipingpu concrete face rockfill dam during the Wenchuan earthquake based on a generalized plasticity model[J]. Computers and Geotechnics, 2013, 49(S0): 111-122.
-
期刊类型引用(1)
1. 杨云洲. 烂泥山水库沥青混凝土心墙坝抗震安全分析. 陕西水利. 2024(08): 15-17 . 
百度学术
其他类型引用(0)
-
其他相关附件
计量
- 文章访问数: 341
- HTML全文浏览量: 37
- PDF下载量: 119
- 被引次数: 1
