基于扫频响应反推黏弹性材料分数阶导数模型

孙 伟， 王 茁， 朱明伟

(东北大学机械工程与自动化学院 沈阳，110819)

引 言

利用黏弹性材料耗能实现对结构件的振动抑制已成为一种重要的被动阻尼减振方法。创建可以有效预估黏弹性复合结构动力学特性的模型是该阻尼减振技术的基础，要获得分析模型，必须首先知道黏弹性材料的力学特性参数。黏弹性阻尼材料的力学特性参数包括储能模量和损耗因子(或耗能模量)，具有频率依赖性的特点，这是黏弹性阻尼材料最独特的力学特点。黏弹性材料的频率依赖性使有效辨识材料的力学特性参数成为一项具有挑战性的研究任务。

目前，研究人员已提出大量的模型来表征黏弹性材料的频率依赖性，典型的模型包括：复模量模型[1]、Golla-Hughes-McTavish模型[2]、滞弹性位移场模型[3]和分数阶导数模型[4]等。这使辨识黏弹性材料参数变为获得上述表征模型的各个特性参数。在上述模型中，分数阶导数模型可以有效表征真实黏弹性材料在时域和频域的力学行为，因而得到了广泛应用[5-7]。

将数值分析和振动实验相结合的反推法是辨识黏弹性材料力学特性参数的一种常用方法，很多学者针对黏弹性复合梁结构进行了反推辨识。Barkanov等[8]基于响应面优化技术，通过实测复合梁的振动衰减响应来辨识用复模量模型表征的黏弹性材料参数。Shi等[9]通过实测复合梁的共振频率和模态损耗因子，利用数值优化法反推黏弹性材料的弹性模量和损耗因子。上述研究最大的缺陷是仅能针对共振频率处的数据点进行参数辨识，而由于数据点数量上的欠缺可能导致辨识的黏弹性参数与真实值存在偏差。文献[10]针对一个自由阻尼层悬臂梁，用力锤激励测得频响函数，采用两步法反推出用分数阶导数模型表征的黏弹性材料参数。该方法已经不局限于仅针对共振频率点进行参数辨识，是黏弹性材料力学特性参数辨识方法上的一个巨大进步。但是，锤击测试涉及激励位置难以准确确定的问题，不精确的激励点可能影响黏弹性材料参数的辨识结果。另外，上述研究在创建的理论模型中只考虑了单一的黏弹性材料阻尼，没有考虑梁夹持端以及空气中存在的阻尼(这里定义为剩余等效黏性阻尼)。这样的模型将导致测量的阻尼值总是大于计算的阻尼值，将影响黏弹性材料参数的辨识结果。

笔者提出一种基于扫频响应来反推黏弹性材料分数阶导数模型参数的方法。该方法具有以下特点：a. 以具有自由阻尼层的悬臂板形试件为对象开展黏弹性材料参数辨识研究；b. 理论模型中同时考虑了材料阻尼和剩余等效黏性阻尼；c. 利用共振点及非共振点处的数据实施分数阶导数模型参数辨识；d. 利用基础激励下扫频得到的频域响应进行参数辨识。研究中建立了基础激励作用下黏弹性复合板有限元方程，给出了求解振动响应的方法。提出面向振动响应、基于灵敏度法的匹配计算来实现黏弹性参数的反推辨识。最后，以贴敷ZN-1型黏弹性材料的悬臂钛板为对象，实现了该黏弹性材料用分数阶导数模型描述的力学特性参数辨识。

1 反推法辨识原理

如图1所示，基于反推法的黏弹性材料参数辨识主要涉及3项技术，分别为振动测试、理论建模及匹配计算。待辨识的黏弹性材料分数阶导数模型可表示为

(1)

图1 基于反推法辨识黏弹性材料分数阶导数模型的原理Fig.1 The procedure of identifying the fractional derivative model of viscoelastic materials

对于不同的聚合物组合以及承受的力学状态的差异，通常黏弹性材料表现出不同的力学特点[10]。这里选择贴敷黏弹性材料的板型试件来辨识黏弹性材料参数，是因为基于梁型试件辨识的结果不一定完全满足板型试件动力学分析的需要。

辨识原理中的振动测试包含两部分：a.获得阻尼处理前悬臂板的各阶模态阻尼比；b.获得基础激励作用下黏弹性阻尼板的频域扫频响应。其中，模态阻尼比将作为剩余等效黏性阻尼引入分析模型，而测得的频域扫频响应主要用于反推辨识黏弹性材料参数。理论分析主要用于求解黏弹性阻尼板的频域振动响应，笔者创建的模型中同时考虑了黏弹性材料阻尼以及剩余等效黏弹性阻尼。

匹配计算是为了使理论模型获得的振动响应值与实验值的偏差最小，进而反推出所要获得的分数阶导数模型参数。这里提出了以振动响应为目标函数、基于灵敏度的匹配计算方法。需要说明的是，匹配计算可同时基于共振及非共振点进行材料参数辨识，但是考虑到共振点信噪比较高，因而应优先选择共振点。此外，选择的共振及非共振频率点总数应大于分数阶导数模型待辨识参数的个数，对于本实例选择的频率点数应大于4。

2 黏弹性复合板振动响应分析

2.1 运动方程的建立及响应求解

如图2所示，当薄板产生弯曲振动时，自由阻尼层会随基体一起运动，在阻尼层内部将产生拉压变形而产生耗能，从而起到减振的作用。笔者认为对于这样的结构，层间剪切效应可以被忽略，而将黏弹性阻尼层与基体视为一体。事实上，文献[10-11]也是按照这种忽略层间效应来对黏弹性自由阻尼层复合结构进行建模的。

设坐标轴恰好位于复合板中性面上，Hv为黏弹性阻尼层厚度，He为金属基体厚度，Hev为基体及黏弹性层中面之间的距离，d为金属基体中面到复合结构中性面的距离。

图2 基础激励下的黏弹性复合悬臂板Fig.2 The cantilever plate attached with viscoelastic FLD under base excitation

对式(1)的分数阶导数模型进行实部及虚部分解，可得到复模量模型，其表达式为

(2)

用此复模量模型进行黏弹性阻尼板的振动特性分析。由于在辨识算法中仅涉及若干频率点对应的振动响应值，频率信息已包含在内，因而在响应的求解过程中不需要考虑黏弹性材料的频率依赖性。

对于处于基础激励作用下的该黏弹性复合板，其频域运动方程为

(-ω2M+iωC+K*)X=F

(3)

其中：M为复合板的质量矩阵；C为由夹持边界及空气产生的剩余等效黏性阻尼矩阵；K*为基体及黏弹性阻尼层共同产生的复刚度矩阵(包含了材料阻尼)；X和F分别为相对响应向量和激振力向量。

X和F表达式为

X=U0-GU

(4a)

(4b)

这里采用模态叠加法求解黏弹性复合板在基础激励作用下的振动响应。求解实模态的特征方程可表示为

(5)

求得的各阶模态振型向量可组成模态振型矩阵φ，用模态振型矩阵φ对式(3)进行解耦，得到一组独立的以模态坐标xNr(r=1,2,…,n)表达的运动方程，即

(6)

式(6)括号内第1项和第3项系数分别为

其中：ηr为第r阶模态损耗因子，来自于黏弹性材料阻尼。

根据振型阻尼假设，式(6)中括号内第2项可以写为

(8)

其中：ξr为第r阶剩余等效模态阻尼比，可按照结构阻尼处理前的各阶模态阻尼比引入。

将式(7)，(8)代入式(6)，可以将方程转换到模态坐标下，整理后得到

(9)

由式(9)可知，在响应的求解中同时考虑了两种阻尼：a. 由黏弹性材料产生的模态损耗因子ηr；b. 由剩余等效黏性阻尼产生的模态阻尼比ξr。

(10)

结构在基础激励作用下振动响应为

(11)

其中：|·|为求模运算。

实际计算时一般不用取所有n阶模态，考虑的阶次m只需大于分析频率范围内结构的最高阶次即可。

2.2 刚度及质量矩阵的生成

对于图2所示的黏弹性自由阻尼复合板，选用四边形板单元(该单元共有4个节点，每个节点有3个自由度)进行模拟。按照以下流程求解刚度矩阵和质量矩阵。

1) 将薄板复合结构划分为若干个单元。

2) 求解单元的刚度及质量矩阵，求解式分别为

De和Dv分别为基体和黏弹性阻尼层的弹性矩阵，表达式为

(13a)

(13b)

其中：μe和μv分别为弹性基体和黏弹性层的泊松比。

复合板的中性面到基体中心的距离d表示为

d=EvHvHev/(EeHe+EvHv)

(14)

3) 在考虑约束条件的基础上组合单元刚度矩阵及单元质量矩阵，最终可形成黏弹性复合板的总刚度矩阵K*和总质量矩阵M。

3 基于灵敏度的匹配计算

参照式(1)，将待辨识的分数阶导数模型的4个参数E0，E∞，τ和α作为设计变量，采用差分法获得灵敏度。以E0为例，第i个频率点振动响应对设计变量E0的灵敏度SE0,i可表示为

(15)

其中：ΔE0为设计变量的变化量(即步长)，这里设为10-4E0；ΔxE0,i为设计变量变化ΔE0时振动响应的变化量(残差)。

同样，分别求出振动响应对E∞，τ和α的灵敏度SE∞,i，Sτ,i和Sα,i。假设在整个辨识中共选择P个频率点，则求得的针对设计变量的灵敏度可组成灵敏度矩阵S，即

(16)

该匹配计算的目标函数为

(17)

其中：Wx，Wb分别为振动响应和设计变量的权重矩阵；Δb为设计变量变化量所组成的向量；Δx为对应各频率点振动响应的残差组成的向量。

Δb和Δx的表达式为

(18a)

(18b)

需要说明的是，式(18)中的残差对应所有设计变量变化时响应的变化量。

得到如下迭代公式

(19)

经过若干次迭代，满足收敛条件后迭代终止，获得对应分数阶导数模型待辨识的4个参数。需要说明的是，考虑到共振点对应的响应值信噪比高，设置权重时应将共振点的权重设为大值，将非共振点权重设为小值。

4 实例研究

以贴敷ZN-1型黏弹性材料的悬臂钛板为例，基于扫频获得频域响应反推分数阶导数模型参数。

4.1 试件及振动特性测试

阻尼处理前后的钛板如图3所示。该薄板的长、宽分别为154.1 mm，109.6 mm，钛板基体的厚度为1.32 m，薄板夹持区长度为20 mm，拧紧力矩为32 N·m。在薄板的一侧贴敷ZN-1型黏弹性材料，厚度为0.5 mm。钛板的材料参数是已知的：储能模量为110.32 GPa，密度为4 420 kg/m3，损耗因子为0.000 7。ZN-1型黏弹性材料的密度可实测获得，具体为789.5 kg/m3，黏弹性材料及钛板的泊松比分别取0.49和0.3，ZN-1型黏弹性材料的分数阶导数模型参数需要辨识获得。

实验现场如图4所示，涉及的仪器设备如表1所示。整个测试过程中均用激光测振仪拾振，参照图4中的坐标，拾振点位置为x=11 mm，y=67 mm。在基于反推法的辨识中，实验系统与理论模型相对应是至关重要的。剔除夹持区后，实验结构为受振动台基础激励的、单面全部贴敷ZN-1型阻尼材料的和由夹具固定的(实践表明，32 N·m的螺栓拧紧力矩可以保证板的底部为完全固定约束)悬臂结构系统，这就保证了实验系统与理论模型之间激励条件、结构形式和边界条件均具有一致性。

图3 钛板实验件Fig.3 Specimens used in the experiment

图4 测试现场图Fig.4 Actual photograph of vibration test

表1 主要测试仪器Tab.1 The instruments used in this test

对阻尼处理前的薄板进行锤击测试，采用半功率带宽法由获得的频响函数辨识出各阶模态阻尼比，测试结果如表2所示。利用振动台，设定包含共振频率的扫频区间，对阻尼处理后的薄板进行分段扫频测试。选择第1，3，5阶共振区的频域响应用于黏弹性材料参数辨识，对应的激励幅度分别为0.25，1和2 g，相应的测试结果如图5所示。

表2阻尼处理前薄板的模态阻尼比

Tab.2Modaldampingratiosofthetitaniumplatewithoutdampingtreatment%

4.2 黏弹性材料的力学特性参数辨识

参照图5所示的扫频响应，选择3个共振频率点以及每个共振频率左右各一个非共振频率点(即3个共振点，6个非共振点)，利用本研究方法进行参数辨识。

图5 实验测试及理论计算得到的频域响应Fig.5 Frequency domain responses obtained by experiment and analysis

整个模型共划分了100个单元，121个节点，利用Matlab编程完成整个计算。参考文献[13]选择ZN-1型黏弹性材料分数阶导数模型的初值，具体为：E0=1.8 MPa，E∞=12.5 MPa，τ=2.42×10-4，α=0.99。计算各频率点处振动响应对分数阶导数模型参数的灵敏度以及振动响应的残差向量。接着，将分数阶导数模型各参数的权重均取0.25，在振动响应的权重矩阵中，共振点处的权重取0.2，非共振点处的权重取0.06。最后，将确定的灵敏度矩阵、振动响应残差向量以及权重矩阵代入式(17)进行迭代计算。设收敛条件为：计算得到的振动响应与相应的实验值偏差小于1%。经过若干次迭代，可反推获得ZN-1型黏弹性材料的分数阶导数模型参数，具体结果如表3所示。例如，以激振频率为横轴， ZN-1型黏弹性材料的储能模量或损耗因子为纵轴，则可对该材料进行频率依赖性频域表征，如图6所示。

表3ZN-1型黏弹性材料的分数阶导数模型参数

Tab.3FractionalderivativemodelparametersofZN-1viscoelasticmaterial

参数E0E∞τα辨识值1．5MPa11．3MPa1．02×10-40．97

图6 ZN-1型黏弹性材料的频率依赖性表征Fig.6 Characterization of ZN-1 viscoelastic material with frequency dependent characteristics

为了说明辨识结果的合理性，将获得的黏弹性材料分数阶导数模型写成复模量的形式并代入到黏弹性阻尼板振动响应分析模型中。黏弹性复合板的频域振动响应需按照各频率点进行求解(计算步长为1 Hz)。首先，针对每个频率值确定黏弹性材料的储能模量和损耗因子；其次，确定总刚度矩阵K*和总质量矩阵M，并进行模态分析；然后，取8阶模态，将表2中的阻尼值作为剩余等效黏性阻尼输入分析模型，按照模态叠加法求解每个频率点对应的响应值；最后，依次类推求解各频率点对应的响应值，绘制与实验频率范围一致的频域响应曲线，如图5所示。提取各指定频率点对应的理论计算与实验测试获得的振动响应值并进行比对，如表4所示。可以看出，理论计算获得的黏弹性阻尼板频域响应值与实测值几乎一致，证明了辨识结果的合理性。

表4实验测试及理论计算获得的振动响应

Tab.4Vibrationresponsesobtainedbyexperimentandanalyticalcalculation

频率点/Hz振动响应/(m·s-1)实验测试理论计算偏差/%56．540．03090．03041．6259．13∗0．23510．23560．2162．030．02440．02471．23363．090．13080．13291．61367．11∗0．48540．48400．28372．210．12350．12591．94642．100．01930．01993．11645．08∗0．06440．06480．62649．360．01300．01353．85

*号为共振点

5 结 论

1) 在同时考虑黏弹性材料阻尼以及剩余等效黏性阻尼的基础上，创建了黏弹性阻尼板振动响应分析模型。以理论计算获得的振动响应与实测值偏差最小为目标函数，研发出基于灵敏度的匹配算法，在同时保证计算精度和计算效率的前提下实现了黏弹性材料分数阶导数模型参数的辨识。

2) 以ZN-1型黏弹性材料为例，用本研究方法辨识出该材料的分数阶导数模型参数。将获得的辨识结果代入到黏弹性阻尼板振动响应分析模型中，并比较计算得到的响应值与实测值，发现两者几乎一致，从而证明了辨识结果的合理性。利用实例中9个频率点(3个共振点，6个非共振点)辨识出的黏弹性材料参数对复合板结构进行动力学分析，完全可以再现实际的黏弹性阻尼板频域振动响应行为。

[1] Martinez-Agirre M, Elejabarrieta M J. Dynamic characterization of high damping viscoelastic materials from vibration test data[J]. Journal of Sound and Vibration, 2011, 330(16): 3930-3943.

[2] McTavish D J, Hughes P C. Modeling of linear viscoelastic space structures[J]. Journal of Vibration and Acoustics, 1993, 115(1): 103-110.

[3] Lesieutre G A, Bianchini E. Time domain modeling of linear viscoelasticity using anelastic displacement fields[J]. Journal of Vibration and Acoustics, 1995, 117(4): 424-430.

[4] Bagley R L, Torvik P J. Fractional calculus in the transient analysis of viscoelastically damped structures[J]. AIAA Journal, 1985, 23(6): 918-925.

[5] Escobedo-Torres J, Ricles J M. The fractional order elastic-viscoelastic equations of motion: formulation and solution methods[J]. Journal of Intelligent Material Systems and Structures, 1998, 9(7): 489-502.

[6] Schmidt A, Gaul L. Finite element formulation of viscoelastic constitutive equations using fractional time derivatives[J]. Nonlinear Dynamics, 2002, 29(1-4): 37-55.

[7] Cortés F, Elejabarrieta M J. Homogenised finite element for transient dynamic analysis of unconstrained layer damping beams involving fractional derivative models[J]. Computational Mechanics, 2007, 40(2): 313-324.

[8] Barkanov E, Skukis E, Petitjean B. Characterisation of viscoelastic layers in sandwich panels via an inverse technique[J]. Journal of Sound and Vibration, 2009, 327(3): 402-412.

[9] Shi Y, Sol H, Hua H. Material parameter identification of sandwich beams by an inverse method[J]. Journal of Sound and Vibration, 2006, 290(3): 1234-1255.

[10] Kim S Y, Lee D H. Identification of fractional-derivative-model parameters of viscoelastic materials from measured FRFs[J]. Journal of Sound and Vibration, 2009, 324(3): 570-586.

[11] Parthasarathy G, Reddy C V R, Ganesan N. Partial coverage of rectangular plates by unconstrained layer damping treatments[J]. Journal of Sound and Vibration, 1985, 102(2): 203-216.

[12] 刘建涛, 杜平安, 黄明镜, 等. 阻尼连续体简谐基础振动的有限元实现方法[J]. 机械工程学报, 2010, 46(1): 109-114.

Liu Jiantao, Du Ping′an, Huang Mingjing, et al. Finite element approach for analyzing continua with damp under harmonic loads from the base[J]. Journal of Mechanical Engineering, 2010, 46(1): 109-114.(in Chinese)

[13] 李恩奇, 盛秀成, 王江枫. ZN-1型黏弹性阻尼材料模型参数修正研究[J]. 强度与环境, 2011, 38(5): 36-41.

Li Enqi, Sheng Xiucheng, Wang Jiangfeng. Parameter update of ZN-1 viscoelastic damping material model[J]. Structure & Environment Engineering, 2011, 38(5): 36-41. (in Chinese)