包裝件在非高斯隨機載荷下響應特征分析方法研究

朱大鵬， 王浩然， 曹興瀟

(1. 蘭州交通大學 交通運輸學院，蘭州 730070； 2. 蘭州交通大學 機電工程學院，蘭州 730070)

包裝件在流通過程中，在車輛運輸環節，長時間受到隨機振動載荷作用。包裝件在隨機振動載荷激勵下，內裝的產品或產品的關鍵部件可能會發生首次穿越破壞或者疲勞損傷。為提高包裝件在隨機振動載荷條件下的可靠性，近年來，國內外研究者圍繞隨機振動環境模擬、包裝件響應分析等，開展了大量研究工作。Lepine等[1]介紹了目前公路運輸中車輛隨機振動載荷的模擬和包裝件振動測試方法，包括非高斯和非平穩隨機振動模擬方法。曾昕等[2]基于運行測試方法，提出一種非平穩指數的表征方法，基于峭度和功率譜密度函數等約束條件，利用貝塔分布隨機數模擬生成非平穩非高斯隨機振動。Lamb等[3]根據采集的公路運輸條件下車輛內的振動數據，分析了道路類型、車輛懸掛參數等對車輛隨機振動的影響。Hosoyama等[4]考慮了隨機振動的非高斯特性，將包裝件建模為單自由度振動系統，研究包裝件動態特性對響應峭度的影響，為優化包裝設計提供了理論依據。Wang等[5]考慮到包裝件的復雜性，基于包裝件中元件線性振動假設，提出了運輸包裝件加速振動測試方法。王志偉等[6]研究了不同包裝件約束、不同振動等級條件下，不同譜型的隨機振動激勵對包裝件加速度響應和動壓力的影響。文獻[7]分析了隨機振動載荷下，無約束包裝件的跳動對包裝件加速度響應非高斯特性的影響。于水源等[8]采用數值模擬方法對包裝件疲勞特性進行了分析。

將包裝件建模為振動系統，分析包裝件在隨機振動條件下的響應特征是分析包裝件振動可靠性或疲勞特性的重要基礎，也是對包裝件進行優化的重要前提。如果隨機振動激勵是平穩高斯的，對于線性包裝件，可采用分析法得到包裝件的響應。對于非線性包裝件，目前可采用FPK(fokker-planck-kolmogorov)法、攝動法、隨機平均法、等效線性化法等分析其隨機響應。由于道路不平[9]、包裝件在運輸過程發生跳動、車輛懸掛特性具有時變性[10]等原因，包裝件在運輸過程中經常受到非高斯隨機振動激勵。近年來，非高斯隨機振動的模擬及非高斯隨機振動條件下包裝件的響應分析已成為運輸包裝領域的研究熱點。蔣瑜等[11-12]基于非高斯隨機振動的PSD，提出了一種新的基于幅值調制和相位重構的算法，生成不同偏斜度和峭度的非高斯隨機振動。Grigoriu[13-14]采用非線性零記憶變換的方法，對高斯分布的隨機振動進行單調變換，實現非高斯隨機振動的模擬。楊喆等[15]結合多項式混沌展開方法和Karhunen-Loeve展開，用標準正態隨機過程的非線性函數表達非高斯隨機過程。

在非高斯隨機振動載荷激勵下，包裝件的響應分析較為困難，即使對于線性振動系統，分析其響應的計算量也較大[16]。對于非線性包裝件，很難采用分析法直接獲得其響應的統計特征參數，目前主要通過蒙特卡洛模擬法分析包裝件的響應統計參數，根據這些統計參數再現包裝件的隨機響應。由于傳統蒙特卡洛法計算量很大，可采用擬蒙特卡洛分析法[17]，通過合理的隨機變量采樣策略，提高傳統蒙特卡洛法的分析效率，提高響應統計參數的收斂速度。雖然采用擬蒙特卡洛法可改進蒙特卡洛分析的效率，但該方法與分析法相比還存在著計算量大、分析效率低的缺點。

本文提出非高斯隨機振動載荷條件下，非線性包裝件的加速度響應統計特征分析方法。首先，本文采用非高斯Karhunen-Loeve展開法[18-20]模擬非高斯隨機振動載荷，將非高斯隨機振動載荷表達為非高斯隨機變量的線性組合。對非線性包裝件的響應進行一階泰勒展開，根據該展開式，分析響應的各階矩統計參數和自相關函數，采用鞍點估計法分析包裝件加速度響應的概率密度函數(probability density function,PDF)和累積分布函數(cumulative distribution function,CDF)，從而可分析包裝件的振動可靠性、疲勞特性等。由于本文方法中，非高斯隨機振動表示為非高斯隨機變量的線性組合，以該非高斯隨機振動的模擬方法為基礎，對非線性包裝件的響應進行的線性化處理(一階泰勒展開)具有誤差小的特點，這與文獻[21]提出的非線性包裝件等效線性化處理方法在本質上是一致的。因此，本文的方法與傳統的基于蒙特卡洛模擬方法相比，具有分析效率高、精度好的優點。

1 非高斯隨機振動的模擬算法

(1)

(2)

向量ξi滿足

E[ξi]=0

(3a)

E[ξiξj]=δij

(3b)

式中，δij為Kronecker-Delta函數。

實際應用中，通常在式(1)中采用有限項表示隨機過程，式(1)可簡化為[22]

(4)

對于呈高斯分布的隨機振動，采用Karhunen-Loeve展開表示該隨機振動時，在式(4)中，ξi為標準正態隨機變量，且滿足式(3a)和式(3b)。對于非高斯隨機振動，則需要對式(4)中的ξi分析，通過合理選擇隨機變量ξi的概率密度函數，實現非高斯隨機振動的數值模擬。確定ξi的算法如下：

(5)

(6)

步驟2根據累積分布函數F，生成M組隨機向量ξi(θm),i=1,2,…,M,m=1,…,N，對其進行標準化處理，使其均值和方差分別為0和1，且滿足式(3a)和式(3b)，將這N組隨機向量ξi(θm)代入式(4)中生成N個非高斯隨機振動樣本

(7)

式中：k為迭代次數；m為模擬的非高斯隨機振動樣本編號

(8)

步驟5用式(9)對ξi進行更新迭代

(9)

步驟6重復步驟2～步驟5，直至模擬的非高斯隨機振動累積分布函數與目標累積分布函數誤差小于給定的閾值，或相鄰迭代步驟中誤差變化很小。

在以上算法步驟中，不僅需要保證模擬的非高斯隨機振動累積分布函數和目標累積分布函數一致，還需要保證在迭代過程中，模擬的非高斯隨機振動的自相關函數與目標自相關函數C(τ)一致。為達到該目標，在每次迭代中，計算出的Karhunen-Loeve展開式中的隨機變量ξi(θm)(i=1,2,…,M,m=1,…,N)應滿足式(3a)和式(3b)，以確保在迭代過程中模擬的非高斯隨機振動的自相關函數C(τ)保持不變。本文采用優化拉丁超立方采樣法[24]，在不改變隨機變量ξi(θm)(i=1,2,…,M,m=1,…,N)的累積分布的基礎上，調整隨機變量ξi的分布位置，減小隨機變量之間相關性，使其滿足式(3a)和式(3b)。首先建立一個N×M矩陣X，在該矩陣中放置生成N個非高斯隨機振動隨機變量ξi(θm)，該矩陣的每列隨機變量符合累積分布函數F，對每列隨機變量，根據其大小排序編號，與矩陣X對應，我們構建一個N×M矩陣R，該矩陣中各元素為矩陣X中各元素的排序編號，即矩陣R各列為1～N的正整數。矩陣R第i列和第j列之間的相關性可由M×M矩陣T定義，該矩陣中的元素Tij為Spearman排序相關系數，由式(10)定義[25]

(10)

由于T是一個正定矩陣，對T進行Cholesky分解運算

T=Q′Q

(11)

式中，Q′為矩陣Q的轉置。對初始的包含隨機變量排列位置的矩陣R進行變換

R′=RQ-1

(12)

根據矩陣R′每列中包含的隨機變量排列位置，調整隨機變量ξi的位置，可大大降低隨機變量ξi之間的相關性。

對隨機變量ξi位置調整后，還需對其標準化，代入Karhunen-Loeve展開式中，式(7)改寫為

(13)

以上方法在隨機變量ξi的分布不變的條件下，調整了其位置，并對其標準化處理，其目的在于令隨機變量ξi滿足式(3a)和式(3b)。根據Phoon等的推導，滿足式(3a)和式(3b)的隨機變量ξi代入式(7)中，模擬的非高斯隨機振動信號的自相關函數與目標自相關函數保持一致。因此，在每次迭代中，模擬的隨機振動的自相關函數保持不變，且與目標自相關函數保持一致，確保了模擬的振動信號頻域特性準確性。

2 包裝件振動響應模擬

將包裝件中的產品看作是一個剛性質量塊，將包裝件建模為支座激勵單自由度振動系統，如圖1所示，其運動方程式為

圖1 單自由度包裝件模型Fig.1 Single degree of freedom package model

(14)

(15)

式中，μξi為ξi的均值。為計算式(14)中的偏微分項，對式(14)中的各項對隨機變量ξi求偏微分

(16)

式中，緩沖材料非線性項對隨機變量ξi的偏微分為

(17)

結合式(16)和式(17)，得

(18)

(19)

(20)

根據式(15)，包裝件加速度響應的均值為

(21)

式中， E[·]為期望值運算，由于在式(21)中

故包裝件響應的均值為

(22)

式中：μξ=[μξ1,μξ2,…,μξM]；μξi為ξi的均值，根據式(15)和式(22)，有

(23)

對式(23)的平方求期望值，可得

(24)

由于隨機變量ξi,(i=1,…,M)的方差為1，且各隨機變量之間不相關，故式(24)中，E[(ξi-μξi)2]=1，且當i≠j時，E[(ξi-μξi)(ξj-μξj)]=0，故根據式(24)可得包裝件響應的方差表達式

(25)

同理，可得包裝件加速度響應的三階統計矩和四階統計矩的表達式

(26)

(27)

包裝件加速度響應的自相關函數為

(28)

對于平穩非高斯隨機振動，式(28)可簡化為單參數函數

(29)

式中，τ=t2-t1。

3 振動響應統計特征分析

根據式(22)、式(25)～式(27)可得包裝系統加速度響應的統計特征值，2020年朱大鵬根據這些統計特征值，結合多項式混沌展開和Karhunen-Loeve展開，模擬加速度響應的時域信號，采用擬蒙特卡洛法分析系統的振動可靠性。該方法雖然可以準確分析出系統可靠性，但計算量大，不適用于復雜系統的可靠性分析，不適用于需要多次分析可靠性的系統優化場合。鞍點估計法[26-27]是分析振動系統可靠性的一種高效準確分析法，在已知式(15)中隨機變量ξi的PDF條件下，鞍點估計可準確估計出包裝系統加速度響應的CDF，從而提供了一種分析響應可靠性的高效方法。但該方法需分析ξi的累積母函數，在式(15)中，ξi經過多次迭代優化后其PDF無法用分析式準確表達，或ξi的累積母函數非常復雜，導致鞍點估計計算量很大，這些因素都限制了鞍點估計法的應用。本文采用基于四階統計矩參數的鞍點估計法[28-29]分析包裝系統的可靠性，避免了分析ξi的累積母函數，具有良好的分析效率和通用性。

(30)

(31)

根據式(30)的定義，可得

(32)

(33a)

(33b)

(33c)

(33d)

(34)

(35)

(36)

(37)

(38)

式中：Φ和φ分別為標準正態分布的累積分布函數和概率密度函數；參數w和v由式(39a)和式(39b)確定

(39a)

(39b)

式中，sgn(ts)=1,0,-1，對應的ts分別大于、等于、小于0。

根據式(36)和式(38)，可得包裝件響應的PDF和CDF，可高效確定包裝件響應的統計特征，包裝件的振動可靠性可由響應的CDF確定。

4 實例分析

本部分采用實例介紹非高斯隨機振動模擬、非線性包裝件振動響應分析方法、振動可靠性分析方法。

圖2 試驗記錄的非高斯振動的時域信號和PSD曲線Fig.2 Recorded field non-Gaussian vibration time domain data and PSD curve

(40)

圖3 的PDF和自相關函數Fig.3 PDF and autocorrelation function of

圖4 模擬的非高斯隨機振動的時域信號和迭代后ξi的PDFFig.4 Simulated time domain non-Gaussian random vibration and PDF of ξi after iterations

表1 包裝件加速度響應統計參數Tab.1 Statistical parameters of package acceleration response

圖5 包裝件加速度響應的時域信號和自相關函數Fig.5 Time domain package acceleration response and the autocorrelation function

圖6 的PDF和CDFFig.6 Analytical PDF and CDF of

5 結 論

本文提出一種平穩非高斯隨機振動載荷下包裝件加速度響應統計特征的高效分析方法，與傳統的基于蒙特卡洛和擬蒙特卡洛分析法相比，本文的方法具有計算效率高、精度好的優點，特別適用于需要進行重復的包裝件加速度響應統計特征分析和可靠性分析的包裝系統優化、包裝件參數優化的場合。論文主要成果如下：

(1) 為確保分析精度，本文將非高斯隨機振動激勵表示為非高斯隨機變量的線性組合，在隨機變量均值處對包裝件加速度響應進行一階泰勒展開估計，由于未對隨機變量進行非線性變換，故采用一階泰勒展開估計的誤差較小，可確保本文方法的準確性。

(2) 本文提出了確定包裝件加速度響應統計特征的分析方法，與蒙特卡洛和擬蒙特卡洛方法相比，具有快速高效的優點。

(3) 確定包裝件響應加速度統計特征需要大量的加速度響應時域樣本，因此需要大量的蒙特卡洛或擬蒙特卡洛分析，造成模擬計算量大、計算效率低，本文依據包裝件加速度響應的前四階統計矩，采用鞍點估計法，可高效準確得到包裝件加速度響應的PDF和CDF，避免了蒙特卡洛或擬蒙特卡洛分析。

(4) 本文將復雜包裝件簡化為單自由度振動系統進行分析，本文基于該簡化模型提出了包裝件振動響應統計特征分析方法，包裝件的簡化模型的準確性對于分析包裝件在運輸過程中的響應、評價包裝件可靠性、對包裝件進行優化設計等至關重要。因此，需建立包裝件等效簡化模型分析方法，根據包裝件振動實驗數據識別包裝件簡化模型彈性特性和阻尼特性的類型，并識別具體的模型參數。該工作是今后本領域的一個重要研究方向。