第 43 卷第 4 期2023 年 8 月
振动、测试与诊断Vol. 43 No. 4
Aug.2023 Journal of Vibration，Measurement & Diagnosis
基于应变响应的结构动态载荷识别方法∗
郑国峰1，2，陈柏先1，陈文1，赵树恩1，肖攀2，刘晓昂3
（1.重庆交通大学机电与车辆工程学院重庆， 400074）
（2.中国汽车工程研究院股份有限公司重庆， 401122）
（3.河北工业大学机械工程学院天津， 300401）
摘要提出了基于结构应变响应的动态载荷识别方法。首先，开展结构动态响应理论分析，推导结构动态响应的应变表达，构建基于应变响应的动态载荷识别模型；其次，基于D‐optimal算法，对应变响应采集所涉及的贴片数量和方位进行优化计算；然后，以方形悬臂薄板和汽车拖曳臂的动态载荷识别为例，分别进行正弦和随机载荷的识别，分析识别激励的误差来源；最后，采用Pearson相关系数和拟合系数，衡量识别
载荷与原始载荷的相似程度和幅值差异情况。结果表明，方形薄板识别激励与原始激励吻合较好，拖曳臂识别激励与原始信号的波形基本保持一致，但幅值有一定的差异。
关键词应变响应；载荷识别；应变采集；D‐optimal算法
中图分类号U467.1；TH139
引言
机械装备部件或结构的动态服役载荷获取可通过直接法和间接法实现。基于结构的传递特性，结合结构动力学反演模型，通过易于采集的结构响应信号来重构零部件动态激励的间接法已成为机械装备零部件服役载荷谱的重要获取方法［1］。
基于结构表面振动响应的激励识别得到国内外学者的广泛研究［2］。工程中，结构动态载荷反演方法主要有时域法、频域法和智能算法。按理论基础不同，频域法分为传递函数求逆［3］、传递函数正则化［4‐5］和逆虚拟激励法等。传递函数求逆和传递函数正则化适用于确定性载荷识别，逆虚拟激励法适用于平稳随机动载荷。时域法分为基于模态空间变换与Duhamel积分展开法［6］、基于状态空间法［7］以及函数逼近法［8‐9］等。智能算法包括遗传算法［10］、神经网络［11］以及改进粒子算法［12］等。针对复杂非线性和不确定结构等问题，智能算法开辟了有效的反演途径。文献［13］研究了
动态载荷反演，在提高算法效率和精度方面有所发展。但针对传感器的优化布置、响应信号的准确采集等影响激励识别的关键因素，鲜有研究报道。
笔者构建了基于结构应变响应的动态载荷识别模型。该方法通过便于测量的应变响应，结合系统的应变模态响应状态，实现动态载荷的识别。对贴片位置、数量和方向进行优化计算，获得用于高精度动态载荷识别的贴片位置和方向。采用Pearson相关系数和拟合系数，衡量识别载荷与原始载荷的相似程度和幅值差异情况，分析了识别过程的误差来源。将该方法运用到工程实际中，对方形悬臂薄板和汽车底盘结构拖曳臂的动态载荷进行识别。识别结果表明：方形薄板识别激励与原始激励吻合较好，拖曳臂识别激励与原始信号的波形基本保持一致，但幅值有一定的差异；导致幅值误差的主要原因是在模态参与因子计算时采用了截断模态。
1 动态载荷识别方法
1.1　结构应变与激励力关系标定　
在工程实际中，对结构简单或受力简单的部件，其服役过程中动态载荷一般采用贴片和标定的方法实现，即在表面粘贴应变片，加载已知力并输出结构表面的应变信号，建立应变与已知加载力的线性关系。利用线性关系将服役过程中采集到的应变信号DOI:10.16450/jki.issn.1004‐6801.2023.04.022
∗ 国家自然科学基金资助项目（52072054）；重庆市自然科学基金面上资助项目（CSTB2022NSCQ‐MSX1266 ）；中国博士后科学基金面上资助项目（2022M713438）；重庆市教委科学技术研究计划资助项目（KJQN2021000713）；重庆交通大学研究生科研创新资助项目（YYK202206）
收稿日期：2021‐07‐09；修回日期：2021‐12‐01
振动、测试与诊断第 43 卷
转换为力信号，即
ε=AF（1）其中：ε为结构表面的应变；A为标定系数矩阵；F为外部动态激励载荷。
动态激励载荷与贴片数不同时，标定矩阵为非方阵，此时应变转换为力信号，需要利用最小二乘法对服役过程中的动态激励进行计算［14‐15］
F=(A T A)-1A Tε（2）测量的应变信号与期望值之间有一定的误差，且每个测量应变信号之间的标准差相对独立。动态激励载荷的计算误差［16‐17］可表示为
Err(F)=σ2(A T A)-1（3）根据式（3）可知，若使动态激励载荷的计算误差越小，则要求(A T A)-1越小，即A T A越大。由于A 与贴片位置有关，若使A T A越大，需要选定合适的应变测点，使标定系数矩阵的行列式最大化，此时该问题归结为D‐optimal优化问题［14］。
1.2　结构动态响应　
结构在运动过程中的振动与响应可利用有限元法将连续结构离散为具有n个自由度的系统，在外部动态激励载荷F作用下构建振动模型
M nn x n+C nn x n+K nn x n=F n（4）其中：M nn，C nn和K nn分别为系统的质量矩阵、阻尼矩阵和刚度矩阵；x n，x n和x n分别为单元的加速度、速度和位移信号。
将系统的位移响应转换到模态坐标系下
x n=ϕnn q n（5）其中：ϕnn为模态矩阵；q n为模态参与因子。
将式（5）代入式（4），得到振动方程在模态坐标系下的表达式为
M nnϕnn q n+C nnϕnn q n+K nnϕnn q n=F n（6）1.3　基于应变响应的动态载荷识别　
基于有限元理论，连续结构体单元内任一点的应变ε与单元节点位移x满足以下关系［16］
ε=Bx（7）其中：B为应变矩阵，是单元形函数对坐标求导后得到的矩阵。
基于结构振动方程，通过式（5）和式（7）获得运动过程中结构的m个单元（m=n/6）的应变响应，即ε
m=B mn x n=B mnϕnn q n=ψmn q n（8）其中：ψmn为系统的应变模态矩阵。
ψmn=B mnϕnn是系统各阶模态响应下的应变表达，为系统固有属性，通过试验或模态分析获得。根据式（8）可知，在某时刻结构上一点的应变响应可表示为该点处各阶应变模态的线性组合。
对于连续系统，当结构形式确定后，应变模态矩阵即可确定。结构的应变响应可通过应变片测量获取。系统在振动过程中的模态参与因子为
q n=(ψT mnψmn)-1ψT mnεm（9）在工程实际中，保留足够的模态阶数能准确识别出外部动态载荷。一般选取较低阶应变模态矩阵ψmn求解模态因子，保留的模态阶次以模态参与因子至少占整体的95%为标准［14］。选取前k阶截断应变模态矩阵ψmk，得到模态参与因子为
q k=(ψT mkψmk)-1ψT mkεm（10）结合应变片采集的应变响应εm，通过式（5）求得系统在前k阶模态参与条件下的位移响应。通过数值求导得到相应速度和加速度响应，并结合式（6）对结构受到的动态载荷进行识别。
1.4　应变测点及方位的优化计算　
1.4.1　应变测量及坐标转换　
计算系统模态参与因子时，要求分别进行系统应变模态矩阵识别和结构表面的应变采集。系统的应变模态矩阵一般通过有限元法获取，结构表面的应变采集通过贴片获取。当布置的测点越多，获得的模态参与因子就越准确。需要考虑结构各阶振动响应，采集应变模态矩阵敏感点，才能准确识别动态加载力信号，并要求应变个数必须大于或等于截断应变模态的阶数，使模态参与因子的计算为超定问题。
由于应变测试只能在结构表面进行，因此采用有限元法对结构进行模态响应分析。获取结构表面的应变响应时，二维壳单元将提取单元形心处的应变，三维实体单元则在表面覆盖壳单元，通过共节点处理后提取壳单元形心处的应变。
系统应变模态矩阵ψmk中的行表示单元个数，列代表模态响应的阶次，元素为应变张量。由于单轴应变片只对轴向的应变敏感，因此随着测量角度的不同，测得的应变也不同。将所有单元在每个阶次下的应变张量进行坐标转换
εx' y' z'=Tεxyz T T（11）其中：T为坐标转换矩阵。
二维壳单元的单元坐标系仅涉及绕z轴的旋转
780
第 4 期
郑国峰，等：基于应变响应的结构动态载荷识别方法
变换，此时坐标转换矩阵可表示为
T =éëêêêùû
úúúúcos θsin θ0-sin θ
cos θ000
1（12）
其中：θ为变换角度，变化范围为0°~180°。
考虑到计算效率等问题，设置变换角度的增量为10°，共18个增量。通过式（11）获取每个模态阶次下所有单元在0°~180°的应变张量。每个单元的应变张量记为εl ij ，其中：上标l =1，2，…，m ，为单元个数；i =1，2，…，k ，为截断模态的阶次；j =0°，10°，…，
170°，为应变的坐标变换角度。在18个坐标转换角度下，所有单元的应变张量组成系统的应变模态集ψ
CS mk ={εl ij }。  1.4.2　应变数量及测试方向优化
根据式（10），利用截断应变模态矩阵计算模态
参与因子，其与期望值有一定的误差。模态参与因子的计算误差表达式为
Err (q k )=σ2(ψ
T mk ψmk )
-1
（13）
其中：
σ为应变测量信号的标准差。模态参与因子的误差越小，要求(正则化损伤识别matlab
ψT mk ψmk )
-1
越
小，即行列式|ψT mk ψmk |越大。其中，截断应变模态
矩阵ψ
mk 在应变模态集ψCS mk ={εl ij }中选择。基于模态参与因子的误差分析，寻贴片位置和贴片方向可归纳为数学问题：在截断应变模态集ψCS mk 中，寻g 个贴片位置（g ≥k ）及其最佳感应方位，使行列式|ψT ψ|最大，即
ìíî
ïïmax ()
||ψ
T ψ      ψ∈{}εl ij    (i =1，2，…，k ；j =0°，10°，…，170°；l =1，2，…，m )s.t.      l =g ，g ≥k
（14）
上述优化问题的具体步骤［18］如下。1） 在应变模态集中寻g 个点的应变粘贴位置和0°方向构成初始矩阵ψ
。
2） 在剩余应变模态集中寻一个点及某个角
度方向，将相应列向量增加到初始矩阵的最后一列，使增加后的矩阵行列式|ψ
+T
ψ+|值最大。其中：ψ+为由g 列矩阵增加一列产生的g +1列矩阵。3） 去掉矩阵ψ+中的一列（不含新增的第g +1
列），使行列式|ψ
-T
ψ-|值最大。其中：ψ-为ψ+中去掉一列后产生的新的g 列矩阵。
4） 循环步骤2和3，在g 个矩阵中不断地增加和
减少列，直至行列式|ψT ψ|最大，得到最佳的贴片位
置和贴片角度。1.5　动态载荷识别流程
动态载荷识别流程主要包含以下步骤。1） 结构有限元分析。在Abaqus 软件中对结构进行网格划分和约束模态分析，获取结构应变模态矩阵ψmn 。利用输出刚度和质量矩阵的插件，导出模型单元及全局的质量、刚度和载荷矩阵，及结构应变模态矩阵，并自动转换成Matlab 矩阵。
2） 确定截断模态阶次。根据模态分析结果，基于保留的模态数量以模态参与因子至少占整体95%的原则，确定截断模态的阶次。
3） 计算最优贴片位置和贴片方向。在Matlab 软件中对应变张量进行坐标转换，形成在不同坐标
转换角度下，所有单元应变张量组成的截断应变模态集ψCS mk 。采用D‐optimal 算法，在截断应变模态
集中同时寻g 个贴片位置和贴片最佳方位，使行列式|ψT mk ψmk |最大。贴片位置和方向确认后，同时也确定结构的截断应变模态矩阵ψ
mk 。4） 截断模态参与因子计算。对结构进行无阻尼瞬态响应分析，根据确定的贴片位置和方位，输出相应的应变数据εm （视为测试应变数据）。根据截断应变模态矩阵，计算截断模态参与因子。
5） 动态载荷识别。根据式（10）得到的截断模态参与因子，结合式（5）计算结构的位移响应，通过求导获得速度和加速度响应，并将结果代入式（4），即可对动态载荷进行识别。
2 平面薄板动态载荷识别的仿真验证
笔者以100 mm×100 mm×2 mm 方形悬臂薄板为例进行动态载荷识别。采用正弦和随机载荷对薄板进行激励，基于薄板模态应变响应，运用笔者提出的动态载荷识别方法对载荷进行识别。
对平面薄板结构进行网格划分，形成20 mm×20 mm 的离散结构单元，得到结构单元数m =25，自由度n=150。首先，进行约束模态分析，获取结构的应变模态矩阵；其次，基于方形薄板的约束模态分析和瞬态响应分析结果，在正弦信号F (t )=200 sin (30t )和随机载荷激励的作用下，得到方形薄板在外部激励
781
振 动、 测试与诊断第 43 卷
作用下前10阶模态参与因子，如表1所示。
在2种激励作用下，前5阶的模态参与因子分别占整体的95.14%和95.27%，符合高于95%的要
求，因此截断模态阶次至少为5阶。2.1　正弦信号的识别结果与分析　
确定截断模态的阶次后，将应变模态在Matlab 中进行坐标转换，采用D‐optimal 优化算法确定贴片的数量和方向。为确定截断模态阶次对贴片数量及动态载荷识别结果的影响，计算分析了截断模态为5~10阶的贴片位置、数量和方向，以及动态载荷识别的情况。
图1为不同截断模态下薄板贴片位置及贴片方位示意图。图中，方形薄板被划分为25个单元，左上角数值表示每个单元的编号，载荷垂直作用于薄板的右下角处，单元中心的短黑线表示应变片。各阶截断模态下，水平短黑线表示应变片粘贴方向与x 轴的夹角为0°，斜短黑线表示应变片粘贴方向与x 轴的夹角为80°。
按照图1的形式对方形薄板进行贴片，采集其应变响应，采样频率为0.1，0.2，0.5，1.0和2.0 kHz 。基于应变响应对动态载荷进行识别，得到薄板在不同阶次截断模态下不同应变采样率的动态载荷识别结果，如图2所示。可见，不同应变采样频率下识别载荷与原始载荷的吻合度不同。
采用Pearson 相关系数衡量识别载荷与原始载荷之间的吻合度，Pearson 相关系数［19］ 定义为
r =
M
∑i =1
M x i
y i
-∑i =1M x i
∑i =1
M
y
i
M
∑i =1
M
x
2
i
-
()
∑i =1
M
x i 2
M
∑i =1
M
y
2i
-
()
∑i =1
M
y i
2
（15）
表1　方形薄板在外部激励作用下前10阶模态参与因子Tab.1　First ten modal participation factors of
square thin plate under external excitation
阶次12345678910
正弦载荷模态参与因子4.81-1.240.110.040.060.030.060.010.02-0.20
百分比/%73.1018.841.670.610.910.460.910.150.303.04
随机载荷
模态参与因子2.46-0.650.060.020.030.010.030.010.01-0.10
百分比/%
72.7819.231.780.590.890.300.890.300.302.96
图1　不同截断模态下薄板贴片位置及贴片方位示意图
Fig.1　Schematic diagram of position and orientation of strain gauge attach to thin plate under different order truncated modes
782
第 4 期
郑国峰，等：基于应变响应的结构动态载荷识别方法
其中：x 为原始载荷信号，其数据记为x i ，i =1，2，…，M ；M 为原始载荷的数据总数；y 为识别载荷，其数据记为y i ，
i =1，2，…，M 。表2为不同采样率下识别载荷与原始载荷Pearson 相关系数。可见，随着采样率的增大，识别载荷与原始载荷的相关系数变大。
对识别载荷与原始载荷进行拟合，设拟合截距为0。拟合得到曲线与y =x 越接近且相关系数越大，则说明2条曲线的波形越相似。拟合系数与1的差异越大，表示信号的幅值差异越大，差异的程度以拟合系数为倍数。不同应变采样率下识别载荷与原始载荷拟合系数如表3所示。
图3为不同截断模态下相关系数变化趋势，根据图3得到以下结论。
1） 应变采样率从0.1 kHz 增加到0.5 kHz 的过程中，相关系数的变化趋势较大；当采样频率达到0.5 kHz 后，相关系数的变化趋势较为缓和，除了采用第9阶截断模态进行识别的载荷以外，其余阶次识别出的动态载荷与原始载荷的吻合度均达到0.98以上。
2） 在相同应变信号采样频率的情况下，第9阶截断模态比其他阶次模态得到的相关系数小。基于薄板有限元模态分析结果可知，在该阶段模态下，有多个应变片安装在振型节点上，导致测量结果较小，从而使载荷识别精度低。
根据图2，3和表2，3得到以下结论。
1） 应变采集信号的采样率越高，识别得到的载荷精度越高，但此时识别载荷的噪声严重，需要对识别载荷进行滤波处理。
2） 根据图2可知，第8阶截断模态对应识别
载
图2 薄板在不同阶次截断模态下不同应变采样率的动态载荷识别结果
Fig.2 Dynamic load identification results of thin plate under different order truncated modes and different strain sampling rates
表2 不同采样率下识别载荷与原始载荷Pearson 相关系数Tab.2 Pearson coefficient between identification load and
original load with different strain sampling rate
截断模态阶次5678910
采样率/ kHz
0.10.911 80.909 70.909 00.911 80.853 40.905 7
0.20.976 20.973 80.973 10.976 00.911 90.969 3
0.50.995 00.992 60.992 20.994 50.938 10.987 4
1.00.997 70.995 30.994 80.997 10.940 60.990 0
2.00.998 40.996 00.995 50.997 80.987 30.990 7
表3 不同应变采样率下识别载荷与原始载荷拟合系数Tab.3 Fit coefficients between identification load and
original load with different strain sampling rate
截断模态阶次5678910
采样率/ kHz
0.117.133 00.635 90.735 30.991 50.220 20.410 6
0.217.624 00.652 90.754 71.022 70.224 20.424 5
0.517.868 00.663 40.769 21.022 10.236 80.419 3
1.017.891 00.663 90.769 81.022 30.237 10.419 5
2.017.893 0
0.663 90.769 91.022 00.234 80.419 2
783

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。