基于混合近似模型的汽車正面碰撞耐撞性優(yōu)化設(shè)計(jì)

王國(guó)春 成艾國(guó) 顧紀(jì)超 宋 凱 鐘志華

湖南大學(xué)汽車車身先進(jìn)設(shè)計(jì)制造國(guó)家重點(diǎn)實(shí)驗(yàn)室，長(zhǎng)沙，410082

0 引言

汽車耐撞性分析是一個(gè)涉及幾何非線性、材料非線性和邊界條件非線性的動(dòng)力學(xué)問(wèn)題[1]。有限元仿真技術(shù)在汽車耐撞性分析中已被廣泛應(yīng)用，可以通過(guò)仿真計(jì)算得到多個(gè)方案的計(jì)算結(jié)果。然而，汽車耐撞性分析問(wèn)題本身依然是一個(gè)黑匣子問(wèn)題，而且由于方案計(jì)算耗時(shí)長(zhǎng)，限制了提交計(jì)算的方案數(shù)量，也無(wú)法利用全局優(yōu)化算法進(jìn)行迭代計(jì)算。近似模型又叫代理模型，它使用由設(shè)計(jì)參數(shù)構(gòu)建的簡(jiǎn)單的函數(shù)來(lái)代替耗時(shí)的黑匣子問(wèn)題，使得這些問(wèn)題的優(yōu)化設(shè)計(jì)成為可能，其最大的特點(diǎn)是計(jì)算速度快。但是，近似模型也存在兩個(gè)缺陷：一是在模擬不同問(wèn)題甚至同一問(wèn)題的不同參數(shù)時(shí)，近似模型的精度差別很大，有時(shí)甚至?xí)玫藉e(cuò)誤的結(jié)果；二是近似模型在模擬較大設(shè)計(jì)空間和較多設(shè)計(jì)變量的問(wèn)題時(shí)局部空間的精度不高，因而阻礙了它們?cè)趦?yōu)化設(shè)計(jì)中的應(yīng)用[2－3]。

Shan等[4]開(kāi)發(fā)了一種可以系統(tǒng)地將設(shè)計(jì)空間縮小到一個(gè)相對(duì)小的區(qū)間的方法；Gu等[5]提出了基于混合代理模型的優(yōu)化方法，并行采用多種近似模型進(jìn)行計(jì)算，對(duì)結(jié)果進(jìn)行比較，逐漸縮小搜索空間以尋求最優(yōu)值。本文參考混合近似模型和縮小空間的思想，采用三種各具特點(diǎn)的近似模型——Kriging模型、徑向基函數(shù)模型和二階多項(xiàng)式響應(yīng)面模型進(jìn)行搜索，同時(shí)采用傳統(tǒng)的應(yīng)用全局優(yōu)化算法優(yōu)化近似模型的方法，綜合考慮耐撞性問(wèn)題的特點(diǎn)，結(jié)合兩種優(yōu)化流程，在三種類型的變量空間搜索最優(yōu)值，在最優(yōu)點(diǎn)進(jìn)行敏感空間的迭代計(jì)算，最終得到耐撞性目標(biāo)——B柱加速度峰值最小的前部零件厚度組合，以較高的效率得到了滿足工程精度和效率的最優(yōu)方案，同時(shí)在一定程度上消除了由于近似模型選擇不當(dāng)而導(dǎo)致結(jié)果錯(cuò)誤的現(xiàn)象。

1 三種近似模型的基本理論

1.1 Kriging模型

Kriging模型采用了插值方法，它擬合的曲面通過(guò)所有的樣本點(diǎn)。從統(tǒng)計(jì)意義上說(shuō)，這是一種從變量相關(guān)性和變異性出發(fā)，在有限區(qū)域內(nèi)對(duì)區(qū)域化變量的取值進(jìn)行無(wú)偏、最優(yōu)估計(jì)的方法[6]。Kriging近似模型的基本理論可簡(jiǎn)述如下[7]。

它由全局模型與局部偏差疊加而成，可表示為

式中，Y（X）為未知的近似模型；f（X）為已知的多項(xiàng)式函數(shù)；Z（X）為均值為零、方差為σ2、協(xié)方差不為零的隨機(jī)過(guò)程。

f（X）提供了設(shè)計(jì)空間的全局近似模型，一般情況下可取為常數(shù)β，而Z（X）則在全局模型的基礎(chǔ)上創(chuàng)建了局部偏差。

1.2 徑向基函數(shù)模型

徑向基函數(shù)（radial basis function，RBF）模型與Kriging模型相似，徑向基函數(shù)也是由一系列使用相同方法建立的函數(shù)構(gòu)造的，只是基函數(shù)不同，其中基函數(shù)可以是線性的也可以是非線性的，不同的基函數(shù)給出的精度不同。本文采用歐拉范數(shù)作為基函數(shù)，此外還有很多常用的基函數(shù)，比如線性函數(shù)、立方函數(shù)、高斯函數(shù)等，同時(shí)使用者可以開(kāi)發(fā)新的基函數(shù)來(lái)滿足特殊需求，一些常用的基函數(shù)及其性能比較可參見(jiàn)文獻(xiàn)[8]。

1.3 最小二乘響應(yīng)面方法

最小二乘響應(yīng)面法（least squares regression，LSR）用一個(gè)確定的顯式函數(shù)（通常采用二次函數(shù)）來(lái)近似表達(dá)結(jié)構(gòu)參數(shù)與結(jié)構(gòu)響應(yīng)之間的關(guān)系，應(yīng)用最小二乘法確定顯式函數(shù)中的待定參數(shù)。

2 拉丁超立方試驗(yàn)設(shè)計(jì)

在近似模型的構(gòu)造過(guò)程中，設(shè)計(jì)點(diǎn)的選擇是很重要的。任意選擇的設(shè)計(jì)點(diǎn)會(huì)導(dǎo)致不精確的模型，試驗(yàn)設(shè)計(jì)的理論可以幫助確定分布合理的設(shè)計(jì)點(diǎn)。拉丁超立方試驗(yàn)設(shè)計(jì)（Latin hypercube design，LHD）被稱為是一種“填充空間的設(shè)計(jì)”，該方法讓每個(gè)因子都有相同的水平，并且對(duì)因素的每個(gè)水平所做的試驗(yàn)次數(shù)相等，以實(shí)現(xiàn)設(shè)計(jì)空間中的均勻抽樣，從而指定用于定義設(shè)計(jì)矩陣的n個(gè)點(diǎn)，是一種研究多因素、采樣大型空間的試驗(yàn)設(shè)計(jì)方法，該方法效率高、采樣均衡性好[9]。如圖1、圖2所示，對(duì)于一個(gè)兩因素4次試驗(yàn)的試驗(yàn)設(shè)計(jì)，利用全因子試驗(yàn)方法只能考察因素的2個(gè)水平，利用拉丁超立方法，則能考察4個(gè)水平，并且因素的每個(gè)水平都做了一次試驗(yàn)，顯然，試驗(yàn)點(diǎn)在空間的排列比全因子法要均勻。

圖1 全因子試驗(yàn)設(shè)計(jì)

圖2 拉丁超立方試驗(yàn)設(shè)計(jì)

3 計(jì)算實(shí)例

3.1 整車正面碰撞有限元模型

整車正面碰撞有限元模型共有112萬(wàn)個(gè)節(jié)點(diǎn)，121萬(wàn)個(gè)單元，按照試驗(yàn)標(biāo)準(zhǔn)配重后的質(zhì)量為1740kg，以13.8m／s的初始速度撞擊剛性墻，碰撞仿真過(guò)程計(jì)算時(shí)間為100ms。本文中需要計(jì)算的有限元模型如圖3所示。

圖3 整車正面碰撞有限元模型

3.2 近似模型擬合與優(yōu)化計(jì)算

在整車正面碰撞中，B柱的加速度峰值越低，車身吸收能量的時(shí)間越長(zhǎng)，對(duì)乘員的損傷越小，汽車就越安全。本例是在保證各前部結(jié)構(gòu)零件質(zhì)量和其吸收的內(nèi)能滿足設(shè)定的約束條件的前提下，實(shí)現(xiàn)耐撞性重要指標(biāo)B柱加速度峰值最小的優(yōu)化設(shè)計(jì)。

優(yōu)化問(wèn)題的數(shù)學(xué)模型如下：

設(shè)計(jì)變量為T1～T10

式中，aB為B柱加速度最大值，要求優(yōu)化計(jì)算求其最小值；m為圖4中15個(gè)部件的質(zhì)量之和，約束質(zhì)量增加不超過(guò)初始方案的5%，為31.08kg，設(shè)為31kg；EI為15個(gè)部件吸收的內(nèi)能之和，允許吸收的內(nèi)能略有下降，設(shè)為69 000J；T1～T10為15個(gè)部件的10種板料厚度（對(duì)稱件共用一個(gè)厚度變量）；T為各部件的初始厚度。

T9根據(jù)實(shí)際項(xiàng)目要求設(shè)定了不同的范圍。所有零件的厚度具體取值范圍、響應(yīng)的初始值及經(jīng)過(guò)有限元驗(yàn)證的優(yōu)化值如表1所示。

表1 變量與響應(yīng)的初始值、上下限、優(yōu)化值

圖4 優(yōu)化厚度的零件

優(yōu)化流程如圖5所示。初始樣本點(diǎn)應(yīng)用LS－DYNA軟件，通過(guò)25臺(tái)PC同時(shí)計(jì)算，每個(gè)正面碰撞有限元模型樣本耗時(shí)24h左右，經(jīng)過(guò)96h左右的連續(xù)計(jì)算，得到全部樣本的仿真結(jié)果。收斂判斷的準(zhǔn)則是：如果目標(biāo)值差值的絕對(duì)值小于0.1g或者目標(biāo)值大于前一輪的目標(biāo)值，都認(rèn)為收斂。

傳統(tǒng)的采用近似模型的優(yōu)化方法是：構(gòu)建近似模型，采用全局優(yōu)化算法進(jìn)行計(jì)算，輸出最優(yōu)解，調(diào)用有限元模型進(jìn)行驗(yàn)證。這種方法的優(yōu)點(diǎn)是能夠快速找到近似模型的全局最優(yōu)解，但是，由于近似模型在設(shè)計(jì)空間較大時(shí)局部空間精度不同，如果全局最優(yōu)解恰好存在于精度不高的設(shè)計(jì)空間，則得到的解就會(huì)存在相對(duì)誤差較大甚至錯(cuò)誤的情況，另外，利用單一模型計(jì)算時(shí)，有時(shí)會(huì)出現(xiàn)某個(gè)響應(yīng)在搜尋到的最優(yōu)點(diǎn)處相對(duì)誤差偏大的現(xiàn)象。

圖5 優(yōu)化流程圖

本文改良了傳統(tǒng)方法，對(duì)三種近似模型進(jìn)行全局優(yōu)化計(jì)算；同時(shí)為避免近似模型在某一空間可能存在相對(duì)誤差偏大的問(wèn)題，在三種類型的變量敏感空間中生成密集的樣本，再分別構(gòu)建三種近似模型進(jìn)行計(jì)算，在所得的樣本中進(jìn)行搜索，選取滿足約束條件的目標(biāo)值最小解；最后選出所有的滿足約束條件的各近似模型的目標(biāo)值最小解，一并進(jìn)行有限元計(jì)算驗(yàn)證，找到有限元模型中的滿足約束的目標(biāo)值最小解，然后在三種敏感空間中重新密集撒點(diǎn)，循環(huán)比較搜尋最優(yōu)點(diǎn)，直至收斂為止。

全局空間指由原有變量范圍構(gòu)成的設(shè)計(jì)空間；局部空間指由滿足約束要求的目標(biāo)值最小的10個(gè)初始有限元模型樣本點(diǎn)的變量范圍構(gòu)成的設(shè)計(jì)空間；敏感空間指由滿足約束的目標(biāo)值最小點(diǎn)的變量上下浮動(dòng)一定比例構(gòu)造的設(shè)計(jì)空間，文中生成了變量分別上下浮動(dòng)2.5%、5%、10%的三種由小到大的敏感空間，以進(jìn)行逐步擴(kuò)大的目標(biāo)搜索。

3.3 優(yōu)化過(guò)程分析

表2～表4列出了第一輪的計(jì)算結(jié)果，其中，序號(hào)是在提交LS－DYNA計(jì)算過(guò)程中的順序，以使三個(gè)表格的方案與驗(yàn)證結(jié)果能夠一一對(duì)應(yīng)，編號(hào)的先后沒(méi)有特殊意義。表2列出的是通過(guò)在三種類型的空間中進(jìn)行密集撒點(diǎn)，應(yīng)用近似模型進(jìn)行計(jì)算，并按照約束條件搜索得到的每種近似模型的目標(biāo)值最小的方案。這里需要指出的是，對(duì)于三個(gè)不同范圍的敏感空間，將按照敏感空間的類型選出每種近似模型的目標(biāo)值最小的方案，如表2中選出的7、8、9號(hào)方案。表3列出的是通過(guò)模擬退火全局優(yōu)化算法對(duì)三種近似模型在全局空間按照約束條件進(jìn)行搜索得到的目標(biāo)值最小的方案。最終將通過(guò)兩種優(yōu)化流程、三種空間類型、三種近似模型技術(shù)選出的12組方案代入有限元模型進(jìn)行驗(yàn)證計(jì)算，得到的結(jié)果如表4所示，可看出，3號(hào)方案為第一輪計(jì)算中目標(biāo)值aB最小的方案，其值為41.17g，3號(hào)方案是在全局空間進(jìn)行密集撒點(diǎn)，通過(guò)RBF近似模型計(jì)算得到的方案，aB的相對(duì)誤差只有0.81%。相對(duì)誤差表示了每個(gè)樣本點(diǎn)的近似模型相對(duì)于物理模型的計(jì)算精度，其表達(dá)式如下[10]：

er為每個(gè)樣本點(diǎn)的相對(duì)誤差。其中，當(dāng)｜y｜≠0時(shí)，ε＝0；當(dāng)｜y｜＝0時(shí)，ε＝0.01。

表2 第一輪近似模型計(jì)算結(jié)果

表3 第一輪近似模型全局優(yōu)化結(jié)果

表4 第一輪計(jì)算結(jié)果驗(yàn)證

由表4可以發(fā)現(xiàn)，同一種近似模型在同一樣本點(diǎn)，不同響應(yīng)的相對(duì)誤差大小不同，在不同樣本點(diǎn)同一響應(yīng)的相對(duì)誤差有時(shí)差別很大。從第一輪12個(gè)樣本點(diǎn)的驗(yàn)證結(jié)果來(lái)看，RBF、Kriging、LSR模型的近似精度依次降低。利用自適應(yīng)模擬退火算法（adaptive simulation annealing，ASA）在全局空間對(duì)LSR近似模型進(jìn)行優(yōu)化計(jì)算的目標(biāo)值最小方案只有32.02g，但是其相對(duì)誤差卻高達(dá)28.64%。在對(duì)RBF、Kriging近似模型進(jìn)行全局優(yōu)化計(jì)算時(shí)其目標(biāo)值分別為39.91g和40.56g，雖然均小于3號(hào)方案近似模型的計(jì)算值，但是其相對(duì)誤差分別為6.72%和4.23%，物理模型的計(jì)算值為42.79g和42.35g，都大于3號(hào)方案物理模型的計(jì)算值。這一點(diǎn)并不能說(shuō)明直接針對(duì)近似模型進(jìn)行優(yōu)化計(jì)算的方法不好，只能說(shuō)明當(dāng)變量處于不同的空間位置時(shí)，其精度有時(shí)會(huì)差別很大，這也就是常說(shuō)的空間依賴性，這也具有一定的偶然性。如果在最優(yōu)點(diǎn)的相對(duì)誤差非常小或者在可接受的范圍內(nèi)，則這種方法的尋優(yōu)效率是最高的。

在完成第一輪計(jì)算得到目標(biāo)值最小的3號(hào)方案后，在其變量的三個(gè)敏感空間重新進(jìn)行密集撒點(diǎn)，按照約束條件搜索得到的三個(gè)敏感空間三種近似模型的計(jì)算結(jié)果如表5所示。選出三種近似模型目標(biāo)值最小的三個(gè)方案代入有限元模型進(jìn)行計(jì)算后，得到第二輪計(jì)算的目標(biāo)值最小的1號(hào)方案，如表6所示，它是在第一輪目標(biāo)值最小方案變量變動(dòng)±2.5%后的敏感空間中由Kriging模型計(jì)算得到的，aB為40.54g，相對(duì)誤差為1.58%，比第一輪減小了0.63g，不滿足收斂條件，繼續(xù)進(jìn)行迭代。

在第二輪得到的目標(biāo)值最小的1號(hào)方案的變量的三個(gè)敏感空間中重新進(jìn)行密集撒點(diǎn)，計(jì)算結(jié)果如表7、表8所示，得到第三輪目標(biāo)值最小方案，它是在第二輪目標(biāo)值最小方案變量變動(dòng)±10%后的敏感空間中由LSR模型計(jì)算得到的。aB為41.99g，相對(duì)誤差為9.09%，比第二輪增大了1.45g，滿足了收斂了條件，因而輸出第二輪的1號(hào)方案作為最優(yōu)結(jié)果。

表5 第二輪近似模型計(jì)算結(jié)果

表6 第二輪計(jì)算結(jié)果驗(yàn)證

表7 第三輪近似模型計(jì)算結(jié)果

表8 第三輪計(jì)算結(jié)果驗(yàn)證

B柱加速度的初始值和最終優(yōu)化值的對(duì)比如圖6所示。

圖6 B柱加速度初始值優(yōu)化結(jié)果對(duì)比

如表1所示，aB由初始值46.34g減小到40.54g，減小了12.52%，有較大程度的降低；m由29.60kg增大到30.178kg，增大了1.95%；EI由70 620J增大到71 046J；兩個(gè)約束條件均滿足初始的約束條件控制要求。

4 結(jié)語(yǔ)

本文結(jié)合拉丁方試驗(yàn)設(shè)計(jì)、三種近似模型技術(shù)、正面碰撞有限元仿真技術(shù)、模擬退火全局優(yōu)化算法、密集撒點(diǎn)優(yōu)化搜索方法、三種空間兩種流程的同步優(yōu)化方法，實(shí)現(xiàn)了正面碰撞耐撞性的優(yōu)化設(shè)計(jì)，較好地改善了耐撞性指標(biāo)，得到了最優(yōu)的板厚組合。這一方法的應(yīng)用將能夠在設(shè)計(jì)階段較好地改善汽車耐撞性設(shè)計(jì)，減少后期的實(shí)驗(yàn)次數(shù)，縮短開(kāi)發(fā)周期。對(duì)于零件的最終厚度組合方案，傳統(tǒng)的方法通常是把優(yōu)化結(jié)果進(jìn)行圓整，保留一位小數(shù)，得到最優(yōu)化的厚度組合。這是不符合實(shí)際工程情況的，因?yàn)椋汗こ讨兴冒宀暮穸仁且粋€(gè)連續(xù)值，厚度是存在公差的；零件在沖壓成形的過(guò)程中存在不同程度的減薄，通常要求減薄率大于5%而小于20%。因此，考慮響應(yīng)對(duì)厚度的敏感性，應(yīng)該綜合考慮以上兩因素，同時(shí)應(yīng)用穩(wěn)健性優(yōu)化方法，才能得到最好的工程解決方案，這一點(diǎn)將在后期的研究中進(jìn)行研究。

[1]Hong Z H.Finite Element Procedures for Contact－impact Problems[M].Cambridge：Oxford University Press，1993.

[2]Myers R H，Montgomery D C.Response Surface Methodology：Process and Product Optimization Using Designed Experiments[M].Toronto：John Wiley ＆Sons，2002.

[3]Chen V C P，Tsui K L，Barton R R，et al.A Review on Design，Modeling and Applications of Computer Experiments[J].IIE Transactions，2006，38：273-291.

[4]Shan S，Wang G G.Space Exploration and Global Optimization for Computationally Intensive Design Problems：a Rough Set Based Approach[J].Structural Multidisciplinary Optimization，2004，28：427-441.

[5]Gu J，Li G Y，Dong Z.Hybrid and Adaptive Metamodel Based Global Optimization[C]／／Proceedings of the ASME 2009International Design Engineering Technical Conferences ＆ Computers and Information in Engineering Conference，IDETC／CIE 2009.San Diego，California，2009：DETC2009－87121.

[6]張崎，李興斯.基于Kriging模型的結(jié)構(gòu)可靠性分析[J].計(jì)算力學(xué)學(xué)報(bào)，2006，23（2）：175-176.

[7]丁彥闖，兆文忠.基于Kriging模型的焊接構(gòu)架抗疲勞優(yōu)化設(shè)計(jì)[J].大連交通大學(xué)學(xué)報(bào)，2008，29（2）：7-11.

[8]Fang H B，Horstemeyer M F.Global Response Approximation with Radial Basis Functions[J].Engineering Optimization，2006，38（4）：407-424.

[9]方開(kāi)泰.均勻試驗(yàn)設(shè)計(jì)的理論、方法和應(yīng)用——?dú)v史回顧[J].數(shù)理統(tǒng)計(jì)與管理，2004，23（3）：69-80.

[10]謝延敏，于滬平，陳軍，等.基于Kriging模型的可靠度計(jì)算[J].上海交通大學(xué)學(xué)報(bào)，2007，41（2）：177-178.