非岩性地基钍基熔盐堆厂房楼层反应谱分析与优化

Analysis and Optimization of Floor Response Spectra of Thorium-based Molten Salt Reactor Buildings Built on Non-Lithological Foundations

摘要

本文提出了匹配反应谱和包络功率谱的人工地震波混合模拟方法；详细给出了水平及竖向地震作用下人工边界5个面等效应力大小及方向；利用EERA软件进行了土层反应计算，得到自由场不同场地深度的地震反应；建立了钍基熔盐堆厂房及土体的有限元模型，进行了考虑土-结构相互作用效应的地震响应分析；提出了竖向楼层反应谱优化方法。结果表明：提出的人工地震波混合模拟方法有较高的目标设计反应谱匹配精度，且能很好地包络目标功率谱密度；相对于水平方向，由于楼板局部振动特性不同，竖直方向楼板反应谱会在楼板局部振动固有频率附近出现峰值，对核设备抗震分析带来不利影响；通过调整结构局部布置，改变楼板局部振动特性，可以改善该部分楼板的竖向反应谱。

Abstract

A hybrid simulation method for artificial seismic waves was developed to match target response spectra while enveloping specified power spectral densities. The equivalent stress magnitudes and orientations on five artificial boundary surfaces were quantified under both horizontal and vertical seismic excitations. Soil response analyses were conducted using EERA software to compute free-field seismic motions at varying depths. A finite element model integrating the thorium-based molten salt reactor （TMSR） building with underlying soil strata was established， and seismic response analyses incorporating soil-structure interaction （SSI） effects were performed. Additionally， a vertical floor response spectrum （VFRS） optimization methodology was proposed. The results demonstrate that： The hybrid simulation method achieves high accuracy in matching target design response spectra and effectively envelops target power spectral densities； Due to distinct local vibrational characteristics of floor slabs， the vertical response spectrum may exhibit amplification peaks near local natural frequencies—unlike horizontal spectra—which adversely impacts seismic assessments of nuclear equipment； Adjusting the structural layout to modify local vibrational characteristics improves the vertical response spectrum of the targeted floor slab sections.

Graphical abstract

关键词

钍基熔盐堆 / 人工地震波混合模拟 / 地震波动输入 / 土-结构相互作用 / 楼层谱优化

Key words

thorium molten salt reactor / mixed simulation of artificial seismic waves / earthquake wave input / soil-structure interaction / floor spectrum optimization

0 引 言

钍基熔盐堆（TMSR）作为我国第四代核能技术的战略部署，具有钍资源高效利用、核废料少、本质安全及高温输出等显著优势，特别契合我国钍资源丰富特点。目前上海建工集团参与的2 MW实验堆已建成，10 MW研究堆选址甘肃民勤红砂岗，其30 m深核岛厂房地基为砾质砂岩非岩性地层（平均剪切波速约600 m/s）^［1］。

楼层谱作为核设施抗震设计的核心依据，在非岩性场地建设中必须考虑土-结构相互作用（SSI）。国际主流采用SASSI^［2］等专用软件进行SSI分析，近年来利用通用有限元三维建模的直接法^［3-4］也得到广泛应用。研究堆抗震分析需突破两大关键技术：①构建匹配设计反应谱且包络目标功率谱的人工时程；②建立包含SSI效应的厂房动力时程计算模型。

核电站抗震设计要求人工地震波需匹配目标设计反应谱（TDRS）并包络目标功率谱密度函数（TPSD）^［1］。现有方法可分为时域、频域和混合三类，均通过迭代优化种子加速度时程的反应谱匹配度来实现。时域方法通过在种子加速度时程上叠加时间历程提升匹配精度，如Nie等^［5］提出的贪婪小波方法（GWM）通过单次迭代添加小波显著提高了计算效率。频域方法基于傅里叶级数生成种子加速度时程，然后迭代修正傅里叶级数以便生成与TDRS匹配的人工加速度时程，虽计算快速但精度较低。混合方法结合两者优势，如Zhang等^［6］和何佳等^［7］分别通过窄带过程模型改进了时程修正效果。在TPSD包络研究方面，王海涛等^［8］基于窄带过程实现了TDRS、TPSD双参数匹配。人工智能技术在人工地震波生成领域的应用日益广泛，如陈晓明等^［9］通过扩散模型实现了双参数匹配的人工地震波合成。

在土-结构相互作用（SSI）的核电厂房动力时程分析中，人工边界设置与土-结构耦合建模是两大核心问题。针对有限土体模型无法反映无限地基辐射阻尼效应的问题，需通过人工边界模拟波动能量耗散，并施加等效输入荷载确保内部波动场与自由场协调。现有边界条件可分为全局边界和局部边界两类，全局边界（如边界元法^［10］、无穷元法^［11］）虽满足无限域场方程和辐射条件^［12］，但计算成本过高，多限于线性分析；局部边界则通过物理效应模拟波动传播，典型代表包括黏性边界及其改进型黏弹性边界^［13］。刘晶波团队^［14］针对传统黏性边界的低频稳定性不足，提出黏弹性边界理论，谷音等^［15］发展出适用于二维/三维的等效黏弹性人工边界单元。在地震输入方面，刘晶波等^［14］将地震动转化为人工边界上的等效荷载，实现了波动的输入。

土-结构耦合建模方法按理论基础可分为三类：集中质量参数法通过弹簧-质量-阻尼系统表征地基动力特性^［17］，子结构法通过子系统耦合求解地震响应^［18］，直接法则将土-结构作为整体建模，兼具非线性分析与三维模拟优势^［3-4］。其中直接法采用板壳单元描述厂房、实体单元构建土体，通过共节点或Coupling方式实现多尺度力传递^{［3，19-20］}，成为核电站精细化分析的主流方法。

现有人工地震波模拟方法在实现设计反应谱（TDRS）与目标功率谱（TPSD）高精度双参数匹配时，往往难以兼顾计算效率与收敛稳定性；同时，多向地震作用下人工边界等效应力的计算方法，以及竖向地震动与楼层反应谱的耦合优化机制尚未完善。为了有效进行10 MW钍基熔盐研究堆厂房的楼面加速度反应谱分析，本文提出匹配反应谱和包络功率谱密度的人工加速度时程（人工地震波）混合模拟方法，水平及竖向地震作用下人工边界等效输入荷载具体计算方法及方向，楼层竖向反应谱优化方法，为钍基熔盐研究堆核设备抗震设计和选型打好基础。

1 人工加速度时程的混合模拟

为了生成匹配TDRS和包络TPSD的人工加速度时程，本文提出由基于频域方法的模拟和基于时域方法的混合模拟组成的时频域（混合）模拟方法。对于模拟人工加速度时程，频域调整算法使用简便但匹配TDRS的精度较低，而时域调整算法（如窄带叠加法）匹配TDRS的精度较高但收敛速度较低。本文提出的频域-时域混合模拟算法兼具时域调整算法匹配精度高和频域调整算法使用便捷及收敛速度快的优点。

1.1　人工加速度时程时频域混合模拟方法

1） 频域模拟方法

频域方法利用近似匹配TDRS的PSD，由傅里叶级数（截断）展开生成平稳加速度时程，对平稳加速度时程进行幅值调制得到种子加速度时程，不断迭代修正傅里叶系数的方法来匹配TDRS^［6］，生成初步匹配TDRS的人工加速度时程。频域方法步骤如下：

（1） 生成种子地震波

令TDRS为S_{\mathrm{T}}\left(\omega \right)，则由PSD和TDRS的近似转换关系，可以得到近似匹配S_{\mathrm{T}}\left(\omega \right)的双边PSD：

G_{\mathrm{z}}\left(\omega \right)=\frac{\zeta }{\mathrm{\pi }\omega }\frac{S_{\mathrm{T}}^{\mathrm{2}}\left(\omega \right)}{\mathrm{l}\mathrm{n}\left(\omega T_{\mathrm{0}}\right)-\mathrm{l}\mathrm{n}\left[-\mathrm{\pi }\mathrm{l}\mathrm{n}\left(\mathrm{1}-r\right)\right]}（1）

式中：G_{\mathrm{z}}\left(\omega \right)为人工加速度时程的双边PSD；\omega为圆频率；\zeta为阻尼比；T_{\mathrm{0}}为人工动时程持续时间；r为超越G_{\mathrm{z}}\left(\omega \right)的概率，一般取r\le \mathrm{0.15}。

由傅里叶级数（截断）展开生成平稳加速度时程：

x^{(\mathrm{0})}\left(t\right)=\sqrt[]{\mathrm{2}}{\displaystyle \sum _{i=\mathrm{0}}^{N-\mathrm{1}}}A\left({\omega }_i\right)\mathrm{c}\mathrm{o}\mathrm{s}\left({\omega }_{i}t+{\phi }_i\right)（2）

式中：N为傅里叶级数截断数；幅值（傅里叶系数）A\left({\omega }_i\right)=\sqrt[]{\mathrm{2}{G}_{\mathrm{z}}\left({\omega }_i\right)\mathrm{\Delta }\omega }；圆频率{\omega }_i=i\mathrm{\Delta }\omega；\mathrm{\Delta }\omega ={\mathrm{2}\mathrm{\pi }/T}_{\mathrm{0}}为频率采样间隔；初相位{\phi }_i(i=\mathrm{1},\cdots ,N)为\left[\mathrm{0,2}\mathrm{\pi }\right]上均匀分布的随机数。

对平稳时程x^{(\mathrm{0})}\left(t\right)进行幅值调制，可以得到种子人工加速度时程：

a^{\left(\mathrm{0}\right)}\left(t\right)=m\left(t\right)x^{(\mathrm{0})}\left(t\right)（3）

式中：m\left(t\right)为地震波强度包络函数，本项目由地震安评报告给出：

（4）

其中的参数值分别为t_{\mathrm{1}}=\mathrm{2.5}\mathrm{s}、t_{\mathrm{2}}=\mathrm{17.5}\mathrm{s}和c=\mathrm{0.1}。

（2） 地震波迭代修正

一般来说，a^{\left(\mathrm{0}\right)}\left(t\right)的RS不能满足TDRS的匹配要求，需要进行迭代逼近，使人工加速度时程匹配TDRS。修正a^{\left(\mathrm{0}\right)}\left(t\right)的频域方法是迭代修正式（2）中的幅值A\left({\omega }_i\right)：

（5）

式中：A^{(l)}\left({\omega }_i\right)、A^{(l+\mathrm{1})}\left({\omega }_i\right)分别为第l、l+\mathrm{1}次迭代时的幅值；S_a^{\left(l\right)}\left({\omega }_j\right)为第l次迭代生成的人工加速度时程a^{\left(l\right)}\left(t\right)的RS；{\omega }_j为频率控制点处的频率；N_{j,l}、N_{j,r}为控制频段\left[{\omega }_{j,l},{\omega }_{j,r}\right]的左右端点号：

{\omega }_{j,l}=\frac{{\omega }_{j-\mathrm{1}}+{\omega }_j}{\mathrm{2}},{\omega }_{j,r}=\frac{{\omega }_j+{\omega }_{j+\mathrm{1}}}{\mathrm{2}}（6）

N_{j,l}=\frac{{\omega }_{j,l}}{\mathrm{\Delta }\omega },N_{j,r}=\frac{{\omega }_{j,r}}{\mathrm{\Delta }\omega }（7）

注意，由式（7）计算的N_{j,l}和N_{j,r}应为正整数。

2） 匹配TPSD方法

为了避免地震动在某些频带处能量不足，需对人工加速度时程进行调整，使人工时程PSD包络80%的TPSD，方法如下：

（1） 对第l轮迭代得到的人工时程a^{\left(l\right)}\left(t\right)进行离散傅里叶变换，可以得到离散频率点{\omega }_i处的傅里叶幅值\left|{\widehat{F}}_g^{\left(l\right)}\left({\omega }_i\right)\right|\left(i=\mathrm{0,1},\cdots ,N-\mathrm{1}\right)，假设强度包络随机过程a^{\left(l\right)}\left(t\right)的边际功率谱与其平稳段功率谱相等，可以近似计算a^{\left(l\right)}\left(t\right)的PSD：

{\widehat{G}}_g^{\left(l\right)}\left({\omega }_i\right)=\frac{T_{\mathrm{0}}}{{\int }_{\mathrm{0}}^{T_{\mathrm{0}}}{m}^{\mathrm{2}}\left(t\right)\mathrm{d}t}\frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }{T}_{\mathrm{0}}}{\left|{\widehat{F}}_g^{\left(l\right)}\left({\omega }_i\right)\right|}^{\mathrm{2}}（8）

（2） 增量PSD计算式：

\mathrm{\Delta }{G}_g^{\left(l\right)}\left({\omega }_i\right)=\mathrm{m}\mathrm{a}\mathrm{x}\left[\mathrm{0.8}{G}_{\mathrm{T}}\left({\omega }_i\right)-{\widehat{G}}_g^{\left(l\right)}\left({\omega }_i\right),\mathrm{0}\right]（9）

式中：G_{\mathrm{T}}\left(\omega \right)为TDRS对应的TPSD。当项目地震安评报告未提供TPSD时，可由下式给出^［1］：

G_{\mathrm{T}}\left(\omega \right)={\alpha }_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{2}}{G}_{\mathrm{0}}\left(\omega \right){\left[\frac{S_{\mathrm{T}}\left(\omega \right)}{S_{\mathrm{0}}\left(\omega \right)}\right]}^{\mathrm{2}}（10）

其中：S_{\mathrm{T}}\left(\omega \right)为TDRS；{\alpha }_{\mathrm{m}\mathrm{a}\mathrm{x}}为S_{\mathrm{T}}\left(\omega \right)的地震影响系数最大值；S_{\mathrm{0}}\left(\omega \right)、G_{\mathrm{0}}\left(\omega \right)分别为美国核管会导则RG1.60水平向标准反应谱以及该标准反应谱对应的（双边）PSD（以{\alpha }_{\mathrm{m}\mathrm{a}\mathrm{x}}=1.0标定）^［1］。

（3） 对\mathrm{\Delta }{G}_g^{\left(l\right)}\left(\omega \right)进行离散逆傅里叶变换，可以得到增量平稳时间历程\mathrm{\Delta }{x}^{\left(l\right)}\left(t\right)，对\mathrm{\Delta }{x}^{\left(l\right)}\left(t\right)进行幅值调制，得到增量人工加速度时程：

\mathrm{\Delta }{a}_g^{(l)}=m\left(t\right)\mathrm{\Delta }{x}^{\left(l\right)}\left(t\right)（11）

式中：m\left(t\right)为地震波幅值包络函数。

（4） 将增量人工加速度时程\mathrm{\Delta }{a}_g^{(l)}叠加到a^{\left(l\right)}\left(t\right)，则得到包络80%TPSD的加速度时程。

3） 时域模拟方法

为得到非常高的反应谱匹配度，在时域中继续修正由频域方法生成的人工加速度时程a^{\text{'}}\left(t\right)。当频域方法频率控制点处加速度时程反应谱与设计反应谱的比值属于［0.9，1.3］范围以后，进入加速度时程的时域调整。本文采用称为窄带叠加法的时域方法进一步修正a^{\text{'}}\left(t\right)。按照国内外规范要求，人工地震动时程需同时匹配多个阻尼比的TDRS。拟合出同时匹配多个阻尼比TDRS的人工地震动时程是比较困难的，特别是还需包络TPSD，利用不同阻尼反应谱谱值在时域中发生时刻可能不同这一性质，通过在时域中叠加不同频率的窄带时程，可以减少不同阻尼反应谱之间相互影响，做到同时匹配3个阻尼比的TDRS^［21］。窄带叠加法的原理如下^［6-9］：

令以a^{\text{'}}\left(t\right)为种子时程经l-\mathrm{1}轮（叠加窄带增量激励过程）修正的时间历程{a^{\text{'}}}^{\left(l\right)}\left(t\right)的RS为S_{a^{\text{'}}}^{\left(l\right)}\left(\omega \right)，\omega为圆频率，则其在频率控制点{\omega }_j处与TDRS的差值为\mathrm{\Delta }{S}^{\left(l\right)}\left({\omega }_j\right)=S_{\mathrm{T}}\left({\omega }_j\right)-S_{a^{\text{'}}}^{\left(l\right)}\left({\omega }_j\right)。令在{a^{\text{'}}}^{\left(l\right)}\left(t\right)激励下有阻尼固有圆频率为{\omega }_j、阻尼比为\zeta的线性振子的绝对加速度反应为a_{a,j}^{\left(l\right)}\left(t\right)，t_{\mathrm{m}\mathrm{a}\mathrm{x},j}\le T_{\mathrm{0}}为\left|a_{a,j}^{\left(l\right)}\left(t\right)\right|出现最大值的时刻，T_{\mathrm{0}}为人工加速度时程的持续时间，则在a_{a,j}^{\left(l\right)}\left(t\right)上叠加窄带过程

\begin{array}{l}\mathrm{\Delta }{a}_{a,j}^{\left(l\right)}\left(t\right)=\mathrm{s}\mathrm{g}\mathrm{n}\left[a_{a,j}^{\left(l\right)}\left(t_{\mathrm{m}\mathrm{a}\mathrm{x},j}\right)\right]\mathrm{\Delta }{S}^{\left(l\right)}\left({\omega }_j\right)\\ \mathrm{e}\mathrm{x}\mathrm{p}\left[-\frac{\mathrm{1}}{\mathrm{2}}{\left(\frac{t-t_{\mathrm{m}\mathrm{a}\mathrm{x}}}{\sigma }\right)}^{\mathrm{2}}\right]\mathrm{c}\mathrm{o}\mathrm{s}{\omega }_j\left(t-t_{\mathrm{m}\mathrm{a}\mathrm{x}}\right)\end{array}（12）

使得总绝对加速度反应绝对值在t=t_{\mathrm{m}\mathrm{a}\mathrm{x}}处最大且等于S_{\mathrm{T}}\left({\omega }_j\right)，其中的\sigma为控制窄带过程\mathrm{\Delta }{a}_{a,j}^{\left(l\right)}\left(t\right)带宽的系数：

\sigma =\frac{\mathrm{\pi }}{{\omega }_{c,j}}\frac{\mathrm{1}}{N_{c,j}}（13）

式中：{\omega }_{c,j}=\mathrm{m}\mathrm{i}\mathrm{n}\left(\frac{{\omega }_j-{\omega }_{j-\mathrm{1}}}{\mathrm{2}},\frac{{\omega }_{j+\mathrm{1}}-{\omega }_j}{\mathrm{2}}\right)，由频率控制点间距确定；参数N_{c,j}由试错法确定为常数3。

对应于增量反应过程\mathrm{\Delta }{a}_{a,j}^{\left(l\right)}\left(t\right)的增量激励过程为

\mathrm{\Delta }{a}^{\text{'}}{}_j^{\left(l\right)}\left(t\right)=\mathrm{R}\mathrm{e}\left\{{\mathrm{ℱ}}^{-\mathrm{1}}\left[\frac{\mathrm{\Delta }{A}_{a,j}^{\left(l\right)}\left(\omega \right)}{H_{a,j}\left(\omega \right)}\right]\right\}（14）

式中：H_{a,j}\left(\omega \right)代表所考虑的线性振子绝对加速度反应的传递函数；\mathrm{\Delta }{A}_{a,j}^{\left(l\right)}\left(\omega \right)=\mathrm{ℱ}\left[\mathrm{\Delta }{a}_{a,j}^{\left(l\right)}\left(t\right)\right]，为\mathrm{\Delta }{a}_{a,j}^{\left(l\right)}\left(t\right)的傅里叶变换；\mathrm{R}\mathrm{e}\left(\cdot \right)代表实部；{\mathrm{ℱ}}^{-\mathrm{1}}\left(\cdot \right)代表傅里叶逆变换。

在第l轮迭代步完成频率控制点{\omega }_j处TDRS匹配的时间历程为a^{\text{'}}{}_{}^{\left(l\right)}\left(t\right)+\mathrm{\Delta }{a}^{\text{'}}{}_j^{\left(l\right)}\left(t\right)。当对所有频率控制点进行窄带增量响应过程处理后，可获得经一次时域调整的时间历程a^{\text{'}}{}_{}^{\left(l+\mathrm{1}\right)}\left(t\right)。

在进行人工加速度时程的模拟时，先由频域方法和TPSD包络方法组成的混合模拟方法，经多次迭代产生初步匹配TDRS和包络TPSD的人工加速度时程，为进一步提高与TDRS的匹配度，再通过由时域方法和TPSD包络方法组成的混合模拟方法，经过多次迭代产生与TDRS的匹配度更高同时又包络TPSD的人工加速度时程。

1.2　地面人工加速度时程

根据地震安评报告^［22］，钍基熔盐研究堆SL-2级（100年超越概率为1%）阻尼比为5%的场地自由地表水平、竖直方向TDRS特征值见表1，谱值与周期之间为线性对数关系。与该TDRS对应的近似PSD和TPSD由公式（1）和（10）计算。

采用本文提出的混合模拟方法产生的匹配100年超越概率1%TDRS和包络TPSD的2条水平、竖向人工加速度时程如图1所示，人工地震波的时间间隔\mathrm{\Delta }t=\mathrm{0.01}\mathrm{s}，总持续时间T_{\mathrm{0}}=\mathrm{35}\mathrm{s}。按照要求^［1，23］，在0.2~33 Hz频率范围内，反应谱频率控制点数不得少于75个；在0.1~50 Hz频率范围内，所需要的控制点数不得少于100个，上述频率点应大体均匀分布在对数坐标下；本文频率控制点分别为111、120个，且满足规范^［1］表3.4.1频率增量的要求。图2和3绘出了人工加速度时程与TDRS的匹配结果和包络80%的TPSD的结果。需要指出的是，为了使钍基熔盐堆厂房楼层反应谱的分析结果偏于保守一些，图3中的加速度时程功率谱为频率控制点处的功率谱数值，而不是核电站设计规范^［1］要求的±0.2f范围内的平均（平滑处理）值，通常基于频率控制点处功率谱差值调整人工加速度时程将具有稍大的能量（相比平滑处理）。对于全部频率控制点，水平和竖向人工加速度时程的RS小于TDRS的控制点各有1个，不多于5个且不低于TDRS的90%；所有控制点RS不超过TDRS的1.3倍^［1］。0.3~24 Hz范围内的频率控制点，2条人工加速度时程的PSD都大于相应的80%的TPSD值，满足规范^［1］要求。

2 自由场不同深度土层地震反应计算

在强震作用下，土体属性一般认为是非线性的，当土体剪应变γ_d增加时，土体的动剪切模量比G_d/G_max（G_d为动剪切模量，G_max为最大剪切模量）呈下降趋势，阻尼比ζ呈上升趋势。需拟合动剪切模量比、阻尼比与剪应变关系的曲线，从而通过一维等效线性化方法得到地震作用下的动剪切模量及阻尼比，并通过场地地震动反应计算得到不同场地深度的加速度、速度、位移及应力时程，作为波动法的输入条件。本文采用EERA软件^［23］进行土层反应计算，由上述拟合得到的地表人工加速度时程，得到自由场不同深度的地震反应。

为确定研究堆岩土力学试验参数，中国科学院武汉岩土力学研究所对10 MW钍基熔盐研究堆场地钻孔岩芯进行了室内力学试验，通过动三轴和共振柱试验得到动剪切模量比、动弹性模量比、阻尼比与剪应变和轴向应变的关系，表2为动剪切模量比、阻尼比随剪应变变化的试验数据^［22］。图4展示了基于表2深度52~58 m土层试验数据，通过EERA软件拟合得到的该深度强风化正长花岗岩（②号土）动剪切模量比、阻尼比与剪应变关系曲线，其余深度土层的关系曲线均可通过相同方法完成拟合。

3 波动法等效应力时程计算

刘晶波^［14］等提出的波动法将地震动视为作用于人工边界上的等效荷载，从而实现地震动的输入。当人工边界底面垂直入射剪切波（SV波）时，SV波为横波，振动方向与波前进方向垂直，产生水平方向振动（水平地震）；人工边界底面垂直入射纵波（P波）时，振动方向与波前进方向平行，产生竖直方向振动（竖向地震）。现对水平、竖向地震作用下人工边界5个面（前后左右底）上等效荷载进行推导。

人工边界是从半无限自由场中取出有限的计算区域，在人工边界内部面输入等效应力时程，准确实现波动输入的条件是在人工边界上施加的等效应力时程荷载应使人工边界上的位移和应力与原自由场相同，即满足下式^［14］：

\begin{array}{l}\delta (x_B,y_B,z_B,t)={\delta }_{\mathrm{0}}(x_B,y_B,z_B,t)\\ \tau (x_B,y_B,z_B,t)={\tau }_{\mathrm{0}}(x_B,y_B,z_B,t)\end{array}（15）

式中：\delta \mathrm{、}\tau和{\delta }_{\mathrm{0}}\mathrm{、}{\tau }_{\mathrm{0}}分别为人工边界内部节点处（坐标为x_B,y_B,z_B，时间t）位移、应力和自由场位移、应力。

3.1　底部垂直入射SV波（水平地震）

如图5所示底面人工边界，设X、Y轴水平，Z垂轴竖直向下，沿X、Y、Z轴位移u\mathrm{、}v\mathrm{、}w。底部垂直入射剪切波（SV），传播方向-Z，土体产生X向水平振动。由于振动方向为X向，因此不存在沿Y向剪应力。假设在人工边界（并联弹簧和阻尼）内部点B施加的应力时程为F_B(t)，采用力学分析中的脱离体概念，在人工边界内部点B处将人工边界与土体截断（图5），图中F_B(t)为人工边界对土体的作用力，因此对于土体微小脱离体\mathrm{\Sigma }{F}_x=\mathrm{0}：

\begin{array}{l}{F}_B(t)=\tau (x_B,y_B,z_B,t)+f_B(t)=\\ {\tau }_{\mathrm{0}}(x_B,y_B,z_B,t)+f_B(t)\end{array}（16）

对于人工边界脱离体，由弹簧和阻尼器构成的物理元件运动方程为

\begin{array}{l}{f}_B(t)=C_B\dot{\delta }(x_B,y_B,z_B,t)+K_B\delta (x_B,y_B,z_B,t)=\\ C_B{\dot{\delta }}_{\mathrm{0}}(x_B,y_B,z_B,t)+K_B{\delta }_{\mathrm{0}}(x_B,y_B,z_B,t)\end{array}（17）

式中：C_B和K_B分别为分布阻尼器的阻尼系数和分布弹簧的刚度系数。将式（17）代入（16）可得^［14］

\begin{array}{l}{F}_B(t)={\tau }_{\mathrm{0}}(x_B,y_B,z_B,t)+\\ C_B{\dot{\delta }}_{\mathrm{0}}(x_B,y_B,z_B,t)+K_B{\delta }_{\mathrm{0}}(x_B,y_B,z_B,t)\end{array}（18）

上式等式右边分为两部分：一部分为自由场应力项{\tau }_{\mathrm{0}}（第一项），它是自由场地震动在人工边界面上的剪应力分布；另外一部分为速度位移项C_B{\dot{\delta }}_{\mathrm{0}}\mathrm{、}{K}_B{\delta }_{\mathrm{0}}（第二、三项），是用来平衡人工边节点速度和位移引起的附加应力。

3.2　底部垂直入射P波（竖向地震）

当P波底部垂直入射时，如图6所示底面人工边界，传播方向-Z，土体产生Z向振动。由\mathrm{\Sigma }{F}_z=\mathrm{0}，同理可得

\begin{array}{l}{F}_B(t)={\sigma }_{z\mathrm{0}}(x_B,y_B,z_B,t)+\\ C_B{\dot{\delta }}_{\mathrm{0}}(x_B,y_B,z_B,t)+K_B{\delta }_{\mathrm{0}}(x_B,y_B,z_B,t)\end{array}（19）

其中，{\sigma }_{z\mathrm{0}}为自由场地震动在人工边界面上的Z向正应力分布。

由于对称（任一竖向平面都是对称面），有水平位移u=\mathrm{0},v=\mathrm{0}及水平轴向应变{\epsilon }_x=\mathrm{0},{\epsilon }_y=\mathrm{0}。由弹性力学物理方程

\left.\begin{array}{l}{\sigma }_x=(\mathrm{2}G+\lambda ){\epsilon }_x+\lambda {\epsilon }_y+\lambda {\epsilon }_z\\ {\sigma }_y=(\mathrm{2}G+\lambda ){\epsilon }_y+\lambda {\epsilon }_x+\lambda {\epsilon }_z\\ {\sigma }_z=(\mathrm{2}G+\lambda ){\epsilon }_z+\lambda {\epsilon }_x+\lambda {\epsilon }_y\end{array}\right\}（20）

可得

{\sigma }_x={\sigma }_y=\frac{\lambda }{\mathrm{2}G+\lambda }{\sigma }_z（21）

式中：第一拉梅常数\lambda =\frac{\mu E}{(\mathrm{1}+\mu )(\mathrm{1}-\mathrm{2}\mu )}；剪切模量G=\frac{\mu }{\mathrm{2}(\mathrm{1}+\mu )}；E为弹性模量；\mu为泊松比。

3.3　人工边界面等效应力大小及方向

由以上分析，表3给出了相应坐标系底部垂直入射SV波（横波，地震动水平X、Y向）和P波（纵波，地震动竖直Z向）时人工边界5个面需要考虑的应力方向，5个面为面向表中所示坐标系所对应的前、后、左、右、底面。表中速度位移项是人工边界位移、速度时程引起弹簧和阻尼的附加应力，自由波场不同深度土层的位移、速度时程由上节EERA软件求出，乘以相应的K_B和C_B后得到弹簧和阻尼的附加应力。地震水平输入应力项：人工边界前后左右底面的应力为剪应力，与自由场该深度水平剪应力大小相等。如表3中所示X向地震输入，振动方向X向，左右面垂直于X轴，如图5所示应力为{\tau }_{xz}，大小与该深度水平剪力{\tau }_{zx}相等，方向分别为Z-、Z+；振动方向为X向，不存在沿Y向剪应力，因此前后面剪应力为0。地震竖直输入应力项：人工边界5个面应力为正应力，底部正应力与自由场该深度正应力大小相等，左右前后面正应力与自由场该深度正应力满足式（21）。其中，X、Y、Z代表三个坐标轴，+ -表示沿坐标轴的正或负方向输入。

4 10 MW钍基熔盐研究堆SSI计算

10 MW钍基熔盐研究堆核岛厂房长86 m，宽80 m，其中地下部分长86 m，宽47 m。地上三层（局部二层），屋面标高为+20 m，地下二层（局部有夹层），地下室底板标高为-27 m。采用剪力墙结构体系，反应堆容器等关键核设备位于地下室，因防辐射、承载力要求，地下室部分剪力墙、楼板厚度较厚。

按照规范^［1］，土-结构相互作用（SSI）分析可采用三类数值模型：集中质量参数模型、有限元模型（直接法）、子结构模型。地基近似为均匀弹性介质时，可采用集中参数模型，模拟出地基的等效弹簧刚度和等效阻尼系数，对地基采用弹簧阻尼器进行简化；非匀质地基宜采用有限元模型。本文采用有限元模型将土-结构作为整体进行精细化建模，研究堆的有限元模型如图7所示，该模型由厂房结构、土体及外围的人工边界组成。

厂房结构采用梁单元和板壳单元，其中梁、柱采用梁单元，楼板、剪力墙采用板壳单元。厂房结构有限元模型可以通过YJK转通用有限元接口程序转换生成inp文件，导入到通用有限元程序生成厂房模型后，利用网格—偏移创建实体层的方法从底板和地下室外墙板壳单元生成共节点土体三维实体单元。YJK里输入的隔墙、设备等荷载，接口程序自动转化为有限元软件的点质量，极大提高了输入效率。

土壤-结构共同作用（SSI）数值模拟分析模型中，地基是一个半无限开放系统，由于计算技术和时间成本的限制，只能从地基中截取结构周围的有限土体，截断边界以黏弹性人工边界模拟真实土体的辐射阻尼。土体模型为一个长616 m、宽330 m、高58 m的长方体，侧边界至厂房同侧边缘的距离大于同一方向厂房尺寸的3倍，底边界距基础底板的距离小于基础埋深的3倍，但底边界土体的动弹性模量E_d已经不小于土体表层动弹性模量的10倍，底部边界可视为刚性边界^［1］，为确保计算模型准确性，所有边界仍采用黏弹性人工边界。人工边界采用三维一致等效黏弹性人工边界^［15］，在土体模型的4个侧面和底面延伸着一层厚度为0.5 m、质量为0的三维实体单元，其等效弹性模量、等效剪切模量、等效泊松比、法向和切向黏弹性人工边界修正系数按文献［15］取值。

进行楼层谱计算时，根据规范^［1］要求，极限安全地震动（100年超越概率1%）对应的阻尼比，钢筋混凝土结构取7%，管道系统取4%，设备一般取3%。根据规范^［24］，考虑混凝土开裂影响，采用混凝土弹性模量的80%进行计算。土体材料参数采用EERA软件一维土柱等效线性化迭代计算获得的动弹性模量、动剪切模量、动阻尼比。在土体三维实体单元与人工边界相邻的面输入第3节求得的等效面力时程（也可转化为该土体三维实体单元的体力时程进行输入），进行SSI动力响应计算时，采用隐式动力分析，Rayleigh（瑞丽）阻尼。不同材料（包括混凝土、设备、土层、黏弹性人工边界）的阻尼均以阻尼比的形式给出，有限元计算时需把阻尼比转换为瑞丽阻尼，其公式为^［25］

\left.\begin{array}{l}\alpha =\mathrm{2}{\omega }_{\mathrm{1}}{\omega }_{\mathrm{2}}\left({\xi }_{\mathrm{1}}{\omega }_{\mathrm{2}}-{\xi }_{\mathrm{2}}{\omega }_{\mathrm{1}}\right)/\left({{\omega }_{\mathrm{2}}}^{\mathrm{2}}-{{\omega }_{\mathrm{1}}}^{\mathrm{2}}\right)\\ \beta =\mathrm{2}\left({\xi }_{\mathrm{2}}{\omega }_{\mathrm{2}}-{\xi }_{\mathrm{1}}{\omega }_{\mathrm{1}}\right)/\left({{\omega }_{\mathrm{2}}}^{\mathrm{2}}-{{\omega }_{\mathrm{1}}}^{\mathrm{2}}\right)\end{array}\right\}（22）

式中：\alpha为质量阻尼系数；\beta为刚度阻尼系数；{\omega }_{\mathrm{1}}~{\omega }_{\mathrm{2}}为选定的固有圆频率范围；{\xi }_{\mathrm{1}}、{\xi }_{\mathrm{2}}为对应{\omega }_{\mathrm{1}}\mathrm{、}{\omega }_{\mathrm{2}}的材料阻尼比（通常假定不同频率对应的阻尼比相同）。确定瑞丽阻尼比的原则是：在{\omega }_{\mathrm{1}}和{\omega }_{\mathrm{2}}处等于材料阻尼比，\left[{\omega }_{\mathrm{1}},{\omega }_{\mathrm{2}}\right]要覆盖结构分析中感兴趣的频段：在该频段范围以内，瑞丽阻尼比略小于材料阻尼比，结构反应将略大于实际响应（偏安全）；在该频段范围以外，瑞丽阻尼比迅速增大，振动响应被有效抑制。本文通过模态分析中质量参与系数来确定{\omega }_{\mathrm{1}}\mathrm{、}{\omega }_{\mathrm{2}}，分别取质量参与系数较大的低阶和高阶振型频率，同时控制\left[{\omega }_{\mathrm{1}},{\omega }_{\mathrm{2}}\right]频段累积质量参与系数\ge \mathrm{90}\mathrm{\%}。

5 10MW钍基熔盐研究堆楼板反应谱优化思路

由于10 MW钍基熔盐堆主要核岛设备都放在地下，主要对100年超越概率1%，阻尼比5%的地下部分楼板反应谱进行初步优化。楼层的响应主要取决于厂房楼板该方向上的自振特性。由于楼板水平刚度大，因此各层楼板的水平响应与该层厂房整体水平响应类似，地下楼层受土体缓冲作用的影响更为明显，反应谱峰值向低频偏移，每层楼板不同位置处的水平反应谱差别不大^［26］。楼板的竖向刚度在跨中和支座差别很大，引起楼板局部振型频率差异，地下楼层每层楼板不同位置处的竖向反应谱有较大差别。

1号点与2号点的位置如图8所示：其中1点为-18 m楼层板跨中点，2点为柱支座点。二者地震反应谱的对比如图9所示：该层这两点水平方向反应谱基本一致，体现出整体水平相应特性；由于土的缓冲作用，放大作用主要集中在1~4 Hz低频范围；由于处于地下，对地震动的放大效应较小，水平谱值也小于地表目标谱。这两点在竖向振动特性上的差异导致竖向反应谱差别较大，板跨中点1在楼板局部振动固有频率14 Hz附近出现峰值11.1 m/s²，已经超过了地表目标谱最大值9.3 m/s²。优化楼板竖向反应谱可通过以下结构措施实现：增加楼板厚度、增设梁构件、在楼板下方布置柱或剪力墙以调整其竖向振动特性。在不影响下层房间使用功能前提下，在1点增加X、Y方向梁（300×800）及柱（400×400），高频竖向反应谱值明显降低，14 Hz附近峰值降低为5.9 m/s²［图9（b）］。

地下楼层为钍基熔盐堆核岛设备集中放置区域，减小楼层谱的放大效应对设备抗震分析尤为重要。地下水平向楼层谱谱型相似，由于处于地下，对地震动的放大效应较小。地下竖向楼板反应谱由于楼板局部振动特性差异，楼板不同位置的竖向响应差异较大，通过对结构局部布置调整改变楼板局部振动特性，可以改善该处竖向反应谱。

6 结 论

本文采用混合模拟方法拟合生成地表人工加速度时程；利用EERA软件进行土层反应计算，得到自由场不同场地深度的地震反应；详细给出了水平、竖向地震作用下人工边界5个面等效应力大小和方向；通过对钍基熔盐研究堆及土体的建模，进行了SSI效应下的地震响应分析，并对地下楼层竖向反应谱进行了初步优化，得到以下结论：

（1） 本文的混合模拟方法采用“频域预调+时域准确匹配”双阶段迭代策略，通过频域初筛减少无效迭代次数，再利用时域叠加窄带时程提高匹配精度，高效生成了满足规范要求的高TDRS匹配精度且包络TPSD的人工地震动时程，为钍基熔盐堆抗震分析提供了基准输入。

（2） 相对于水平方向，地下部分竖直方向楼层反应谱由于楼板局部振动特性不同，楼板不同位置的反应谱差异较大，板中位置受楼板局部振动影响，会在楼板局部振动固有频率附近出现峰值，对核设备抗震分析带来不利影响。

（3） 通过结构局部调整（如增加板厚、增设梁柱墙等），增加竖向刚度，改变楼板局部振动特性，可削减板中竖向反应谱峰值，是改善该部分楼板竖向反应谱的有效手段。本文通过板跨中增设梁柱，跨中点竖向反应谱减小为原值的0.53倍，大幅度优化了竖向反应谱。

参考文献

[1] 中国地震局.核电厂抗震设计标准：GB 50267—2019［S］.北京：中国计划出版社，2019.

[2] China Earthquake Administration.Standard for seismic design of nuclear power plants：GB 50267—2019［S］.Beijin：China Planning Press，2019.（in Chinese）

[3] LYSMER J，OSTADAN F，CHIN C.A system for analysis of soil-structure interaction［D］.Berkeley，California：University of California，1999.

[4] 刘艺诚，王晓，王晓艳，土壤-结构相互作用的TMSR-LF1厂房楼层反应谱分析［J］.核技术，2021，44（7）：070602.

[5] LIU Yicheng，WANG Xiao，WANG Xiaoyan，et al.TMSR-LF1 plant floor response spectrum analysis of soil-structure interaction［J］.Nuclear Techniques，2021，44（7）：070602.（in Chinese）

[6] WANG X X，ZHOU Q，ZHU K X，et al.Analysis of seismic soil-structure interaction for a nuclear power plant （HTR-10）［J］.Science and Technology of Nuclear Installations，2017，2017（1）：2358403.

[7] NIE J，GRAIZER V，SEBER D.A greedy algorithm for wavelet-based time domain response spectrum matching［J］.Nuclear Engineering and Design，2023，410（8）：1.1-1.17.

[8] ZHANG Yushan，ZHAO Fengxin.Artificial ground motion compatible with specified peak ground displacement and target multi-damping response spectra［J］.Nuclear Engineering and Design，2010，240（10）：2571-2578.

[9] 何佳，王海涛.用于窄带叠加人工时程模拟的一种窄带构造算法［J］.核动力工程，2012，33（3）：69-73.

[10] HE Jia，WANG Haitao.A narrowband construction algorithm for narrowband overlay artificial time history simulation［J］.Nuclear Power Engineering，2012，33（3）：69-73.（in Chinese）

[11] 王海涛，何树延，张征明，同时满足反应谱匹配和功率谱密度包络要求的人工时程生成算法［J］.核动力工程，2009，30（1）：112-115.

[12] WANG Haitao，HE Shuyan，ZHANG Zhengming，et al.Artificial time history generation algorithm that simultaneously meets the requirements of response spectrum matching and power spectral density envelope［J］.Nuclear Power Engineering，2009，30（1）：112-115.（in Chinese）

[13] CHEN X M，LÜ F H，ZHANG J D，et al.Artificial seismic waves generation for complex matching conditions based on diffusion model［J］.Soil Dynamics and Earthquake Engineering，2025，192（5）：109290.

[14] BESKOS D E.Boundary element methods in dynamic analysis：Part II （1986-1996）［J］.Applied Mechanics Reviews，1997，50（3）：149-197.

[15] 王志亮，阳栋.地震载荷下核岛动力响应数值模拟［J］.岩土力学，2013，34（S2）：400-406.

[16] WANG Zhiliang，YANG Dong.Numerical simulation of dynamic response of nuclear island under earthquake load［J］.Geotechnical Mechanics，2013，34（S2）：400-406.（in Chinese）

[17] 廖振鹏.工程波动理论导论［M］.2版.北京：科学出版社，2002.

[18] LIAO Zhenpeng.Introduction to engineering wave theory［M］.Second Edition.Beijing：Science Press，2002.（in Chinese）

[19] LYSMER J，KUHLEMEYER R L.Finite dynamic model for infinite media［J］.Journal of the Engineering Mechanics Division，1969，95（4）：859-878.

[20] 刘晶波，吕彦东.结构-地基动力相互作用问题分析的一种直接方法［J］.土木工程学报，1998，31（3）：55-64.

[21] LIU Jingbo，LÜ Yandong.A direct method for analyzing the dynamic interaction between structure and foundation［J］.Journal of Civil Engineering，1998，31（3）：55-64.（in Chinese）

[22] 谷音，刘晶波，杜义欣.三维一致黏弹性人工边界及等效黏弹性边界单元［J］.工程力学，2007，24（12）：31-37.

[23] GU Yin，LIU Jingbo，DU Yixin.Three dimensional consistent viscoelastic artificial boundary and equivalent viscoelastic boundary element［J］.Engineering Mechanics，2007，24（12）：31-37.（in Chinese）

[24] 隋翊.超设计基准荷载作用下核岛厂房动力响应特性的精细化研究［D］.大连：大连理工大学，2022.

[25] SUI Yi.Refined study on dynamic response characteristics of nuclear island building under beyond design benchmark load［D］.Dalian：Dalian University of Technology，2022.（in Chinese）

[26] WOLF J P，SOMAINI D R.Approximate dynamic model of embedded foundation in time domain［J］.Earthquake Engineering & Structural Dynamics，1986，14（5）：683-703.

[27] LYSMER J，TABATABAIE-RAISSI M，TAJIRIAN F，et al.SASSI：A system for analysis of soil-structure interaction［R］.Berkeley：University of California，Berkeley，1981.

[28] WANG W D，LI H W，WANG J X.Progressive collapse analysis of concrete-filled steel tubular column to steel beam connections using multi-scale model［J］.Structures，2017，9（2）：123-133.

[29] 石永久，王萌，王元清.基于多尺度模型的钢框架抗震性能分析［J］.工程力学，2011，28（12）：20-26.

[30] SHI Yongjiu，WANG Meng，WANG Yuanqing.Seismic performance analysis of steel frames based on multi-scale models［J］.Engineering Mechanics，2011，28 （12）：20-26.（in Chinese）

[31] 侯春林，李小军，潘蓉.不同法规要求的核电厂设计地震动时程［J］.核安全，2012（4）：46-50.

[32] HOU Chunlin，LI Xiaojun，PAN Rong.Design seismic time history of nuclear power plants under different regulatory requirements［J］.Nuclear Safety，2012（4）：46-50.（in Chinese）

[33] 中国地震局地质研究所.小型模块化钍基熔盐堆研究设施项目地震安全性评价专题报告［R］.武汉：中国地震局地质研究所，2022.10.

[34] Geological Research Institute of China Earthquake Administration.Special report on seismic safety evaluation of small modular thorium molten salt reactor research facility project［R］.Wuhan：Geological Research Institute of China Earthquake Administration，2022.10.（in Chinese）

[35] BARDET J P，ICHII K，LIN C H.EERA—A computer program for Equivalent-linear Earthquake site Response Analyses of layered soil deposits［R］.2000.

[36] FEMA.FEMA356，Prestandard and Commentary for the Seismic Rehabilitation of Buildings［S］.Washington，D.C：Federal Emergency Management Agency，2000.

[37] 张洪才.ANSYS 14.0理论解析与工程应用实例［M］.北京：机械工业出版社，2013.

[38] ZHANG Hongcai.ANSYS 14.0 theoretical analysis and engineering application examples［M］.Beijing：Machinery Industry Press，2013.（in Chinese）

[39] 朱凯昕，王晓欣，史力，三向地震载荷下HTR-10厂房土壤-结构相互作用分析［J］.原子能科学技术，2018，52（10）：1817-1824.

[40] ZHU Kaixin，WANG Xiaoxin，SHI Li，et al.Analysis of soil-structure interaction in HTR-10 plant under three directional seismic load［J］.Atomic Energy Science and Technology，2018，52（10）：1817-1824.（in Chinese）

基金资助

国家重点研发计划项目(2023YFC3806100)

上海建工集团基础研究项目(24JCYJ-05)