彈性連桿機(jī)構(gòu)的漸變動力學(xué)及可靠性靈敏度

劉曉靜

（山西工程職業(yè)學(xué)院汽車與軌道交通工程系，山西 太原 030009）

1 前言

早期，運動機(jī)構(gòu)轉(zhuǎn)速較低、剛性較大，因此可以忽略機(jī)構(gòu)彈性變形對其工作性能的影響。然而隨著現(xiàn)代工業(yè)的發(fā)展，機(jī)構(gòu)速度高、構(gòu)件質(zhì)量輕的趨勢使得機(jī)構(gòu)呈現(xiàn)柔性化，導(dǎo)致彈性機(jī)構(gòu)的彈性變形增大，使機(jī)構(gòu)實際的運動偏離了預(yù)期運動。與此同時，彈性變形會產(chǎn)生較大的動應(yīng)力，進(jìn)而導(dǎo)致機(jī)構(gòu)疲勞破壞。

因而研究彈性機(jī)構(gòu)的漸變動力學(xué)及可靠性靈敏度就顯得尤為重要。文獻(xiàn)［1］最早應(yīng)用偽剛體模型方法進(jìn)行分析，對彈性四桿機(jī)構(gòu)進(jìn)行運動綜合。文獻(xiàn)［2］應(yīng)用有限元分析方法探索了機(jī)構(gòu)的各階固有頻率均值和轉(zhuǎn)速在低階下的關(guān)系。文獻(xiàn)［3］則提出了機(jī)械動態(tài)與漸變可靠性理論，基于首次穿越理論的可靠性分析研究現(xiàn)狀進(jìn)行了詳細(xì)的闡述。但目前國內(nèi)外計算動應(yīng)力的方法比較復(fù)雜［4］，存在計算工作量大，工作效率低的問題。1989年，文獻(xiàn)［5］將Kriging理論推廣到確定性計算機(jī)實驗的設(shè)計和分析領(lǐng)域，并給出了一種較實用的Kriging 算法（或Kriging 模型）。Kriging 模型根據(jù)相關(guān)函數(shù)使其具有局部估計的特點，是一種保證方差最小的無偏估計的模型，因此它受到近似確定計算機(jī)模型越來越多的青睞［6］。因此，文中采用Kriging替代模型，使其具有較好的擬合度，同時為后續(xù)計算可靠性靈敏度較大程度地降低計算量。

這里以兩節(jié)點八自由度梁單元為基礎(chǔ)，建立彈性連桿機(jī)構(gòu)的動力學(xué)模型，獲得具有漸變參數(shù)的系統(tǒng)動力學(xué)方程，進(jìn)而建立Kriging近似替代模型。最后通過算例，考慮構(gòu)件彈性變形求解其動力學(xué)特性和可靠性靈敏度，使計算結(jié)果更加接近工程實際。利用響應(yīng)面近似替代模型可以較大程度地提高計算效率，通過修正構(gòu)件參數(shù)可以優(yōu)化彈性連桿機(jī)構(gòu)的可靠性。

2 彈性連桿機(jī)構(gòu)動力學(xué)建模

2.1 梁單元運動學(xué)關(guān)系

選取兩節(jié)點八自由度等截面梁為研究對象，單元內(nèi)的縱向位移（SW）和軸向位移（SV）分別滿足五次Hermit多項式和線性分布，任意一點的縱向位移（SW）和軸向位移（SV）都可通過“型函數(shù)”和單元坐標(biāo)系中的廣義坐標(biāo)表示為：

式中：u1，u5—節(jié)點軸向位移；u2，u6—節(jié)點縱向位移；u3，u7—節(jié)點彈性轉(zhuǎn)角；u4，u8—節(jié)點曲率；?i—型函數(shù)。

這里在忽略剛?cè)狁詈享椀那疤嵯拢芯苛簡卧\動學(xué)關(guān)系［7］。假定機(jī)構(gòu)的剛體（加）速度和彈性（加）速度的矢量和為實際（加）速度。

結(jié)合運動關(guān)系及式（1），可以得到梁單元軸線上坐標(biāo)x處的速度為：

2.2 梁單元微分方程

為了得到梁單元的動力學(xué)方程，這里需要引入單元的動能、變形能和Lagrange方程［8］。

根據(jù)材料力學(xué)和式（2），分別得到單元動能和變形能的矩陣形式為：

將式（3）、式（4）代入式（5）中，整理得到梁單元運動微分方程式為：

根據(jù)機(jī)構(gòu)運動學(xué)假定，式（6）可改寫為：

2.3 系統(tǒng)微分方程

為使單元運動方程組集成系統(tǒng)運動方程，這里引入一組整體坐標(biāo)系中的單元廣義坐標(biāo)［10］。

它與單元坐標(biāo)系中的廣義坐標(biāo)列陣關(guān)系為：

式中：R—坐標(biāo)轉(zhuǎn)換矩陣。

將式（9）代入式（7）中，得到整體坐標(biāo)系中的梁單元運動微分方程為：

定義系統(tǒng)廣義坐標(biāo)列陣為：U={U1，...，UN}T

式中：N—彈性力學(xué)分析的自由度數(shù)。

坐標(biāo)協(xié)調(diào)矩陣Bi可將單元編號Uie和系統(tǒng)編號U實現(xiàn)轉(zhuǎn)換，表示為：

將式（11）代入式（10）得到以系統(tǒng)坐標(biāo)為未知量的單元運動微分方程，再將各單元對應(yīng)坐標(biāo)系下的運動微分方程進(jìn)行疊加，得到系統(tǒng)運動方程為：

當(dāng)考慮阻尼時，系統(tǒng)運行方程為：

其中，阻尼矩陣C常取比例阻尼形式，即：

式中：α、β—材料的比例系數(shù)，一般由前兩階固有頻率確定。

3 彈性連桿機(jī)構(gòu)漸變動力學(xué)特性

3.1 固有頻率和特性

由于阻尼對機(jī)構(gòu)的固有頻率影響非常小，根據(jù)無阻尼自由振動方程，可以得到機(jī)構(gòu)的固有頻率［11］和模態(tài)。

根據(jù)機(jī)構(gòu)的自由振動得到下列等式：

上式存在非零解的條件是系數(shù)矩陣行列式為0，即：由此得到特征值λ=ω2，將其代入式（15）中得到對應(yīng)的特征向量（模態(tài)）X。其中X()

i對應(yīng)特征值，i= 1，2，3，…，n。

3.2 動態(tài)應(yīng)力應(yīng)變分析

梁單元任意一點的彎曲應(yīng)變?yōu)椋?/p>

式中：ξ—在y方向上的坐標(biāo)；S''W—單元縱向位移對于坐標(biāo)x的二次導(dǎo)數(shù)。

梁單元中，線性單元的拉伸應(yīng)變［12］為：

式中：單元軸向和縱向位移對于坐標(biāo)軸變化率取平方引起的附加拉壓應(yīng)變。

同理求得拉壓應(yīng)變后，梁單元任一點的應(yīng)變是彎曲應(yīng)變和拉壓應(yīng)變的疊加，由此得到總體應(yīng)變?yōu)椋?/p>

矩陣；U通過u=RiBiU得到。

此過程中沒有考慮剪切應(yīng)變的影響，因為一般情況下這種影響是微小的，可以忽略。

根據(jù)σ=Eε，當(dāng)柔順桿件彎曲時，任意截面外邊緣的彎曲應(yīng)力為：

任意截面的拉壓應(yīng)力為：

式（22）中可以看出單元中拉壓應(yīng)力為常值。因此單元中外邊緣任意一點的應(yīng)力為：

所有截面中最大動應(yīng)力的絕對值為：

任意一個位置的動應(yīng)力可由式（24）得到，按公式計算可得不同時間不同位置的單元動應(yīng)力。

4 彈性連桿機(jī)構(gòu)的可靠性靈敏度

4.1 Kriging模型

假設(shè)給定m個樣本點及其響應(yīng)，y表示響應(yīng)的預(yù)測值，F(xiàn)(β：，l，x)表示回歸模型，zl(x)表示隨機(jī)函數(shù)（隨機(jī)過程），則響應(yīng)的預(yù)測值模型［13］為：

隨機(jī)部分zl(x)滿足以下統(tǒng)計特征：

定義隨機(jī)部分任意樣本點間的相關(guān)方程為：

待測點與樣本點的相關(guān)向量中的元素為：任意預(yù)測點x處的響應(yīng)可通過樣本點線性組合得到，即：

為了保證擬合的無偏差性，預(yù)測誤差均值應(yīng)為0、方差為極小值，進(jìn)而得到系數(shù)C。

因為Fβ?Y，則模擬的最小二乘估計為：

記殘差為Rγ*=Y-Fβ，可得：

由式（31）可知，只要給定f(x)，r(x)便可得到任意待測點的預(yù)測響應(yīng)。

預(yù)測誤差的方差為：

式中：y(x)—響應(yīng)的真實值。

為保證方差最小，需要計算不確定的參數(shù)σ和相關(guān)系數(shù)R中的θ，則其對數(shù)似然函數(shù)為：

ln(P(x))

則問題轉(zhuǎn)化為求解最大化無約束問題：

4.2 Kriging響應(yīng)面構(gòu)建中常用的試驗設(shè)計方法

響應(yīng)面的構(gòu)造，首先是初始樣本點的選取，可以通過隨機(jī)采樣獲取，也可以通過優(yōu)化試驗設(shè)計進(jìn)行采樣，最終使初始樣本點在樣本空間內(nèi)均勻分布，以便獲得良好的近似Kriging替代模型。為了獲得良好的數(shù)據(jù)，常用的試驗設(shè)計方法有：全因子試驗設(shè)計、部分因子試驗設(shè)計、中心復(fù)合試驗設(shè)計、拉丁方設(shè)計和序列探索實驗設(shè)計方法等。在眾多試驗方法中，拉丁超立方（LHD）是使用較為廣泛的方法，下面對拉丁超立方進(jìn)行一個簡單介紹。

假設(shè)在n維空間內(nèi)抽取m個樣本點，拉丁超立方抽（LDS）的步驟為：

（1）將n維中每一維劃分為互不重疊的m個區(qū)間，使得每個區(qū)間有相同的概率（通常考慮到均勻分布，每個區(qū)間的長度相同）。

（2）在不同維不同區(qū)間中隨機(jī)抽取一個點。

（3）再從每個區(qū)間里隨機(jī)抽?。?）中選取的任一維點，并將它們組成向量，得到初始樣本點。

4.3 擬合優(yōu)度（決定系數(shù)）

響應(yīng)面擬合結(jié)果可以通過擬合優(yōu)度（決定系數(shù)）來判定。

假定Y是m個樣本點的真實響應(yīng)值，Y是真實響應(yīng)值的均值，y是樣本點通過Kriging模型得到的預(yù)測響應(yīng)值?？傠x平方和（SST）與殘差平方和（SSE）分別為：

則擬合優(yōu)度：

4.4 失效概率

假定結(jié)構(gòu)的功能函數(shù)為Z=gX(X)，隨機(jī)變量X的概率密度函數(shù)為fX(x)。按fX(x)對X進(jìn)行隨機(jī)抽樣，所得的樣本代入功能函數(shù)Z=gX(X)，Z＜0則記為失效一次。設(shè)模擬次數(shù)為N，其中Z＜0的次數(shù)為nf次。根據(jù)概率論總的大數(shù)據(jù)定律中的Bernoulli定理［14］知，N次獨立試驗中的失效概率 收斂于該事件的概率nf N，pf則失效概率［14］估計值為：

結(jié)構(gòu)的失效概率可以表示為：

式中：I(x)—x的指示函數(shù)。規(guī)定x＜0時，I(x) = 1；x≥0，I(x) =0。根據(jù)式（39），X的第i個樣本值為xi，則pf的估計值為：

4.5 可靠性靈敏度

靈敏度［15］為失效概率pf對基本變量xi的分布參數(shù)θ(k)x i(i=1，2，...，n；k= 1，2，...，m)，其 中mi為第i個變量xi分布參數(shù)總數(shù)的偏導(dǎo)數(shù)，將失效概率的積分對分布參數(shù)θ(k)x i求導(dǎo)數(shù)，則得到可靠性靈敏度：

5 算例

5.1 彈性連桿機(jī)構(gòu)[17]動力學(xué)分析

平面彈性連桿機(jī)構(gòu)的主要參數(shù)，主動桿的轉(zhuǎn)速為300r/min，如表1所示。將彈性四桿機(jī)構(gòu)的模型結(jié)構(gòu)，每個部件劃分為4個單元，每個單元節(jié)點有4個節(jié)點廣義坐標(biāo)，如圖1所示。

圖1 彈性四桿機(jī)構(gòu)模型結(jié)構(gòu)Fig.1 The Module of Elastic Four-Bar Mechanism

表1 四桿機(jī)構(gòu)參數(shù)Tab.1 Parameters of the Four-Bar Mechanism

在該轉(zhuǎn)速下系統(tǒng)的前三階固有頻率，如圖2所示。

圖2 系統(tǒng)前三階固有頻率

表示主動桿的轉(zhuǎn)速分別為300r/min、600r/min和1200r/min不同時刻機(jī)構(gòu)對應(yīng)的最大應(yīng)力，如圖3所示。隨著轉(zhuǎn)速的增加，機(jī)構(gòu)的最大動應(yīng)力隨之增加。

圖3 系統(tǒng)最大動應(yīng)力Fig.3 The System Maximum Dynamic Stress

5.2 彈性連桿機(jī)構(gòu)的可靠性靈敏度

首先通過拉丁超立方［18］抽取50個樣本點，根據(jù)動力學(xué)特性得到結(jié)構(gòu)真實響應(yīng)值。根據(jù)樣本點及其響應(yīng)建立Kriging近似替代模型。另取誤差檢驗參考點100個。Kriging近似替代模型的預(yù)測值和真實值的擬合效果的對比，從圖中可知有較好的擬合結(jié)果，同時得到擬合優(yōu)度（決定系數(shù)）為：R2= 0.9852，如圖4所示。

圖4 Kriging模型擬合效果Fig.4 Model Fitting Effect of Kriging

主動桿轉(zhuǎn)速為300r/min，在Kriging模型的基礎(chǔ)上，利用拉丁超立方抽樣方式進(jìn)行抽樣，設(shè)定抽樣模擬次數(shù)為5000次，最大應(yīng)力隨著模擬次數(shù)的增加的變化曲線，如圖5所示。

圖5 基于Kriging模型最大應(yīng)力隨抽樣模擬次數(shù)的變化曲線Fig.5 The Maximum Stress Changes with the Increase of the Simulation Curve

由圖5可以看出：隨著抽樣模擬次數(shù)的增加，最大應(yīng)力抽樣結(jié)果趨于平穩(wěn)收斂，在最大應(yīng)力為261.06MPa處收斂，該圖表明近似替代模型可以準(zhǔn)確計算機(jī)構(gòu)實際可靠度。

彈性連桿機(jī)構(gòu)桿件尺寸存在隨機(jī)性，具體參數(shù)，如表2所示。選?。?3σ，3σ］為樣本采集范圍。

表2 機(jī)構(gòu)變量參數(shù)表Tab.2 Variable Parameter List

主動桿件轉(zhuǎn)速為300r/min，根據(jù)表2變量參數(shù)表利用拉丁超立方抽樣1000次計算曲柄中點彈性縱向位移隨著主動桿件轉(zhuǎn)動的失效概率，最大的失效概率為0.427 出現(xiàn)在（0～60）°內(nèi)，如圖6所示。

圖6 曲柄中點彈性位移失效概率Fig.6 The Failure of Crank Midpoint

曲柄中點彈性縱向位移對于各桿件參數(shù)的均值可靠性靈敏度和標(biāo)準(zhǔn)差靈敏度，如圖7、圖8所示。由圖中可以看出在同一時刻，不同參數(shù)的均值和標(biāo)準(zhǔn)差對失效概率會產(chǎn)生不同的影響，即靈敏度值有正有負(fù)。在不同時刻，相同參數(shù)的均值和標(biāo)準(zhǔn)差對失效概率產(chǎn)生相同的影響，但是對失效概率的影響程度即靈敏度值會隨著主動桿件位置的改變而發(fā)生相應(yīng)的改變，系統(tǒng)失效概率最大時各參數(shù)的靈敏度達(dá)到極值。

圖7 曲柄中點彈性位移各參數(shù)的均值靈敏度Fig.7 Mean Sensitivity of Each Crank Midpoint Elastic Displacement Parameters

圖8 曲柄中點彈性位移各參數(shù)的標(biāo)準(zhǔn)差靈敏度Fig.8 Standard Deviation Sensitivity of Each Crank Midpoint Elastic Displacement Parameters

靈敏度直方圖，如圖9、圖10所示。

圖9 均值靈敏度直方圖Fig.9 Histogram of Mean Sensitivity

圖10 標(biāo)準(zhǔn)差靈敏度直方圖Fig.10 Histogram of Standard Deviation Sensitivity

從圖中可以得到桿件L2的均值和標(biāo)準(zhǔn)差對系統(tǒng)可靠性都具有最大的影響，L2的均值靈敏度值為?Pf/?μ2=-25.694，標(biāo)準(zhǔn)差靈敏度值為?Pf/?σ2=-4.89648，由此可見靈敏度對于可靠性的影響是消極的。

而均值產(chǎn)生積極影響的是參數(shù)L1，標(biāo)準(zhǔn)差中產(chǎn)生積極影響的是參數(shù)L4。

6 結(jié)束語

（1）這里以梁單元為基礎(chǔ)，利用廣義坐標(biāo)和Lagrange方程建立單元運動微分方程。通過坐標(biāo)轉(zhuǎn)換矩陣和協(xié)調(diào)矩陣，單元運動微分方程組集得到系統(tǒng)運動微分方程。利用動力學(xué)模型對彈性連桿機(jī)構(gòu)的漸變動力學(xué)特性進(jìn)行分析。

（2）利用拉丁超立方抽樣，對比真實值和Kriging響應(yīng)面近似替代模型得到的預(yù)測值，分析可知響應(yīng)近似替代模型具有良好的擬合度。

因此在計算彈性連桿機(jī)構(gòu)的漸變動力學(xué)和可靠性靈敏度時，基于近似替代模型，可以較大程度地提高計算效率。

（3）以平面彈性連桿機(jī)構(gòu)為例，進(jìn)一步證實Kriging近似替代模型具有良好的擬合度。隨著主動桿轉(zhuǎn)速的提高，機(jī)構(gòu)任意時刻的最大動應(yīng)力隨之增大，因此控制機(jī)構(gòu)的動應(yīng)力可以通過控制主動桿件的轉(zhuǎn)速實現(xiàn)。

（4）基于Kriging替代模型，平面彈性四桿機(jī)構(gòu)曲柄中點的彈性縱向位移的失效概率最大值出現(xiàn)在（0～60）°內(nèi)，最大的失效概率為0.427。

桿件L2的均值和標(biāo)準(zhǔn)差對系統(tǒng)可靠性均具有最大且消極的影響。而只有桿件L1的均值及L4的標(biāo)準(zhǔn)差對系統(tǒng)的可靠性產(chǎn)生積極影響。