采用包絡(luò)-準則法優(yōu)化考慮瞬態(tài)傳熱的薄板結(jié)構(gòu)

孫 國， 閻 琨， 蔡賢輝， 程耿東

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

1 引 言

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

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

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

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

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

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

圖1 薄板Fig.1 Thin plate

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

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

滿足材料體積約束條件：

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 數(shù)值算例分析

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

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

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

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

圖2 四邊簡支板Fig.2 Simply-supported rectangular plate

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

為了校核優(yōu)化設(shè)計，在給定200s的時間段內(nèi)，計算初始模型、第一次優(yōu)化設(shè)計（OPT-1）及第二次優(yōu)化設(shè)計（OPT-2）在瞬態(tài)傳熱過程下的位移和應(yīng)力響應(yīng)歷程。圖5給出了傳熱中間時刻t＝50s時兩種優(yōu)化設(shè)計的應(yīng)力分布，使用相同的標尺對比應(yīng)力。可以看出，OPT-1設(shè)計在溫度荷載作用下，存在高應(yīng)力區(qū)域最大Mises應(yīng)力650MPa；OPT-2設(shè)計有效降低了上述區(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）

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

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

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

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

4 結(jié) 論

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