朱大鹏，魏 洁

(兰州交通大学 交通运输学院，兰州 730070)

包装件在流通过程中，在装卸、运输、储存等环节，会受到各种外界载荷的作用，如跌落冲击、随机振动、静压力载荷等。其中跌落冲击载荷的作用较为剧烈，很容易引起包装件的损坏，因此，国内外研究者在跌落冲击载荷条件下对包装件的动态响应、破损边界曲线、缓冲材料的保护能力等开展了一系列研究[1-7]。在运输过程中，包装件长时间受到随机振动载荷的激励。在随机振动条件下，包装件中脆弱部件的加速度响应可能会超出产品的脆值，从而造成包装件损坏。因此，很有必要研究包装件在运输过程中的随机振动载荷作用下的振动可靠性，可靠性分析为优化包装设计和产品设计，提高包装件和产品的可靠性，减少包装件在流通过程中的损坏等提供了理论基础和优化依据。目前，研究包装件振动可靠性的文献较少，文献[8-10]应用一阶可靠性法(First Order Reliability Method,FORM)和镜像激励法分析了包装件的振动可靠性，但未考虑包装件模型中参数的不确定性。

由于不确定性不可避免地存在于包装件模型的各参数中，如缓冲材料动态特性的不确定性、脆弱部件和产品主体之间连接特性的不确定性，这使得包装件成为一个参数不确定的系统。在随机振动激励下，由于包装件参数不确定性的影响，包装件振动可靠性也是不确定的，因此，很有必要研究包装件参数不确定对振动可靠性变化的影响，并对各参数进行灵敏度分析，以具体量化各不确定参数的变化对可靠性变化的影响程度，为包装件和产品的优化提供依据。为了对由参数不确定性所引起的系统响应的不确定度进行合理的量化，人们提出了各种分析和模拟方法，其中原始蒙特卡洛法是一种通用、准确的方法[11]，但该方法需要进行大量的数值模拟，收敛速度慢，对于振动可靠性问题，由于失效概率通常较小，为确保准确性，需生成大量的数据进行分析，导致分析速度很慢，这也限制了该方法的应用，该方法目前通常用于验证其他方法的准确性。目前对系统响应进行不确定度量化的理论分析方法主要包括Karhunen-Loeve展开法[12]、正交级数展开法[13]、优化线性估计法[14]、多项式混沌展开法[15]等，其中，多项式混沌展开法在不确定度量化中得到了广泛的应用。该方法将随机变量或随机过程表示为标准正态随机变量的正交多项式与确定性参数乘积相叠加的形式，多项式混沌展开可准确表达具有一定统计特性或分布特性的随机变量或随机过程，目前该方法广泛应用于不确定系统的响应分析和量化[16-17]、不确定参数识别[18-19]、灵敏度分析等[20]。灵敏度分析用于量化不确定参数的变化对系统响应变化的影响程度，是进行系统优化的一个重要步骤。目前灵敏度分析主要包括三大类：局部灵敏度分析、全局灵敏度分析、区域灵敏度分析。局部灵敏度分析用于评估不确定参数在名义值处的变化对系统响应的影响程度[21]，而全局灵敏度分析研究在整个不确定参数范围内参数变化对响应变化的影响[20-22]，区域灵敏度分析用于评估不确定参数处于不同分布区间时，其变化对系统响应变化的影响[23]。

本文考虑包装件中参数的不确定性，采用多项式混沌展开量化由于参数不确定性引起的随机变化的可靠性指标，由于采用分析法研究不确定参数对可靠性变化的影响较为复杂，本文采用非嵌入式方法[24]对包装件中不确定参数均匀合理采样，应用内积法求解混沌多项式的系数，获得随机变化的可靠性指标的多项式混沌展开表达式。为确保准确性，采用原始蒙特卡洛法对该分析方法进行验证。基于不确定参数的多项式混沌展开表达，分析包装件的全局灵敏度指数，量化分析包装件中各参数变化对包装件可靠性变化的影响程度。

1 参数确定的包装件振动可靠性分析

1.1 随机振动激励的表示

(1)

(2)

Σ=Φ·Λ·ΦT

(3)

式中:矩阵Φ=[Φ1,Φ2,…,ΦN]是关于特征向量Φi(i=1,2,…,N)的正交矩阵，上标T表示矩阵的转置，矩阵Λ=diag[λ1,λ2,…,λN]是关于特征向量的对角矩阵，则零均值随机振动激励可利用Karhunen-Loeve展开式表示为[26-27]：

(4)

(5)

的最小的r值即为主成分的个数。应用式(2)～(5)，可将具有一定谱特征的随机振动表达到标准正态变量空间u中。

1.2 包装件的振动可靠性分析

a(t)·u

(6)

式中，a(t)=[a1(t),a2(t),…,an(t)]T，

(7)

根据式(6)，定义包装件中关键部件在随机振动激励下的极限状态方程：

g(u,Gc,t)=Gc-a(t)·u

(8)

式中:Gc为包装件的脆值。对于线性包装件，式(8)所示的极限状态方程为r维标准正态空间中的一个超平面，包装件在随机振动激励作用下，若g(u,Gc,t)>0，则包装件安全，反之，若g(u,Gc,t)≤0，则包装件发生首次穿越损坏。应用一阶可靠性方法，线性包装件在式(4)所表达的随机振动激励下，损坏概率为：

(9)

式中:Ф[·]为标准正态累积分布函数，参数β称为可靠性指标，其定义为坐标原点到g(u,Gc,t)所示的超平面的距离：

(10)

2 参数不确定性的包装件振动可靠性分析

2.1 不确定参数的表示

多项式混沌展开是分析不确定系统的有效工具，对于任意有限方差的不确定参数θ，应用多项式混沌展开法可表示为[28-29]：

(11)

式中:m为随机变量的个数，φi是向量ξ=[ξ1,ξ2,…,ξm]T的函数，其中，ξi(i=1,2…,m)为标准正态随机变量，φi类型由θ的分布类型确定，如高斯分布的变量θ由Hermite多项式表示，Gamma分布的变量θ由Laguerre多项式表示，Beta分布的变量θ由Jacobi多项式表示等[30]，参数a为确定性参数。式(11)可简写为：

(12)

2.2 包装件可靠性指标的表示

由于缓冲材料特性存在着不确定性，包装件中产品的生产、装配等环节也不可避免地存在着不确定因素，导致包装件模型中的参数通常是不确定的。这些参数不确定性会影响到包装件的振动可靠性，因此包装件在随机振动作用下的可靠性也是不确定的，表征包装件振动可靠性的参数是式(10)所示的可靠性指标β，即β是一个随机变量，应用多项式混沌展开可将β表示为：

(13)

式中:函数Ψj(η)(j=1,2,…,∞)是关于向量η的正交多项式，参数bj(j=1,2,…,∞)是确定性参数，η=[η1,η2,…,ηm]为包含各不确定参数的向量，m为不确定参数的个数。由于不确定参数η不一定是标准正态变量，为构建可靠性指标β的多项式混沌展开，需要应用等概率变换原则[31]在随机变量ηi和标准正态随机变量ξi之间构建一一映射关系：

i=1,2,…,m

(14)

式中:Fi是不确定参数ηi的累积分布函数，Φ是标准正态变量ξi的累积分布函数，上标-1表示逆累积分布。因此，不确定参数向量η可用标准正态随机向量ξ表示：

η=T(ξ)

(15)

式中:T(·)表示等概率转换[31]。为减少分析计算量，通常采用有限的M项正交多项式近似表达随机变量，式(13)可写为：

(16)

式中:Ψ(ξ)=[Ψ0(ξ),Ψ1(ξ),…,ΨM-1(ξ)]，b=[b1,b2,…,bM-1]T。根据多项式混沌展开的收敛性分析[32]，通常用三阶或四阶多项式混沌展开即可准确表达一个有限方差不确定变量，如果阶次p超过四阶，则对分析结果影响不大，则式(16)中M可由下式确定[32]：

(17)

2.3 求解多项式混沌展开的系数

在求解不确定参数的多项式混沌表达式中的系数b时，目前广泛采用的方法有两类，嵌入式和非嵌入式。其中嵌入式方法基于Galerkin法，通过求解确定性方程可得多项式混沌展开的系数，对于包含多个随机变量的系统，应用该方法确定随机变量的多项式混沌展开表达式的系数较为复杂。本文采用非嵌入方法合理选择随机参数的采样点，根据这些采样点，应用最优化方法或者内积法求解系数b。其中，最优化方法应用下式求解：

b=argminE[(β(η)-Ψ(ξ)b)2]

(18)

式中:E[·]表示期望值。应用式(18)求解b时，为提高准确性，通常要求采样点数n至少为M的三倍，在采样点中，考虑到随机参数的分布特征，通常在参数均值附近的采样点较多，在概率分布的尾部采样点较少，应用总体最小二乘法分析后，可靠性指标的分析值在其概率密度分布的尾部会出现较大的误差。概率分布的尾部值代表了随机变量数值较大和较小的部分，对于分析包装件振动可靠性非常重要，为提高分析精度，本文采用内积法求解系数b。式(16)两边同乘向量Ψj(ξ)，考虑到混沌多项式的正交性，可得下式：

(19)

式中:〈·,·〉表示总体平均，其运算为Hilbert空间的向量内积：

(20)

式中:W(ξ)为积分过程中的权值函数，其定义为：

(21)

根据式(21)和式(19)可求得式(16)中可靠性指数的多项式混沌展开式的系数bj。其中式(19)中，分母的内积可采用理论分析的方法获得，分子的内积采用数值分析的方法获得。求得多项式混沌展开的系数后，包装件振动可靠性指标可用下式进行估计：

(22)

考虑到Ψ的正交特性，可靠性指标的均值和方差为：

(23)

(24)

3 可靠性指标的全局灵敏度分析

本文采用Sobol法分析可靠性指标的全局灵敏度，对可靠性指标β进行Sobol分解[34]，将其表达为：

(25)

则β的方差可表达为：

(26)

根据文献[20]中的推导，方差D可表示为：

(27)

其中

1≤i1

(28)

全局灵敏度指数定义为：

(29)

根据式(27)和(29)，全局灵敏度指标满足关系式：

(30)

将式(22)所示可靠性指数进行Sobol分解：

(31)

为确保式(31)和式(25)等效，必须准确定义集合I，如在式(31)中，用Ii可选出所有的单变量项，用Ii,j可选出所有的双变量项，用Ii,j,k可选出所有三变量项。I定义为：

Ii1,i2,…,is=

(32)

(33)

在多个不确定因素i1,i2,…,is共同作用下，可靠性指标β的高阶灵敏度指标为：

(34)

由以上分析可知，在获得了随机变量的多项式混沌展开表达式后，通过合理选择多项式混沌展开式中的对应项，可方便地求出各随机参数单独作用和多个随机参数共同作用时包装件可靠性指标的Sobol指标，可量化分析随机参数变化对可靠性指标变化的影响。

4 实例分析

对于非线性包装件，为简化分析，通常采用一个等效的线性包装件来替代非线性包装件，因此，本文采用二自由度线性系统模拟包装件，如图1所示。

图1 包装件简化模型

本文采用EPE泡沫作缓冲材料，厚度为30 mm，密度为21 kg/m3，用质量块-泡沫系统模拟包装件，质量块重1.1 kg，用冲击锤激励质量块，记录系统响应，实验系统如图2所示。用自由响应数据分析系统的固有频率ω1和阻尼比ζ1[35]，由于EPE泡沫特性的不确定性，识别出的系统固有频率和阻尼比也呈随机变化，本文用100块缓冲材料试样识别出的系统的ω1和ζ1的随机分布情况如图3所示。用极大似然估计法分析ω1和ζ1的概率密度分布情况，分析不同概率分布和真实的数据分布之间关于累积分布函数的误差，分析结果如图4所示，图4(a)和(b)分别给出了ω1和ζ1的实验数据实际概率和不同的概率分布方程曲线的对比图，并给出了不同概率分布条件下的累积分布误差。由图4可以看出，ω1符合logistic分布，ζ1符合loglogistic分布。由于产品主体和关键部件之间的连接特性参数ω2和ζ2很难直接用实验的方法获得，为考虑ω2和ζ2的随机变化对包装件振动可靠性的影响，本文假定ω2和ζ2呈正态分布，假定ω2均值为500 rad/s，ζ2的均值为0.1。包装件模型中的不确定参数的基本情况，如表1所示。

图2 测试ω1和ζ1的实验装置

(a)

(a)

表1 包装件模型中不确定参数的基本情况

由于包装件模型中参数值的随机变化，包装件在随机振动激励下的振动可靠性指标β也对应地是一个随机变量，由于随机参数的数量为4，在多项式混沌展开式中，假定阶次为4，根据式(17)，用多项式混沌展开时，展开式中包含70项，根据拉丁超立方合理采样原则，应用非嵌入法分析可靠性指标时，采样数目至少为210个以上，本文采样点数选用500个。应用非嵌入法分析包装件可靠性指标的具体步骤如下：

(3)包装件运输过程中的随机振动，如图5所示为卡车在一段公路行驶时记录设备记录的随机振动的PSD曲线，车辆行驶速度为60 km/h，车辆载重为20 t。根据PSD曲线，根据式(4)，将该随机振动表示在标准正态随机变量空间；

图5 记录的随机振动PSD曲线

(4)根据ηi具体值，可求出包装件的加速度单位脉冲响应函数为ha。应用式(6)、(7)、(10)计算包装件的可靠性指标，计算出的可靠性指标的概率分布情况如图6所示，分析可靠性指标的概率分布类型可知，可靠性指标基本呈正态分布，累积概率分布误差为1.17%；

图6 用非嵌入法分析的可靠性指标的分布

(5)用多项式混沌展开式表达随机可靠性指标：根据可靠性指标的分布类型，由于可靠性指标呈正态分布，故选用正交Hermite多项式作为式(16)中标准正态随机变量的多项式类型，根据式(19)～(21)计算多项式的系数。

用多项式混沌展开法分析的可靠性指标的概率密度方程曲线如图7所示，本文采用原始蒙特卡洛法(N=105)验证了分析结果的准确性，从图7中可以看出，本文采用的分析方法具有良好的准确性，与原始蒙特卡洛法相比，该方法大大提高了计算效率。

图7 包装件振动可靠性指标的概率密度分布

本文用多项式混沌展开准确表达包装件振动可靠度指标的概率分布情况，因此，可利用式(33)和(34)计算可靠性指标β对各参数的灵敏度指标，分析结果如图8所示。从图中可以看出，缓冲材料的阻尼的单参数灵敏度最大，产品的脆弱部件和产品之间的阻尼比的单参数灵敏度最小，各参数的总体灵敏度比单参数灵敏度均有明显的提高，这表明各个随机参数共同作用时所引起的可靠性指标的变化，相比随机参数单独作用所引起的可靠性指标的变化，有明显的提高，这也表明，在参数不确定性的影响下，包装件振动可靠度指标的变化主要是由各不确定参数共同作用所引起的。

图8 可靠性指标对各随机参数的灵敏度

5 总结与展望

包装件在流通过程中，由于随机振动的激励作用，包装件中产品可能发生加速度首次穿越损坏，本文用可靠性指标表征包装件的瞬时失效概率。由于包装件的参数是不确定的，如缓冲材料的动态特性、产品主体和产品脆弱部件之间连接特性，这些参数不确定性会引起包装件振动可靠性指标的不确定性，造成了包装件的可靠性不是一个确定的值，而是具有一定概率分布特性的随机变量。本文研究在模型参数变化的情况下，包装件可靠性指标的变化情况，研究可靠性指标相对于各模型参数的灵敏度。

本文应用Karhunen-Loeve展开法将具有一定谱特征的平稳随机振动表达在标准正态随机变量空间中，应用FORM法计算包装件的振动可靠性指标。考虑包装件中缓冲材料的动态特性参数的不确定性参数ω1和ζ1，以及产品主体和脆弱部件连接特性的不确定性参数ω2和ζ2，应用等概率变换的原则，将这些不确定参数等效转换到四维标准正态随机变量空间中，采用非嵌入法分析在这些不确定的包装件参数变化时振动可靠性指标的变化情况。为在确保分析精度的条件下减少分析工作量，在四维标准正态空间中，采用拉丁超立方采样原则，将采样的标准正态随机变量转换为包装件参数，分析得到对应的振动可靠性指标。根据振动可靠性指标样本的分布情况合理选择多项式混沌展开式中的正交多项式类型，应用非嵌入法得到该展开式中的系数。在准确获得可靠性指标的多项式混沌展开式后，应用Sobol法分析可靠性指标的全局灵敏度。根据本文的分析，在分析包装件振动可靠性时，考虑到不确定因素的影响，可得以下结论：

(1)考虑到缓冲材料特性的不确定性，包装件中产品主体和脆弱部件之间连接特性的不确定性，包装件在随机振动作用下的振动可靠性是一个不确定变量，基本呈正态分布，在包装件设计、优化和振动可靠性分析时，需考虑振动可靠性的不确定性。

(2)包装件可靠性指标的波动主要是由多个包装件不确定参数共同作用引起的，如果要减小振动可靠性的不确定性，需同时减小各参数的波动。

(3)对于单个包装件参数，根据灵敏度分析结果，缓冲材料阻尼的不确定性对包装件可靠性指标的波动影响最大，其次是缓冲材料的弹性特性的不确定性，而产品主体和脆弱部件之间的阻尼的不确定性对振动可靠性的波动影响最小，其影响几乎可以忽略。

在实际应用中，为优化包装件参数，提高包装件振动可靠性，今后可在以下两个方面开展深入研究：① 实际的包装件的结构、特性较为复杂，本文用二自由度集中参数模型替代真实的包装件，可以总体性地了解缓冲材料特性的不确定性、产品主体和脆弱部件之间连接特性的不确定性对包装件振动可靠性波动的影响，为更加清晰了解包装件真实振动可靠性及参数敏感度，需构建能够真实反映包装件动态特性的有限元模型，并在模型中分析振动可靠性和参数敏感度；② 考虑到随机振动条件下的包装件的振动可靠性，需构建包装件在随机振动条件下的设计和优化方法，提高包装件在随机振动条件下的可靠性，或在振动可靠性约束条件下，实现包装成本最优。