非线性包装系统跌落冲击动力学响应分析的NHB方法

杜兴丹，陈安军，2

(1.江南大学 机械工程学院,江苏 无锡 214122；2.国家轻工业包装制品质量监督检测中心,江苏 无锡 214122)

1968年牛顿首次提出破损边界理论[1]，为缓冲包装理论的发展奠定了重要思想基础,但该理论基于单自由度线性系统。包装工程领域中大多数缓冲包装材料为非线性，其中典型非线性缓冲材料可描述为三次型、正切型、双曲正切型等。缓冲材料的强非线性给理论分析带来较大困难，故非线性缓冲包装系统的动力学分析与评价通常采用数值分析方法[2-3]，但数值分析方法难以获得动力学响应的解析表达，缺少明确的物理意义，给相关参数对系统响应影响分析带来不便。近年来，非线性问题近似分析方法得到发展，如同伦摄动法（HPM）[4-5]、変分迭代法(VIM)[6-7]、何氏频率-振幅关系（HFAF）[8-9]、能量平衡法（EBM）[10-11]等，并广泛应用于非线性工程问题。

针对非线性保守系统，作为牛顿线性化方法和谐波平衡法的有机结合，一种新的非线性分析方法——牛顿谐波平衡法[12-15]（NHB）被提出，该法不仅可以避免使用谐波平衡法时求解复杂非线性代数方程组的问题，且不依赖于小参数，对强非线性振动问题分析可获得满意的精度。虽包装工程产品的跌落冲击动力学问题与振动问题的控制方程相同，但一般振动问题的分析，更多的关注系统的振动频率、周期以及周期解，而包装系统动力学分析与包装系统的缓冲设计需重点关注跌落冲击过程系统响应的位移峰值、加速度峰值以及跌落冲击持续时间等重要参数，且两者初始条件不同，为使NHB方法应用于非线性包装系统跌落冲击问题，需讨论满足系统初始条件的牛顿谐波平衡法的分析逼近解。

1 NHB的基本思想

针对无阻尼非线性缓冲包装系统，其一般形式跌落冲击动力学方程及初始条件可表示为

其中：f(x)是关于位移x的非线性奇函数，引入新变量τ=ωt，方程式（1）改写为

其中：（’）表示对变量τ求导，ω表示系统频率。应用牛顿线性化方法[13-15]，位移和频率分别表示为

其中：Δx1和Δω12分别是x1和ω12的微小增量。为了简化表述，令ω2=Ω，方程式（5）改写为

式（4）和式（6）代入方程式（3），有

将方程式（7）线性化，且不计高阶小量，有

式中：fx(x1)表示f(x1)对x的1阶偏导。为了满足初始条件表达式（2），令初始逼近解为

式（9）代入方程式（8），有

其中：Δx1是一个关于τ的以2π为周期的周期函数，满足初始条件

为了获得1阶近似解，令

将式（12）代入方程式（10），经整理并令sin(τ)前系数为0，则1阶近似解

为了获得2阶近似解，令

将式（15）代入方程式（10），令sin(τ)和sin(3τ)的系数分别为0，可得关于c和ΔΩ1的线性方程组，解该方程组获得参数c和ΔΩ1表达式，则二阶近似解

2 NHB在跌落冲击中的应用

为验证NHB方法的有效性，以三次非线性缓冲包装系统为研究对象，获得跌落冲击动力学响应的解析表达式，并与変分迭代法[16-17]和龙格-库塔法比较，验证NHB近似解析解的精度。

2.1 模型与方程

对无阻尼三次非线性包装系统，跌落冲击动力学模型如图1示，跌落冲击动力学方程及初始条件可表示为

式中：m0为产品质量，x0为压缩变形量，k0是弹性系数，r为非线性系数，g为重力加速度，h为跌落高度，?为产品加速度。引入参数无量纲位移及时间参数分别为无量纲动力学方程及初始条件改写为

引入新变量τ=ωt1，方程式（20）改写为

将方程式（4）和式（6）代入到方程式（22）中，并经线性化处理，不计高阶小量，有

为满足跌落冲击初始条件，方程式（9）为初始逼近解，代入方程式（23）得

图1 三次非线性包装系统的跌落冲击模型

注意到

将式（25）代入方程式（24）整理得

2.2 1阶近似解

为获得一阶近似解，将式（12）代入方程式（26）有

令sin(τ)的系数为0，得1阶近似解

2.3 2阶近似解

为了获得2阶近似解，将式（15）代入方程式（26），有

令sin(τ)及sin(3τ)的系数为0，得

将式（28）代入式（31）和式（32），解上述方程组得

则2阶近似解为

跌落冲击动力学方程式（22）的位移、速度、加速度2阶解析表达式分别为

跌落冲击持续时间为

为了求解未知参数A，联立无量纲初始条件式（21）、速度表达式（37）和频率表达式（34）有

由式（39）、式（40）、及式（41）知，系统跌落冲击时间、无量纲位移最大值与加速度最大值取决于参数A。对给定系统，当跌落高度或无量纲跌落冲击速度已知，参数A可由式（42）式确定，且参数A取决于系统跌落高度或无量纲跌落冲击速度，进而通过式（39）、式（40）及式（41）等解析表达式确定无量纲位移、加速度最大值及跌落冲击时间，分析系统跌落高度与相关参数对跌落冲击响应的影响。

3 算例分析

为了评价运用NHB法求解跌落冲击问题近似解的精度，参考文献[16-17]算例，相关参数：m=10kg，h=0.6m，k0=600Ncm-1和r=72Ncm-3，应用R-K、NHB和VIM[16-17]法获得系统跌落冲击的无量纲位移、加速度响应分别如图2、图3所示。

图2 跌落冲击无量纲位移响应

图3 跌落冲击无量纲加速度响应

为进一步探讨跌落高度对NHB方法的精度影响，分别取跌落高度0.2 m、0.4 m、0.6 m、0.8 m、1.0 m，通过R-K、NHB和VIM法获得的跌落冲击无量纲最大位移、最大加速度以及跌落冲击持续时间对比及相对误差分析如表1、表2和表3所示。

由图2、图3和表1、表2、表3知：

(1)相同条件下，NHB获得的2阶近似解析解与龙格库塔法结果最为接近。通过位移、加速度响应分析可知，NHB2阶近似解精度优于VIM近似解，VIM近似解精度优于NHB1阶近似解。跌落冲击时间相对误差分析表明，NHB2阶近似解精度优于1阶近似解，NHB1阶近似解精度优于VIM1阶近似解。

(2)随跌落高度的增加，NHB1阶、NHB2阶解析解获得的产品最大位移、最大加速度以及跌落冲击时间相对误差均出现增大，但在包装工程应用中其跌落高度一般不超过1 m，在此条件下NHB2阶近似解最大相对误差控制在2.5%以内，可满足工程设计的要求。

表1 不同跌落高度无量纲位移响应最大值比较

表2 不同跌落高度无量纲加速度响应最大值比较

表3 不同跌落高度跌落冲击无量纲持续时间比较

4 结语

针对非线性包装系统跌落冲击动力学问题，应用NHB方法，获得系统跌落冲击位移响应、加速度响应及跌落冲击持续时间的1阶和2阶近似解析解，且跌落冲击位移、加速度最大值以及跌落冲击持续时间为代数方程，形式简单，物理意义明确，易于分析系统相关参数对跌落冲击响应的影响。与龙格-库塔法比较，结果表明由NHB获得的2阶近似解具有满意的精度，且其精度明显优于変分迭代法。随高度增加，尽管解析解精度有所降低，但当跌落高度在一定范围之内，解析解精度仍能够满足工程设计需求。