采用包絡-準則法優化考慮瞬態傳熱的薄板結構

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

（大連理工大學 工業裝備結構分析國家重點實驗室 工程力學系，大連116024）

1 引 言

在工程結構中，薄板結構廣泛應用于航空航天、機械和交通等領域，往往起著承載和隔斷空間的作用，因此，彈性薄板結構分析及最優設計受到持續的關注。程耿東等［1－3］以薄板的柔順性為目標函數，由變分原理出發，推導出優化準則，構造了求解板的優化厚度分布函數的迭代格式，以提高薄板結構的剛度。

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

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

2 瞬態溫度場下薄板結構優化設計

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

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

圖1 薄板Fig.1 Thin plate

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

在單純外力作用且無初始應力的條件下，式（1）給出的柔順性正比于結構彈性應變能，如式（2）所示。

滿足材料體積約束條件：

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

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

引入指標函數g（x，y）：

求解關于厚度h（x，y）的方程，優化條件可重新寫為

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

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

薄板彎曲最大應力出現在上下表面處，以上表面為例，

指標g（x，y）可以由表面最大應力表示，該值和板內各點的應變能密度有關，這一特點使得將厚度修改式（6）應用于受到機械力和穩態溫度荷載聯合作用的薄板的柔順性優化設計中也很有效。

2.2 瞬態溫度場下薄板結構優化設計

當考慮瞬態熱-結構耦合問題時，假定結構受到的外力和熱荷載是一個典型的準穩態過程，溫度場按導熱方程瞬態求解，應力場取不同時刻溫度場，按穩態求解，忽略結構變形對溫度場的影響。

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

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

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

該包絡-準則方法包括兩個環節的優化設計。首先，針對外力荷載q，進行一個結構柔順性的優化設計；以這一設計為基礎，通過計算多個時刻指標函數g并取其包絡，對其進行修正迭代，以消除瞬態溫度場作用下的局部高應力。其具體步驟和有關公式如下。

（1）首先利用式（5，6）進行外力q作用下的優化設計，為方便表述，將這個優化環節稱為OPT-1，其結果作為進一步考慮瞬態傳熱優化設計的初值。為了與溫度荷載下的計算結果區分，將式（5）改寫為下標q的形式。

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

式中 下標qT代表載荷q和溫度T聯合作用，取gqT在多個時刻ti的值，并構造包絡函數GqT為

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

（3）構造局部修正的權重函數W（x，y）

利用GqT計算對應的應變能密度SqT，

將SqT最大值歸一化，構造表征應變能密度水平的權重函數W（x，y），

（4）通過對比兩種狀態下的指標函數gq（x，y）和G（x，y）W（x，y），實現瞬態過程中高應力區的局部修正，修正后的指標變量為

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

3 數值算例分析

如上文所述，包絡-準則法優化包括OPT-1和OPT-2兩個環節，OPT-2采用OPT-1設計結果作為初值。算例的主優化程序通過Matlab編程實現，調用ANSYS軟件執行有限元模型修改和分析任務并輸出計算結果數據文件，Matlab主優化程序進而讀取這些數據文件，執行OPT-1或OPT-2優化計算，輸出更新后的厚度數據文件，以供下一次迭代計算采用。

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

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

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

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

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

為了校核優化設計，在給定200s的時間段內，計算初始模型、第一次優化設計（OPT-1）及第二次優化設計（OPT-2）在瞬態傳熱過程下的位移和應力響應歷程。圖5給出了傳熱中間時刻t＝50s時兩種優化設計的應力分布，使用相同的標尺對比應力。可以看出，OPT-1設計在溫度荷載作用下，存在高應力區域最大Mises應力650MPa；OPT-2設計有效降低了上述區域的應力幅值，最大Mises應力降為282.5MPa。

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

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

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

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

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

圖7 不同時刻Mises應力最大值Fig.7 Maximum Mises stress via time

4 結 論

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