設(shè)計(jì)相關(guān)動載荷作用下的水下耐壓結(jié)構(gòu)拓?fù)鋬?yōu)化

詹 可，蔣垣騰，趙 敏

(上海交通大學(xué) 船舶海洋與建筑工程學(xué)院 海洋工程國家重點(diǎn)實(shí)驗(yàn)室，上海 200240)

潛水器是進(jìn)行深?？碧胶秃Ｑ罂茖W(xué)研究的重要裝備，也是一個國家科學(xué)技術(shù)水平的重要體現(xiàn)。耐壓殼作為潛水器的關(guān)鍵結(jié)構(gòu)，其設(shè)計(jì)對潛水器的安全性和空間利用率等性能具有重要影響，在保障設(shè)備工作和人員安全方面起著重要的作用。目前學(xué)者們往往按照傳統(tǒng)的金屬材料耐壓結(jié)構(gòu)設(shè)計(jì)方法開展水下耐壓殼結(jié)構(gòu)設(shè)計(jì)，即將靜水壓力作為環(huán)境載荷作用于耐壓結(jié)構(gòu)上，然后開展安全性設(shè)計(jì)。但是，水下耐壓結(jié)構(gòu)也會面臨一種特別的沖擊載荷——內(nèi)爆。當(dāng)水下耐壓結(jié)構(gòu)不能承受外部水壓而被壓潰塌陷時，流場靜水壓力轉(zhuǎn)化為流體動能，水流壓縮結(jié)構(gòu)至最小限度時，會發(fā)生水錘型的沖擊，水流動能轉(zhuǎn)化為沖擊波壓力對周圍結(jié)構(gòu)造成破壞，這種沖擊載荷就是內(nèi)爆載荷。單體結(jié)構(gòu)內(nèi)爆誘發(fā)多體結(jié)構(gòu)內(nèi)爆的現(xiàn)象就是殉爆。美國于2009年研制成功的11 000米級“海神號”無人潛水器搭載了大量陶瓷耐壓結(jié)構(gòu)，其結(jié)構(gòu)如圖1所示[1]，該潛水器在探索位于新西蘭的世界第二深海溝克馬德克海溝時下潛至9 990 m處由于單體陶瓷耐壓結(jié)構(gòu)內(nèi)爆引發(fā)多體殉爆，導(dǎo)致整個潛水器損毀，所有設(shè)備幾乎成為碎片。

圖1 “海神號”上陶瓷耐壓結(jié)構(gòu)Fig. 1 Ceramic pressure hulls on hybrid remotely operated vehicle Nereus

美國海軍[2]2004年發(fā)布的水下無人潛水器(UUV)計(jì)劃要求潛水器的設(shè)計(jì)必須嚴(yán)格考慮耐壓支承結(jié)構(gòu)內(nèi)爆所產(chǎn)生的影響。Turner[3]使用中空玻璃耐壓結(jié)構(gòu)在6.996 MPa的壓力筒中進(jìn)行水下內(nèi)爆試驗(yàn)，得到了沖擊波壓力曲線；陳鋒華和趙敏[4]對該耐壓結(jié)構(gòu)進(jìn)行內(nèi)爆模擬，取得了與試驗(yàn)值吻合的數(shù)值計(jì)算結(jié)果，并分析了流場壓力波動特性，如圖2所示，圖2 (a) 為測試玻璃球體及監(jiān)測裝置，圖2 (b) 為數(shù)值模擬7 MPa下內(nèi)爆后的壓力曲線。

圖2 測試玻璃球體及監(jiān)測裝置和數(shù)值模擬7 MPa下內(nèi)爆后的壓力曲線Fig. 2 Test stand with glass sphere, blade and three pressure sensors numerical simulation of pressure curve after implosion at 7 MPa

從圖2(b)可以看到，內(nèi)爆后產(chǎn)生的壓力大小隨著時間在不斷的發(fā)生變化，由于存在多次壓縮與反彈過程，壓力出現(xiàn)多個周期性峰值，對于周圍耐壓結(jié)構(gòu)而言，環(huán)境載荷不再保持為靜載。因此，在進(jìn)行耐壓殼優(yōu)化設(shè)計(jì)時，要考慮動態(tài)載荷即載荷隨時間變化對優(yōu)化結(jié)果的影響。圖2(b)中內(nèi)爆后產(chǎn)生的壓力存在多個周期性峰值，從載荷的類別上，可以將內(nèi)爆后產(chǎn)生的壓力以最高峰值為界分為兩個階段，從防護(hù)的角度研究水下耐壓結(jié)構(gòu)優(yōu)化設(shè)計(jì)，由于單體內(nèi)爆瞬間產(chǎn)生極大水錘型沖擊，第一階段可通過設(shè)置耐壓結(jié)構(gòu)防護(hù)裝置，防止陶瓷耐壓結(jié)構(gòu)單體內(nèi)爆導(dǎo)致的多體殉爆。第二階段可將周期性壓力變化簡化成簡諧載荷的形式，開展設(shè)計(jì)相關(guān)動載荷作用下的水下耐壓結(jié)構(gòu)拓?fù)鋬?yōu)化理論及方法研究。圖2(b)右上角為峰值后的壓力曲線放大圖，可以看到其載荷幅值波動在2 MPa以內(nèi)，其載荷變化形式類似于簡諧載荷。從橫軸的時間尺度來看，周期約為1 ms，頻率約為6 000 rad/s，即動載荷的頻率較高。在多耐壓結(jié)構(gòu)共存的情況下，若某一水下耐壓結(jié)構(gòu)內(nèi)爆產(chǎn)生周期性動態(tài)載荷，將會對其他耐壓結(jié)構(gòu)造成影響，僅考慮靜載的耐壓結(jié)構(gòu)傳統(tǒng)優(yōu)化設(shè)計(jì)方法在該工況下難以適用。因此，研究設(shè)計(jì)相關(guān)動載荷作用下的水下耐壓結(jié)構(gòu)優(yōu)化設(shè)計(jì)，是具有重要的理論意義和工程價值的。

目前，動力學(xué)相關(guān)的拓?fù)鋬?yōu)化主要集中在最大化結(jié)構(gòu)特征頻率設(shè)計(jì)[5-7]，以及以頻率為約束的優(yōu)化問題[8]。但是，關(guān)于結(jié)構(gòu)動響應(yīng)比如以結(jié)構(gòu)的動柔順度為目標(biāo)的研究較少。Ma等[9]首次將動柔順度的概念引入到結(jié)構(gòu)動力學(xué)拓?fù)鋬?yōu)化中，并以動柔順度為目標(biāo)，利用均勻化方法研究結(jié)構(gòu)動力學(xué)問題。Olhoff和Du[10-11]提出了頻率漸變的IF方法和GIF方法研究動柔順度相關(guān)問題，并且采用質(zhì)量修正的插值模型去避免虛假模態(tài)問題，其結(jié)果表明該方法簡單有效，求得的頻率值準(zhǔn)確。Liu等[12]給出了簡諧載荷下不同列式的拓?fù)鋬?yōu)化問題，用ESO算法給出了相應(yīng)的以動柔順度為目標(biāo)的結(jié)構(gòu)動力學(xué)優(yōu)化設(shè)計(jì)，得到了較為理想的結(jié)果。Silva等[13]研究了不同因素對穩(wěn)態(tài)強(qiáng)迫振動下的單材料結(jié)構(gòu)拓?fù)鋬?yōu)化的影響。以上研究均是關(guān)于固定載荷的優(yōu)化設(shè)計(jì)問題，然而本文物理背景下研究的問題是設(shè)計(jì)相關(guān)動載荷作用下的優(yōu)化問題，關(guān)于這類問題的研究相對較少，Olhoff和Du[10]和張暉等[14]利用其開發(fā)的邊界搜索算法來研究此類問題。張暉等[14]進(jìn)行了受內(nèi)壓的容器和水箱的優(yōu)化設(shè)計(jì)，其采用RAMP(rational approximation of material properties)插值模型來避免虛假模態(tài)，并對結(jié)構(gòu)拓?fù)涞淖兓o出了定性的解釋。但是，上述文獻(xiàn)并未將設(shè)計(jì)相關(guān)動載荷作用下的優(yōu)化問題引入到水下耐壓結(jié)構(gòu)設(shè)計(jì)問題中，并且對于得到的結(jié)構(gòu)沒有給出更為深入的解釋。

本文研究聚焦于設(shè)計(jì)相關(guān)動載荷下的水下耐壓結(jié)構(gòu)拓?fù)鋬?yōu)化問題，其中設(shè)計(jì)相關(guān)載荷的邊界搜索是區(qū)別于固定載荷動力學(xué)問題的關(guān)鍵環(huán)節(jié)。目前關(guān)于加載面搜索算法的研究主要分為兩類：一類是通過對節(jié)點(diǎn)密度的處理來搜索加載面，Du與Olhoff[15]提出利用貝塞爾樣條曲線代表加載邊界，由節(jié)點(diǎn)相對密度值插值得到的等值點(diǎn)來作為樣條曲線控制點(diǎn)的等值線方法和Zhang等[16]提出的單元節(jié)點(diǎn)加載法均屬于此類，該類方法的優(yōu)點(diǎn)是易于實(shí)現(xiàn)，而且避免了載荷靈敏度分析；另一類通過構(gòu)造多物理場來跟蹤壓力載荷的變化，Chen和Kikuchi[17]、Bourdin和Chambolle[18]分別提出的通過流體流動來識別壓力載荷表面的技術(shù)。Kumar等[19]利用巖石力學(xué)中的Darcy定律來處理設(shè)計(jì)相關(guān)載荷。Wang和Qian[20]提出的基于密度梯度的方法來跟蹤加載面。多物理場的優(yōu)勢在于壓力直接作用于兩相交界面上，無需進(jìn)行加載面搜索，但其有限元建模相對復(fù)雜。除上述兩類搜索算法外，也有一些學(xué)者進(jìn)行了水平集方法下的加載面搜索算法研究，郭旭和趙康[21]，Jiang和Zhao[22]利用水平集演化技術(shù)，提出了一種設(shè)計(jì)相關(guān)載荷作用邊界的搜索方法，通過適當(dāng)?shù)臄?shù)學(xué)變換，方便的處理施加在結(jié)構(gòu)上的拓?fù)湎嚓P(guān)載荷，避免了在密度法中繁瑣的邊界提取工作。

采用Ibhadode等[23]提出的BILE(boundary identification-load evolution)模型作為邊界搜索算法，這一邊界識別和載荷演化的模型通過設(shè)置單元偽密度閾值來搜索加載面，該閾值在每次迭代中通過特定的方式遞增以進(jìn)行邊界識別。同時，引入一個參數(shù)來定義兩個邊界識別步驟之間載荷演化的迭代次數(shù)，也可以控制優(yōu)化的速度。BILE模型屬于第一類加載面搜索算法，無需進(jìn)行載荷的靈敏度分析。因此，該模型易于應(yīng)用，而且在80～100次迭代步后得到的數(shù)值結(jié)果都具有不錯的計(jì)算速度和可行性[23]。

通過將內(nèi)爆后耐壓結(jié)構(gòu)所受載荷簡化為周期性載荷，研究設(shè)計(jì)相關(guān)動載荷作用下水下耐壓結(jié)構(gòu)的拓?fù)鋬?yōu)化理論及方法，采用SIMP(solid isotropic microstructures with penalization)方法，選擇Ibhadode等[23]提出的BILE模型作為邊界搜索算法，并結(jié)合杜建鑌[24]提出的修正的材料插值模型來避免局部模態(tài)現(xiàn)象，采用靜柔度權(quán)重因子的拓?fù)鋬?yōu)化數(shù)學(xué)模型來解決高頻下的水下耐壓結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)，給出在不同頻率下的水下耐壓結(jié)構(gòu)最優(yōu)構(gòu)型，指導(dǎo)耐壓結(jié)構(gòu)的概念設(shè)計(jì)，為結(jié)構(gòu)設(shè)計(jì)人員提供思路。

1 問題提出

1.1 有限元分析

考慮受簡諧變化的壓力設(shè)計(jì)相關(guān)動載荷拓?fù)鋬?yōu)化問題，用有限元方法離散連續(xù)體結(jié)構(gòu)時，結(jié)構(gòu)在隨時間變化的力作用下的動力響應(yīng)平衡方程可寫成：

(1)

不考慮阻尼C=0，且在諧響應(yīng)的情況下可以得到:

(K-w2M)u=F

(2)

其中，M,K分別代表了結(jié)構(gòu)的總質(zhì)量矩陣和剛度矩陣，利用有限元方法可得：

(3)

其中，n是設(shè)計(jì)域中單元的總數(shù)。mi和ki是單元質(zhì)量和剛度矩陣，通過式(4)、(5)計(jì)算：

(4)

(5)

其中，V是單元的體積，ρ和D分別是質(zhì)量密度和本構(gòu)矩陣，N和B表示單元形狀函數(shù)和應(yīng)變位移矩陣，mc為單元的一致質(zhì)量矩陣，ml為單元的集中質(zhì)量矩陣。本文采用形如式(4)中的組合形式來減少計(jì)算誤差。

1.2 拓?fù)鋬?yōu)化數(shù)學(xué)模型

基于密度法的連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化問題，將設(shè)計(jì)變量松弛為0～1之間任意值的連續(xù)變量優(yōu)化問題。由于這樣的松弛方式將導(dǎo)致最優(yōu)結(jié)構(gòu)中出現(xiàn)中間密度的材料，為了減少灰度單元，常采用懲罰因子來對單元偽密度進(jìn)行懲罰。最常采用的是SIMP方法，設(shè)給定設(shè)計(jì)區(qū)域?yàn)棣?，用N個單元離散，每個有限單元的相對密度為：xi(i=1,2……N)。對于各項(xiàng)同性材料，由SIMP模型可以得到單元的彈性模量：

(6)

式中：E0為固體材料的彈性模量；p為懲罰系數(shù)，在文中取3；Emin是為了避免總體剛度矩陣奇異性而引入的非零常數(shù)，取0.001。

本研究給定體積約束下最小化結(jié)構(gòu)的動柔順度問題，當(dāng)使用標(biāo)準(zhǔn)的SIMP法進(jìn)行動力學(xué)優(yōu)化時，存在局部模態(tài)(也稱虛假模態(tài))現(xiàn)象。局部模態(tài)現(xiàn)象是指在優(yōu)化過程中，設(shè)計(jì)域內(nèi)會出現(xiàn)一些低密度區(qū)域，由于這些區(qū)域的剛度與質(zhì)量比遠(yuǎn)低于其他位置，會導(dǎo)致最后求解的振動模態(tài)出現(xiàn)劇烈的振蕩。因此，在動力學(xué)拓?fù)鋬?yōu)化中需要對局部模態(tài)進(jìn)行處理。Lazarus和Hagiwara[25]建議采用0.1作為低密度區(qū)的門檻值，Pedersen[6]提出一種修正的剛度插值模型；杜建鑌[24]借鑒Pedersen[6]提出的插值形式的思路，提出對低密度區(qū)的質(zhì)量插值模型進(jìn)行修正，采用杜建鑌[24]提出的質(zhì)量修正模型：

me=xeρe,ρ≥0.1

(7)

me=xeqρe,ρ<0.1

(8)

其中，me為單元的質(zhì)量，xe為單元的相對密度，ρe為單元的物理密度,q是進(jìn)行單元質(zhì)量插值時的懲罰因子，文中q取為 6。因此，優(yōu)化問題的離散數(shù)學(xué)模型為：

(9)

當(dāng)頻率接近或高于一階固有頻率時，采用上述優(yōu)化模型無法得到可靠的結(jié)果(3.1節(jié)會給出相關(guān)算例來說明)，需引入靜柔度約束，或目標(biāo)函數(shù)中引入含靜柔度的權(quán)重因子[13]，采用第二種方式，其數(shù)學(xué)模型為：

(10)

其中，F(xiàn)為設(shè)計(jì)相關(guān)動載荷的載荷幅值向量，U為位移幅值向量，K、M分別為整體剛度和整體質(zhì)量矩陣，w為激勵頻率，Cd為結(jié)構(gòu)動柔順度，Cs為靜柔度，計(jì)算KU=F得到，V為實(shí)際材料體積，V0為設(shè)計(jì)域體積，f為給定材料體分比，n為設(shè)計(jì)域中單元的個數(shù)。

1.3 靈敏度分析

為了解決式(9)中列出的優(yōu)化問題，需要目標(biāo)函數(shù)對設(shè)計(jì)變量的導(dǎo)數(shù)。采用優(yōu)化準(zhǔn)則法對優(yōu)化問題進(jìn)行求解。目標(biāo)函數(shù)對設(shè)計(jì)變量的求導(dǎo)為：

(11)

根據(jù)KdU=F，可得：

(12)

(13)

代入式(11)可得：

(14)

對于BILE模型(第2節(jié)中將給出具體的搜索方法)，單元密度被用于定義加載面，壓力載荷直接作用于加載節(jié)點(diǎn)，每個等效力在每次迭代中大致保持相同的大小。在兩次迭代之間，由于加載面的長度和加載節(jié)點(diǎn)的數(shù)量同時增加或減少，每個等效力的大小都不會有顯著變化。因此，在整個優(yōu)化過程中，每個等效節(jié)點(diǎn)力的大小保持不變[23]。載荷對設(shè)計(jì)變量的伴隨靈敏度為0，則靈敏度為：

(15)

即：

(16)

2 BILE模型

在設(shè)計(jì)相關(guān)載荷作用下結(jié)構(gòu)拓?fù)鋬?yōu)化的研究難點(diǎn)在于如何確定設(shè)計(jì)相關(guān)載荷作用的加載面。采用Ibhadode等[23]提出的BILE模型作為加載面搜索算法，因?yàn)槠涫諗克俣瓤?，能夠?0～100個迭代步給出與其他算法相似的結(jié)構(gòu)，且本文為應(yīng)用研究，BILE模型能夠更快給出一個較為準(zhǔn)確的結(jié)果，有效提高工作效率。

該方法的基本流程為：

1)確定每次迭代的邊界節(jié)點(diǎn)。根據(jù)單元密度閾值以及每個節(jié)點(diǎn)相鄰單元的平均密度值確定邊界節(jié)點(diǎn)如圖3(a)和3(b)，第一次迭代時由于各單元密度值相等，邊界節(jié)點(diǎn)為設(shè)計(jì)域四周節(jié)點(diǎn)。而迭代過程中，某一節(jié)點(diǎn)n的選取公式為:

(17)

其中，xn,i為節(jié)點(diǎn)周圍的4個單元，xsn,l為當(dāng)前節(jié)點(diǎn)鄰接的4個節(jié)點(diǎn)Sn的四周單元，若當(dāng)前節(jié)點(diǎn)四周的單元密度符合公式的選取規(guī)則，則當(dāng)前節(jié)點(diǎn)被標(biāo)記為邊界節(jié)點(diǎn)，某一次的選取如圖3(b)所示。

2)從生成的邊界節(jié)點(diǎn)中選擇加載節(jié)點(diǎn)[如圖3(c)]形成加載面。在BILE模型中，選取最外側(cè)邊界點(diǎn)作為加載節(jié)點(diǎn)，具體是否選取某一邊的加載節(jié)點(diǎn)取決于初始邊界條件和邊界條件與加載節(jié)點(diǎn)關(guān)系。

圖3 節(jié)點(diǎn)n的鄰接節(jié)點(diǎn)以及某一次優(yōu)化過程中的邊界節(jié)點(diǎn)和載荷節(jié)點(diǎn)Fig. 3 Adjacent node of node n, boundary node in a certain optimization process, load node in a certain optimization process

3)施加等效節(jié)點(diǎn)力。對于一個加載節(jié)點(diǎn)，該點(diǎn)上力的角度由連接兩個最近的加載節(jié)點(diǎn)到垂直正方向的直線的法線逆時針來計(jì)算。

由于是通過單元偽密度閾值并人為設(shè)定一定規(guī)則來搜索壓力加載面，不必求解力的伴隨靈敏度項(xiàng)，這樣不僅簡化了目標(biāo)函數(shù)靈敏度的求解過程，而且也加快求解的速度，提高了整體的計(jì)算效率。

3 數(shù)值算例

由于實(shí)際結(jié)構(gòu)均存在阻尼，阻尼對高頻響應(yīng)有很強(qiáng)的抑制作用，對低頻響應(yīng)的影響較小[24]，這里主要關(guān)注小于或略高于結(jié)構(gòu)一階固有頻率的低頻共振優(yōu)化問題，研究不同設(shè)計(jì)域和約束以及頻率條件對結(jié)構(gòu)拓?fù)涞挠绊憽?/p>

Synergy HT酶標(biāo)儀（美國BioTeK公司）；Mini-PROTEAN Tetra蛋白電泳儀、Trans-Blot SD半干轉(zhuǎn)膜系統(tǒng)（美國Bio-rad公司）；Direct-Q超純水儀（美國Millipore公司）；5417R高速冷凍離心機(jī)（德國Eppendorf公司）；IKA T18 basic勻漿器（德國ULTRATURRAX公司）；AX70顯微照相系統(tǒng)（日本Olympus公司）；ImageQuant LAS 4000全自動圖像分析系統(tǒng)（美國GE公司）。

采用優(yōu)化準(zhǔn)則法進(jìn)行設(shè)計(jì)變量的更新，利用MATLAB中的eigs函數(shù)求解結(jié)構(gòu)頻率，其求解對稱矩陣采用Lanczos算法。在本文算例中，Sigmund[26]、Sigmund和Maute[27]提出的靈敏度過濾方法被用來消除計(jì)算中的棋盤格，網(wǎng)格依賴性等數(shù)值不穩(wěn)定現(xiàn)象。

3.1 受內(nèi)壓的拱形結(jié)構(gòu)

作為承受壓力載荷的結(jié)構(gòu)優(yōu)化問題中的經(jīng)典算例，選擇內(nèi)部受壓的拱形結(jié)構(gòu)優(yōu)化設(shè)計(jì)作為第一個算例，用以驗(yàn)證搜索算法在動態(tài)問題上的適用性。初始的設(shè)計(jì)域及約束條件等如圖4所示。動載荷幅值為1.0，實(shí)體材料的材料屬性均采用無量綱化的取值，彈性模量取值為1.0，泊松比取值為0.3，質(zhì)量密度ρ=1.0×10-6，本節(jié)算例的上述參數(shù)取值相同。設(shè)計(jì)區(qū)域離散成 40×20個正方形單元，且最大允許的材料體積占初始總體積的體積分?jǐn)?shù)為50%。

圖4 初始設(shè)計(jì)域及邊界條件Fig. 4 Design domain and boundary conditions

此前，已有許多學(xué)者進(jìn)行了設(shè)計(jì)相關(guān)載荷下的拓?fù)鋬?yōu)化研究。圖5表明本文優(yōu)化結(jié)果的拓?fù)湫问脚cZhang等[16]研究以及Xia等[28]研究的結(jié)果類似。表1開展了本文優(yōu)化結(jié)果的柔度與Xia等[28]研究和Zhang等[16]研究結(jié)果的比較，當(dāng)楊氏模量采用1.0，網(wǎng)格為100×50時，本文柔度值為10.02，小于表1中Xia等[28]研究的結(jié)果；當(dāng)楊氏模量采用100，網(wǎng)格為40×20時，本文的柔度值為0.07，略小于表1中Zhang等[28]研究的結(jié)果。本文的方法獲得了較小柔度的拓?fù)浣Y(jié)構(gòu)，優(yōu)于參考文獻(xiàn)的結(jié)果。

圖5 本文結(jié)果與不同文獻(xiàn)結(jié)果的對比Fig. 5 Comparison of optimization results between this paper and other references

表1 設(shè)計(jì)相關(guān)載荷下拓?fù)鋬?yōu)化的數(shù)值示例對比：受內(nèi)壓拱形結(jié)構(gòu)Tab. 1 Comparison of numerical examples of topology optimization under design-dependent loads: arch structure subjected to internal pressure

圖6給出了采用式(9)的優(yōu)化模型時隨激勵頻率變化目標(biāo)函數(shù)以及拓?fù)湫问降淖兓?。可以看到，?dāng)激勵頻率低于500 rad/s時，隨著激勵頻率的增大，雖然優(yōu)化結(jié)果的動柔度略有增大，但拓?fù)浣Y(jié)構(gòu)形式并無明顯差異。當(dāng)激勵頻率為200 rad/s時，本文優(yōu)化結(jié)果與張暉等[14]的研究結(jié)果對比如圖7(a)和7(c)所示，拓?fù)湫问綆缀跻恢隆ｋS著激勵頻率增大到600 rad/s時，由于頻率過大，求解并不收斂，中間結(jié)果如圖6所示。圖7(b)為激勵頻率800 rad/s時的中間結(jié)果，雖然與圖7 (d)張暉等[14]的研究結(jié)果類似，但繼續(xù)優(yōu)化后求解并不收斂，最終結(jié)果發(fā)散；張暉等[14]在其論文中已經(jīng)說明其結(jié)果不滿足體積約束，非可行解。由此可見，傳統(tǒng)的以單一動柔度為目標(biāo)的拓?fù)鋬?yōu)化數(shù)學(xué)模型存在缺陷。

圖6 不同激勵頻率下的優(yōu)化結(jié)果對比曲線Fig. 6 Comparison curve of optimization results under different excitation frequencies

圖7 不同激勵頻率下本文優(yōu)化結(jié)果與參考文獻(xiàn)[14]的對比Fig. 7 Comparison of optimization results between this paper and reference [14] with different excitation frequencies

因此，引入了如式(10)的目標(biāo)函數(shù)——含靜柔度權(quán)重因子的優(yōu)化模型來計(jì)算高頻下的結(jié)果，其有效性已在Silva等[13]研究中得到證明。該方法中η的選取對結(jié)果影響較大。表2給出了激勵頻率為800 rad/s時，不同η下優(yōu)化結(jié)果，其中nlter為迭代步數(shù)。

表2 激勵頻率為800 rad/s時不同的η對結(jié)果的影響Tab. 2 The effect of different η on the result when excitation fequency is 800 rad/s

當(dāng)η大于0.6時，拱形結(jié)構(gòu)中會存在一些灰度單元，且結(jié)果并不對稱；隨著η的減小，優(yōu)化結(jié)果呈現(xiàn)“0～1”分布，且η為0.4時結(jié)構(gòu)的一階固有頻率為808.2 rad/s，比靜載時的結(jié)構(gòu)頻率750.5 rad/s高；因此，在利用式(10)的優(yōu)化模型進(jìn)行受內(nèi)壓的拱形結(jié)構(gòu)優(yōu)化時，建議η在0.2～0.4取值。

3.2 設(shè)計(jì)相關(guān)動載荷下的水下耐壓結(jié)構(gòu)平面應(yīng)變模型優(yōu)化設(shè)計(jì)

圖8 設(shè)計(jì)域與邊界條件Fig. 8 The design domain and boundary conditions

圖9和圖10分別為不同激勵頻率下的優(yōu)化結(jié)果以及激勵頻率200 rad/s時的優(yōu)化過程。

圖9 六種激勵頻率下的最優(yōu)拓?fù)銯ig. 9 Optimal topology under six excitation frequencies

圖10 設(shè)計(jì)問題的優(yōu)化過程Fig. 10 Optimization history for the design problem

圖9(a)中為靜水壓(激勵頻率為0 rad/s)下的優(yōu)化結(jié)果，與以往學(xué)者[30-31]所得到靜水壓下的結(jié)果基本一致，其一階固有頻率為418.1 rad/s。圖9(b)～(f)中隨著激勵頻率的增大，由于算法振蕩，導(dǎo)致優(yōu)化結(jié)構(gòu)的拓?fù)錁?gòu)型在局部邊界處灰度單元減少，但優(yōu)化結(jié)構(gòu)的整體構(gòu)型基本保持不變。在設(shè)計(jì)域中無固定端等強(qiáng)約束時，圖8所示的初始設(shè)計(jì)域在低于靜水壓下結(jié)構(gòu)的一階固有頻率的不同激勵頻率下，圓環(huán)均為最優(yōu)的結(jié)構(gòu)形式，進(jìn)一步驗(yàn)證了圓柱殼作為耐壓結(jié)構(gòu)的合理性。當(dāng)激勵頻率高于靜水壓下結(jié)構(gòu)的一階固有頻率(418.1 rad/s)時，采用式(10)的模型來進(jìn)行優(yōu)化。激勵頻率為500 rad/s時不同η下的結(jié)果對比如表3。

表3 激勵頻率為500 rad/s時η對結(jié)果的影響Tab. 3 The effect of η on the result when excitation fequency is 500 rad/s

從表4可以發(fā)現(xiàn)，當(dāng)η為0.80時，表3中優(yōu)化結(jié)果的一階固有頻率較圓環(huán)型耐壓結(jié)構(gòu)小，且并不呈現(xiàn)出中心對稱的構(gòu)型，不可取；隨著η減小到0.40，優(yōu)化結(jié)果呈現(xiàn)較為清晰的“0～1”分布，且表4中顯示結(jié)構(gòu)的一階固有頻率普遍高于圓環(huán)型耐壓結(jié)構(gòu)的一階固有頻率，但大于0.30時，結(jié)構(gòu)對稱性沒有小于0.20時好，在利用式(10)的優(yōu)化模型進(jìn)行設(shè)計(jì)相關(guān)動載荷下的水下耐壓結(jié)構(gòu)平面應(yīng)變模型優(yōu)化設(shè)計(jì)時，建議η在0.05～0.20內(nèi)取值。

表4 激勵頻率為500 rad/s時不同η下結(jié)果的第一階固有頻率Tab. 4 Fundamental frequency of results under different η when excitation fequency is 500 rad/s

3.3 設(shè)計(jì)相關(guān)動載荷下的水下耐壓結(jié)構(gòu)軸對稱模型優(yōu)化設(shè)計(jì)

除了上述開展耐壓結(jié)構(gòu)平面應(yīng)變模型的研究外，一些學(xué)者[30-31]進(jìn)行了靜載環(huán)境軸對稱下水下耐壓結(jié)構(gòu)的設(shè)計(jì)研究。在前人的基礎(chǔ)上，開展設(shè)計(jì)相關(guān)動載荷作用下水下耐壓結(jié)構(gòu)的拓?fù)鋬?yōu)化研究，其設(shè)計(jì)域與邊界條件等如圖11(a)所示?？紤]到圖11(a)中的模型及邊界條件的對稱性，取 1/4 模型如圖11(b)所示作為設(shè)計(jì)域，并離散成120×48個正方形單元，且最大允許的材料體積占初始總體積的體積分?jǐn)?shù)為35%。

圖11 設(shè)計(jì)域與邊界條件Fig. 11 The design domain and boundary condition

圖12(a)中結(jié)構(gòu)的第一階固有頻率值為1 253.8 rad/s。從圖12(b)、12(c)及圖13優(yōu)化結(jié)果可以看到，隨著激勵頻率的增大，固定約束處的材料聚集明顯，整體材料分布有靠近約束的趨勢。

圖12 三種激勵頻率下的1/4軸對稱模型最優(yōu)拓?fù)銯ig. 12 Optimal topology under three excitation frequencies of 1/4 axisymmetric model

圖13 設(shè)計(jì)問題的優(yōu)化過程Fig. 13 Optimization history for the design problem

定義固定端約束兩側(cè)5列單元密度和為約束處面積，表5及圖14反映出隨著激勵頻率增大，約束處面積不斷增大，這與杜建鑌[24]所做的流道設(shè)計(jì)中隨著激勵頻率變大，頂部稍薄，支撐處有所加強(qiáng)有相似之處。為了解釋這種現(xiàn)象，考慮不同頻率下最終優(yōu)化構(gòu)型的第一階模態(tài)剛度與模態(tài)質(zhì)量比值的變化：

圖14 不同激勵頻率下固定約束處體積Fig. 14 The volume of the fixed restraint under different excitation frequencies

表5 不同頻率下約束處體積Tab. 5 Restricted volume at different frequencies

(18)

其中，wn為結(jié)構(gòu)的圓頻率，ur為結(jié)構(gòu)的第r階模態(tài)向量，K、M分別為結(jié)構(gòu)總的剛度矩陣和質(zhì)量矩陣。由于頻率不同，結(jié)構(gòu)的第r階模態(tài)剛度和第r階模態(tài)質(zhì)量均在發(fā)生變化，無法很好地反映材料變化對結(jié)構(gòu)模態(tài)剛度和結(jié)構(gòu)模態(tài)質(zhì)量的影響。因此，利用振動力學(xué)模態(tài)質(zhì)量正則化方法[32]，可得相應(yīng)的正則模態(tài)剛度為：

(19)

通過比較正則模態(tài)剛度的變化，可以比較材料分布變化對結(jié)構(gòu)一階固有頻率的影響。圖15給出了不同激勵頻率下優(yōu)化結(jié)構(gòu)的第一階正則模態(tài)剛度變化情況。隨著激勵頻率的增大，優(yōu)化結(jié)構(gòu)的正則模態(tài)剛度也在不斷增大。因此，拓?fù)湫问缴系牟煌?，?shí)際上反映的是優(yōu)化過程中材料的增加和去除對結(jié)構(gòu)一階模態(tài)剛度與一階模態(tài)質(zhì)量比值的影響，若某位置材料的增加或去除能使其比值變大，則能夠提高結(jié)構(gòu)的一階固有頻率。圖16優(yōu)化結(jié)果的一階固有頻率隨激勵頻率的變化情況也驗(yàn)證了上述觀點(diǎn)。因此，對于多球交接的水下耐壓結(jié)構(gòu)形式，在多球交接的位置采用環(huán)向加強(qiáng)肋可以提高結(jié)構(gòu)的動力學(xué)特性，使耐壓殼能夠適應(yīng)更加復(fù)雜的水下環(huán)境，提高潛水器抗沖擊性能。

圖15 優(yōu)化結(jié)果一階正則模態(tài)剛度隨激勵頻率變化Fig. 15 The first-order canonical mode stiffness varies with the excitation frequency

圖16 優(yōu)化結(jié)果一階固有頻率隨激勵頻率變化Fig. 16 The first-order frequency varies with the excitation frequency

當(dāng)激勵頻率高于靜水壓下結(jié)構(gòu)的第一階固有頻率(1 253.8 rad/s)時，采用式(10)的數(shù)學(xué)模型來進(jìn)行優(yōu)化。表6給出了激勵頻率為1 300 rad/s時不同η的拓?fù)鋬?yōu)化結(jié)果。

表6 激勵頻率為1 300 rad/s時不同η對結(jié)果的影響Tab. 6 The effect of different η on the result when excitation fequency is 1 300 rad/s

當(dāng)η大于0.20時，優(yōu)化結(jié)果出現(xiàn)較多的灰度單元，結(jié)果不可取；隨著η的減小，優(yōu)化結(jié)果呈現(xiàn)較為清晰的“0～1”分布。因此，在利用式(10)的優(yōu)化模型進(jìn)行設(shè)計(jì)相關(guān)動載荷作用下的水下耐壓結(jié)構(gòu)軸對稱模型優(yōu)化設(shè)計(jì)時，建議η在0.01～0.15內(nèi)取值。

為了進(jìn)一步說明優(yōu)化結(jié)構(gòu)的應(yīng)用可行性，圖17(a)、17(b)分別給出了MIT 水下潛水器設(shè)計(jì)團(tuán)隊(duì)設(shè)計(jì)的一種多球交接耐壓殼體模型[33]及張建等[34]提出的多蛋交接耐壓殼三蛋形殼截面，本文得到的多球交接耐壓結(jié)構(gòu)如圖17(c)所示?？梢姡疚闹械玫降亩嗲蚪唤幽蛪航Y(jié)構(gòu)球殼交接處與圖17(a)、17(b)非常相似，在連接處會進(jìn)行局部加強(qiáng)，結(jié)構(gòu)一階固有頻率更高，動力學(xué)性能更好，說明本文對于水下耐壓結(jié)構(gòu)的探索具有一定的工程應(yīng)用價值。

圖17 考慮設(shè)計(jì)相關(guān)動載荷得到的優(yōu)化結(jié)果與其他文獻(xiàn)結(jié)果的比較Fig. 17 Comparison of optimization results between this paper under design-dependent dynamic loading and other references

4 結(jié) 語

通過將內(nèi)爆后耐壓結(jié)構(gòu)第二階段所受高頻率的周期性變化載荷簡化為簡諧載荷，基于Ibhadode等提出的BILE加載面搜索算法，研究了設(shè)計(jì)相關(guān)動載荷作用下的水下耐壓結(jié)構(gòu)拓?fù)鋬?yōu)化研究及方法。同時，采用了靜柔度權(quán)重因子的拓?fù)鋬?yōu)化數(shù)學(xué)模型來解決高頻下的水下耐壓結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)，數(shù)值算例驗(yàn)證算法有效性。

著重研究了設(shè)計(jì)相關(guān)動載荷作用下的水下耐壓結(jié)構(gòu)拓?fù)鋬?yōu)化問題，探索了水下耐壓結(jié)構(gòu)新形式。低頻時，圓環(huán)型耐壓結(jié)構(gòu)與靜水壓下相似，但對于多球交接耐壓結(jié)構(gòu)，隨著頻率增加，材料往約束處聚集以提升結(jié)構(gòu)的動力學(xué)性能；接近或高于一階頻率時，圓環(huán)型耐壓結(jié)構(gòu)與多球交接型耐壓結(jié)構(gòu)形式均與靜水壓下的結(jié)果存在明顯差異，多球交接型耐壓結(jié)構(gòu)與現(xiàn)有的耐壓結(jié)構(gòu)概念設(shè)計(jì)結(jié)果相似。最終，探究了其在工程實(shí)際中的應(yīng)用，對新型水下耐壓結(jié)構(gòu)的概念設(shè)計(jì)具有積極的意義。

本文的研究只是對設(shè)計(jì)相關(guān)動載荷作用的水下耐壓結(jié)構(gòu)優(yōu)化設(shè)計(jì)的初步探究，在未來的工作中，將致力于該方法的三維拓展及屈曲約束下的水下耐壓結(jié)構(gòu)優(yōu)化設(shè)計(jì)以及其對應(yīng)的三維拓?fù)鋬?yōu)化問題。