CN101916242B - 非传统阻尼矩阵识别方法 - Google Patents

Abstract

本发明为解决环境激励下现有模态测试技术中只有低阶的模态信息能够获得的现状，克服现有阻尼矩阵识别方法在一定程度上需要相对高阶模态信息的问题，提出新的能够处理三维空间问题的基于单元矩阵分解的非传统阻尼体系阻尼矩阵识别方法，这种方法所需要的实测振型可以是复振型，不需要质量归一化，能更好的适应目前只有有限的低阶模态能够相对准确测得的现状，与基于输出响应的模态参数识别技术结合，可以进行环境激励下的海洋平台结构非传统阻尼矩阵的识别。

Description

非传统阻尼矩阵识别方法

技术领域

本发明涉及对海洋工程结构物的非传统阻尼矩阵识别方法，特别涉及一种仅利用少数低阶实测复模态的非传统阻尼矩阵识别方法。

背景技术

海洋平台等大型工程结构在其服役期间都会表现出某种程度的能量耗散，只要涉及结构的能量耗散，阻尼在结构的响应过程中就起着至关重要的作用。传统意义上讲，这种能量耗散归因于结构固有的阻尼，材料、连接与非结构性构件的内部摩擦是影响结构阻尼的重要因素。为便于数学处理，通常假定阻尼矩阵是质量矩阵与刚度矩阵的线性组合，即常用的“比例阻尼”模型，然而这种模型将导致系统的频率和振型皆为实数。实际上，真实的振动系统并不一定表现出这一性质。非传统阻尼模型才是更普遍的现象，较比例阻尼模型更符合工程问题的本质特性。该模型的特点是阻尼矩阵无法被结构的无阻尼振型解耦，在对结构进行动力学分析的时候会导致复振型。目前工程上比较广泛接受的一种非传统阻尼模型是1992年由美国国家地震研究中心的Liang博士等人提出的。Liang等人认为系统的阻尼矩阵可以表达为一个对角阵和比例阻尼矩阵之和，然而对角阵选取的合理与否直接影响度该模型的精度。

Lancaster(1961)提出了利用结构的特征值和特征向量直接计算结构的质量矩阵、刚度矩阵及阻尼矩阵的方法。这种方法只适用于粘滞阻尼体系，并且需要一种特别的关于质量矩阵与阻尼矩阵的归一方法。就如Lancaster说的那样，该方法的不足之处在于需要对振型归一化，但是现有的测试方法还无法获得归一化后的振型，从而限制了该方法的进一步应用。

Hasselman(1972)研究了两种阻尼形式，即比例阻尼与非比例阻尼形式。该方法通过确定加速度信号的一致部分与积分部分的相位差，进而构造阻尼矩阵的非对角线元素。但是，该方法只有在测试数据无噪声时才能成立。

Beliveau(1976)利用结构的固有频率、阻尼比、振型及相位角进行阻尼矩阵的识别。该方法用到了贝叶斯与牛顿拉普拉斯理论，通过逐步迭代获得阻尼矩阵。对每一个特征向量，该方法要求解n阶线性方程，从而降低了该方法的有效性。同时，该方法虽然可考虑参数的不确定性，但并不意味着能获得比较理想的结果。

Ibrahim(1983)假定已知结构的有限元模型以及实测复模态。在此基础上将振型进行质量归一化。该方法可应用于较多自由度的体系，但是该方法只能对实测位置处单元的阻尼矩阵进行识别。

Fabunmi，Chang与Vorwald(1988)提出了在频域内利用强迫振动信息进行阻尼矩阵识别的技术。其前提是已知结构的质量矩阵、刚度矩阵以及频率响应。虽然应用这种方法获得的阻尼矩阵能够在一定程度上重现实测数据，但是往往很难获得与真实阻尼矩阵完全一致的阻尼估计矩阵。

Minas与Inman(1991)假定质量矩阵与刚度矩阵从有限元模型获得，特征值与特征向量由试验测试获得。求解阻尼体系的特征方程重可获得结构的阻尼矩阵。但是，该方法仅限于对称的正定阻尼矩阵形式。

Chen，Ju与Tsuei(1996)用频域信息估计结构的阻尼矩阵。该方法虽然可考虑噪声影响，但必须在多个频率处进行求解。

Gaylard(1996)用质量做为权函数提出了阻尼矩阵识别的确定性方法。该方法是一种时域方法，计算工作量大，同时要用到卷积。他们的算例表明，在阻尼矩阵识别的迭代过程中，用识别后的质量矩阵代替真实的质量矩阵能显著降低阻尼矩阵识别时的鲁棒性。从另一个算例也可以看出，该方法只有在瑞利阻尼时才能获得较好的结果。

Srikantha(2007)对现有频域阻尼识别方法进行了总结，将现有方法分为三类：矩阵方法，模态方法及改进方法，并用四个不同的算例进行对文中所列方法的优缺点进行了对比。

George(2009)在前期研究的基础上将模态阻尼识别模型从传统阻尼形式推广到了非传统阻尼形式。但他们的模型针对的是线性框架结构，而且其识别结果为模态阻尼比。

上述方法一方面需要较高阶的实测振型，另一方面对实测振型的归一化、变量矩阵的选择也提出了较高的要求，而这些条件又是实际工程应用中无法完全提供的，从而限制了进一步工程应用。

发明内容

本发明提出一种仅利用少数低阶实测复模态的非传统阻尼矩阵识别方法，该方法所涉及的变量矩阵形式、数值具有通用性，并且实测模态无需质量归一化处理。

为解决上述技术问题，本发明非传统阻尼矩阵识别方法，其特征在于，包括如下步骤：

A、建立结构有限元数值模型，获取拟识别结构的质量矩阵M_t与刚度矩阵K_t，以结构的阻尼矩阵C_t为待识别项；

B、选取标准基向量ω_i，其中i代表结构的模态阶次，并将质量矩阵M_t、刚度矩阵K_t及标准基向量ω_i数据存储入专用存储器中；

C、利用传感器获取结构动力响应信号；

D、利用模态参数识别技术得到模态参数λ_j、Φ_j，同时将其存储入所述专用存储器中，其中，λ_j为识拟别结构的第j阶特征值，Φ_j为拟识别结构的第j阶特征向量；

E、确定初始单元矩阵分布形态及初始值C_e0；

F、阻尼矩阵识别：

F1、从所述专用存储器中读取上述步骤B、D中存储的数据，即ω_i、M_t、K_t、λ_j、Φ_j；

F2、构造线性方程组，将C_t引入方程组中，通过求解该方程组获得新的估计阻尼矩阵

使

尽可能的逼近C_t；

所述F2步骤中包括如下步骤：

F21、结构振动信号及模态参数提取

利用传感器测试实际结构，并提取模态频率f、模态振型Φ及模态阻尼系数ξ，其中

F22、构造线性方程组： ${\mathrm{\ω}}_i^T{K}_t{\mathrm{\Φ}}_j+{\mathrm{\λ}}_j{\mathrm{\ω}}_i^T{C}_t{\mathrm{\Φ}}_j+{\mathrm{\λ}}_j^2{\mathrm{\ω}}_i^T{M}_t{\mathrm{\Φ}}_j=0;$

F23、建立非传统阻尼矩阵C_t与初始阻尼矩阵C₀的关系，

$C_0=\underset{i=1}{\overset{\mathrm{Ne}}{\mathrm{\Σ}}}{C}_j$

则

$C_t=C_0+\underset{l=1}{\overset{\mathrm{Ne}}{\mathrm{\Σ}}}\underset{k=1}{\overset{N_C}{\mathrm{\Σ}}}{\mathrm{\α}}_p{C}_{l,k}$

其中Ne为系统的单元数量；α_p为修正系数；并且p＝1，2，…，Ne×N_C，N_C为子矩阵C_l，k的个数；C_l，k为第l个单元的第k个在整体坐标系下的子矩阵；

F23、根据上述F21及F22步骤获取估计阻尼矩阵

$\hat{C}=C_0+\underset{l=1}{\overset{\mathrm{Ne}}{\mathrm{\Σ}}}\underset{k=1}{\overset{N_C}{\mathrm{\Σ}}}{\mathrm{\α}}_p{C}_{l,k}$

其中α_p为

中第p个元素。

进一步地，所述E步骤中，当C_e0的分布形态完全未知，三维结构可假定为12×12的无非零元素的矩阵，二维、一维结构分别假定为12×12、6×6的无非零元素的矩阵；当C_e0的分布形态已知或部分已知时，只需关心非零元素。

进一步地，所述E步骤中，确定C_e0初始值时，当有参考数据时取参考值，当无参考数据时可假定初始值为常数。

进一步地，所述B步骤中i＜＝结构的总自由度数目。

与现有技术相比，本发明的优点和积极效果是：

1)本发明可识别非传统阻尼形式的阻尼矩阵，较现有技术而言，所需实测模态阶数更少，并且实测模态可以是复数的形式，与现有某些技术的明显区别在于无需实数化的近似处理；

2)本发明以初始单元阻尼矩阵C_e0作为变量矩阵，较现有技术而言，该变量矩阵的分布形态、数值可根据实际情况变化。即使在C_e0完全未知的情况下，可选用常规的矩阵(如单位阵等)来替代，具有通用性；

3)本发明亦可识别结构的阻尼系数，与现有技术相比，本发明可以实现阻尼系数的跳跃式识别，如当只测得结构的第三阶模态而无第一、第二阶模态时，应用本发明可正确识别出实测的结构第三阶阻尼系数。同时，随着实测模态阶数的增多，更高阶的未测出的阻尼系数其识别精度也会提高；

总之，对于如海洋平台等大型结构，借助动力测试数据获得该结构的所有模态是不可能实现的，尤其对于环境激励下结构的动力响应测试问题，而本发明中所需要的数据为无需质量归一化的少数低阶复模态，因此，可以进行环境激励下的海洋平台结构非传统阻尼矩阵的识别，具有实际应用价值。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作一简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1：本发明的导管架式海洋平台结构有限元模型图；

图2：导管架式海洋平台结构的真实阻尼矩阵；

图3：利用前2阶实测模态识别的阻尼矩阵与真实阻尼矩阵对比图；

图4：利用前3阶实测模态识别的阻尼矩阵与真实阻尼矩阵对比图；

表1：导管架式海洋平台结构的前10阶阻尼比；

表2：利用前2阶实测模态识别的前10阶识别阻尼比与真实阻尼比；

表3：利用前3阶实测模态识别的前10阶识别阻尼比与真实阻尼比。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明的创新之处在本算例中的体现为：1)平台结构单元阻尼矩阵分布形态完全未知；2)只有低阶模态信息测得，本例取前三阶；3)实测振型无需归一化、配对等问题，无需实数化近似处理。

下面结合附图和具体实施方式对本发明作进一步详细的说明。

一，具体算法

a)建立结构有限元数值模型，获取拟识别结构的质量矩阵M_t与刚度矩阵K_t，以结构的阻尼矩阵C_t为待识别项；

b)选取标准基向量ω_i，i代表结构的模态阶次，并且i＜＝结构的总自由度数目，如对于一个5自由度的结构， ${\mathrm{\ω}}_3=\left\{\begin{array}{c}0\\ 0\\ 1\\ 0\\ 0\end{array}\right\};$ 然后将质量矩阵M_t、刚度矩阵K_t及标准基向量ω_i数据存储入专用存储器中。这样处理的优势在于：标准基向量ω_i来源于数学概念，形式、数值相对固定，后期的方程构造更具普适性；

c)获取结构动力响应信号，其或是加速度，和/或是速度，和/或是位移。

d)利用模态参数识别技术得到模态参数λ_j、Φ_j，其中λ_j、Φ_j均为复数形式，无需进行简化处理，同时将其存储入所述专用存储器中；

e)变量矩阵选取：以初始单元阻尼矩阵C_e0作为变量矩阵。

①C_e0分布形态确定

当C_e0的分布形态完全未知，三维结构可假定为12×12的无非零元素的矩阵，二维、一维结构分别假定为12×12、6×6的无非零元素的矩阵；当C_e0的分布形态已知或部分已知时，只需关心非零元素，相对于分布形态完全未知的情况，需修正的未知数数目得到降低，但可看做分布形态完全未知情况的特例。

②C_e0初始值的确定

当有参考数据时取参考值，该参考值可来源于传统阻尼模型；当无参考数据时可假定初始值为常数(如取为1)。

这样处理的优势在于：①C_e0的分布形态可根据实际情况变化，能够更好的反映结构的实际阻尼分布特性；②变量矩阵的选取更具通用性，可应用现有研究资料提供一个初始单元阻尼矩阵C_e0，即使无此条件时亦可采用常数阵来替代，工程应用前景好。

f)阻尼矩阵识别：

①从专用存储器中读取上述步骤a)、b)中存储的数据，即M_t、K_t、ω_i；

②结构振动信号及模态参数提取

借助传感器测试实际结构，利用模态参数识别方法提取模态频率f、模态振型Φ及模态阻尼系数ξ，其中

③构造线性方程组： ${\mathrm{\ω}}_i^T{K}_t{\mathrm{\Φ}}_j+{\mathrm{\λ}}_j{\mathrm{\ω}}_i^T{C}_t{\mathrm{\Φ}}_j+{\mathrm{\λ}}_j^2{\mathrm{\ω}}_i^T{M}_t{\mathrm{\Φ}}_j=0,$ 假设可构造Nm个，其中N_m＝N_i×N_j，并且N_j为所取标准基向量的个数，N_i为实测模态的阶数。

④建立非传统阻尼矩阵C_t与初始阻尼矩阵C₀的关系，即非传统阻尼矩阵C_t认为是对初始阻尼矩阵C₀修正的结果，

$C_0=\underset{i=1}{\overset{\mathrm{Ne}}{\mathrm{\Σ}}}{C}_j$

则

$C_t=C_0+\underset{l=1}{\overset{\mathrm{Ne}}{\mathrm{\Σ}}}\underset{k=1}{\overset{N_C}{\mathrm{\Σ}}}{\mathrm{\α}}_p{C}_{l,k}$

其中Ne为系统的单元数量；α_p为修正系数；并且p＝1，2，…，Ne×N_C，N_C为子矩阵C_l，k的个数；C_l，k为第l个单元的第k个在整体坐标系下的子矩阵。

⑤将步骤④带入步骤③，得

其中

与

⑥步骤⑤进一步整理为

其中

⑦步骤⑥写成矩阵的形式

其中

为N_m×(Ne×n)复矩阵；α为Ne×n的列向量；为Nm维的列向量。

⑧复数矩阵转变为实数矩阵，写成统一的形式：

分别用Re(z)与Im(z)记为复数z的实部和虚部，则可将步骤⑥转变为

Gc＝d

其中

且

c＝{α}

⑨G矩阵的分解

$G={\mathrm{U\ΣV}}^T=U_p{\mathrm{\Σ}}_p{V}_p^T$

其中，U为维数2N_m×2N_m的矩阵并且满足U^TU＝I；v为维数N_C×N_C的矩阵并且满足V^TV＝I；∑为维数2N_m×N_C的对角矩阵并且其对角线元素为奇异值。这些奇异值通常是按降序的顺序排列的，并且这些奇异值中有时有值为0的元素。并且

$\mathrm{\Σ}=\left[\begin{array}{cc}{\mathrm{\Σ}}_P& 0\\ 0& 0\end{array}\right]$

其中，∑_p为p×p的对角阵。式中，U_p与V_p分别为矩阵U与V的前p列。

修正系数的求解

其中，

为广义逆矩阵。

⑩估计阻尼矩阵

的获取

$\hat{C}=C_0+\underset{l=1}{\overset{\mathrm{Ne}}{\mathrm{\Σ}}}\underset{k=1}{\overset{N_C}{\mathrm{\Σ}}}{\mathrm{\α}}_p{C}_{l,k}$

其中α_p为中第p个元素。

二、三维海洋平台应用实施例

1、以四腿导管架平台模型进行研究，见图1。该模型共有40个管单元组成，其中4根平台腿离散为20个单元。所用材料的杨氏模量为2.1×10¹¹Pa，泊松比为0.3，密度为7860Kg/m³，即单位长度的质量为9.825Kg/m。

为构造非传统阻尼矩阵，首先假定该平台模型C_n＝10^-5K_n，然后令真实阻尼矩阵，即非传统阻尼矩阵为C_t＝γ_nC_n，并且参数γ_n服从均值为0方差为4.5的高斯分布。真实阻尼矩阵C_t见图2，导管架式海洋平台结构的前10阶阻尼比见表1。

表1导管架式海洋平台结构的前10阶阻尼比

阶次	真实阻尼比
		1	0.00042895
2	0.00040097
		3	0.00055237
4	0.0007118
		5	0.0008844
6	0.00082969
		7	0.00101
8	0.0010351
		9	0.0018523
10	0.0026275

2、初始阻尼矩阵分布形态及初始值

因初始单元阻尼矩阵未知，假定其具有与单元刚度矩阵一致的分布形态，又因其初始值也未知，故假定各元素均为常数1。因此，各单元初始阻尼矩阵为

$c_j=\left[\begin{array}{cccccccccccc}1& 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 1& 0& 1& 0& 0& 0& 1\\ 0& 0& 1& 0& 1& 0& 0& 0& 1& 0& 1& 0\\ 0& 0& 0& 1& 0& 0& 0& 0& 0& 1& 0& 0\\ 0& 0& 1& 0& 1& 0& 0& 0& 1& 0& 1& 0\\ 0& 1& 0& 0& 0& 1& 0& 1& 0& 0& 0& 1\\ 1& 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 1& 0& 1& 0& 0& 0& 1\\ 0& 0& 1& 0& 1& 0& 0& 0& 1& 0& 1& 0\\ 0& 0& 0& 1& 0& 0& 0& 0& 0& 1& 0& 0\\ 0& 0& 1& 0& 1& 0& 0& 0& 1& 0& 1& 0\\ 0& 1& 0& 0& 0& 1& 0& 1& 0& 0& 0& 1\end{array}\right]$

3、标准基向量的选取

平台模型共离散为40个单元，即Ne＝40；总自由度数为120，理论上标准基向量可取120个，具体数目可根据需修正未知数的个数确定，本例取120个。

4、阻尼矩阵识别

假定只有前两阶模态能够测得，即Ni＝2，识别后的阻尼矩阵与真实阻尼矩阵对比见图3，阻尼比的对比情况见表2。

表2利用前2阶实测模态识别的前10阶识别阻尼比与真实阻尼比

阶次	真实阻尼比	估算阻尼比
			1	0.00042895	0.00042895
2	0.00040097	0.00040097
			3	0.00055237	0.00050104
4	0.0007118	0.000536
			5	0.0008844	-0.00094747
6	0.00082969	0.0011229
			7	0.00101	0.00068559

8	0.0010351	-0.00080283
			9	0.0018523	0.0045632
10	0.0026275	0.0063219

表2说明在初始单元阻尼矩阵未知且只有前2阶模态信息时，应用本发明后实测模态的前两阶阻尼比均可得到很好的估计。图3表明，估算的阻尼矩阵

与真实阻尼矩阵C_t在分布形态上是一致的，保证了真实阻尼矩阵的非传统阻尼分布特性，同时

在数值上也比较接近C_t，证明非传统阻尼矩阵的识别结果也比较好。

当第3阶模态亦可测得时，即前三阶模态均可测得时，本工况主要研究随着模态阶次的增加，本发明的识别精度是否能够得到提高。应用本发明，识别后的阻尼矩阵与真实阻尼矩阵对比见图4，阻尼比的对比情况见表3。

表3利用前3阶实测模态识别的前10阶识别阻尼比与真实阻尼比

阶次	真实阻尼比	估算阻尼比
			1	0.00042895	0.00042895
2	0.00040097	0.00040097
			3	0.00055237	0.00055237
4	0.0007118	0.00071183
			5	0.0008844	0.0008836
6	0.00082969	0.00082943
			7	0.00101	0.0010094
8	0.0010351	0.0010344
			9	0.0018523	0.0018458
10	0.0026275	0.0026218

表3表明，本发明不仅可以准确估计实测模态阻尼系数，同时，随着实测模态阶次的增加，未测出的模态阻尼系数其识别精度也会提高。图4也说明非传统阻尼矩阵整体识别精度得到提高。

总之，算例充分证明本发明在只有少数低阶实测复模态时即可比较准确的识别出结构的非传统阻尼分布特性，并且能够以矩阵的形式表示。

最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。

Claims (4)

1.非传统阻尼矩阵识别方法，其特征在于，包括如下步骤：

A、建立结构有限元数值模型，获取拟识别结构的质量矩阵M_t与刚度矩阵K_t，以结构的阻尼矩阵C_t为待识别项；

B、选取标准基向量ω_i，其中i代表结构的模态阶次，并将质量矩阵M_t、刚度矩阵K_t及标准基向量ω_i数据存储入专用存储器中；

C、利用传感器获取结构动力响应信号；

D、利用模态参数识别技术得到模态参数λ_j、Φ_j，同时将其存储入所述专用存储器中，其中，λ_j为识拟别结构的第j阶特征值，Φ_j为拟识别结构的第j阶特征向量；

E、确定初始单元矩阵分布形态及初始值C_e0；

F、阻尼矩阵识别：

F1、从所述专用存储器中读取上述步骤B、D中存储的数据，即ω_i、M_t、K_t、λ_j、Φ_j；

F2、构造线性方程组，将C_t引入方程组中，通过求解该方程组获得新的估计阻尼矩阵 

使 尽可能的逼近C_t；

所述F2步骤中包括如下步骤：

F21、结构振动信号及模态参数提取

利用传感器测试实际结构，并提取模态频率f、模态振型Φ及模态阻尼系数ξ，其中 

F22、构造线性方程组：

F23、建立非传统阻尼矩阵C_t与初始阻尼矩阵C₀的关系， 

则

其中Ne为系统的单元数量；α_p为修正系数；并且p＝1，2，…，Ne×N_C，N_C为子矩阵

的个数；

为第                                                  个单元的第k个在整体坐标系下的子矩阵；

F23、根据上述F21及F22步骤获取估计阻尼矩阵 

其中α_p为 中第p个元素。

2.根据权利要求1所述的非传统阻尼矩阵识别方法，其特征在于，所述E步骤中，当C_e0的分布形态完全未知，三维结构可假定为12×12的无非零元素的矩阵，二维、一维结构分别假定为12×12、6×6的无非零元素的矩阵；当C_e0的分布形态已知或部分已知时，只需关心非零元素。

3.根据权利要求2所述的非传统阻尼矩阵识别方法，其特征在于，所述E步骤中，确定C_e0初始值时，当有参考数据时取参考值，当无参考数据时可假定初始值为常数。

4.根据权利要求3所述的非传统阻尼矩阵识别方法，其特征在于，所述B步骤中i＜＝结构的总自由度数目。 

Images