張曉頔， 郝志明， 李一堃， 劉遠(yuǎn)東

(1.中國工程物理研究院總體工程研究所， 四川 綿陽 621900； 2.成都大學(xué)高等研究院， 成都 610106)

連接結(jié)構(gòu)在航空航天、軍事裝備和核電等領(lǐng)域中均有廣泛應(yīng)用[1-2]，大型裝備如航天器、軍用車輛、反應(yīng)堆壓力容器等是由不同的零部件通過螺栓預(yù)緊連接、鉚釘連接、楔環(huán)連接等多種連接方式連接而成。在結(jié)構(gòu)系統(tǒng)的拓?fù)鋬?yōu)化設(shè)計(jì)中，認(rèn)為連接結(jié)構(gòu)的接觸表面不發(fā)生相對運(yùn)動，將其簡化為線性模型分析，計(jì)算效率更高。實(shí)際上連接結(jié)構(gòu)的剛度和阻尼具有不連續(xù)性、遲滯性等非線性特性，分布于結(jié)構(gòu)系統(tǒng)不同位置的剛度不同，對整體結(jié)構(gòu)的力學(xué)特性的影響也不同[3]。由于連接接觸界面在外載荷作用下發(fā)生微觀滑移現(xiàn)象甚至宏觀滑移現(xiàn)象，引起能量耗散并產(chǎn)生干摩擦阻尼，占結(jié)構(gòu)整體阻尼的90%；同時連接接觸界面的滑移還會引起結(jié)構(gòu)局部剛度非線性變化。因此，在對復(fù)雜工程結(jié)構(gòu)開展拓?fù)鋬?yōu)化研究時，考慮其連接部分的非線性具有重要的意義。

結(jié)構(gòu)拓?fù)鋬?yōu)化研究經(jīng)歷了從離散體結(jié)構(gòu)到連續(xù)體結(jié)構(gòu)的發(fā)展過程。從早期由Michell的桁架理論和Dorn等人提出的基結(jié)構(gòu)法建立起離散結(jié)構(gòu)的拓?fù)鋬?yōu)化方法，發(fā)展到由Bends?e和Kikuchi提出的結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)均勻化方法，標(biāo)志著連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)的誕生。Xie[4]又提出了一種連續(xù)體拓?fù)鋬?yōu)化方法——漸進(jìn)結(jié)構(gòu)優(yōu)化方法(Evolutionary Structural Optimization，ESO)，與均勻化方法相比，它不產(chǎn)生0-1之間的中間密度，拓?fù)湫问角逦?，也避免了多變量的?shù)學(xué)規(guī)劃求解，這種方法已成功應(yīng)用于大量工程優(yōu)化問題的求解。目前關(guān)于線彈性材料和小變形的結(jié)構(gòu)拓?fù)鋬?yōu)化研究較為充分，但考慮結(jié)構(gòu)非線性行為的研究仍相對較少[5]。Ryu等[6]率先進(jìn)行了非線性響應(yīng)的設(shè)計(jì)靈敏度分析，開辟了非線性結(jié)構(gòu)優(yōu)化研究的先例。Buhl等[7]研究了幾何非線性結(jié)構(gòu)的拓?fù)鋬?yōu)化，指出非線性和線性情況下優(yōu)化結(jié)果的區(qū)別。Yuge等[8]將材料非線性引入結(jié)構(gòu)拓?fù)鋬?yōu)化，改進(jìn)彈塑性均勻化方法，提高了計(jì)算效率。Jung等[9]在優(yōu)化中同時考慮了幾何非線性和材料非線性，進(jìn)行了數(shù)值驗(yàn)證?？紤]界面非線性的結(jié)構(gòu)優(yōu)化研究文獻(xiàn)很少，僅在尺寸優(yōu)化和形狀優(yōu)化方面有少量工作[10-11]，尚未見界面非線性問題的拓?fù)鋬?yōu)化研究報(bào)道。

開展含螺栓連接結(jié)構(gòu)的非線性拓?fù)鋬?yōu)化設(shè)計(jì)，需模擬連接界面的接觸行為。傳統(tǒng)摩擦模型[12]如庫倫摩擦模型在有限元分析中主要采用接觸算法，計(jì)算較為耗時且不易收斂，難以應(yīng)用于拓?fù)鋬?yōu)化分析。因此，尋求一種更加合理的模型描述連接界面的接觸行為十分必要。連接接觸本構(gòu)模型[13]是一類能描述摩擦接觸非線性遲滯行為的模型，如Iwan模型、Bouc-Wen模型、Valanis模型和剪切層模型等。其中Iwan模型[14-15]由一系列理想彈塑性元件組成，最初用以描述金屬材料的彈塑性力學(xué)行為。Segalman[16]最早采用并聯(lián)-串聯(lián)形式的Iwan模型描述連接接觸非線性力學(xué)行為，提出了一個四參數(shù)Iwan模型，能夠準(zhǔn)確反映微觀滑移能量耗散的冪次關(guān)系，但不能描述發(fā)生宏觀滑移后的殘余剛度現(xiàn)象。在四參數(shù)Iwan模型的基礎(chǔ)上，Li等[17-18]提出了一個含截?cái)鄡缏煞植己碗p脈沖的非均勻密度函數(shù)的六參數(shù)Iwan模型，可以同時描述微觀滑移階段能量耗散冪次關(guān)系和宏觀滑移階段殘余剛度現(xiàn)象。此外，該模型引入有限元分析，不需要采用接觸算法，計(jì)算易于收斂[19-20]。

本文旨在建立一種含界面強(qiáng)非線性的靜力學(xué)拓?fù)鋬?yōu)化方法。采用六參數(shù)Iwan模型描述含螺栓連接的非線性力學(xué)行為，ESO方法結(jié)合結(jié)構(gòu)非線性分析，用ANSYS軟件的參數(shù)化設(shè)計(jì)語言(APDL語言)編寫命令流文件，形成一套優(yōu)化算法，應(yīng)用于含螺栓連接圓筒結(jié)構(gòu)的拓?fù)鋬?yōu)化設(shè)計(jì)。

1 螺栓連接的非線性力學(xué)模型

Iwan模型由若干Jenkins單元(彈簧-滑塊)構(gòu)成。圖1所示為并聯(lián)-串聯(lián)型Iwan模型，其中ki為彈簧單元的剛度，fi為滑塊單元的屈服力，F(xiàn)和u為外力和位移。令每個Jenkins單元的剛度均為k，位移載荷為x，單元屈服力f以密度函數(shù)r(f) 表示，可以得到Iwan模型表達(dá)式如下：

(1)

將式(1)對x求二階導(dǎo)數(shù)后整理得

(2)

再令φ=f/k,ρ(φ)=k2r(f)，則式(1)可改寫為

(3)

(4)

式中ρ(φ)為密度函數(shù)。Li等[17-18]提出一個含截?cái)鄡缏煞植己碗p脈沖的非均勻密度函數(shù)，如圖2所示，其表達(dá)式為

(5)

將式(5)代入式(3)中積分可以得到六參數(shù)Iwan模型的力-位移曲線關(guān)系，其示意圖見圖3。O-A段力-位移關(guān)系呈線性；A-B段為微觀滑移階段，力-位移關(guān)系呈非線性；從B點(diǎn)開始模型進(jìn)入宏觀滑移階段，力-位移關(guān)系再次呈線性。

將六參數(shù)Iwan模型應(yīng)用于數(shù)值計(jì)算時，將其離散為n+2個并聯(lián)的Jenkins單元，如圖4所示[20]。前n個Jenkins單元J1～Jn描述微觀滑移階段連接結(jié)構(gòu)剛度的變化；第n+1個Jenkins單元Jn+1描述宏觀滑移時刻剛度的變化量，其剛度為K2；第n+2個Jenkins單元Jn+2為描述宏觀滑移后接觸界面殘余剛度現(xiàn)象的彈簧單元，其剛度為K∞。

2 考慮接觸非線性的漸進(jìn)結(jié)構(gòu)優(yōu)化

一個理想的結(jié)構(gòu)，其每一部分的應(yīng)力應(yīng)該接近于相同的安全水平[21]。而受力分析表明結(jié)構(gòu)中的應(yīng)力分布并不均勻，可以認(rèn)為應(yīng)力水平較高的區(qū)域是結(jié)構(gòu)破壞的主要區(qū)域，而在應(yīng)力水平較低的區(qū)域材料處于沒有充分利用的狀態(tài)。以單元的等效應(yīng)力作為刪除判據(jù)，將低應(yīng)力的材料逐漸刪除，更新結(jié)構(gòu)設(shè)計(jì)，使得優(yōu)化后的結(jié)構(gòu)應(yīng)力水平更加均勻，材料構(gòu)成趨于合理，同時也提高了材料的利用率。對于各向同性材料可采用von Mises應(yīng)力作為單元刪除判據(jù)。

(6)

式中Sij為應(yīng)力偏張量，且i=j。

(7)

則認(rèn)為該單元處于低應(yīng)力狀態(tài)，是無效的或低效的單元，對結(jié)構(gòu)整體性能貢獻(xiàn)最小，應(yīng)從結(jié)構(gòu)中刪除。

保持當(dāng)前刪除率不變，重復(fù)刪除單元，直至再無單元滿足式(7)，即優(yōu)化結(jié)果對于當(dāng)前刪除率RR已達(dá)到穩(wěn)定狀態(tài)。此時引入一個進(jìn)化率ER，將刪除率修改為

RR=RR+ER

(8)

以新的刪除率重復(fù)刪除和更新刪除率過程，直到結(jié)構(gòu)滿足給定的約束條件為止。

含螺栓連接結(jié)構(gòu)的拓?fù)鋬?yōu)化的約束條件應(yīng)包括兩個方面：1)結(jié)構(gòu)的強(qiáng)度要求，即結(jié)構(gòu)的應(yīng)力不超過許用應(yīng)力；2)結(jié)構(gòu)的連接剛度要求，即一般不允許螺栓連接界面產(chǎn)生宏觀滑移，造成螺栓受剪切載荷作用。以結(jié)構(gòu)的質(zhì)量最小為目標(biāo)函數(shù)，結(jié)構(gòu)的最大等效應(yīng)力滿足要求、連接不產(chǎn)生宏觀滑移為約束條件，以設(shè)計(jì)區(qū)域材料單元的存在狀態(tài)為設(shè)計(jì)變量，建立如下拓?fù)鋬?yōu)化數(shù)學(xué)模型：

(9)

基于應(yīng)力水平的含螺栓連接結(jié)構(gòu)的非線性漸進(jìn)結(jié)構(gòu)優(yōu)化方法的步驟如下：

1) 建立含螺栓連接的組合結(jié)構(gòu)非線性有限元模型。利用六參數(shù)Iwan模型描述螺栓連接界面的非線性特性；定義載荷及邊界約束條件；

2) 確定優(yōu)化準(zhǔn)則。本文采用應(yīng)力準(zhǔn)則，如式(7)；

3) 定義初始刪除率RR和進(jìn)化率ER；明確結(jié)構(gòu)的許用應(yīng)力和螺栓連接處相對位移(滑移量)允許值；

4) 采用增量法分步施加載荷，對當(dāng)前離散結(jié)構(gòu)進(jìn)行非線性靜力學(xué)求解；依據(jù)當(dāng)前刪除率RR，將滿足式(7)的低應(yīng)力單元組成集合Ωσ，刪除集合中的單元；

5) 保持當(dāng)前刪除率不變，重復(fù)步驟4)，直至再無單元滿足優(yōu)化準(zhǔn)則，即優(yōu)化結(jié)果對前刪除率RR已達(dá)到穩(wěn)定狀態(tài)。此時引入進(jìn)化率ER，按式(8)對刪除率進(jìn)行修改；

6) 以新的刪除率重復(fù)步驟4)～5)，直到單元最大等效應(yīng)力達(dá)到許用應(yīng)力或螺栓連接處的相對位移達(dá)到允許值，優(yōu)化完成。

3 算例

3.1 模型的建立

圖5(a)所示為一含螺栓連接的圓筒組合結(jié)構(gòu)，由圓筒和底座兩個部分組成。底座上半部分為一內(nèi)徑120 mm、外徑370 mm、厚度10 mm的圓盤，其上表面與圓筒連為一體，下半部分為直徑370 mm、厚度10 mm的圓板。三對單螺栓連接組件以120°間隔均勻排布在直徑為300 mm的圓周上，將底座上下部分連接起來。螺栓預(yù)緊力為7 117 N(圖5中未畫出連接螺栓)；圓筒高為335 mm，內(nèi)徑為230mm，外徑為234 mm，厚為2 mm。結(jié)構(gòu)材料為鋼材，彈性模量=200 GPa，泊松比ν=0.3，密度ρ=7 850 kg/m3，許用應(yīng)力σ*=105 MPa。

有限元建模采用SOLID185單元，以六面體網(wǎng)格將結(jié)構(gòu)劃分為31 412個單元，如圖5(b)所示。

文獻(xiàn)[23]對預(yù)緊力為7 117 N時的螺栓連接實(shí)驗(yàn)結(jié)果進(jìn)行了參數(shù)辨識，具體值見表1。將表1的參數(shù)代入六參數(shù)Iwan模型，獲得如圖3所示的單個螺栓連接的力-相對位移關(guān)系(骨線方程)。圖3B點(diǎn)為宏觀滑移的起始點(diǎn)，對應(yīng)的力為4 430 N，相對位移φ2=1.02×10-5m。若取不發(fā)生宏觀滑移的安全系數(shù)[n]=1.2，則允許的最大相對滑移量Smax≤0.85×10-5m。

表1 預(yù)緊力7 117 N的參數(shù)辨識結(jié)果

3.2 拓?fù)鋬?yōu)化

在ANSYS有限元分析中，每對螺栓連接組件之間均以六參數(shù)Iwan模型描述，將每對連接離散為13個COMBIN40單元和1個COMBIN14單元(圖4)，各單元的單元剛度和屈服力的選取參見文獻(xiàn)[23]，APDL命令流如下。

/prep7

!定義單元類型!

et,1,40

et,2,14

!定義單元的剛度和屈服力!

R,1,1.31E+07,,,,2,,

R,2,1.31E+07,,,,8,,

R,3,1.31E+07,,,,15,,

R,4,1.31E+07,,,,24,,

R,5,1.31E+07,,,,35,,

R,6,1.31e+07,,,,46,,

R,7,1.31e+07,,,,57,,

R,8,1.31e+07,,,,70,,

R,9,1.31e+07,,,,83,,

R,10,1.31e+07,,,,97,,

R,11,1.31e+07,,,,111,,

R,12,1.31e+07,,,,126,,

R,13,3.48e8,,,,3549.6,,

R,14,2.065e7,,,,,,

!將1-13號單元設(shè)置為類型1!

*do,rin,1,13

type,1,

real,rin

e,m,n

*enddo

!將14號單元設(shè)置為類型2!

type,2,

real,14,

e,m,n

(說明：命令流中m和n分別為每對螺栓連接的兩個接觸面中心節(jié)點(diǎn)的編號。)

約束底座下端面，對筒體上端面施加7.5 MPa的均勻壓力，并選定初始刪除率RR=0.2%，進(jìn)化率ER=0.1%。利用編寫的命令流輸入文件，完成非線性求解、優(yōu)化循環(huán)及后處理等過程[22]。

優(yōu)化后的圓筒結(jié)構(gòu)如圖6所示，最大等效應(yīng)力為105 MPa(對應(yīng)刪除率為60%)。最終優(yōu)化構(gòu)型中肋條的上部沿圓周均勻分布，承受加載在筒體上端面的均布載荷；下部與螺栓連接結(jié)構(gòu)的排布相對應(yīng)，呈120°分別集中，、向下方的連接部分傳力。很明顯，優(yōu)化后得到拓?fù)錁?gòu)型很好地反映了結(jié)構(gòu)的受力和傳力特征。

表2中列出了隨刪除率的增大，筒體的單元最大von Mises應(yīng)力的變化。從表2可以看出，筒體的最大等效應(yīng)力隨刪除率增加而逐漸增大，當(dāng)刪除率達(dá)60%時，筒體最大等效應(yīng)力達(dá)到許用應(yīng)力，此時連接單元的相對滑移量Smax=8.12×10-6m，小于允許的最大相對滑移量8.5×10-6m，即連接處沒有產(chǎn)生宏觀滑移。

表2 筒體應(yīng)力隨刪除率的變化

4 結(jié)論

本文采用六參數(shù)Iwan模型描述螺栓連接的非線性力學(xué)行為，再基于漸進(jìn)結(jié)構(gòu)優(yōu)化方法和非線性有限元分析，建立含螺栓連接圓筒結(jié)構(gòu)的非線性靜力學(xué)拓?fù)鋬?yōu)化方法。利用ANSYS軟件中的APDL語言編寫命令流輸入文件，實(shí)現(xiàn)優(yōu)化流程。通過含螺栓連接圓筒組合結(jié)構(gòu)的優(yōu)化算例，驗(yàn)證了含界面強(qiáng)非線性的靜力學(xué)拓?fù)鋬?yōu)化方法可以應(yīng)用于含螺栓連接復(fù)雜工程結(jié)構(gòu)的靜力學(xué)拓?fù)鋬?yōu)化設(shè)計(jì)。