數據驅動的曲面構件形狀?拓撲協同優化方法

高天賀，田闊，黃蕾，張澍，李增聰

大連理工大學 工業裝備結構分析優化與CAE 軟件全國重點實驗室 工程力學系，大連 116024

曲面構件是組成航空航天結構的關鍵零部件［1-3］?；鸺械恼髡?，飛機機翼、尾噴管以及密封艙壁板、艙門等都是曲面構件。曲面構件在運載火箭中占干重比80%以上，在飛機中占干重比50%以上［4］。隨著我國航空航天事業的不斷發展，對航空航天結構的壽命及可靠性提出更高的要求［5］，所以研究曲面構件的輕量化設計方法具有重要意義。

近年來，拓撲優化方法憑借其不依賴指定的初始結構的優勢，已逐漸成為曲面構件概念設計階段的重要手段［6］。曲面構件的拓撲優化方法主要包括固體各向同性材料懲罰法（Solid Isotro?pic Material Penalized，SIMP）［7］、水 平 集 法（Level Set，LS）［8］、可變形孔洞法（Moving Mor?phable Void，MMV）［9］等。其中SIMP 法以單元密度為設計變量，具有收斂性好、靈敏度分析簡單等優點，但應用于復雜曲面構件時卻容易產生棋盤格或灰度單元等問題，通常需要引入一些過濾技術［10］或B 樣條參數化方法［11］。相比于SIMP法，LS 方法用一個高維度的水平集隱式表達拓撲結構邊界，產生的孔洞或筋條更加清晰。但是，LS 方法收斂性較差，利用該方法分析非線性問題的靈敏度較為困難。近年來，相關學者將LS 方法與等幾何分析相結合［12］，提高了拓撲邊界應力的計算精度，解決了應力約束拓撲優化問題應力奇異的挑戰。除上述方法外，張衛紅團隊［13］還提出一種幾何背景網格法，用于航天三維曲面薄壁殼的加筋設計。Feng 等［11］針對曲殼結構提出了一種有效的B 樣條參數化方法，與基于密度的方法相比，結果更清晰且無棋盤格現象。Zhou 等［14］參考蠕蟲的幾何形態，提出了基于仿生B-樣條偏移的加筋設計方法。Dong 等［15］基于自然界分支系統生長機理，提出了面向薄壁結構抗屈曲的最佳加筋布局自適應生長方法。Hu 等［16］對自適應增長方法進行了改進，針對加筋角度、位置和長度等參數進行了尺寸優化，提高了結構的固有頻率。綜上，目前學者們已針對曲面構件的拓撲優化方法開展了大量的研究，相比傳統依賴人工經驗的設計方法，拓撲優化往往更容易產生創新性設計構型。

然而，由于曲面構件的某些性能受曲面形狀的影響很大，當曲面構件形狀固定時，僅依靠拓撲優化方法難以保證獲得的材料分布是最優的［17］。目前，學者們針對曲面構件形狀-拓撲協同優化問題開展了一定研究。最早的形狀-拓撲協同優化方法，是在拓撲優化后，利用CAD 方法［18］或一些光滑的隱式曲線［19］對孔洞形狀進行局部優化。進一步，為將壁面形狀與拓撲優化協同考慮，相關學者提出了一種形狀與拓撲交替進行，直到優化問題收斂的分步迭代優化方法［20］。該方法在優化過程中會改變曲面形狀，為拓撲優化提供了更為靈活的設計域。然而，由于形狀和拓撲變量被單獨處理，形狀與拓撲變量的相互作用無法得到充分的考慮。近年來，學者們開始將SIMP 法［21-22］等拓撲優化方法與形狀設計變量的靈敏度分析方法相結合，通過推導靈敏度公式［23］，進行形狀與拓撲的耦合求解。然而對于復雜問題，靈敏度的推導較為困難，將耦合求解方法應用于實際工程還存在較大難度。如何在考慮形狀-拓撲相互作用的同時，避免復雜的設計變量靈敏度推導，使其能較好的應用于工程實際中的復雜問題，是形狀-拓撲協同優化的研究挑戰。

近年來，數據驅動方法受到國內外學者的廣泛關注，為解決形狀-拓撲的協同優化問題提供了新的途徑。目前已有許多新穎的數據驅動的形狀優化［24］和拓撲優化［25］方法，然而如何應用數據驅動思想，將形狀優化與拓撲優化協同考慮，還有待進一步研究。針對數據驅動的建模，相關學者先后提出了網格映射［26］和網格變形［27］，并成功應用于回轉曲面的拓撲優化模型重構［28］。該方法為參數化建立有限元網格模型，實現形狀-拓撲協同優化提供了可能。數據驅動方法可以利用一定數目的樣本，擬合具有較高非線性程度的代理模型，并利用代理模型代替真實計算，能有效提高優化求解的效率［29-31］。Meng 等［32］構 建 了Kriging 自 適 應 代 理模型，并開展了多學科協同優化設計，顯著提高了高度非線性系統協同優化的精度和效率。因此，基于數據驅動的網格變形建模及代理模型優化方法，提出一種數據驅動的曲面構件形狀-拓撲協同優化方法，旨在拓寬航空航天曲面構件設計空間，為航空航天曲面構件提供更為輕量化的創新構型。

本文主要包括3 部分，第1 部分介紹數據驅動的曲面構件形狀-拓撲協同優化方法；第2 部分開展簡支梁和航天器艙門算例研究，以驗證所提出方法的有效性；第3 部分為結論與展望。

1 數據驅動的形狀?拓撲協同優化方法

首先介紹形狀方程-網格變形驅動的形狀優化方法，該方法通過用數學方程中的少量參數描述控制點的位移，將形狀設計變量用少量的參數表示。然后，介紹考慮沖壓約束的SIMP 拓撲優化方法，針對每一個網格變形的設計域開展拓撲優化設計，獲得優化結果。最后，介紹形狀-拓撲協同優化的技術路線，包括離線、在線、更新3 個階段。

1.1 形狀方程?網格變形驅動的形狀優化方法

傳統的基于節點的形狀優化方法往往面臨著設計變量規模較大、網格更新困難和容易收斂到不滿足工藝要求的構型設計等問題［33］。因此，提出一種以形狀方程描述形狀設計變量、以網格變形為驅動的形狀優化方法。該方法能夠減少形狀設計變量個數且避免形狀改變后網格重新剖分，便于獲得滿足真實需求的設計構型。如圖1 所示，該方法主要包括以下步驟：

圖1 形狀方程-網格變形驅動的形狀優化方法Fig.1 Shape equation-mesh deformation-driven shape optimization method

步驟1：建立有限元模型，劃分形狀優化的可設計域和不可設計域。在剖分好的網格模型中均勻選取控制點，以保證網格變形后網格邊界的光滑性。

步驟2：建立描述壁面或截面形狀的數學方程，根據數學方程求解控制點網格變形后的位置坐標。

步驟3：基于徑向基函數神經網絡，擬合控制點變形前后的兩個坐標數據集合，獲得二者之間的映射關系。

步驟4：基于訓練好的映射關系，對原網格模型所有節點坐標進行預測，生成新的網格模型。為保證網格質量，對變形后的網格進行拉普拉斯網格光順化處理，獲得邊界平滑、網格質量較高的有限元網格模型［27］。

經過上述步驟，即可完成網格變形驅動的參數化建模過程。需要說明的是，上述方法將數學方程描述壁面形狀和網格變形改變壁面形狀相結合，僅需要提供數學方程的若干參數，便可自動的進行控制點的移動、映射關系的訓練，生成可用于有限元分析的網格模型。

1.2 固體各向同性材料懲罰拓撲優化方法

SIMP 法［34］作為一種經典的拓撲優化方法，具有收斂性好、靈敏度分析簡單等優點。因此，以SIMP 法為例，開展形狀-拓撲協同優化的研究。以結構的應變能最小為優化目標，體積為約束，優化問題為

式中：ρi為單元i的偽密度值；N為模型設計域的單元總數；為避免剛度陣奇異，δ取0.001；F為各單元的等效節點力；U為各單元節點位移；K為整體剛度矩陣；V0為結構初始體積；Vf為約束體分比。

曲面加筋是一種常見的結構形式，為獲得符合工藝約束的拓撲結果，往往需要約束沖壓方向每一列材料的偽密度相等，來施加沖壓約束［35］。以圖2 為例，假設沖壓方向共四層網格，則其沖壓約束表達式為

圖2 沖壓約束效果示意圖Fig.2 Schematic diagram of effect of stamping constraint

沖壓方向網格單元偽密度在優化過程中保持相等。

通過求解上述優化問題，獲得優化后的偽密度場，即可獲得結構的拓撲優化構型。

1.3 數據驅動的形狀?拓撲協同優化方法

一種數據驅動的形狀-拓撲協同優化方法如圖3 所示，所提出方法主要包括以下3 個階段：

圖3 數據驅動的形狀-拓撲協同優化方法流程圖Fig.3 Flow chart of data-driven shape-topology collab?orative optimization method

1） 離線階段。離線階段需要利用少量的樣本點建立全局高精度的代理模型。首先，基于拉丁超立方采樣方法［36］在設計空間進行采樣獲得形狀變量樣本（即形狀方程中的參數組合），使樣本點具有較好的空間覆蓋性，避免構建的代理模型局部坍塌［37-40］。然后，將形狀變量樣本代入形狀方程，利用網格變形技術參數化、自動化地獲得樣本點數據對應形狀的網格模型?；赟IMP法對獲得的網格模型進行拓撲優化，獲得拓撲優化后的響應結果（如應變能、最大位移、頻率等）。訓練徑向基函數（Radial Basis Function，RBF）代理模型［41］，獲得形狀設計變量（數學方程中的若干參數）與拓撲優化響應結果之間的映射關系。其中，代理模型的精度指標［1］為留一驗證的R2和RRMSE，其表達式分別為

式中：k為所有樣本點的個數；yi為第i個樣本點響應的真實值；y?i為第i個樣本點響應的預測值；yˉ為所有樣本響應的平均值。R2越接近1，RRMSE越接近0，代理模型的預測精度越高。

在實際應用中，可根據實際問題的難易程度及計算資源確定合適的訓練樣本數目。需要指出的是，當邊界和載荷條件發生變化時，代理模型需要重新訓練。為保證在線階段優化能夠穩定快速收斂，離線階段建立的代理模型需要保證具有較高的全局精度，如R2大于0.9。如果精度過低，可以通過自適應序列加點策略（如CVVoronoi 方 法［38］、EI 加 點 方 法［39］等）進 行 更 新 和逐步提升精度。

2） 在線階段?；谌謱炈惴▍f方差矩陣自適應進化策略（Covariance Matrix Adaptive Evolution Strategy，CMA-ES）［42］和RBF 代理模型進行優化。該階段基于離線階段訓練好的代理模型開展高效的優化，無需進行拓撲優化設計。由于代理模型根據輸入獲得輸出幾乎無耗時，所以可在優化過程中使用充分多的種群數目和迭代次數。CMA-ES 優化結束后，以CMAES 優化得到的最優結果為起點，利用序列二次規劃 算 法（Sequential Quadratic Programming，SQP）［43］對RBF 進 行 再 次 局 部 尋 優，確 保 真 正收斂。

3） 更新階段。為避免程序無限制的運行下去導致RBF 代理模型在局部密集加點，甚至產生矩陣奇異等問題。根據韓忠華等［44-45］以及Fang等［46］的研究成果，利用最大迭代次數限制更新次數。最大迭代次數可根據問題難易程度或代理模型精度設置。所提出方法利用形狀方程驅動的網格變形技術降低了設計變量維度和優化難度，并通過在離線階段構建具有較高精度的代理模型、在線階段通過CMA-ES 全局搜索結合SQP 局部搜索的模式，保證優化算法的穩定收斂。因此通常不需要設置過大的迭代次數就足夠保證優化收斂。在線階段結束后，判斷是否達到最大迭代次數，如果未達到最大迭代次數，則將在線階段得到的優化解真實計算后加入原數據集，并重新訓練代理模型以達到提高代理模型精度的目的。代理模型更新結束后，基于重新訓練的代理模型重新進行在線階段的優化，直至達到最大迭代次數。更新階段可通過加點更新提升代理模型的局部精度，自適應的保證優化結果的穩定，減少拉丁超立方采樣隨機性對優化結果的影響。

核心創新點為提出了一種數據驅動的形狀-拓撲協同優化框架，通過建立形狀變量與拓撲優化響應的代理模型實現形狀與拓撲的協同優化。其中，離線階段的代理模型的采樣方法以及在線計算的優化算法可以根據需要進行替換，并不影響整個優化框架。

2 算例研究

2.1 算例一：簡支梁算例

首先，基于一個簡支梁算例來驗證所提出方法的有效性。簡支梁的尺寸和載荷邊界如圖4所示，考慮載荷及邊界的對稱性，對梁進行二分之一簡化。取梁的幾何中心為坐標原點，x軸為梁的跨度方向。梁上端面承受f=0.1 MPa 均布載荷，采用六面體單元對半跨簡支梁進行網格劃分，單元數目為1 200 個。簡支梁的彈性模量和泊松比分別為E= 67 GPa 和υ= 0.3，材料密度為ρ= 2.700 g/cm3。

圖4 均布壓強載荷作用下的簡支梁示意圖Fig.4 Schematic diagram of a simply supported beam under sinusoidal pressure load

梁的形狀變化如圖5 所示，要求梁沿跨中截面保持左右對稱，半跨梁的外輪廓必須為直線。梁跨中截面可以上下移動，記移動量為a，沿y軸方向為正，梁跨中截面高度可以增加或減小，記改變量為2b，高度增加為正。在半跨梁的上下端面均勻選取控制點，控制點坐標可由a和b表示：

圖5 簡支梁形狀優化示意圖Fig.5 Schematic diagram of shape optimization for the simply supported beam

式中：坐標原點位于變形前梁跨中截面的幾何中心；yU和yL分別為梁上下端面控制點的縱坐標。

形狀-拓撲優化以控制形狀的兩個參數a和b為設計變量，以拓撲優化后應變能最小為目標，約束拓撲優化后模型質量小于1.4 g。拓撲優化以應變能最小為目標，約束體分比小于0.4。

離線階段中，利用拉丁超立方采樣方法在形狀設計空間內抽取15 個樣本點，并利用網格變形技術和拓撲優化方法進行建模和計算，獲得樣本點的應變能和質量數據，并利用RBF 代理模型擬合設計變量與應變能、質量的映射關系。其中應變能的預測精度R2為0.956 5，質量的預測精度R2大于0.999 9。使用RBF 的原因是其訓練收斂速度快，耗時少，且具有較強的逼近能力［47］。根據測試對比，利用高斯過程代理模型預測應變能的精度R2為0.949 9，深度神經網絡的預測精度R2為0.818 4。進行一次留一驗證耗時分別為：RBF 耗時0.001 114 s；高斯過程耗時2.007 636 s；深度神經網絡耗時2 216.64 s?？梢园l現RBF 的訓練速度和精度均較高。然后進入在線階段的優化，CMA-ES 最大迭代次數設置為100，CMA-ES 結束后進行SQP 的繼續尋優，最大迭代次數為20。將最終優化結果真實計算后進入更新階段并更新代理模型，代理模型更新次數設為5 次。優化迭代過程及優化目標（Objective，Obj）的迭代曲線如圖6 所示，隨著代理模型的更新迭代，形狀-拓撲協同優化的目標函數逐漸減小并最終收斂，代理模型優化解的形狀和拓撲結果也趨于收斂，最終獲得優化后的形狀與孔洞構型。

圖6 簡支梁形狀-拓撲優化目標函數迭代過程Fig.6 Iterative process of objective function in shapetopology optimization of simply supported beam

進一步，將所提出方法與初始設計（僅拓撲優化）的優化結果進行對比。為保證公平對比在有限計算資源下不同優化方法的尋優能力，本文將計算時間作為控制變量，通過控制樣本點數目一致保證各方法計算時間一致。通過對比相同計算耗時下的優化結果，驗證所提出方法在有限計算資源下具有較高的尋優能力。對比結果如表1 所示，優化得到的形狀及拓撲結果如圖7 及圖8 所示，其中所提出方法（形狀-拓撲協同優化設計Ⅰ）優化得到的應變能為0.019 9 mJ（1 mJ=10?3J），與初始設計相比應變能降低20.08%，說明設計域形狀的改變為結構設計帶來更大的優化空間。此外，也與形狀-拓撲協同優化設計Ⅱ（19 代理模型樣本+1 精細校核）和Ⅲ（不構建代理模型，直接CMA-ES 優化）進行了對比，結果顯示所提出方法在相同計算資源下可獲得力學性能更優的設計。其中精細校核與優化加點的區別在于，優化加點需要將真實計算后的優化解加入訓練集并更新代理模型，而精細校核只需要真實計算優化解，不需要更新代理模型。由于達到相同的尋優效果所提出方法用時最少，也驗證了所提出方法的計算效率較高。此外，本節也嘗試定量對比無代理模型方法（形狀-拓撲協同優化設計Ⅴ）與所提出方法的計算效率。然而，無代理模型方法迭代100次（所提出方法耗時的5倍）仍與形狀-拓撲優化Ⅰ的優化結果有9%的差距。說明所提出方法的優化效率是無代理模型的5 倍以上。形狀設計變量上下限優化結果如圖7 及表1 所示。此外，還測試了設計域增大后相同材料用量約束的單一拓撲優化，其中擴大后設計域跨中寬度為120 mm，足夠包絡所提出方法設計空間的所有形狀。優化結果如圖7 （c）所示，優化構型不清晰，優化后應變能為0.040 7 mJ，相比所提出方法優化結果增大了1 倍以上。說明拓撲優化設計域并不是越大越好，直接擴大拓撲優化設計域并不總會得到與形狀-拓撲協同優化類似的結果，進一步驗證了所提出方法的必要性。

表1 所提出方法與其他方法簡支梁優化結果對比Table 1 Comparison between results of simply supported beam from proposed method and other methods

圖7 簡支梁上下限形狀及擴大設計域后拓撲結果對比Fig.7 Comparison of topology optimization results be?tween upper and lower bounds shapes and ex?tended design domain of simply supported beams

圖8 所提出方法與其他方法優化后簡支梁形狀和拓撲結果Fig.8 Shape and topology results of simply supported beam optimized by the proposed method and other methods

最后，本節也嘗試測試不同采樣方法對優化結果及空間覆蓋性的影響?？臻g覆蓋性利用代理模型在驗證集上的精度表示。因此，首先利用全因子（Full Factorial， FF）試驗設計［48］在設計空間均勻抽取100 個樣本點作為驗證集。然后，利用FF法均勻抽取15 個樣本點進行代理模型構建，最終優化結果為形狀-拓撲優化Ⅳ。經驗證，FF 法樣本點構建的代理模型在驗證集上的R2為0.974 6，LHS 的R2為0.993 9，說明FF 抽樣法的空間覆蓋度較差。在該問題中，兩種采樣方法的優化結果應變能誤差為1.5%，說明算法較為穩定，改變采樣方法，并不會對最終優化結果造成較大的影響。

2.2 算例二：航天器艙門算例

為進一步驗證所提出方法的有效性，本節針對航天器艙門開展算例研究。艙門的簡化模型及截面尺寸如圖9 所示，艙門門體安裝在門框上，艙門半徑為560 mm，包括蒙皮在內的設計域厚度為30 mm，門體中心半徑100 mm 的圓面在優化過程中保持水平。門體外表面承受p= 0.09 MPa 均布壓強作用，門框底部設置固支邊界。根據對稱性取四分之一模型進行簡化，對簡化后模型截面設置對稱邊界條件。采用六面體單元對艙門結構進行網格劃分，單元數目為57 780 個。艙門結構材料屬性如表2 所示。

表2 艙門結構材料屬性Table 2 Material properties of cabin door

圖9 艙門載荷邊界及形狀變化示意圖Fig.9 Schematic diagram of cabin door load boundary and shape change

艙門門體截面控制點y坐標可用超橢圓方程表示，根據模型的尺寸，可以構造出門體截面的數學方程：式中：c和n為艙門算例的兩個形狀設計變量；c為超橢圓的短半軸長度，可以表示艙門截面邊界的內凹程度，c越大，艙門截面邊界內凹程度越大，當c接近0 時，艙門截面邊界形狀趨于水平直線；n為超橢圓的次數，n為1 時艙門截面邊界為直線，n為2 時為標準橢圓。通過改變c和n，可實現艙門截面形狀的改變。

形狀-拓撲優化通過優化形狀設計變量，使拓撲優化后的應變能最小。拓撲優化以應變能最小為目標，以優化后的體積小于等于2 022.47 cm3（初始設計形狀設計域體積的30%）為約束，保證優化方法最終構型體積相同。為保證拓撲優化結果滿足工藝制造要求，對拓撲優化問題施加沖壓約束，沖壓方向為艙門門體的面外方向。此外，施加尺寸約束，最小尺寸為20 mm，最大尺寸為60 mm，最小間隙設為200 mm。

離線階段中，初始樣本點個數為30 個，RBF代 理 模 型 預 測 精 度R2為0.931 0，RRMSE 為0.674 3。在線優化CMA-ES 設置最大迭代次數為100，SQP 最大迭代次數為20。代理模型更新次數設置為15 次，優化迭代過程及優化目標（Obj）的迭代曲線如圖10 所示，最終可獲得清晰的設計構型。

圖10 艙門形狀-拓撲協同優化目標函數迭代過程Fig.10 Iterative process of objective function in shapetopology optimization of cabin door

形狀設計變量上下限優化結果如圖11 及表3 所示，上下限形狀優化結果分別為15 040 mJ和14 038 mJ。

表3 所提出方法與其他方法艙門優化結果對比Table 3 Comparison between results of cabin door from proposed method and other methods

圖11 艙門上下限形狀及擴大設計域后拓撲結果對比Fig.11 Comparison of topology optimization results be?tween upper and lower bounds shapes and ex?tended design domain of cabin doors

補充設計域足夠包絡住形狀設計空間的單一拓撲優化結果，如圖11（c）所示。最終優化結果如表3 所示，所提出方法優化結果（形狀-拓撲協同優化設計Ⅰ）為7 437.38 mJ，對比相同質量下初始設計（僅拓撲優化），應變能結果降低了37.93%，驗證了形狀-拓撲協同優化的有效性。擴大后的設計域上端面水平，下端面斜率為1，門體中心厚度為120 mm，優化后應變能為29 981.2 mJ。由于只進行拓撲優化，忽略了蒙皮形狀對拓撲優化結果的影響，該優化結果相比所提出方法優化結果增大了4倍以上。同時，將形狀-拓撲協同優化設計Ⅰ與形狀-拓撲協同優化設計Ⅱ（無代理模型更新）、Ⅲ（不構建代理模型，直接優化）進行對比，結果顯示所提出方法在相同計算資源下可獲得力學性能更優的設計，與算例一結論完全相同。由圖12 可知，所提出方法優化結果形狀規則且筋條清晰，由于沖壓約束在拓撲優化中的使用，使拓撲優化產生的筋條與蒙皮保持垂直，符合工藝要求。

圖12 所提出方法與其他方法艙門優化后形狀和拓撲結果Fig.12 Shape and topology results of cabin door optimized by the proposed method and other methods

3 結 論

通過網格變形技術和代理模型方法，實現了壁面形狀和拓撲的協同優化。最終形成了離線-在線-更新三階段的優化路線。面向工程中簡支梁和航天器艙門結構開展了算例研究，結果表明，相比于固定形狀的拓撲優化結果，形狀-拓撲協同優化的應變能可分別降低20.08%和37.93%，驗證了所提出方法的有效性。此外，還將所提出方法與無代理模型的形狀-拓撲優化、無代理模型更新階段的形狀-拓撲優化進行了對比，結果表明，在相同計算資源下，相比無代理模型方法，所提出方法優化的應變能分別減小18.11%和4.87%，驗證了代理模型可提高相同計算資源下的尋優能力；相比無代理模型更新階段的方法，所提出方法的應變能分別減小3.86%和1.09%，驗證了代理模型更新階段可保證最終優化收斂。

在后續研究中，將探索設計域內包含不可變結構功能性孔洞的網格變形方法，實現更為智能的參數化建模。此外，還將研究代理模型預測精度對尋優能力的影響，嘗試利用代理模型精度自適應平衡采樣數目與代理模型更新次數，并擴展到更為復雜的載荷工況和結構形式，進一步提高方法對復雜問題的適用性。