孫 國(guó)， 閻 琨， 蔡賢輝， 程耿東

（大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室 工程力學(xué)系，大連116024）

1 引 言

在工程結(jié)構(gòu)中，薄板結(jié)構(gòu)廣泛應(yīng)用于航空航天、機(jī)械和交通等領(lǐng)域，往往起著承載和隔斷空間的作用，因此，彈性薄板結(jié)構(gòu)分析及最優(yōu)設(shè)計(jì)受到持續(xù)的關(guān)注。程耿東等［1－3］以薄板的柔順性為目標(biāo)函數(shù)，由變分原理出發(fā)，推導(dǎo)出優(yōu)化準(zhǔn)則，構(gòu)造了求解板的優(yōu)化厚度分布函數(shù)的迭代格式，以提高薄板結(jié)構(gòu)的剛度。

在重要的工業(yè)裝備結(jié)構(gòu)中，受到高溫的熱結(jié)構(gòu)失效往往會(huì)造成嚴(yán)重后果，其設(shè)計(jì)成為核心技術(shù)。即使是穩(wěn)定的溫度場(chǎng)，熱荷載作用于結(jié)構(gòu)時(shí)，一方面產(chǎn)生溫度應(yīng)力和溫度變形，另一方面降低材料的很多力學(xué)性能，給結(jié)構(gòu)設(shè)計(jì)帶來(lái)很多困難，某些條件下增加材料用量使得結(jié)構(gòu)內(nèi)的熱應(yīng)力更大，而這和沒(méi)有溫度荷載的情況正好相反。熱結(jié)構(gòu)設(shè)計(jì)面臨的這些特殊問(wèn)題需要研究特定的設(shè)計(jì)方法。Rao等［4］基于溫度邊界用最優(yōu)準(zhǔn)則法設(shè)計(jì)，Adelman等［5］進(jìn)行了熱彈性的滿應(yīng)力設(shè)計(jì)，Haftka等［6］通過(guò)近似方法進(jìn)行了耦合熱結(jié)構(gòu)設(shè)計(jì)，Tortorelli等［7］分析了熱彈性問(wèn)題的結(jié)構(gòu)響應(yīng)靈敏度等。近年多為以控制結(jié)構(gòu)應(yīng)力和位移為目標(biāo)的拓?fù)鋬?yōu)化研究。左孔天等［8］采用多目標(biāo)優(yōu)化方法進(jìn)行傳熱結(jié)構(gòu)的拓?fù)鋬?yōu)化設(shè)計(jì)。Pedersen等［9－10］在結(jié)構(gòu)柔順性優(yōu)化之后，增加一個(gè)啟發(fā)式的應(yīng)變能密度的均勻化過(guò)程，以處理局部應(yīng)力過(guò)高的問(wèn)題。張衛(wèi)紅等［11］用柔順性和彈性應(yīng)變能兩種指標(biāo)進(jìn)行熱結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)，給出數(shù)值算例并指出應(yīng)變能指標(biāo)對(duì)于降低應(yīng)力水平更為有效。倪曉琴等［12］推導(dǎo)了實(shí)心彈性薄板結(jié)構(gòu)在溫度荷載作用下的優(yōu)化準(zhǔn)則，有效地降低了板殼結(jié)構(gòu)在穩(wěn)態(tài)溫度場(chǎng)下的熱變形。

需要考慮瞬態(tài)傳熱過(guò)程的薄板結(jié)構(gòu)優(yōu)化設(shè)計(jì)更加困難，因?yàn)闇囟群奢d不僅與時(shí)間有關(guān)，也與模型厚度相關(guān)，而在優(yōu)化過(guò)程中厚度分布會(huì)不斷修改，意味著不再滿足前期一些優(yōu)化問(wèn)題假定的均勻穩(wěn)態(tài)溫度場(chǎng)分布，更為嚴(yán)重的是，如果疊加在結(jié)構(gòu)上的溫度應(yīng)力變化范圍很大，穩(wěn)態(tài)優(yōu)化設(shè)計(jì)結(jié)果在瞬態(tài)傳熱過(guò)程中可能出現(xiàn)較大的局部應(yīng)力，需要建立非穩(wěn)態(tài)溫度場(chǎng)下的優(yōu)化方法。Kok等［13］利用響應(yīng)面方法進(jìn)行瞬態(tài)熱應(yīng)力問(wèn)題的優(yōu)化求解；Turteltaub［14］提出了瞬態(tài)傳熱下對(duì)材料特性進(jìn)行拓?fù)鋬?yōu)化；Molaei Najafabadi等［15］進(jìn)行了考慮瞬態(tài)傳熱過(guò)程的功能梯度材料的參數(shù)優(yōu)化設(shè)計(jì)。本文對(duì)已有的穩(wěn)態(tài)溫度場(chǎng)下薄板的優(yōu)化準(zhǔn)則加以局部修正，以處理瞬態(tài)傳熱過(guò)程中的薄板結(jié)構(gòu)優(yōu)化問(wèn)題。

2 瞬態(tài)溫度場(chǎng)下薄板結(jié)構(gòu)優(yōu)化設(shè)計(jì)

2.1 只有外力作用下彈性薄板的最小柔順性設(shè)計(jì)

圖1為平面薄板結(jié)構(gòu)，在平板中面上建立坐標(biāo)面xOy，z軸沿中面法向，板厚h＝h（x，y）。

圖1 薄板Fig.1 Thin plate

如果考慮薄板承受垂直于中面的機(jī)械荷載q的作用，通常將目標(biāo)函數(shù)選取為結(jié)構(gòu)的柔順性：

在單純外力作用且無(wú)初始應(yīng)力的條件下，式（1）給出的柔順性正比于結(jié)構(gòu)彈性應(yīng)變能，如式（2）所示。

滿足材料體積約束條件：

設(shè)計(jì)變量h（x，y）除了滿足積分約束，還有最大最小限制hmin≤h（x，y）≤hmax。采用變分法可以得到該優(yōu)化問(wèn)題的最優(yōu)條件如式（4）所示，具體推導(dǎo)不再贅述［1］。

式中 Mx，My和Mxy為在截面單位寬度上的彎矩及扭矩，ν是泊松比，λ＊是拉格朗日乘子。

引入指標(biāo)函數(shù)g（x，y）：

求解關(guān)于厚度h（x，y）的方程，優(yōu)化條件可重新寫(xiě)為

λ＊為優(yōu)化問(wèn)題的拉格朗日乘子，式（6）不能直接求出h（x，y），可以采用迭代形式求解。根據(jù)體積約束條件，可以解出拉格朗日乘子λ＊。

式中V為材料總體積，Ωu，Ωmin和Ωmax分別代表厚度取中間值、最小值及最大值的區(qū)域。

薄板彎曲最大應(yīng)力出現(xiàn)在上下表面處，以上表面為例，

指標(biāo)g（x，y）可以由表面最大應(yīng)力表示，該值和板內(nèi)各點(diǎn)的應(yīng)變能密度有關(guān)，這一特點(diǎn)使得將厚度修改式（6）應(yīng)用于受到機(jī)械力和穩(wěn)態(tài)溫度荷載聯(lián)合作用的薄板的柔順性優(yōu)化設(shè)計(jì)中也很有效。

2.2 瞬態(tài)溫度場(chǎng)下薄板結(jié)構(gòu)優(yōu)化設(shè)計(jì)

當(dāng)考慮瞬態(tài)熱-結(jié)構(gòu)耦合問(wèn)題時(shí)，假定結(jié)構(gòu)受到的外力和熱荷載是一個(gè)典型的準(zhǔn)穩(wěn)態(tài)過(guò)程，溫度場(chǎng)按導(dǎo)熱方程瞬態(tài)求解，應(yīng)力場(chǎng)取不同時(shí)刻溫度場(chǎng)，按穩(wěn)態(tài)求解，忽略結(jié)構(gòu)變形對(duì)溫度場(chǎng)的影響。

從優(yōu)化設(shè)計(jì)的角度來(lái)看，瞬態(tài)溫度場(chǎng)T（x，y，z，t）的變化特性使得Mx，My和Mxy同樣是隨時(shí)間變化的，則穩(wěn)態(tài)溫度場(chǎng)下的優(yōu)化準(zhǔn)則難以適用，直接建立非穩(wěn)態(tài)溫度場(chǎng)下的優(yōu)化準(zhǔn)則同樣十分困難。為了處理這種復(fù)雜分布的T（x，y，z，t），將一個(gè)考慮瞬態(tài)傳熱復(fù)雜溫度荷載下的彈性薄板優(yōu)化問(wèn)題分成兩個(gè)子問(wèn)題，即控制撓度w的優(yōu)化問(wèn)題和控制局部熱應(yīng)力較高的問(wèn)題。

基于已有的薄板優(yōu)化有關(guān)研究工作，在單獨(dú)外力q作用下，基于柔順性的優(yōu)化結(jié)果可以提高結(jié)構(gòu)抗彎剛度，進(jìn)而有效降低機(jī)械荷載q作用下的撓度w，溫度荷載對(duì)于撓度w的影響主要源于溫度彎矩MT的影響（面內(nèi)熱膨脹對(duì)于撓度w的影響可以忽略），這個(gè)基于柔順性的優(yōu)化結(jié)果仍然可以控制q和MT共同作用下的撓度。

為解決傳熱過(guò)程中出現(xiàn)的局部熱應(yīng)力較高的問(wèn)題，本文給出一種包絡(luò)-準(zhǔn)則方法，對(duì)瞬態(tài)傳熱過(guò)程中高應(yīng)力區(qū)域進(jìn)行局部修改處理。考慮構(gòu)造表征傳熱過(guò)程中多個(gè)時(shí)刻下g（x，y）最大值的包絡(luò)函數(shù)G（x，y）和局部修正的權(quán)重函數(shù)W（x，y）。構(gòu)造包絡(luò)函數(shù)時(shí)間點(diǎn)的具體取值需要兼顧計(jì)算效率與包絡(luò)函數(shù)的精度。針對(duì)一種常見(jiàn)的薄板結(jié)構(gòu)單側(cè)對(duì)流傳熱的瞬態(tài)過(guò)程，為了避免在后期接近熱平衡的時(shí)間段內(nèi)配置過(guò)多時(shí)間點(diǎn)，經(jīng)驗(yàn)性地采用逐漸增大間隔時(shí)間點(diǎn)計(jì)算包絡(luò)函數(shù)?？紤]到具體選點(diǎn)方案對(duì)包絡(luò)函數(shù)計(jì)算結(jié)果的精度是有影響的，針對(duì)初始設(shè)計(jì)模型，將擬采用的時(shí)間選點(diǎn)方案與相對(duì)密集選點(diǎn)方案計(jì)算結(jié)果相比較，以保證能有足夠的精度。為實(shí)現(xiàn)局部修正，對(duì)于高應(yīng)力區(qū)域，令W（x，y）近似等于1，其他位置則相對(duì)較小，修正后的指標(biāo)量G（x，y）W（x，y）相當(dāng)于在局部高應(yīng)力區(qū)域采用了一個(gè)類似的指標(biāo)量G（x，y），而其他應(yīng)力較低的區(qū)域則不做修正。經(jīng)過(guò)一定次數(shù)迭代，局部應(yīng)力集中現(xiàn)象逐漸得到緩解。另一方面，進(jìn)一步的迭代會(huì)使修正范圍增大，由于包括了瞬態(tài)過(guò)程多個(gè)時(shí)刻的影響，迭代過(guò)程的收斂性和有效性不能嚴(yán)格保證，較大的修正范圍可能導(dǎo)致設(shè)計(jì)結(jié)果出現(xiàn)劣化，因此這個(gè)局部修正迭代過(guò)程需要檢查最大應(yīng)力變化，并采用最優(yōu)設(shè)計(jì)結(jié)果。

該包絡(luò)-準(zhǔn)則方法包括兩個(gè)環(huán)節(jié)的優(yōu)化設(shè)計(jì)。首先，針對(duì)外力荷載q，進(jìn)行一個(gè)結(jié)構(gòu)柔順性的優(yōu)化設(shè)計(jì)；以這一設(shè)計(jì)為基礎(chǔ)，通過(guò)計(jì)算多個(gè)時(shí)刻指標(biāo)函數(shù)g并取其包絡(luò)，對(duì)其進(jìn)行修正迭代，以消除瞬態(tài)溫度場(chǎng)作用下的局部高應(yīng)力。其具體步驟和有關(guān)公式如下。

（1）首先利用式（5，6）進(jìn)行外力q作用下的優(yōu)化設(shè)計(jì)，為方便表述，將這個(gè)優(yōu)化環(huán)節(jié)稱為OPT-1，其結(jié)果作為進(jìn)一步考慮瞬態(tài)傳熱優(yōu)化設(shè)計(jì)的初值。為了與溫度荷載下的計(jì)算結(jié)果區(qū)分，將式（5）改寫(xiě)為下標(biāo)q的形式。

（2）為了在設(shè)計(jì)中考慮瞬態(tài)傳熱過(guò)程的影響，定義瞬態(tài)傳熱過(guò)程t時(shí)刻的指標(biāo)函數(shù)gqT（x，y，t）值為

式中 下標(biāo)qT代表載荷q和溫度T聯(lián)合作用，取gqT在多個(gè)時(shí)刻ti的值，并構(gòu)造包絡(luò)函數(shù)GqT為

式中ti（i＝1，2，…，n）代表瞬態(tài)不同時(shí)刻。

（3）構(gòu)造局部修正的權(quán)重函數(shù)W（x，y）

利用GqT計(jì)算對(duì)應(yīng)的應(yīng)變能密度SqT，

將SqT最大值歸一化，構(gòu)造表征應(yīng)變能密度水平的權(quán)重函數(shù)W（x，y），

（4）通過(guò)對(duì)比兩種狀態(tài)下的指標(biāo)函數(shù)gq（x，y）和G（x，y）W（x，y），實(shí)現(xiàn)瞬態(tài)過(guò)程中高應(yīng)力區(qū)的局部修正，修正后的指標(biāo)變量為

通過(guò)與式（6，7）類似的迭代計(jì)算，求解厚度h（x，y）。為方便表述，將這個(gè)優(yōu)化環(huán)節(jié)稱為OPT-2。

3 數(shù)值算例分析

如上文所述，包絡(luò)-準(zhǔn)則法優(yōu)化包括OPT-1和OPT-2兩個(gè)環(huán)節(jié)，OPT-2采用OPT-1設(shè)計(jì)結(jié)果作為初值。算例的主優(yōu)化程序通過(guò)Matlab編程實(shí)現(xiàn)，調(diào)用ANSYS軟件執(zhí)行有限元模型修改和分析任務(wù)并輸出計(jì)算結(jié)果數(shù)據(jù)文件，Matlab主優(yōu)化程序進(jìn)而讀取這些數(shù)據(jù)文件，執(zhí)行OPT-1或OPT-2優(yōu)化計(jì)算，輸出更新后的厚度數(shù)據(jù)文件，以供下一次迭代計(jì)算采用。

在OPT-2涉及的瞬態(tài)傳熱分析中，必須考慮沿厚度方向瞬態(tài)溫度變化的影響，瞬態(tài)溫度場(chǎng)分析采用在平面和厚度方向上具有傳導(dǎo)熱能力的層狀殼單元Shell-131，為4節(jié)點(diǎn)單元，該單元適用于殼單元的穩(wěn)態(tài)或瞬態(tài)熱分析，計(jì)算生成的溫度可以傳遞至殼單元進(jìn)行熱應(yīng)力求解。應(yīng)力分析采用Shell-181單元，為4節(jié)點(diǎn)單元，每個(gè)結(jié)點(diǎn)具有6個(gè)自由度，適合對(duì)薄殼到具有一定厚度的殼體結(jié)構(gòu)進(jìn)行分析。上述兩種單元?dú)んw厚度和積分點(diǎn)數(shù)量的設(shè)定可以在橫截面的定義中完成。本文將shdll-131單元和shell-181單元都設(shè)定為3層，每層積分點(diǎn)數(shù)為3個(gè)，單元沿厚度共有7個(gè)積分點(diǎn)（因?qū)娱g兩處積分點(diǎn)重合），較多的積分點(diǎn)有助于保證分層計(jì)算的精度，但同時(shí)導(dǎo)致計(jì)算量的增加。將4個(gè)層面上的溫度計(jì)算結(jié)果傳遞到shell-181單元，構(gòu)成一個(gè)沿厚度的3層線性溫度函數(shù)。

采用的材料特性為，彈性模量E＝200GPa，泊松比ν＝0.31，熱膨脹系數(shù)α＝1.3E-5，導(dǎo)熱系數(shù)k＝15.0W／m·K，比熱容c＝450J／kgK，密度ρ＝8900kg／m3。

如圖2所示為四邊簡(jiǎn)支板，a＝b＝1.0m，根據(jù)對(duì)稱性特征，取其1／4進(jìn)行分析，均勻劃分為2500個(gè)單元。給定下表面壓力P＝100kPa，下表面對(duì)流換熱系數(shù)hg＝250W／m2·K，環(huán)境溫度TW＝500℃，外側(cè)表面為絕熱邊界條件。結(jié)構(gòu)初始溫度T0＝0℃，壁板結(jié)構(gòu)初始設(shè)計(jì)的厚度h0＝8mm，設(shè)計(jì)變量上下限為hmin＝2mm和hmax＝20mm?？紤]到初始階段溫度變化較快，同時(shí)避免過(guò)大的計(jì)算負(fù)擔(dān)，采用時(shí)間點(diǎn)取方案ti＝｛0，10，20，50，100，150，200｝（s），另采用較密集方案ti＝｛0，10，20，…，200｝（s）作為參照進(jìn)行對(duì)比，G（x，y）計(jì)算結(jié)果分別記作單元向量形式和，以范數(shù)之比表征兩者相對(duì)誤差，‖G1－G0‖2／‖G0‖2＝0.66%，‖G1－G0‖∞／‖G0‖∞＝1.3%，可見(jiàn)兩種選點(diǎn)方案計(jì)算結(jié)果較為接近。

圖2 四邊簡(jiǎn)支板Fig.2 Simply-supported rectangular plate

經(jīng)過(guò)20次迭代后，OPT-2過(guò)程最大應(yīng)力迭代變化如圖3所示，在第10次迭代得到最優(yōu)設(shè)計(jì)結(jié)果。兩次優(yōu)化厚度分布如圖4所示，可以看出兩個(gè)優(yōu)化設(shè)計(jì)在整體布局上具有相似性。

為了校核優(yōu)化設(shè)計(jì)，在給定200s的時(shí)間段內(nèi)，計(jì)算初始模型、第一次優(yōu)化設(shè)計(jì)（OPT-1）及第二次優(yōu)化設(shè)計(jì)（OPT-2）在瞬態(tài)傳熱過(guò)程下的位移和應(yīng)力響應(yīng)歷程。圖5給出了傳熱中間時(shí)刻t＝50s時(shí)兩種優(yōu)化設(shè)計(jì)的應(yīng)力分布，使用相同的標(biāo)尺對(duì)比應(yīng)力?？梢钥闯觯琌PT-1設(shè)計(jì)在溫度荷載作用下，存在高應(yīng)力區(qū)域最大Mises應(yīng)力650MPa；OPT-2設(shè)計(jì)有效降低了上述區(qū)域的應(yīng)力幅值，最大Mises應(yīng)力降為282.5MPa。

圖3 Mises應(yīng)力最大值迭代歷程（OPT-2）Fig.3 Maximum Mises stress via iterations（OPT-2）

圖4 優(yōu)化后厚度分布（單位：m）Fig.4 Optimized thickness distribution（unit：m）

在考慮的時(shí)間段內(nèi)，結(jié)構(gòu)中心位置的z向撓度w如圖6所示，等效應(yīng)力最大值如圖7所示。從圖6可以看出，在考慮的時(shí)間段內(nèi)，四邊簡(jiǎn)支板中心位置撓度w隨時(shí)間發(fā)生變化，相對(duì)于原始均勻厚度模型中心位置撓度值44.1mm，OPT-1和 OPT-2設(shè)計(jì)都具有較小的中心撓度w，分別為10.49mm和13.81mm。從圖7可以看出，在所考慮的200s時(shí)間段內(nèi)，最大Mises應(yīng)力同樣呈較強(qiáng)的隨時(shí)間變化特性，原始均勻厚度模型的Mises應(yīng)力最大值為535MPa，由于均勻厚度，其變化相對(duì)平緩，OPT-1對(duì)應(yīng)的 Mises應(yīng)力最大值為683MPa，OPT-2對(duì)應(yīng)的 Mises應(yīng)力最大值為330.4MPa，可見(jiàn)OPT-2相較于OPT-1，200s瞬態(tài)傳熱過(guò)程中的熱應(yīng)力幅值最大值下降了51.6%。

圖5 t＝50s時(shí)Mises應(yīng)力（單位：MPa）Fig.5 Mises stress at time t＝50s（unit：MPa）

圖6 不同時(shí)刻中心位置撓度wFig.6 Deflection of center position wvia time

圖7 不同時(shí)刻Mises應(yīng)力最大值Fig.7 Maximum Mises stress via time

4 結(jié) 論

針對(duì)瞬態(tài)傳熱過(guò)程中存在的溫度場(chǎng)不均衡以及局部應(yīng)力水平較高問(wèn)題，采用包絡(luò)-準(zhǔn)則優(yōu)化方法進(jìn)行處理。第一次優(yōu)化設(shè)計(jì)（OPT-1）使橫向位移得到控制，但是應(yīng)力水平增高較大；經(jīng)第二次優(yōu)化設(shè)計(jì)（OPT-2），相對(duì)于第一次優(yōu)化，位移增量相對(duì)較小，而應(yīng)力水平得到有效降低。經(jīng)采用的包絡(luò)-準(zhǔn)則優(yōu)化方法優(yōu)化后，在控制撓度的同時(shí)，也使得結(jié)構(gòu)局部較高的應(yīng)力水平明顯降低，結(jié)構(gòu)剛度和強(qiáng)度特性都得到明顯改善，顯示了處理此類問(wèn)題的有效性。