冲击现象广泛存在于产品生产、运输和工作过程中。为了提高产品的抗冲击性能，通常采用冲击试验机在地面模拟作用在产品上的真实冲击环境，并形成考核、评定产品冲击环境适应性的一系列标准^[1-2]。半正弦波形发生器作为冲击试验中的关键部件，可以通过改变其尺寸、材料属性来获得相应的半正弦冲击脉冲^[3-4]。因此，为了高效地设计半正弦波形发生器，有必要准确建立半正弦波形发生器的动力学模型和相应的参数反辨识方法。

针对半正弦波形发生器的动力学模型，学者们做了很多探索。文献[5]考虑了半正弦波形发生器的线性刚度并将冲击试验系统简化为单自由度动力学系统进行分析，该模型简单实用，但误差很大；文献[6]采用Arruda-Boyce模型作为半正弦波形发生器的材料本构模型，并利用有限元软件进行正交仿真试验，再通过支持向量机拟合得到描述半正弦波形发生器动力学行为的代理模型。该方法考虑了半正弦波形发生器的硬度、直径、厚度等材料和结构特性，具有一定的通用性，但Arruda-Boyce模型中的相关参数是通过静态拉伸试验得到的，而半正弦波形发生器实际工作在冲击压缩环境下，并且Arruda-Boyce模型没有考虑其黏弹性特性；文献[7]考虑到橡胶材料的硬非线性，选用三次型材料模型建立了半正弦波形发生器的超弹性动力学模型，该模型简单实用，但没有给出确定模型中基本参数的方法，并且也没有考虑冲击脉冲的波形形状变化；文献[8]借鉴橡胶隔振器的静态单轴拉伸和压缩应力-应变曲线的拟合方程设计出了半正弦波形发生器的超弹性模型，该模型考虑了半正弦波形发生器的面积和厚度因素，但没有考虑形状因素，同时因为该文献主要考虑低速冲击工况，工作应变率较低，所以忽略了波形发生器的黏弹性特性；文献[9]给出了半正弦波形发生器静刚度、动刚度、冲击刚度之间的经验关系，同时还利用有限元软件开发了半正弦波形发生器冲击响应的虚拟仿真系统，但其本质仍然是基于Mooney-Rivlin模型的橡胶超弹性模型，并且试验波形与计算波形偏差较大。文献[10]采用达芬方程和Pan-Yang模型相结合的方法构建了描述半正弦波形发生器的非线性动力学模型，该模型展现了较高的波形预测精度，但模型黏弹性项中的待辨识系数缺乏物理意义，且待辨识系数较多，计算量较大，辨识效率较低。

总结发现，之前关于半正弦波形发生器动力学模型研究的主要问题有：缺乏物理意义明确的黏弹性项，没有考虑波形发生器的形状和尺寸因素，模型参数获取困难，不关注理论计算波形与试验波形的整体波形差别。本文首先基于有阻尼达芬方程^[11]建立了半正弦波形发生器的非线性动力学模型，该模型具有形式简单、每一项都有物理意义的优点。同时借鉴圆柱形橡胶隔振器的冲击刚度公式，将半正弦波形发生器的尺寸因素和形状因素考虑进模型之中。然后基于量子遗传算法(quantum genetic algorithm, QGA)，以计算所得冲击波形和试验所得冲击波形的拟合度最高为优化目标，对所建立的非线性动力学模型的相关参数进行反辨识研究。最后以某型半正弦波形发生器为例，通过在跌落式冲击试验机上试验，对上述模型和参数标定方法进行验证，试验结果表明：所建立的半正弦波形发生器的非线性动力学模型具备较高预测精度，其中峰值预测误差和脉宽预测误差均在5%以下，在跌落高度1 m以下时计算波形和试验波形的拟合误差小于15%；提出的参数反辨识方法稳定性好、收敛速度和精度高。本文的研究内容有助于精确预测和模拟半正弦冲击试验波形，也有助于根据试验要求高效合理地设计半正弦波形发生器，同时对于设计、调试、标定冲击试验机也有一定的参考价值。

1 非线性动力学模型

如图 1所示，半正弦波形发生器一般采用橡胶(或聚氨酯)材质并设计为圆柱形。针对橡胶材料本构模型，研究发现橡胶材料同时兼具超弹性和黏弹性的特点^[12]。橡胶产品受激励后的响应力为非线性弹性力与黏弹性力的叠加^[13]。

1.1 非线性弹性力建模

文献[3, 7-8]指出半正弦波形发生器在冲击激励作用下具有硬非线性特性。同时文献[14]在研究橡胶隔振器的冲击响应特性时指出：橡胶隔振器的变形速率对超弹性部分影响很小。综合上述结论，并考虑到达芬方程可以较好地描述硬非线性动力学系统，本文建立的非线性弹性力方程为：

(1)

式中, K(x)为半正弦波形发生器非线性弹性力; k为线性刚度系数; β为硬非线性刚度系数; x(t)为t时刻波形发生器在冲击工况下的变形, 简记为x。

为了引入半正弦波形发生器的材料因素、尺寸因素和形状因素, 总结文献[15]中圆柱形橡胶隔振器的静刚度经验公式, 并乘以一定的修正系数, 得

(2)

式中, α₁为修正系数, 可由试验得到; 为初始冲击速度; n为半正弦波形发生器承载面积和自由表面面积之比; H为半正弦波形发生器的邵氏硬度。

文献[18]指出:半正弦波形发生器的冲击刚度大约是其静刚度的3倍左右。因此

(3)

考虑到半正弦波形发生器的变形量x一般在毫米量级, 代入(3)式可得非线性刚度系数大约是线性刚度系数的10⁴倍, 因此可以假设非线性刚度系数为

(4)

式中, α₂为修正系数, 可由试验辨识得到。

1.2 黏弹性力建模

在半正弦波形发生器动力学特性的研究过程中, 虽然很少关注其黏弹性特性, 但针对橡胶隔振器的黏弹性特性, 却有不少研究。文献[17]指出Maxwell模型和Kelvin-Voigt模型在广义上完全等效; 文献[18]指出橡胶材料的黏弹性特性会导致其动态特性与静态特性有所不同, 主要表现为应变滞后于应力, 且随应变率的增大而显著变化; 文献[19]给出了橡胶材料在大变形、高应变率下的动态压缩性能, 但没有给出通用的计算模型; 文献[20]考虑了应变率的影响, 并给出了橡胶隔振器在振动工况下的超-黏弹性模型, 但该模型待标定参数太多(共11个), 并且个别参数的物理意义并不明确。本文综合考虑计算成本和参数物理意义, 采用线性黏弹性模型来描述半正弦波形发生器的阻尼系数, 则可以得到

(5)

式中, 为黏弹性力; c_eq为半正弦波形发生器的当量黏性阻尼, 在不同的冲击工况下会有波动, 具体可由试验辨识。为t时刻半正弦波形发生器在冲击工况下的变形速率, 下文简记为。

为保证模型的无量纲性, 参考单自由度系统弹簧质量系统阻尼比的概念^[11]将阻尼系数改写为

(6)

式中，ξ₁为无量纲的阻尼比, 对于天然橡胶, 其值在0~1之间; m是冲击物体的质量; ξ₁可以通过试验辨识得到, 此外ξ₁与α₁, α₂类似, 也起着修正系数的作用, 但其有直接的物理意义。

1.3 建立非线性动力学模型

综上, 建立计算半正弦波形发生器冲击脉冲的非线性动力学模型为

(7)

2 冲击响应测试方法

如图 2所示, 本文采用跌落式冲击试验机测试半正弦波形发生器的冲击响应。该型冲击试验机通过气压驱动将台面提升到不同的高度, 然后自由释放来获得不同的冲击速度。配合相应的采集设备, 可以直接采集到每次冲击的加速度信号, 并将其作为半正弦波形发生器的冲击响应。

3 冲击响应计算方法

为了精确计算出半正弦波形发生器的冲击响应, 需要结合上述跌落式冲击试验机的相关特性建立冲击响应动力学方程。文献[3, 21]均将该冲击响应试验等效为如图 3所示的单自由度无阻尼自由振动系统。图 3中，m为台面质量，h为台面跌落高度，k为半正弦波形发生器的线性刚度，x为t时刻台面的位移。规定台面与波形发生器接触时刻为冲击起始时刻, 分离时刻为冲击结束时刻, 因此台面位移与波形发生器的变形量一致。综上可得该系统的动力学方程为。

文献[10]指出, 上述单自由度模型忽略了诸多关键因素:①底座、阻尼器和缓冲气囊组成的减振系统的影响; ②台面自身的柔性; ③台面从跌落到与波形发生器接触过程中的速度损失; ④台面和底座自身重力对半正弦波形发生器和缓冲气囊变形量的影响; ⑤半正弦波形发生器自身刚度和阻尼的非线性特性。综合考虑上述因素, 本文建立的冲击响应动力学试验系统模型如图 4所示。

图 4中，m₁为台面质量，m₂为底座质量，h为台面跌落高度，x₁, x₂分别为t时刻台面和底座的绝对位移, 同样规定台面与波形发生器接触时刻为冲击起始时刻, 分离时刻为冲击结束时刻, 并忽略波形发生器自身质量, 从而x₁-x₂为波形发生器的变形量; 结合(1)式可得半正弦波形发生器的非线性弹性力为为波形发生器的阻尼系数; F₂(x₂)为缓冲气囊变形引起的恢复力; c₂为阻尼器的阻尼系数。借鉴文献[10]的建模方法, 建立该系统的冲击动力学方程为

(8)

式中，μ为台面从跌落到碰撞波形发生器过程中的速度损失系数, 取值范围为(0.7, 1)^{[10, 21]}。

至于台面的自身柔性因素对所使用的冲击试验机的台面进行模态分析可知该台面的固有频率在803.74 Hz以上, 因此相对于固有频率70 Hz左右的半正弦波形发生器而言可以等效为刚体。

结合相关初始条件, 可以采用龙格库塔方法求解非线性动力学方程组(8)得到半正弦波形发生器的冲击响应脉冲。采用龙格库塔方法的原因有:1)龙格库塔方法具有极强的求解常微分方程的能力; 2)龙格库塔方法可以直接计算出半正弦波形发生器的时域波形, 从而避免了转换带来的各种误差。

4 基于QGA的参数辨识方法

(8) 式中的待辨识参数有:m₁, k, β, c₁, m₂, F₂, c₂, h, μ。其中m₁, m₂和F₂的数据可以从冲击试验机的使用说明书中直接得到, 此外冲击试验机台面的提升高度h也可精确测得。当确定半正弦波形发生器之后, 可以精确测得其直径、厚度和硬度, 并代入(2)式计算, 则只须辨识出α₁即可计算出k。同理, 只须辨识出α₂即可计算出β。由(6)式可知如果已知ξ₁即可计算出c₁。借鉴(6)式, c₂也可以获得类似的形式, 即, 其中k₂=dF₂/dx₂, 为缓冲气囊的刚度。因为F₂可知, 所以只需须识出ξ₂即可计算出c₂。综上, 待辨识参数共5项:α₁, α₂, ξ₁, ξ₂, μ, 并且上述待辨识参数均是无量纲量。

4.1 构造优化模型

本文以待辨识参数为优化变量, 以预测波形和实测波形的确定系数最小为优化目标, 通过将半正弦波形发生器的冲击响应计算模型转化为优化模型来辨识模型中的未知参数。显然, 半正弦波形发生器的冲击响应是与参数α₁, α₂, ξ₁, ξ₂, μ相关的函数, 不妨记为Y(α₁, α₂, ξ₁, ξ₂, μ)。

在建立优化模型之前, 须要确定上述5项待辨识参数的约束范围。已知ξ₁的取值范围为(0, 1);μ的取值范围为(0.70, 1);ξ₂为阻尼器的阻尼比, 因为阻尼器的阻尼值较大, 时常出现过阻尼情况, 所以将其取值范围暂定为(0, 2);α₁和α₂的取值范围可以通过若干次预试验暂时确定为(0, 4)。综上, 可得如下优化模型

(9)

式中, n为冲击响应信号的维数; y_i为计算所得冲击响应信号; 为相同条件下试验所得冲击响应信号; 为计算所得冲击响应信号各采样点加速度值的平均值; R²为确定系数, 表征了计算波形和试验波形的近似程度, R²越小, 二者近似程度越高。

4.2 QGA描述

本文采用QGA求解(9)式所示的优化模型。QGA继承了标准遗传算法不依赖梯度, 全局搜索能力强的优点, 并且通过将量子的态矢量表达引入遗传编码, 利用量子逻辑门实现染色体的演化, 从而一条量子状态的染色体可以描述多种系统状态, 有效解决了传统遗传算法收敛速度慢、易陷入局部收敛的问题^[22]。QGA的结构示意图如图 5所示^[22]，图 5中:j为当前进化代数; N为群体数; f(P(t_j))(j=1, 2, …, N)为适应度函数, 与R²一致, 函数中的P(t_j)为运行到第j代的种群。

5 试验验证 5.1 相关试验参数确定

随机选取一个半正弦波形发生器, 测得其承载面积和自由表面面积并计算得n=0.71;测得其硬度H=50。通过冲击试验机铭牌读取其台面质量、底座质量分别为m₁=120 kg, m₂=2 000 kg; 查阅缓冲气囊的恢复力-位移曲线, 并拟合得缓冲气囊刚度方程为k₂(x₂)=-4.8×10^-3x₂+27.6×10^-6x₂²。

5.2 结果

将5.1节中的相关数据代入优化模型, 并结合试验进行参数辨识计算。以台面跌落高度300 mm为例, 基于经典单自由度模型的计算波形和试验波形对比如图 6所示, 显然计算误差较大。

设置12组不同的台面跌落高度进行试验, 并基于提出的非线性动力学模型进行计算, 对比计算结果和试验结果如图 7a)~l)所示。注意:每组跌落高度对应的计算结果为运行QGA 5次并取其中确定系数最小的一次的结果; 各组跌落高度对应的计算波形和试验波形的相关数据如表 1所示; 如图 7a)所示, 表 1中波形脉宽是按JEDEC标准在峰值加速度10%处测得的; 波形确定系数、脉宽误差、峰值误差随跌落高度变化曲线如图 8所示; 为考察算法收敛性和稳定性, 以跌落高度300 mm为例, 运行QGA 10次, 记录各次辨识出的参数值及其标准差如表 2所示, 并绘出前5次的进化过程如图 9所示; 不同跌落高度下辨识出的参数值如表 3所示; 辨识出的所有参数值随跌落高度变化曲线如图 10所示。

6 讨论

由图 6可知，经典单自由度模型的计算误差较大。同时亦可发现，相对于标准半正弦波，半正弦波形发生器冲击响应脉冲具有“尖峰瘦腰”的特性^[10]，并且波形上升段与下降段不对称。文献[8]指出：二者分别是由半正弦波形发生器的超弹性和黏弹性引起的。

由图 7可知：采用本文所建立的非线性动力学模型计算得到的冲击脉冲与试验冲击脉冲的近似程度较高。整体计算精度虽略低于文献[10]的结果，但由图 9可见该方法有着较高的收敛精度和速度，进化速度明显优于文献[10]的结果。

由表 1和图 8可知：随着跌落高度增加(即冲击速度增加)，诸多误差指标也随之增大。造成误差指标不断增大的原因主要是随着冲击速度增加，橡胶材料的硬非线性特性会持续增加，超出了达芬方程的描述能力。但峰值计算误差和脉宽计算误差均在5%以内，而波形确定系数，在跌落高度为1 m以下时均小于15%。如需进一步增大冲击速度，可以考虑其他非线性描述能力更强的模型，如正切型材料模型、双曲正切型材料模型等^[7]。

由表 2可知：本文设计的基于QGA的参数辨识方法辨识出的诸参数值的标准差较低，10次计算中只有1次确定系数稍有波动，即算法稳定性较好。

由表 3和图 10可知：①随着跌落高度增加，α₁(即半正弦波形发生器的线性刚度)有增大趋势。这是半正弦波形发生器硬特性的表现，该趋势也与图 7中试验波形的脉宽随跌落高度增加而减小相对应。文献[23]在研究橡胶隔振器时，通过试验也得到了相同的趋势，但文献[3, 7]通过经典单自由度模型得出半正弦波形发生器刚度与跌落高度无关的结论显然与试验结果不符；②随着跌落高度增加，α₂(即半正弦波形发生器的非线性刚度)也有增大趋势，在本文设计的跌落高度范围内为(4)式所示的非线性刚度的1.3~3.2倍。这是半正弦波形发生器硬非线性的表现，该趋势与图 7中试验波形的“尖峰瘦腰”特性随跌落高度增加而越明显相对应，并且其变化范围也与文献[10]中“半正弦波形发生器的冲击刚度大约是其静刚度的1.95~4.4倍”的结论类似；③随着跌落高度增加，ξ₁(半正弦波形发生器的阻尼比)有减小趋势，在本文设计的跌落高度范围内，ξ₁的变化范围为0.009~0.053。冲击阻尼比较小的原因是冲击载荷作用时间短，在冲击作用之前，阻尼力还来不及从结构中消耗较多能量，显然冲击速度越小，阻尼比越大，冲击速度越大，阻尼比越小，本文辨识出的半正弦波形发生器的冲击阻尼比的变化趋势与之契合。文献[3, 7-8, 21]正是因为冲击阻尼比较小而直接忽略了半正弦波形发生器的阻尼，文献[23]在测定橡胶隔振器冲击刚度时，指出橡胶隔振器的冲击阻尼比大约在0.04左右，该结论与本文的辨识结果类似。④随着跌落高度增加，ξ₂(阻尼器的阻尼比)有增大趋势，在本文设计的跌落高度范围内，ξ₂的变化范围为0.741~1.306。阻尼器的作用本身就是提供阻尼，因此设计值较大是合理的，在本试验中甚至出现过阻尼情况，该趋势与文献[24]中的试验结论契合；⑤随着跌落高度增加，μ(初始速度损失系数)有减小趋势，该结论与文献[25]的结论一致。

7 结论

1) 以冲击试验用半正弦波形发生器的宏观动力学模型为研究对象，基于有阻尼达芬方程建立了半正弦波形发生器的非线性动力学模型，具有形式简单、物理意义明确、预测精度高的优点；同时结合QGA对该动力学模型中的诸参数进行了辨识研究。

2) 借鉴圆柱形橡胶隔振器的静刚度和冲击刚度公式，将波形发生器的尺寸和形状因素考虑进模型之中，同时也提高了算法的收敛速度和精度。

3) 试验结果表明：基于本文提出的半正弦波形发生器非线性动力学模型计算出的加速度波形和试验所得的加速度波形的拟合误差在跌落高度1 m以下时小于15%，而峰值加速度相对误差和脉宽相对误差均在5%以内；本文提出的基于QGA的参数辨识方法具有较好的稳定性，较高的收敛速度和较高的精度。针对辨识出的诸参数值的变化范围和趋势的合理性进行了论证，结果显示二者均符合一定的物理规律，并和前人的研究成果相对应。

4) 本文所提出的半正弦波形发生器非线性动力学模型和参数辨识方法直接基于时域分析获得，避免了诸多转换带来的误差；较好地解决了之前半正弦波形发生器动力学模型忽略了黏弹性特性、没考虑形状和尺寸因素、模型参数获取困难、不关注计算波形与试验波形形状差别的问题；有助于精确预测冲击波形和高效地设计半正弦波形发生器；对于设计、调试、标定冲击试验机也有帮助。