非能动安全壳冷却系统(passive containment cooling system, PCS)是第三代核电厂区别于传统核电厂的显著特点。由于液体的晃动频率远低于结构的自振频率, 地震动中长周期分量对冷却水箱的影响无法忽视。此外, 高置冷却水箱的设置对地震作用下核电厂辅助厂房关键设备层楼层谱的影响, 也要重点考虑。因此, 需要研究水箱结构的动力特性, 从而更准确地进行核岛结构的动力反应分析。

在核电厂非能动冷却水箱动力特性研究方面, 赵春风等^[1-3]考虑流固耦合作用, 研究了AP1000核电厂屏蔽厂房重力水箱在地震动作用下的动力响应和重力水箱的减震效果。党俊杰等^[4]根据AP1000非能动冷却水箱及屏蔽厂房的固有频率建立与其固有频率近似的圆柱水箱试验模型, 测得试验模型的固有频率和水的晃动频率。曾晓佳等^[5]通过1/40的模型试验, 测量不同液深时顶盖位置的冲击力, 验证了液体晃动冲击力计算公式的可行性。

在水箱动力特性试验研究方面, 目前的试验手段一般是通过激光位移传感器进行液体晃动位移的测量, 得到其晃动频率。王为等^[6]采用位移传感器测得了液体波高的时程曲线, 通过对数据的处理, 得到晃动频率和阻尼比。王剑楠等^[7]研制了船舶液舱晃荡模拟系统, 通过激光位移传感器来监测晃荡过程液体运动情况, 记录液体晃动波高。此外, 李松等^[8]通过激光扫描采集液面的速度响应, 对速度信号进行FFT分析, 得到水平放置的半充液圆柱容器的液体晃动低阶频率值。在阻尼比计算方面, 只有一些几何形状比较规则的容器才有解析解或近似解, 一般采用试验方法。

晃动基频和阻尼比是液体晃动动力特性的关键参数, 也是基于Housner模型进行液体晃动分析模型简化处理的关键参数。相关规范中^[9-10]基于相关学者的研究给出了规则圆柱水箱液体晃动频率的简化计算公式, 但对于超规范非规则环形圆柱水箱, 并没有进行详细介绍。晃动阻尼比方面, 虽然规范也根据相关研究^[11]给出了建议值, 但是对于特殊结构, 也缺乏相关的研究和验证。

本文建立水箱结构缩尺模型, 通过振动台试验测试水箱结构内壁不同高度处动水压和波高衰减数据, 计算获得液体的晃动频率和阻尼比, 并与理论计算结果进行对比, 验证方法的合理性。

1 模型振动台试验测试 1.1 模型介绍

根据实际结构参数, 按照几何相似比1/16制作试验模型。为了满足水箱结构刚性假定和波高测试需要, 水箱侧壁和底板采用微粒混凝土, 用镀锌铁丝模拟钢筋, 水箱顶板采用钢化玻璃。试验模型剖面及各部分尺寸如图 1所示, 模型结构各部分材料参数如表 1所示。

针对液体晃动, 现有试验研究采用的测试设备一般是激光位移传感器或扫描式激光测振仪。一方面, 高精度的激光位移传感器或相关采集系统成本高；另一方面, 为满足激光测试要求, 需要在液面使用有机溶剂(香蕉水)溶解有色物质来增加激光反射, 该方法有一定的危险性, 并且该方法在测试波高时对飞溅液体无法捕捉, 液体表面的有机溶剂会在一定程度影响液体自身的晃动特性。因此, 本文试验采用水压力传感器进行动水压测量, 并进行晃动频率的识别。

在水箱结构内壁X截面、Y截面分别沿高度布置两组共16个水压力传感器。同时, 在X截面和Y截面水箱内壁布置刻度尺, 钢化玻璃顶板相应位置布置摄像头, 记录液体自由晃动段的波高衰减数据, 传感器布置情况见图 2。

1.2 试验过程

通过振动台施加两组幅值分别为0.05 g和0.10 g的单向、三向谐波, 记录关键位置的动水压数据和波高衰减情况, 进而求得液体的晃动频率和阻尼比。试验中输入的归一化谐波信号如图 3所示。

1.3 试验数据分析 1.3.1 晃动频率测试与识别

利用动水压数据计算液体的晃动频率基于数据信号的Fourier变换。Fourier变换的典型用途是将信号分解成振幅分量和频率分量：

$ F\left( \omega \right) = \int_{ - \infty }^{ + \infty } {f\left( t \right){{\rm{e}}^{ - {\rm{i}}\omega t}}{\rm{d}}t} $

(1)

式中：f(t)是随时间t变化的数据信号, F(ω)是f(t)的频谱函数, |F(ω)|是f(t)的振幅频谱。

试验过程中, 在外部激励结束后, 液体仍会继续保持晃动并逐步衰减, 此时自由晃动时的动水压参数只反应水箱中液体的固有属性, 与外荷载无关。居荣初等给出了如图 4所示的简单圆柱形储液罐侧壁动水压力的计算公式^[12]：

$ \begin{array}{*{20}{c}} {p\left( {a,\theta ,z,t} \right) = - \rho a\sin \theta \cdot }\\ {\left[ {\sum\limits_{n = 1}^\infty {{{\ddot q}_n}\left( t \right)\frac{{{d_n}\cosh \left( {{\sigma _n}\frac{{z + h}}{a}} \right)}}{{\cosh \left( {{\sigma _n}\frac{h}{a}} \right)}}} + {{\ddot x}_0}\left( t \right)} \right]} \end{array} $

(2)

$ {d_n} = \frac{2}{{\sigma _n^2 - 1}} $

(3)

式中：ρ为液体密度, a为储液罐半径, θ为圆心与侧壁某点连线与Y轴夹角, h为液体高度, $\ddot x$₀(t)为外荷载加速度, σ_n为J′₁(σ_n)=0的根, J′₁为一阶贝塞尔函数的导数, q(t)为液体振动的广义坐标, 满足下列方程：

$ {{\ddot q}_n}\left( t \right) + v{{\dot q}_n}\left( t \right) + \omega _n^2q\left( t \right) = - {{\ddot x}_0}\left( t \right) $

(4)

式中：ν为流体阻尼系数(量纲为1/s)。方程4与单自由度质点的强迫振动方程具有相同的形式, 和单自由度系统一样, ω_n可视为液体的晃动频率。液体自由振动阶段, 储液罐侧壁动水压力为

$ \begin{array}{*{20}{c}} {p\left( {a,\theta ,z,t} \right) = - \rho a\sin \theta \cdot }\\ {\sum\limits_{n = 1}^\infty {{{\ddot q}_n}\left( t \right)\frac{{{d_n}\cosh \left( {{\sigma _n}\frac{{z + h}}{a}} \right)}}{{\cosh \left( {{\sigma _n}\frac{h}{a}} \right)}}} } \end{array} $

(5)

从方程5可以看出, 侧壁动水压力可以用液体振动的广义加速度表示, 类比由加速度Fourier变换得到自振频率的方法, 通过处理动水压数据, 求得液体晃动频率。取自由晃动段的动水压数据, 如图 5(a)所示, 并对该自由晃动段数据进行Fourier变换, 如图 5(b)所示, 进而识别液体各阶晃动频率。

通过动水压数据变换计算得到的液体前4阶晃动频率统计结果见表 2。通过表中数据可以看出, 试验测得的前四阶晃动频率变异系数均小于5%, 测试结果离散性小, 实验数据合理、可靠。

1.3.2 阻尼比识别

液体的晃动阻尼主要来自于水箱与液体之间的相互摩擦和液体之间的粘滞作用, 数值大小通常与四个因素有关^[13-14]：箱体壁面黏性阻尼作用、自由液面处的黏性耗散、液体内部黏性耗散作用和毛细作用。由于参数的复杂性, 晃动阻尼比一般通过试验方法获得。本文中, 试验水箱模型基频大于30 Hz, 认为该水箱为刚性。对于规则圆柱刚性水箱, 内部为低黏性液体以及无隔板时, 文献[10]给出晃动阻尼比建议值为0.5%。对非规则结构, 规范并没有给出建议值。

利用对数衰减法计算阻尼比。计算公式如下

$ \delta = \ln \frac{{{\mu _i}}}{{{\mu _{i + 1}}}} = \frac{{2{\rm{ \mathsf{ π} }}\xi }}{{\sqrt {1 - \xi } }} $

(6)

式中：δ表示对数衰减率, μ_i、μ_i+1分别表示任意两相邻振动的峰值, ξ表示阻尼比。

由于液体晃动阻尼很小, 自由振动衰减很慢, 为获得更高精度, 采用相隔j周的振动峰值比来计算阻尼比。由此可得小阻尼体系阻尼比的近似计算公式：

$ \xi \approx \frac{1}{{2{\rm{ \mathsf{ π} }}j}}\ln \frac{{{\mu _i}}}{{{\mu _{i + j}}}} $

(7)

选取40个视频样本, 利用对数衰减法计算阻尼比, 计算结果如图 6所示。通过多组波高衰减数据计算得到的阻尼比变异系数小于5%, 测试结果的离散性较小, 实验数据合理、可靠。试验测试得到的阻尼比均值为0.513 8%, 与规范建议值0.5%吻合良好。通过对比可以看出, 刚性水箱液体晃动阻尼比数值小, 水箱形状影响并不敏感。

2 模型晃动频率数值计算

在求解流体与结构耦合作用(FSI)方面, 有限元软件Adina具有较强的求解能力和可靠性。本文利用Adina软件建立试验水箱结构数值模型, 计算液体的晃动频率, 选用3D-solid单元模拟水箱混凝土部分, 3D-shell单元模拟钢化玻璃顶板, 3D potential-based fluid单元模拟流体。有限元模型如图 7所示。

考虑到结构对称性原因, 除去重复的晃动模态, 与试验测得的前4阶晃动频率进行比较。从表 3中可以看出, 除了第三阶晃动频率误差较大外, 第一、二、四阶晃动频率误差均小于1%, 通过试验结果与数值分析结果的比较, 验证了数值分析方法的合理性, 也验证了试验测试方法的准确性。

3 等效模型理论分析

非能动冷却水箱截面非规则, 其动力参数求解复杂, 而与其形状近似的环形圆柱水箱内液体晃动的非粘性流体动力方程有闭型解析解, 可将其等效简化为环形圆柱水箱来进行晃动频率的近似计算。简化示意图见图 8, 图中R_o是水箱外半径, R_i是水箱内半径, h是冷却水箱实际液面高度, h_a是等效简化成环形圆柱水箱后液面高度。在简化过程中, 保证水箱外半径R_o、水箱内半径R_i、水的体积V不变, 通过简化前后体积不变计算出等效液面高度h_a。

环形圆柱水箱晃动频率的理论公式如下^[15]

$ {f_i} = \frac{{\sqrt {\left( {g/{R_o}} \right){\xi _i}\tanh {k_i}} }}{{2{\rm{ \mathsf{ π} }}}} $

(8)

式中：g是重力加速度；ξ_i是方程J′₁(ξ_i)Y′₁(kξ_i)-J′₁(kξ_i)Y′₁(ξ_i)=0的第i次根, J₁、Y₁分别是第一类和第二类一阶贝塞尔函数, J′₁、Y′₁分别是J₁、Y₁的一阶导数, k=R_i/R_o；k_i=(h_a/R_o)ξ_i。

根据图 1中的模型参数, 利用式(8) 求得试验模型液体晃动一阶频率：

$ k = {R_i}/{R_o} = 0.392 $

(9)

$ {\xi _1} = 1.4710 $

(10)

$ {k_1} = \left( {{h_a}/{R_o}} \right){\xi _1} = 0.734 $

(11)

$ {f_1} = \frac{{\sqrt {\left( {g/{R_o}} \right){\xi _1}\tanh {k_1}} }}{{2{\rm{ \mathsf{ π} }}}} = 0.5072\;{\rm{Hz}} $

(12)

将简化公式计算得到的晃动基频与试验均值进行比较, 如表 4所示。理论计算结果稍大于试验均值, 误差小于3%, 验证了采用该方法简化计算液体晃动频率合理、有效。

4 结论

1) 本文提出的利用动水压自由振动段时程曲线Fourier变换识别液体晃动频率的方法合理、有效, 通过与数值分析结果的对比, 验证了该方法具有较高的准确性。

2) 本文提出的方法可以在测量动水压这一关键的液体晃动动力学参数的同时, 识别液体晃动频率, 避免液面大幅波动或飞溅对激光位移传感器测试造成的误差, 在保证基频测试准确度的情况下, 可以识别出多阶液体晃动频率。

3) 根据波高的衰减关系计算得到液体晃动阻尼比, 尽管试验模型为非规则环形圆柱水箱, 但是晃动阻尼比与相关规范建议值0.5%吻合良好, 说明参数受截面形状影响不明显。

4) 试验测试得到的晃动基频与相关规范、文献中提出的等效简化计算方法吻合良好, 可以采用该方法来简化计算。