王国春 成艾国 顾纪超 宋 凯 钟志华

湖南大学汽车车身先进设计制造国家重点实验室，长沙，410082

0 引言

汽车耐撞性分析是一个涉及几何非线性、材料非线性和边界条件非线性的动力学问题[1]。有限元仿真技术在汽车耐撞性分析中已被广泛应用，可以通过仿真计算得到多个方案的计算结果。然而，汽车耐撞性分析问题本身依然是一个黑匣子问题，而且由于方案计算耗时长，限制了提交计算的方案数量，也无法利用全局优化算法进行迭代计算。近似模型又叫代理模型，它使用由设计参数构建的简单的函数来代替耗时的黑匣子问题，使得这些问题的优化设计成为可能，其最大的特点是计算速度快。但是，近似模型也存在两个缺陷：一是在模拟不同问题甚至同一问题的不同参数时，近似模型的精度差别很大，有时甚至会得到错误的结果；二是近似模型在模拟较大设计空间和较多设计变量的问题时局部空间的精度不高，因而阻碍了它们在优化设计中的应用[2－3]。

Shan等[4]开发了一种可以系统地将设计空间缩小到一个相对小的区间的方法；Gu等[5]提出了基于混合代理模型的优化方法，并行采用多种近似模型进行计算，对结果进行比较，逐渐缩小搜索空间以寻求最优值。本文参考混合近似模型和缩小空间的思想，采用三种各具特点的近似模型——Kriging模型、径向基函数模型和二阶多项式响应面模型进行搜索，同时采用传统的应用全局优化算法优化近似模型的方法，综合考虑耐撞性问题的特点，结合两种优化流程，在三种类型的变量空间搜索最优值，在最优点进行敏感空间的迭代计算，最终得到耐撞性目标——B柱加速度峰值最小的前部零件厚度组合，以较高的效率得到了满足工程精度和效率的最优方案，同时在一定程度上消除了由于近似模型选择不当而导致结果错误的现象。

1 三种近似模型的基本理论

1.1 Kriging模型

Kriging模型采用了插值方法，它拟合的曲面通过所有的样本点。从统计意义上说，这是一种从变量相关性和变异性出发，在有限区域内对区域化变量的取值进行无偏、最优估计的方法[6]。Kriging近似模型的基本理论可简述如下[7]。

它由全局模型与局部偏差叠加而成，可表示为

式中，Y（X）为未知的近似模型；f（X）为已知的多项式函数；Z（X）为均值为零、方差为σ2、协方差不为零的随机过程。

f（X）提供了设计空间的全局近似模型，一般情况下可取为常数β，而Z（X）则在全局模型的基础上创建了局部偏差。

1.2 径向基函数模型

径向基函数（radial basis function，RBF）模型与Kriging模型相似，径向基函数也是由一系列使用相同方法建立的函数构造的，只是基函数不同，其中基函数可以是线性的也可以是非线性的，不同的基函数给出的精度不同。本文采用欧拉范数作为基函数，此外还有很多常用的基函数，比如线性函数、立方函数、高斯函数等，同时使用者可以开发新的基函数来满足特殊需求，一些常用的基函数及其性能比较可参见文献[8]。

1.3 最小二乘响应面方法

最小二乘响应面法（least squares regression，LSR）用一个确定的显式函数（通常采用二次函数）来近似表达结构参数与结构响应之间的关系，应用最小二乘法确定显式函数中的待定参数。

2 拉丁超立方试验设计

在近似模型的构造过程中，设计点的选择是很重要的。任意选择的设计点会导致不精确的模型，试验设计的理论可以帮助确定分布合理的设计点。拉丁超立方试验设计（Latin hypercube design，LHD）被称为是一种“填充空间的设计”，该方法让每个因子都有相同的水平，并且对因素的每个水平所做的试验次数相等，以实现设计空间中的均匀抽样，从而指定用于定义设计矩阵的n个点，是一种研究多因素、采样大型空间的试验设计方法，该方法效率高、采样均衡性好[9]。如图1、图2所示，对于一个两因素4次试验的试验设计，利用全因子试验方法只能考察因素的2个水平，利用拉丁超立方法，则能考察4个水平，并且因素的每个水平都做了一次试验，显然，试验点在空间的排列比全因子法要均匀。

图1 全因子试验设计

图2 拉丁超立方试验设计

3 计算实例

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

整车正面碰撞有限元模型共有112万个节点，121万个单元，按照试验标准配重后的质量为1740kg，以13.8m／s的初始速度撞击刚性墙，碰撞仿真过程计算时间为100ms。本文中需要计算的有限元模型如图3所示。

图3 整车正面碰撞有限元模型

3.2 近似模型拟合与优化计算

在整车正面碰撞中，B柱的加速度峰值越低，车身吸收能量的时间越长，对乘员的损伤越小，汽车就越安全。本例是在保证各前部结构零件质量和其吸收的内能满足设定的约束条件的前提下，实现耐撞性重要指标B柱加速度峰值最小的优化设计。

优化问题的数学模型如下：

设计变量为T1～T10

式中，aB为B柱加速度最大值，要求优化计算求其最小值；m为图4中15个部件的质量之和，约束质量增加不超过初始方案的5%，为31.08kg，设为31kg；EI为15个部件吸收的内能之和，允许吸收的内能略有下降，设为69 000J；T1～T10为15个部件的10种板料厚度（对称件共用一个厚度变量）；T为各部件的初始厚度。

T9根据实际项目要求设定了不同的范围。所有零件的厚度具体取值范围、响应的初始值及经过有限元验证的优化值如表1所示。

表1 变量与响应的初始值、上下限、优化值

图4 优化厚度的零件

优化流程如图5所示。初始样本点应用LS－DYNA软件，通过25台PC同时计算，每个正面碰撞有限元模型样本耗时24h左右，经过96h左右的连续计算，得到全部样本的仿真结果。收敛判断的准则是：如果目标值差值的绝对值小于0.1g或者目标值大于前一轮的目标值，都认为收敛。

传统的采用近似模型的优化方法是：构建近似模型，采用全局优化算法进行计算，输出最优解，调用有限元模型进行验证。这种方法的优点是能够快速找到近似模型的全局最优解，但是，由于近似模型在设计空间较大时局部空间精度不同，如果全局最优解恰好存在于精度不高的设计空间，则得到的解就会存在相对误差较大甚至错误的情况，另外，利用单一模型计算时，有时会出现某个响应在搜寻到的最优点处相对误差偏大的现象。

图5 优化流程图

本文改良了传统方法，对三种近似模型进行全局优化计算；同时为避免近似模型在某一空间可能存在相对误差偏大的问题，在三种类型的变量敏感空间中生成密集的样本，再分别构建三种近似模型进行计算，在所得的样本中进行搜索，选取满足约束条件的目标值最小解；最后选出所有的满足约束条件的各近似模型的目标值最小解，一并进行有限元计算验证，找到有限元模型中的满足约束的目标值最小解，然后在三种敏感空间中重新密集撒点，循环比较搜寻最优点，直至收敛为止。

全局空间指由原有变量范围构成的设计空间；局部空间指由满足约束要求的目标值最小的10个初始有限元模型样本点的变量范围构成的设计空间；敏感空间指由满足约束的目标值最小点的变量上下浮动一定比例构造的设计空间，文中生成了变量分别上下浮动2.5%、5%、10%的三种由小到大的敏感空间，以进行逐步扩大的目标搜索。

3.3 优化过程分析

表2～表4列出了第一轮的计算结果，其中，序号是在提交LS－DYNA计算过程中的顺序，以使三个表格的方案与验证结果能够一一对应，编号的先后没有特殊意义。表2列出的是通过在三种类型的空间中进行密集撒点，应用近似模型进行计算，并按照约束条件搜索得到的每种近似模型的目标值最小的方案。这里需要指出的是，对于三个不同范围的敏感空间，将按照敏感空间的类型选出每种近似模型的目标值最小的方案，如表2中选出的7、8、9号方案。表3列出的是通过模拟退火全局优化算法对三种近似模型在全局空间按照约束条件进行搜索得到的目标值最小的方案。最终将通过两种优化流程、三种空间类型、三种近似模型技术选出的12组方案代入有限元模型进行验证计算，得到的结果如表4所示，可看出，3号方案为第一轮计算中目标值aB最小的方案，其值为41.17g，3号方案是在全局空间进行密集撒点，通过RBF近似模型计算得到的方案，aB的相对误差只有0.81%。相对误差表示了每个样本点的近似模型相对于物理模型的计算精度，其表达式如下[10]：

er为每个样本点的相对误差。其中，当｜y｜≠0时，ε＝0；当｜y｜＝0时，ε＝0.01。

表2 第一轮近似模型计算结果

表3 第一轮近似模型全局优化结果

表4 第一轮计算结果验证

由表4可以发现，同一种近似模型在同一样本点，不同响应的相对误差大小不同，在不同样本点同一响应的相对误差有时差别很大。从第一轮12个样本点的验证结果来看，RBF、Kriging、LSR模型的近似精度依次降低。利用自适应模拟退火算法（adaptive simulation annealing，ASA）在全局空间对LSR近似模型进行优化计算的目标值最小方案只有32.02g，但是其相对误差却高达28.64%。在对RBF、Kriging近似模型进行全局优化计算时其目标值分别为39.91g和40.56g，虽然均小于3号方案近似模型的计算值，但是其相对误差分别为6.72%和4.23%，物理模型的计算值为42.79g和42.35g，都大于3号方案物理模型的计算值。这一点并不能说明直接针对近似模型进行优化计算的方法不好，只能说明当变量处于不同的空间位置时，其精度有时会差别很大，这也就是常说的空间依赖性，这也具有一定的偶然性。如果在最优点的相对误差非常小或者在可接受的范围内，则这种方法的寻优效率是最高的。

在完成第一轮计算得到目标值最小的3号方案后，在其变量的三个敏感空间重新进行密集撒点，按照约束条件搜索得到的三个敏感空间三种近似模型的计算结果如表5所示。选出三种近似模型目标值最小的三个方案代入有限元模型进行计算后，得到第二轮计算的目标值最小的1号方案，如表6所示，它是在第一轮目标值最小方案变量变动±2.5%后的敏感空间中由Kriging模型计算得到的，aB为40.54g，相对误差为1.58%，比第一轮减小了0.63g，不满足收敛条件，继续进行迭代。

在第二轮得到的目标值最小的1号方案的变量的三个敏感空间中重新进行密集撒点，计算结果如表7、表8所示，得到第三轮目标值最小方案，它是在第二轮目标值最小方案变量变动±10%后的敏感空间中由LSR模型计算得到的。aB为41.99g，相对误差为9.09%，比第二轮增大了1.45g，满足了收敛了条件，因而输出第二轮的1号方案作为最优结果。

表5 第二轮近似模型计算结果

表6 第二轮计算结果验证

表7 第三轮近似模型计算结果

表8 第三轮计算结果验证

B柱加速度的初始值和最终优化值的对比如图6所示。

图6 B柱加速度初始值优化结果对比

如表1所示，aB由初始值46.34g减小到40.54g，减小了12.52%，有较大程度的降低；m由29.60kg增大到30.178kg，增大了1.95%；EI由70 620J增大到71 046J；两个约束条件均满足初始的约束条件控制要求。

4 结语

本文结合拉丁方试验设计、三种近似模型技术、正面碰撞有限元仿真技术、模拟退火全局优化算法、密集撒点优化搜索方法、三种空间两种流程的同步优化方法，实现了正面碰撞耐撞性的优化设计，较好地改善了耐撞性指标，得到了最优的板厚组合。这一方法的应用将能够在设计阶段较好地改善汽车耐撞性设计，减少后期的实验次数，缩短开发周期。对于零件的最终厚度组合方案，传统的方法通常是把优化结果进行圆整，保留一位小数，得到最优化的厚度组合。这是不符合实际工程情况的，因为：工程中所用板材厚度是一个连续值，厚度是存在公差的；零件在冲压成形的过程中存在不同程度的减薄，通常要求减薄率大于5%而小于20%。因此，考虑响应对厚度的敏感性，应该综合考虑以上两因素，同时应用稳健性优化方法，才能得到最好的工程解决方案，这一点将在后期的研究中进行研究。

[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模型的结构可靠性分析[J].计算力学学报，2006，23（2）：175-176.

[7]丁彦闯，兆文忠.基于Kriging模型的焊接构架抗疲劳优化设计[J].大连交通大学学报，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]方开泰.均匀试验设计的理论、方法和应用——历史回顾[J].数理统计与管理，2004，23（3）：69-80.

[10]谢延敏，于沪平，陈军，等.基于Kriging模型的可靠度计算[J].上海交通大学学报，2007，41（2）：177-178.