核电站防波堤地震动力响应及破坏机理分析

杨 勋，楼云锋，余克勤，葛鸿辉，金先龙

(1.上海交通大学 机械系统与振动国家重点实验室，上海 200240;2.上海交通大学 机械与动力工程学院，上海 200240;3.上海核工程研究设计院，上海 200233)

福岛核事故发生之后，核电站设施在地震中的运行安全问题受到了世界范围内的关注。防波堤是核电站海防工程的重要组成部分，用以保护整个核电厂区免受海潮波浪影响，地震时一旦发生破坏，造成的危害不堪设想。因此，深入开展核电站防波堤的抗震性能研究，准确把握其在地震中的动力特性及破坏机理具有重要意义［1－2］。

目前防波堤地震响应的有限元仿真多数采用基于等价粘弹性的材料本构进行等价线性分析，该方法在参数选择和应用方面积累了丰富的试验资料和工程经验。但这种方法存在较明显的缺点:不能考虑应力路径的影响;不能计算永久变形;在大应变时计算误差较大等。基于弹塑性本构模型的真非线性分析方法能够较好地模拟土体的实际响应特性，可以考虑堆石料在高围压下的剪切破坏，模拟复杂加载路径下堆石料的变形特性，并且能够直接计算出结构的残余变形，在理论上比等价线性方法更为合理［3－4］。

强震作用下防波堤表面混凝土面板、空心块、镇脚等护面结构与堆石体可能出现局部脱离、大滑移等接触非线性问题。以往研究中多是通过建立接触单元模拟接触面的力学行为，但此种方法局限于小变形假设，无法考虑沿接触面发生的大的滑移、分离运动。基于罚函数的动态主从接触法，将问题转化为无约束泛函极值问题，适用于分析接触面间可能发生的大滑移以及局部脱离问题［5－6］。

本文利用LS-DYNA有限元仿真软件，应用动力弹塑性分析方法，考虑结构间的动态接触作用，对核电站防波堤地震响应特性进行了仿真分析，重点讨论了防波堤挡浪墙地震位移响应、残余变形以及防波堤剪应变变化等情况，同时对地震动卓越频率、地震动峰值(PGA)以及护坡长度参数变化对防波堤地震响应的影响规律进行了分析，从而更好地对地震作用下防波堤的抗震安全进行准确评价，为合理的抗震措施提供理论基础。

1 仿真模型的建立

1.1 有限元模型及边界条件

本文的研究对象为某核电站厂区的斜坡式防波堤，该防波堤全长1 818.9 m，堤顶挡浪墙前后高程分别为11.4 m和12.7 m，根据工程的结构特点，将防波堤－地基体系的变形问题简化为平面应变问题对其进行数值仿真，其标准断面如图1所示。

图1 防波堤结构标准断面Fig.1 Cross-section of breakwater

为了克服模型边界对计算结果的影响，需要在土体四周施加人工边界来模拟无限区域。本文采用刘晶波等［7］提出的粘弹性人工边界对LS-DYNA自带的粘性边界进行改进，具体方法为在已经建好的有限元模型边界上沿法向延伸一层厚度相等、边界固定的边界单元，边界单元材料参数的计算公式为:

式中:E'为等效边界单元的输入弹性模量;E为与人工边界节点相连介质的弹性模量;rb为散射波源到人工边界的距离;h为等效实体单元沿人工边界法向的厚度。

本文在护面与堆石间建立动态主从接触模型。将护面结构作为主接触面，堆石体节点作为从接触面。每一时步进行护面体节点和堆石体接触面的穿透判断，无穿透不作处理，穿透则在护面体节点与堆石体接触面间引入法向接触力:

式中:m为接触穿透量;ni为接触点i处主接触面的外法线单位矢量;ki为接触刚度系数。在护面体节点上直接施加该法向接触力，而在堆石体接触面则施加其反力。接触面间的摩擦力采用库仑摩擦模型。

图2所示为防波堤有限元模型，模型共包括3 344个单元。模型底部和两侧为粘弹性人工边界，顶部为自由边界。

图2 防波堤有限元模型Fig.2 Finite element model of breakwater

1.2 计算材料参数

动力有限元弹塑性分析中常用的岩土材料屈服准则有Drucker-Parger(D－P)准则与Mohr-Coulomb(MC)准则。与M-C准则相比，本文在仿真中采用的D-P准则屈服面不存在尖角，避免了数值计算上的困难，其形式为:

式中:I1、J2分别为应力张量的第一不变量和偏应力张量的第二不变量;α、k为与岩土材料内摩擦角φ和粘聚力C有关的常数;不同的α和k在π平面上代表不同的圆，如图3所示。

图3 各屈服准则在π平面上的曲线Fig.3 The yield surface on the deviator plane

以往研究表明［8］，LS-DYNA中默认使用的DP1屈服准则在分析岩土工程问题时存在比较大的误差，而对于岩土工程平面应变问题，DP4屈服准则具有较高的计算精度，因此在应用中要对屈服准则做出修正。

设 c1、φ1为已知的初始粘聚力和内摩擦角，c2、φ2为采用平面应变DP4屈服准则的粘聚力和内摩擦角，即转换后的参数值。

对于关联流动法则下平面应变莫尔－库伦内切DP4准则，其α和k值可表示为:

对于莫尔－库伦外角外接圆DP1准则，其α和k值可表示为:

令 α1= α2，k1=k2，联立即可求得:

此时的c2、φ2即为采用平面 DP4屈服准则时在LS-DYNA中输入的粘聚力和内摩擦角。经过转后的土体计算参数如表1所示。挡浪墙及护面为混凝土材料，其材料参数如表2所示。

表1 岩土体材料参数Tab.1 Material parameters of rock and soil

表2 混凝土材料参数Tab.2 Material Parameters of concrete

1.3 模态及阻尼

在地震响应分析中，阻尼对地震响应起明显的衰减作用，但由于可以得到的阻尼资料比较有限，要想精确的确定结构阻尼是非常困难的，实际分析中，通常采用Rayleigh阻尼，其表达式为:

式中:α、β为待定系数;M为质量矩阵;K为刚度矩阵。利用结构振动的前两阶自然频率ω1、ω2及阻尼比ξ可以求得待定系数［9－10］:

地基土体的阻尼比通常在2% －5%;护面结构的阻尼比通常在2% －10%。以往研究表明［11］，弹塑性材料发生塑性流动时可以耗散大部分能量，因此在弹塑性动力分析中可以取相对小的阻尼比，本文取地基土体和护面结构的阻尼比分别为2%和5%，由模态分析可得 ω1=5.082、ω2=7.724，故可求得地基土体材料的 α、β 为0.122 和0.003;护面体材料的 α、β 为 0.305和 0.007 5。

1.4 初始地应力

天然土体中存在着由自重产生的初始应力，在进行动力分析之前首先要进行静力计算，以获得准确的初始地应力场，本文中静力计算分为三步:

(1)设置各种材料为弹性模型，设置干密度，使模型在重力作用下达到平衡。

(2)将堆石及地基土体材料设置为D－P模型，并赋予表1中所示的真实塑性参数，重新达到平衡。

(3)施加静水压力，模型重新达到平衡，得到初始地应力场。

2 计算结果分析

本文结合《核电厂抗震设计规范》［12］及《水运工程抗震设计规范》［13］中相关的规定，采用极限安全地震震动输入对防波堤结构进行地震响应仿真。规范中规定极限安全地震动水平加速度峰值不得低于0.15 g，竖直设计加速度峰值采用水平加速度峰值的2/3。本文对场地地震安全性评价报告给出的50年超越概率10%的基岩加速度时程曲线进行调整，图4所示为水平地震动加速度时程曲线，地震波峰值加速度为0.175 g，卓越频率为 0.7 Hz。

图4 水平向地震动激励Fig.4 The lateral acceleration time-history curve

防波堤的主要作用是保护核电厂区免受海潮波浪影响，挡浪墙高程是其发挥作用的关键设计参数，因此挡浪墙在地震中的位移响应及残余变形要予以重点关注。图5所示为挡浪墙前、后侧挡板顶部节点位移时程，从图中可以看出，竖直方向上的地震动对挡浪墙水平位移影响很小，两种地震输入下，挡浪墙前、后侧挡板顶点水平位移变化规律及幅值基本相同。双向地震动输入下，前、后侧挡板顶点水平位移响应峰值分别为0.321 m和0.3 m。与水平位移相比，前、后挡板的竖直位移幅值差别明显，两种输入工况下的竖向位移变化规律也不同。在相同输入工况下，后侧挡板的竖向位移大于前侧的，这是因为前侧堆石的滑移量大于外侧堆石的滑移量。双向输入工况下的挡板竖向位移明显大于单向输入工况下的结果，前、后侧挡板的竖直位移响应峰值分别为0.492 m和0.611 m。

图5 挡浪墙前、后侧挡板顶部节点位移时程Fig.5 Displacement time-history curves of wave barrier

表3给出了挡浪墙前、后侧挡板顶部节点的震后残余变形值。挡浪墙水平向残余变形对竖直向地震动输入不敏感，两种地震动输入下，前、后挡板水平残余变形量只分别相差5.39%和5.45%。竖直向残余变形受竖直方向上地震动影响较大，前、后挡板竖直向残余变形量分别相差11.31%和15.50%，因此对防波堤地震响应分析时，不能忽略竖向地震动输入的影响。双向地震动输入下，挡浪墙前、后侧挡板的竖直向残余变形量分别为－0.492 m和－0.611 m，挡浪墙高程的减小将影响防波堤的防浪性能。

表3 挡浪墙残余变形Tab.3 Permanent displacement of wave barrier

防波堤的破坏现象可宏观地描述为:在地震作用下，地基软弱土层首先发生较明显的沉降及滑移变形，并逐渐引起堤顶挡浪墙出现沉降，内、外侧护坡结构发生明显滑移，防波堤两侧坡脚处地面隆起。图6所示为双向地震输入地震结束时防波堤的水平位移，可以看出地震结束时，内、外海侧护坡分别向两侧产生了0.826 m和0.485 m的滑移。图7所示为防波堤竖直位移，可以看出受到护坡结构滑移及竖直向地震动的共同影响，防波堤顶部产生了明显的沉降，最大值为0.716 m，发生在挡浪墙后挡板外侧。

图6 防波堤水平位移云图Fig.6 The lateral displacement contour of breakwater

图7 防波堤竖直位移云图Fig.7 The vertical displacement contour of breakwater

上述防波堤破坏形式主要是由地基土体及堆石体内部积累了较大的塑性变形引起的，图8所示为双向输入地震结束时防波堤结构的塑性剪应变云图，由图可知地基软弱土层的塑性剪应变较大，最大剪应变达到17.28%。当防波堤内部的塑性剪切应变积累到一定程度时，结构将沿图中红色实线所示的潜在最危险滑移面发生失稳破坏。因此在防波堤选址及施工中要对地基软土层尤其是浅层软土予以重点关注。

图8 防波堤塑性剪应变云图Fig.8 The plastic shear strain contour of breakwater

3 影响因素分析

防波堤震损的主要模式为挡浪墙滑移及沉降，为了进一步研究防波堤在地震作用下的挡浪墙残余变形，分别对不同地震动卓越频率、地震动峰值及外侧护坡长度进行变参数分析。

3.1 地震动卓越频率的影响

为研究输入地震动卓越频率对挡浪墙残余变形的影响，本文从太平洋地震工程研究中心数据库中选取了五条不同卓越频率的真实地震波，其卓越频率分别为 0.25 Hz、0.70 Hz、1.46 Hz、2.41 Hz 和 3.10 Hz，并将其加速度峰值均调整为0.2 g，其傅里叶谱如图9所示。

图9 地震波傅里叶谱Fig.9 Fourier spectrum of seismic wave

图10所示为不同地震动卓越频率下挡浪墙前、后侧挡板的水平及竖直残余变形变化情况。由图可知，在输入地震动峰值相同的情况下，结构基频附近，输入地震动卓越频率对结构响应有显著的影响;当输入地震动卓越频率接近系统基频时，挡浪墙残余变形达到最大;输入地震动卓越频率远离系统基频时，挡浪墙残余变形减小，并且在相对高频区域，输入地震动卓越频率对残余变形的影响逐渐减小。

图10 地震动频率对残余变形的影响Fig.10 Variation of horizontal and vertical permanent displacements with remarkable frequence of ground motion

3.2 地震动峰值(PGA)的影响

本文选择卓越频率为0.7 Hz的地震波作为输入，其他参数不变，只改变输入地震动峰值，分别为0.1 g、0.15 g、0.2 g、0.25 g 和 0.3 g。

图11所示为不同地震动峰值输入下的挡浪墙残余变形。由图可知，挡浪墙水平及竖直方向上的残余变形随着输入地震动峰值的增大而增大，且竖直方向上的残余变形对地震动峰值的变化更为敏感。

3.3 护坡长度的影响

其他参数不变，只改变防波堤内海侧护坡长度，分别为7.5 m、12.5 m、17.5 m 和 22.5 m，挡浪墙残余变形如图12所示。从图中可以看出，随着护坡长度的增大，挡浪墙水平及竖直方向上的残余变形均有所减小，因此在设计过程中可以考虑适当增加防波堤外侧护坡的长度，以提高其抗震性能。

图11 地震峰值对残余变形的影响Fig.11 Variation of horizontal and vertical permanent displacements with peak value of ground motion

图12 护坡长度对残余变形的影响Fig.12 Variation of horizontal and vertical permanent displacements with length of slope protection

4 结论

本文运用LS-DYNA软件，采用动力弹塑性本构模型对某核电站防波堤进行了地震动力响应分析，所得结论如下:

(1)强震作用下，防波堤地基软弱土层发生较明显的沉降及滑移变形，并逐渐引起堤顶挡浪墙出现沉降，内、外侧护坡结构发生明显滑移，防波堤两侧坡脚处地面隆起。

(2)地震过程中地基软弱土层累积的塑性剪应变较大，是导致防波堤上部变形的主要原因，在防波堤选址及施工中要对地基软土层尤其是浅层软土予以重点关注。竖直方向地震动对防波堤挡浪墙位移响应及震后残余变形影响显著，在防波堤抗震评价及稳定性分析中不能忽略。

(3)地震动卓越频率对挡浪墙震后残余变形影响显著，建议在对防波堤进行抗震安全计算时要尽量选取卓越周期接近结构基频的地震波。挡浪墙震后残余变形随地震动输入峰值增大而增大，且竖直方向上的残余变形对地震动峰值的变化更为敏感。挡浪墙震后残余变形随护坡长度增大而减小，在设计过程中可以考虑适当增加防波堤内海侧护坡的长度，以提高其抗震性能。

［1］李炎保，蒋学炼，刘 任.防波堤损坏特点与其成因的关系［J］.海洋工程，2006，24(2):130－138.LI Yan-bao，JIANG Xue-lian，LIU Ren.Discussion on the relationship between characteristics and the reasons of breakwater failures［J］.The Ocean Engineering，2006，24(2):130－138.

［2］刘立平.防波堤的震损机理与抗震设防［J］.水运工程，1992，2:12 －15.LIU Li-ping.Failure mechanism and anti-seismic measures of breakwater［J］.Port and Water Engineering，1992，2:12－15.

［3］Kramer S L.Geotechnical earthquake engineering［M］.Upper Saddle River，New Jersey:Prentice Hall，1996.

［4］Wang X，Wang LB.Dynamic analysis of a water-soil-pore water coupling system ［J］.Comput Struct，2007，85(11 －14):1020－1031.

［5］Yuksel Y，Cetin K O，Ozguven O，et al.Seismic response of a rubble mound breakwater in Turkey［J］.Proceedings of the Institution of Civil Engineers(Martime Engineering，2004，157(4):151－161.

［6］Cihan K，Yuksel Y.Deformation of rubble-mound breakwaters under cyclic loads［J］.Coastal Engineering，2011，58:528 －539.

［7］刘晶波，李 彬.三维黏弹性静－动力统一人工边界［J］.中国科学(E 辑)，2005，35(9):966 －980.LIU Jing-bo，LI Bin.A unified viscous-spring artificial boundary for 3 － D static and dynamic applications［J］.Science in China(Series E)，2005，35(9):966 －980.

［8］邓楚键，何国杰，郑颖人.基于M－C准则的D－P系列准则在岩土工程中的应用研究［J］.岩土工程学报，2006，28(6):735－739.DENG Chu-jian，HE Guo-jie，ZHENG Ying-ren.Studies on Durcker-Prager yield criterions based on M-C yield criterion and application in geotechnical engineering［J］.Chinese Journal of Geotechnical Engineering，2006，28(6):735 －739.

［9］Dickenson S E，Yang D S.Seismically-induced deformations of caisson retaining walls in improved soils［J］.Proc Geotechnol EarthquakeEngSoilDyn，1998，75(2):1071－1082.

［10］Alyami M，Wilkinson S M，Rouainia M.Numerical analysis of deformation behavior of quay walls under earthquake loading［J］.Soil Dyn Earthquake Eng，2009，29:525 － 536.

［11］Arablouei A，Gharabaghi A R M，Ghalandarzadeh A，et al.Effects ofseawater-structure-soilinteraction on seismic performance of caisson-type quay wall［J］.Computers and Structures，2011，89:2439 －2459.

［12］GB50267－97，核电厂抗震设计规范［S］.中国国家地震局，1997.

［13］JTS 146－2012，水运工程抗震设计规范［S］.中华人民共和国交通运输部，2012.