核电厂结构竖向隔震理论及设计应用

韩世界，王长虹

(上海大学 土木工程系，上海 200444)

核电作为一种清洁、高效、优质的绿色能源，用以替代煤、石油等高污染性传统能源，为世界上各个国家所大力发展。由于核燃料的高放射性，核电厂房一旦经受地震灾害造成核泄露，则随之带来的生命伤亡、环境污染和经济损失将是难以估量[1]。

核电工程的结构抗震能力是保障核能安全应用的重要内容[2]。至今，地震灾害已引起多起严重的核电工程安全事故。如2011年日本东北部海域发生里氏9.0级大地震，引发巨大海啸，使得福岛核电厂多个机组发生停堆。强震使得核电厂外部电网中断，应急柴油发动机也因为海啸丧失功能，冷却功能失效，导致内部燃料过热熔毁，发生爆炸造成核泄漏。福岛核电厂辐射物质泄漏最终定性为7级核事故，这是迄今为止人类核能发展史上最严重的一次安全事故[3]。2011年美国东海岸发生5.8级地震，弗吉尼亚州震中附近的12个核电厂均有震感。在该次地震中，北安娜核电厂核电机组实现了停堆，核电厂丧失了场外电源，4台应急柴油机全部启动，厂房水平方向和竖向振动的加速度均超出规范要求。反应堆厂房内部的墙体出现了裂缝，放射性废物储存罐也发生了滑移[4]。

隔震结构是世界上广泛使用的控制地震响应的技术。目前，核电厂结构隔震研究的热点为三维隔震理论。Fujita等[5]提出一种三维隔震系统构想，采用碟形弹簧作为竖向隔震器，橡胶支座作为水平隔震器。Somaki等[6]开发了三维隔震系统，由叠层橡胶支座用于水平方向隔震，碟形弹簧用于竖向隔震，并进行了足尺试验，验证其力学性能。Whittaker等[7]开发了一种计算核电厂隔震性能的程序，嵌入了开源代码OpenSees。Najafijozani等[8]探讨了多种不同的核电厂结构竖向隔震系统。

20世纪80年代初期，我国开始进行建筑结构水平隔震技术的研究，如摩擦滑移支座隔震技术和叠层铅芯橡胶支座隔震技术。熊世树[9]在组合碟形弹簧层间填充粘弹性材料，增大阻尼性能，并和铅芯橡胶支座组合成为三维隔震支座。陈兆涛等[10]采用水平隔震铅芯橡胶支座和竖向隔震液压油缸组合成为新的三维隔震支座。

我国属于多地震国家，基本烈度大于7度的面积约占国土面积的1/3，大于或等于6度的面积达到国土的60%。另外，我国目前已经进入核电工程快速发展时期，在高地震区或一般地震区发展核电工程将成为可能[11]。王涛等[12-13]开始核电厂结构的三维隔震技术的理论研究，并进行了几何缩尺比例为1/15的振动台试验。魏陆顺等[14]进行了核电厂结构三维隔震研究，提出了一种抗摇摆装置。

竖向振动和摇摆效应控制是核电厂结构三维隔震理论的核心问题，但是目前缺乏对两者解耦问题的研究。本文首先将介绍我国自主创新设计的CAP1400型核电厂[15]结构和三维隔震支座参数。建立核电厂结构竖向振动和摇摆效应长(短)轴方向的双自由度梁-弹簧模型，讨论竖向隔震支座刚度分布对竖向振动和摇摆效应的影响。并采用大型有限元软件ANSYS对核电厂结构进行三维数值模拟验证。分析结果以满足预定的竖向振动加速度、位移和摇摆效应的限值为目标，确定了最优竖向隔震支座数量和分布位置。

1 设计方法与参数

如图1所示，采用半逆的方法计算最优总刚度，即先假设一个竖向总刚度，通过动力学理论计算得出最优支座位置，最终校准支座布置数量和总刚度。沿CAP1400型核电厂长(短)轴方向建立双自由度的梁-弹簧模型。通过刚度凝聚的方法，集中为3个等刚度支座群，逐步改变中间支座群的位置，寻找控制竖向振动和摇摆效应的最优支座布置。

图1 半逆法的设计流程

1.1 CAP1400型核电厂参数

如图2所示，CAP1400型核电厂由一座核岛和两座裙房组成。核电厂结构全长为92.10 m，宽为60.47 m，高为87.75 m，占地面积为461.59 m2。质心高为19.13 m，距离图2(d)中左边缘和下边缘分别为45.84 m和27.03 m。

(a)三维图

CAP1400型核电厂是在AP1000型核电厂的基础上完成了引进、消化、吸收和再创新的自主化历程，总重为20.6万吨，上部结构可视为刚体。在夏祖讽[16]提出的AP1000型核电厂的抗震设计参数基础上，竖向隔震设计目标如表1所示。

表1 核电厂结构竖向隔震设计目标

1.2 隔震支座参数

如图3所示，根据刘文光等[17]发明的倾斜旋转型三维隔震支座进行隔震设计。装置由上部水平铅芯橡胶支座和下部的倾斜旋转的铅芯橡胶支座组成，中部通过转动钢板连接。

(a)

倾斜旋转型三维隔震支座具有良好的竖向隔震刚度和阻尼耗能特性，竖向隔震参数如表2所示。

表2 三维隔震支座的竖向力学参数

图4 三维隔震支座竖向力-位移关系图

2 理论模型

将CAP1400型核电厂视为一个带有转动惯量的单质点，通过刚性结构与底板梁联结，在竖向地震作用下产生上下振动和摇摆效应。简化竖向隔震支座为弹簧，底板梁通过弹簧支承于刚性地基上。

2.1 动力学模型

动力学模型如图5所示。其中：L为核电厂的长(短)轴长度；M为核电厂的集中质量；J为核电厂的转动惯量；H为核电厂重心处距底板梁的高度；x0为核电厂重心距底板梁中心水平距离，其中长轴方向为0.213 m，短轴方向为5.031 m；Ki(i=1,2,3)为三个竖向隔震支座群的凝聚刚度;x1为中间支座位置距底板梁中心的水平距离；C为阻尼系数；u为平动(竖向)自由度；θ为摇摆(转动)自由度。

图5 CAP1400型核电厂的梁-弹簧模型

根据达朗贝尔原理，核电厂结构单质点双自由度模型的振动方程为：

(1)

(2)

2.2 动力学参数

梁-弹簧模型绕梁中点转动，假设竖向位移向上为正，摇摆效应逆时钟转动为正。

当质点有单位竖向振动(u=1)或摇摆效应时(θ=1)，杆端的弹簧力如图6所示。

如图6(a)所示，通过竖向静力平衡条件得到：

(3)

如图6(b)所示，通过弯矩平衡条件得到：

(a)梁-弹簧模型单位竖向振动

(4)

通过式(3)和式(4)得到刚度矩阵为：

(5)

同理，当梁-弹簧有单位竖向振动或摇摆效应加速度时，由静力平衡条件得到的质量矩阵为：

(6)

根据质量矩阵推算的阻尼矩阵为：

[C]=2πζ[M]

(7)

式中:ζ为梁-弹簧模型的阻尼系数，即核电厂隔震系统等效阻尼比，可取ζ=0.2[18]。

2.3 计算方法

将位移{U}按振型分解为:

{U}=[A]{q}

(8)

将式(8)代入式(1)并利用正交条件,可得到2个独立的广义坐标方程：

(9)

式中:ωi和ξi分别为i阶振型的圆频率和阻尼比；i=1,2；γi为振型参与系数，可表示为：

(10)

3 理论分析结果

根据核电厂场地的多样性与复杂性，选择4条不同周期的天然波：Iwate波、Lomap波、Sansimeo波和El-Centro波。根据美国核监管委员会颁布的监管指导RG1.6反应谱标准，拟合得到2条人工波。地震波加速度反应谱如图7所示，周期覆盖了0.01～10 s的范围。

图7 地震波加速度反应谱

根据CAP1400型核电厂结构和隔震支座力学参数，采用MATLAB编制计算程序，地震波激励振幅统一取为0.6 g，计算得到梁-弹簧模型的竖向振动、摇摆效应与隔震支座总刚度K的关系曲线。

3.1 竖向总刚度

采用半逆法结合核电厂的竖向位移和加速度响应计算最优的总刚度。

图8(a)给出了在不同的竖向刚度下，梁-弹簧模型竖向加速度的响应曲线。图8(b)为图8(a)中方形线框的放大区域。该区域红色虚线范围内的刚度值不但满足表1的设计要求，并且支座布置数量较少。竖向总刚度取值为1.0×104～1.4×104kN/mm。

(a)梁-弹簧模型竖向加速度与刚度关系曲线

图9(a)给出了在不同的竖向总刚度下，梁-弹簧模型摇摆加速度响应的关系曲线。图9(b)为方形线框内放大区域，弹簧刚度选取范围为1.0×104～1.4×104kN/mm仍能满足表1的设计要求。

(a)梁-弹簧模型摇摆加速度与刚度的关系曲线

图10给出了在不同的竖向总刚度下，梁-弹簧模型竖向位移响应的关系曲线。综合考虑图7的取值范围和表1的设计要求，梁-弹簧模型的竖向总刚度仍然选取为：1.0×104～1.4×104kN/mm。

图10 梁-弹簧模型竖向位移与刚度的关系曲线

图11给出了在不同的竖向总刚度下，梁-弹簧模型摇摆角度响应的关系曲线。

图11 梁-弹簧模型摇摆角度与刚度关系曲线

通过以上关系曲线的分析，梁-弹簧模型的竖向总刚度取值范围为1.0×104～1.4×104kN/mm，将满足表1设计要求的限值。在随后的计算中，选取刚度平均值1.231×104kN/mm作为设计值。

3.2 沿长轴的支座合理布置

图12给出了在中间支座群不同的布置位置上，梁-弹簧模型竖向加速度响应的关系曲线。图中水平短直线为表1设计要求的限值。

图12 梁-弹簧模型最大竖向加速度

图13给出了在中间支座群不同的位置上，梁-弹簧模型摇摆加速度响应的关系曲线。通过分析发现当中间支座群布置在底板梁中心位置时，摇摆效应最小。

图13 梁-弹簧模型最大摇摆加速度

图14给出了在中间支座群不同的位置上，梁-弹簧模型竖向位移响应的关系曲线。图中水平短直线为表1设计要求的限值。

图14 梁-弹簧模型最大竖向位移

图15给出了在中间支座群不同的布置位置上，梁-弹簧模型摇摆角度响应的关系曲线。通过分析可以发现当中间支座群布置在底板梁中心位置时，摇摆效应最小。

图15 梁-弹簧模型最大摇摆角度

3.3 沿短轴的支座合理布置

同理，CAP1400型核电厂隔振体系沿短轴方向也可分成3个支座群，通过计算得到的竖向、摇摆的加速度和位移响应相对于长轴的动力响应小。

图16为中间支座布置位置与模型摇摆加速度响应的关系曲线，结合表1控制核电厂结构摇摆效应，x1选取竖直短直线范围内的距离时摇摆效应最小，值域为0～5.031 m。

图16 梁-弹簧模型短边方向最大摇摆加速度

3.4 隔震支座布置

结合梁-弹簧模型长、短边支座合理布置位置与核CAP1400型核电厂工程参数，给出了如图17所示的隔震支座平面布置图。

图17 隔震支座平面布置图

平面布置图共采用361个三维旋转型隔震支座，分为4个支座群，隔震支座的竖向力学参数如表2所示。图中实线交点为支座布置点，间距为3 m。沿长轴方向，左、右两端支座群均布置7×7+8×8个隔震支座。核岛部分支座群采用环向布置115个隔震支座，裙房部分支座群布置4×9个隔震支座。长、短轴虚线交点为核岛重心所在位置，4个支座群的刚度中心布置在长、短轴虚线交点的下方2 m处，如虚线圆圈所示。

4 数值模型验证

采用大型有限元软件ANSYS建立CAP1400型核电厂有限元模型。上部核电厂结构采用实体单元，隔震支座采用弹簧单元，通过硬点布设支座位置。

图18 CAP1400型核电厂三维数值分析模型图

根据CAP1400型核电厂和隔震支座工程参数，在6条振幅为0.6 g的地震波激励下，计算得到核电厂结构重心处理论与数值分析时程曲线。

以El-Centro波激励为例，核电厂结构的竖向位移时程曲线与数值分析得到的时程曲线如图19所示。理论解与数值解的频率高度相似；理论的竖向位移幅值为0.018 m，数值分析得到竖向位移幅值为0.019 m，相对误差为8.3%，且均满足表1的设计要求。

图19 核电厂结构竖向位移的理论与数值分析时程曲线

由于文章的篇幅限制，将核电厂在6条地震波激励下的竖向位移、摇摆角度、竖向加速度和摇摆加速度的响应列入表3。计算得到的数值解和理论解的频率相似。对比其幅值响应，在6条地震波的激励下，通过支座的合理布置，核电厂的竖向响应均满足隔震设计要求。

表3 不同地震波激励下核电厂的竖向响应

5 结 论

为了控制核电厂在地震作用下的竖向振动和摇摆效应，我国第三代CAP1400型核电厂拟采用竖向隔震技术。通过竖向隔震的研究，包括隔震方案的设计、计算模型的选择以及隔震结构动力响应分析，得到以下三点结论。

(1)通过建立双自由度的梁-弹簧模型，进行解耦分析，得到了核电厂结构的竖向振动、摇摆效应与隔震支座刚度的关系曲线。

(2)通过有限元计算结果和理论计算结果对比，两者接近，验证了梁-弹簧理论模型的有效性。

(3)通过研究核电厂的竖向抗震要求，提出了一种控制支座布置和数量，降低核电厂摇摆效应的半逆法设计理论。

附录

%Author:Hanshijie

%E-mail:haoshijie@shu.edu.cn

%Date:2020/01/02

clear

clc

s=1 200; %Cycle index

kTotal=zeros(1,s); %Total stiffness of nuclear power plant

pVAcc=zeros(4,s);pRAcc=zeros(4,s);%The acceleration response of the model under different stiffness pVAcc:Vertical acceleration

pVDis=zeros(4,s);pRDis=zeros(4,s);%The displacemen response of the model under different stiffness pRDis:Rotational displacement

x0=0.213;x1=0.1;h=20;l=92;

for j=1:s

if j<=1 000

m=206;J=206*1 010.7;k1=6.77+6.77*2*(j);k2=6.77+6.77*2*(j);k3=6.77+6.77*2*(j);

else

m=206;J=206*1 010.7;k1=2*6.77*10^3+6.77*200*(j-1 000);k2=2*6.77*10^3+6.77*200*(j-1 000);k3=2*6.77*10^3+6.77*200*(j-1 000);

end

kTotal(j)=k1*3;

m11=m;m12=m*x0;m21=m12;m22=m*x0^2+J+m*h^2;

m=[m11,m12;m21,m22];

k11=k1+k2+k3;k12=(k2-k1+k3*2*x1/l)*l/2;k21=k12;k22=(k1+k2+k3*4*x1^2/l^2)*l^2/4;

k=[k11,k12;k21,k22];

%Modal analysis method

cn=2; %Model degree of freedom

[x,d]=eig(k,m);

d=diag(sqrt(d));%Solve for the circular frequency of the structure

% The circular frequencies of the structure are arranged in the order from small to large, and the order of their size is recorded

[d,indexf]=sort(d);

T=2*pi./d;

x=x(:,indexf);

f=1./T;

%Find the vibration mode participation coefficient

zhcan=zeros(cn,1);

for i=1:cn

x(:,i)=x(:,i)/x(cn,i);

zhcan(i)=x(:,i)'*m*ones(cn,1)/(x(:,i)'*m*x(:,i));

end

PCoeff=zhcan';

%Loading seismic wave

Lo1=load('Iwate.txt');

Lo1=Lo1/max(abs(Lo1));

n1=length(Lo1);

Lo2=load('LOMAP.txt');

Lo2=Lo2/max(abs(Lo2));

n2=length(Lo2);

Lo3=load('Sansimeo.txt');

Lo3=Lo3/max(abs(Lo3));

n3=length(Lo3);

Lo4=load('ELCENTRO.txt');

Lo4=Lo4/max(abs(Lo4));

n4=length(Lo4);

Lo5=load('New1x.txt');

Lo5=Lo5/max(abs(Lo5));

n5=length(Lo5);

Lo6=load('New2x.txt');

Lo6=Lo6/max(abs(Lo6));

%The dynamic solution of the single-degree of freedom model

%The newmakebate is solving function of structural dynamics equation

M1=x(:,1)'*m*x(:,1);K1=x(:,1)'*k*x(:,1);

M2=x(:,2)'*m*x(:,2);K2=x(:,2)'*k*x(:,2);

c1=2*0.2*M1*d(1);c2=2*0.2*M2*d(2);

m0=x(:,1)'*m*[0.6*9.8;0];

m00=x(:,2)'*m*[0.6*9.8;0];

[u_1,v_1,aa_1]=newmakebate(Lo1,n1,m0,M1,K1,c1);

[u2_1,v2_1,aa2_1]=newmakebate(Lo1,n1,m00,M2,K2,c2);

[u_2,v_2,aa_2]=newmakebate(Lo2,n2,m0,M1,K1,c1);

[u2_2,v2_2,aa2_2]=newmakebate(Lo2,n2,m00,M2,K2,c2);

[u_3,v_3,aa_3]=newmakebate(Lo3,n3,m0,M1,K1,c1);

[u2_3,v2_3,aa2_3]=newmakebate(Lo3,n3,m00,M2,K2,c2);

[u_4,v_4,aa_4]=newmakebate(Lo4,n4,m0,M1,K1,c1);

[u2_4,v2_4,aa2_4]=newmakebate(Lo4,n4,m00,M2,K2,c2);

[u_5,v_5,aa_5]=newmakebate(lo5,n5,m0,M1,K1,c1);

[u2_5,v2_5,aa2_5]=newmakebate(lo5,n5,m00,M2,K2,c2);

[u_6,v_6,aa_6]=newmakebate(lo6,n6,m0,M1,K1,c1);

[u2_6,v2_6,aa2_6]=newmakebate(lo6,n6,m00,M2,K2,c2);

%Synthesis of acceleration

Acc1=zeros(2,n1);

Acc1(1,:)=PCoeff(1)*aa_1;Acc1(2,:)=PCoeff(2)*aa2_1;

Acc1=x*Acc1;

pVAcc(1,j)=max(Acc1(1,:));

pRAcc(1,j)=max(Acc1(2,:));

Acc2=zeros(2,n2);

Acc2(1,:)=PCoeff(1)*aa_2;Acc2(2,:)=PCoeff(2)*aa2_2;

Acc2=x*Acc2;

pVAcc(2,j)=max(Acc2(1,:));

pRAcc(2,j)=max(Acc2(2,:));

Acc3=zeros(2,n3);

Acc3(1,:)=PCoeff(1)*aa_3;Acc3(2,:)=PCoeff(2)*aa2_3;

Acc3=x*Acc3;

pVAcc(3,j)=max(Acc3(1,:));

pRAcc(3,j)=max(Acc3(2,:));

Acc4=zeros(2,n4);

Acc4(1,:)=PCoeff(1)*aa_4;Acc4(2,:)=PCoeff(2)*aa2_4;

Acc4=x*Acc4;

pVAcc(4,j)=max(Acc4(1,:));

pRAcc(4,j)=max(Acc4(2,:));

Acc5=zeros(2,n5);

Acc5(1,:)=y(1)*aa_5;Acc5(2,:)=y(2)*aa2_5;

Acc5=x*Acc5;

pVAcc(5,j)=max(Acc5(1,:));

pRAcc(5,j)=max(Acc5(2,:));

Acc6=zeros(2,n6);

Acc6(1,:)=y(1)*aa_6;Acc6(2,:)=y(2)*aa2_6;

Acc6=x*Acc6;

pVAcc(6,j)=max(Acc6(1,:));

pRAcc(6,j)=max(Acc6(2,:));

%Synthesis of displacement

Dis1=zeros(2,n1);

Dis1(1,:)=y(1)*u_1;Dis1(2,:)=y(2)*u2_1;

Dis1=x*Dis1;

pVDis(1,j)=max(Dis1(1,:));

pRDis(1,j)=max(Dis1(2,:));

Dis2=zeros(2,n2);

Dis2(1,:)=y(1)*u_2;Dis2(2,:)=y(2)*u2_2;

Dis2=x*Dis2;

pVDis(2,j)=max(Dis2(1,:));

pRDis(2,j)=max(Dis2(2,:));

Dis3=zeros(2,n3);

Dis3(1,:)=y(1)*u_3;Dis3(2,:)=y(2)*u2_3;

Dis3=x*Dis3;

pVDis(3,j)=max(Dis3(1,:));

pRDis(3,j)=max(Dis3(2,:));

Dis4=zeros(2,n4);

Dis4(1,:)=y(1)*u_4;Dis4(2,:)=y(2)*u2_4;

Dis4=x*Dis4;

pVDis(4,j)=max(Dis4(1,:));

pRDis(4,j)=max(Dis4(2,:));

Dis5=zeros(2,n5);

Dis5(1,:)=y(1)*u_5;Dis5(2,:)=y(2)*u2_5;

Dis5=x*Dis5;

pVDis(5,j)=max(Dis5(1,:));

pRDis(5,j)=max(Dis5(2,:));

Dis6=zeros(2,n6);

Dis6(1,:)=y(1)*u_6;Dis6(2,:)=y(2)*u2_6;

Dis6=x*Dis6;

pVDis(6,j)=max(Dis6(1,:));

pRDis(6,j)=max(Dis6(2,:));

end

以上均为matlab程序，请复制在matlab中运行；

其他程序请访问：https://pan.baidu.com/s/1P3ijj55_MytTIpzWouz3eg

提取码：xj36