詹 可，蔣垣騰，趙 敏

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

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

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

美國(guó)海軍[2]2004年發(fā)布的水下無(wú)人潛水器(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]對(duì)該耐壓結(jié)構(gòu)進(jìn)行內(nèi)爆模擬，取得了與試驗(yàn)值吻合的數(shù)值計(jì)算結(jié)果，并分析了流場(chǎng)壓力波動(dòng)特性，如圖2所示，圖2 (a) 為測(cè)試玻璃球體及監(jiān)測(cè)裝置，圖2 (b) 為數(shù)值模擬7 MPa下內(nèi)爆后的壓力曲線。

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

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

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

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

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

1.1 有限元分析

考慮受簡(jiǎn)諧變化的壓力設(shè)計(jì)相關(guān)動(dòng)載荷拓?fù)鋬?yōu)化問(wèn)題，用有限元方法離散連續(xù)體結(jié)構(gòu)時(shí)，結(jié)構(gòu)在隨時(shí)間變化的力作用下的動(dòng)力響應(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ì)量和剛度矩陣，通過(guò)式(4)、(5)計(jì)算：

(4)

(5)

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

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

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

(6)

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

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

me=xeρe,ρ≥0.1

(7)

me=xeqρe,ρ<0.1

(8)

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

(9)

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

(10)

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

1.3 靈敏度分析

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

(11)

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

(12)

(13)

代入式(11)可得：

(14)

對(duì)于BILE模型(第2節(jié)中將給出具體的搜索方法)，單元密度被用于定義加載面，壓力載荷直接作用于加載節(jié)點(diǎn)，每個(gè)等效力在每次迭代中大致保持相同的大小。在兩次迭代之間，由于加載面的長(zhǎng)度和加載節(jié)點(diǎn)的數(shù)量同時(shí)增加或減少，每個(gè)等效力的大小都不會(huì)有顯著變化。因此，在整個(gè)優(yōu)化過(guò)程中，每個(gè)等效節(jié)點(diǎn)力的大小保持不變[23]。載荷對(duì)設(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個(gè)迭代步給出與其他算法相似的結(jié)構(gòu)，且本文為應(yīng)用研究，BILE模型能夠更快給出一個(gè)較為準(zhǔn)確的結(jié)果，有效提高工作效率。

該方法的基本流程為：

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

(17)

其中，xn,i為節(jié)點(diǎn)周圍的4個(gè)單元，xsn,l為當(dāng)前節(jié)點(diǎn)鄰接的4個(gè)節(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)化過(guò)程中的邊界節(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)力。對(duì)于一個(gè)加載節(jié)點(diǎn)，該點(diǎn)上力的角度由連接兩個(gè)最近的加載節(jié)點(diǎn)到垂直正方向的直線的法線逆時(shí)針來(lái)計(jì)算。

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

3 數(shù)值算例

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

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

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

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

作為承受壓力載荷的結(jié)構(gòu)優(yōu)化問(wèn)題中的經(jīng)典算例，選擇內(nèi)部受壓的拱形結(jié)構(gòu)優(yōu)化設(shè)計(jì)作為第一個(gè)算例，用以驗(yàn)證搜索算法在動(dòng)態(tài)問(wèn)題上的適用性。初始的設(shè)計(jì)域及約束條件等如圖4所示。動(dòng)載荷幅值為1.0，實(shí)體材料的材料屬性均采用無(wú)量綱化的取值，彈性模量取值為1.0，泊松比取值為0.3，質(zhì)量密度ρ=1.0×10-6，本節(jié)算例的上述參數(shù)取值相同。設(shè)計(jì)區(qū)域離散成 40×20個(gè)正方形單元，且最大允許的材料體積占初始總體積的體積分?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開(kāi)展了本文優(yōu)化結(jié)果的柔度與Xia等[28]研究和Zhang等[16]研究結(jié)果的比較，當(dāng)楊氏模量采用1.0，網(wǎng)格為100×50時(shí)，本文柔度值為10.02，小于表1中Xia等[28]研究的結(jié)果；當(dāng)楊氏模量采用100，網(wǎng)格為40×20時(shí)，本文的柔度值為0.07，略小于表1中Zhang等[28]研究的結(jié)果。本文的方法獲得了較小柔度的拓?fù)浣Y(jié)構(gòu)，優(yōu)于參考文獻(xiàn)的結(jié)果。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

除了上述開(kāi)展耐壓結(jié)構(gòu)平面應(yīng)變模型的研究外，一些學(xué)者[30-31]進(jìn)行了靜載環(huán)境軸對(duì)稱下水下耐壓結(jié)構(gòu)的設(shè)計(jì)研究。在前人的基礎(chǔ)上，開(kāi)展設(shè)計(jì)相關(guān)動(dòng)載荷作用下水下耐壓結(jié)構(gòu)的拓?fù)鋬?yōu)化研究，其設(shè)計(jì)域與邊界條件等如圖11(a)所示?？紤]到圖11(a)中的模型及邊界條件的對(duì)稱性，取 1/4 模型如圖11(b)所示作為設(shè)計(jì)域，并離散成120×48個(gè)正方形單元，且最大允許的材料體積占初始總體積的體積分?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é)果可以看到，隨著激勵(lì)頻率的增大，固定約束處的材料聚集明顯，整體材料分布有靠近約束的趨勢(shì)。

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

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

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

圖14 不同激勵(lì)頻率下固定約束處體積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ā)生變化，無(wú)法很好地反映材料變化對(duì)結(jié)構(gòu)模態(tài)剛度和結(jié)構(gòu)模態(tài)質(zhì)量的影響。因此，利用振動(dòng)力學(xué)模態(tài)質(zhì)量正則化方法[32]，可得相應(yīng)的正則模態(tài)剛度為：

(19)

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

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

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

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

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

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

為了進(jìn)一步說(shuō)明優(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)所示?？梢?jiàn)，本文中得到的多球交接耐壓結(jié)構(gòu)球殼交接處與圖17(a)、17(b)非常相似，在連接處會(huì)進(jìn)行局部加強(qiáng)，結(jié)構(gòu)一階固有頻率更高，動(dòng)力學(xué)性能更好，說(shuō)明本文對(duì)于水下耐壓結(jié)構(gòu)的探索具有一定的工程應(yīng)用價(jià)值。

圖17 考慮設(shè)計(jì)相關(guān)動(dòng)載荷得到的優(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é) 語(yǔ)

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

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

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