基于通用生成函数法的CFRP防撞梁碰撞可靠性设计优化*

刘展鹏 周金宇

（1.盐城工学院，盐城 224000；2.金陵科技学院，南京 211169）

主题词：碳纤维增强复合材料 通用生成函数法 耐撞性 可靠性设计优化 非均匀聚类Kriging模型

1 前言

防撞梁在汽车的安全性设计中有着重要的地位，汽车防撞梁可以吸收碰撞时的能量，给予车身缓冲，对汽车的主体结构起到防护作用。汽车的碰撞过程是高度非线性的复杂过程，设计变量、随机变量与功能函数、目标函数之间不存在显式的函数关系式，常用的方法是构建功能函数和目标函数的代理模型，从而对汽车耐撞性与安全性进行优化。众多学者以此展开研究，仲伟东等[1]结合响应面模型和多目标遗传算法对防撞梁进行优化设计。王爽等[2]采用灰色关联及熵权法对碳纤维增强复合材料（Carbon Fibre Reinforced Plastics，CFRP）防撞梁进行优化设计。谢晖等[3]采用响应面法和自适应模拟退火（Adaptive Simulated Annealing，ASA）算法对防撞梁进行设计优化。张鑫等[4]采用Hyperkriging模型和遗传算法对防撞梁厚度进行优化。陈静等[5]通过对复合材料防撞梁厚度进行多目标优化，使得防撞梁质量明显减轻。Belingardi[6]从材料、结构两方面对汽车前防撞梁进行了优化，改善了前防撞梁的吸能效果。曹立波等[7]采用中心复合试验设计方法和自适应响应面法对防撞梁进行优化设计。童小伟等[8]采用Hyperkriging 方法结合遗传算法对碳纤维增强复合材料防撞梁进行优化设计。王庆等[9]采用Kriging 代理模型方法和遗传算法对碳纤维增强树脂复合材料防撞梁结构进行优化。Duan[10]等将响应面分析法（Response Surface Methodology，RSM）和非支配遗传算法（Non-dominated Sorting Genetic Algorithm-II，NSGA-Ⅱ）相结合，对复合材料防撞梁进行优化设计。Jalauddin[11]等研究不同材料的防撞梁的耐撞性和轻量化程度，最后得出碳纤维增强树脂复合材料的耐撞性与轻量化更好。Liu[12]等结合Kriging近似模型和多目标优化算法对防撞梁进行优化设计。J.Hilmann等[13]使用遗传算法对防撞梁进行低速碰撞研究并对防撞梁结构进行优化设计。Simon 等[14]使用LS-DYNA 对质量相同的钢制和铝合金制防撞梁进行低速碰撞仿真分析，结果表明铝合金防撞梁的比吸能更高。上述文献主要采用传统的优化算法进行确定性的优化设计，但在实际工程中存在大量的不确定性，所以对复合材料防撞梁进行碰撞结构可靠性分析十分必要。

本文采用二次样条曲线对防撞梁轮廓进行描述，轮廓的形状主要由各个控制点横坐标进行控制，故选取轮廓曲线控制点横坐标为设计变量，同时选取厚度与铺层角度作为设计变量，以此建立数学模型。由于目标函数及功能函数为隐式函数，故使用有限元建模与代理模型法拟合出显式关系式。最后在可靠性分析环节引入通用生成函数（Universal Generating Function，UGF）法，使用遗传算法进行优化。

2 CFRP防撞梁的碰撞RBDO数学模型

根据国家标准GB 17354—1998，以防撞梁轮廓曲线控制点横坐标xi（i=1,2）、厚度t以及复合材料铺层角度θj为设计参数，以防撞梁碰撞力峰值Fmax作为约束条件，以最大化比吸能ESEA为设计目标，构建防撞梁的碰撞可靠性优化设计（Reliability-Based Design Optimization，RBDO）RBDO数学模型。防撞梁的截面尺寸如图1所示。

图1 防撞梁截面

2.1 建立RBDO数学模型

为了使数学模型的表达更加清晰，铺层角度值采用公式进行转化，将角度值用0，1，2，3来表示。铺层角度θj可由公式表示为：

式中，xj+3的取值范围为{0，1，2，3}。

通过式（1）即可转换成铺层角度的取值。铺层角度如图2所示。

图2 铺层角度

复合材料的弹性模量和剪切模量是非常重要的参数，将弹性模量E1、E2和剪切模量G1、G2设定为随机参数。假设各随机参数相互独立，其中E1、E2服从正态分布，G1、G2服从随机分布，参数信息如表1所示。从耐撞性角度考虑，需要防撞梁的总吸能E尽可能大，碰撞力峰值Fmax尽可能小；从轻量化角度考虑，需要防撞梁质量m尽可能小。因此，为了同时满足防撞梁耐撞性以及轻量化的要求，将比吸能ESEA作为优化目标，以碰撞力峰值Fmax作为约束，构建数学模型如下：

表1 随机参数信息表

式中，d为确定性设计变量；P=(E1,E2,G1,G2)为随机参数向量；Pr{G(d,P)≤0} 为结构失效概率；Φ()为标准正态分布函数，许用可靠度指标βT=2；x3为防撞梁厚度。

2.2 耐撞性与轻量化评价

通常将吸能量H、比吸能ESEA、碰撞力峰值Fmax、平均碰撞力Fave、防撞梁侵入量D等参数作为结构耐撞性的评价参数。其中吸能量H是在碰撞过程中，防撞梁变形所吸收的能量，数值越大，说明防撞梁吸能效果越好；比吸能ESEA为防撞梁单位质量的吸能量，为总吸能量与质量的比值，优化比吸能可以同时优化耐撞性和实现轻量化，具有重大的意义。碰撞力峰值Fmax是碰撞过程中防撞梁结构所受的最大冲击力，数值越小，说明防撞梁所受的冲击力越小，安全性能越强。

3 CFRP防撞梁碰撞工况有限元仿真

3.1 CFRP防撞梁有限元建模

由于汽车CFRP防撞梁的轮廓曲线为左右对称的自由曲线，可用二次B样条（B-Spline）曲线进行描述，用控制点来改变轮廓曲线的形状。对B-Spline曲线的描述由MATLAB编程实现，二次均匀B-Spline曲线可以表示为：

式中，u为描述曲线的控制点横坐标;n为控制点数量；a、b分别为防撞梁长度的上、下限；Pi为B-Spline 曲线控制点坐标；Ni,p(u)为B-Spline曲线函数，其中i为样条曲线基函数的个数，p为样条曲线的次数，其表达式为：

基于上述公式由B-Spline 曲线描述的CFRP 防撞梁轮廓如图3所示。由图3可知，OA段由直线构成，AC段曲线是由点B控制的二次均匀B-Spline 曲线，BD段曲线是由点C控制的二次均匀B-Spline 曲线，CE段曲线是由D点控制的二次均匀B-Spline 曲线。B、D点的横坐标为设计变量中的曲线关键点横坐标。

图3 组合B-Spline曲线

3.2 CFRP防撞梁碰撞工况仿真

当前，防撞梁正面碰撞工况仿真模型有2种建模方式：一是将防撞梁进行固定，赋予刚性墙一定的速度和质量，使刚性墙撞击防撞梁；二是将刚性墙固定不动，使防撞梁以一定的速度进行撞击。由于第2种方案更加贴近实际的碰撞情况，故本文采用第2种方案进行仿真。

将防撞梁以4 km/h的速度撞击固定的刚性墙，防撞梁与刚性墙的接触形式为自动面面接触，将静摩擦因数设定为0.3，动摩擦因数设定为0.2时模型不会出现穿透现象[15]。

利用MATLAB 与Ansys 进行联合仿真，首先在MATLAB中输入防撞梁的结构参数和材料参数，采用二次B-Spline曲线完成对防撞梁轮廓曲线的描述，将轮廓曲线的坐标以txt格式输入Ansys进行有限元建模。由于Ansys库中没有所需的材料模型，而LS-DYNA的材料库较为齐全，所以需在LS-DYNA 中进行建模。建模完成后，在LS-PrePost中进行碰撞仿真，得到目标数据后再通过txt格式输出到MATLAB 中进行计算，最终得到碰撞力和比吸能数据。具体有限元碰撞仿真过程如图4所示。

图4 有限元碰撞仿真过程

4 试验设计方法与克里金代理模型

本节通过拉丁超立方试验法选取样本点，再通过上节介绍的CFRP 防撞梁的碰撞工况有限元仿真获取样本点的响应值，之后便可通过克里金（Kriging）代理模型来建立目标函数以及功能函数的代理模型。

4.1 试验设计

汽车碰撞问题为高度非线性问题，单次的碰撞仿真会消耗大量的时间，优化设计需要多次迭代并调用仿真数据。解决该问题的通常做法是选取适当的试验设计方法，在样本空间中生成足够的样本，通过这些样本得到近似模型，完成优化设计，而非碰撞仿真。试验设计的选择标准为：生成的样本足够多且均匀。最理想的试验设计是以较少的样本点得到较为精准的预测结果，多数学者在处理碰撞仿真等非线性问题时，常采用拉丁超立方法。

进行n次抽样，拉丁超立方法将m个随机变量分别等分为n个区间，样本空间即被分为n个m维区间。对于每个变量，可以保证n次随机抽样分别落在各小区间，样本点可等概率地分布到整个随机空间内。通过拉丁超立方法生成了120个数据，同时根据复合材料铺层的铺设原则，应避免将同一铺层角度的铺层连续放置，一般不超过3层，以避免出现应力集中和内部微裂纹，根据此原则对样本点进行删除。由于篇幅有限仅展示部分结果，如表2所示。

表2 拉丁超立方试验设计表

4.2 建立Kriging代理模型

在实际工程问题中，通常需要使用有限元的方式对工程进行分析，但会缺少显式的极限状态函数以及目标函数，造成可靠度计算的困难。为解决该问题，需引入代理模型法，常见的代理模型有响应面模型、Kriging 代理模型、人工神经网络模型等。其中Kriging近似模型有较好的全局拟合精度，且考虑了局部偏差，可以较为精准地拟合变量值与响应值之间的数学关系，使得优化设计的计算成本降低。半参数化的Kriging 模型不需要建立一个显式的数学模型，相较于参数化模型，如响应面模型，半参数化的Kriging 模型的应用更加灵活。

Kriging 近似模型由一个参数模型和一个非参数随机过程构成：

式中，F(K)为参数模型，形式为多项式回归方程；Z(K)为模型局部偏差的随机项，其均值为0，方差、协方差不为0。

5 RBDO分析

基于前文所得到的数学模型和隐式函数对CFRP防撞梁进行可靠性设计优化。传统的可靠性设计优化问题的本质为两层嵌套的双循环优化问题，分别为内层可靠性分析循环和外层优化循环。可靠性分析方法主要分为3 类：近似解析法、数值模拟法和间接代理法。近似解析法包括一阶可靠性法、二阶可靠性法等，但该类方法具有固有局限性，在进行可靠性分析时精度偏低，当随机变量为非正态、功能函数高度非线性时，无法得到可行解。数值模拟法的典型方法有蒙特卡罗模拟法，虽然该方法适用范围广、精度高，但在处理小失效概率问题时计算成本极高，无法解决大型工程结构中的可靠性分析问题。基于以上考虑，本文将通用生成函数法引入防撞梁碰撞可靠性分析，在随机变量非正态、功能函数高度非线性的情况下得到更加精准的结果。

5.1 基于UGF的可靠性分析

对于连续变量S，假设其累积分布函数与概率密度函数分别为FS(s)和fS(s),将该变量在其定义域（smin，smax）内均匀离散成m个点，而每个离散点Si对应的概率值Pi为：

式中，δ=(smax-smin)/m为离散步长。

通过式（6）计算可得离散数据集{(Si,Pi)|i=1,2,..,m},从而定义连续型随机变量S的UGF为：

式中，US(z)为随机变量S的UGF；z为离散后随机变量组成的多项式。

对于拥有n维连续型随机向量S的工程结构，需要在可靠性分析时根据式（6）、式（7）获得S各分量的UGF，记为：

式中，指数项siji为随机变量Si的第ji个状态值；Piji为其对应的概率值。

设G(S)为工程结构的功能函数，则G(S)≥0 时结构可靠，G(S)＜0 时结构失效。为了获得结构的总体UGF，需要对每个变量的UGF进行复合运算，表达式为：

式中，UG(z)为针对功能函数G(S)的结构UGF；⊗G为复合算子。

式（9）可进一步化简为：

式中，M为随机变量离散状态组合总数。

最后对UGF 的系数项求和，即可得到最终可靠度R：

式中，ψ()为条件求和算子；I()为示性函数，当Gi＜0时取0，否则取1。

5.2 UGF-直接映射方法

RBDO 问题是一个双层嵌套问题，针对RBDO 的求解效率问题，当前研究已提出了解耦算法，其思路是将RBDO的嵌套解耦成一系列确定性设计优化和可靠性分析组成的迭代过程。其代表方法是序列优化与可靠性评估（Sequential Optimization and Reliability Assessment，SORA）方法，通过引入偏移向量使最可能失效点（Most Probable Point，MPP）落在可行域内，由于引入了偏移向量，称为间接映射。采用UGF-直接映射方法[16]，可通过响应面模型来拟合设计变量与可靠度指标的函数关系式，避开MPP 且不需要偏移向量。由于SORA 法对初始偏移向量较为敏感，若不能确定初始偏移向量，则会导致优化结果不稳定，同时，基于MPP 来求解初始偏移向量存在原理性误差，不适用于非正态变量和非线性极限状态问题，故采用UGF-直接映射法可提高求解效率。

UGF-直接映射法的原理为通过构建响应面模型来拟合设计点和可靠度指标之间的函数关系，生成的函数定义为指标函数，将指标函数带入RBDO 的概率约束，将含有概率约束的不确定优化转化为确定性优化，转化后的UGF-直接映射法的可靠性优化设计数学模型可表示为

式中，(d,P)为指标函数；许用可靠度指标βT=2。

5.3 非均匀聚类技术

UGF可靠性分析法将随机变量离散化后，需将离散后的状态项进行状态组合，但在实际工程中随机变量较多，会产生过多的状态组合，使功能函数的评估次数随之增多，导致可靠性分析求解效率极低，此时可通过非均匀聚类技术解决此问题。

非均匀聚类算法对传统K-均值（K-Means）聚类算法的聚类中心更新策略进行了改进，提出了簇心概率加权方法。在可靠性分析中，极限状态空间附近的点对结构的失效或可靠性分析更加敏感，另外，概率密度大的点对可靠性分析的结果影响也较大，基于这两点对数据点进行加权，使得数据点在极限状态曲面附近的同时也要向概率密度大的区域集中，可以使需要的设计点数量减少且提高可靠性分析的精度。

簇心概率加权方法的具体操作为：为了考虑概率密度的影响，先对坐标点进行坐标概率加权，即在划分高维空间坐标点对应的簇时，将同一簇的坐标点乘以对应的取值概率，得到坐标向量后进行求和，所得的向量记作E1。将各坐标点的取值概率求和，所得的值为P；为了向极限状态空间附近靠近，则对坐标点的功能函数值加权，功能函数值加权系数j为：

然后将同一簇的坐标点乘以每个坐标点对应的功能函数值加权系数，将得到的值进行累加，所得的向量记作E2。将各功能函数值加权系数进行求和，记作g；最后对中心点进行概率加权更新：

式中，G为各坐标点对应的功能函数值；B为簇心概率加权系数；f为组合系数，通常取0.3。

具体流程如图5所示，以本文4个随机变量为例，将各变量离散成20个状态项。通过非均匀聚类算法的两步聚类，将160 000 个状态项缩减至8 000 个，在可靠性分析环节中，只需对全空间聚类所得到的8 000个状态项进行功能函数评估，相比较不使用非均匀聚类方法的160 000 次功能函数评估，效率大幅提升。同时，在进行高维聚类时，通常会用到大量的多重循环语句，导致代码运行过慢，针对这一问题，本文采用矢量化编程技术来避免大量循环，从而提高代码运行效率。

图5 非均匀聚类流程

6 CFRP防撞梁碰撞RBDO流程

在CFRP 防撞梁的RBDO 问题中，设计变量为控制点横坐标、防撞梁厚度以及复合材料铺层角度。其中控制点横坐标为连续型设计变量，厚度和铺层角度为离散型设计变量，且厚度和铺层角度要求为整数，针对上述的混合整数非线性优化问题，传统的优化算法无法得到结果，本文拟采用遗传算法进行优化。

6.1 优化流程

CFRP防撞梁碰撞结构可靠性优化设计具体步骤如下：

a.通过拉丁超立方设计方法对所有设计变量进行试验设计，生成样本点；

b.基于生成的样本点，同时结合MATLAB与Ansys的相互调用进行防撞梁的有限元碰撞仿真；

c.基于生成的样本点，通过MATLAB 与Ansys 计算相应的目标值；

d.基于计算出的目标值，通过Matlab 中的Kriging工具箱，即可构建有关目标函数及功能函数的Kriging代理模型，并验证其精度，若不满足要求则补充样本点并回到步骤c，若满足要求则进行下一步；

e.确定初始设计点D0，并通过UGF 法求解对应可靠度指标；

f.构建可靠度指标响应面。

g.使用遗传算法计算当前设计点Dk，采用非均匀聚类技术并结合UGF法计算对应的可靠度指标；

h.验证结果是否收敛，若不收敛，则将初始设计点和当前设计点合并，返回步骤f，若结果收敛，则输出最优结果。

优化流程如图6所示。

图6 优化流程

6.2 结果对比

将传统矩算法、蒙特卡罗法、UGF 法进行比较：算法1 为基于FORM 的双循环法；算法2 为基于FORM 的SORA法；算法3为UGF-直接映射法，可靠性分析由UGF法计算；算法4的可靠性分析由蒙特卡罗法计算；算法5为UGF-直接映射法的变形，优化流程中使用了非均匀聚类技术。算法1与算法2属于传统矩算法，算法4属于蒙特卡罗法，算法3和算法5属于UGF法。相对误差ε为：

式中，β为目标设计点的可靠度指标；βT为许用可靠度指标，取值为2。

不同算法优化结果对比如表3所示。

表3 不同算法优化结果对比

由表3可知，当极限状态函数高度非线性及随机变量非正态时，基于FORM 的双循环法、基于FORM 的SORA法无法收敛，UGF-直接映射法、MCS-直接映射法以及基于非均匀聚类算法的UGF-直接映射法的可靠度指标均满足许用可靠度指标要求。由于每次迭代过程中，均需对功能函数进行调用计算，所以可将功能函数调用次数作为衡量算法效率的指标。由此可得，基于非均匀聚类算法的UGF-直接映射法效率最高，普通UGF-直接映射法效率次之，而MCS-直接映射法的效率最低。

7 结束语

本文对CFRP 防撞梁进行了可靠性优化设计。对传统矩方法、SORA 法、UGF-直接映射法、MCS-直接映射法、结合非均匀聚类法的UGF-直接映射法进行了对比。结果表明，传统矩方法无法收敛；MCS-直接映射法求解精度最高，但计算成本过高；UGF-直接映射法可以稳定收敛，且求解精度较高，计算成本较小。结合非均匀聚类法的UGF-直接映射法的精度略低，但求解效率得到了极大提升，在应对更多随机变量的RBDO问题时更加适用。将初始设计点代入目标函数后可得初始目标值为81.22 J/kg，经过算法5 的优化后得到的目标为103.64 J/kg，增幅为21.6%，达到预期优化效果。关于样件试制与试验的问题，刘成龙[17]通过对复合材料样件进行拉伸试验，将计算试样的可靠度与理论仿真的结果进行对比，结果显示，理论结果和试验结果仅存在6.8%的误差，可间接验证UGF 的合理性和有效性。但关于碰撞吸能试验因条件不足，直接验证存在难度，需后续加以深度研究。