包裝件隨機振動時域響應和PSD 快速獲取方法研究

薛陽，宋衛生

（河南牧業經濟學院，河南 鄭州 450046）

在公路運輸環節，貨物（包裝件）受到的激勵主要是隨機振動，來自運輸工具（車輛）底板。要深入地分析包裝件或包裝件內易損零部件的振動行為，需要得到包裝件或易損件的振動響應和其功率譜密度(PSD)。目前有兩種方法可以實現，分別是從試驗和仿真分析的角度出發。

試驗方法是將加速度傳感器固定在包裝件指定位置上，記錄運輸環節包裝件的加速度時域信號，然后將采集到的時域信號轉換成加速度功率譜密度。也可以將隨機振動記錄儀固定在車輛底板指定位置上，采集車輛底板的加速度時域信號，轉換成加速度功率譜密度后在隨機振動試驗臺上復現，而后利用隨機振動試驗臺獲取包裝件的振動響應和功率譜密度[1]。試驗方法操作簡單，精度高，但耗時長，成本較高。

仿真分析方法是首先建立路面不平度模型，生成路面對車輛的激勵；建立車輛的多自由度振動模型，得到車輛底板即包裝件底部的振動響應，該響應即為包裝件受到的振動激勵；根據實際情況建立包裝件的單自由度/多自由度或者非線性振動模型，最后得到其振動響應和功率譜密度[2]。仿真分析方法耗時非常短，精度相對較高，但對操作者的理論水平和數值計算程序編寫能力要求較高。基于仿真分析方法的以上優缺點，本文意圖清晰地梳理整個設計過程，并結合相應程序進行介紹。

1 路面不平度模型構建

1.1 路面頻域模型

國家標準GB 7031 建議采用路面頻域模型即功率譜密度來描述路面不平度，如下所示：

式中，和分別是參考空間頻率和空間頻率，單位是，，是譜密度指數，一般設置為2，和分別表示對應空間頻率下的路面位移譜密度，單位為，也稱為路面不平度系數。

為了方便地研究車輛與包裝件的振動特性，常常需要將路面不平度的頻域模型轉換成時域模型，其中常用的方法是濾波白噪聲生成法[3]。

1.2 路面時域模型

利用濾波白噪聲生成法可以建立路面不平度的時域模型，如下所示：

式中，是路面受到的隨機激勵，為車輛速度，單位是m/s，是零均值白噪聲過程，其方差，和是對應路面等級下的常數[4]。

基于上述分析，可基于Matlab 軟件的Simulink仿真模塊建立路面對四輪車輛輸入的時域仿真模型[5]，如圖1 所示。

圖1 路面對四輪車輛輸入的時域Simulink 仿真模型

2 四輪車輛振動模型構建

基于拉格朗日方程，可以建立多自由度運輸車輛的振動方程[2,4]：

其中，分別是系統的質量、剛度和阻尼矩陣，是車輛的位移響應列向量，是輪胎剛度矩陣，是路面對車輛前后四輪隨機激勵的列向量。

為了方便地使用Simulink 仿真模塊對上述模型進行求解，可以用狀態空間表達式對振動方程進行轉換，進而得到車輛底板的時域響應。假設車輛速度為15m/s，行駛在C 級路面上，仿真時間設置為50s，則四輪車輛底板的位移響應、加速度響應如圖2 所示，加速度均方根值為0.1739g。加速度功率譜密度可以使用Welch 平均周期圖法來進行繪制，在Matlab 軟件中的調用函數是pwelch( )，計算得到的加速度功率譜總均方根值為0.1739g。

圖2 四輪車輛底板的位移和加速度響應

3 包裝件振動模型構建

在上一部分，已經介紹如何獲得四輪車輛底板的時域響應。建立好包裝件振動模型，將四輪車輛底板的響應作為包裝件的輸入，可通過計算得到包裝件的振動響應。如果包裝件是單自由度/多自由度線性振動模型，可以用上一部分的狀態空間表達式對包裝件的振動方程進行處理。

如果包裝件是非線性振動系統，比如包裝件內緩沖材料是三次非線性或者包裝件內易損零部件需要看作是簡支梁式/懸臂梁式結構的情況，筆者建議包裝件的模型求解部分使用自寫程序的方式來解決。因為Simulink 仿真模型自帶的state-space 在continues模塊組下只能描述線性系統。接下來以二自由度非線性包裝件為研究對象介紹后續程序如何設計。

二自由度包裝件系統模型如圖3 所示，其中m1是易損零部件的質量，m2為除易損件外包裝件主體的質量。假設易損件與主體連接部分的彈性力是線性，包裝件主體與車輛底板之間緩沖材料的彈性力為三次非線性，表達式為F(x)=k2x+rx3，k2是初始彈性系數，r 是非線性常數，c1和c2是相應部分的阻尼，u 是四輪車輛底板的隨機位移激勵，已經在上一部分完成求解。對于此二自由度振動問題，可以以易損件靜平衡位置處為坐標原點，建立坐標系x1，以包裝件主體靜平衡位置處為坐標原點，建立坐標系x2，兩個坐標系都以向上為正方向，可建立系統的動力學方程：

圖3 二自由度包裝件系統模型

接下來可以用龍格庫塔數值算法對式4 的動力學方程進行求解，需要注意的是，由于系統的外部隨機激勵u和不能用一個關于時間的函數表示，因此需要對u和進行插值處理，找到任意時刻t 對應的值。程序設計如下所示：

y0=[0;0;0;0];%易損件與包裝件主體的初位移和初速度

tspan = u(:,1);%計算時間跨度

[TT, Y] = ode45(@ (t, y) odefun_2(t, y, k1,k2, c1, c2, m1, m2, r, u, du), tspan, y0);

function dydt = odefun_2(t, y, k1, k2, c1,c2, m1, m2, r, u, du)

dydt = zeros(4,1);

dydt(1) = y(3);

dydt(2) = y(4);

%對位移激勵做一個插值

Ut = interp1(u(:,1), u(:,2) , t, ‘spline’);

dut = interp1(u(:,1), du, t, ‘spline’);

dydt(3) = (k1/m1)＊(y(2) - y(1)) + (c1/m1)＊(y(4)- y(3));

dydt(4) = (k2/m2)＊(ut - y(2)) + r＊(ut -y(2))^3 + ...

(c2/m2)＊(dut - y(4)) - (k1/m2)＊(y(2) - y(1)) -(c1/m2)＊(y(4) - y(3));

end

運行上述代碼后，可以得到包裝件主體以及易損件的位移響應，如圖4 所示。

圖4 包裝件主體和易損件的位移響應

易損件和包裝件主體的加速度響應可通過下面代碼計算，結果如圖5 所示。計算各自的加速度均方根值分別是0.0037g 和0.0042g。

圖5 包裝件主體和易損件的加速度響應

%輸出每個時刻易損件的加速度ddXX1

ddXX1 = (k1/m1)＊(Y(:,2) - Y(:,1)) + (c1/m1)＊(Y(:,4) - Y(:,3));

ddXX1 = ddXX1/g;

disp ([‘易損件的加速度均方根值為 ‘ num2str(rms(ddXX1)) ‘ g’]);

dut = interp1 (u(:,1), du, TT, ‘spline’);

ut = interp1 (u(:,1), u(:,2) , TT, ‘spline’);

%輸出每個時刻包裝件的加速度 ddXX2

ddXX2 = (k2/m2)＊(ut - Y(:,2)) + r＊(ut -Y(:,2)).^3 + ...

(c2/m2)＊(dut - Y(:,4)) - (k1/m2)＊(Y(:,2) -Y(:,1) ) - (c1/m2)＊(Y(:,4) - Y(:,3));

ddXX2= ddXX2/g;

disp([‘包裝件加速度均方根值為 ‘ num2str(rms( ddXX2)) ‘ g’ ]);

使用平均周期圖法可以計算易損件的加速度功率譜密度，具體代碼如下。計算結果如圖6 所示。計算得到的加速度功率譜密度總均方根值為0.0038g，與車輛底板相比，隨機振動強度大幅減弱，說明緩沖包裝可以有效地保護易損零部件。

圖6 易損件的加速度功率譜密度

xn = ddXX1;

nfft = 1024;

noverlap = 20; %數據無重疊

window = boxcar(100); %矩形窗

[pwelchx,Fxx] = pwelch (xn, window, 20,nfft, Fs, ‘psd’);

plot(Fxx,pwelchx);

P S D_G r m s_p w e l c h x = C a l_P S D_Grms(Fxx,pwelchx); %計算功率譜密度的總均方根值

disp([‘平均周期法PSD 的總均方根值為 ‘num2str(PSD_Grms_pwelchx) ‘ g’]);

xlabel(‘頻率(Hz)’); ylabel(‘易損件加速度功率譜密度(g^2/Hz)’);

set(gca,’yscale’,’log’);%取對數實現縱坐標不均勻分布set(gca,’xscale’,’log’);%取對數實現縱坐標不均勻分布

xlim([0 200])

4 結語

本文較為詳細地從仿真角度介紹了如何快速計算公路運輸過程包裝件隨機振動時域響應和功率譜密度的方法，并且給出了相關程序。上述程序適用于非線性包裝件振動模型，對于易損件為簡支梁或懸臂梁結構，可以運用有限差分法對模型進行求解。由于篇幅有限，筆者在這里不過多贅述。希望本文可以為相關研究者提供一些參考，同時為包裝工程本科生學習運輸包裝課程時提供一些幫助，加深對包裝件隨機振動的認識。