包裝件隨機(jī)振動(dòng)時(shí)域響應(yīng)和PSD 快速獲取方法研究

薛陽，宋衛(wèi)生

（河南牧業(yè)經(jīng)濟(jì)學(xué)院，河南 鄭州 450046）

在公路運(yùn)輸環(huán)節(jié)，貨物（包裝件）受到的激勵(lì)主要是隨機(jī)振動(dòng)，來自運(yùn)輸工具（車輛）底板。要深入地分析包裝件或包裝件內(nèi)易損零部件的振動(dòng)行為，需要得到包裝件或易損件的振動(dòng)響應(yīng)和其功率譜密度(PSD)。目前有兩種方法可以實(shí)現(xiàn)，分別是從試驗(yàn)和仿真分析的角度出發(fā)。

試驗(yàn)方法是將加速度傳感器固定在包裝件指定位置上，記錄運(yùn)輸環(huán)節(jié)包裝件的加速度時(shí)域信號(hào)，然后將采集到的時(shí)域信號(hào)轉(zhuǎn)換成加速度功率譜密度。也可以將隨機(jī)振動(dòng)記錄儀固定在車輛底板指定位置上，采集車輛底板的加速度時(shí)域信號(hào)，轉(zhuǎn)換成加速度功率譜密度后在隨機(jī)振動(dòng)試驗(yàn)臺(tái)上復(fù)現(xiàn)，而后利用隨機(jī)振動(dòng)試驗(yàn)臺(tái)獲取包裝件的振動(dòng)響應(yīng)和功率譜密度[1]。試驗(yàn)方法操作簡單，精度高，但耗時(shí)長，成本較高。

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

1 路面不平度模型構(gòu)建

1.1 路面頻域模型

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

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

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

1.2 路面時(shí)域模型

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

式中，是路面受到的隨機(jī)激勵(lì)，為車輛速度，單位是m/s，是零均值白噪聲過程，其方差，和是對(duì)應(yīng)路面等級(jí)下的常數(shù)[4]。

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

圖1 路面對(duì)四輪車輛輸入的時(shí)域Simulink 仿真模型

2 四輪車輛振動(dòng)模型構(gòu)建

基于拉格朗日方程，可以建立多自由度運(yùn)輸車輛的振動(dòng)方程[2,4]：

其中，分別是系統(tǒng)的質(zhì)量、剛度和阻尼矩陣，是車輛的位移響應(yīng)列向量，是輪胎剛度矩陣，是路面對(duì)車輛前后四輪隨機(jī)激勵(lì)的列向量。

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

圖2 四輪車輛底板的位移和加速度響應(yīng)

3 包裝件振動(dòng)模型構(gòu)建

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

如果包裝件是非線性振動(dòng)系統(tǒng)，比如包裝件內(nèi)緩沖材料是三次非線性或者包裝件內(nèi)易損零部件需要看作是簡支梁式/懸臂梁式結(jié)構(gòu)的情況，筆者建議包裝件的模型求解部分使用自寫程序的方式來解決。因?yàn)镾imulink 仿真模型自帶的state-space 在continues模塊組下只能描述線性系統(tǒng)。接下來以二自由度非線性包裝件為研究對(duì)象介紹后續(xù)程序如何設(shè)計(jì)。

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

圖3 二自由度包裝件系統(tǒng)模型

接下來可以用龍格庫塔數(shù)值算法對(duì)式4 的動(dòng)力學(xué)方程進(jìn)行求解，需要注意的是，由于系統(tǒng)的外部隨機(jī)激勵(lì)u和不能用一個(gè)關(guān)于時(shí)間的函數(shù)表示，因此需要對(duì)u和進(jìn)行插值處理，找到任意時(shí)刻t 對(duì)應(yīng)的值。程序設(shè)計(jì)如下所示：

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

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

[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);

%對(duì)位移激勵(lì)做一個(gè)插值

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

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

圖4 包裝件主體和易損件的位移響應(yīng)

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

圖5 包裝件主體和易損件的加速度響應(yīng)

%輸出每個(gè)時(shí)刻易損件的加速度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’);

%輸出每個(gè)時(shí)刻包裝件的加速度 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’ ]);

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

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

xn = ddXX1;

nfft = 1024;

noverlap = 20; %數(shù)據(jù)無重疊

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); %計(jì)算功率譜密度的總均方根值

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

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

set(gca,’yscale’,’log’);%取對(duì)數(shù)實(shí)現(xiàn)縱坐標(biāo)不均勻分布set(gca,’xscale’,’log’);%取對(duì)數(shù)實(shí)現(xiàn)縱坐標(biāo)不均勻分布

xlim([0 200])

4 結(jié)語

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