基于通用生成函數法的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 的合理性和有效性。但關于碰撞吸能試驗因條件不足，直接驗證存在難度，需后續加以深度研究。