基于參數化水平集方法的水下耐壓結構的拓撲優化

蔣垣騰，趙敏

（1.上海交通大學海洋工程國家重點實驗室,上海 200240；2.高新船舶與深海開發裝備協同創新中心,上海 200240）

0 引 言

隨著潛水器的工作環境向著更大的深度發展，其耐壓結構起到極大的作用。潛水器的耐壓結構用來裝載電子元件以及檢測設備，為確保這些設備不會因為受到海水的壓力與腐蝕而損壞，其耐壓結構需要有足夠的強度和穩定性。但是，其結構的質量增加同樣造成總體性能的下降，因此，潛水器的耐壓結構優化的目的是在保證結構強度與穩定性的條件下使結構的質量最小，并且使其結構的浮力與重力之比最大。

在以往的耐壓結構設計中，常見的潛水器耐壓結構有球形、橢球形、圓筒形等。采取的優化手段往往是通過尺寸優化或者簡單的形狀優化，尺寸優化方法主要是以結構的厚度、半徑以及加強筋等的尺寸大小為設計變量進行的以輕量化等為目標函數的優化設計，形狀優化則是在給定的形狀下，以與設計結構物理參量無關的形狀參數作為設計變量進行的結構優化設計，這類優化方式得到的結構在該形狀下可以滿足最優化條件。因此，究竟如何選取潛水器耐壓結構的最優拓撲形式，即材料空間分布的最佳拓撲形式，有必要在常規的設計經驗基礎上，采用一種新的設計方法進行深入探究。目前，越來越多的設計者將目光投向拓撲優化，拓撲優化方法是一種根據給定的負載情況、約束條件和性能指標，在給定的區域內對材料分布進行優化的數學方法。在拓撲優化方法下，以結構的形狀作為設計變量，通過這種優化手段，尋求得到一個滿足設計要求的最佳拓撲形狀。這一方法很大程度上將打破工程設計人員的思維定式，在工程實踐中具有一定的指導意義。

在水下耐壓結構的優化設計中，施加在結構上的載荷與傳統拓撲優化設計中的固定載荷不同，該類靜水壓力載荷隨著結構輪廓的變化而發生改變，因此，此類載荷又被稱作設計相關載荷[1]。關于在設計相關載荷作用下的拓撲優化問題，其關鍵問題在于受力邊界的搜索。

在以往的研究中，在受力邊界的搜索方法上，國內外許多研究者提出了多種解決方案。2000 年，由Hammar和Olhoff[2]開創了設計相關壓力載荷拓撲優化的先河，采用基于SIMP模型的變密度方法，在固體各向同性材料的灰度密度分布范圍內利用貝塞爾樣條曲線代表加載邊界，樣條曲線的控制點是由節點相對密度插值得到的等值點；2007 年，Sigmund 和Clausen[3]提出了一種混合位移壓力有限元公式，其中使用重疊壓力變量將空白單元定義為靜水不可壓縮流體；2008年，Zhang 等[4]提出了新的邊界搜索算法，通過預設搜索的起點與終點，并規劃壓力的加載路線，直接將壓力加載于結構單元邊界上，并避免了載荷的靈敏度分析；2009年，Zheng 等[5]引入了偽等勢函數來確定結構的邊界，采用伴隨法進行靈敏度分析；2012年，Lee 等[6]同樣直接連接單元邊界上的等密度值點作為壓力作用邊界，但是在設計模型中引入初始空白區域，壓力則是由空白區域加載到設計區域邊界上。從而，無需預先定義壓力作用曲線的端點，這種方法的優點在內部受壓和多壓力艙的優化設計中得到了充分的體現。并且，利用高斯積分的方法將等效節點力表達為節點密度值和密度閾值的表達式，并提出了載荷靈敏度計算的解析方法。

也有研究者針對水下耐壓結構設計這一工程問題，給出了解決方法。2006 年，茍鵬等[7]提出一種以柔度最小化為目標的優化算法OC-LE 法，用以解決在設計相關載荷作用下的結構優化問題，并給出外部邊界受靜水壓力作用下的結構設計算例，證明該方法在水下耐壓結構設計中的可行性；2015年，王存福等[8]提出一種基于圖像分割技術的加載面搜索算法，利用Li 等[9]提出的DRLSE 模型實現圖像邊界探測，為了避免載荷靈敏度分析，運用了新的優化算法模型，即認為加載面在連續的一定迭代步內保持不變，并且，著重討論了水下耐壓結構的設計問題；2017年，戴揚等[10]討論了采用不同復合材料的水下耐壓結構設計問題，從復合材料的鋪層方式以及角度的變化等因素，探究了復合材料在水下耐壓結構設計的前景以及應用意義。

上文提到的文獻都是采用密度法對在設計相關載荷下的結構拓撲優化設計或水下耐壓結構設計進行的研究，密度法雖然具有較強的拓撲表現能力以及成熟的優化算法，但是由于密度法得到的結構邊界過于模糊，導致在計算過程中有大量時間花費在邊界搜索上，因此降低了計算效率。但是，水平集方法在優化過程中將時刻提供輪廓清晰的結構邊界，這在處理表面發生多物理場相互作用的問題時具有優勢，并縮短計算時間，大大提高計算效率，在工程實踐中，這將縮短設計周期，具有一定的應用價值。

由于上述優勢所在，水平集方法吸引了許多學者進行研究。2005年，趙康等[11]提出一種利用水平集演化技術求解設計相關載荷作用下結構拓撲問題的數值方法，由于結構的邊界可以用零水平集加以描述，文中利用適當的數學變換，方便地處理施加在結構上的設計相關載荷，避免了以往算法中繁復的邊界提取工作以及為了處理設計相關載荷所采取的特殊技巧；2015 年，Xia 等[12]分別利用兩個水平集函數及其對應的水平集，來確定結構的壓力邊界以及進行水平集方法的拓撲優化，并且，為了保證優化結構的最優性，文中采用了一種周長正則化處理的方法；2019年，Picelli等[13]采用一種偽單元材料插值的固定有限網格，文中邊界搜索策略是采用由拉普拉斯方程控制的流體場來計算靜水壓力場加載線彈性結構，并將固定網格的偽單元材料法與等效節點力相結合。

在以往的研究中，學者普遍采用傳統水平集方法應用于設計相關載荷下的拓撲優化設計中，但是，傳統水平集方法在處理該問題上存在一些缺陷：第一，傳統水平集方法的偏微分方程的解法主要采用有限差分法，有限差分方法不適應不規則網格；第二，由于有限差分法對迭代步長的限制，因而導致尋優過程過長；第三，在結構的輪廓演化過程中，反復的重新初始化操作影響了計算效率。此外，在不考慮拓撲導數的條件下，傳統水平集方法的優化結果依賴于初始結構輪廓選取，即傳統水平集方法無法自行為設計結構開孔。為了提高計算效率，并提高計算結果的可信性，本文采用基于徑向基函數的參數化水平集方法。參數化水平集方法由Wang 等[14]在2006 年首次提出；2018 年，Wei 等[15]提供了以柔度最小為目標的88行參數化水平集方法下的MATLAB 程序，該方法繼承了傳統水平集方法的優點：如結構輪廓光滑且清晰，無需重新劃分網格等，將該方法應用在設計相關載荷作用下的問題可以解決傳統水平集方法中在該問題上諸如需要重新初始化等的一些弊端，最為重要的是這一方法將提高耐壓結構優化設計的效率。

本文將采用參數化水平集方法來解決水下耐壓結構優化設計問題，研究在體積約束下不同邊界約束條件下的以水下耐壓結構柔度最小化為目標的拓撲優化設計。由于水平集方法中結構的邊界信息隱含在水平集函數中，因此，本文提出一種邊界掃描的方法，簡便地提取壓力邊界信息，同時提出與參數化水平集方法匹配的拉格朗日算子更新策略，實現水下耐壓結構的拓撲優化設計，并通過四個經典水下耐壓結構設計的數值算例以證明本文方法的有效性。

1 參數化水平集方法

參數化水平集方法是指采用徑向基函數插值得到的水平集函數的方法，其中，徑向基函數的插值精度和插值過程中的計算效率與所選取的徑向基函數的種類有關，徑向基函數一般按種類分為全局徑向基函數(GSRBF)和緊支徑向基函數(CSRBF)[15]。

本文將采用的徑向基函數為MQ（Multiquadric）樣條曲線，這種徑向基函數屬于全局徑向基函數，這種曲線用公式表示為

式中，D為設計區域，xi為設計區域的一個樣本點，c為一個常數。一般而言，采用MQ 樣條曲線插值的精度由樣本點的選取個數與c的大小決定，本文中c取10-4。在水平集方法中，水平集函數正是一個相對結構輪廓而言的大型高維度標量場，隨時間推演而變化，并且沒有具體的函數表達形式。因此，將一組徑向基函數分布在對應的n個樣本點處，對水平集函數進行擬合，表示為

采用徑向基函數對水平集函數進行插值擬合，將Hamilton-Jacobi 方程從偏微分方程簡化成為一個常微分方程，將對計算產生積極影響，主要表現為以下幾點：

首先，采用徑向基函數插值后得到的水平集函數保證了其光滑性和連續性要求，并在其零水平集的附近能保持符號距離函數的特性，因此這種方法不需要在優化過程中反復地重新初始化，大幅減輕計算壓力并提高計算效率。

其次，步長不再受限于CFL 條件，步長的選取相較于傳統水平集方法更大，并得到合理的結構。這將會加快目標函數向最優解收斂的時間，并減少得到最優解所需的迭代步數。同樣，這一因素將會大幅提高計算效率。

此外，在傳統水平集方法中，初始設計結構的選取對最終設計結構的影響較大，而在基于徑向基函數插值的參數化水平集方法中，優化結構對初始設計結構的選取依賴較小，即在這種方法下，結構在優化過程中可以自行生成孔洞，這使得優化結果更具有普遍性意義。

2 問題描述

首先，本文所討論的問題為考慮靜水壓力載荷下的結構拓撲優化，為了簡化所涉及的優化問題，本文在二維空間內討論，并且討論的問題均認為處于平面應力狀態下。在設計相關載荷作用下的優化問題中，結構輪廓邊界可以劃分為三部分：

式中：ΓN為Neumann 邊界，即壓力所施加的邊界；ΓD為Dirichlet 邊界，即固定位移邊界；Γ0為自由邊界，在Γ0邊界上位移不固定且沒有施加外力。在優化過程中，這三種邊界將隨著結構輪廓的演變而變化。

結構的力平衡條件以及邊界約束條件如式（7）所示：

式中，u為結構的位移場，f為結構內部的體積力，為了將問題簡化，本文將不考慮體積力，即令f=0，p0（p0≤0）為Neumann邊界上的壓力。

本文目標函數為結構的柔度，結構的柔度等于外力在結構作用點上作的功，在設計相關載荷作用下，目標函數可以寫作

以上是本文優化問題的數學模型，設計變量為結構的形狀。由于優化算法要基于梯度的下降方向，因此，本文引入形狀導數的概念，關于形狀導數的具體定義與定理，已有多篇文獻[12,16]給出，本文將不再贅述。

3 靈敏度分析

本文將采用形狀導數的理論得到Hamilton-Jacobi方程的速度場，對于優化問題的目標函數，存在兩種約束：結構平衡方程以及邊界條件和設計約束，即體積約束。所以，下面將逐一求解在兩種約束下的速度場，即先求解在結構平衡方程以及邊界條件約束下的速度場，再在此基礎上，求解在體積約束下的速度場。

根據平衡方程，并通過高斯散度定理，可以得到平衡方程的弱形式為

4 體積約束處理

體積約束通過確定上文所提到的拉格朗日算子來實現，此前，有多篇論文討論過拉格朗日乘子的選取問題，例如：采用固定拉格朗日乘子[16]，通過每一迭代步中得到的速度場求得拉格朗日乘子[11]以及采用二分法[14]求解等，這些方法得到的拉格朗日算子應用于基于徑向基函數插值的水平集方法表現相對保守，收斂過程較長；針對載荷數量較大時，如設計相關載荷條件下，不能保證使目標函數在收斂過程中過渡平緩。因此，本文采用以下的拉格朗日算子的更新方法：

5 邊界搜索方法以及等效靜水壓力載荷

為了處理設計相關載荷，需要確定載荷作用的邊界，本文的搜索策略為：沿外載荷作用方向逐行或逐列進行定向掃描，在掃描方向上依次對每一個節點在與掃描方向垂直的方向進行側向掃描，直到節點的?≥0 后停止掃描，并標記為主標記點。之后再把與其鄰近的單元節點標記為副標記點。只有在主副標記點上，靜水壓力系數p0(x)不為0。壓力邊界搜索實現的基本思想如圖1所示。

圖1 靜水壓力邊界搜索的實現方式Fig.1 Realization of hydrostatic pressure boundary search

本文提出的壓力邊界搜索方法可以有效解決有結構遮蔽區域的問題，得到的壓力邊界將更加準確。

為了準確且方便地將設計相關載荷施加在結構的壓力邊界上，將討論設計相關載荷的計算方法，通過有限元分析中等效節點力的方式以及線積分轉化為體積分的公式，可以得到

6 數值算例

本文的算法總結如下：

（1）初始化水平集函數，定義結構輪廓以及結構內外區域；

（2）輸入一組徑向基函數并對水平集函數進行擬合，從而得到這組徑向基函數的插值系數；

（3）輸入材料參數以及結構的邊界條件；

（4）進入迭代步，檢測壓力邊界并計算壓力大小，形成有限元分析中的力向量F；

（5）組裝當前結構的總剛度矩陣K，通過有限元分析得到位移場U；

（6）計算當前結構的形狀導數以及拉格朗日算子，得到Hamilton-Jacobi方程的速度場；

（7）演化水平集函數，更新徑向基函數前的插值系數；

（8）判斷收斂性，比較目標函數在相鄰5 個迭代步中的相對誤差是否在合理范圍內，若不滿足收斂性條件，則返回到（4）。

本文通過具體數值算例來表明上述方法處理設計相關載荷作用下結構拓撲優化問題的有效性，數值算例為平面應力拓撲優化問題，楊氏模量的取值為1，泊松比的取值為0.3，編程語言為Python。

6.1 受外壓拱形結構設計

本文將先通過水平集方法中的經典算例，初步證明方法的有效性。受外壓的拱形耐壓結構模型的設計區域、位移邊界條件以及初始壓力載荷條件如圖2 所示，設計區域為（6 m×4 m）的矩形。結構在左右兩端部附近受到固定約束，其約束寬度為0.3 m，設計材料體積為不超過設計區域的30%。設計區域的網格劃分為60×40，受靜水壓力的大小為1。

圖2 受外壓拱形結構設計模型Fig.2 Design model of arched structure subjected to external pressure

受外壓的拱形耐壓結構設計的最終結果如圖3 所示，最終柔度為81.56。這一算例作為經典算例出現在中外多篇水平集方法的研究論文[11-13]中，本文得到的結構形式與以往論文相類似，算例的計算模型與趙康等[11]文中相似，趙康等對比了在初始設計域設置空洞與否以及在不同網格尺寸下的結果，得出剛度最大的結構為初始設置空洞的結果，其柔度為92.32。拓撲優化設計作為概念設計階段的主要工具，其主要目的是為工程設計提供方案，因此在概念設計中得到的結果將極大程度影響結構的性能，這將會更好地為工程實踐服務。

圖3 本文最終結果與趙康等[11]得到的結果對比Fig.3 Comparison between results in this paper and those obtained by Zhao et al[11]

6.2 靜水壓力下的結構設計

耐壓結構在水下工程中起到至關重要的作用，尤其是在潛水器的設計過程中。在實際工程中，耐壓結構的設計需要考慮諸多方面的因素，為了方便研究，本文只在靜水壓力下分為兩種情況進行討論，在固定約束下受靜水壓力的耐壓結構設計和自由約束下受靜水壓力的耐壓結構設計，這兩種情況在工程中具有應用意義，可應用于潛水器的封頭以及平行中體的形狀設計中。

6.2.1 在固定約束下受靜水壓力的結構設計

在固定約束下受靜水壓力的耐壓結構設計的設計區域、位移邊界條件以及初始壓力載荷條件如圖4 所示，設計區域為一個正方形（4 m×4 m），其固定約束點位于正方形中線距中心點1.5 m 處，設計約束為體積約束，最終結構的體積不超過設計區域的20%。且為了計算的簡便性，本算例取用1/4 模型，設計區域的網格劃分為40×40，受靜水壓力的大小為1。

圖4 在固定約束下受靜水壓力的結構設計模型和1/4模型Fig.4 Design model and 1/4 model of underwater pressure structure under fixed constraints

1/4 模型的結構優化結果如圖5 所示，最終結構的柔度為16.41，目前本算例沒有出現在以往的水平集方法的研究文章中。因此，對比在密度法下得到的結構，本文算例結果與Zheng 等[5]、Li 等[17]得到的結果類似。Zheng等[5]在文中沒有給出具體的柔度值，因此，將本文的柔度與Li等[17]得到的結果相比較，在其論文中，楊氏模量取100，體積約束取設計區域的30%。在算例中提供了目標函數柔度在優化過程中的收斂曲線，再經過換算，在本文體積約束的條件下，本文結果的柔度依然小于Li 等[17]得到的結果（見圖6）。

圖5 在固定約束下受靜水壓力的結構設計的中間結果以及最終結果Fig.5 Intermediate results and final results of underwater pressure structure design subjected to hydrostatic pressure under fixed constraints

圖6 本文結果與Zheng[5]等得到的結果對比Fig.6 Comparison between results in this paper and those obtained by Zheng et al[5]

6.2.2 自由約束下受靜水壓力下的結構設計

自由約束下受靜水壓力的耐壓結構設計的設計區域、位移邊界條件以及初始壓力載荷條件如圖7所示，設計區域同樣為一個正方形（4 m×4 m），本算例中在正方形的兩條對稱軸上設計約束，限制對稱軸上的區域只能在中心線方向移動，約束寬度為0.9 m，設計約束為體積約束，最終結構的體積不超過設計區域的10%。同樣，為了計算的簡便性，本算例將取用1/4模型，設計區域的網格劃分為40×40，受靜水壓力的大小為1。

圖7 自由約束下受靜水壓力的結構設計模型與1/4模型Fig.7 Design model and 1/4 model of underwater pressure structure under free constraints

1/4 模型的結構優化結果如圖8 所示，最終結構的柔度為13.13，目前這一算例也沒有出現在以往水平集方法的研究文章中。因此，對比在密度法下得到的結構，本文算例結果與王存福等[8]的論文中得到的結果完全相似，在王存福等人的論文中，取材料的楊氏模量為210 GPa，泊松比為0.3，靜水壓力大小取2 MPa，體積約束與本文相同，王存福等文中得到的結構最終柔度為250.15 J，無量綱化后的最終柔度值為13.15，本文結構的柔度略小于王存福等文中的柔度，得到的結構趨近于圓環狀（見圖9），這一結果與實際設計中采用的形狀相同。

圖8 自由約束下受靜水壓力的結構設計的中間結果以及最終結果Fig.8 Intermediate results and final results of underwater pressure structure design subjected to hydrostatic pressure under free constraints

圖9 本文結果與王存福等[8]得到的結果對比Fig.9 Comparison between results in this paper and those obtained by Wang et al[8]

6.3 空心水下耐壓結構優化設計

在潛水器耐壓結構的設計中，一般需要提供一定的儲物空間，因此空心耐壓結構的優化設計更加具有工程意義。與靜水壓力下的耐壓結構設計不同，上述討論給出了潛水器橫剖面的優化形狀，在本節將討論潛水器的縱剖面的優化形狀。空心耐壓結構的設計模型中，在其矩形設計區域（20 m×8 m）中預留一塊矩形（16 m×4 m）空白區域，空心耐壓結構設計模型的設計區域、位移邊界條件以及初始壓力載荷條件如圖10 所示，作用在設計模型上的約束在圖10 中進行了標注，約束類型分為完全固定和縱向固定兩種情況。設計材料體積為不超過設計區域的30%，為了降低計算壓力，本算例將選用1/4模型，設計區域的網格劃分為100×40，受靜水壓力的大小為1（見圖11）。

圖10 空心耐壓結構的設計模型Fig.10 Design model of hollow underwater pressure structure

圖11 空心耐壓結構的1/4模型Fig.11 1/4 model of hollow underwater pressure structure

1/4 模型的結構優化結果如圖12 中所示，其中圖12（a）是約束處為完全固定約束的，最終結構的柔度為80.32，圖12（b）是約束處為縱向固定約束的，最終結構的柔度為94.88。目前這一算例也沒有出現在以往的水平集方法的研究文章中。因此，與用密度法下得到的結構進行對比，本文結果與以往論文中得到的多球交接結構在形式上類似，用Li等[17]的算例對比本文，其楊氏模量以及受力大小分別為7.0×1010Pa、10 000 N/m，其體積約束取設計區域的30%，得到的結構柔度值為0.5877 和0.6771,兩個算例換算為本文單位結果分別為83.96 和96.72，本文結果與Li 等[17]得到的結構柔度值相比較小，因此，證明了本文結果的有效性。圖13 中比較了本文結果與Li 等[17]的結果，圖14 為MIT 水下潛水設計團隊[17]提出的多球交接耐壓結構模型，本文結果的結構形狀與其類似，這一概念設計的提出預示著這種結構形式具有工程參考價值，且本文方法將會引導工程設計的實踐。

圖12 空心耐壓結構設計的1/4模型的最終結果Fig.12 Final results of 1/4 model of hollow underwater pressure structure

圖13 本文最終結果與Li等[17]得到的結果對比Fig 13 Comparison of the results in this paper and those obtained by Li et al[17]

圖14 MIT提出的多球交接耐壓結構模型[18]Fig.14 Multiple intersecting spherical pressure hull model proposed by MIT[18]

7 結 語

本文研究了考慮靜水壓力下水下耐壓結構的拓撲優化設計問題，該問題是一個典型的設計相關載荷作用下的拓撲優化問題。由于傳統水平集方法在該問題中的低效性，因而選取基于徑向基函數插值的水平集方法，并利用參數化水平集方法可以得到結構輪廓清晰邊界的特性，提出了一種可掃描有遮蔽區域的邊界搜索方法。該方法利用了參數化水平集方法中邊界信息隱含于水平集函數的特性，通過逐列或逐行掃描的方式得到了結構受靜水壓力的邊界，可以有效并完整地得到結構的受力邊界，在此基礎上可以有效解決靜水壓力下水下耐壓結構的拓撲優化設計問題。

本文重點討論了水下耐壓結構設計以及在殼體內部局部加強對耐壓結構剛度的影響。水下耐壓結構拓撲優化設計是以柔度最小化為目標，以結構體積為設計約束，結合改進的拉格朗日算子，在平穩的優化過程中，得到了在不同邊界約束條件下水下耐壓結構的優化設計結果，文中的邊界約束條件是對諸如加強筋以及球殼交接處等的簡化。優化結果表明，該方法可以處理諸如多球交接耐壓結構等存在較為復雜邊界下的水下耐壓結構設計問題，經過對比可知，該方法與以往研究得到的結構形式相吻合，且得到的結構相較以往論文得到的結構具有較高的剛性。在未來新型水下耐壓結構的探索過程中，該方法可以為工程設計人員的設計工作提供參考依據，具有重要的工程應用價值。

綜上，本文提出的邊界搜索方法在二維平面應力問題中簡易且有效，為了進一步與工程實際情況接軌，接下來將針對在三維空間中的水下耐壓結構設計問題、以及探究考慮應力約束下的水下耐壓結構設計問題進行研究。