基于全局動態(tài)應(yīng)力解空間譜單元插值的關(guān)鍵時間點識別

張艷崗　毛虎平　蘇鐵熊　李　坤　王　軍

中北大學(xué),太原,030051

?

基于全局動態(tài)應(yīng)力解空間譜單元插值的關(guān)鍵時間點識別

張艷崗毛虎平蘇鐵熊李坤王軍

中北大學(xué),太原,030051

摘要：針對結(jié)構(gòu)動態(tài)響應(yīng)優(yōu)化中動態(tài)分析的復(fù)雜性與高耗時性問題,提出了基于全局動態(tài)應(yīng)力解空間譜單元插值的關(guān)鍵點識別方法，找到結(jié)構(gòu)動態(tài)響應(yīng)下最危險的時刻。首先利用模態(tài)疊加法，獲得結(jié)構(gòu)的模態(tài)應(yīng)力分布，并計算全局動態(tài)應(yīng)力解空間，然后利用譜單元離散動態(tài)應(yīng)力絕對極大值點曲線，采用Lagrange插值并調(diào)用區(qū)域細分全局優(yōu)化求解器，找到全局動態(tài)應(yīng)力的極大值與極小值，即關(guān)鍵時間點。將該方法應(yīng)用到124桿平面桁架和均布力與集中力共同作用的結(jié)構(gòu)關(guān)鍵時間點識別問題中，識別結(jié)果說明了該方法的可行性和有效性。

關(guān)鍵詞：解空間譜單元離散；全局動態(tài)應(yīng)力；Lagrange插值；關(guān)鍵時間點

0引言

工程機械結(jié)構(gòu)幾乎都在動態(tài)載荷作用下工作，其各種性能均表現(xiàn)為時間的函數(shù)。為了使機械結(jié)構(gòu)的動態(tài)性能達到極致，確保其在動載環(huán)境下可靠工作,目前最行之有效的辦法是進行結(jié)構(gòu)動態(tài)響應(yīng)優(yōu)化設(shè)計。然而,在動態(tài)響應(yīng)優(yōu)化過程中，優(yōu)化器需要反復(fù)調(diào)用與設(shè)計變量相關(guān)的目標函數(shù)和約束函數(shù)，最直接的方法是通過結(jié)構(gòu)動態(tài)分析獲得相關(guān)數(shù)據(jù)，然而結(jié)構(gòu)瞬態(tài)動力學(xué)分析非常復(fù)雜和耗時[1-2],因此動態(tài)響應(yīng)優(yōu)化設(shè)計存在很大困難，有時是不太可能實現(xiàn)的。

動態(tài)優(yōu)化設(shè)計的參數(shù)變化使得應(yīng)力解空間在整個時間域和結(jié)構(gòu)空間域上變化,即不論是目標函數(shù)還是約束函數(shù)在隨時間變化的同時，其空間位置也在發(fā)生變化。與此同時，結(jié)構(gòu)靜態(tài)優(yōu)化技術(shù)已經(jīng)趨于成熟，在此背景下將動態(tài)優(yōu)化問題轉(zhuǎn)化為靜態(tài)優(yōu)化問題成為了可能。既然要轉(zhuǎn)化，那么首先要解決的問題是在哪個時刻轉(zhuǎn)化，即等效轉(zhuǎn)化為結(jié)構(gòu)在哪個時刻的動態(tài)響應(yīng)，這是動態(tài)優(yōu)化轉(zhuǎn)化等效靜態(tài)優(yōu)化的關(guān)鍵一步。文獻[3]研究了3種識別關(guān)鍵時間點的方法:第一種方法采用自適應(yīng)搜索技術(shù)，在3點或5點擬合二次函數(shù)，然后運用優(yōu)化方法得到極值點，該方法適合于響應(yīng)比較緩慢的情況；第二種方法利用最小二乘拋物線樣條函數(shù)來近似結(jié)構(gòu)的動態(tài)響應(yīng)，該方法可以有效近似噪聲動態(tài)響應(yīng)；第三種方法是超峰值法，采用包括設(shè)計變量倒數(shù)的一階泰勒展開來近似約束函數(shù)。此外,文獻[3]中的關(guān)鍵時間點是針對直接求解動態(tài)優(yōu)化問題提出的，僅對與時間相關(guān)的約束函數(shù)進行了近似，并利用近似函數(shù)獲得關(guān)鍵時間點。這種方法滿足了與時間相關(guān)的約束函數(shù)，認為時間域上的所有時間點都滿足約束，而沒有考慮目標函數(shù)通常也是與時間相關(guān)的函數(shù)。在結(jié)構(gòu)動態(tài)優(yōu)化過程中，對目標函數(shù)的求解也需要進行動態(tài)響應(yīng)分析，因此沒有從根本上解決問題。

文獻[4]采用動態(tài)位移響應(yīng)對時間的微分等于零來確定關(guān)鍵時間點，然而微分求解比較困難，并沒有詳細研究求解方法；而且，即便已經(jīng)確定關(guān)鍵時間點，即動態(tài)位移響應(yīng)的極值點為關(guān)鍵時間點，雖然從一定意義上可以說明結(jié)構(gòu)在動態(tài)載荷作用下位移響應(yīng)達到極值的時刻最危險，然而位移達到極值的時刻應(yīng)力不一定也達到極值時刻[5]；再者目前的研究文獻中僅僅以結(jié)構(gòu)某一個節(jié)點的動態(tài)位移絕對值最大為結(jié)構(gòu)關(guān)鍵時間點,沒有從結(jié)構(gòu)全局考慮最危險的時刻。因此,將某一節(jié)點動態(tài)位移響應(yīng)達到極值的時刻作為結(jié)構(gòu)最危險的時刻來進行結(jié)構(gòu)設(shè)計，結(jié)構(gòu)可能欠設(shè)計或過設(shè)計。

課題組在前期研究中發(fā)現(xiàn),利用譜單元離散插值一個未知函數(shù)的精度非常高[6-7],并在系統(tǒng)動力學(xué)分析中得到應(yīng)用驗證,可是對于零部件結(jié)構(gòu)的進一步設(shè)計需要借助有限元軟件來實現(xiàn)。筆者前期在研究動態(tài)優(yōu)化時，提出了基于動態(tài)應(yīng)力解空間譜元離散的關(guān)鍵時間點識別方法[8],旨在找到結(jié)構(gòu)受動態(tài)載荷作用最危險時刻，為下一步等效靜態(tài)轉(zhuǎn)化提供轉(zhuǎn)化時間點;筆者在進一步研究中發(fā)現(xiàn)結(jié)構(gòu)最危險的時刻并不是某一個單元的應(yīng)力達到極值的時刻,而應(yīng)該是結(jié)構(gòu)全局動態(tài)應(yīng)力達到極值的時刻。全局動態(tài)應(yīng)力是結(jié)構(gòu)在動態(tài)載荷的作用下，結(jié)構(gòu)所有單元中的應(yīng)力最大值,并且最大應(yīng)力隨時間變化，最大應(yīng)力發(fā)生的單元也隨著時間變化。因此，本文從結(jié)構(gòu)全局動態(tài)應(yīng)力的角度來確定關(guān)鍵時間點，提出全局動態(tài)應(yīng)力解空間譜單元插值的關(guān)鍵時間點識別方法。

1譜單元插值

對未知函數(shù)的插值,學(xué)者們一直不斷探索新的方法以追求高精度,但往往由于被插值函數(shù)的復(fù)雜性，插值效果并不理想。研究中發(fā)現(xiàn)在正交多項式的零點插值時,能夠獲得高精度插值函數(shù)[9]。

(1)

Legendre正交多項式是Sturm-Liouville form方程的解：

[(1-ξ2)P′(ξ)]′+n(n+1)P(ξ)=0

(2)

式中,n為Legendre正交多項式的階數(shù)。

它的權(quán)函數(shù)ω=1。k階Legendre多項式可以定義如下：

(3)

由于Legendre多項式的零點不包括區(qū)間端點,因此引入Lobbato多項式。Lobbato多項式是通過Legendre多項式的微分定義的,即

(4)

其滿足正交特性:

(5)

式中,δi j為Kronecher函數(shù)。

GLL點定義如下：

(6)

式(6)的解為GLL點，其中,N為插值次數(shù)。

求解式(6)可以獲得GLL點坐標值,16次插值GLL點的分布如圖1所示。對龍格函數(shù)分別進行平均離散插值和GLL點離散插值，可以清楚地看到GLL點離散插值的高精度優(yōu)勢,見圖2。

圖1　16次插值GLL點分布

圖2　平均離散與GLL 點離散插值精度比較

2全局動態(tài)應(yīng)力

承受動態(tài)載荷的結(jié)構(gòu)響應(yīng)都是時間的函數(shù)，動態(tài)應(yīng)力也不例外。對于成熟的有限元法來說，結(jié)構(gòu)動態(tài)應(yīng)力不僅是時間的函數(shù)還是位置的函數(shù)，即結(jié)構(gòu)的每一個空間位置點的動態(tài)應(yīng)力均是時間的函數(shù)，變化趨勢基本一致，但由于受動態(tài)載荷的結(jié)構(gòu)是振動的，結(jié)構(gòu)每一空間位置的動態(tài)應(yīng)力達到極值的時刻不一定相同。本文在前期研究基礎(chǔ)上提出更加合理的基于全局動態(tài)應(yīng)力解空間關(guān)鍵時間點識別方法。

結(jié)構(gòu)在動態(tài)載荷作用下，其任一點i的動態(tài)應(yīng)力滿足[10]:

(7)

全局動態(tài)應(yīng)力對應(yīng)的關(guān)鍵時間點可以表示為

tCTP(|σ|G)=arg max|σ(i)(t)|

(8)

3關(guān)鍵時間點識別

結(jié)構(gòu)動態(tài)優(yōu)化需要反復(fù)進行動態(tài)分析，每迭代一次需要n+1次動態(tài)分析[11]。因此,基于等效原理將動態(tài)載荷轉(zhuǎn)化為靜態(tài)載荷，利用結(jié)構(gòu)靜態(tài)優(yōu)化代替動態(tài)優(yōu)化，可以有效提高動態(tài)優(yōu)化效率。本文提出基于全局動態(tài)應(yīng)力解空間譜單元插值的關(guān)鍵點識別方法，其流程見圖3，具體步驟如下：

圖3　基于全局動態(tài)應(yīng)力的關(guān)鍵時間點識別流程

(1)建立有限元動力模型。運動微分方程可表達為

(9)

M=ΦTmΦC=ΦTc Φ

K=ΦTkΦF=ΦTf(t)

(2)動力學(xué)分析。采用模態(tài)疊加法,先求解式(9)，獲得q,然后得到結(jié)構(gòu)的節(jié)點位移x=qΦ,再由公式σk(t)=DBkqkΦ獲得第k階模態(tài)應(yīng)力,通過式(7)得到結(jié)構(gòu)動態(tài)應(yīng)力,最后應(yīng)用式(8)獲得結(jié)構(gòu)全局動態(tài)應(yīng)力。

(3)全局動態(tài)應(yīng)力解空間譜單元離散。利用式(2)～式(6)將全局動態(tài)應(yīng)力進行解空間譜單元離散，并計算GLL點對應(yīng)的結(jié)構(gòu)全局動態(tài)應(yīng)力的值|σ|G(tGLL)。

(4)Lagrange插值結(jié)構(gòu)全局動態(tài)應(yīng)力。Lagrange插值函數(shù)為

(10)

其中,Φk+1(ξ)=(ξ-ξ1)(ξ-ξ2)…(ξ-ξk+1)為GLL點對應(yīng)的全局動態(tài)應(yīng)力插值基函數(shù),GLL點對應(yīng)的全局動態(tài)應(yīng)力值等于步驟(3)中獲得的值；Lk(ξ)為全局動態(tài)應(yīng)力在GLL點上Lagrange插值函數(shù)。

(5)識別時間關(guān)鍵點。由于結(jié)構(gòu)動態(tài)應(yīng)力不僅是時間的函數(shù)還是位置的函數(shù)，即結(jié)構(gòu)的每一個空間位置點的動態(tài)應(yīng)力均是時間的函數(shù)，因此全局動態(tài)應(yīng)力解空間函數(shù)特性極為復(fù)雜，一般為非線性多峰值函數(shù)。結(jié)合前期優(yōu)化設(shè)計經(jīng)驗，文中采用區(qū)域細分算法(DIRECT)計算結(jié)構(gòu)全局動態(tài)應(yīng)力的絕對極大值點。DIRECT算法是一種全局優(yōu)化算法，該算法既不需要確定目標函數(shù)的具體方程，同時也不要求相關(guān)的導(dǎo)數(shù)信息，而是自動根據(jù)優(yōu)化迭代過程中的每個采樣點處的目標函數(shù)估值和超矩形的區(qū)域細分情況選擇搜索方向，因而非常適合于黑盒函數(shù)的仿真優(yōu)化，而且它能夠保證收斂到全局最優(yōu)點。采用的基于全局動態(tài)應(yīng)力關(guān)鍵時間點識別的數(shù)學(xué)模型如下：

(11)

其中，M為全局動態(tài)應(yīng)力解空間離散的單元數(shù)；k為每一個單元離散GLL點數(shù)；|σ(i)(tk)|max為第i個單元的第k個GLL點對應(yīng)的全局動態(tài)應(yīng)力；t0為仿真開始時間,通常為0;t1為仿真結(jié)束時間。

4應(yīng)用驗證

4.1124桿平面桁架

(a)124桿平面桁架

(b)動態(tài)載荷圖4　124桿平面桁架及動態(tài)載荷

如圖4所示,該平面桁架結(jié)構(gòu)有49個鉸鏈,94個自由度。彈性模量E=207 GPa,泊松比ν=0.3,密度ρ=7850 kg/m3,桿的截面積為64.5 mm2。動態(tài)載荷為半正弦函數(shù)(圖4)。在節(jié)點1、20、19、18、17、16、15的x正方向作用同樣大的動態(tài)載荷,在節(jié)點1、2、3、4、5的y負方向也作用同樣大的動態(tài)載荷。

圖5所示為124桿平面桁架全局動態(tài)應(yīng)力譜單元離散插值結(jié)果，從圖中可以看出結(jié)構(gòu)的全局動態(tài)應(yīng)力非常復(fù)雜，是高度非線性多峰值函數(shù)。譜單元Lagrange插值全局動態(tài)應(yīng)力精度非常高，應(yīng)用DIRECT全局優(yōu)化器準確地找到了關(guān)鍵時間點。

圖5　124桿平面桁架全局動態(tài)應(yīng)力

圖6　單元15的軸向應(yīng)力時間歷程

圖7　節(jié)點15的x方向位移時間歷程

如果從某一節(jié)點的動態(tài)位移或者動態(tài)應(yīng)力出發(fā)，尋找關(guān)鍵時間點，當然可以找到對應(yīng)的“關(guān)鍵時間點”,見圖6、圖7。將不同識別方法得到的關(guān)鍵時間點匯總至表1，可以發(fā)現(xiàn)：基于位移識別方法獲得的兩點為t1=0.096 698 411 s,t2=0.192 408 258 s，t1時刻節(jié)點15x方向正位移達到最大,t2時刻節(jié)點15x方向負位移達到最大；基于局部應(yīng)力識別方法獲得的兩點為t1=0.106 521 572 s,t2=0.192 519 952 s,t1時刻單元15拉應(yīng)力達到最大,t2時刻單元15壓應(yīng)力達到最大;基于全局動態(tài)應(yīng)力識別方法也獲得兩點，即t1=0.102 063 859 s，t2=0.192 058 523 s,t1時刻全局某一點拉應(yīng)力達到最大,t2時刻全局某一點壓應(yīng)力達到最大。位移關(guān)鍵點與局部應(yīng)力關(guān)鍵點不一致,最大相差0.009 823 161 s,也就是說位移達到最大時,局部應(yīng)力還沒有達到最大，存在一個時間差；而全局應(yīng)力關(guān)鍵點與局部應(yīng)力關(guān)鍵點也不一致，最大相差0.004 457 713 s,說明局部應(yīng)力達到最大時，全局應(yīng)力還沒有達到最大，存在一個時間差。

表1　124桿平面桁架關(guān)鍵時間點識別結(jié)果

4.2均布力與集中力作用的平面桁架

圖8所示為均勻分布力與集中力共同作用的平面桁架結(jié)構(gòu)。彈性模量E=207 GPa,密度ρ=7850 kg/m3,豎直桿橫截面積為0.03 m2,水平桿橫截面積為0.0068 m2。半正弦均勻分布力F1(t)作用在水平桿上,半正弦集中力F2(t)作用在節(jié)點2和節(jié)點3處。

圖8　均布力與集中力作用的平面桁架

圖9所示為應(yīng)用本文方法得到的桁架全局動態(tài)應(yīng)力及關(guān)鍵時間點，從圖中可以看出全局動態(tài)應(yīng)力不僅呈非線性多峰值狀態(tài)，而且似乎沒有規(guī)律可循；可以肯定的是無論全局動態(tài)應(yīng)力解空間多么復(fù)雜，譜單元Lagrange插值總是體現(xiàn)其優(yōu)越性，在此基礎(chǔ)上找到關(guān)鍵時間點。

圖9　均布力與集中力作用的平面桁架全局動態(tài)應(yīng)力

圖10　均布力與集中力作用的平面桁架單元10的動態(tài)應(yīng)力

圖11　均布力與集中力作用的平面桁架節(jié)點12的位移響應(yīng)

圖10、圖11所示分別是從局部動態(tài)應(yīng)力與位移的角度來識別關(guān)鍵時間點的結(jié)果。將不同識別方法得到的關(guān)鍵時間點匯總至表2，可以發(fā)現(xiàn)3種方法總是有時間差的。另外基于全局動態(tài)應(yīng)力的關(guān)鍵時間點識別方法還可獲得結(jié)構(gòu)最大應(yīng)力值,這為結(jié)構(gòu)動態(tài)設(shè)計從另一個角度提供參考。當然，如果以位移或局部應(yīng)力為依據(jù)識別關(guān)鍵時間點，其對應(yīng)時刻不是結(jié)構(gòu)最大應(yīng)力的時刻，這樣可能造成結(jié)構(gòu)欠設(shè)計。

表2　均布力與集中力作用的

5結(jié)論

(1)結(jié)構(gòu)動態(tài)應(yīng)力最大的時刻與動態(tài)位移最大的時刻存在一個時間差，以動態(tài)應(yīng)力最大時刻為結(jié)構(gòu)關(guān)鍵時間點，雖然相比以位移最大時刻為結(jié)構(gòu)關(guān)鍵時間點合理,然而結(jié)構(gòu)某一點應(yīng)力達到最大時，全局動態(tài)應(yīng)力并不一定達到最大，因此以全局動態(tài)應(yīng)力最大時刻為結(jié)構(gòu)關(guān)鍵時間點最合理。

(2)應(yīng)用譜單元離散插值全局動態(tài)應(yīng)力解空間,可獲得高精度的插值函數(shù)。結(jié)構(gòu)全局動態(tài)應(yīng)力呈現(xiàn)極度非線性多峰值，采用在GLL點處Lagrange插值獲得高精度的全局動態(tài)響應(yīng)解空間近似函數(shù)。

(3)通過全局優(yōu)化算法對譜單元離散插值函數(shù)進行尋優(yōu)，獲得滿意的關(guān)鍵時間點。應(yīng)用局部優(yōu)化算法對譜單元離散插值函數(shù)尋找最大值點，往往會由于初始點的不同而找到不同的局部最優(yōu)點，區(qū)域細分全局優(yōu)化算法可以找到全局最大值與全局最小值，對應(yīng)的就是全局應(yīng)力最大拉應(yīng)力與全局應(yīng)力最大壓應(yīng)力。

參考文獻:

[1]James M L.Vibration of Mechanical and Structural Systems[M].2nd ed. New York:Harper Collins,1994.

[2]Clough R W.Dynamics of Structures[M].New York:McGraw-Hill,1993.

[3]Grandhi R V, Haftka R T, Watson L T. Design-oriented Identification of Critical Times in Transient Response [J]. AIAA Journal, 1986, 24(4):649-656.

[4]Choi W S, Park G J. Transformation of Dynamic Loads into Equivalent Static Loads Based on Modal Analysis [J]. Int. J. Numer. Meth. Engng.,1999,46: 29-43.

[5]張艷崗.基于關(guān)鍵時間點的能量等效靜態(tài)載荷法及結(jié)構(gòu)動態(tài)響應(yīng)優(yōu)化研究[D].太原：中北大學(xué),2014.

[6]Zhang Yangang, Su Tiexiong. Dynamic Analysis of Diesel Engine Piston Based on Time Spectral Element Method [J]. Applied Mechanics and Materials，2013，415：565-568.

[7]毛虎平，吳義忠，陳立平. 基于時間譜元法的動態(tài)響應(yīng)優(yōu)化[J]. 機械工程學(xué)報，2010，46(16)：79-87.

Mao Huping, Wu Yizhong, Chen Liping. Dynamic Response Optimization Based on Time Spectral Element Method[J].Journal of Mechanical Engineering, 2010，46(16)：79-87.

[8]張艷崗,蘇鐵熊，毛虎平，等.動態(tài)應(yīng)力解空間譜元離散的關(guān)鍵時間點識別方法[J].機械工程學(xué)報,2014,50(3):82-87.

Zhang Yangang,Su Tiexiong,Mao Huping,et al.Critical Time Points Identification Method for Solution Space of Dynamic Stress Based on Spectral Element [J].Journal of Mechanical Engineering,2014,50(3):82-87.

[9]Pozrikidis C.Introduction to Finite and Spectral Element Methods Using Matlab[M].New York:Chapman and Hall/CRC,2005．

[10]張汝清,殷學(xué)綱,董明.計算結(jié)構(gòu)動力學(xué) [M].重慶: 重慶大學(xué)出版社, 1987.

[11]毛虎平, 吳義忠, 陳立平. 基于多領(lǐng)域仿真的SQP并行優(yōu)化算法[J]. 中國機械工程, 2009, 20(15): 1823-1829.

Mao Huping,Wu Yizhong,Chen Liping.SQP Parallel Optimization Algorithm Based on Multi-domain Simulation[J].China Mechanical Engineering,2009,20(15):1823-1829.

(編輯袁興玲)

Critical Time Point Identification for Global Dynamic Stress Based on Spectral Element Interpolation of Solution Space

Zhang YangangMao HupingSu TiexiongLi KunWang Jun

North University of China,Taiyuan,030051

Abstract:According to the complexity and high-cost characteristics of dynamic analyses in the structural dynamic optimization processes,a critical time points identification method for solution space of global dynamic stress based spectral element interpolation was proposed.The most dangerous time points of the structure dynamic response were found.The modal stress distribution was obtained by using the modal superposition method,and then, the solution space domain of the global dynamic stress was computed.The Lagrange interpolation techniques were applied in the solution space domain to get the high precision function of the global dynamic stress. The absolute maximum points, that is, the critical time points,were found by executing the global DIRECT optimizer. Last, two examples of 124-member plane truss and the structure subjected to uniform and concentrated force were used to illustrate the feasibility and validity of the proposed method.

Key words:spectral element discretion of solution space;global dynamic stress;Lagrange interpolation;critical time point

作者簡介：張艷崗,男,1981年生。中北大學(xué)機械與動力工程學(xué)院講師、博士。研究方向為結(jié)構(gòu)動力學(xué)建模與仿真、結(jié)構(gòu)動態(tài)響應(yīng)優(yōu)化。毛虎平,男,1974年生。中北大學(xué)機械與動力工程學(xué)院副教授、博士。蘇鐵熊,男,1963年生。中北大學(xué)機電工程學(xué)院教授、博士研究生導(dǎo)師。李坤,男,1981年生。中北大學(xué)機械與動力工程學(xué)院講師、博士。王軍，男, 1979年生。中北大學(xué)機械與動力工程學(xué)院博士生。

中圖分類號：TH122

DOI：10.3969/j.issn.1004-132X.2016.05.021

基金項目：國家自然科學(xué)基金資助項目(51275489)

收稿日期：2014-12-10