三维隔震单层球面网壳结构频率相关变阻尼半主动控制

2022-03-08石运东

地震工程与工程振动 2022年1期

石运东，韦鑫，丁阳，王涛

（1.中国地震局工程力学研究所，中国地震局地震工程与工程振动重点实验室，黑龙江哈尔滨 150080；2.天津大学建筑工程学院，天津 300350）

引言

大跨空间结构由于跨度较大，在三维方向周期均较长，地震作用下其竖向响应与水平向响应往往是同量级的［1-4］。近年来学者们研发提出的三维隔震技术，可同时降低大跨空间结构的三维动力响应［4-10］。然而，由于地震动的不确定性，长周期地震时地震动的低频成分可能接近大跨空间结构的隔震频率并引发共振现象，降低三维隔震效果。如陈兆涛等［9-10］在开展大跨空间结构三维隔震支座振动台试验时，地震动竖向分量主频接近竖向隔震频率引起共振，导致响应放大。

为降低大跨空间隔震结构在共振时的动力响应，可在隔震层增设阻尼器。然而，目前常用的被动形式阻尼器设计参数固定，无法确保在共振与非共振等不同情况时的最优控制［10］。半主动控制结合了主动控制和被动控制的特点，能根据外界荷载和结构自身的振动反应实时调整控制策略，在大跨空间结构研究中受到较多关注。如倪莉等［11］结合网壳结构高次超静定的特点提出了将半主动控制器做成变刚度变阻尼可控制杆件代替网壳结构中某些杆件对网壳结构进行减震控制的策略。瞿伟廉等［12］对装有ER/MR 智能阻尼器的空间网壳结构的地震反应进行了计算。李忠献等［13］通过磁流变阻尼器对大跨平板网架结构的多维地震反应进行震动控制。张路［14］研究了体育馆三角桁架屋盖结构利用磁流变阻尼器作为立柱与桁架斜撑时的减震效果等。以上研究多通过在非隔震结构中使用智能阻尼器替换结构杆件或者添加支撑的方法进行时域内的半主动控制。相比于非隔震结构，隔震结构的响应主要由隔震层特性决定，且受到地震动频谱特征的影响显著。既有研究结果表明，对于隔震结构，地震动频谱特征不同时隔震层最优设计参数不同［15-16］。非隔震结构中基于时域的控制方法在地震动频谱特征变化时需要间接根据结构响应的大小调整控制力，较难实现大跨空间隔震结构在不同频率激励时的优化控制。

文中以K8型单层球面网壳结构为对象，从频域角度提出一种地震动主频相关加速度阈值变阻尼半主动控制方法，研究其在三维隔震大跨空间结构中的控制效果。基于ABAQUS 软件的DLOAD 子程序，对三维隔震单层球面网壳结构进行半主动控制模拟与时程分析，验证半主动控制算法的控制效果。

1 隔震层阻尼与地震动主频对结构响应影响

1.1 单自由度体系阻尼对传递率影响规律

对三维隔震结构在单方向以单自由度结构进行简化分析。单自由度体系受迫振动时，定义结构绝对加速度响应幅值与输入绝对加速度幅值的比值为加速度传递率TRa［17］：

定义结构相对位移响应幅值与输入绝对加速度幅值的比值为位移传递率TRd［17］：

式中：f为激励频率；fn为结构固有频率；ζ为结构阻尼比。图1给出了加速度传递率TRa以及位移传递率TRd曲线。由图可知，当f/fn＜1.414时，阻尼比越大，加速度与位移传递率越低，定义f/fn比值区间（0，1.414）为共振区。当f/fn＞1.414 时，阻尼比越小，加速度传递率越低，定义f/fn比值区间（1.414，+∞）为非共振区。为降低结构加速度，在共振区应选取较大阻尼比，而在非共振区应选取较小阻尼比。在非共振区，虽然阻尼比越小，位移传递率越高，但此时的高频地震动引起的位移值较小，不会引起隔震支座安全问题。

图1 不同阻尼比时单自由度体系响应传递率Fig.1 Response transmissibility for SDOF system with different damping ratio

1.2 基于地震动主频的最优阻尼比

实际地震动含有较宽频带的频率成分，难以直接确定其对应的优化阻尼比。文中通过对地震动时程进行快速傅里叶变换（FFT），计算频谱图中峰值处所对应的频率作为地震动主频。并定义，当地震动主频与结构频率比处于共振区时，地震动称为长周期地震动；而处于非共振区时称为短周期地震动。

图2 给出了来自于同一台站El Centro Array #9 所记录的两条地震动Mtn（Borrego Mtn 地震，1940 年）与Imperial-01（Imperial Valley-01地震，1938年）水平方向分量的加速度时程与FFT 频谱，两条地震动均调幅至9度罕遇地震620 Gal。由图中可知，Mtn主频为0.60 Hz，Imperial-01主频为6.70 Hz。假设某单自由度体系自振频率fn=0.74 Hz，则根据前述定义，Mtn属于长周期地震动，Imperial-01属于短周期地震动。

图2 加速度时程与FFT频谱Fig.2 Acceleration time history and FFT spectrum

对此单自由度结构进行动力时程分析，阻尼比ζ分别取0.1 与0.707。结果表明，阻尼比0.1 与0.707 时，在Mtn 地震作用下，加速度分别为9.61 m/s2与4.42 m/s2，位移分别为0.44 m 与0.14 m，大阻尼比效果更好；在Imperial-01 地震作用下，加速度分别为1.02 m/s2与1.68 m/s2，位移分别为0.05 m 与0.02 m，小阻尼比效果更好。利用地震动主频作为指标判断优化阻尼比，时程分析结果与传递函数分析结果吻合。

2 地震动主频相关加速度阈值变阻尼半主动控制

2.1 地震动主频相关两级变阻尼半主动控制

为了实现实时最优阻尼控制，文中提出一种地震动主频相关两级变阻尼半主动控制方法。两级阻尼的确定主要考虑阻尼设备所能提供阻尼能力及结构振动控制的目标要求。半主动控制基本准则如式（3）所示：

式中：ζ为半主动控制实时阻尼比；ζhigh为共振区使用较高阻尼比；ζlow为共振区外使用较低阻尼比。

2.2 加速度阈值变阻尼优化

两级变阻尼半主动控制基本准则存在一个缺陷：在大小阻尼比切换时，阻尼力突变。此现象的根本原因在于两级可变阻尼差别较大。为弥补这种算法的缺陷，提出一种加速度阈值阻尼优化方法。

考虑一单自由度隔震结构，其阻尼力FD计算公式如式（4）：

式中：m为结构质量；ωn为隔震结构圆频率；ζ为隔震层阻尼比；v为隔震结构相对地面速度。定义阻尼力突变产生的结构加速度为附加加速度，计算如式（5）：

式中：a为附加加速度；FD（t）为当前t时刻的阻尼力；FD（t+Δt）为t+Δt时刻阻尼力。当Δt足够小时，t时刻与t+Δt时刻结构速度接近，假设其均为v（t）。若可容许附加加速度为amax，结合式（4）以及式（5）得t+Δt时刻阻尼比计算式（6）：

式中：ζ（t）为当前t时刻的隔震层阻尼比；ζ（t+Δt）为t+Δt时刻隔震层阻尼比。式中，加减号根据频率比所处区段决定。若t时刻频率比处于共振区内，t+Δt时刻频率比处于非共振区内，则隔震层阻尼比应由大阻尼状态调整为小阻尼状态，取减号；反之，隔震层阻尼比由小阻尼状态调整为大阻尼状态取加号。

2.3 长短周期地震动分区识别

文中地震动长短周期的定义由地震动主频与隔震结构主频相对关系而定。第1.2节进行地震动主频计算时，实际是对已知完整地震动数据进行整体分析。而实测时，仅能针对实时所测部分数据进行逐步分析。文中使用基于快速傅里叶变换（FFT）的主频检测方法。

该方法采集结构附近地面的地震动加速度，对于采样时长Δt内的数据进行FFT分析，获取所采样数据的主频信息。而在下一时刻（+1/Samp，Samp 为采样频率），FFT 所分析数据进行一个采样数据的更新，继续计算新的主频信息，以此形成主频信息的时程。从该方法可以看出，实际所测地震动主频存在时间滞后效应，即由采集数据获得的主频f不能真实反应实时的频率，而是检测时间点之前Δt时间内的历史数据的频率。当地震动频率特征随时间变化较快时，所测主频可能出现较大误差。因此，为提高信号的实时性，Δt不能过长。式（7）给出了采样时长Δt与所测频率分辨率Δ f的相互关系：

由式（7）可知，频率分辨率Δ f与采样时长Δt成反比。例如，采样时间Δt=0.5 s，则频率分辨率Δ f为2 Hz，即所测得的地震动主频f只能反馈得到结果0、2、4 Hz 等频率。当需要获得准确频率信息时，应尽量增加采样时长Δt［15］，这与前述提高信号实时性的要求相悖。对于文中所提出的两级变阻尼算法，实际无需对于地震动主频的精确检测，而只需准确识别地震动主频与长短周期地震动分界线（f=1.414fn）的相互关系。如图3 所示，对于所能识别出的第i个频率值为iΔf，根据实际应用中四舍五入的计算方法，其对应的地震动真实频率为区间［（i-0.5）Δ f，（i+0.5）Δ f］。例如，图中所示某时刻地震动数据的真实主频为f1，假设按照前述定义其属于长周期地震动范围（f1＜1.414fn）。对其无需精确识别为f1，而只要其所处iΔf对应区间位于长短周期分界线左侧即可，即满足：

图3 快速傅里叶变换分辨率设置Fig.3 Setting for resolution of FFT analysis

虽然无需对主频进行精确检测，但过短的采样时长Δt可能无法有效检测长周期成分。图4 给出了前述Mtn 地震动与Imperial-03 地震动（Imperial Valley-03 地震，1951年，El Centro Array#9台站）X方向分量的主频时程，图中虚线给出了整条地震动时程分析所得主频。对于长周期地震动Mtn，在Δt=0.48 s（i=0，fn=0.74 Hz）时，由于采样时长过短，导致地震动主频检测结果受到短时频率影响过大，未能准确地捕捉地震动的低频信息，其对于地震动加速度时程中（如t=10 s左右）长周期频率成分的控制可能出现错误。对于短周期地震动Imperial-03，在Δt=10.02 s（i=10，fn=0.74 Hz）时，由于采样时长过长，历史信息影响实时频率的判断。如在时程后半部分受到历史数据影响，将应该是低频的部分误判为高频。

图4 不同FFT分析设置时地震动主频测试结果Fig.4 Dominant frequency of ground motions with different settings of FFT analysis

不同于常用时域内的半主动控制方法，地震动主频相关加速度阈值变阻尼半主动控制通过频域的快速傅里叶变换分析地震动激励的频率特征，然后依据地震动激励的频率与结构自身频率的相互关系，“实时”调整半主动控制所提供给隔震结构的阻尼，从而达到在不同频率地震动激励时的优化控制。

3 基于结构有限元模型的半主动控制模拟方法

前述所提出的半主动控制方法是基于单自由度结构。为验证所提出控制算法在三维隔震单层球面网壳结构中的应用效果，有必要进行整体结构的响应分析。既有半主动控制算法研究通常借助MATLAB、MATLAB/Simulink 等软件完成。但由于三维隔震单层球面网壳结构模型建模复杂，在MATLAB 等软件中不易实现。本节结合ABAQUS 有限元软件，提出一种基于结构有限元模型的半主动控制模拟方法。实际工程结构中结构响应无需利用有限元计算，可直接通过传感器反馈获得。

ABAQUS 软件的Connector 单元一般可用来模拟特定构件（如隔震与阻尼元件）的刚度与阻尼。文中通过取消激活Connector 单元阻尼参数，将Connector 单元的常阻尼属性替换成变阻尼属性，并通过用户自定义子程序DLOAD施加荷载的方法实现对结构半主动控制阻尼力的施加。

图5给出了DLOAD模拟方法工作流程图。该方法借助子程序URDFIL实现对于结构模型参数的采集与输出。DLOAD 通过对于输出的结构模型参数的读取调用，实现控制力的计算。在对结构模型施加控制力时，DLOAD 子程序需要确定荷载施加位置与施加方向。为实现三维隔震单层球面网壳结构多点阻尼力控制，在共享的子程序中根据不同荷载点位的三维坐标以条件语句的形式计算N个不同点位的阻尼力输出F1…FN。在ABAQUS计算的每一个荷载步，不同荷载点位以该位置的坐标信息，自动寻找子程序中相应的阻尼力值规定并进行输出。另外，对于同一点位的多个方向，控制力各不相同，而DLOAD子程序对于同一点位只能输出一个相同幅值的力（方向可不同），因此需在理论荷载施加位置相近的点位构造多个荷载施加位置（尽量靠近理论荷载位置，减小引起的附加弯矩等），从而实现多点多维控制［18］。

图5 DLOAD模拟方法工作流程图Fig.5 Flowchart of the DLOAD simulation method

4 半主动控制效果分析

4.1 有限元模型及地震动

对一个90 m跨度的K8型单层球面网壳进行时程分析，以验证提出的频率相关加速度阈值变阻尼控制对三维隔震单层球面网壳结构动力响应的控制效果。单层球面网壳的矢跨比为1/4，屋面荷载取1.5 kN/m2，结构参数如表1所示，计算模型如图6所示。在有限元软件ABAQUS中，建立无隔震与三维隔震的单层球面网壳模型。上部结构杆件采用B31梁单元模拟，材料本构采用理想弹塑性模型，所有杆件屈服强度为235 MPa。分析时上部结构采用Rayleigh 阻尼，阻尼比取0.03。经过分析，无隔震时上部结构水平与竖向基本周期分别为0.30 s和0.25 s。三维隔震水平与竖向隔震设计周期分别取上部结构水平与竖向基本周期的4倍。支座数目共32个，阻尼器共计96个，隔震支座和阻尼器均采用Connector单元模拟。

表1 单层球面网壳结构参数Table 1 Parameters of single-layer spherical reticulated shell

图6 单层球面网壳计算模型Fig.6 Single-layer reticulated spherical shell model

选取了12 条具有不同频率特征的地震动进行有限元计算，信息如表2所示。根据地震动主频以及结构频率，选择X向地震动分量进行长短周期区分。根据对长、短周期地震动定义，表中1～6号地震记录为长周期地震动，而7～12号地震记录为短周期地震动。将X向最大加速度调幅至620 Gal，Y、Z方向最大加速度按照X、Y、Z这3个方向1∶0.85∶0.65比例调幅，以计算单层球面网壳在9度罕遇地震作用下的响应。

表2 地震记录Table 2 Earthquake records

通过模态分析，三维隔震结构在X、Y、Z这3 个方向的主频fn分别为0.74 Hz（1.35 s）、0.74 Hz（1.35 s）、0.92 Hz（1.09 s）。根据式（8），水平X、Y方向选择i=2，竖直Z方向选择i=3，计算可得X、Y、Z这3 个方向通过FFT分析地震动主频的频率分辨率Δ f分别为0.42、0.42、0.37 Hz，所对应采样时长Δt分别为2.39、2.39、2.69 s。定义工况C1、C7、NODO、ACCTDO 分别表示隔震层3 个方向采用较小阻尼比0.1 的被动控制、采用较大阻尼比0.707的被动控制、采用频率相关但无加速度阈值优化的变阻尼半主动控制（no damping optimization，NODO）以及采用频率相关加速度阈值变阻尼半主动控制（acceleration threshold damping optimization，ACCTDO）。ACCTDO计算时阻尼力变化引起的容许附加加速度取amax=0.05 m/s2。

文中对于结构的加速度与轴力2个指标进行评估。定义指标J1为结构顶点加速度，J2为结构杆件最大轴力。为了评估控制效果，将隔震层三向阻尼比采用较小值0.1 以及较大值0.707 的2 个计算结果中响应较小的值称为最优被动控制结果JOPC，并将以上4 类工况的响应与最优被动控制结果进行对比，以响应放大系数η（η=J/JOPC）进行评估。

4.2 DLOAD控制模拟方法准确性验证

使用DLOAD 模拟方法对粘滞阻尼器提供的阻尼力进行等效替换模拟，与常规使用Connector 单元模拟阻尼力的结果进行对比。选取表2 中Loma 与Hollister 地震动，将原波输入至网壳结构进行时程分析。图7给出了Loma地震时两种工况的结构顶点加速度时程。由图可知，两者的模拟结果吻合度高。以Connector 模拟方法结果为标准，在Loma地震作用下，DLOAD 方法三向加速度峰值模拟误差分别为1.9%、1.9%、3.7%，隔震层位移峰值模拟误差分别为0%、0%、0.1%。在Hollister 地震作用下，DLOAD 三向加速度峰值模拟误差分别为4.1%、3.5%、5.3%。隔震层位移峰值模拟误差分别为0.2%、0%、0.8%。由结果可知，DLOAD 模拟方法与Connector模拟方法的结果拟合良好，可用于三维隔震单层球面网壳结构的半主动控制模拟。此方法也适用于其它类型结构的半主动以及主动控制研究。

图7 Loma地震时结构顶点加速度对比Fig.7 Vertex acceleration under Loma motion

4.3 节点加速度控制效果分析

表3 给出了不同地震动作用下单层球面网壳顶点的最大加速度，图8 为4 种工况X向与Z向加速度响应放大系数对比。其中，前6条地震动为长周期，后6条地震动为短周期。根据模拟结果，在被动控制时，恒定阻尼比难以在不同周期时保持最优的控制效果。长周期地震动时，C1 工况（阻尼比0.1）结构三向加速度放大系数（J1/J1OPC）平均值分别为1.93，1.76，1.51；短周期地震动时，C7工况（阻尼比0.707）结构三向加速度放大系数平均值分别为1.42，1.45，1.69。以上2 类情况3 个方向的加速度放大系数平均值为1.63。由表3 与图8可以看出，频率相关加速度阈值变阻尼半主动控制ACCTDO 工况时，结构加速度响应基本与最优被动控制结果接近。根据计算结果，ACCTDO 工况时所有地震动作用下结构不同方向的加速度响应放大系数平均值为1.01，可以显著提高不同频率地震动时的控制效果。

图8 顶点加速度放大系数Fig.8 Amplification coefficient of vertex acceleration

表3 顶点最大加速度响应Table 3 Peak acceleration of the vertex

图9 给出了Imperial-03 地震动时的顶点Z方向的加速度响应。该地震动在开始阶段频率较高，而在后期阶段频率较低。采用被动控制时，开始阶段C1 工况加速度低于C7 工况加速度，后期阶段C7 工况加速度低于C1 工况加速度。而使用ACCTDO 进行控制时，阻尼自动调整优化，使得开始高频振动阶段的响应接近于C1工况结果，而后期低频振动阶段的响应接近于C7工况结果。因此，所提出的频率相关加速度阈值变阻尼半主动控制方法可以有效分析地震动频率，并及时做出阻尼的优化调整，可以获得与预先已知地震动而进行最优被动控制相近甚至更低的结果。适应于不同频率特性的地震动作用下的结构震动控制。

图9 Imperial-03波不同工况顶点Z向加速度时程Fig.9 Time history of Z acceleration of Imperial-03 under different conditions

在部分地震动，特别是长周期地震动时，如果直接采用高低阻尼切换而不进行加速度阈值优化（NODO工况），结构的加速度会有显著放大。如Loma 地震时，使用加速度阈值阻尼优化时X向加速度为4.66 m/s2，而不进行优化时X向加速度为12.16 m/s2，相对提高161%。图10给出了Loma地震时顶点X方向加速度时程以及半主动控制阻尼时程。图中存在较明显阻尼差异。如t=5.19 s时，根据识别的地震动主频，阻尼需由0.1切换为0.707。不进行阻尼优化时阻尼比突然增加，导致了图10（a）的加速度突增现象。而根据加速度阈值优化时，阻尼比则缓慢地变化，避免了结构受力以及加速度的突然变化。因此，文中提出的加速度阈值阻尼优化方法可以有效避免两级阻尼高低切换导致的阻尼力突变以及过大的附加加速度。

图10 Loma地震时顶点X方向加速度时程及半主动控制阻尼比时程曲线Fig.10 Time history of X acceleration of vertex and controlled damping ratio under Loma

4.4 杆件内力控制效果分析

表4给出了不同地震动作用下单层球面网壳杆件的最大轴力（包括静力作用下轴力）。图11为4种工况时环杆与斜杆轴力响应放大系数对比。与加速度计算结果类似，在被动控制时，恒定阻尼比难以在不同频率特征地震动时保持最优的轴力控制效果。长周期地震动时，C1工况（阻尼比0.1）结构环杆、肋杆以及斜杆轴力放大系数（J2/J2OPC）平均值分别为1.44，1.73，1.67；短周期地震动时，C7 工况（阻尼比0.707）结构环杆、肋杆以及斜杆轴力放大系数平均值分别为1.70，1.68，1.64。以上2 类情况3 个方向的轴力放大系数平均值为1.64。由图11可以看出，频率相关加速度阈值变阻尼半主动控制ACCTDO 工况时，环杆、肋杆以及斜杆轴力响应基本与最优被动控制结果接近。根据计算结果，ACCTDO工况时所有地震动作用下环杆、肋杆以及斜杆轴力放大系数平均值为1.01。采用基于加速度阈值的变阻尼控制方法对于结构轴力的控制改善效果与加速度控制改善效果接近。