時間譜元法在動態響應優化中的應用

毛虎平，吳義忠，李建軍，王應紅

（1．中北大學機電工程學院，山西 太原 430074；2．華中科技大學國家CAD支撐軟件工程技術研究中心，湖北 武漢4 30074；3．中北大學理學院，山西 太原 030051；4．上海汽車集團股份有限公司乘用車公司，上海 201804）

引 言

機械結構幾乎都在動載荷作用下運行，機器的各種動態性能是與時間相關的函數。欲使機器的性能最優，動態優化是很有必要的。故而必須處理由于動態載荷作用機器的各種動態響應，包括目標函數與約束函數。在時間域上，較精確地求出相關響應，并且要求滿足時間函數約束，而且不能漏掉任何響應，特別是響應極值；而應用很小步長計算動態響應，一定程度上可以減小時間函數約束的可能。在優化過程中，每一次迭代都要重新計算時間函數目標函數、時間函數約束，因此，要求求解動態響應必須高精度而且高效。

有限差分法是比較流行的求解動態響應微分方程的方法。應用很小時間步長計算二階微分方程，自適應調節時間步長，直到獲得要求的精度。時間步長直接影響計算的穩定性和準確性，這極大地影響了計算效率。特別是當結構受脈沖載荷時，差分法很難獲得滿意的結果。文獻［1］應用譜元法的h型精細方案對受脈沖的動態系統分析求解，在沒有增加單元數量的基礎上獲得了譜收斂精度。

文獻［2］中，提出將設計變量和位移響應、設計變量和位移響應與速度響應、設計變量和位移響應與速度響應加速度響應分別作為優化變量3種方案，應用有限差分法近似和時間相關的約束，運動微分方程被作為等式約束處理。對于動態響應最大值最小問題，文獻［3］提出了直接處理目標函數的方法，解決了迭代過程中最大響應震蕩造成收斂困難的問題。文獻［4］應用時間有限元法將動態響應微分方程轉化為代數方程，然后進行了二階靈敏度分析，為動態響應優化奠定基礎。文獻［5］中，應用了時間譜元法離散時間微分方程，GLL點法和關鍵點法處理與時間相關的約束，然而GLL點法不穩定，對單元數和插值次數比較敏感，關鍵點法由于每一次迭代的最大值的位置發生變化，因此收斂速度非常慢甚至是發散。

20世紀70年開始研究承受瞬態載荷的結構優化。2006年，Kang等首次把動態響應優化作為優化的一個分支［6］。在其優化中，設計變量每迭代一次，時間響應需要計算一次，同時約束必須在整個時間區間內被滿足。處理時間約束有幾種方法。一種是只在響應的全局最大值處滿足約束；另一種更加穩定的方法是使在更小時間步長上滿足約束，保證中間點處的約束失效不可能發生。在準靜態方法中，應用了這種方法使多個等靜態載荷滿足約束而不是動態載荷［7］。在這些方法中，約束數大大增加。因為在優化迭代過程中，要計算這些約束的靈敏度，這樣優化耗費也在增加。再一種更加有效的方法是只在響應的局部極值點處處理約束。這樣可以減少約束的數量和靈敏度的次數［8］。

本文提出關鍵點及其相鄰的GLL點法處理與時間相關的約束，在最大值附近的一些點能夠被包括在約束集，克服當前迭代的局部最大值可能在下一次迭代不是局部最大值的問題。

1 時間譜元法

1984年Patera提出了譜元法［9］，在流體動力學中得到應用，其融合了有限元法的處理邊界、結構的靈活性與譜方法的快算收斂性。在要求相同精度的情況下，譜元法能夠采用較少的單元減小了計算開銷。在一個單元內，將時間離散為與GLL多項式零點相對應的網格點，在這些點上進行Lagrange插值。

從理論上分析，在一定點數上插值，當這些點是對應的正交多項式的零點時獲得的插值精度最高［10］。如圖1所示，均勻點分布和GLL點分布分別近似Runge函數式。

圖1 GLL點和平均點Lagrange插值Runge函數Fig．1 Lagrange interpolation in GLL points and average points for Runge function

在有限元法中，隨著插值次數增加，近似的數值誤差極大或近似失敗；而GLL插值有明顯的優勢。GLL點在標準區間分布如圖2所示。

時間譜元法的主要步驟：（Ⅰ）首先轉化二階微分方程組為一階微分方程組，應用Bubnov－Galerkin原理轉化為等效的積分形式，代入單元積分表達式，得到單元譜元方程。（Ⅱ）將時間段劃分成若干單元，每個單元由若干個節點組成。（Ⅲ）用Legendre

圖2 GLL點分布（p表示單元內插值節點數，在［－1，1］區間上）Fig．2 GLL point distribution（The prepresents the number of interpolation nodes．Intervels are［－1，－1］）

正交多項式作為基函數，將每一個單元的近似解表達為基函數的線性組合。（Ⅳ）將所有單元的譜元方程按一定規則合成其總體譜元方程。（Ⅴ）解時間的總體譜元方程，得到其全局近似解。

1．1 動態響應微分方程及其形式轉化

動態響應微分方程可表達為

式中M為質量矩陣，C為阻尼矩陣，K為剛度矩陣，F為時間函數載荷向量，x為位移向量為加速度向量為速度向量。初始條件v0；其中，M，C，K不隨時間變化；為時間的函數，F為時間的任意函數，t∈［t0，tn］。

1．2 單元劃分

圖3 時間譜元法離散時間Fig．3 Divided the time domain temporal SEM

如圖3所示，將仿真時間t∈［t0，tn］分割為n個互不相交的單元，即［t0，t1］，［t1，t2］，［t2，t3］，…，［tn－1，tn］，對每單元配置若干個點。在配置點時，可以根據動態載荷的變化情況，在載荷突變的地方配置多一些點，在載荷平緩的地方可以配置少一些點。這樣能保證在動態載荷的所有地方都能有足夠的精度。

1．3 單元分析

區間分段以后，將每一個時間段劃分單元，為提高近似解的精度，可以在單元區間上再增加插值點來增加近似解的階次，這在有限元里被稱為高階元，即hp方法。

在時間譜元法中，為了達到譜方法的收斂速度，并且增加數值計算的精度，單元內插值點應放置在特殊的配置點上，它們是由N＋1個GLL點構成的，而基函數由有限項正交多項式的和來表示，構成單元上任意配置點的形函數，如此可在有限的插值點上獲得指數級的收斂速度。

在每一個單元上，將定義在插值點上的形函數通過Lagrange插值得到近似函數。換句話說，任意一個定義在參考單元上的函數?x（j）（ξ），都可以用下式近似

式中（ξ）為j單元的k次Lagrange多項式，ξk為定義在［－1，1］上的 GLL點，x（j）（ξk）為單元j上未知節點在GLL點的值。

有兩種類型正交多項式可用于時間譜元法中：分別為Chebyshev和Legendre正交多項式。這里只討論Legendre展開，勒讓德多項式是Sturm－Liouville Form方程的解它的權函數ω＝1。k階Legendre多項式可以定義如下

由于Legendre多項式的零點不包括區間端點，因此引入Lobbato多項式。Lobbato多項式是通過Legendre多項式的微分定義的［10］

將式（3）展開，其中一個方程可以表達為

式中As為耦合矩陣。

將每一個狀態變量在給定時間段內離散化，近似為m次Lagrange多項式（4），在每一個單元上應用Bubnov－Galerkin法使得插值誤差最小［11］，得出

將每一個單元用矩陣形式表示為

1．4 集成總體時間譜元方程并求解

總體合成的任務是將所有單元譜元方程按序疊加，得到總體譜元方程，即線性方程組。對于個狀態變量，通過耦合矩陣的張量叉乘得到全部狀態變量的全局組裝式

化簡式（11）得到

2 動態響應優化

機械結構動態響應優化設計問題的目標函數為選擇設計變量向量X（x1，x2，…，xk）T，以使系統在受到瞬變載荷作用下，在給定時間區間［0，T］內滿足瞬時動態性能（振幅或相對位移極限、動應力、動應變破壞極限等）及設計變量允許變化范圍等約束條件，并使系統某些關心位置或坐標的最大動態響應（位移、速度、加速度）在某種意義或準則下達到最優。其數學模型可表示為

式中X為系統的設計變量向量，由系統的幾何參數和物理參數組成；z（t）為系統的狀態變量向量，由系統運動狀態的廣義坐標如位移、速度、加速度等組成，其要滿足機械系統運動規律的運動方程，即動態特性的數學描述（t）＝p［X，t，z（t），F（t）］；J為目標函數；hi，gj分別表示等式約束與不等式約束，它們在所有時間點（?t∈［0，T］）上都要滿足。

處理與時間相關的約束有兩種方法（見文獻［12］）：一種是約束在所有GLL點得到滿足；另一種是約束在每一個單元的絕對值極值點得到滿足。第一種方法要求GLL點之間的距離盡量的小，這樣可以減少漏掉GLL之間的絕對值極值點得不到滿足的可能性，因此約束數量比較龐大，同時譜元法的求解開銷也較大；第二種方法可以用滿足精度要求的較少譜單元求解運動微分方程，在每一個單元內，對其高次Lagrange函數進行一維搜索找到單元絕對值極值點，當然單元絕對值極值點是震蕩的，因此每迭代一步，單元絕對值極值點要重新計算。文獻［5］中說明了GLL點法處理與時間相關的約束，非常不穩定，但單元數和插值點數達到一定大時，才能找到最優值，否則找不到最優值；而關鍵點法，當前迭代的局部最大值可能在下一次迭代不是局部最大值，這樣收斂速度變慢甚至不收斂。本文提出將所有單元響應絕對最大值以及其相鄰的2個GLL點作為約束，在最大值附近的一些點能夠被包括在約束集，如圖4所示。劃分 3 個單 元，即［t0，t1］，［t1，t2］，［t2，t3］，關鍵點分別為p1，p2，p3，其相鄰的 GLL點分別為p11，p12，p21，p22，p31，p32。

圖4 單元絕對極值點及其相鄰的兩個GLL點作為約束Fig．4 The constraints are the absolute extreme value and the two adjacent GLL points

這樣，關鍵點集PC可以采用下式表示

式中tmin／max＝argminL（t），tGLL0為tmin／max的左端附近GLL點，tGLL1為tmin／max的右端附近 GLL點。從式（14）中可以很清楚地看出本文與作者前期研究的不同，關鍵點法僅僅考慮tmin／max點，而沒有考慮左右端的點。

3 算例分析

在各種各樣的機械系統（如，減振器、武器后座機構、飛機起落架等）中，系統的主質量在該系統達到穩態之前可能作大幅度振動，當激勵力的頻率接近于主系統的固有頻率時，系統可能發生破壞，對激勵的瞬時動態響應必須約束。對于這種問題，可以用抽象模型的優化來解決。

3．1 線性兩自由度減震器最優設計

已知兩自由度系統，參數如圖5所示，m1＝4．534kg，m2＝0．453 4kg，k1＝1 749．03N／m，激勵頻率是主質量固有頻率Ωn的1．2倍，即ω＝23．57。求阻尼器的剛度系數k2和阻尼系數c使得主質量的最大位移響應最小，同時滿足振動空間約束條件減振器位移響應和主質量位移響應之差不能超過3倍的主質量最大響應，以及穩態約束條件。初始條件均為零。

圖5 減振器Fig．5 Shock absorber

引入符號：

那么，優化模型表示為

根據文獻［12］，引入人工設計變量x5，可以把式（18）轉化為

式中a＝0．030 48m；x2st，x4st分別表示主質量和阻尼器的穩態位移響應，表達式為

表1和2為與文獻［13］優化結果比較，可以看出本文提出的關鍵點及其相鄰的GLL點法處理與時間相關的約束非常穩定，對不同的單元和插值數都找到了最優值，而且迭代次數比較少，最大的13次；而在作者的前期研究（文獻［5］）中關鍵點法處理與時間相關的約束不穩定，比如當Nel＝4，m＝5以及當Nel＝10，m＝6時，沒有找到最優值，并且迭代次數是關鍵點及其相鄰的GLL點法的兩倍多。兩自由度問題是最簡單的多自由度問題，這樣就能說明本文研究結果更具有實際意義。

其實，從耗時來分析，還可以給SQP提供解析梯度函數或者簡單的數值梯度函數，這樣函數的評估次數能明顯減少［14］；SQP算法采用差分法計算梯度，每次迭代需要多次進行仿真函數評估（評估次數等于設計變量維數的兩倍乘以實例數＋1），優化求解很費時［15］；對于基于譜元法的動態響應優化，每評估一次函數值就要運用譜元法更新一次，很明顯，單元數越多及其插值次數越高，耗時越多，減少函數評估次數，自然耗時就能下降。從文獻［12］可知，作者研究的單自由度的緩沖器優化結果波動的原因是是否有網格點在響應的絕對值最大點處，文獻［12］第10頁圖4說明了此問題。因此這里只比較了關鍵點法與本文提出的方法。表3是參考文獻中的數據和本文的數據，說明本文方法的正確性和可行性。

表1 關鍵點及其相鄰的GLL點法數值實驗結果Tab．1 Numerical experimental results of the critical points and adjacent GLL points

表2 關鍵點法數值實驗結果Tab．2 Numerical experimental results of the critical points

圖6 主質量動態響應Fig．6 The dynamic response of main mass

表3 文獻數據與本文方法優化結果數據的比較Tab．3 Comparing of the literature and this paper

3．2 汽車懸掛系統動態優化設計

汽車系統可近似簡化為圖7所示的5自由度的動力學模型。圖中m1為駕駛員及其座位的質量，它由彈簧k1和阻尼器c1支撐在車體上。汽車車體的質量、前后車軸和車輪的質量分別為m2，m4和m5。汽車車體由連于前后車軸的彈簧k2和k3及阻尼器c2和c3所支撐。k4，k5和c4，c5表示輪胎的剛度和阻尼系數。車體對其質量中心的轉動慣量用I表示。L表示前后輪距長度。f1（t）表示由于汽車行駛道路表面起伏不平所激起的前后輪的位移函數。yi（t）為5個自由度的廣義坐標。

圖7 5自由度汽車簡化模型Fig．7 Five DOFs automobile simplified model

（1）選擇設計變量

本優化問題僅取模型中與懸架系統有關的彈簧系數ki和阻尼系統ci（i＝1，2，3）作為設計變量。即

（2）建立目標函數

以司機和乘客座位的瞬變動態最小對其機械系統的結構進行動態響應優化設計，即在汽車懸掛系統的結構設計時，要求汽車在各種車速和路面條件下，使得司機和乘客座位的最大加速度響應最小。可表示為

（3）建立動態響應優化數學模型

式中 動態響應微分方程組為等式約束，θ1為車體與司機座位間的相對位移約束，θ2為車體與前輪之間的相對位移約束，θ3為車體與后輪之間的相對位移約束，θ4為路面與前輪之間的相對位移約束，θ5為路面與后輪之間的相對位移約束。

本優化問題僅取模型中與懸架系統有關的彈簧系數ki和阻尼系統ci（i＝1，2，3）作為設計變量，以司機和乘客座位的瞬變動態響應最小為目標，以系統狀態方程和振幅為約束。在如圖8所示路面條件下行駛。

圖8 路面輪廓Fig．8 Pavement surface outline

汽車速度v＝11．43m／s。路面激勵頻率為ωi＝πv／Li，L1＝9．144m，L2＝3．657 6m，后輪到達前輪未知的時間間隔為tσ＝0．266 7s，其中路面起伏不平所激起前后輪的垂直位移函數f1（t），其值為

由算例1可以看出本文提出的方法處理與時間有關的約束的多自由度動態問題的優越性，因此算例2僅僅采用本文提出的方法對一組數據進行計算，結果如表4所示。單元數Nel＝20，插值點數p＝6時，就得到比文獻更好的最優值，迭代僅僅51次；而關鍵點法不僅結果沒有本文提出的方法好而且迭代次數為72次，如圖9所示。

表4 本文方法與關鍵點法數值實驗結果Tab．4 Numerical experimental results of the method for this paper and the critical points method

圖9 座位的動態響應迭代過程（Nel＝20，p＝6）Fig．9 The iterative process of seat＇s dynamic response（Nel＝20，p＝6）

4 結 論

機械系統動態響應優化設計有很好的應用前景［11］。動態響應必須滿足依賴于時間的微分方程。為了獲得最優解且滿足與時間有關的約束，要求獲得整個時間上的響應。本文應用譜元法計算整個時間內的響應，將微分方程組轉化為代數方程組；應用序列二次規劃（SQP）優化算法嵌套黃金分割法進行優化計算。

提出了關鍵點及其相鄰的GLL點法處理與時間有關的約束的處理，與關鍵點方法相比不僅提高了優化的穩定性，即優化過程對時間譜單元數和單元插值次數不敏感，無論劃分多少個單元，多少次插值，都能獲得最優解，而且迭代次數提高近一倍；這是由于動態響應優化中，時間函數目標與時間函數約束的極值點隨著參數的變化是震蕩的，而關鍵點法中只考慮關鍵點（極值點），這樣造成優化不穩定，有時得不到最優解，有時可以獲得，即使獲得最優解，其迭代次數也很大，效率較低；本文提出的方法解決了此問題。在本文方法的時間譜元法中，提出單元離散的疏密大小與插值次數看作動態載荷變化程度的函數，即在動態載荷變化急劇的區域，單元離散較密且插值次數高；而在動態載荷變化平緩的區域，單元離散較稀且插值次數低。通過分析2自由度減振器設計問題和5自由度汽車懸掛系統設計問題，驗證了本文提出方法的可行性和優越性。

二自由度減振器設計問題代表最簡單的多自由度動力學設計問題，對于機械零部件的動態優化設計問題以及動載荷彈性分布參數系統的優化設計，可以參考本文的方法解決，比如在固定端承受振動輸入的矩形變截面梁的動態優化設計，在不同邊界條件下承受垂直平面均布瞬態動載荷作用的彈性梁的動態優化設計等。

［1］ Kurdi M H，Beran P S．Spectral element method in time for rapidly actuated systems［J］．Journal of Computational Physics，2007，227（3）：1 809—1 835．

［2］ Wang Q，Jasbir S A．Alternative formulations for transient dynamic response optimization［A］．46th AIAA／ASME／ASCE／AHS／ASC Structures，Structural Dynamics ＆ Materials Conference［C］．Austin，Texas，2005．

［3］ Choi D H，Park H S，Kim M S．A direct treatment of min－max dynamic response optimization problems［J］．AIAA－93－1352－CP．

［4］ Park S，Kapania R K，Kim S J．Nonlinear transient response and second－order sensitivity using time finite element method［J］．AIAA J，1999，37（5）：613—622．

［5］ 毛虎平，吳義忠，陳立平．基于時間譜元法的動態響應優化［J］．機械工程學報，2010，46（16）：79—87．

Mao Huping，Wu Yizhong，Chen Liping．Optimization of dynamic respond based on temporal spectral element method［J］．Chinese Journal of Mechanical Engineering，2010，46（16）：79—87．

［6］ Kang B S，Park G J，Arora J S．A review of optimization of structures subjected to transient loads［J］．Struct．Multidisc．Optim．，2006，31（2）：81—95．

［7］ Choi W S，Park G J．Structural optimization using equivalent static loads at all time intervals［J］．Compute Methods Appl．Mech．Eng．，2002，191：2 077—2 094．

［8］ Grandhi R V，Haftka R T，Watson L T．Design oriented identification of critical times in transient response［J］．AIAA Pap．，1984：649—656．

［9］ Patera A T．A spectral element method for fluid dynamics laminar flow in a channel expansion［J］．Journal of Computational Physics，1984，54：468—488．

［10］Pozrikidis C．Introduction to Finite and Spectral Element Methods Using Matlab［M］．Chapman and Hall／CRC，2005．

［11］Seymour V P．On the Legendre－Gauss－Lobatto points and weights［J］．Journal of Scientific Computing，1999，14（4）：347—355．

［12］Kurdi M H，Beran P S．Optimization of dynamic response using temporal spectral element method［A］．46th AIAA Aerospace Sciences Meeting and Exhibit［C］．Reno，Nevada，2008．

［13］E J豪格，J S阿羅拉，著．實用最優設計—機械系統與結構系統［M］．郁永熙，丁慧梁，譯．北京：科學出版社，1985．

Haug E J，Arora J S．The Practical and Optimal Design［M］．Beijing：Science Publishing，1985．

［14］薛定宇，陳陽泉，著．高等應用數學問題的 MATLAB求解［M］．北京：清華大學出版社，2004．

Xue Dingyu，Chen Yangquan． Advanced Applied Mathematical Problem Solutions with MATLAB［M］．Beijing：Tsinghua University Press，2004．

［15］毛虎平，吳義忠，陳立平．基于多領域仿真的SQP并行優化算法［J］．中國機械工程，2009，20（15）：1 823—1 829．

Mao Huping，Wu Yizhong，Chen Liping．SQP parallel optimization algorithm based on multi－domain simulation［J］．China Mechanic Engineering，2009，20（15）：1 823—1 829．

［16］Hsieh C C，Arora J S．Design sensitivity analysis and optimization of dynamic response［J］．Computer Methods in Applied Mechanics and Engineering，1984，43：195—219．

［17］Hsieh C C，Arora J S．A hybrid formulation for treatment of point－wise state variable constraints in dynamic response optimization［J］．Computer Methods in Applied Mechanics and Engineering，1985，48：171—189．

［18］Paeng J K，Arora J S．Dynamic response optimization of mechanical system with multiplier methods［J］．ASME Journal of Mechanism，Transmission，and Automation in Design，1989，111：73—80．

［19］Choi D H，Park H S，Kim M S．A direct treatment of min－max dynamic response optimization problems［J］．AIAA－93－1352－CP．