基于序列響應(yīng)面方法的高強(qiáng)度鋼材料參數(shù)反求

張 勇 陸 勇

華僑大學(xué)，廈門，361021

0 引言

汽車碰撞安全性是汽車設(shè)計(jì)過程中主要考慮的因素之一［1］，隨著計(jì)算機(jī)硬件技術(shù)及數(shù)值仿真技術(shù)的快速發(fā)展，基于數(shù)值仿真技術(shù)的汽車安全性設(shè)計(jì)正逐漸替代傳統(tǒng)的物理整車碰撞安全性實(shí)驗(yàn)。但是，整車碰撞安全性的仿真精度在很大程度上取決于所選材料的本構(gòu)方程，而本構(gòu)方程中的材料參數(shù)，如彈性模量、硬化指數(shù)、應(yīng)變率、應(yīng)變速度導(dǎo)致的溫度效應(yīng)等將直接影響到仿真結(jié)果的精確性［2］，因此，準(zhǔn)確獲得材料參數(shù)對(duì)提高碰撞安全仿真結(jié)果的可靠性有著顯著的科學(xué)意義和實(shí)際應(yīng)用價(jià)值。

此外，準(zhǔn)靜態(tài)工況下進(jìn)行的單向拉伸試驗(yàn)得到的材料參數(shù)，往往難于反映材料在大變形、大應(yīng)變、材料快速變形過程中的發(fā)熱等對(duì)材料性能的影響。材料參數(shù)的諸多不確定性使得碰撞安全性仿真的結(jié)果時(shí)常與實(shí)驗(yàn)結(jié)果具有較大的誤差。因此，本文采用有限元技術(shù)與數(shù)值優(yōu)化相結(jié)合的方法來反求材料的各類參數(shù)。該方法采用優(yōu)化方法進(jìn)行迭代尋優(yōu)，從而調(diào)整實(shí)驗(yàn)曲線點(diǎn)與仿真曲線點(diǎn)的逼近程度。但是對(duì)于碰撞安全這類涉及大位移、大變形的高度非線性仿真問題，當(dāng)數(shù)值優(yōu)化方法需要成千次調(diào)用正問題的計(jì)算時(shí)，計(jì)算成本將成幾何級(jí)數(shù)的增長(zhǎng)，導(dǎo)致優(yōu)化設(shè)計(jì)不可行，此外，基于整個(gè)設(shè)計(jì)空間構(gòu)造的近似模型的精度常常較低，導(dǎo)致基于該模型的優(yōu)化解的精度差，且不易找到全局最優(yōu)點(diǎn)［3］。因此，筆者采用逐次逼近的序列響應(yīng)面方法，使其能在較小的設(shè)計(jì)空間中構(gòu)建出高精度的近似模型，從而較好地解決上述問題。

基于物理實(shí)驗(yàn)與數(shù)值仿真的幾何與材料參數(shù)的識(shí)別研究，已取得一些成果［4－13］，但是，基于近似模型方法的材料參數(shù)反求，特別對(duì)于高強(qiáng)度鋼板DP800的材料參數(shù)反求的文獻(xiàn)仍然鮮有報(bào)道。因此，筆者將序列響應(yīng)面近似模型方法與小種群遺傳算法相結(jié)合，在基于DP800制作的薄壁吸能管碰撞試驗(yàn)數(shù)據(jù)的基礎(chǔ)上，基于真實(shí)的試驗(yàn)條件建立了相關(guān)的有限元仿真模型，從而來開展針對(duì)高強(qiáng)度鋼板DP800的各項(xiàng)材料性能參數(shù)的反求。

1 序列響應(yīng)面方法

1.1 響應(yīng)面模型的構(gòu)建

響應(yīng)面方法是一種將試驗(yàn)設(shè)計(jì)與數(shù)理統(tǒng)計(jì)相結(jié)合來建立經(jīng)驗(yàn)?zāi)Ｐ偷囊环N優(yōu)化方法［3－4］，它采用擬合與插值的方法分析試驗(yàn)設(shè)計(jì)點(diǎn)，從而構(gòu)建出目標(biāo)函數(shù)與約束函數(shù)的響應(yīng)面近似模型，并以該模型來替代物理有限元模型進(jìn)行尋優(yōu)并預(yù)測(cè)非試驗(yàn)設(shè)計(jì)點(diǎn)的響應(yīng)值。

實(shí)際應(yīng)用過程中，一般采用二階多項(xiàng)式響應(yīng)面模型，其函數(shù)關(guān)系為

式中，j為設(shè)計(jì)變量的個(gè)數(shù)；β0、βi、βii、βil為待定系數(shù)。

當(dāng)選取了n個(gè)試驗(yàn)設(shè)計(jì)點(diǎn)時(shí)，響應(yīng)面函數(shù)可以通過矩陣形式表示為

式中，ε為誤差矢量。

回歸系數(shù)β采用最小二乘擬合方法確定，即

1.2 響應(yīng)面模型的更新與逐次逼近

設(shè)計(jì)空間的大小對(duì)響應(yīng)面近似模型的精度具有重要的影響，本文采用在優(yōu)化過程的不同階段逐步更新與縮放設(shè)計(jì)空間的方法，使優(yōu)化設(shè)計(jì)問題的求解域從初始設(shè)計(jì)空間中一個(gè)較大的設(shè)計(jì)域轉(zhuǎn)化成一個(gè)較小的信賴域。與在整個(gè)設(shè)計(jì)空間構(gòu)建的近似模型相比，較小的信賴域內(nèi)構(gòu)建的近似模型的精度更高。

設(shè)經(jīng)過第k次優(yōu)化設(shè)計(jì)后，設(shè)計(jì)空間將變?yōu)榈趉＋1次信賴域，而第k＋1次信賴域的中心即為第k次信賴域中的優(yōu)化設(shè)計(jì)點(diǎn)x［k］，此次信賴域的縮放和移動(dòng)關(guān)系可用下式表示［3］：

為了保證信賴域更新過程中各變量的縮放比例保持一致，取λ作為所有變量的部分因子，也稱之為變換因子，其表達(dá)式為

此時(shí)，第i個(gè)設(shè)計(jì)變量在第k＋1次設(shè)計(jì)子空間的上下界為

在優(yōu)化設(shè)計(jì)過程中，整個(gè)信賴域的移動(dòng)和縮放如圖1所示。

圖1 信賴域的變化示意圖

通過設(shè)計(jì)空間的逐步移動(dòng)和縮放，優(yōu)化問題將收斂在一個(gè)較小的信賴域中進(jìn)行，得到精度較高的響應(yīng)面近似模型與最優(yōu)設(shè)計(jì)點(diǎn)，該方法能夠有效地消除噪聲，改善優(yōu)化算法的收斂性能。

1.3 D最優(yōu)試驗(yàn)設(shè)計(jì)

在整個(gè)設(shè)計(jì)空間選取有限數(shù)量的樣本點(diǎn)，使之盡可能地反映設(shè)計(jì)空間的特性，即稱為試驗(yàn)設(shè)計(jì)。試驗(yàn)設(shè)計(jì)是構(gòu)建近似模型的基礎(chǔ)，設(shè)計(jì)樣本點(diǎn)的選取是否合適對(duì)后續(xù)的響應(yīng)面近似模型的構(gòu)建起著非常重要的作用，它將直接影響構(gòu)建的近似模型的精度，故試驗(yàn)設(shè)計(jì)方法的選取也顯得尤為重要。因此，如何在設(shè)計(jì)空間內(nèi)最優(yōu)布置試驗(yàn)設(shè)計(jì)點(diǎn)，也是當(dāng)前試驗(yàn)設(shè)計(jì)理論研究的熱點(diǎn)之一。本文采用了Draper等［14］提出的D最優(yōu)試驗(yàn)設(shè)計(jì)方法來獲取采樣點(diǎn)。

2 DP800高強(qiáng)度鋼材料參數(shù)的反求

2.1 DP800薄壁吸能梁的試驗(yàn)與反求仿真模型

薄壁吸能梁的碰撞試驗(yàn)示意圖以及它的幾何截面參數(shù)分別如圖2與圖3所示。其中，試件夾持在夾具中，其夾持長(zhǎng)度L2＝100mm，碰撞吸能自由端的長(zhǎng)度L1＝310mm，碰撞質(zhì)量臺(tái)車的質(zhì)量為600kg，它以15m／s的速度正面碰撞DP800吸能梁，碰撞過程中通過力與加速度傳感器采集碰撞力與減速度［4］。

圖3中，薄壁梁的截面寬度b1＝60mm，截面高度b2＝60mm，界面的圓角半徑R1＝3mm，R2＝2mm，其帽形梁的焊接法蘭長(zhǎng)度為25mm，高強(qiáng)度鋼板DP800的板料厚度為1.2mm。碰撞試驗(yàn)后得到的碰撞力及其變形量等相關(guān)試驗(yàn)數(shù)據(jù)源于參考文獻(xiàn)［4］。

圖2 薄壁梁碰撞吸能示意圖

圖3 薄壁梁的幾何截面圖

為了模擬整個(gè)碰撞過程并反求DP800的材料參數(shù)，筆者按照試驗(yàn)條件建立了其有限元模型，如圖4所示，其碰撞仿真過程在非線性有限元軟件LS－DYNA中進(jìn)行。薄壁梁在高速碰撞過程中的快速塑性變形使得高強(qiáng)度鋼DP800的溫度迅速上升，導(dǎo)致材料產(chǎn)生溫度軟化效應(yīng)，同時(shí)，高速碰撞過程中，高強(qiáng)度鋼DP800的應(yīng)變率變化范圍也對(duì)仿真精度具有較大的影響。因此，本文采用綜合考慮了溫度軟化效應(yīng)、應(yīng)變強(qiáng)化效應(yīng)、應(yīng)變率強(qiáng)化效應(yīng)等影響的Johnson–Cook（J－C）本構(gòu)材料模型，并以它作為碰撞仿真模型中高強(qiáng)度鋼DP800的本構(gòu)模型。

圖4 薄壁梁的有限元仿真圖

J－C模型是一種動(dòng)態(tài)的本構(gòu)模型，J—C模型的等效流動(dòng)應(yīng)力隨溫度、應(yīng)變、應(yīng)變率的變化而變化，其具體本構(gòu)表達(dá)式為

2.2 DP800材料參數(shù)的反求方法

本文采用小種群遺傳算法對(duì)高強(qiáng)度DP800的材料參數(shù)進(jìn)行反求，但是，如果遺傳算法直接針對(duì)圖4的有限元模型進(jìn)行尋優(yōu)，上千次的有限元模型的計(jì)算需要被執(zhí)行，從而導(dǎo)致反求效率極其低下。因此，筆者采用序列響應(yīng)面方法，在各次信賴域中構(gòu)建設(shè)計(jì)目標(biāo)的近似模型，而使得小種群遺傳算法在近似模型上尋優(yōu)，從而來提高整個(gè)參數(shù)反求過程的計(jì)算效率。反求過程即是使仿真曲線與該碰撞力曲線相逼近的過程，其逼近示意如圖5所示。

圖5 仿真與實(shí)驗(yàn)曲線逼近示意圖

由圖5可知，當(dāng)實(shí)驗(yàn)曲線上的插值點(diǎn)較多時(shí)，要使得仿真曲線上的對(duì)應(yīng)點(diǎn)與其逼近，將會(huì)調(diào)用大量的有限元計(jì)算，這就需要在仿真曲線與實(shí)驗(yàn)曲線的對(duì)應(yīng)插值點(diǎn)之間都構(gòu)建一個(gè)近似模型來解決該問題。由此可設(shè)立反求過程中的目標(biāo)函數(shù)為

式中，K為插值曲線的點(diǎn)的個(gè)數(shù)；x為設(shè)計(jì)變量，x＝A，B，C，m，n；Wi為權(quán)重系數(shù)；fi（x）為響應(yīng)面近似模型的計(jì)算值；Ti為實(shí)驗(yàn)測(cè)試點(diǎn)的值；Si為殘差比例系數(shù)；ei（x）為殘差值。

因此，DP800高強(qiáng)度鋼板的材料參數(shù)反求問題即轉(zhuǎn)化為以材料參數(shù)A、B、C、m、n為設(shè)計(jì)變量，以式（9）為設(shè)計(jì)目標(biāo)的優(yōu)化設(shè)計(jì)問題。其數(shù)學(xué)模型可以表示為

DP800薄壁梁材料參數(shù)反求過程可以表示為：首先對(duì)DP800梁進(jìn)行碰撞試驗(yàn)，得到碰撞力與變形長(zhǎng)度的試驗(yàn)曲線；然后，在整個(gè)設(shè)計(jì)空間中采用D最優(yōu)試驗(yàn)設(shè)計(jì)采樣樣本點(diǎn)，從而構(gòu)建DP800薄壁梁的各個(gè)試驗(yàn)曲線點(diǎn)的近似模型，并按式（9）形成目標(biāo)函數(shù)；最終，通過在小種群遺傳算法尋優(yōu)及信賴域的移動(dòng)和縮放來不斷更新近似模型，直到優(yōu)化收斂為止。DP800高強(qiáng)度鋼材料參數(shù)反求的設(shè)計(jì)流程如圖6所示。

圖6 DP800材料參數(shù)反求流程圖

2.3 DP800材料參數(shù)的反求結(jié)果及分析

基于序列響應(yīng)面方法在信賴域內(nèi)構(gòu)建了如式（10）所示的目標(biāo)函數(shù)，由于小種群遺傳算法是通過交叉、變異等方式進(jìn)行尋優(yōu)的，不需要梯度信息，且在處理設(shè)計(jì)響應(yīng)非線性程度高、不連續(xù)問題時(shí)，具有較好的全局最優(yōu)搜尋能力，它也可以避免計(jì)算過早的收斂，并且在最優(yōu)區(qū)域附近比普通遺傳算法有更好的收斂性［13，15］，因此，本文采用小種群遺傳算法迭代優(yōu)化。經(jīng)過8次設(shè)計(jì)空間的更新優(yōu)化設(shè)計(jì)后，目標(biāo)函數(shù)F收斂。DP800高強(qiáng)度鋼板的材料參數(shù)A、B、C、m、n的初始值、8次優(yōu)化迭代過程中的取值、設(shè)計(jì)目標(biāo)F的初始值及最優(yōu)值如表1所示。

由表1分析可知，優(yōu)化設(shè)計(jì)過程經(jīng)過8次響應(yīng)面近似模型的更新與優(yōu)化，設(shè)計(jì)目標(biāo)F的收斂曲線如圖7所示。材料參數(shù)A、B、C、m、n的信賴域的空間也逐漸縮小，材料參數(shù)反求過程共調(diào)用了81次DP800薄壁梁有限元計(jì)算模型，與直接采用傳統(tǒng)優(yōu)化算法的反求方法相比，序列響應(yīng)面方法所需調(diào)用正問題的次數(shù)極大地減少，從而大大提高了DP800高強(qiáng)度鋼板材料參數(shù)反求的效率。采用反求得到的材料參數(shù)A、B、C、m、n對(duì)薄壁梁進(jìn)行有限元仿真計(jì)算，從而得到仿真運(yùn)行的碰撞力與變形量曲線，它與試驗(yàn)曲線的對(duì)比如圖8所示。

由圖8的試驗(yàn)曲線與通過反求得到的材料參數(shù)A、B、C、m、n計(jì)算曲線的對(duì)比可以看出，兩條曲線整體趨勢(shì)上基本一致，各個(gè)峰值的位置也基本相同，從而驗(yàn)證了基于序列響應(yīng)面的材料參數(shù)反求方法以及基于該方法得到的材料參數(shù)值的可行性與正確性。

表1 材料參數(shù)及設(shè)計(jì)目標(biāo)的優(yōu)化歷程

3 結(jié)束語

本文提出了將序列響應(yīng)面方法與小種群遺傳算法相結(jié)合的材料參數(shù)反求方法，將材料參數(shù)反求問題轉(zhuǎn)變?yōu)橐粌?yōu)化問題，通過信賴域的縮放和更新來不斷地提高近似模型的精度，使整個(gè)反求過程的效率得到大大提高，通過對(duì)比反求得到的材料參數(shù)的碰撞力與變形量仿真曲線與試驗(yàn)曲線可知，該方法反求的材料參數(shù)具有較高的精度，能夠較好地達(dá)到工程應(yīng)用的要求。同時(shí)，該方法能有效地減少材料參數(shù)獲取過程中的試驗(yàn)次數(shù)，節(jié)約成本。

圖7 目標(biāo)函數(shù)收斂曲線

圖8 試驗(yàn)與反求參數(shù)仿真曲線對(duì)比圖

［1］ 鐘志華，張維剛，曹立波，等.汽車碰撞安全技術(shù)［M］.北京：機(jī)械工業(yè)出版社，2003.

［2］ 曹銀鋒，李光耀，鐘志華.金屬成型材料參數(shù)的反求技術(shù)［J］.計(jì)算力學(xué)學(xué)報(bào)，2004，21（3）：292－296.

［3］ Kurtaran H，Eskandarian A，Marzouguiel D，et al.Crashworthiness Design Optimization Using Successive Response Surface Approximation［J］.Computational Mechanics，2002，29（4／5）：409－421.

［4］ Tarigopula V，Langseth M，Hopperstad O，et al.Axial Crushing of Thin－walled High－strength Steel Sections［J］.International Journal of Impact Engineering，2006，32（5）：847－882.

［5］ 鄭剛，李光耀，孫光永.基于近似模型的拉延筋幾何參數(shù)反求［J］.中國(guó)機(jī)械工程，2006，17（19）：1988－1992.

［6］ Nagel G M，Thambiratnam D P.A Numerical Study on the Impact Response and Energy Absorption of Tapered Thin－ walled Tubes［J］.International Journal of Mechanical Sciences，2004，46（2）：201－216.

［7］ Abedrabbo N，Mayer R，Thompson A，et al.Crash Response of Advanced High－strength Steel Tubes：Experiment and Model［J］.International Journal of Impact Engineering，2009，36（8）：1044－1057.

［8］ Han Haipeng，Taheri F，Pegg N.Quasi－static and Dynamic Crushing Behaviors of Aluminum and Steel Tubes with a Cutout［J］.Thin－walled Structures，2007，45（3）：283－300.

［9］ 譚飛，韓旭.基于代理模型的功能梯度梁的材料特性參數(shù)反求［J］.復(fù)合材料學(xué)報(bào)，2008，25（5）：176－180.

［10］ 劉迪輝，汪晨，李光耀.薄壁鋼管材料參數(shù)反求［J］.中國(guó)機(jī)械工程，2008，19（6）：688－690.

［11］ Khalfallah H，Bel Hadj Salah H，Dogui A.Anisotropic Parameter Identification Using Inhomogeneous Tensile Test［J］.European Journal of Mechanics A／Solids，2002，21（6）：927－942.

［12］ Gelin J C，Ghouati O.An Inverse Method for Material Parameters Estimation in the Inelastic［J］.Computational Mechanics，1995，16（3）：143－150.

［13］ 韓利芬，高暉，李光耀，等.神經(jīng)網(wǎng)絡(luò)與遺傳算法在拉延筋參數(shù)反求中的應(yīng)用［J］.機(jī)械工程學(xué)報(bào)，2005，41（5）：172－176.

［14］ Redhe M，F(xiàn)orsberge J，Janssone T，et al.Using the Response Surface Methodology and the D－optimality Criterion in Crashworthiness Related Problems［J］.Structural and Multidisciplinary Optimization，2002，24（3）：185－194.

［15］ Liu G R，Han X.Computational Inverse Techniques in Nondestructive Evaluation［M］.Boca Raton：CRC Press，2003.