基于ANSYS平台的非岩性地基条件下核岛厂房结构三维地震响应分析

尹训强 赵文燕 王桂萱

（大连大学土木工程技术研究与开发中心，大连116622）

0 引 言

核电是高效型、安全型能源，然而由于国内不同地区地质条件差异较大，以及近些年来我国核电的飞速发展导致核电厂址选址问题日益突出，因此，内陆非岩性厂址已成为目前核电站建设的必然选择［1-2］。其中，无限地基和近场地基非线性特征对核电结构地震响应的影响成为核电厂选址的关键技术问题。

考虑土-结构相互作用效应（SSI）的方法主要有比例边界有限元法（SBFEM）［3］、阻尼抽取法（DSEM）［4］等子结构法以及人工边界法［5］等直接法。与其他方法相比，人工边界法通过在半无限地基的外边界设置反射与折射等能量传播方式，能够有效吸收反射波动能量，同时可方便考虑桩土相互作用和近场地基非线性的影响，其中，由于黏性边界［6］易于操作，概念清晰，在工程中应用较广泛，同时也是核电抗震规范推荐方法。

近年来，非岩性地基［7］场地条件下核岛厂房地震响应分析也备受国内外学者关注，并取得了很多有价值的成果。王桂萱，盛超等［8］建立了考虑近场非线性和土-结构相互作用的核岛厂房三维地震响应分析模型，研究非岩性地基对结构响应的影响。邹德高，隋翊等［9］在考虑地基土的非线性动力特性，采用人工边界和等效结点荷载的地震动输入方法的基础上，考虑核岛结构与深厚覆盖土层地基的相互作用，研究结构动力响应的变化规律。尹训强，薛志强等［10］采用等效线性和黏性人工边界模拟近场地基特征，对土质地基条件下的核岛厂房结构进行了拟三维地震分析；然而，拟三维分析需在两垂直平面内简化结构模型，也不能实现三方向上地震波的输入，不能真实反映土体运动对变形特性参数的影响。

针对以上问题，本文基于等价线性法和黏性边界的相关理论，开发了三维等价线性单元和黏性人工边界单元，并结合UPFs二次开发特点嵌入至通用有限元软件ANSYS中，建立了非岩性地基-核岛厂房结构的动力相互作用分析模型，进而通过数值算例及实际工程对该计算模型进行合理性及实用性的验证分析，并探讨无限地基与近场非线性特征对核电工程结构地震响应的影响。

1 非岩性地基条件下核岛厂房结构地震响应分析模型

如图1所示，核岛厂房结构-非岩性地基相互作用分析模型由近场地基与模拟无限地基的黏性人工边界模型所构成。如前所述，近场非岩性地基的非线性动力特性和不均匀性，以及无限地基辐射阻尼是需要考虑的影响因素。

图1 非岩性地基条件下的相互作用系统分析模型Fig.1 Structure-foundation interaction model of nuclear island factory building under the condition of non-lithology foundation

1.1 近场地基非线性模拟

在强震作用下，非岩性地基条件下土体的非线性动力特征显著，直接非线性时域方法由于本构模型的复杂性以及数值计算的稳定性在核电工程中并未得到广泛应用。美国土木工程师协会（ASCE）建议的等价线性法［11］概念明确且优势突出，是不可替代的主流方法。

等价线性法［12-13］通过线性迭代求解非线性问题，关键思想是若动力分析迭代前后单元应变变化显著，则利用每个单元的等效剪应变γeff=0.65γmax，参照所赋予材料的D-γ与G-γ关系曲线重新得出修正剪切模量G和阻尼比D，进行下一次动力分析，直到应变变化满足收敛要求。具体计算过程如图2所示。

图2 等价线性解析流程图Fig.2 Equivalent linear flow chart

1.2 黏性人工边界模型及地震动输入

由于近场地基由等价线性法来描述，每次迭代的动力分析中，土体则认为是线性材料。因此，无限地基辐射阻尼的影响可由黏性人工边界来模拟。

黏性人工边界模型［14］的核心是在地基的边界处布置阻尼器吸收边界处的反射波波动能量（图1），并通过施加节点等效荷载模拟截断外边界处实际波场的应力边界条件。在三维计算模型中，输入地震动转化为等效荷载的具体计算公式为

式中：F为施加在人工边界节点处的等效荷载；下标n表示当前边界面的外法线方向；下标tj表示对应的切线方向；Ct与Cn分别为黏性人工边界节点处的切向阻尼系数和法向阻尼系数；u̇f为自由波场条件下边界处节点的速度；σf为自由条件下边界处的正应力；τf为相应的切应力。

对于非均质近场地基的地震动输入，可先进行自由场分析得到外边界处的节点响应，而后利用式（1）进行各节点处的等效荷载计算，在具体数值实现中，可建立土柱有限元模型进行自由场分析［15］。

2 基于ANSYS平台的地震响应分析模型开发

2.1 三维等价线性单元

传统的有限元软件中，数值迭代往往是在某一荷载步中完成，无法实现线性动力分析的循环迭代。为了在ANSYS中基于等价线性法描述非岩性地基的动力特性，必须使地基单元能够描述D-γ与G-γ关系曲线且能实现线性迭代过程。本文提出一种地基侧三维等价线性单元，该单元形式是对已有等参元材料模型的修改。

在用户可编程特性（UPFs）中，用户可修改的UserElem.F是创建用户单元的接口子程序。通过该接口程序，可从标准ANSYS程序中获取原地基单元的源数据，包括节点号、节点坐标、弹性模量、密度、积分点等，该子程序接口中可同时创建多个单元，不同的单元类型由KEYOPT关键选项实现转换。等价线性单元与标准单元的不同关键在于每次线性动力迭代分析后的刚度阵和阻尼阵的变化，该处需要通过UPFs提供的自定义函数来实现阻尼比和剪切模量更新后实时数据的传递，之后便可计算单元的刚度阵以及阻尼阵等输出数据返回给标准程序，保证下一次动力分析时为修正的刚度阵及阻尼阵。该用户单元子程序的编制流程如图3所示，程序的正文部分与标准等参元的编制大致相同，仅为剪切模量G和阻尼比D的差异及误差分析的添加，需要注意的是，第2步中需要接口子程序User01.f实现实时数据的更新。

图3 等价线性用户单元子程序的编制流程Fig.3 The compiling sketch of interface subroutine

2.2 黏性人工边界单元

以往的文献中，人工边界的开发往往是将物理元件参数乘以每个节点的控制面积，然后在每个节点施加相应的单元类型，这样的开发形式导致前处理及控制面积的计算均较为繁琐。本文提出在原地基网格的基础上创建一种三维黏性人工边界面单元，该单元形式是对外已有边界节点的修正，网格划分与原始网格一致，易于操作。

如图4所示为所开发单元的创建过程，通过节点1-4所围成的区域实现人工边界的施加。首先，计算该单元的面积并得到每个节点的控制区域；然后，可由物理元件系数得到边界节点在该单元区域内对阻尼阵的贡献并将其团聚到相应节点，从而可以得到阻尼对角矩阵，并将其添加到整体的计算矩阵之中。依照图4所示的坐标系，边界单元阻尼矩阵的表达式如下：

图4 三维黏性边界单元示意图Fig.4 Sketch of three-dimensional viscoelastic boundary element

式中：CBN为所施加阻尼器的切向系数；CBT为法向系数；Ai（i=4）为边界节点控制区域的面积。

另外，在单元子程序编制中，可方便获取节点的坐标以及每一荷载步的响应（位移、速度以及加速度），从而可依据式（1）实现节点等效荷载的施加。

2.3 时域计算模型的建立

非岩性地基条件下的核电结构地震响应时域模型在ANSYS的建立，三维等价线性单元及黏性人工边界单元用户子程序的编制与连接是关键技术。子程序的编制是对标准ANSYS的功能添加，需编译连接成功后才能实现新的单元类型模块的嵌入，进而结合成熟的有限元软件前后处理及求解器，实现本文模型的建立，详细操作过程如表1所示。

表1 所开发时域模型的建立流程详解Table 1 Flowchart of the establishment of proposed model

3 数值验证

本节以某拟建内陆核电场址的地基资料为基础，进行自由场动力分析，并分别采用所开发的计算模型与商用软件SuperFLUSH/2D进行对比，验证本文提出模型的合理性及精确性。

3.1 材料参数

依据地质勘察资料，该厂址地基水平层状分布较为显著，为简化计算，假定为水平成层地基。各分层的材料属性如表2所示。

表2 地基计算参数Table 2 Calculation parameters of soil

图5给出由试验得到的各土层对应的G-γ和D-γ曲线。

图5 各土层对应的G-γ和D-γ曲线Fig.5 G-γand D-γcurves corresponding to each layer

3.2 有限元模型

基于ANSYS的三维土柱自由场有限元模型如图6（a）所示，底部为半无限地基，四周为自由边界，选取截断范围为4.0 m×4.0 m×47.7 m，各土层的深度由表1所示每层土厚度可知。基于SuperFLUSH/2D的二维有限元模型如图6（b）所示，利用该模型计算时，两个水平方向的响应通过二次动力计算获得。在计算分析中，在基岩处的网格最大高度为3.0 m，在软弱土层局部网格最大高度为1.0 m，满足影响频率范围内简谐波传播对网格尺寸的要求。

图6 土柱有限元模型Fig.6 Finite element model of soil pile

3.3 输入地震动

根据实际场地地震安评报告，计算中所施加的地震动加速度时程曲线如图7所示。水平X向

图7 输入的加速度时程曲线Fig.7 Enter acceleration time curve

不同方向随地基深度的加速度幅值，不难看出，如图9所示，新开发模型的计算结果与SuperFLUSH/2D的计算结果在三个方向上总体趋势基本一致，数值上也差别较小。可进一步验证的加速度峰值为0.105g，水平Y向的加速度峰值为0.104g，竖直Z向的加速度峰值为0.124g，时间步长为0.01 s，总持时为40 s。

图9 加速度峰值沿高程的分布图Fig.9 The distribution of Peak acceleration along the height

3.4 计算结果分析

通过对不同高程处的加速度反应谱及加速度幅值随深度的变化规律来考察所开发模型。为方便比较，将所开发模型命名为SANSVE。

图8所示为两种计算模型在不同高程（顶部、中部和底部）处不同方向的加速度反应谱对比，不难发现，新开发模型SANSVE的计算结果与商用软件SuperFLUSH/2D的计算结果整体趋势相同，基本完全吻合，从数值上来看也差别很小，在水平和竖直方向三个高程的反应谱峰值相差最大分别0.000 212g、0.000 189g和0.000 261g。

图8 不同高程处加速度反应谱对比Fig.8 Comparison of acceleration response spectra at different elevations

所开发计算模型的合理有效性。

4 工程应用

以该拟建核电厂所采用AP1000堆型为研究对象，根据水平成层地基的分布特征，进行非岩性地基条件下的核岛厂房结构三维地震响应分析。

4.1 计算模型

核岛厂房结构为质点-梁模型，具体参数可参见文献［8］，地基参数取值如表2所示，而各土层对应的G-γ和D-γ曲线如图5所示。地基模拟范围为从筏板边缘向四周延伸筏板长度的1倍，深度方向则自地表面至玄武岩层顶面取47.7 m。如图10所示所建立核岛厂房结构-地基系统有限元模型和核岛厂房质点-梁模量（带筏板）的局部模型，整体有限元模型共有281 277个单元、253 650个节点。

4.2 计算工况

为更好地研究土体的非线性特征和无限地基辐射阻尼效应，选取三种不同的计算工况进行对比分析，具体如表3所示。

表3 计算工况Table 3 Calculation conditions

限于篇幅，选取包络性节点的计算结果作比较，所选节点位置如图10（a）所示，分别为标高+333.13 m处的PCS水箱顶部（节点48）、标高+169.00 m处的内部结构顶部（节点210）。

图10 核岛厂房结构-地基系统有限元计算模型Fig.10 Structure-foundation system finite element calculation model of nuclear power plant site

4.3 计算结果分析

图11、图12所示为节点48和节点210的加速度反应谱曲线，由于节点位置的不同，两节点的响应分析会有较明显的差异，因此该分析的结果具有很好的代表性。通过对比工况1和工况2可以看出：采用黏性人工边界模拟无限地基辐射阻尼后，核岛的加速度峰值明显降低，尽管两者的加速度反应谱曲线趋势大致相同，但是加速度的量值却降低了，水平方向反应谱加速度峰值最大降低率约为31.2%，竖直方向反应谱加速度峰值最大降低率约为28.6%。

另外，由表4可知，工况1和工况2虽然出现峰值的频率基本相同，但两者的加速度峰值相差很大。这主要是因为刚性地基模型不能吸收散射波的波动能量，对模拟辐射阻尼效应具有一定的局限性。而黏性边界模型可以很好地吸收散射波的逸散能量，从而使核岛结构的响应减小，更加趋于安全。

表4 加速度反应谱的峰值加速度对比Table 4 Peak acceleration comparison of response spectrum

考虑土体非线性特征对核岛结构响应的影响，从图11、图12可以看出，工况2和工况3之间存在差值，较之工况2不考虑土体的非线性特征的影响，工况3中考虑土体非线性效应影响后，加速度峰值显著降低，其峰值频率也移至较低频率。在水平和竖直两个方向的加速度幅值最大分别下降0.3g和0.4g，峰值频率最大分别偏移了0.9 Hz和4.5 Hz。这主要是因为考虑土体非线性特征后，地基刚度变小，而阻尼比会增大，从而导致地基主频减小，反应谱峰值也随之减小。

图11 不同工况条件下节点48处加速度反应谱对比Fig.11 Comparison of acceleration reflection spectra under different working conditions of node 48

图12 不同工况条件下节点210的加速度反应谱对比Fig.12 Comparison of acceleration reflection spectra under different working conditions of node 210

由工况1和工况3的对比可知，考虑黏性边界和近场非线性时，加速度反应谱有着非常明显的变化，水平向加速度幅值最大下降了0.7g，峰值频率最大偏移了0.9 Hz。竖直向加速度幅值最大下降了0.7g，峰值频率最大偏移了4.5 Hz。反应谱峰值最大减小幅度约为46.7%，减小幅值约为0.7g，且加速度峰值频率向较低频偏移。进一步说明，在核岛厂房结构软土地基抗震安全性分析中，考虑无限地基辐射阻尼效应及土体非线性的必要性。

5 结 论

（1）基于UPFs的二次开发特点，在ANSYS进行非岩性地基条件下核岛厂房结构三维地震响应分析计算模型的开发具有可行性，经算例验证，该模型表现出良好的合理性及精确性，可操作性强，具有较高的工程实用价值。

（2）非岩性地基对核岛厂房结构的地震响应有较明显的影响，在抗震安全评价中考虑近场地基非线性和远场无限地基辐射阻尼是十分必要的。

（3）基于通用有限元软件的开发，可在满足专业要求的前提下，最大限度地发挥软件的优势，利用其丰富的单元类型及强大的求解能力，解决更复杂的实际工程问题。

（4）所开发模型具有直接法的优点，可对不符合抗震设防要求的地基进行加固处理及优化数值模拟，考虑桩-土-结构的动力相互作用，进一步对地基处理后的方案进行安全性评价。

致谢 感谢同济大学建筑设计研究院(集团)有限公司陈曦高级工程师在本项目中对作者的悉心指导和帮助。