基于COMSOL多物理场仿真计算下水利大坝静、动力特性分析研究

刘 茵

(东莞市水务工程质量安全监督站，广东 东莞 523109)

水利工程中不仅静力荷载下安全稳定性是工程设计师们持续关注的问题[1-3]，同样针对性解决水工建筑抗震设防乃是重要的研究课题[4-6]，结合静、动力分析，确保水工结构安全设计是工程运营的必要保障。刘婷婷[7]、赵建坤等[8]、李松平等[9]引入物理模型试验理论，设计有原型尺寸室内试验模型，设计开展大坝、溢洪道等水利设施的渗流试验、静力破坏实验等，为工程设计提供重要试验参数。当然物理模型试验亦可在动力分析中使用，利用振动台试验设备，设计开展坝体动荷载下模拟破坏过程，研究坝体失稳过程动力响应特征，为实际抗震设计提供参考[10，11]。模型试验成本较高，使用效率不高，而数值计算作为一种高效率手段，可借助仿真计算研究大坝、水闸、溢洪道或其他水工建筑的渗流场、静力特性、动力响应等，为工程设计方案对比优化、工程安全性检验等提供重要佐证[12-14]。本文根据水利枢纽工程中坝体基本工况，设计开展静力工况与动力荷载下水利大坝的应力、位移等参数分析，评价水工建筑的安全稳定性，为工程设计加固等提供计算依据。

1 工程概况

1.1 工程资料

地区为提升水资源利用效率，考虑拟建一水电站枢纽工程，该工程包括有上游蓄水库、拦水大坝、溢洪道、引水隧洞及其他水利设施，承担着地区防洪、蓄水、发电等作用，可降低区域枯水季缺水率10%以上，可供水量超过30万 m3/d，上游水库正常蓄水库容为1.5亿 m3，水库含沙量为7.8 kg/m3，年可发电量为90亿 kW·h，其安全稳定性乃是区域内水利安全的重要保障。溢洪道设定堰顶高程为945 m，具有多孔式水闸作为流量设施，每个孔直径均为2.2 m，最大泄流量设计为5 800 m3/s，水闸采用预应力闸墩作为结构加固措施，墩厚直径为1.2 m，预应力锚索最大张拉吨位为6 580 kN，确保水闸受水力冲击动力下闸室结构静、动力特性满足安全设计要求；另溢洪道泄流下游设计有挡土边墙作为水利冲刷控制设施，最大位移不超过2 mm，结构张拉应力稳定，所连接的消能池深度为1.4 m，采用变坡坎设计，减弱水力侵蚀作用。拦水大坝顶高程为948 m，坝高最高处为125 m，轴线长度为792 m，设计以止水面板作为坝身防渗系统，另由于坝基所在场地含有黏土层，沉降变形较大，因而设计有防渗墙与垫层，其中防渗墙底部位于基岩层，深度为3.2 m，厚度为0.6 m，垫层铺设厚度为0.5 m，有效降低黏土层沉降变形，坝身最大渗透坡降不超过0.22，渗透稳定性较佳。根据工程场地勘察，该枢纽工程设计为Ⅶ级烈度，防洪标准按照洪峰流量5 650 m3/s设计，拦水大坝采用拱坝形式设计，分左、右拱形，厚度最大分别为52 m、42 m，曲率半径最大分别为233 m、135 m，基岩层为半风化灰岩，钻孔资料显示完整性较佳，室内测试强度超过55 MPa，覆盖层软弱层中以黏土层为设计着重关注土层，其他土层含水量适中，沉降变形较小，可塑性较佳，颗粒粒径以0.25～4.75 mm为主。该拦水大坝在静、动力荷载下的稳定性乃是该枢纽工程的重点考虑因素，工程设计部门考虑利用COMSOL多物理场仿真平台开展水利大坝静、动力分析，此对预判或评价水利工程安全有效运营具有重要意义。

图1 大坝有限元模型

1.2 工程建模

利用COMSOL多物理场仿真平台建立拱坝有限元模型，并进行精密网格划分[15，16]，共获得网格单元125 633个，节点数98 246个，所建立模型如图1所示。静力荷载下坝体模型顶部为单向约束，而模型底部为多向约束，所有荷载均直接施加在模型边界上。动力荷载下采用反应谱叠加法进行求解，施加的反应谱为工程场地可靠性加速度标准谱，峰值加速度为2.5 m/s2，其反应谱如图2所示，动力荷载下另还包括有动水压力Pw，其计算式如式(1)所示。另为分析方便，本文计算模型中X、Y、Z正向分别为坝体右岸向、顺水流向、模型向上方向，模型计算范围包括有三向分别为760 m、462 m、155 m。

(1)

式中：α指动水压力计算系数；ρw指液体密度；H0指水位深度；h指动水压力计算深度。

图2 加速度反应谱

2 大坝静力特性分析

本文根据静力荷载下施加应力不同，共划分出5种工况，其中A、B、C工况均无温度热应力与泥沙沉降应力两种荷载，此三种工况分别对应水库正常蓄水期、完建期及死水位，三种工况上、下游水位分别946m/719.5m、940m/704m、922m/704m，另有D工况包括有温度热应力荷载，E工况包括有泥沙沉降应力荷载，D、E工况中水位与正常蓄水期一致，分别对5种工况下大坝静力特性展开计算分析。

2.1 应力特征

根据对五种工况应力特征计算，获得各工况中拉、压应力变化特征，如图3所示。从图中可看出，五种工况中拉应力最大为B工况，达6.78 MPa，而D、E工况中最大拉应力基本相近，相比B工况分别降低了78.9%、79%，C工况中拉应力最低，仅为B工况最大拉应力的16.6%，A工况最大拉应力较之B工况差距最小，其减少幅度仅为16.2%。分析认为，完建期水利大坝张拉应力最大，而正常蓄水期后，由于水位上升，对坝体产生静力弯矩，此亦为其拉应力产生来源，且正常蓄水期拉应力的产生与完建期基本相近；温度热荷载与泥沙沉降应力两者施加在计算模型中后，产生的张拉应力效果基本一致，表明两者荷载对大坝张拉应力的影响基本一致，在计算张拉应力时可作为一类考虑。当水位处于最低时，结构张拉应力较低，且不超过结构材料允许拉应力设计值。

图3 拉、压应力变化特征

与拉应力呈对比的是，压应力在A～C工况中基本接近，差距幅度不超过0.4%，均稳定在9.6 MPa，D、E工况中最大压应力高于A～C工况，分别为14.85 MPa、13.75 MPa，相比A～C工况稳定压应力要增大了54.7%、43.2%，即温度热荷载与泥沙沉降应力对结构压应力具有显著影响，且为正向促进效应。分析表明，大坝在运行过程中不可避免会出现温度热荷载不均的现象，坝体上下变形发生不一致，导致产生剪力弯矩，此一定程度上促生了坝体上最大压应力的增高现象；而泥沙沉降应力的施加，最直接的结果即是增大了坝体自重，最终表现压应力增大的结果。从最大压应力位于区域可知，A～C工况中最大压应力均位于坝体下游面坝肩处，此主要为上游水压力乃是压应力产生重要来源，故而表现在下游面上压应力较大；D、E工况中最大压应力分别位于坝中心轴线处，两者分布区域基本类似，但E工况更靠近坝址。

图4 各向最大位移值变化特征

2.2 位移分析

根据各工况中各向位移计算，获得不同工况下各向最大位移值变化特征，如图4所示。从图中可知，各工况中除D工况外，其他工况X向位移方向均指向坝体左岸，量值上以B工况最大，达35.4 mm，A、E工况最大位移较之前者分别降低了19.5%、13.4%，从X向位移特征可知，D工况中温度热荷载的施加，可导致坝体上位移方向的改变，且其量值亦较低，相比B、C、E工况中X向位移分别降低了80.6%、73%、77.6%。在Y向位移中，各工况位移方向均为一致，顺水流向发展，最大位移乃是D工况，达45.68 mm，而位移最低为C工况，仅为12.56 mm，且A～C工况中Y向位移随上游水位增大而递增，其中A、B工况中Y向位移较之C工况分别增长了1.83倍、72.5%；笔者认为，Y向位移与上游水位密切相关，亦与温度荷载、泥沙沉降应力等附加荷载有关，当坝体负荷愈大，则坝体Y向位移愈显著。各工况中Z向位移方向均为向下沉降，最大值为E工况，达21.1 mm，而D工况中Z向位移为最小，仅为E工况中的31.2%，另一方面A～C工况中Z向位移基本接近，稳定在14.66 mm，上游水位变化，对A～C工况中Z向位移无显著影响，而温度热荷载对Z向沉降位移具有削弱作用，不可忽视泥沙沉降应力可增大坝体沉降，相比正常蓄水期沉降值增大了43.9%，因而坝体沉降应考虑全面，减少由于大坝沉降而引起的工程失稳可能性。

3 大坝动力特性分析

为分析大坝动力特性差异性，对计算模型分别施加三个不同峰值加速度反应谱[17]，其他参数保持一致，三个峰值加速度分别为2.5 m/s2(1#方案)、5m/s2(2#方案)、7.5 m/s2(3#方案)，以此计算地震动荷载差异性下动力响应特征。

3.1 模态分析

根据模态分析计算获得各阶次大坝自振频率变化特征，如图5所示。根据图中自振频率变化可知，各峰值加速度计算方案中自振频率随阶次均为递增关系，在1#方案中第1阶次自振频率为1.26 Hz，而第5阶次、10阶次、15阶次自振频率较之前者分别递增了57.1%、165.1%、287.3%，在1#方案中平均每阶次自振频率增幅为7.7%。而在2#方案中，第5阶次、10阶次、15阶次自振频率较之第1阶次分别递增了66.6%、171.8%、289.1%，当递增1阶次，自振频率可增长9.7%；同样在3#方案中每阶次的增幅为11.2%，分析表明峰值加速度愈大，各阶次自振频率增长幅度愈快，但各计算方案中振型特征基本类似，1～5阶均为平动，而6～10阶次为对称性振型，10～15阶次为反对称振型，16～20阶次保持为组合振型特征。

图5 各阶次大坝自振频率变化特征

3.2 动力响应特征

为分析三个计算方案下坝体动力响应特征，设计针对坝体上、下游面上各特征点应力、位移响应值展开计算，各特征点分布如图6所示。

图6 大坝各特征点分布图

根据对各特征点各向最大位移计算，获得不同方案下各向位移对比特征，如图7所示。从图中可看出，峰值加速度与各向最大位移均为正相关关系，峰值加速度愈大，则各向最大位移愈大，在X向中1#方案第4个特征点位移为2.3 mm，而在2#、3#方案中同一特征点最大位移较前者分别增大了1.05倍、3.11倍，在各特征点中， 3#方案相比1#、2#方案的X向位移差距幅度分别为2.7倍～3.4倍、63.4%～1.14倍，峰值加速度增大2.5 m/s2，各特征点中X向位移平均增幅为1.95倍。Y、Z向位移水平与峰值加速度关系与X向类似，而Y、Z向位移在峰值加速度2.5 m/s2增长量下，平均增幅分别为62.8%、1.89倍，表明X向位移受峰值加速度影响更为敏感。

图7 不同方案下各向位移对比特征

图8 各特征点拉应力受峰值加速度影响变化特征

同理，根据计算可获得各特征点最大拉应力受峰值加速度影响变化特征，如图8所示。当动荷载加速度峰值愈大，则拉应力响应值愈高，且均超过结构材料安全允许值，在第5个特征点中，1#方案最大拉应力为2.7 MPa，而2#、3#方案相比之分别增大了18.5%、66.7%，在所有20个特征点中，2#方案相比1#方案的递增幅度平均为51.1%，而3#方案相比2#方案的递增幅度平均为1.93倍，峰值加速度增大2.5 m/s2，大坝拉应力平均提高了1.1倍。由此可见，地震动荷载的等级对大坝等工程抗震设防具有重要参考价值。

4 结语

本文主要得到以下几点结论：

(1)静力荷载下，拉应力最大工况为完建期，达6.78 MPa，死水位工况拉应力最低，此两工况中压应力基本接近，稳定在9.6 MPa，温度热荷载与泥沙沉降应力两工况中的拉应力效果相近，但对压应力具有正向促进效应。

(2)静力分析下X向位移以完建期最大，指向左岸；Y向位移最大为温度热荷载工况，达45.68 mm，Z向位移以泥沙沉降应力工况最大，无温度荷载与泥沙沉降荷载的三个工况中Z向位移相近，稳定在14.66 mm。

(3)模态分析中自振频率随阶次均为递增，峰值加速度递增，自振频率愈大，1#、2#、3#方案中每阶次自振频率递增幅度分别为7.7%、9.7%、11.2%。

(4)峰值加速度与各向最大位移、拉应力均为正相关关系，峰值加速度增大2.5 m/s2，各特征点中X、Y、Z向位移平均增幅分别为1.95倍、62.8%、1.89倍，而拉应力平均提高了1.1倍。