引言

水工结构的基础是1个半无限体，其真实的边界条件是基础位移以及应力在无穷远处不受上部结构的影响。在土石坝数值分析中，存在对仿真模型无限域或半无限域的边界条件的简化和截取计算范围的问题，通常采取的办法是人为截取一定范围宽度的地基边界，将无限域变成有限域。特别是对于三维模型而言，边界截取范围太小，影响计算精度；截取范围太大，计算量大且在很大程度上浪费了时间和资源。旁多水电站工程处于海拔4000m以上的西藏高原复杂地质条件地区，考虑在地震荷载作用下坝址区域基础为深覆盖层，防渗墙采用悬挂式的特殊结构条件下，土石坝结构动力特性和抗震性能受边界条件的影响情况亟待研究。目前，国内外尚未见有专门针对高原地区复杂地质条件下采用无限元法对土石坝深厚覆盖层悬挂式防渗墙特殊结构下抗震性能影响的研究。

20世纪中后期，计算机与数值分析计算的出现和迅速发展，为解决无限域或半无限域类工程分析提供了强有力的工具和方法，从而为使用有限单元法和无限元法求解非线性工程实际问题成为可能。Ungless(1973)提出无限元的概念，Bettess等(1977)在研究流体波动时提出了无限元的应用，Chow等(1981)提出谐振无限固体单元，将无限元的研究引入固体中波的传播分析，Zhang等(1987)将这一方法应用于地基-结构动力相互作用分析中。20世纪90年代后，国外研究大都集中于无限元理论基础分析之上(Khalili等，1997；Noorzaie等，2002；Yerli等，2003)，仅有少部分学者将无限元的理论引用于地基-结构相互作用(SSI)实际模型数值求解中。周小溪等(2013)通过对深覆盖层土石坝地震响应特征分析，采用人工截取地基范围和施加边界条件，认为地震荷载对防渗墙无明显影响。Han等(2016)通过现场测定深覆盖层工程力学性质和压力分布规律，研究发现通过实验室很难确定力学性质和分布规律。西藏实际工程因地质条件复杂，覆盖层深、防渗墙为悬挂式。上述文献中大部分的分析模型为一维或二维，不能反映实际工程的真实受力情况，计算仿真三维模型结构在一定程度上存在缺陷，在反映实际工程地震动作用下的结构动力特性和抗震性能可能存在较大误差。仅有部分学者建立模型为三维，其分析的实际工程地基地质条件简单。随着科研工作的进一步深入，高原地区复杂地质条件下工程建设与抗震安全越来越受到专家及学者的关注，特别是西藏高原地区基础为深覆盖层的情况。赵崇斌(2010)改进了Chow等(1981)的二维无限元，提出了动力映射无限单元；向前(2004)和孙萍等(2007)采用无限元法研究了波在无限介质中能量逐渐弥散的现象；燕柳斌(2004)、陈跃庆等(2006)利用无限元方法模拟结构与地基相互作用以及不同土性地基动力相互作用具有相似的规律；董光辉(2011)、温立峰等(2015)和余翔等(2018)通过对深覆盖土石坝非线性动力分析，认为在地震动响应计算中应考虑辐射阻尼的影响，尽量改进边界条件和计算方法。杨正权等(2016)对超深覆盖层土石坝的边界处理方法进行了深入研究，将底部边界取到基岩面、侧向为3—5倍坝高，并采用固定边界的处理方式。

本文针对西藏高原地区旁多水电站土石坝实际工程，考虑坝址区域地基为深覆盖层情况，使用有限单元法与无限元法进行非线性动力分析，分析了土石坝在有限单元法和采用映射无限元离散远场地基下的加速度、动位移等动力响应之间的差异，讨论了2种方法在高原地区复杂地质条件下对分析深覆盖层悬挂式防渗墙土石坝结构动力特性和抗震性能的影响。

1.   模型及参数

西藏旁多水电站地处拉萨河流域中游，坝址位于西藏自治区林周县旁多乡下游1.5km，工程位于海拔4000m以上，坝基覆盖层深度达200m，整个枢纽坝址区具有深覆盖层、高地震烈度等特点。西藏位于印度洋板块和亚欧板块交界处，因印度洋板块的挤压抬升作用，导致西藏地区海拔越来越高，由于该地区自古就是板块运动活跃的地方，因此地震发生频繁。坝基覆盖层由含混合土碎（块）石、冲积卵石混合土（Q₄^al）、冲积卵石混合土（Q₃^al）和冰水积卵石混合土（Q₂^fgl）组成。大坝结构为碾压式沥青混凝土心墙砂砾石坝，沥青混凝土心墙顶部厚度0.70m，底部最大厚度1.00m（坝高的1/70），最大坝高72.30m，坝顶长1052.00m，坝顶宽12.00m，大坝最大断面见图 1。

图 1  坝体剖面尺寸

Figure 1.  Detail size parameters of the dam body

下载: 全尺寸图片 幻灯片

1.1   计算模型

计算工况为正常蓄水期加地震动分别计算有限单元法3倍坝高、5倍坝高和无限元法3倍坝高，坝体、坝基均采用八节点六面体单元进行离散，有限单元类型为三维实体单元（C3D8），无限单元类型为三维无限元实体单元（CIN3D8），对坝体分别进行有限单元法和无限元法计算分析，三维计算模型共55483个单元，60778个节点。上、下游各取3倍坝高，左右坝肩向外各取1.5倍坝高，建基面以下分别取3倍、5倍坝高。有限单元法底部输入地震动，侧向约束垂直于各面的自由度；无限元法底部输入地震动，侧向自由。三维模型计算网格图见图 2。

图 2  三维计算网格模型

Figure 2.  Calculation mesh of 3D model

下载: 全尺寸图片 幻灯片

1.2   材料参数

目前，土石坝的地震动响应分析分为等效线性和非线性两大类。后者往往适用于一维场地反应分析，在土石坝三维动力分析中应用较少且参数难以准确确定，且受计算时间效率的限制。本文通过前者迭代的方法来近似地反映土体的非线性，材料等效线性模型计算参数见表 1。在土石坝粘弹性边界条件中，土的剪切模量G和阻尼比λ是剪切应变$\overline {{\gamma _{\rm{d}}}} $的函数，引用沈珠江等（1996）提出的公式：

表 1  材料等效线性模型参数

Table 1.  Parameters of material equivalent linear model

材料	k1	k2	n	λmax
坝壳料	20.43	2336	0.268	0.19
心墙	21.06	1106	0.556	0.25
地基	15.66	450	0.500	0.20

下载: 导出CSV 

| 显示表格

式中：P_a为大气压强；${{\sigma '}_3}$是围压；k₁、k₂和n是由试验确定的材料参数；$\overline {{\gamma _{\rm{d}}}} $是归一化的剪应变，可根据地震过程中的最大动剪应变${\gamma _{{\rm{dmax}}}}$按（3）式计算：

通过对公式（1）—（3）进行迭代，动剪切模量比和阻尼比随之更新，对迭代计算中最大剪应变随深度的分布图观察迭代效果，一般3—4次可达到收敛要求。

2.   三维映射无限元原理

采用单向映射无限元模拟地基无限域问题，在基岩耦合面输入三向地震动，大坝近域部分地基用有限元离散，远场地基用映射无限元离散，八节点三维映射无限元模型如图 3所示。节点1—4与有限元节点相耦合，节点5—8为无限元中间节点，其余4个节点在无穷远处。

图 3  三维八节点映射无限元

Figure 3.  3D infinite element of 8-node mapping

下载: 全尺寸图片 幻灯片

半空间无限域共有5个面需要用到单向映射无限元，节点编号顺序不同所对应的形函数和映射函数不同，其中形函数N_i为：

单向映射无限元映射函数M_i表示为：

单元刚度矩阵按式（8）计算：

式中，ξ、η、ζ为局部坐标；[D]为弹性矩阵；[B]为刚度矩阵，[B]=[B₁，B₂，B₃，B₄，B₅，B₆，B₇，B₈]。

3.   动力响应分析

大坝场址区基本烈度为Ⅷ度，参考坝址工程地勘资料及区域地质条件，地震动输入采用迁安波，时程采样间隔0.01s，持时为15s，峰值为1.32m/s²。横河向、顺河向和竖直向输入加速度之比为3^:3^:2，计算中输入地震波如图 4所示。通过静力计算提取初始应力，假定最大剪切模量初次迭代数值，结合公式（1）、（2）和（3）对最大剪切模量进行计算、提取，不同材料区域的参数为等效线性模型参数。为简化迭代计算过程，采用FORTRAN编程对最大动剪切模量比和阻尼比进行计算。从各迭代计算中最大剪应变随高度的分布情况（图 5）看出，经过4次迭代，第3次和第4次迭代基本吻合，可以利用第4次迭代结果获得的动剪切模量比和阻尼比进行正式的计算。

图 4  地震加速度时程

Figure 4.  Acceleration time history of earthquake

下载: 全尺寸图片 幻灯片

图 5  各迭代计算中最大剪应变随高度的分布

Figure 5.  Distribution of maximum shear strain with depth

下载: 全尺寸图片 幻灯片

3.1   动位移响应

图 6为最大动位移沿最大断面中轴线的分布，由图可知，顺河向最大动位移随着基岩厚度增加而减小，有限单元法（3倍）、（5倍）和无限元法（3倍）最大动位移分别为4.03cm、3.88cm和3.71cm。通过数据分析可见，远场地基采用无限元离散下的顺河向最大动位移最小；竖向最大动位移最大响应在坝底处，有限单元法（3倍）曲线拟合和另外两者在坝体1/3以上相差较大，无限元法（3倍）和有限单元法（5倍）拟合系数相对更高。此外，通过数据分析可知在蓄水作用的影响下，上游面板堆石体内顺河向最大动位移在同一高程处比下游堆石体内大；随着坝体高程的增加顺河向动位移和坝轴向动位移增大。经对3个几何模型的计算结果进行分析比较，结果表明：考虑地基远场弹性作用对最大动位移的影响是不可忽略的，特别是针对复杂地质条件深覆盖层情况，在相同坝基深度情况下，无限元法较有限单元法更具有一定的合理性和精确性。

图 6  最大动位移沿坝高的分布

Figure 6.  Max-dynamic displacement distribution of dam

下载: 全尺寸图片 幻灯片

图 7为坝顶典型节点动位移时程曲线。顺河向有限单元法（3倍）、（5倍）和无限元法（3倍）动位移最大幅值分别为0.95cm、1.16cm、1.02cm；顺河向和竖向动位移响应时程曲线无明显差异。动位移幅值较大的时间均出现在3—6.5s，但由于阻尼作用，坝顶动位移时程曲线与输入地震波相比均呈现滞后效应。由图 7（b）可知，顺河向动位移有限单元法（3倍）时程曲线幅值明显较高，有限单元法（5倍）和无限元法（3倍）动位移时程曲线幅值小、拟合程度较高。

图 7  坝顶动位移时程曲线

Figure 7.  Time history of dynamic displacement at top of dam

下载: 全尺寸图片 幻灯片

3.2   加速度响应

图 8为坝顶典型节点加速度响应时程曲线。顺河向有限单元法（3倍）、（5倍）和无限元法（3倍）坝顶加速度最大幅值分别为2.06m/s²、2.40m/s²和2.52m/s²，放大系数分别为1.56、1.82和1.84；坝轴向加速度最大幅值分别为1.85m/s²、2.15m/s²和1.87m/s²。通过数据分析可知，无限元法（3倍）与有限单元法（5倍）加速度时程曲线幅值拟合系数较高，且加速度响应时程曲线幅值（不包括最大幅值）均小于有限单元法（3倍）。人工截取地基范围和施加边界条件存在一定的缺陷，而无限元法弥补了地基弹性对地震动影响的问题。加速度幅值较大的时间均出现在3—6.5s，但由于阻尼作用，坝顶顺河向加速度时程曲线与输入地震动相比均呈现滞后的现象。顺河向和竖向加速度均随坝高的增加呈递增趋势，主要以顺河向加速度为主。

图 8  坝顶加速度响应时程曲线

Figure 8.  Time history of peak acceleration at top of dam

下载: 全尺寸图片 幻灯片

图 9为最大加速度放大系数沿坝体最大断面中轴线高度的分布。顺河向最大加速度有限单元法（3倍）、（5倍）和无限元法（3倍）最大值分别为2.20m/s²、2.40m/s²和2.52m/s²，放大系数分别为1.67、1.81和1.91。从坝轴向最大加速度放大系数分布看，无限元法（3倍）和有限单元法（5倍）中下部有一定的差异，但并不明显。最大加速度随坝体高程的增加而增大，最大加速度响应发生在坝顶靠近左岸处。从顺河向放大系数看，无限元法（3倍）放大系数沿坝高的分布与有限单元法（5倍）拟合程度较高，有限单元法（3倍）坝底处放大系数比另外两者大，但坝顶处相反。从竖向放大系数分布上看，无限元法（3倍）与有限单元法（5倍）差异不明显，但有限单元法（3倍）分布曲线相对另外两者曲线低，计算精度低。

图 9  最大加速度放大系数沿坝高的分布

Figure 9.  Maximum acceleration distribution of dam

下载: 全尺寸图片 幻灯片

4.   结论

因西藏高原旁多水电站处于地质条件复杂地区，考虑坝址区域基础为深覆盖层和悬挂式防渗墙的特殊结构条件，通过映射无限元模拟远场地基的手段，对西藏高原地区实际工程在深厚覆盖层坝体-坝基相互作用下的地震动响应进行了分析，探讨了无限元法在高原地区复杂地质条件下土石坝结构动力特性和抗震性能分析中的应用。计算结果表明：采用3倍坝基的无限元法和采用5倍坝基的有限单元法得到的计算结果基本一致。因此，在相同坝基深度情况下，无限元法能够更真实反应坝体-坝基相互作用下结构动力特性和抗震性能，表明在西藏高原地区深覆盖层悬挂式防渗墙中采用无限元法对土石坝非线性动力分析切实可行，为今后西藏高原地区实际工程中分析具有相似结构特点的大坝结构动力特性和抗震性能提供了分析方法和可行的分析手段。