《1 前言》

1 前言

一般依据结构模型计算和预测的结构响应与实测值往往差别较大。造成差别的因素很多, 主要是在建模过程中对结构的材料特性、联结形式以及外部作用进行了理想化。同时, 在结构建造过程中存在一定的施工误差, 使得依据结构设计图纸而建立的有限元模型与实际结构在构件尺寸甚至结构形态上存在一定的偏差。另一方面, 对实际结构进行测试的过程中, 测试方法和手段还不够完善, 测试数据在采集、传输和存储过程中容易受到噪声污染。为了提高模型的计算精度, 必须对建立的模型进行修正。模型修正是采取一定的方法和手段, 使建立的结构模型更精确地代表实际结构, 尽可能地减小计算和测试结果间的偏差。模型修正过程如图1所示 [1,2]

《图1》

图1模型修正过程示意图

图1模型修正过程示意图  

Fig.1 Sketch map of model updating

早在1965年, Guyan发表了模型缩聚方面的文章, 开创了模型修正的研究工作 [3]。经过几十年的发展, 现有的模型修正理论主要有两大类:Lagrange乘子法和罚函数法 (penalty function method) 。根据约束方程的不同, Lagrange乘子法又可分为:参考基准法 (reference basis methods) 、混合矩阵法 (matrix mixing) 、特征结构匹配法 (eigenstructure assignment methods) 、特征值取逆法 (inverse eigenvalue methods) 等。现有的这些模型修正方法, 由于依据的基本测量信息和修正对象的不同, 各种修正方法间差异显著, 但就具体的模型修正过程来看, 模型修正方法总体上可以分为矩阵型和参数型两类模型修正方法。关于这些经典的模型修正技术, 文献[4]中有详细的评述。

这些研究工作为实现复杂结构模型的修正工作奠定了重要基础, 但正如综述文献[5,6]在比较了众多的修正方法后所评述的那样, 现阶段要实现复杂结构模型的修正还很困难, 而困难的根本原因是模型修正理论自身, 复杂结构模型的修正理论需要进一步的研究和发展。大型结构模型修正中存在的关键问题是测量信息有限、修正目标众多以及模型修正前后的关联性, 试图通过单一的修正过程来完全解决这些关键性问题无疑十分困难, 把复杂的修正过程来分步实现成为可能 [7]。基于此, 提出多重子步模型修正方法 (MSMUM, multi-stage model updating method) , 以期为实际工程结构的模型修正提供新的方法。

《2 多重子步模型修正方法》

2 多重子步模型修正方法

为了分步实现模型修正过程, 在结构建模过程中, 借助子结构方法, 使得模型在整体上包含有限个子结构, 同时在子结构内部可以详细模拟实际结构的几何和空间构形, 以满足工程结构计算的需求。具体的模型修正实施步骤为:

Step 1 在结构整体层次上针对各个子结构进行误差定位。由于结构整体包含的子结构数目有限, 结构整体层次上误差定位过程中的未知量数目较少, 便于计算求解。

Step 2 针对误差定位后的子结构进行修正。为了修正误差定位后的子结构, 必须增加测量信息, 当测量信息完备时, 计算和测量自由度间完全匹配, 则可以实现子结构内部构件参数的修正。结构用有限数目的子结构进行划分时, 子结构所对应的是具有相当自由度的超级单元。当超级单元的内部自由度较大时, 对子结构直接进行参数修正仍然存在修正未知量过多以及自由度匹配困难的问题。此时, 在子结构内部进行类似于Step 1的误差定位过程, 进一步缩小误差修正范围。实现子结构内部的误差定位时, 可以在子结构内部嵌套子结构, 即建立多重子结构模型, 针对下一层的各个子结构进行误差定位。

Step 3 当Step 2误差定位后的子结构 (或多重子结构) 规模降低到一定程度时, 可以实现计算和测量自由度间的匹配, 方便修正构件参数, 完成误差定位子结构模型的修正。

Step 4 针对结构整体层次上误差定位后的一些子结构, 分别进行Step 2和Step 3的修正过程, 最终实现整体结构模型的修正。

上述模型修正过程见图2。

《图2》

图2多重子步模型修正过程示意图

图2多重子步模型修正过程示意图  

Fig.2 Process of multi-stage model updating

修正过程可以针对包含多重子结构的结构模型在不同层次的子结构上进行误差定位和参数修正, 具有多重修正的意义。同时, 修正各个子结构的过程相互独立, 可以对具有误差的子结构逐步修正, 因此上述模型修正过程称为多重子步模型修正方法 (MSMUM) 。

《2.1误差定位》

2.1误差定位

对包含有限数目子结构的整体模型进行误差定位仅相当于小规模结构模型的修正过程, 常规的损伤识别技术 (此时把损伤看成模型的广义误差) 和经典的矩阵型模型修正方法均可用来实现误差定位。借用矩阵型模型修正方法的基本思想, 通过修正各单元 (或超级单元) 的刚度参与系数来实现误差定位。

假设模型修正时依据的基础数据为动力测试数据, 构造目标函数为

J=i=1kΚϕ^i-λ^iΜϕ^i2(1)

其中λ^i=ω¯^i2为振动系统测试的特征值, K为结构整体刚度矩阵, M为结构整体质量矩阵, ϕ^i为测试的第i阶模态, k为测试的模态阶数。

假设系统的质量矩阵没有误差, 只有刚度阵具有误差, 且修正后模型的刚度阵可以用初始模型中各单元 (或超级单元) 的刚度阵线性表出 [8], 即

Κ=i=1nαiΚ0i(2)

式中, αi为第i个单元的刚度矩阵修正系数;K0i为原始模型中第i个单元的刚度矩阵;n为整体结构模型的单元数。

把式 (2) 代入式 (1) 有

J=j=1kQjX-λ^jRjΙn2(3)Qj=[Κ01ϕ^j,Κ02ϕ^j,,Κ0nϕ^j],Rj=[Μ01ϕ^j,Μ02ϕ^j,Μ0nϕ^j]X=[α1,α2,,αn]Τ,Ιn=[1,1,,1]1×nΤ

式 (3) 是关于X的纯二次函数, 使其取极小值, 由最小二乘法可得

JX=0(4)

j=1k(QjΤQjX-λ^jQjΤRjΙn)=0(5)

Η=j=1kQjΤQj,B=j=1kλ^jQjΤRjΙn, 上式简写为

ΗX=B(6)

求解式 (6) , 可得到待修正系数X=[α1, α2, …, αn]T的解答, 进而实现误差定位。

《2.2参数修正》

2.2参数修正

对于包含误差的子结构, 进行相应的参数修正在规模上相当于小规模结构模型的参数修正, 但有根本的区别。修正小规模的模型时, 依据整体尺度上的测量信息对结构整体模型进行修正;而子结构内部单元参数的修正过程中, 所依赖的测量信息是结构整体尺度上的测量信息, 而子结构是结构的局部。因此, 如何由整体尺度上的测量信息来实现局部模型的参数修正。

考察包含若干子结构的结构模型, 如图3所示, 对于每个子结构来说, 在参与整体结构计算时要进行内部自由度凝聚, 仅保留边界结点 (或称外部结点) 的自由度。

《图3》

图3子结构示意图

图3子结构示意图  

Fig.3 FE model with substructure

设通过适当的结点编号, 使子结构的刚度矩阵以及相应的结点位移和载荷可以分块表示如下:

[ΚbbΚbiΚibΚii]{abai}={ΡbΡi}(7)

式中ab为边界结点位移向量;ai为内部结点位移向量。

在静力凝聚的过程中, 一般通过高斯-约当消去法进行消元运算得到 [9]

S˜[ΚbbΚbiΚibΚii]=[Κ¯bb0biΚ¯ibΙii](8)

其中S˜为一系列初等矩阵的乘积, 是可逆的;0bib×i阶零矩阵;Iiii×i阶单位矩阵;Κ¯ib=Κii-1ΚibΚ¯bb=Κbb-ΚbiΚii-1Κib。子结构参与整体结构计算时的单元刚度 (即缩聚后超级单元的刚度) 为Κ¯bb

通过矩阵的分块运算, 可以得到

S˜=[Ιbb-ΚbiΚii-10ibΚii-1](9)

对子结构中的第i个单元有

fi=kssi(10)

式中ki为第i个单元的刚度矩阵;fi为第i个单元的结点力矢量;si为第i个单元的结点位移矢量。

如果子结构中共有n个单元, 如式 (10) 所示的单元结点力—结点位移关系式共有n个。把这n个未经整合的关系式写为

(f1f2fn)Τ=[diag(k1k2kn)](s1s2sn)Τ(11)

上式中把各单元的刚度矩阵ki作为矩阵主对角线上的子矩阵列入。把上式简记为

f=kss(12)

经过整合的整体刚度矩阵为K, 则结构整体的结点平衡方程为

Κδ=F(13)

式中 δ为结构的结点位移矢量, 包括结构的全部自由度;F为结点载荷矢量。

设结构的位移矢量δ与式 (12) 中的列阵s之间存在以下关系:

s=Aδ(14)

其中A为变换矩阵。

由虚功原理可以给出 [10]

Κ=AΤksA(15)

对结构中的第i个单元, 设与该单元相关的待修正的参数为Pi, 则有

ki=Τipi(16)

考虑式 (11) 和式 (12) , 可以得到

ks=diag(k1k2kn)=[diag(Τ1Τ2Τn)][diag(Ρ1Ρ2Ρn)](17)

把上式简记为

ks=Τp(18)

其中T为与刚度有关的参数关联矩阵;P为待修正的参数矩阵。

把式 (18) 代入式 (15) , 可得

Κ=AΤΤpA(19)

由式 (19) 及式 (9) , 对子结构有

S˜AΤΤpA=[Κ¯bb0biΚ¯ibΙii](20)

把上式两边左乘矩阵[0ibIii]、右乘矩阵[ab0i]T

Κ¯ibab=[0ibΙii]S˜AΤΤpA[ab0i](21)

由式 (7) 和式 (8) , 在子结构内部无结点载荷时,

Κ¯ibab=-ai(22)

由式 (9) 和式 (22) 可得

ai=-[0ibΚii-1]AΤΤpA[ab0i](23)

利用式 (23) 可以依据边界结点的响应来计算子结构内部结点的响应。当子结构内部构件的参数产生一定量的摄动δp时, 由上式可得内部结点响应的改变量为

δai=-[0ibΚii-1]AΤΤδpA[ab0i](24)

对误差定位后的子结构, 由于内部构件的参数存在一定的误差, 其内部结点响应的测量值与依据式 (23) 给出的计算值总有一定的差别, 记为Δai。设实际测量值为a*i, 则

Δai=ai*-ai(25)

δaiai, 代入式 (24) 可得

Δai=-[0ibΚii-1]AΤΤΔpA[ab0i](26)

求解式 (26) , 即可给出子结构内部构件的参数修正量Δp

在上面进行子结构内部构件参数修正的过程中, 所依据的边界结点和内部结点的测量信息均是整体结构尺度上的测量信息, 因此, 多重子步模型修正方法也可以实现一些学者提出的逐段模型修正的设想。

《3 算例》

3 算例

为了验证上面提出的多重子步模型修正过程, 考察如图4所示悬臂梁, 假设梁只发生平面变形。把悬臂梁划分为10个单元 (图4中的①, ②, …, ⑩) , 相应的结点数为11个 (图4中的1, 2, …, 11) 。采用两结点平面梁单元进行划分, 各结点的自由度为3, 整个结构的自由度为30。要利用矩阵型修正理论直接对悬臂梁整体进行误差定位至少需要完整的30个自由度的测量数据。现把悬臂梁用3个子结构进行划分 (如图4中3个虚线椭圆所示的子结构Ⅰ, Ⅱ, Ⅲ) , 每个子结构中有3个单元, 这样在包含子结构的划分方案下整体结构只有12个自由度, 整体结构的自由度显著降低。

《图4》

图4包含3个子结构的悬臂梁模型

图4包含3个子结构的悬臂梁模型  

Fig.4 Cantilever model with 3 substructures

设计了误差模式A (E=0.9E0, E=0.8E0, E=0.6E0) 和B (E=0.75E0) , 其中E=0.9E0表示子结构Ⅰ中的3个单元的弹性模量均为初始值的0.9倍, E, E的意义依此类推;E表示单元⑨的弹性模量。修正前初始模型中各单元弹性模量均取为E0。误差模式A中各子结构内部单元发生一致性误差, 而误差模式B中发生非一致性误差。

用初始模型以及在误差模式A和B下的实际模型计算的频率值如表1所示, 模型A表示误差模式A下对应的真实模型, 模型B为误差模式B下对应的真实模型。


  

表1两类误差模式下计算的前10阶频率值  

Table 1 The first 10 computed frequencies using two kinds of error modes

 

《图5》

表1 两类误差模式下计算的前10阶频率值

先采用刚度参与法进行误差定位。在误差定位过程中选取不同数目的模态进行定位, 考虑到动态测量中测试转动自由度分量非常困难, 只采用各阶模态中平动自由度分量进行分析, 由式 (6) 计算的各单元 (超级单元) 刚度参与系数随采用的模态数目的变化如图5所示。可以看出, 对误差模式A, 采用前4阶模态进行定位时, 即可比较精确地给出各单元的刚度参与系数。但对于误差模式B, 直到采用了前8阶模态进行定位时误差单元的刚度参与系数才与实际数值吻合较好, 这主要是因为子结构中的部分单元具有误差时, 所对应的超级单元刚度无法通过初始模型中的单元刚度矩阵乘以倍数来精确得到, 但从误差定位的角度, 只要能判断出刚度参与系数非1, 即可实现相应的误差定位。从图5可以看出, 选取2阶到3阶模态即可达到误差定位的目的, 这在实际的模态测量时较容易达到。

《图6》

图5刚度参与系数与选取的模态数目关系图

图5刚度参与系数与选取的模态数目关系图  

Fig.5 Stiffness coefficients versus modal numbers

根据式 (26) , 在实现误差定位后的子结构进行修正时, 每个子结构中待修正的参数为3个, 而进行整体一次性修正参数为10个, 这样每子步进行修正的参数数目显著减少, 降低了方程求解的规模。图6给出了误差模式A和B参数修正的结果, Target为相应单元参数的目标值, 1Mode为仅用一阶模态的修正结果, 2Mode为采用二阶模态的修正结果, 后者比前者的结果好。但随着采用的模态数目的增加, 相应的计算量将显著增加, 同时测量误差的影响将增大, 对实际结构应该选取合适数目的模态来进行修正, 适当模态数目的选取有待进一步研究。

《图7》

图6两种误差模式下的参数修正结果

图6两种误差模式下的参数修正结果  

Fig.6 Updated parameters under two error modes

从上面的算例可以看出, 多重子步模型修正方法能有效地实现误差定位和参数修正。用该方法对润扬长江大桥南汊悬索桥的裸塔有限元模型进行了修正。润扬长江大桥南汊悬索桥桥塔净高219 m (图7a) , 对于这样高耸的大型结构无法用常规的人工激励振动的方法进行测试, 而采用了环境激励振动测试的方法。测试系统选用11只941B型伺服式超低频加速度传感器, 其中一只传感器作为参考点, 其余10只传感器按图7b所示方式布置, 测量方向为水平方向, 分为顺桥向和横桥向两部分进行。信号采集与分析采用AZ-Cras振动与动态信号采集与分析系统, 采用环境激励下的模态参数识别方法。

《图8》

图7润扬长江大桥南汊悬索桥桥塔和传感器位置

图7润扬长江大桥南汊悬索桥桥塔和传感器位置  

Fig.7 Runyang bridge tower and the sensors location

对润扬悬索桥桥塔建立了有限元模型。考虑到桥塔属于高耸结构, 截面尺寸相对于自身纵向、横向的尺寸较小, 采用空间二结点等参梁单元来建立其有限元模型。桥塔是钢筋混凝土结构, 其各部位的钢筋配比率以及施工过程的养护条件有所差异, 各部位相应的材料参数可能差别较大, 在建立模型时采用了4组不同的材料 (M1, …, M4) 。另外考虑到塔架的各段截面尺寸有所变化, 采用了6组不同的截面参数 (S1, …, S6) , 根据桥塔的自然区段把结构划分为9个子结构 (SE1, …, SE6, SE1′, …, SE3′) , 其中SE1和SE1′, SE2和SE2′以及SE3和SE3′是对称的子结构, 参数两两相同, 桥塔的子结构划分和参数分布情况如图8所示。

《图9》

图8模型中子结构划分及参数分布情况

图8模型中子结构划分及参数分布情况  

Fig.8 Substructures and related parameter distribution in the tower model

在给定每一区段的截面参数时仔细按照结构施工图纸, 采用大型商用软件ABAQUS建立不同区段的局部有限元模型。依据局部模型, 计算出相应区段的截面参数, 各区段的截面参数是可靠的, 不需要进一步修正。建模过程中, 同时考虑了数学离散精度的影响。通过数值试验, 选定计算收敛性良好的单元疏密程度, 模型共有262个梁单元。在建立桥塔的初始模型时, 由于无法预先精确了解各区段的材料参数, 采用塔架浇注混凝土类型C50的统计材料参数作为各段的初始材料参数, 暂不考虑钢筋及其配比率的影响。

用多重子步模型修正方法对建立的桥塔有限元模型进行了修正。在误差定位的过程中, 仅采用了6只在节点处传感器测量的信息。误差主要集中在子结构SE2, SE2′, SE3和SE3′。考虑到对称性, 参数修正时只对SE2, SE3进行修正。子结构SE2中有52个单元, SE3中有24个单元, 实际进行参数修正的单元共76个。由于传感器数目的限制, 加之现场进行较大规模测试十分困难, 进行子结构内部单元的参数修正所需要增加的测量信息通过对上述10个传感器的测量值进行样条插值得到。

依据修正后模型的计算结果见表2, 第3—5列是前期根据最小秩修正方法修正后模型的计算结果 [11], 第6—8列是根据所提出的多重子步修正后的计算结果。在采用最小秩修正方法进行修正时需要对262个单元同时进行修正, 计算规模是多重子步模型修正过程的3.5倍。对比的结果可见, 用笔者提出的方法修正结果与实测值更为接近, 修正后模型计算的模态与实际测量模态的相关性更高 (MAC值更大) , 表明MSMUM是有效可靠的。


  

表2修正模型计算的结果与测试结果比较  

Table 2 Comparison of computed results via the updated model with measured ones

 

《图10》

表2 修正模型计算的结果与测试结果比较

《4 结论》

4 结论

多重子步模型修正的新方法, 先进行误差定位, 再实现参数修正, 把复杂的模型修正过程分步实现。利用子结构建模技术, 在结构整体上只包含数目非常有限的子结构所对应的超级单元, 利用较少的测量信息即可实现误差定位。通过建立子结构内部结点响应与待修正参数间的关系, 进一步对具有误差的子结构内部单元的参数进行修正。数值算例表明, 多重子步模型修正方法在进行误差定位时需要的测量信息有限, 对包含误差的局部模型进行参数修正时计算规模显著降低, 更有利于实际工程结构模型的修正。对润扬大桥南汊悬索桥实际桥塔模型的修正工作表明, 与前期利用其他修正方法修正的结果比较, 该方法减少了待修正参数的数目, 约为最小秩修正方法的29% (76/262≈29%) , 显著降低了修正过程中的计算规模, 且修正后的结果更为合理。