自复位摩擦耗能支撑结构多尺度模型更新数值混合模拟方法

王 涛，李 勐，孟丽岩，许国山，刘家秀，谢婧怡

(1.黑龙江科技大学 建筑工程学院，哈尔滨 150022；2.哈尔滨工业大学 土木工程学院，哈尔滨 150090)

混合试验是一种结合物理子结构试验加载和数值模拟技术的试验手段，可以有效评价结构的地震响应，在土木工程领域得到了广泛的应用[1-4]。日本学者Hakuno[5]最早提出混合试验的概念，将结构中重要的部分作为物理结构进行试验加载，数值子结构假定有可靠的模型在计算机中模拟。近几年来，混合试验已在数值积分算法[6-7]、时滞补偿[8-10]、模型更新[11-14]等关键方面取得了很大的进展。混合试验的数值部分建模是否准确，将影响混合试验精度，模型更新是解决这一问题的有效手段。模型更新混合试验方法将模型更新技术与混合试验方法相结合，利用试验子结构观测数据在线识别试验子结构模型参数，实时更新数值子结构模型参数，以提高数值模型精度。近几年来，众多研究学者针对此方法进行了深入研究。Ou[15]对在线模型更新的实时混合仿真进行了精度分析，通过从物理子结构中提取信息来提高保真度，将模型更新算法集成到实时混合仿真中，采用两个不同复杂度的Bouc-Wen模型作为目标模型，并评价该方法的性能。陈凡[16]提出基于AUKF(adaptive unscented Kalman filter，AUKF)的框架结构离线模型更新混合试验方法，在传统UKF的基础上加入方差自适应模块，能够减轻初始参数设定对参数识别结果的影响。对以Bouc-Wen为恢复力模型的防屈曲支撑进行低周反复加载数值模拟，通过AUKF算法进行Bouc-Wen模型参数识别，验证了所提出算法的准确性。陈再现[17]提出基于均匀设计的模型更新混合模拟方法，在试验子结构的基础上，通过均匀设计的手段构造验算子结构样本空间，得到最符合试验子结构的数值模型用于模型参数更新，有效提升识别效率。

然而，模型更新混合试验的误差来源主要有两方面，其中包括物理子结构边界条件的非完整性和本构模型误差。当物理子结构边界自由度较多时，受加载设备、试验场地以及试验经费等多方面的限制，其边界条件往往难以完全实现，而边界条件的缺失势必改变结构的受力状态，影响结构性能评估的准确性。Wu[18]提出了数值混合模拟方法，并将其应用在防屈曲支撑-钢筋混凝土框架结构中，解决了边界条件对模拟结果的影响，提高了混合试验的模拟精度。此外，当原型结构由多种材料所构成，如配备有耗能减震构件的新型结构，在以往的混合试验中仅仅对主要结构进行材料本构模型参数识别更新，从而忽略其他不同材料所构成的重要构件的模型参数识别更新，这样处理势必会导致在试验过程中数值子结构的反应与原型结构的真实反应存在较大差异，利用本构参数存在较大误差的数值子结构将难以准确的模拟原型结构在地震作用下的结构响应，同时也难以为物理子结构提供准确的边界条件[19]。

鉴于此，为了解决复杂结构数值模型精度，本文提出基于容积卡尔曼滤波器(cubature Kalman filter，CKF)的多尺度模型更新数值混合模拟方法，针对二层带有自复位摩擦耗能支撑(self-centering friction energy dissipation braces，SCED)结构钢框架进行数值混合模拟，来验证所提方法的有效性。

1 CKF模型更新算法

1.1 CKF模型更新算法步骤

CKF模型更新算法需建立如下系统的状态方程和观测方程

Xk=f(Xk-1,uk-1)+Vk-1

(1)

Yk=h(Xk,uk)+Wk

(2)

式中：u为输入向量；V为过程噪声，假定过程噪声向量是均值为0的高斯白噪声，其协方差矩阵为R；W为观测噪声，假定观测噪声向量是均值为0的高斯白噪声，其协方差矩阵为Q；X为系统状态向量；Y为系统观测向量；k为算法运行步数，k=1,2,…N。

求出算法第k-1步容积点

(3)

式中：i为容积点个数，i=1,2,…2n；n为状态的维数；

Sk-1为运行算法过程中k-1步误差协方差分解的一个下三角矩阵，其表达式为

(4)

ζ为容积点集，其表达式为

(5)

(6)

式中,[e]i为第i个容积点。

(7)

(8)

(9)

(10)

重采样

(11)

(12)

(13)

(14)

(15)

(16)

式中,Kk为运行算法第k步的卡尔曼增益。

校正状态量

(17)

式中,yk为第k步的试验观测值。

校正协方差阵

(18)

CKF模型更新算法流程图如图1所示。

图1 CKF模型更新算法流程图Fig.1 CKF model update algorithm flow chart

2 物理子结构数值模型

通过第1章可以看出，建立模型更新算法的状态方程和观测方程需要借助物理子结构数值模型方程得以实现。因此，物理子结构的模型方程尤为重要，本文主要采用两种层次的本构模型来进行数值模拟，其中包括：构件层次的本构模型和材料层次的本构模型。

2.1 构件层次的本构模型

为验证多尺度模型更新数值混合模拟方法的有效性，以SCED为研究对象，采用基于改进的旗形模型(modified flag-shaped model,MFS)[20]的构件层次本构模型模拟SCED结构的滞回性能，对MFS模型进行参数在线识别，SCED构造图与MFS模型力与位移变化关系如图2所示。

图2中：P为外部荷载；T0为预应力；F为能量耗散器中的静摩擦力；Fs为内摩擦装置的激活力。

图2 自复位摩擦耗能支撑构造图Fig.2 Self-centering energy dissipation braces structure

MFS模型由线模型、双线性模型、具有可滑移弹塑性模型组成，自复位摩擦耗能支撑机理如图3所示。MFS模型初始阶段刚度K为线性模型刚度K1、双线性模型弹性刚度K2、具有可滑移弹塑性模型的第一刚度K3之和，即K=K1+K2+K3。

(a) 线性部分

其中K1、K2、K3由式(19)～(21)给出

K1=A·K

(19)

K2=(1-A)·(1-Q)·K

(20)

K3=(1-A)·Q·K

(21)

式中,A为激活刚度比，其表达式为

A=(1-Q)·b

(22)

Q为耗能率，是能量耗散的重要指标，其表达式为

Q=K3/(K2+K3)

(23)

当摩擦装置启动后，支撑刚度减小至预应力筋刚度K1，卸载初始阶段的卸载刚度等于K1+K3，当力到达激活点b时，刚度变成K1+K2。

MFS的恢复力F可表达为

F=K1x+K2R(x)+K3z

(24)

式中,R(x)为由双线性弹性模型提供的恢复力，kN；

R(x)具体表达式可由式(25)给出

(25)

式中：H(x)为海维赛德阶跃函数；b为激活位移，mm；z为模型滞回位移，mm；z微分表达式由式(26)给出

(26)

2.2 材料层次的本构模型

钢框架非线性梁柱单元采用材料层次本构模型建模，本文选择OpenSees软件中具有各向同性应变硬化的单轴Giuffre-Menegotto-Pinto钢材料模型，此模型由Filippou[21]所提出，OpenSees将此材料命名为Steel02，其应力与应变关系如图4所示。

图4 钢材应力应变关系图Fig.4 Stress-strain diagram of steel

图中：R0、R1、R2为控制弹性到塑性发展的参数；ξ、ξ1、ξ2为方程系数；b为第二刚度系数。

钢材本构模型可由如下非线性函数表示

σs=g(fy,Es,b1,R0,R1,R2,a1,a2,a3,a4)

(27)

式中：σs为钢材任意时刻的应力；g为历史变量相关的非线性函数；fy为钢材的屈服强度；Es为弹性模量；b1为第二刚度系数；a1～a4为控制同性强化的系数。

3 数值混合模拟验证

3.1 数值混合模拟原理

在模型更新过程中，模型参数的识别方法主要分为在线识别和离线识别两类[22-23]。为了提高复杂结构模型更新数值混合模拟精度，本文主要应用在线参数识别方法，验证笔者所提方法的有效性。基于多尺度模型更新混合试验方法的基本思想是：通过对整体结构加载，可以获得两种层面物理子结构的准确加载命令，并能够考虑试验对象与结构其它部分之间的耦合作用影响；物理子结构当前步恢复力不返回给整体结构，仅用于在线估计模型参数，通过估计的模型参数更新整体结构模型参数进而得到结构结构反应，从而评价整体结构抗震性能。

该数值混合模拟方法包括四个模块：整体结构分析模块、物理子结构试验模块、材料模型参数识别模块、构件模型加载及构件模型参数识别模块。首先，根据试验条件确定物理子结构，通过作动器实现边界条件。采用OpenSees有限元软件进行整体结构加载分析，MATLAB软件进行参数识别，OpenSees有限元软件和MATLAB软件之间的数据通讯基于TCP/IP协议的Socket通讯技术实现[24]。数值混合模拟的原理示意图如图5所示。

图5 多尺度模型更新数值混合模拟原理

3.2 多尺度模型更新数值混合模拟方法步骤

步骤(1)对于整体结构，将二层框架部分在OpenSees端离散化为数值模型，SCED在MATLAB端利用模型方程建立数值模型，同时在整体结构中选择用于提供模型参数更新部分的物理子结构，材料层次选用结构底层钢框架柱作为物理子结构，构件层次选用结构底层SCED作为物理子结构。

步骤(6)以最新识别得到的钢框架柱材料本构模型参数及SCED模型参数更新所对应整体结构数值模型的材料本构模型及SCED模型参数。

步骤(7)重复步骤(3)～(7)，直至试验结束。

3.3 参数设置

为了验证本文所提出自复位摩擦耗能支撑结构多尺度模型更新数值混合模拟方法的可行性，针对二层带有SCED结构的二层钢框架进行数值混合模拟，设定钢框架底柱与基础为刚接，该框架层高为3.6 m，跨度为6.0 m。本文为方便建模，梁柱截面均采用H型钢，材料参数为(HW300×300×10×15)。钢框架的梁和柱采用基于位移的梁柱单元(element disp Beam Column)，每个单元选取3个Gauss-Lobatto积分点，钢材截面为纤维截面。在对截面进行纤维单元定义时，为了保证试验精度且快速完成计算，分别沿局部坐标系的两个方向对截面进行划分，沿长度方向划分出20个子区域，沿宽度方向划分出2个子区域。钢材纤维的本构模型选择具有各向同性应变硬化的单轴 Giuffre-Menegotto-Pinto(Steel02)钢材料对象。对于物理子结构的选取，在子结构试验中通常将子结构的边界建立在反弯点处，以避免模拟转动自由度。因此，本文除将二层框架结构下层SCED作为物理子结构外，还将底层钢框架柱的一半作为物理子结构。

本次数值混合模拟二层框架结构及结构底层钢框架柱采用OpenSees模拟，SCED采用MATLAB模拟。本次试验仿真假定：试验子结构与数值子结构采用相同模型、试验子结构先于数值子结构进入非线性、试验子结构和更新的数值子结构加载路径差异对结构无影响。框架结构模型参数如下所示：每层结构质量为Mn1=Mn2=2 000 t，框架结构的层间水平刚度刚度为Kn1=Kn2=80 000 kN/m，框架结构的层间阻尼系数为Cn1=Cn2=1 550 kN/(m/s)，此处n1、n2分别为数值结构的一层和二层，支撑与楼面夹角均为28.81°，每层SCED结构恢复力模型采用MFS模型。地震作用选取EI Centro(1940,NS)地震波，地震动峰值加速度为400 cm/s2，加载制度如图6所示。

图6 EI Centro地震动加载Fig.6 EI Centro loading system

本次模拟采用CKF方法识别物理子结构的模型参数，识别物理子结构的模型参数包括：下层支撑滞变位移、加载位移及MFS模型关键参数b、K、Q、钢框架柱屈服强度fy、弹性模量E、钢材硬化系数b1，物理子结构的模型参数值如表1所示。利用识别关键参数值去更新数值子结构模型参数，通过模型参数识别结果及整体结构的震后响应验证多尺度模型更新数值混合模拟的可行性。本次识别过程中，物理子结构的识别分为两个部分，第一部分为底层钢框架柱识别部分，第二部分为底层MFS模型识别部分，此处需要明确的是，两次识别在整个模拟过程中同时进行同时结束，以确保试验准确性。

表1 模型参数Tab.1 Model parameters

(28)

式中,x为加载位移。

y=x4·x·A+(1-A)·(1-x5)·x4·R(x)+

(1-A)·x5·x4·x1+VQ1

(29)

(30)

(31)

VQ1=10-26×I5

(32)

式中,I5为5×5的单位矩阵。

WR1=10-10×I1

(33)

式中,I1为1×1的单位矩阵。

VQ2=0.4×I3

(34)

式中,I3为3×3的单位矩阵。

WR2=0.4×I1

(35)

(36)

(36)

(37)

CKF识别模型参数结果如图7所示。

(a) 滞变位移识别值

使用CKF方法得到的支撑恢复力如图8所示。

(a) 上层支撑恢复力

使用CKF方法得到的支撑滞回曲线及整体结构滞回曲线如图9所示。

(a) 下层支撑滞回曲线

由图7、图8、图9可以看出，使用CKF方法得到的模型参数识别值及滞回曲线与真实值吻合程度较高，收敛至真实值的速度较快，且与真实值误差较小。从图8中可以发现，相比真实值模拟得到的支撑恢复力，基于多尺度模型更新数值混合模拟方法所得到的支撑恢复力与采用参数真实值得到的支撑恢复力吻合良好；虽然在初期由于参数识别过程的波动导致基于多尺度模型更新数值混合模拟方法所得到的支撑恢复力与采用参数真实值进行非线性分析得到的支撑恢复力存在一定差异，由于参数识别属于反问题求解，参数识别值初期都会经历一定的波动，之后逐渐收敛到一个稳态值。从耗能角度分析，CKF方法得到的下层支撑滞回曲线、上层支撑滞回曲线、整体结构滞回曲线与真实值的相对误差分别为5.69%、7.29%、6.47%，相对误差较小，可以更好地模拟结构真实响应。

为了能定量评价算法识别精度，定义独立仿真的相对误差(relative error,RE)为

(38)

模型参数识别值的相对误差如图10。

结合图7、图10可以看出，在参数识别过程中的前10 s，参数识别值并没有达到稳定阶段，参数识别值波动较为明显，导致参数识别值相对误差较大，这是因为每次运行的观测噪声是随机的，从本质上讲CKF算法是一种随机算法，每一次识别的结果都存在差异，并且算法初始参数的设置对识别结果也有一定的影响，往往试验者常根据经验估计算法初始参数值，从而导致参数识别前期产生较大的波动。此外，模型参数敏感性对识别结果也有一定的影响，敏感性小的参数识别起来相对难度更大。随着时间的增长，多尺度模型更新数值混合模拟方法得到的模型参数识别值相对误差逐渐降低，最终达到平稳阶段。MFS模型参数b、K、Q的平均相对误差分别为18.75%、4.38%、13.69%，钢材料模型参数fy、E、b1的平均相对误差分别为11.23%、5.26%、19.56%，可见多尺度模型更新混合模拟方法具有较好的识别精度。

图10 参数相对误差图Fig.10 Relative error diagram of parameters

4 对比验证

为了能够更好说明多尺度模型更新数值混合模拟方法的优势，分别单一对钢框架模型参数及单一对构件模型进行模型更新，数值模拟过程中单纯考虑更新一种层面上的模型参数，整体结构滞回曲线如图11所示。通过对二层整体钢框架结构的耗能、残余变形、顶层相对位移及最大层间位移角分别进行分析，以验证多尺度模型更新数值混合模拟方法的可行性及有效性。

由图11(a)可以看出，单独进行钢材料模型更新对整体结构滞回曲线的结果影响很大，不足以反映结构真实响应。为了定量分析整体结构的耗能能力，采用结构的累计耗能值作为衡量指标，即采用EI Centro地震动加载产生滞回环的面积作为衡量的指标。从耗能的角度分析，单纯对钢材料模型更新得到的滞回曲线与真实值耗能相对误差为34.35%，而多尺度模型更新数值混合模拟方法得到的滞回曲线与真实值耗能相对误差为6.47%，多尺度模型更新数值混合模拟方法相对误差降低了27.88%。

(a) 钢材料模型更新

由图11(b)可以看出，单独进行MFS模型更新对整体结构滞回曲线的结果影响很大，与结构真实加载响应存在较大差异。从耗能的角度分析，MFS模型更新与真实值耗能相对误差为29.32%，而多尺度模型更新数值混合模拟方法与真实值耗能相对误差为6.47%，多尺度模型更新数值混合模拟方法相对误差降低了22.85%。

通过单独对钢材料模型更新、MFS模型更新以及多尺度模型更新数值混合模拟三种方式，分别计算整体结构的顶层相对位移，验证所提出方法的有效性，结构顶层相对位移如图12所示。

(a) 顶层相对位移

由图12(a)可以看出，单独对钢材料模型更新及单独对MFS模型更新所得到的顶层最大相对位移分别为0.120 m和0.125 m，产生的最大层间位移角分别为1/60和1/57.6，而多尺度模型更新数值模拟方法得到的顶层最大相对位移仅为0.095 m，产生的最大层间位移角为1/75.79，顶层最大相对位移与真实值的相对误差为2.96%，最大层间位移角与真实值的相对误差为5.6%，与单独对钢材料模型更新和单独对MFS模型更新所得到的相对误差相比，最大相对位移相对误差分别降低了20.41%和25.55%，最大层间位移角分别降低了28.23%和30.02%。

由图12(b)可以看出，单独对钢材料模型更新及单独对MFS模型更新所得到的结构残余变形分别为5.72×10-5m和5.95×10-5m，结构残余变形与真实值相对误差分别为25.44%和30.48%，而多尺度模型更新数值混合模拟方法得到的结构残余变形为4.63×10-5m，和真实残余变形相对误差为1.53%，与单独对钢材料模型更新及单独对MFS模型更新所得到的相对误差相比，分别降低了23.91%和28.95%。相对误差图如图13所示。

图13 相对误差图Fig.13 Relative error diagram

图13为三种不同模型更新方法对整体结构的滞回曲线耗能量、残余变形、顶层相对位移及最大层间位移角的相对误差图。由图13可知，多尺度模型更新数值混合模拟方法与钢材料模型更新方法及MFS模型更新方法相比，多尺度模型更新数值混合模拟精度更高。

5 结 论

本文提出自复位摩擦耗能支撑结构多尺度模型更新数值混合模拟方法，以二层带有SCED的钢结构框架为例，通过使用CKF算法进行多尺度模型更新数值混合模拟。得出以下结论：

(1)基于CKF的多尺度模型更新数值混合模拟方法解决了大型复杂结构的模型更新，提高了含有多种不同类型构件结构的模型更新数值混合模拟精度，扩展了数值混合模拟的工程应用能力。

(2)与单一材料模型更新传统数值混合模拟方法相比，基于CKF的多尺度模型更新数值混合模拟方法极大程度降低了耗能、残余变形、顶层相对位移、最大层间位移角的相对误差，验证了多尺度模型更新数值混合模拟方法的有效性，具有广阔的应用前景。