基于三维比例边界有限元的高拱坝地震响应模拟

Seismic response simulation of high arch dams based on 3D scaled boundary finite element method

摘要

针对传统有限元法在高拱坝抗震分析中存在的网格依赖性强及横缝接触模拟困难等问题，引入比例边界有限元法（SBFEM），构建了高拱坝-地基系统的非线性地震响应分析模型。该模型结合八叉树网格划分技术与拉格朗日乘子法，实现了对大坝横缝接触效应的精确模拟。以某300 m级高拱坝为例，对比分析了SBFEM与ABAQUS软件的计算结果，验证了模型的有效性，并研究了横缝张开效应对大坝位移及应力分布的影响。结果表明：SBFEM能够有效模拟高拱坝的动力特性；横缝的存在显著改变了大坝的受力状态，使部分拱向应力向梁向转移。

Abstract

Aiming at the problems of strong mesh dependency and difficulty in simulating transverse joint contact in the seismic analysis of high arch dams using the traditional finite element method， the scaled boundary finite element method （SBFEM） was introduced， and a nonlinear seismic response analysis model for the high arch dam-foundation system was constructed. By combining octree meshing technology and the Lagrange multiplier method， accurate simulation of the contact effects at the dam's transverse joints was achieved by the model. Taking a 300-meter-class high arch dam as an example， the computational results of SBFEM and ABAQUS software were compared and analyzed， and the effectiveness of the model was verified. In addition， the influence of the transverse joint opening effect on the dam's displacement and stress distribution was studied. The results indicated that the dynamic characteristics of high arch dams could be effectively simulated by SBFEM； The stress state of the dam was significantly altered by the presence of transverse joints， causing part of the arch-direction load to be transferred to the cantilever direction.

Graphical abstract

关键词

比例边界有限元法（SBFEM） / 高拱坝 / 地震响应 / 横缝 / 八叉树网格 / 拉格朗日乘子

Key words

scaled boundary finite element method （SBFEM） / high arch dam / seismic response / transverse joint / Octree mesh / Lagrange multiplier

0 引言

我国西南地区水能资源丰富，已建成或在建多座300 m级特高拱坝，这些工程多处于高烈度地震区，其抗震安全性问题备受学术界与工程界的关注^［1-2］。为适应地基变形及防止混凝土收缩开裂，高拱坝通常沿拱向设置若干横缝。虽然横缝在静力作用下通常保持闭合，但在强震作用下可能发生反复地张开、闭合及滑移，这种非线性接触行为将显著改变结构的动力特性与应力分布，甚至诱发局部破坏。因此，在数值模拟中精确考虑横缝的动态接触效应^［3］，是获取高拱坝真实地震响应、科学评估其抗震安全性的前提。

尽管有限元法（FEM）在高拱坝抗震分析中应用广泛，但在处理包含复杂几何边界及横缝接触界面的大坝-地基系统时，其前处理网格剖分往往极为耗时且面临应力奇异性模拟的挑战^［4-5］。相比之下，SONG等^［6］提出的比例边界有限元法（SBFEM）为解决此类难题提供了全新的途径。作为一种半解析数值方法，SBFEM仅需在环向进行有限元离散，径向采用解析解，这使其能够精确捕捉裂纹尖端或角点处的应力奇异性，且支持构造任意面数的多面体单元。这一特性使得SBFEM与八叉树（Octree）等自动网格生成技术完美结合，实现复杂结构网格划分的高度自动化^［7-8］。

近年来，SBFEM在断裂力学^［9-10］、无限域波动^［11］和不确定性^［12］等问题中取得了长足进展，但将其应用于高拱坝这类大型复杂三维结构的非线性接触动力分析的相关研究尚不多见。鉴于此，本研究基于三维SBFEM理论框架，结合高效的八叉树网格生成技术，构建了高拱坝-地基系统的整体分析模型。研究中引入拉格朗日（Lagrange）乘子法来严格施加横缝界面的非线性接触约束，实现了对横缝张开、闭合和滑移行为的精确模拟。以某300 m级高拱坝为例，通过与商业有限元软件ABAQUS计算结果的对比，验证了该方法的有效性，并探讨了横缝局部非线性对大坝整体地震响应和损伤演化的影响机制，为高拱坝抗震安全评价提供了一条新的数值模拟途径。

1 三维比例边界有限元理论与方法

1.1 三维比例边界有限元基本方程

在比例边界坐标系（\xi ,\eta ,\zeta）中，\xi为径向坐标，\eta、\zeta为环向坐标。在比例边界有限元法中需要设置比例中心，比例中心对子域内任意一点可见，在边界上进行离散化，在多面体上的每个面都要离散成三角形和四边形，如图1所示。子域内任意一点b的笛卡儿坐标（ x， y， z ），可由\xi、\eta和\zeta表示：

\mathit{x}(\xi ,\eta ,\zeta )=\xi \mathit{N}(\eta ,\zeta ){\mathit{x}}_{\mathrm{b}}，(1)

\mathit{y}(\xi ,\eta ,\zeta )=\xi \mathit{N}(\eta ,\zeta ){\mathit{y}}_{\mathrm{b}}，(2)

\mathit{z}(\xi ,\eta ,\zeta )=\xi \mathit{N}(\eta ,\zeta ){\mathit{z}}_{\mathrm{b}}，(3)

式中：\mathit{N}(\eta ,\zeta )为二维形函数；{\mathit{x}}_{\mathrm{b}}、{\mathit{y}}_{\mathrm{b}}和{\mathit{z}}_{\mathrm{b}}为边界上结点坐标集。

在\xi =\mathrm{1}的单元边界上，采用类似有限元的位移插值方法，线单元的位移场可通过形函数\mathit{N}(\eta ,\zeta )对节点位移进行插值求解。假设单元对应径向常数\xi的环向位移都采用相同的形函数\mathit{N}(\eta ,\zeta )插值，设位移解形式：

\mathit{u}(\xi ,\eta ,\zeta )=\mathit{N}(\eta ,\zeta )\mathit{u}(\xi )，(4)

式中：\mathit{u}(\xi )为径向上各点位移。

三维比例边界有限元关于位移的控制方程可采用加权余量法或虚功原理^［8，13］推导：

{\mathit{E}}_{\mathrm{0}}{\xi }^{\mathrm{2}}\mathit{u}(\xi ){,}_{\xi \xi }+\left(\mathrm{2}{\mathit{E}}_{\mathrm{0}}+{\mathit{E}}_{\mathrm{1}}^{\mathrm{T}}-{\mathit{E}}_{\mathrm{1}}\right)\xi \mathit{u}(\xi ){,}_{\xi }+({\mathit{E}}_{\mathrm{1}}^{\mathrm{T}}-{\mathit{E}}_{\mathrm{2}})\mathit{u}(\xi )=\mathrm{0}，(5)

式中：{\mathit{E}}_{\mathrm{0}}、{\mathit{E}}_{\mathrm{1}}和{\mathit{E}}_{\mathrm{2}}为取决于边界的几何系数矩阵。

求解该控制方程，可推导得到环向位移解析解，并代入式（4）中即可得到单元内任一点的位移和相应的应力，具体计算详见文献［8］和［13］。

1.2 基于SBFEM八叉树网格-FEM六面体网格耦合

高拱坝是典型的薄壳结构，坝体厚度相对于其跨度和高度较小，且受力状态复杂，为精确捕捉坝体内部的应力梯度及横缝的开合状态，必须采用较密的网格进行离散^［14］。大坝地基模拟范围通常延伸至数倍坝高，若全域采用均匀细网格，计算自由度将呈指数级增长，严重制约计算效率。为了解决坝体“薄结构、密网格”与地基“大体积、疏网格”之间的矛盾，兼顾计算精度与效率，本研究采用SBFEM八叉树-FEM六面体耦合网格建模策略。

针对坝体区域几何尺度小且需精细分析的特点，采用成熟的有限元（FEM）六面体网格进行精细离散。在地基区域，采用SBFEM多面体单元，利用八叉树算法对地基进行递归剖分，生成疏密过渡的六面体网格，并允许存在悬挂节点（Hanging Nodes）^［15］。SBFEM将每个八叉树方格视为一个多面体子域，能够直接处理悬挂节点带来的不连续性，无需进行复杂的过渡网格加密。在坝体与地基交界面处，通过节点自由度映射实现位移协调，保证能量传递的连续性。

2 比例边界有限元横缝动接触模型

横缝的动力接触问题本质上是一个含不等式约束的边界非线性问题。本研究采用Lagrange乘子法来精确满足接触面的非嵌入条件^［16］。该方法是求解接触约束最小化问题的主要方法之一，通过引入Lagrange乘子将接触问题转化为无约束问题求解。接触状态有3种：张开、黏结和滑移。用Lagrange乘子{\lambda }_{\mathrm{N}}和{\lambda }_{\mathrm{T}}分别表示发生接触时，接触面上的法向和切向接触力。法向Lagrange乘子场{\lambda }_{\mathrm{N}}表示当接触发生时迫使u_{\mathrm{N}}（法向间隙）等于0（裂纹面接触不嵌入）的压力，切向Lagrange乘子场{\lambda }_{\mathrm{T}}表示在黏结区迫使满足黏结条件的切向力。

在满足库伦摩擦定律条件下，分离、黏结和滑移3种接触状态^［17］，可以通过应用Lagrange乘子法添加摩擦接触条件得到以下3个方程。

1） 接触不嵌入条件

{\int }_{S_{\mathrm{c}}}{u}_{\mathrm{N}}\cdot \delta {\lambda }_{\mathrm{N}}\mathrm{d}S=\mathrm{0},\forall \delta {\lambda }_{\mathrm{N}}，(6)

式中：u_{\mathrm{N}}为接触点对法向相对距离，S_{\mathrm{c}}为接触界面。

2） 切向黏结条件

{\int }_{S_{\mathrm{T}}}\delta u_{\mathrm{T}}{\lambda }_{\mathrm{T}}\mathrm{d}S=\mathrm{0},|\begin{array}{cc}\forall \mathrm{\delta }{\mathrm{\lambda }}_{\mathrm{N}}& \left({\mathrm{\lambda }}_{\mathrm{T}}\le \mathrm{\mu }‖{\mathrm{\lambda }}_{\mathrm{N}}‖\right)\end{array}，(7)

式中：u_{\mathrm{T}}为接触点对切向相对距离。

在切向满足滑动接触状态时，切向运动不再受约束，而切向力满足：

{\lambda }_{\mathrm{T}}=-\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}(u_{\mathrm{T}})\mu \left|{\lambda }_{\mathrm{N}}\right|，(8)

式中：\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}符号函数，\mu为接触面摩擦系数。

3） 内外力（接触力作为一种外力）应满足系统的虚功原理，即：

\begin{array}{l}\delta W^{\mathrm{i}\mathrm{n}\mathrm{t}}(\delta u,u)-\delta W^{\mathrm{e}\mathrm{x}\mathrm{t}}(\delta u,u)-\\ {\int }_{S_{\mathrm{c}}}\delta u_{\mathrm{N}}{\lambda }_{\mathrm{N}}\mathrm{d}S-{\int }_{S_{\mathrm{T}}}\delta u_{\mathrm{T}}{\lambda }_{\mathrm{T}}\mathrm{d}S-\\ {\int }_{S_{\mathrm{T}}}\delta u_{\mathrm{T}}\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}(u_{\mathrm{T}})\mu {\lambda }_{\mathrm{N}}\mathrm{d}S=\mathrm{0}\end{array}。（9）

比例边界有限元中，横缝面接触分别在两个单元的边界上（\xi =\mathrm{1}），采用插值得到的Lagrange乘子，根据虚功方程式（9）以及约束方程式（8）可得到裂纹面满足滑移接触的动接触控制方程：

\left[\begin{array}{cc}\mathit{M}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}\end{array}\right]\left[\begin{array}{c}{\ddot{\mathit{u}}}_{\mathrm{0}}\\ \mathrm{0}\end{array}\right]+\left[\begin{array}{cc}\mathit{C}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}\end{array}\right]\left[\begin{array}{c}{\dot{\mathit{u}}}_{\mathrm{0}}\\ \mathrm{0}\end{array}\right]+\left[\begin{array}{l}\mathit{K}{\mathit{C}}_{\mathrm{N}}^{\mathrm{T}}-\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}({\mathit{u}}_{\mathrm{T}})\mu \cdot {\mathit{C}}_{\mathrm{T}}^{\mathrm{T}}\\ {\mathit{C}}_{\mathrm{N}}\mathrm{0}\end{array}\right]\left[\begin{array}{c}{\mathit{u}}_{\mathrm{0}}\\ {\mathit{\lambda }}_{\mathrm{N}}\end{array}\right]=\left[\begin{array}{c}\mathit{F}\\ \mathrm{0}\end{array}\right]，(10)

式中：\mathit{M}为质量矩阵，{\ddot{\mathit{u}}}_{\mathrm{0}}为节点加速度，{\dot{\mathit{u}}}_{\mathrm{0}}为节点速度，{\mathit{C}}_{\mathrm{N}}^{\mathrm{T}} 和 {\mathit{C}}_{\mathrm{T}}^{\mathrm{T}}分别为法向和切向接触积分系数矩阵。

采用Newmark隐式算法对控制方程进行时间积分。

3 基于比例边界有限元的高拱坝-地基系统模型

3.1 工程概况与几何模型

以澜沧江流域某300 m级双曲拱坝为研究对象。大坝最大坝高294.5 m，坝顶高程1245 m。坝顶中心线弧长892.8 m，厚高比约为0.25，属于典型的特高薄拱坝。工程抗震设防烈度为9度。模型坐标系定义为：原点位于坝底拱冠处，X轴指向右岸，Y轴顺河向下游，Z轴竖直向上。为了消除边界截断对地震波反射的影响，地基模拟范围在各个方向延伸至2倍坝高，总尺寸约为1500 m×1500 m×1000 m。

根据该工程的勘察设计资料，表1为该坝的材料参数，动力分析时坝体混凝土弹性模量取为静弹模的1.5倍，阻尼采用瑞利阻尼，阻尼比按照规范取5%。坝基简化为均质各向同性的弹性材料，采用无质量地基来减小因地基截断边界产生的反射波对坝体的作用，以此来近似模拟无限地基作用。

3.2 网格划分

近场为有限元网格和远场为比例边界有限元八叉树网格，整体模型单元总数为24 328，节点总数为30 216（其中坝体单元数为5248，坝体节点数为7040）。根据研究，取7条或11条横缝模拟时，横缝最大缝宽值基本等同于25条和40条横缝模拟的结果，综合考虑后本次计算选择7条横缝进行模拟，从坝体上游面从左至右分别记为J-1、J-2、J-3、J-4、J-5、J-6和J-7。图2给出了高拱坝-地基SBFEM整体网格图、坝体网格图和横缝具体位置示意图。

3.3 计算结果与分析

为验证SBFEM接触模型的准确性，首先在空库及正常蓄水位工况下进行线弹性地震响应计算，并将结果与商业软件ABAQUS的计算结果进行比对。计算中采用El-Centro地震波，水平地震加速度峰值为0.344 g，图3给出了顺河向、横河向及竖向的加速度时程曲线。

图4和图5分别为空库时拱冠处上游顺河向位移时程曲线及J1缝处张开-闭合曲线，结果显示，SBFEM计算所得的拱冠处最大顺河向位移为0.102 m，ABAQUS计算结果为0.098 m，相对误差小于5%。在横缝张开度方面，左岸J1缝的最大张开度本文解为2.8 cm，ABAQUS解为2.75 cm。两者的位移时程曲线及应力云图分布吻合度极高，证明了本文提出的SBFEM-FEM耦合模型在处理接触非线性及波传播方面的精度满足工程要求。

在此基础上，计算正常高水位（水深为287 m）与地震共同作用的线弹性地震响应，图6为拱冠处的顺河向位移时程响应曲线比较，从图6中可看出，静动力响应的变化规律基本一致。图7（a）和图7（b）为两种方法的上游面最大主应力云图比较，可以看出，应力云图的分布规律基本相同。

3.4 横缝对高拱坝非线性地震响应的影响

为了研究横缝对高拱坝非线性地震响应的影响，基于SBFEM分别对坝体带有7条横缝以及无横缝两种情况进行了计算模拟，因为低水位时对横缝张开不利，本次计算荷载组合考虑了低水位（水深228 m）与地震的共同作用。

图8给出了坝体上游面拱冠处分别在有无横缝时的位移响应，从图8中可看出，无横缝情况的初始位移为0.334 m，小于有横缝情况下的初始位移0.524 m，这是由于横缝张开使得径向位移变大，张开时结构处于分割的块体，抵抗变形的能力小于无横缝下的完整结构。考虑横缝情况下，位移在3.78 s处取得最大值，为0.729 m，在5.54 s处取得最小值，为0.322 m；而不考虑横缝情况下，位移在2.48 s处取得最大值，为0.689 m，在5.40 s处取得最小值，为-0.011 m。

图9给出了有无横缝模型两种情况下的顺河向位移取得最大值时的位移云图。在缝两侧位移明显表现出不连续性，且横缝拱坝的顺河向位移明显要大于无横缝的整体拱坝。

图10给出了顺河向位移取得最大值时的拱向应力云图。从图10可看出，由于横缝张开效应，带横缝拱坝上部拱向应力（压）明显变小。图11为顺河向位移取得最大值时的竖向应力云图。图11表明，横缝模型的梁向拉应力明显增大，明显表现出一部分荷载由拱向转化为梁向承担。

4 结束语

本研究充分发挥比例边界有限元法（SBFEM）半解析特性与网格划分灵活的优势，构建了高拱坝-地基系统非线性地震响应分析模型，并通过拉格朗日乘子法实现了对横缝接触约束的精确模拟，数值算例验证了本模型计算结果与ABAQUS计算结果的高度一致性。分析结果揭示了横缝对大坝受力状态的动态调节机制，即考虑横缝张开后，坝体拱向应力显著释放，部分拱向应力转由梁向应力承担，从而造成梁向应力明显增大。

为聚焦SBFEM在处理复杂边界接触问题上的算法优势，本文暂未引入坝体材料的本构非线性特征，后续研究将在此基础上，引入能够反映混凝土力学特性的损伤或塑性本构模型，以更全面地评估高拱坝的抗震性能。

参考文献

[1] 陈厚群， 吴胜兴， 党发宁. 高拱坝抗震安全［M］. 北京： 中国电力出版社， 2012.

[2] CHEN Houqun， WU Shengxing， DANG Faning. Seismic safety of high arch dams［M］. Beijing： China Electric Power Press， 2012．

[3] 孙立国， 李斌， 江守燕， 等. 考虑坝-基-库水相互作用的拱坝底缝开裂模拟［J］. 三峡大学学报（自然科学版）， 2024， 46（6）： 25-31.

[4] SUN Liguo， LI Bin， JIANG Shouyan， et al. Nonlinear cracking simulation of a concrete arch dam considering dam-foundation-reservoir interaction［J］. Journal of China Three Gorges University （Natural Sciences）， 2024， 46（6）： 25-31.

[5] 江德军， 黄会宝， 巨淑君， 等. 泸定Ms6. 8级地震大岗山拱坝横缝及阻尼器工作性态分析［J］. 岩土力学， 2024， 45（5）： 1527-1539.

[6] JIANG Dejun， HUANG Huibao， JU Shujun， et al. Analysis of working behavior of transverse joints and dampers of Dagangshan arch dam during Luding Ms6.8 earthquake［J］. Rock and Soil Mechanics， 2024， 45（5）： 1527-1539.

[7] 杨建， 谢伟， 张志伟. 基于多面体单元的三维ES-FEM理论研究［J］. 西北工业大学学报， 2021， 39（4）： 747-752.

[8] YANG Jian， XIE Wei， ZHANG Zhiwei. Exploring three-dimensional edge-based smoothed finite element method based on polyhedral mesh to study elastic mechanics problems［J］. Journal of Northwestern Polytechnical University， 2021， 39（4）： 747-752.

[9] 石东洋， 周钱. BBM-Burgers方程的非协调有限元方法的超收敛分析［J］. 信阳师范学院学报（自然科学版）， 2024， 37（2）： 182-189.

[10] SHI Dongyang， ZHOU Qian. Superconvergence analysis of nonconforming finite element method for BBM-Burgers equation［J］. Journal of Xinyang Normal University（Natural Science Edition）， 2024， 37（2）： 182-189.

[11] SONG Chongmin， WOLF J P. The scaled boundary finite-element method： Alias consistent infinitesimal finite-element cell method for elastodynamics［J］. Computer Methods in Applied Mechanics and Engineering， 1997， 147（3/4）： 329-355.

[12] 章鹏， 杜成斌， 赵文虎. 基于图像八叉树的三维比例边界有限元多面体网格生成算法［J］. 河海大学学报（自然科学版）， 2024， 52（1）： 46-54.

[13] ZHANG Peng， DU Chengbin， ZHAO Wenhu. Study on mesh generation algorithm of polyhedron in 3D scaled boundary finite element method based on image Octree［J］. Journal of Hohai University（Natural Sciences）， 2024， 52（1）： 46-54.

[14] 陈灯红， 梁羽翔， 刘云辉， 等. 耦合SBFEM-FEM的高拱坝系统地震响应分析［J］. 水力发电学报， 2025， 44（4）： 72-84.

[15] CHEN Denghong， LIANG Yuxiang， LIU Yunhui， et al. Seismic response analysis of high arch dam system with coupled SBFEM-FEM［J］. Journal of Hydroelectric Engineering， 2025， 44（4）： 72-84.

[16] 蒋新新， 钟红， 李云途， 等. 重力坝地震断裂的多边形比例边界有限元模型研究［J］. 水利学报， 2024， 55（1）： 115-125.

[17] JIANG Xinxin， ZHONG Hong， LI Yuntu， et al. Dynamic fracture modelling of concrete gravity dams by polygon scaled boundary finite element method under earthquake［J］. Journal of Hydraulic Engineering， 2024， 55（1）： 115-125.

[18] SONG Chongmin， OOI E T， NATARAJAN S. A review of the scaled boundary finite element method for two-dimensional linear elastic fracture mechanics［J］. Engineering Fracture Mechanics， 2018， 187： 45-73.

[19] 许贺， 邹德高， 孔宪京. 基于FEM-SBFEM的坝-库水动力耦合简化分析方法［J］. 工程力学， 2019， 36（12）： 37-43.

[20] XU He， ZOU Degao， KONG Xianjing. A simplified dam-reservoir dynamic coupling analysis method based on FEM-SBFEM［J］. Engineering Mechanics， 2019， 36（12）： 37-43.

[21] SHEN Xiaowei， DU Chengbin， JIANG Shouyan， et al. Multivariate uncertainty analysis of fracture problems through model order reduction accelerated SBFEM［J］. Applied Mathematical Modelling， 2024， 125（Part A）： 218-240.

[22] SONG Chongmin. The scaled boundary finite element method： Introduction to theory and implementation［M］. Hoboken： John Wiley & Sons， 2018.

[23] 张建伟， 刘鹏飞， 王涛， 等. 考虑横缝接触耦合坝体材料非线性下的拱坝损伤开裂分析［J］. 水力发电， 2020， 46（3）： 66-72.

[24] ZHANG Jianwei， LIU Pengfei， WANG Tao， et al. Damage cracking analysis of arch dam considering nonlinear contact between transverse joints and material nonlinearity［J］. Water Power， 2020， 46（3）： 66-72.

[25] SAPUTRA A A， EISENTRÄGER S， GRAVENKAMP H， et al. Three-dimensional image-based numerical homogenisation using octree meshes［J］. Computers & Structures， 2020， 237： 106263.

[26] ZHANG Peng， DU Chengbin， ZHAO Wenhu， et al. Dynamic crack face contact and propagation simulation based on the scaled boundary finite element method［J］. Computer Methods in Applied Mechanics and Engineering， 2021， 385： 114044.

[27] 杜成斌， 李云， 杜宁宇， 等. 横缝对高拱坝系统非线性地震响应的影响分析［J］. 三峡大学学报（自然科学版）， 2023， 45（5）： 38-44.

[28] DU Chengbin， LI Yun， DU Ningyu， et al. Effect of transverse joints on nonlinear seismic response of high arch dam system［J］. Journal of China Three Gorges University（Natural Sciences）， 2023， 45（5）： 38-44.

基金资助

国家自然科学基金项目(12202179)

南京工程学院大学生创新训练项目(202511276040)