柔性体的有限元模型修正技术研究

由上图可以看出第 5、第 6 个节点以及第 10、第 11 个节点的残余力远大于零，所以可
-6-
http://www.paper.edu.cn
以判断这四个节点存在着损伤。第 5、第 6 个节点对应梁的第 5 单元，第 10、第 11 个节点 对应梁的第 10 单元。至此，我们可以确定第 5 单元和第 10 单元存在损伤或者建模误差。 然后利用本文所提出的方法求解损伤程度，经过上述的化简方法，20 个单元只含有 40 个方程，每个单元会求出两个损伤系数。求得各个单元的损伤系数见表 3 所示。
2．基于残余力的结构损伤识别
对于发生损伤的有限元模型，其特征方程为：
* (K * − λ * j M)φ j = 0
(1) (2)
*
* *
K * = K − ∆K
其中 M 表示质量矩阵，并假设损伤后结构的质量不发生变化， K, K 分别表示损伤前 后系统的刚度矩阵， ∆K 表示损伤前后刚度矩阵的变化，它们均为 n 阶方阵。 λ j , φ j 分别表 示损伤后第 j 个特征值和振型。
⎡ c j1 ⎤ ⎢c ⎥ j2 cj = ⎢ ⎥ ⎢ # ⎥ ⎢ ⎥ ⎢ ⎣c j q ⎥ ⎦
其中 c ji = b js + b jd ，其中 q 为节点个数
(7)
3．改进最小范数法
为了求解损伤的程度，通常都会对如下形式的方程利用最小范数法进行求解：
α = S +b j
(8)
其中上标“+”表示广义逆。通常情况下，矩阵 S 包含了振型 φ d j 的信息，所以振型的测 量噪声会对 S 的求广义逆产生很大的影响，导致计算结果出现很大的误差，常常会导致解 的不稳定， 也是最小范数法不能成功用于工程实践中的根本原因之一。 下面对利用最小范数 法求解上述问题的方法进行改进。 首先对单元刚度矩阵 K i 进行 QR 分解，
(4)
记向量 b j = (K − λ j M)φ ，方程(4)可改写为：
* *
∆K φ j* = b j
上式即为残余力方程， b j 被称为第 j 个残余力向量。 (5)可写成如下形式：
(5)
⎧b * ⎫ ⎡∆K 1T ⎤ j1 ⎪ * ⎪ ⎢ T⎥ b j2 ⎪ ⎢∆K 2 ⎥φ * = ⎪ ⎨ ⎬ j ⎢ # ⎥ # ⎪ ⎪ ⎢ T⎥ ⎪b * ⎪ ⎢∆K n ⎦ ⎥ ⎣ ⎩ jn ⎭
⎡α 1 ⎢ α1 ⎢ ⎢ % ⎢ α1 ⎢ ⎢ α2 ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣
+ = C nm *n b j
α2
%
α2 αm αm
⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ Dφ jd nm*nm ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ % ⎥ αm ⎥ ⎦ nm*nm
(11)
-3-
http://www.paper.edu.cn
4．数值算例
仿真实验对象为一根均匀铝制梁（如图 1 所示） ，密度 2900kg/m3。构件总长 60CM，宽 度 2CM，厚度 2CM，将梁划分为 20 单元。
单元1
单元2
单元20
图 1 悬臂梁
表 1 节点和单元标号之间的关系
单元
对应节点号
单元
对应节点号
1 2 3 4 5 6 7 8 9 10
1 1，2 2，3 3，4 4，5 5，6 6，7 7，8 8，9 9，10
d nm*nm
中会产生很多为零的项，那么在方程中所对应的行将没有意义，所以
d nm*nm
可以将其省略，即将 C 相应的列划去，同时将 Dφ j
以及 b j 中相应的行划去。这将使
改变后的矩阵 C 的列数大大降低。同时，利用损伤识别技术可以减少待修正的单元个数， 从而能够求得更加精确的广义逆，得到更加精确的修正结果。 改进的最小范数法仅需要计算矩阵 C 的广义逆，而矩阵 C 是不依赖于模态信息的，所 以测量噪音不对矩阵 C 的求逆产生影响，并且其物理意义明确，这个优点将使得改进的最 小范数法能够获得比较稳定和精确的识别结果。
11 12 13 14 15 16 17 18 19 20
10，11 11，12 12，13 13，14 14，15 15，16 16，17 17，18 18，19 19，20
-4-
http://www.paper.edu.cn
算例 1： 假设第 5 单元结构存在建模误差或者损伤， 损伤后的杨氏模量变为损伤前的 0.8 倍。即 E5 = 0.8 E5 。下图为利用残余力向量进行损伤识别得到的节点上的残余力。
-5-
http://www.paper.edu.cn 表 2 利用改进最小范数法计算结果
单元
损伤系数
-0.00000000000105 -0.000000000002045 0.000000000004397 0.000000000001083 0.000000000008541 0.000000000003702 0.000000000014068 0.000000000005502 0.199999999988215 0.200000000005876 -0.000000000001764 0.000000000010936 -0.000000000002404 0.000000000029788 -0.00000000000615 -0.000000000050027 -0.000000000009485 -0.000000000021455 -0.000000000027635 -0.000000000019483
http://www.paper.edu.cn
柔性体的有限元模型修正技术研究1
李伟明，洪嘉振
上海交通大学工程力学系，上海 (200240)
E-mail：lwming2001@sjtu.edu.cn
摘 要：本文基于残余力向量的损伤识别方法，提出了一种改进的最小范数法，利用 QR 分 解技术， 在求逆过程中消除掉含有噪声的模态项。 通过严格的公式推导证明该方法只需要少 量的模态信息就可得到满意的模型修正与损伤识别结果。 关键词：模型修正；损伤识别；残余力；噪声 中图分类号：O342； V214
1
本课题得到国家自然科学基金（项目编号：10772113）及高等学校博士学科点专项科研基金（项目编号：
20040248013）的资助。 -1-
http://www.paper.edu.cn
代入到(1)得：
* (K − ∆K − λ * j M)φ j = 0
(3)
整理可得：
* ∆K φ j* = (K − λ * j M)φ
*
1 2 3 4 5 6 7 8 9 10
11 12 13 14 15 16 17 18 19 20
算例 2 ：假设第 5 单元与第 10 单元发存在建模误差或者损伤，即 E5 = 0.8 E5 ，
* E10 = 1.1E10 。下图为利用残余力向量进行损伤识别得到的节点上的残Fra Baidu bibliotek力。
图 2 节点上的残余力
1．引言
在柔性多体系统动力学中， 由于结构的复杂性， 通常需要利用有限元以及模态截断技术 来对其进行仿真分析，这就需要能够首先掌握构件真实的物理参数，例如杨氏模量，密度， 尺寸以及阻尼参数等。因为即使是很小的误差，也会对仿真结果带来很大的影响。工程实际 中精确测量材料的物理参数是很难实现的， 经常发现由有限元模型计算的特征值和特征向量 等动态特性参数和试验所得对应的固有频率和振型等模态参数相差甚远。 一个有限元模型只 有在被证明具有一定的可信度时， 才能被应用于柔性多体系统动力学仿真以及进行结构的主 动控制等， 而实验模态模型可以用于验证有限元模型的可靠性。 这就需要利用模型修正技术 对有限元建模进行修正[1~3]。 目前我们所研究的模型修正方法绝大多数不能准确找到模型中究竟哪里存在误差， 许多 不需修改的物理参数也得到了修改， 即使修正后的模态参数能够与真实值相对应， 但是它不 符合实际物理意义，这样就不能够从根本上解决问题。而且在实际中，经过实验测得的模态 参数都是含有噪声的， 几乎所有的修正方法对含有噪声项较大的模型进行修正都会遇到计算 结果不稳定的问题。 模型建模有误差，实际上也就相当于有限元模型与实验模型相比发生“损伤”，所以本文 利用模型修正与损伤识别方法的相似性，首先找到存在误差的部位[4~5]，减少待修正参数的 个数，然后再利用一种改进的最小范数法进行修正，利用 QR 分解技术，经过严格的数学推 导，在求逆过程中消除掉含有噪声的模态项，从而实现准确的误差定位，降低求解过程中噪 声的影响。
∆K = ∑ K i = α 1Q1 R1 + α 2 Q2 R2 + " + α m Qm Rm
i =1 m
αi
⎤ ⎥ ⎥ R i ⎥ % ⎥ α i ⎦ n*n
(9)
注意到上式中含有 n 个（n 为系统自由度数）相同的变量 α i 。
= [Q1
Q2
⎡α 1 ⎢ α1 ⎢ ⎢ % ⎢ α1 ⎢ ⎢ α2 ⎢ ⎢ " Qm ]n*nm ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣
α2
%
α2 αm αm
⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎡ R1 ⎤ ⎥ ⎢R ⎥ ⎥ ⎢ 2⎥ ⎥ ⎢ # ⎥ ⎥ ⎢ ⎥ ⎥ ⎣ Rm ⎦ nm*1 ⎥ ⎥ ⎥ ⎥ % ⎥ αm ⎥ ⎦ nm*nm
= C ΛD
(10) 其中
C = [Q1
Q2
" Qm ]n*nm
⎡ R1 ⎤ ⎢R ⎥ D=⎢ 2⎥ ⎢ # ⎥ ⎢ ⎥ ⎣ R m ⎦ nm*1
K i = Qi R i
其中 Q i 和 R i 只与 K i 有关，系统单元刚度矩阵 K i 确定，那么 Q i 和 R i 即唯一确定。 由 ∆K i = αK i ，可知单元刚度的改变量 ∆K i 改写成：
-2-
http://www.paper.edu.cn
⎡α i ⎢ ∆K i = α i Qi Ri = Qi ⎢ ⎢ ⎢ ⎣
单元
损伤系数
0.000000000014843 -0.000000000023109 -0.000000000000774 -0.000000000034966 -0.000000000004486 -0.000000000121614 -0.000000000006206 0.000000000054556 -0.000000000005029 0.000000000021858 -0.000000000169297 0.000000000014809 0.000000000012062 0.000000000013132 0.000000000008905 0.000000000016321 0.000000000005789 -0.000000000005236 0.000000000009275 0.00000000000399
(6)
其中 ∆K i (i = 1 ~ n) 表示第 i 个自由度刚度改变列向量，b ji 为 b j 中的第 i 个元素。 由(6) 可以看出， ∆K i 与 b ji 一一对应。由残余力向量中的非零元素就可以判断出发生损伤的自由 度，而由损伤自由度就可以判定出损伤单元。 为了使用上的方便并提高损伤识别的精度， 可以定义单元节点残余力向量来判别损伤位 置。
*
图 1 节点上的残余力
由图 1 可以看出第 5、第 6 个节点的残余力远大于零，所以可以判断这两个节点存在着 损伤。第 5、第 6 个节点对应梁的第 5 单元，至此，我们可以确定第 5 单元存在损伤。 利用本文所提出的方法，经过上述的化简方法，20 个单元只含有 40 个方程，每个单元 会求出两个损伤系数。求得各个单元的损伤系数见表 2 所示。
表 3 利用改进最小范数法计算结果
单元
损伤系数
-0.000000000001752 -0.000000000001065 0.000000000003 -0.000000000000202 0.000000000004378 0.000000000001934 0.000000000010145 0.000000000002889 0.199999999999871 0.200000000003102 0.000000000002548 0.000000000004175 0.000000000001727 0.000000000008561 0.000000000001426 -0.000000000006989 -0.000000000001866 -0.000000000002091 -0.100000000007586 -0.100000000002939
从上式可以看出， Λ 为（n*m）*（n*m）维，但是 Λ 却只含有 m 个未知数，也就是说 如果结构有 n 个自由度， 那么对于每一个损伤单元， 会同时出现 n 个含有同一个损伤系数的 方程。理论上说，这 n 个方程求得的解应该完全相等。 矩阵 C 为 n* （n*m）维，如果单元数目过多， 那么 C 在求逆过程中将会遇到很大问题。 但是注意到 Dφ j