基于CAS 法的核电工程异形水箱分布质量计算模型

李建波，杨 波，李志远，丁志新

(1.大连理工大学海岸及近海工程国家重点实验室，大连 116024；2.大连理工大学建设工程学部水利工程学院工程抗震研究所，大连 116024；3.中广核工程有限公司核电安全监控技术与装备国家重点实验室，深圳 518172)

核电工程中存在多种储液结构，如AP1000、华龙1 号顶部的非能动冷却水箱，具有形状复杂、储液量大等特点。在地震激励下，液体发生震荡与容器壁之间发生复杂的相互作用，直接影响核岛结构的动力特性和地震反应。因此液体和储液器壁之间的流-固耦合效应是核电结构抗震中需重点关注的内容。

现有核电结构的抗震分析，主要有两种手段：一是基于理论推导的等效力学模型；二是基于有限元的数值分析。核电结构的抗震分析中更关注的是地震激励下液体晃动对结构产生的相互作用，液体自身的复杂运动在抗震设计中通常可以忽略。因此，为了避免进行复杂的流体晃动计算，核电工程抗震规范中推荐使用Housner 附加质量模型。附加质量法是将液体等效为质量点，并以质量点运动产生的惯性力等效液体对器壁的动水压力，从而复杂的流-固耦合分析解耦为等效质量点的计算和质量点-结构的动力响应分析两个问题，极大地提升了计算效率。这方面的开创性工作由GRAHAM 和RODRIGUEZ[1]完成。他们首先基于线性势流理论首先得到了矩形水箱液体晃动等效模型的精确公式，但无穷级数表达式不便于工程应用。随后，HOUSNER[2-3]基于物理直观分析，对文献[1]中的级数公式进行了简化，得到了脉冲质量和晃动质量的近似简化公式。Housner 模型公式简单且能够较准确地模拟液体晃动中占主要作用的低阶模态，因此在土木、水利工程中得到了广泛使用。GOYAL 和CHOPRA[4]给出了圆形和矩形等不同截面进水塔结构动水压力附加质量的经验公式，为进水塔结构抗震安全设计的初步阶段提供了技术手段。FARID 和GENDELMAN[5]将圆柱形水箱中的液体等效为弹簧质量模型，研究了油罐临界区运动自由度与机械应力之间的数学关系，并确定了储罐的失效模式和关键位置。ZOU 和WANG[6]基于弹簧-质量模型提出了一种新的水平激励下流-固耦合晃动力学模型，并建立了晃动系统的耦合动力学方程，将计算结果与已有文献结果进行比较，验证了提出模型的合理性。刘焕忠等[7]基于石川岛播磨重工业株式会社对储液罐的一系列动力试验[8]结果，给出了储液罐附加质量的经验公式。与日本电气协会火力发电所的耐震设计指南[9]中建议的附加质量经验公式相比，文献[7, 9]的经验公式得到附加质量在分布和数值上存在很大差异。这表明，储液容器的抗震特性还缺乏一定的共识，现有的附加质量模型难以满足工程抗震要求。

以上附加质量模型大都是依据圆柱形容器和矩形容器假设进行的推导，直接应用到形状复杂的核电储液容器将对核电结构抗震设计带来较大的不确定性。为更准确地模拟液体的晃动特性，宝鑫等[10]在Housner 附加质量模型的基础上给出了适用于环形容器的附加质量经验公式，并从结构振动频率、变形和内力三个方面比较了该方法与工程界常用的 60% 和 100% 均布附加质量法的精度。李小军等[11]和雷墉[12]给出了底部为锥形的圆柱形PCS 水箱附加质量的经验公式，并基于ADINA 软件讨论了简化公式的模拟效果。但这些经验公式仍主要针对环形或锥形等单一形状，对于核电工程中更为一般的复杂形状储液容器的适用性仍存在不确定性。

常采用的数值方法包括：光滑粒子流动力学法(SPH)[13-18]、耦合欧拉-拉格朗日法(CEL)[19-21]和任意拉格朗日-欧拉算法(ALE)[22-24]等。宝鑫等[25]借助Fluid80 流体单元研究了土-结构相互作用对储液结构动力反应的影响。由于问题的复杂性，模拟液体与储液结构的动力相互作用需要精细的网格和极短的时间步长，这对计算机的存储和计算性能提出非常高的要求。近年来声学-结构耦合(CAS)方法被用于流-固耦合的分析。王丕光等[26]推导了椭圆柱体地震动水压力的解析解并利用CAS方法进行了验证。RAWAT 等[27]采用CAS 方法，研究了三维矩形水箱和圆柱形水箱在地震激励下的动力响应，推导了液体运动控制方程，并从晃动波高、动水压力分布等方面进行了综合评价。但是针对核电项目，安全壳和地基的有限元模型要百万自由度，直接采用CAS 法考虑流-固耦合效应，需要采用显示积分方法，步长一般为10-5s 量级，对于几十秒长的地震动，其计算量远超工程能接受的范围。因此，在保证能够满足工程需求的前提下，将流体和结构的动力分析进行解耦从而提高流-固耦合分析的计算效率，对于核电结构动力分析具有重要的工程实践价值。

为解决以上问题，本文基于CAS 方法提出一种流-固耦合分析的分布质量模型，兼顾计算效率和计算精度，适于核电工程中异形水箱的抗震分析。基于CAS 方法从动水压强的剥离出了脉冲动水压强，进而得到等效的分布质量模型，在后续分析中将以附加质量产生的惯性作用等效液体对结构的动水压力，实现流-固耦合分析。在上述两步计算中，计算量都要显著小于直接的流-固耦合分析。通过比较AP1000 安全壳的CAS 模型和附加质量模型的结构响应、计算时间，验证了提出模型的合理性和有效性。本文研究内容可为核电工程中流-固耦合抗震分析提供有益参考。

1 计算理论

本文重点关注CAS 方法计算得到的动水压强，从而进行流-固耦合分析的解耦，参考RAWAT[27]中有限元理论部分进行推导。

1.1 液体运动控制方程

假设储液容器的液体不可压缩、无粘、无旋，容器壁为刚性。对于理想液体，流体运动可以采用一个速度势函数φ(x,y,z,t)来表示，从而简化液体波动方程的求解。三维的液体晃动，其控制方程可以用拉普拉斯方程来表示：

式中， ∇2为三维拉普拉斯算子。液体晃动的加速度为ai=∂2φ/∂i2(i=x,y,z)。若不考虑质量力，用液体压力p表示液体运动的控制方程为：

式中：c为流体中的声速，c=；K为体积模量； ρ为密度。需要指出，当p为声压，c为声波时，式(2)即为声场波动方程。对于小幅晃动问题，仅考虑液体线性范围内的晃动，则液体动水压力p可表示为：

式中， ρL为液体的密度。方程的边界条件主要从2 个方面进行考虑，液体和器壁接触面B1和自由液面B2。

液体与储液容器接触的边界B1如图1 所示。在B1处的边界条件可以表示为：

式中：n为向外指向流体区域的法线；an(x,y,z,t)为指向液体区域沿外法线方向的加速度分量。液体自由液面为B2，假定流体相对于静止所时液面产生一个小幅晃动波高h，这个条件被称为线性化表面波条件，其中包含了重力波。则由伯努利方程可得z=HL处的动水压力为：

式中：g为重力加速度；HL为液体深度。自由液面小幅晃动h引起的动水压力可表示为p=ρLgh。将液体晃动简化为等效力学模型时，液体运动可被分解为同器壁运动一致的脉冲分量和自由振动的对流分量。当不考虑液体自由振动的对流分量时，在水平激励作用下液体与器壁共同运动，在自由液面不会产生竖向运动，则z=HL处的动水压力为：

由此可知，在ABAQUS 软件中通过设置不同的边界条件可将液体动水压力进行有效分解。取式(6)作为边界条件时，即设置自由液面压力为0，则计算得到脉冲分量产生的动水压强；若自由液面不设置边界条件式(6)，计算得到的则为包括脉冲分量和对流分量的总动水压强。

1.2 流-固耦合解耦分析

抗震设计中的核电结构一般包括地基、安全壳、PCS 水箱、液体四部分。本文提出的流-固耦合分析方法包括两个步骤：第一，针对PCS 水箱和液体，采用CAS 方法计算简谐荷载激励下的脉冲动水压力，从而得到分布附加质量；第二，将附加质量施加到地基、安全壳、PCS 水箱的模型中，进行抗震分析即可。计算流程图见图2。下文针对以上步骤进行详细说明。

图2 解耦思路示意图Fig.2 Schematic diagram of decoupling method

在步骤一中，可以先将PCS 水箱以外的单元删除，只保留PCS 相关单元，并建立三维声学单元，在水箱底部输入任意正弦荷载，得到容器各节点脉冲动水压强P(t)。假设单元节点控制的面积上脉冲动水压力均匀分布，则动水压力需要乘以该节点的控制面积，可表示为：

式中：Fp为节点处的动水压力；t为时刻；Anode为该节点的等效控制面积。如图3 所示，灰色阴影部分表示节点E的等效控制面积，可表示为：

图3 节点等效面积示意图Fig.3 Schematic diagram of node equivalent area

式中：Ai为与该节点关联单元的面积；di为关联单元节点的个数。

核电结构中，储液容器多为具有钢衬的混凝土结构，满足刚性假定，器壁节点的加速度与激励相同。与液体的脉冲分量比较，动水压力中的对流分量占比通常很小，并且对流分量的振动相位与器壁的运动相位不同，实际上可能起到减震的作用。因此，在计算等效质量时仅考虑脉冲动水压力。当液体运动只考虑脉冲分量时，液体和侧壁加速度保持一致。此时动水压强均由脉冲质量产生，则节点处的动水压力又可以表示为：

式中，a(t)为外部加速度激励。由于只考虑脉冲分量，Fp(t)与a(t)为同相位。为了减小计算误差，脉冲质量取二者的最大值进行相除得到：

基于CAS 方法计算的分布质量不同于简单的将水体总质量取一定的百分比进行平均分布[28]。由CAS 法得到的脉冲动水压力在离散意义上更加准确。计算等效质量时，仅增加了单元节点控制的面积上脉冲动水压力均匀分布的假定，没有引入更多不确定性，因此得到的等效质量在数值和分布上都具有更高的精度。

2 分布质量分析

2.1 CAS 计算模型验证

为验证基于CAS 方法的分布质量模型的精度和效率，采用提出的模型分析矩形储液容器在正弦激励下的动水压力时程、脉冲动水压力产生的作用力以及水体对器壁产生的总作用力，并与Housner 分布质量模型[2]的结果进行比较。

二维矩形水箱，长L和高h均为2 m，液体高HL为1.5 m，单元尺寸为0.1 m。底部加速度激励(t)=sin(4t)，激励持时为5 s，步长 Δt为0.01 s。提取距离自由液面0.1 m 和0.9 m 两个位置动水压强时程，如图4 所示。从图4 中可以看出，提出模型的计算结果与Housner 分布质量模型计算结果吻合性非常好。

图4 距自由液面不同深度处动水压力时程曲线Fig.4 Time-history curve of hydrodynamic pressure at different depths from free liquid level

为了进一步验证本文采用的基于CAS 方法计算的模型的正确性，还提取了容器底部脉冲动水压力产生的作用力以及水体对器壁产生的总作用力进行比较。计算结果如图5 所示，液体晃动产生的脉冲反力时程同外部激励相位保持一致；同时可以看出，Housner 分布质量模型计算结果与本文提出模型计算得到的结果几乎完全吻合，证明了本文方法的合理性和准确性。

图5 基底反力时程曲线比较Fig.5 Comparison of substrate reaction time history curve

2.2 激励对分布质量的影响

为了验证附加质量与激励的关系，取半径为3 m，高为6 m，水深为4.5 m 的圆柱形水箱作为分析对象，分别作用5 种不同的正弦激励。正弦激励如表1 所示，激励1～激励3 幅值相同，圆频率不同；激励3～激励5 幅值不同，圆频率相同。基于CAS 分布质量模型，计算各激励作用下的脉冲附加质量，沿高程分布如图6 所示。从计算结果可以看出， 5 种不同的激励作用下，沿容器高程的分布质量基本一致，表明分布质量是水箱中液体晃动的固有特征，仅与水箱形状、液体密度、液面高度等因素相关，与外荷载无关。因此基于CAS 计算分布质量时，施加形式简单且短持时的正弦荷载激励即可，从而达到简化计算，提高计算效率的效果。

图6 不同激励作用下附加质量沿高程分布图Fig.6 Elevation distribution of the additional mass under different excitations

3 AP1000 模型分析

3.1 AP1000 数值模型建立

3.1.1 模型建立

以第三代先进非能动反应堆AP1000 安全壳为例，进一步表明提出模型的良好精度和广泛适用性。建立了AP1000 安全壳整体模型，为了计算的效率和后处理的简便，核电厂模型仅考虑顶部重力水箱和安全壳两部分，通风口、挡板等构造物在建模时忽略。AP1000 安全壳几何特征[29]见图7，PCS 水箱内部直径Rin、外部直径Rout、高度Hw和厚度Ww分别为 10.668 m、27.13 m、11.8 m 和0.6 m，水箱底部倾斜角度为15°；屏蔽厂房直径Rs、高度Hs和厚度Ws分别为 44.2 m、71 m 和0.92 m。如图8 所示，AP1000 安全壳整体模型采用壳单元进行模拟，共计4766 个单元，最大单元尺寸为4.29 m，其中PCS 水箱模型共计2688 个单元，最大单元尺寸为1.33 m。

图7 AP1000 安全壳几何尺寸Fig.7 AP1000 containment geometry

图8 AP1000 安全壳有限元模型Fig.8 AP1000 containment finite element model

作为对照，同时建立流-固耦合整体模型，液体深度HL为9.8 m，采用声学单元进行模拟，共计6400 个单元，最大单元尺寸为1.33 m。

3.1.2 材料参数和结构阻尼

安全壳结构材料整体采用C40 混凝土建模，材料密度为2700 kg/m³，弹性模量为30 GPa，泊松比为0.3；液体部分采用水进行模拟，密度为1000 kg/m3，剪切模量为1.96 GPa。

在ABAQUS 中进行模态分析时是采用无阻尼系统，但在动力时程分析时，需要对结构设置一定阻尼。本文选取Rayleigh 阻尼，阻尼比为ξ=5%。

3.1.3 地震动的选择

美国ASCE4-16 规范与GB 50267-2019 规范在选取设计地震动方法上相同，取实测地震动或通过设计谱生成人工波，本文选用根据RG 1.60设计谱进行人工合成加速度时程，人工地震波持续时间为28 s，步长为0.01 s，加速度峰值为0.3g，如图9 所示。在计算中仅考虑X方向的地震波。

图9 RG1.60 人工地震波加速度时程Fig.9 RG1.60 artificial acceleration time history

3.2 PCS 水箱分布质量计算

按照图3 中给出的计算思路，完成对AP1000安全壳的抗震分析。首先计算得到PCS 水箱侧壁的质量分布，图10 是PCS 水箱的示意图。图11所示的是水箱内外侧壁各水箱沿X方向切割开后展开的分布质量包络云图。

图10 PCS 水箱示意图Fig.10 View of the PCS tank

图11 PCS 水箱侧壁分布质量分布云图 /kgFig.11 Mass distribution of the side wall of the PCS water tank

从图11 中可以看出，在X方向激励作用下，内外侧壁的分布质量关于Y方向表现较强的对称性，但并不是严格意义上的对称，大小上仍具有一定的区别。对于水箱外壁，分布质量在沿X(-X)方向取得最大且从侧壁顶部向底部逐渐增大，在底部取得最大；分布质量从X(-X)方向往Y(-Y)轴方向呈弧线逐渐减小，在Y(-Y)轴时取得最小值为0。对于水箱内壁，分布质量在沿X(-X)方向取得最大，且从侧壁顶部向底部逐渐减小，在底部取得最小；分布质量仍表现为从X(-X)方向往Y(-Y)轴方向呈弧线逐渐减小，在Y(-Y)轴时取得最小值为0。表明在单向激励作用下，内外侧壁分布质量沿激励方向向两侧逐渐减小，在垂直于激励方向取得最小值。

3.3 AP1000 安全壳动力计算

将3.2 节中计算得到的PCS 水箱侧壁各节点的分布质量值以附加质量点的形式施加到AP1000安全壳顶部的水箱侧壁对应节点进行动力分析，从而实现流-固耦合的解耦分析。分别从计算时间、安全壳响应等多个方面对本模型(后称附加质量模型)和直接采用CAS 方法进行流-固耦合分析(后称流-固耦合模型)计算的结果进行了比较。

3.3.1 计算效率

流-固耦合模型只能采用显式计算方法进行计算；而附加质量模型采用隐式计算方法进行计算。为了比较两种方法的计算效率，充分保证计算时间的可比性，以上的计算分析均在同一台计算机中进行。计算机的配置如下：内存为64 G，CPU 主频为3.6 GHz，6 核12 线程。安全壳网格划分相同，计算所需要的时间如表2 所示。

表2 两种方法计算时间比较Table 2 Comparison of the calculation time between the two methods

流-固耦合模型的计算时间要比附加质量模型高出2 个数量级。分析其主要原因有以下3 点：流-固耦合模型在计算时液体网格同时参与计算，导致计算量大大增加；在处理液体和器壁之间接触关系，判断边界条件上需花费大量时间；相关文献[30 - 31]指出，采用ABAQUS 软件进行显示动力计算时，考虑结构瑞利阻尼会对计算时间和计算稳定性造成影响。因此本文提出的附加质量模型在计算效率上远优于流-固耦合模型，大大提升了求解效率。

3.3.2 基底反力和倾覆力矩

表3 中给出了流-固耦合模型和附加质量模型在水平地震激励下，安全壳底部的基底反力和倾覆力矩的最大值以及计算误差。从表中可以看出，附加质量模型计算结果均大于流-固耦合模型计算结果，基底反力最值的误差约为10.7%，倾覆力矩最值的误差约为18.6%。证明本文基于CAS 的解耦的附加质量模型计算结果偏保守，对工程而言具有一定的安全裕度。

表3 AP1000 底部基底反力和倾覆力矩比较Table 3 Comparison of AP1000 base reaction and overturning moment

图12 比较的是流-固耦合模型基底反力时程和对流分量产生的基底反力时程。对流分量产生的基底反力通过总的基底反力时程减去脉冲分量产生的基底反力时程得到。从图中可以看出，对流分量产生的基底反力在流-固耦合及附加质量模型计算得到的基底反力中占比较小，对基底反力贡献不大。

图12 两种模型基底反力时程和对流分量基底反力时程比较Fig.12 Comparison of base reaction time histories and convective component base reaction time histories between the two models

表3 以及图12 的5 s～10 s 段，流-固耦合模型基底反力计算结果明显小于附加质量模型。其原因是，流-固耦合模型考虑了对流分量，对流分量一般与脉冲分量不同步，在一定程度上可以抵消脉冲分量对侧壁的冲击，相当于提供了一种流体阻尼，从而使结构的基底反力小于附加质量模型的计算结果。从而可以证明本方法的假设，即在安全壳结构抗震分析时，仅考虑脉冲分量作用而忽略对流分量即可满足工程设计需要。

3.3.3 安全壳不同高程加速度时程和反应谱

图13 给出了安全壳外壁不同高程处的加速度时程曲线。沿安全壳侧壁从顶部依次往下取了4 个观测点，编号为N1～N4，如图8 所示。从图13 可以看出，流-固耦合模型和附加质量模型计算得到的不同高程处加速度时程基本吻合。加速度随着高程的增加而增加，在安全壳顶部N1观测点处取得最大值。

图13 观测点N1～N4 加速度时程曲线Fig.13 Time-history curve of acceleration at observation points N1-N4

图14 给出了两种模型在观测点N1～N4的楼层加速度反应谱的比较曲线。从图14 中可以看出，在安全壳顶部观测点N1处反应谱峰值最大，且随着高程的降低呈现减小趋势；附加质量模型各观测点的反应谱值稍大于流-固耦合模型计算得到的结果。观测点N1～N4反应谱峰值误差分别为4.82%、10.09%、8.48%和6.51%，说明本文通过流-固耦合解耦得到的附加质量模型计算得到的楼层加速度反应谱结果偏保守。对工程设计而言，其结果是有利的。

图14 观测点N1～N4 楼层加速度反应谱曲线Fig.14 Response spectrum of floor acceleration at observation point N1-N4

3.3.4 安全壳不同高程加速度峰值

提取了安全壳沿激励方向所有高程位置节点的加速度最大值，提取点位置如图15 所示。不同高程处最大加速度和高程H之间的关系如图16所示，从图16 可以看出，两种方法计算得到的加速度最大值均随着高程的增加逐渐增加，在顶部时最大。附加质量模型计算得到的结果整体来说均大于流-固耦合模型计算所得的结果，仅在高程为38.7 m～47.2 m 的范围内略小于流-固耦合模型。因此总体而言，附加质量模型的计算结果偏保守。

图15 加速度最大值提取节点示意图Fig.15 Node path for extracting maximum acceleration

图16 加速度峰值随高程的变化Fig.16 Variation of the peak acceleration with elevation

3.3.5 安全壳应力包络图

图17(a)、图17(b)所示分别为流-固耦合模型和附加质量模型计算得到的安全壳最大Mises 应力包络分布云图。从图中可以看出，两种方法计算得到的包络云图分布基本一致。在地震激励作用下，上部PCS 水箱上应力值整体较小，在PCS 水箱和壳体交接处较大且分布较为对称；下部安全壳壳体应力值整体偏大，在高程约为60 m 的位置处较小，呈对称分布。

图17 RG1.60 地震波作用下安全壳最大应力包络云图Fig.17 Maximum Mises stress envelope diagram under RG1.60

两种方法最大Mises 应力值均出现在沿激励作用方向上安全壳侧壁转角位置，即图中引线标注部分。其中流-固耦合模型计算值为11.42 MPa，附加质量模型计算结果为10.74 MPa，误差约为5.95%。若采用本文提出的模型计算所得应力进行配筋设计，可以通过乘以裕度系数的方法来进行考虑，即可保证计算结果的保守性。

4 结论

本文提出的模型基于CAS 方法实现了流-固耦合问题的解耦分析，提高了核电工程储液容器的抗震分析效率。利用CAS 方法对水箱形状适应性强的特性，提高了基于Housner 简化模型的流-固耦合分析的精度。分别采用CAS 模型和本文提出的基于CAS 的分布质量模型对正弦激励下简单矩形水箱和地震激励作用AP1000 安全壳进行了动力分析，比较了两种方法的计算效率和计算精度等内容。主要结论如下：

(1) 采用本模型进行核电结构流-固耦合分析时，首先基于CAS 方法施加简谐荷载得到水箱附加质量，再将附加质量施加到核电结构整体模型中，进行抗震分析，大大提升了计算效率。

(2) CAS 模型适用于任意形状水箱的流-固耦合分析。通过设置边界条件可以计算得到各节点位置动水压强的脉冲分量，计算结果与推导的Housner 分布模型的解析解吻合较好。

(3) 针对AP1000 安全壳的抗震分析，本文提出的基于CAS 的分布质量模型计算效率明显高于直接进行流-固耦合分析的CAS 模型；计算得到的最大基底反力、倾覆力矩、不同高程楼层加速度反应谱均大于CAS 模型的计算结果，结果偏于保守。两种模型计算得到的最大应力包络云图基本一致。