基于COMSOL Multiphysical计算下水库土石坝动力响应特征研究

王 霄，谷 静，齐春舫，李 磊

(淮安市水利勘测设计研究院有限公司，江苏 淮安 223005)

1 概述

水利工程中动力抗震特性的研究很有必要，对提升水工结构在区域断裂带上的安全性具有较大帮助。研究水利设施的动力抗震特性与动力响应息息相关，因而主要借助研究水工结构的动力响应特征来开展相关抗震设计[1- 3]。王海波等[4]、程新俊等[5]、杜玉涛等[6]采用振动台试验手段，研究了水工建筑模型在地震台试验条件下动力破坏状态，为水利结构抗震设计提供重要依据。从现场监测角度考虑，以现场监测数据开展地震动力响应关联性研究实质上更有意义[7- 8]，或借助微震等监测仪器[9]，均是较好研究水工结构动力响应的手段。作为水工建筑的动力响应，利用仿真计算平台，施加动荷载至水工模型，可计算获得在不同运营工况下水工结构的动力响应解，为实际水利工程抗震设计提供重要参考[10- 12]。本文利用COMSOL Multiphysical仿真平台，研究了水库土石坝在人工合成地震波下动力响应特征，为实际分析抗震特性提供参考。

2 工程计算概况

2.1 工程资料

某水库位于华中地区，北侧高程最大，库容152万m3，承担着区域内蓄水、防洪、水资源调度的重要作用，水资源供应主要面向农业用水与工业用水，建设有长度为65km的输水渠道，水渠工程全长均采用格宾石笼作为衬砌形式，减少水流损耗，进而提升枯水期水资源灌溉效率的稳定性；工业用水主要通过水库输水隧洞设施传送至泵站，进而调度至各个用水项目。该水库年可供工业用水量超过50万m3，隧洞渠首流量设计为0.6m3/s。为提升水库工程运营安全状态，考虑对水库工程进行检修加固，而该水库工程包括有土石坝体、泄洪闸、消能池等水工设施，其中坝顶高程为56.5m，顶宽度为5m，迎水侧与背水侧坡度均设置为1/3，采用分层堆筑，确保坝体沉降变形值处于安全允许值范围内，坝体全场为1400m。坝身采用混凝土浇筑形成网格，结合植被护坡形式，提升土石坝边坡稳定性，这样在设计之初不仅坝坡安全性得到考虑，坝体防渗同样被重点关注。采用55cm厚度的防渗墙插入基岩层1.5m，并在坝趾、坝肩等重点区域设置有土工膜等防渗材料，在坝基处铺设有20cm厚的碎石垫层，增强坝体稳定性与防渗性。根据已有设计资料，该坝体在静力荷载运营状态下稳定性均较好，且最大渗透坡降不超过0.2，不论是静力稳定性或是渗流稳定性，水库土石坝体均处于安全正常运营。工程设计部门主要针对性解决在动力荷载影响下水库土石坝动力响应特征，探讨水库土石坝当前动力抗震特性。

为确保仿真计算准确性，应对水库土石坝所在工程现场进行地质踏勘分析。由现场踏勘得知，该工程表面覆盖有第四系种植土，松散性较大，颗粒粒径以1mm左右为主，磨圆度较差，厚度为0.3m；在基岩层上覆存在有一层壤土，含砂量较大，承载力适中，输水渠道衬砌结构持力层即为该土层，含水量较小，干燥状态下渗透性较好；基岩层包括有花岗闪长岩与辉绿岩，其中花岗闪长岩按照风化程度划分出两个岩层，上覆岩层乃是全风化状态，厚度为2.5m，颗粒粒径以粗颗粒为主；强风化花岗闪长岩局部夹有破碎带，夹层主要以砂砾石等为主；辉绿岩层厚度为8.7m，强度较大，室内测试单轴抗压强度超过55MPa，现场钻孔取样亦可看出该岩样表面无显著性孔隙，室内测试静水压力下渗透率最低仅为10-21m2，孔隙度最低仅为1%。土石坝体现场岩土层分布如图1所示。

图1 岩土层剖面图

2.2 模型建立与约束荷载

根据该土石坝实际岩土层以及堆筑材料，采用COMSOL仿真平台建立计算模型。为方便计算，本模型计算动力影响范围为基础深度20m、坝体轴线长度10m，以六面体作为基本单元网格[13- 14]，每个单元网格尺寸与输入动荷载频率有关，其关系式如(1)式所示[15- 16]，模型中网格最大尺寸为边长1.5m，经单元网格划分，共获得8286个单元体，节点数12892个，所建立的计算模型如图2所示。仿真平台中所输入的岩土体计算参数均参照室内土工试验实测值，其中壤土渗透系数取值为8.5×10-5cm/s。

图2 土石坝计算模型

(1)

式中，Δl—单元网格尺寸；λ—地震波波长。

动荷载约束则采用人工合成地震波的方式施加至模型中，以土石坝所在场地实际特征周期参数及人工设计的反应谱加速度峰值等，按照式(2)计算获得人工合成地震波反应谱计算曲线，如图3所示。但由于反应谱加速度峰值与地震烈度有关，当取不同加速度峰值时，所获得反应谱势必具有显著差异，因而本文分别设计3个峰值加速度值，分别为0.1g(A方案)、0.2g(B方案)、0.4g(C方案)，得到的人工合成地震波反应谱如图4所示，得到各研究方案下动力荷载速度时程与功率谱曲线如图5所示。模型阻尼参照已有工程资料设定为0.05；另本文计算方案中加速度均为水平方向，竖向加速度峰值按照66.7%的水平加速度进行换算。

图3 人工合成地震波

图4 人工合成地震波反应谱

图5 动力荷载速度时程与功率谱曲线

(2)

式中，βT—加速度反应谱计算值；βmax—加速度峰值；Tg—特征周期；T—持续时间。

3 土石坝加速度响应特征

经仿真计算获得不同动荷载峰值加速度研究方案下的土石坝加速度响应特征，如图6所示。从图6中可看出，土石坝在各方案下的水平加速度动力响应最大值均出现在坝顶，其中A方案最大水平加速度为0.28m/s2，坝底处水平加速度乃是该研究方案下最低值，其中A方案坝顶与坝底加速度最大差幅可达71.4%，而在B、C方案中差幅分别为

图6 土石坝水平加速度响应特征

56.3%、62.8%，三个方案中坝顶与坝底加速度影响差距基本稳定在50%～70%；从整体加速度响应分布特征可知，坝底至坝顶，加速度响应分布呈阶梯递增态势，笔者认为，在动力地震波作用下，坝体结构受地震荷载的放大效应影响，在坝顶处集中受到较大地震破坏，其乃是地震破坏中首当其冲的危险结构区域，故为确保土石坝结构抗震安全性，应着重加固坝体结构。

同理可获得竖向加速度在土石坝上分布变化曲线，如图7所示。从图7中可知，坝体上不同特征部位的竖向加速度响应特征具有一定差异，其中以下游坝坡加速度放大倍数最大，在A方案中，坝体相对高度0.4m时，下游坝坡的竖向加速度放大倍数为2.65，而在相同高度处的坝轴线、上游坝坡竖向加速度放大倍数相比前者分别降低了9.4%、24.5%，即下游坝坡受到的地震荷载放大破坏效应更剧烈；另一方面，当输入动荷载加速度峰值增大至0.2g(B方案)，前述相同位置处的幅度差异变为9.5%、24.4%，表明输入地震荷载加速度峰值并不影响坝体特征部位处加速度放大倍数幅度差异。从三个方案整体加速度响应值可看出，地震荷载输入加速度峰值愈大，则坝体特征部位处加速度放大效应有所减少，当均为坝体相对高度0.6m时，输入的地震荷载峰值加速度为0.1g，即A方案时，下游坝坡加速度响应值放大倍数为2.8，而B、C方案相应位置的加速度响应值放大倍数相比前者分别降低了21.4%、37.5%。再一个加速度响应值放大倍数与坝体相对高度关系可知，以接近坝顶处区域的放大倍数增长速率为最大，相对高度0.8～1m区域内，高度增大0.2m时，A方案下游坝坡平均增长幅度为6.3%，而在坝体至坝中部区域增长速率较前者显著，坝底至坝中部相对高度0.6m时，高度增大0.2m，平均增长幅度可达116.6%；由此表明，竖向加速响应值放大倍数在坝体高度上具有增长差异性，接近坝体区域处为典型破坏区域，其内因即是放大倍数与增长效应均达到饱满状态，增长速率较小。

图7 土石坝竖向加速度响应特征

4 土石坝动力响应特征

4.1 动应力响应特征

根据仿真计算获得土石坝动力响应特征，包括有坝体动应力与动位移，图8为三个研究方案下坝体动应力响应的最大主应力分布特征。从最大主应力等值线分布可看出，压应力乃是坝体内部主导动应力，坝体底部压应力最大，其中三个方案中最大压应力均为0.6MPa左右，但B方案最大压应力出现在持续时间20s时；当持续时间变化时，最大主应力分布特征无显著变化，不论输入地震动加速度何值，变化均较小，表明地震动荷载对坝体的最大主应力影响较小。

图8 坝体动应力响应最大主应力分布特征

相比最大主应力特征，图9为三个研究方案下典型持续时间的最小主应力分布特征。从图9中可看出，在动荷载后期具有一定拉应力外，其他特征部位处的的地震动荷载均以压应力为主，最大压应力均为0.2MPa左右，但C方案中相同持续时间时最大压应力相比A、B方案下均有一定程度降幅，分别为0.19MPa、0.16MPa；分析认为，坝体拉应力主要出现在动荷载施加后期，该坝体全过程主要以受压为主，最大拉应力量值较小，处于堆筑料拉应力安全允许范围内。从各方案不同持续时间下最小主应力分布特征亦可知，A、B方案的最小主应力基本均处于不变状态，随持续时间无显著性变化，表明地震动荷载在较低地震烈度(较小峰值加速度输入值)时，地震动荷载对最小主应力同样影响较小，当加速度峰值输入为0.4g时，最小主应力随持续时间有一定变化，表明高地震烈度中不同持续时间下最小主应力最大值会降低，拉应力会上涨。

图9 坝体动应力响应最小主应力分布特征

4.2 动位移特征

与动应力计算同理，获得水库土石坝水平动位移特征，图10为典型持续时间下水平动位移等值线分布图。当持续时间增大，水平动位移递增，A方案中持续时间5s时最大水平动位移为6.3cm，而持续时间为20s的最大水平动位移相比前者增大了98.4%；从坝体相对高度方向上，坝底至坝顶呈递增态势。对比三个研究方案可知，相同持续时间下地震动荷载峰值加速度愈大，则水平动位移愈大，同为持续时间5s时，B、C方案最大水平动位移相比之下A方案增大了39.7%、98.4%，表明地震荷载造成的破坏烈度更大，则坝体水平方向位移愈大，且方向逐渐由初始指向上游，在后续持续时间(20s)时指向下游。

图10 坝体水平动位移分布特征

同理可获得坝体竖向动位移分布特征，如图11所示，相比水平位移，竖向位移均显著低于前者，B方案中持续时间15s时的最大竖向动位移为5.2cm，而同一时间下的最大水平动位移乃是前者的3.2倍。对比同一方案中各持续时间下竖向动位移发展趋势可知，竖向位移量值在逐渐增大，且由沉降方向逐渐蔓延至垂直向上方向，在C方案中持续时间5s时的最大竖向动位移为8.5cm，方向为向下沉降，而持续时间增大至15、20s后的竖向动位移最大量值相比前者分别增大了2.4、4.1倍，方向均为向上，即坝体已产生变形破坏，分析表明坝体竖向方向上动位移响应破坏主要出现在持续时间后期，前期主要以沉降变形为主。从不同方案下竖向动位移对比可知，当输入地震动荷载峰值加速度愈大，则坝体承受的变形位移量值愈大，A方案中持续时间20s时的最大竖向位移为8.5cm(向下)，而B方案中最大竖向位移相比A方案增大了111.8%(向下)，C方案相比B方案又增大了94.4%(向上)。

图11 坝体竖向动位移分布特征

5 结论

(1)坝体水平加速度最大响应值位于坝顶，坝顶与坝底加速度差距基本稳定在50%～70%；下游坝坡竖向加速度响应值最大，峰值加速度愈大，则加速度放大倍数愈小。

(2)坝体最大主应力以受压为主，最大压应力为0.6MPa左右，最小主应力中拉应力出现在地震荷载后期。

(3)水平动位移随持续时间与坝体高度为递增，地震荷载峰值加速度愈大，则水平动位移愈大，地震荷载后期，水平位移方向为指向下游。

(4)竖向动位移显著低于水平动位移，同一时间下的B方案最大水平动位移乃是竖向位移的3.2倍；地震荷载峰值加速度愈大，竖向动位移量值愈大。