流熱力耦合的高性能結(jié)構(gòu)拓撲優(yōu)化設(shè)計方法

摘要：拓撲優(yōu)化和增材制造技術(shù)的快速發(fā)展為高性能復(fù)雜裝備提供了高效的產(chǎn)品設(shè)計和制造方法。目前高性能結(jié)構(gòu)拓撲優(yōu)化只考慮熱力耦合或者流熱耦合的設(shè)計，且大多局限于簡單結(jié)構(gòu)，未考慮流熱力三場共同作用下的設(shè)計，限制了結(jié)構(gòu)性能的提升。針對流熱力多物理場工況下的高性能復(fù)雜結(jié)構(gòu)設(shè)計這一挑戰(zhàn)，提出了一種流熱力耦合拓撲優(yōu)化方法，以提高復(fù)雜結(jié)構(gòu)的承溫能力。首先引入流場、溫度場和結(jié)構(gòu)位移場的控制方程，對計算域的流固材料進行統(tǒng)一表征；然后以最小化平均溫度為目標(biāo)，以流動能量耗散和結(jié)構(gòu)柔度為約束，建立了流熱力耦合的拓撲優(yōu)化模型，并結(jié)合變分法和拉格朗日函數(shù)開展了設(shè)計變量的靈敏度分析；最后將所建立的拓撲優(yōu)化模型應(yīng)用于渦輪的結(jié)構(gòu)設(shè)計，得到了散熱性能良好、流道分布合理的可增材制造結(jié)構(gòu)。

關(guān)鍵詞：拓撲優(yōu)化；變密度法；多物理場；高性能結(jié)構(gòu)

中圖分類號：TH122；TP391.7

DOI：10.3969/j.issn.1004132X.2024.03.011

0引言

高性能結(jié)構(gòu)設(shè)計一直以來都是非常具有挑戰(zhàn)性的工作。作為一種先進的結(jié)構(gòu)設(shè)計方法，拓撲優(yōu)化（topologyoptimization，TO）因其無限的設(shè)計自由度和正向設(shè)計的理念而成為高性能結(jié)構(gòu)的理想設(shè)計工具［1-2］。近年來，增材制造（additivemanufacturing，AM）因其可以實現(xiàn)傳統(tǒng)工藝難以生產(chǎn)的復(fù)雜幾何結(jié)構(gòu)而引起了廣泛關(guān)注［3-4］。相比于車削、鍛造等傳統(tǒng)減材加工工藝，增材制造通過逐層打印的方式完成由計算機輔助設(shè)計（computeraideddesign，CAD）模型到三維結(jié)構(gòu)的制造，極大地縮短了產(chǎn)品開發(fā)周期，減少了材料成本。拓撲優(yōu)化與增材制造的結(jié)合為高性能復(fù)雜裝備提供了更高效的產(chǎn)品設(shè)計和制造方法，具有顯著的研究意義和社會價值［5］。

拓撲優(yōu)化最初是針對結(jié)構(gòu)力學(xué)問題而開發(fā)的，現(xiàn)已廣泛應(yīng)用于由單一物理場控制的一系列工程問題，例如最小化結(jié)構(gòu)柔度設(shè)計［6］、最佳熱傳導(dǎo)結(jié)構(gòu)的設(shè)計［7-8］、流體流動設(shè)計［9］、波傳播、聲學(xué)設(shè)計［10-11］和生物醫(yī)學(xué)設(shè)計［12-13］等。在結(jié)構(gòu)柔度設(shè)計方面，F(xiàn)ERRARI等［14］給出了針對結(jié)構(gòu)柔度拓撲優(yōu)化的緊湊型MATLAB代碼，討論了加速策略的使用，大大減少了總體計算時間。WU等［15］提出一種針對shell-infill結(jié)構(gòu)的柔度拓撲優(yōu)化方法，該方法能夠生成同時具有外部殼體和內(nèi)部填充材料的多孔結(jié)構(gòu)。WANG等［16］結(jié)合柔度比例濾波器提出了一種新的密度變量更新方案，避免了獲取靈敏度信息時涉及的數(shù)值推導(dǎo)和計算問題，提高了收斂速度。KRIEGESMANN等［17］考慮載荷和空間變化的彈性模量的不確定性，提出了一種穩(wěn)健的拓撲優(yōu)化方法，該方法利用一階二階矩（firstordersecondmoment，F(xiàn)OSM）方法來估計柔度的平均值和方差，降低了計算成本。在熱力耦合結(jié)構(gòu)設(shè)計方面，RODRIGUES等［18］首先使用均勻化方法建立了最小化結(jié)構(gòu)柔度的熱彈性拓撲優(yōu)化模型。YANG等［19］研究了機械力和熱負載共同作用下的動態(tài)柔度拓撲優(yōu)化。高效散熱裝置［20］、特定溫度條件下的最大剛度結(jié)構(gòu)［21］以及受溫度約束的結(jié)構(gòu)［22］等設(shè)計先后被提出。TAKEZAWA等［23］設(shè)計了強度和熱傳導(dǎo)約束下的熱彈性結(jié)構(gòu)。DEATON等［24］提出了一種基于應(yīng)力的p范數(shù)聚合函數(shù)的熱結(jié)構(gòu)設(shè)計。ZHENG等［25］考慮到材料屬性、載荷等相關(guān)的混合不確定性，實現(xiàn)了穩(wěn)健性熱力耦合拓撲優(yōu)化設(shè)計。

受限于流固耦合邊界的識別與表征以及流動問題的大規(guī)模計算量，涉及耦合三維流場的數(shù)值計算非常具有挑戰(zhàn)性。在流固耦合拓撲優(yōu)化設(shè)計方面，YOON［26］研究了考慮流固邊界變化的拓撲優(yōu)化問題，將結(jié)構(gòu)平衡方程和流體N-S方程同時擴展到了流固域，避免了邊界處數(shù)據(jù)的傳輸及迭代求解。PICELLI等［27］提出了一種針對流固耦合系統(tǒng)頻率響應(yīng)的拓撲優(yōu)化方法，該方法可應(yīng)用于汽車和飛機乘客艙的降噪以及水下結(jié)構(gòu)的振動控制等。LI等［28］使用反應(yīng)擴散方程來更新水平集函數(shù)，并利用貼體自適應(yīng)網(wǎng)格進行流固邊界捕捉，研究了三維流固耦合問題。在流熱耦合設(shè)計方面，YOON等［29］在N-S方程中引入了摩擦力項，并忽略了流體內(nèi)部的浮力和粘黏耗散，以最小化熱順應(yīng)性為目標(biāo)設(shè)計了利用強制對流傳熱的二維散熱結(jié)構(gòu)。YU等［30］和YAJI等［31］分別使用移動組件法（movingmorphablecomponent，MMC）和水平集法研究了流動傳熱拓撲優(yōu)化問題。ZHAO等［32］研究了考慮傳熱和結(jié)構(gòu)承載的冷卻直管二維截面設(shè)計問題，并提出了一種由固體、液體和孔洞組成的多相材料拓撲優(yōu)化方法。

上述高性能結(jié)構(gòu)拓撲優(yōu)化工作只考慮熱力耦合、流固耦合及流熱耦合的設(shè)計，且大多研究都局限于簡單幾何模型下的小規(guī)模問題，未考慮流熱力三場共同作用下的設(shè)計。然而，工程上以渦輪機械為代表的很多復(fù)雜裝備處于流動傳熱結(jié)構(gòu)應(yīng)力共同作用下的工況，因此，為實現(xiàn)此類復(fù)雜裝備的高性能結(jié)構(gòu)設(shè)計，研究流熱力耦合的拓撲優(yōu)化更具有實際意義。本文基于變密度法，提出一種高性能結(jié)構(gòu)流熱力耦合拓撲優(yōu)化方法，以提升復(fù)雜結(jié)構(gòu)的性能。首先采用Brinkman懲罰方式對流體域和固體域進行建模，引入流場、溫度場和結(jié)構(gòu)位移場的控制方程。其次以平均溫度為目標(biāo)，流動能量耗散和結(jié)構(gòu)柔度為約束，構(gòu)建了流熱力耦合拓撲優(yōu)化模型。為實現(xiàn)單元設(shè)計變量的靈敏度計算，通過引入弱形式的拉格朗日函數(shù)和變分法推導(dǎo)了連續(xù)形式的伴隨方程。使用基于有限體積法離散的OpenFOAM軟件對控制方程和伴隨方程進行求解，拓撲優(yōu)化問題的求解策略為移動漸近線法（methodofmovingasymptotes，MMA），并利用兩種數(shù)值過濾算法穩(wěn)定優(yōu)化歷程。最后將建立的拓撲優(yōu)化模型應(yīng)用于渦輪的結(jié)構(gòu)設(shè)計，得到了散熱性能良好、流道分布合理的可增材制造結(jié)構(gòu)。需要強調(diào)的是，本文的研究著重于復(fù)雜結(jié)構(gòu)流熱力耦合拓撲優(yōu)化算法的理論及其實現(xiàn)，具體的增材制造過程及可制造性等方面的內(nèi)容不在本文的研究范疇。

1流熱力多物理場控制方程

考慮圖1所示的多物理場問題模型，整體計算域Ω由固體域Ωs和流體域Ωf組成，計算域邊界包含流體入口及出口邊界Γin、Γout，固定溫度邊界ΓT0，熱通量邊界Γq，固定位移邊界ΓD0，應(yīng)力邊界Γf，其余邊界為應(yīng)力自由的絕熱邊界。該問題中的各物理場分別由N-S方程、能量方程和線彈性平衡方程進行描述，各物理場之間的耦合情況如下：

（1）流熱耦合。主要考慮流場對溫度場的影響，即流體的速度變化通過對流和擴散過程影響溫度場的分布，忽略溫度場對流速場的影響。

（2）流力（流固）耦合。主要考慮設(shè)計域中固體材料部分對流速場的阻礙作用，忽略流速場引起的結(jié)構(gòu)變形。

（3）流熱力耦合。主要考慮彈性變形引起的結(jié)構(gòu)應(yīng)力，忽略溫度變化產(chǎn)生的熱應(yīng)力。流固材料屬性均為各向同性，且不隨溫度發(fā)生變化。

1.1流動控制方程

本文將流體視為穩(wěn)態(tài)不可壓縮充分發(fā)展的層流狀態(tài)，流動控制方程由N-S方程、連續(xù)性方程及相應(yīng)邊界條件表示：

1.2能量方程

計算域的溫度場由能量方程進行控制，本文采用穩(wěn)態(tài)熱傳導(dǎo)方程作為能量控制方程，能量方程通過對流項與N-S方程進行弱耦合。能量方程及其邊界條件定義如下：

1.3線彈性平衡方程

計算域固體材料的位移場由線彈性平衡方程進行控制，張量形式的線彈性平衡可表示為

其中，D代表結(jié)構(gòu)位移場；fb為固體體積力項；I為單位矩陣；在固定位移邊界上指定位移D0；在應(yīng)力邊界上指定面力fN；其余邊界為自由邊界。μ和λ為材料的拉梅常數(shù)，它們與彈性模量E和泊松比ν的關(guān)系分別為

2流熱力拓撲優(yōu)化模型

2.1優(yōu)化模型

本文使用基于密度法的框架進行結(jié)構(gòu)的流熱力耦合拓撲優(yōu)化設(shè)計，其基本思想是在設(shè)計域Ω內(nèi)用設(shè)計變量場γ來描述域內(nèi)流固兩相材料的分布：γ=1代表流體材料，γ=0代表固體材料。在優(yōu)化過程中，兩相材料的物理屬性由標(biāo)量場γ插值描述，從而將設(shè)計變量引入多物理場控制方程。

散熱性能和降溫效果通常是影響受熱結(jié)構(gòu)綜合性能的重要因素，是本文優(yōu)化的主要目標(biāo)。本文以結(jié)構(gòu)的平均溫度作為目標(biāo)函數(shù)來衡量降溫散熱效果，平均溫度目標(biāo)函數(shù)Ψ的定義如下：

Ψ=∫ΩTdΩ/∫ΩdΩ（6）

一般情況下，作為冷卻介質(zhì)的流體通過冷卻通道的壓力損失越大表明結(jié)構(gòu)的散熱性能越好，但實際上流體的流動損失受到流動效率要求的限制。此外，在涉及流場的優(yōu)化設(shè)計時若不進行流場能量耗散的約束往往會產(chǎn)生大量懸浮固體，不具備制造條件，因此，本文對流場的能量耗散進行約束，流動能量耗散約束函數(shù)定義如下：

設(shè)計域內(nèi)流體材料的分布在提高結(jié)構(gòu)散熱性能的同時必定會影響結(jié)構(gòu)的剛度，結(jié)構(gòu)剛度不足會導(dǎo)致結(jié)構(gòu)承載性能的降低，因此，在優(yōu)化時引入結(jié)構(gòu)柔度約束函數(shù)以保證結(jié)構(gòu)剛度。柔度約束函數(shù)定義如下：

式中，為流體材料體積分數(shù)上限；為流動能量耗散約束上限；為結(jié)構(gòu)柔度約束上限。

2.2材料插值方案

在基于達西滲流理論的流體拓撲優(yōu)化框架下，BORRVALL等［9］首先采取基于有理數(shù)近似的插值策略（rationalapproximationofmaterialproperties，RAMP），在流體拓撲優(yōu)化問題中得到了很好的應(yīng)用。因此，本文采取RAMP策略對各物理場屬性進行插值描述。

流體和固體材料的流動屬性通過對式（2）中的材料流阻系數(shù)插值進行表征，插值格式為

其中，αmax為材料最大流阻系數(shù)，即固體材料的流阻系數(shù)；qα為控制懲罰強度的常數(shù)，本文取0.01；當(dāng)設(shè)計變量γ=1時，α（γ）=0，流動控制方程退化為N-S方程；當(dāng)γ=0時，α（γ）=αmax，此時流動阻力達到最大值，流速趨近于0，表示此處為固體材料。

流體和固體材料的傳熱屬性通過對材料熱導(dǎo)率插值進行表征，插值格式為

式中，kf、ks分別為流體區(qū)域和固體區(qū)域的熱導(dǎo)率；qk的作用與qα一致，本文取qk=0.01。

流體和固體材料的力學(xué)屬性通過對材料的彈性模量插值進行表征，插值格式為

其中，Emax表示固體材料的彈性模量，Emin表示流體材料即孔洞的彈性模量，通常取Emin=0.001來避免數(shù)值求解過程中出現(xiàn)的矩陣奇異問題。qE的作用與qα和qk一致，本文取qE=0.01。

3靈敏度計算方法

通過計算目標(biāo)函數(shù)對設(shè)計變量的靈敏度從而更新設(shè)計變量的分布是求解拓撲優(yōu)化問題的主要求解策略。由于本文設(shè)計變量γ的變化會引起狀態(tài)場量（速度、溫度、位移）的變化，而狀態(tài)場量對設(shè)計變量的導(dǎo)數(shù)隱含在控制方程中無法直接獲得，這給靈敏度分析帶來了困難。針對這一難點，本文采用基于拉格朗日乘子法和連續(xù)伴隨法的思想來求設(shè)計變量的靈敏度。與常見的離散伴隨法敏度分析相比，連續(xù)伴隨法需要手動推導(dǎo)連續(xù)形式的伴隨方程，進而求解伴隨方程獲得伴隨變量，最終靈敏度可由伴隨變量經(jīng)過簡單的計算獲得。

首先引入弱形式下的拉格朗日方程，即在整個求解域上將狀態(tài)方程與拉格朗日乘子（伴隨變量）相乘后增廣至目標(biāo)函數(shù)（或約束函數(shù)）上。拉格朗日函數(shù)L定義如下：

其中，Φ為待求靈敏度的目標(biāo)或約束函數(shù)；θ為拉格朗日乘子，即伴隨變量；u^、p^、T^、D^分別為伴隨速度、伴隨壓力、伴隨溫度和伴隨位移；R代表狀態(tài)控制方程，即式（1）～式（4）。狀態(tài)控制方程始終滿足R≡0，因此對于某個目標(biāo)或者約束函數(shù)來說始終有L≡Φ，求解目標(biāo)或者約束函數(shù)的靈敏度即可轉(zhuǎn)換為求解式（13）的靈敏度。令式（13）對狀態(tài)變量的變分為0可求得一組伴隨變量。通過引入伴隨變量，靈敏度分析時只需對式（13）中與設(shè)計變量相關(guān)的項直接求導(dǎo)，而不必考慮其余狀態(tài)量的中間效應(yīng)。

3.1能量耗散靈敏度

求能量耗散約束函數(shù)J的靈敏度時，式（13）中Φ=J。由于能量方程與流動方程是解耦求解的，且能量耗散約束與變量無關(guān)，故能量耗散約束的伴隨變量θ=（u^，p^，0，0）。為獲得能量耗散約束對設(shè)計變量的靈敏度，需引入拉格朗日函數(shù)的總變分：

為使任意擾動量δp和δu下，上式均成立，必須滿足其中的積分項為0。令上式中的體積分項為0，得到伴隨方程；令上式中的邊界積分項為0，得到伴隨邊界條件。最終得到的伴隨N-S方程及邊界條件：

需要指出的是，在求解伴隨方程的過程中認為原始方程的量是已知的，即在求解伴隨變量（u^，p^）時，原始場量（u，p）是已知量。至此，將式（21）的解u^i代入式（16）即可得到能量耗散約束的靈敏度。

3.2平均溫度靈敏度

目標(biāo)函數(shù)平均溫度Ψ的靈敏度推導(dǎo)過程與3.1節(jié)一致，需要注意的是該函數(shù)同時與u、p、T相關(guān)，因此需要求解一組伴隨N-S方程和伴隨能量方程以獲得伴隨變量（u^，p^，T^）。平均溫度的靈敏度表達式為

3.3柔度靈敏度

柔度約束函數(shù)C的靈敏度表達式為

4優(yōu)化求解策略

4.1優(yōu)化算法流程

本文使用基于有限體積法離散的OpenFOAM軟件作為多物理場問題和拓撲優(yōu)化問題的求解器，拓撲優(yōu)化求解算法采用移動漸近線法。優(yōu)化算法的流程圖見圖2，具體執(zhí)行步驟如下：

（1）對設(shè)計變量進行初始化，根據(jù)式（10）～式（12）更新材料的流阻系數(shù)、熱導(dǎo)率和彈性模量。

（2）分別求解式（1）、式（3）、式（4）得到當(dāng)前迭代步下的流速場u、溫度場T和位移場D，進而計算式（6）～式（8）得到當(dāng)前迭代步下的目標(biāo)函數(shù)和各約束函數(shù)的值，在求解溫度場T時使用了流速場u的數(shù)據(jù)。

（3）求解式（21）得到能量耗散約束函數(shù)的伴隨變量，求解式（23）、式（24）得到平均溫度目標(biāo)函數(shù)的伴隨變量。根據(jù)式（16）、式（22）和式（25）分別計算能量耗散約束、平均溫度目標(biāo)和結(jié)構(gòu)柔度約束的靈敏度。

（4）根據(jù)得到的靈敏度信息調(diào)用MMA算法更新設(shè)計變量。

（5）判斷當(dāng)前迭代步下設(shè)計變量是否完成收斂或迭代次數(shù)是否達到上限，若滿足條件則拓撲優(yōu)化流程結(jié)束，否則使用當(dāng)前迭代步更新的設(shè)計變量重新更新材料屬性，繼續(xù)迭代求解直至滿足程序結(jié)束條件。

4.2數(shù)值過濾方法

基于密度法的拓撲優(yōu)化流程一般存在棋盤格效應(yīng)、中間密度單元過多等數(shù)值不穩(wěn)定現(xiàn)象。為抑制上述數(shù)值不穩(wěn)定現(xiàn)象，本文使用基于Helmholtz型偏微分方程的PDE過濾［35］和基于雙曲正切函數(shù)的投影［36］對密度場進行過濾，經(jīng)兩層過濾后可得到邊界清晰的拓撲結(jié)果。

首先對設(shè)計變量場γ使用PDE過濾以消除棋盤格效應(yīng)，PDE過濾函數(shù)如下：

5數(shù)值算例

以燃氣輪機為代表的渦輪機械是典型的多物理場工況作用下的復(fù)雜結(jié)構(gòu)。本文以渦輪盤體和渦輪葉片作為研究對象，采用提出的流熱力拓撲優(yōu)化模型對其進行結(jié)構(gòu)設(shè)計。冷卻氣流由渦輪軸心位置處引入渦輪內(nèi)部，通過內(nèi)部冷卻通道截面的對流換熱帶走溫度，從而提高渦輪的承溫能力。渦輪材料及冷卻空氣的相關(guān)屬性見表1。

5.1渦輪盤體設(shè)計

給定渦輪盤體幾何模型如圖5所示，尺寸單

位為mm。為減少計算量，取渦輪盤體的一個旋轉(zhuǎn)周期作為計算域，渦輪盤體計算域及邊界情況如圖6a所示。為保證優(yōu)化設(shè)計時空氣流動能夠充分發(fā)展，入口和出口處都有一段延長導(dǎo)管。中間區(qū)域為拓撲優(yōu)化設(shè)計域，入口和出口處為純流動狀態(tài)的流體材料。計算域離散為6×105個六面體網(wǎng)格，計算時采用24核并行計算。

計算模型邊界條件設(shè)置如下：指定冷卻通道入口處的法向速度為0.05m/s，溫度為20℃，冷卻通道出口處為零壓力絕熱邊界，其余壁面為絕熱無滑移邊界。熱載荷由體熱源引入，給定熱源密度為10kW/m3；力載荷由葉片表面均布壓力引入，給定壓力大小為0.1MPa，同時固定轉(zhuǎn)軸處位移為0。

為了得到流熱力耦合拓撲優(yōu)化各項約束的參考值，首先對圖6b所示參考設(shè)計進行模擬。該結(jié)構(gòu)的流體材料體積分數(shù)為1.5%，經(jīng)仿真分析，圖6b所示結(jié)構(gòu)的能量耗散J0=1.0×10-5kg·m2/s3，結(jié)構(gòu)柔度C0=6×10-10J，平均溫度為130.7℃。為提高結(jié)構(gòu)的承溫能力，在參考設(shè)計的基礎(chǔ)上，設(shè)置許用流體材料分數(shù)為10%，能量耗散約束為=1.5J0，結(jié)構(gòu)柔度約束為=1.1C0，設(shè)計變量的初始值為0.5。拓撲優(yōu)化收斂準(zhǔn)則設(shè)置為最大迭代次數(shù)400或設(shè)計變量的增量‖γk+1－γk‖lt;10－3。

由于設(shè)計流場的多物理場問題存在極強的非線性，為保證數(shù)值計算穩(wěn)定，將初始最大流阻系數(shù)αmax設(shè)置為108，之后隨迭代次數(shù)的增加而逐漸增大至1011。同時，在優(yōu)化前期對能量耗散約束和結(jié)構(gòu)柔度約束也進行放松：在初始迭代階段，設(shè)置能量耗散約束上限為10J0，結(jié)構(gòu)柔度約束上限為2C0。隨著優(yōu)化過程的進行，控制這兩個約束上限分別線性減小至和。渦輪盤體的拓撲優(yōu)化迭代過程如圖7所示，目標(biāo)及約束迭代曲線如圖8所示。

在優(yōu)化前期，流阻系數(shù)較小，能量耗散值也較小；同時由于流體材料的大面積存在，結(jié)構(gòu)柔度較高；由于此時對固體材料的流速懲罰不足，固體材料內(nèi)部仍存在不合實際的對流換熱，故結(jié)構(gòu)的平均溫度較低。隨著優(yōu)化過程的進行，流阻系數(shù)逐漸增大，能量耗散約束上限和結(jié)構(gòu)柔度約束上限逐漸減小，流體材料由分散狀態(tài)逐漸收斂為完整的管道；同時固體材料內(nèi)部的流速逐漸懲罰為零，對流換熱只能發(fā)生在流體管道的表面，而且固體材料的質(zhì)量熱容遠小于流體材料的質(zhì)量熱容，因此結(jié)構(gòu)的平均溫度有所升高。優(yōu)化迭代至300步左右時，設(shè)計域域內(nèi)流體形成了貫通入口和出口的管道，此時管道內(nèi)流體形成連續(xù)的流動，流場能量耗散瞬時減小。最終結(jié)構(gòu)拓撲變化逐漸穩(wěn)定，平均溫度為120℃。圖9為盤體優(yōu)化前后的溫度場云圖，可以看出優(yōu)化后的溫度場整體上低于參考設(shè)計值，渦輪盤體最終設(shè)計結(jié)果如圖10所示。表2給出了渦輪盤體優(yōu)化前后的性能對比，與參考設(shè)計相比，拓撲優(yōu)化后結(jié)構(gòu)的平均溫度降低了8.2%。

5.2渦輪葉片設(shè)計

渦輪葉片是將燃氣熱能轉(zhuǎn)化為渦輪旋轉(zhuǎn)動能的核心部件，其工作溫度已接近甚至超過了材料的熔點，高效的渦輪葉片冷卻通道是保證渦輪穩(wěn)定工作的關(guān)鍵。本文對圖11a所示渦輪葉片模型進行設(shè)計，該模型含有一個入口和一個出口，計算域離散為3×105個六面體網(wǎng)格，計算時采用24核并行計算。

渦輪葉片計算模型邊界條件設(shè)置如下：指定冷卻通道入口處的法向速度為0.05m/s，溫度為20℃，冷卻通道出口處為零壓力絕熱邊界，其余壁面為絕熱無滑移邊界。熱載荷由體熱源引入，給定熱源密度為10kW/m3；力載荷由葉片表面均布壓力引入，給定壓力大小為0.1MPa，同時固定葉根處位移為0。

首先對圖11b所示參考設(shè)計進行模擬，以獲得流熱力耦合拓撲優(yōu)化各項約束的參考值。參考設(shè)計的平均溫度為69.7℃，流體材料體積分數(shù)為1%，流動能量耗散J0=5×10-7kg·m2/s3，結(jié)構(gòu)柔度C0=2×10-10J。由于渦輪葉片處受熱較多且承受的力學(xué)載荷相對較少，故可以適當(dāng)增加流體材料占比。在參考設(shè)計的基礎(chǔ)上設(shè)置許用流體材料分數(shù)為40%，能量耗散約束為=12J0，結(jié)構(gòu)柔度約束為=1.7C0，設(shè)計變量的初始值為0.5。拓撲優(yōu)化收斂準(zhǔn)則設(shè)置為最大迭代次數(shù)400或設(shè)計變量的增量‖γk+1－γk‖lt;10－3。本算例采取與5.1節(jié)一致的放松策略，將初始最大流阻系數(shù)αmax設(shè)置為108，之后隨迭代次數(shù)的增加而逐漸增大至1011。同時，在初始迭代階段，

設(shè)置能量耗散約束上限為25J0，結(jié)構(gòu)柔度約束上限為3C0。隨著優(yōu)化過程的進行，控制這兩個約束上限分別線性減小至和。

圖12和圖13給出了渦輪葉片優(yōu)化設(shè)計的迭代過程和目標(biāo)及約束函數(shù)的迭代曲線。從優(yōu)化結(jié)果可以看出，當(dāng)流體材料許用體積和流動能量耗散許用值增加后，生成的最優(yōu)構(gòu)型更加復(fù)雜。由于大量增加了流體材料用量，即使在優(yōu)化初期對能量耗散和結(jié)構(gòu)柔度約束進行了大量放松，仍難以找到穩(wěn)定的優(yōu)化方向，優(yōu)化前期出現(xiàn)了少許振蕩。隨著流阻系數(shù)的增加，能量耗散逐步增大；能量耗散約束上限及結(jié)構(gòu)柔度約束上限隨迭代步數(shù)的增加逐漸縮緊到許用值，最終在許用能量耗散和結(jié)構(gòu)柔度約束下，平均溫度目標(biāo)收斂至54.5℃。圖14給出了葉片優(yōu)化前后的溫度場分布情況，可以看出優(yōu)化后葉片的溫度明顯低于參考設(shè)計，葉片流道截面形狀如圖15所示。表3給出了葉片優(yōu)化前后的性能對比，與參考設(shè)計相比，拓撲優(yōu)化后葉片結(jié)構(gòu)的平均溫度降低了21.5%。

6結(jié)論

（1）本文針對流熱力多物理場工況下的高性能復(fù)雜結(jié)構(gòu)設(shè)計這一挑戰(zhàn)，提出了一種復(fù)雜結(jié)構(gòu)流熱力耦合拓撲優(yōu)化方法，為高性能結(jié)構(gòu)的設(shè)計和制造提供了新的方案。

（2）基于變分法和拉格朗日乘子法建立了多物理場拓撲優(yōu)化問題的連續(xù)伴隨敏度分析框架，推導(dǎo)了伴隨狀態(tài)方程及其邊界條件。使用移動漸近線法對拓撲優(yōu)化問題進行求解，結(jié)合Helmholtz型偏微分方程過濾和雙曲正切方程投影，改善了拓撲優(yōu)化過程的不穩(wěn)定現(xiàn)象，得到邊界清晰的可增材制造拓撲結(jié)構(gòu)。

（3）將所提的流熱力拓撲優(yōu)化模型應(yīng)用于實際工程裝備的設(shè)計，分別對渦輪盤體和渦輪葉片進行結(jié)構(gòu)設(shè)計，得到了散熱性能良好、流道分布合理的高性能結(jié)構(gòu)。通過比較發(fā)現(xiàn)，本文方法可以獲得比參考設(shè)計承溫能力更好的結(jié)構(gòu)。因此，考慮流熱力耦合的拓撲優(yōu)化設(shè)計有其實際必要性。

參考文獻：

［1］LIUJikai，GAYNORAT，CHENShikui，etal.CurrentandFutureTrendsinTopologyOptimizationforAdditiveManufacturing［J］.StructuralandMultidisciplinaryOptimization，2018，57（6）：2457-2483.

［2］劉成穎，譚鋒，王立平，等.面向機床整機動態(tài)性能的立柱結(jié)構(gòu)優(yōu)化設(shè)計研究［J］.機械工程學(xué)報，2016，52（3）：161-168.

LIUChengying，TANFeng，WANGLiping，etal.ResearchonOptimizationofColumnStructureDesignforDynamicPerformanceofMachineTool［J］.JournalofMechanicalEngineering，2016，52（3）：161-168.

［3］DUPLESSISA，BROECKHOVENC，YADROITSAVAI，etal.BeautifulandFunctional：aReviewofBiomimeticDesigninAdditiveManufacturing［J］.AdditiveManufacturing，2019，27：408-427.

［4］劉偉，李能，周標(biāo)，等.復(fù)雜結(jié)構(gòu)與高性能材料增材制造技術(shù)進展［J］.機械工程學(xué)報，2019，55（20）：128-151.

LIUWei，LINeng，ZHOUBiao，etal.ProgressinAdditiveManufacturingonComplexStructuresandHigh-performanceMaterials［J］.JournalofMechanicalEngineering，2019，55（20）：128-151.

［5］朱繼宏，周涵，王創(chuàng)，等.面向增材制造的拓撲優(yōu)化技術(shù)發(fā)展現(xiàn)狀與未來［J］.航空制造技術(shù)，2020，63（10）：24-38.

ZHUJihong，ZHOUHan，WANGChuang，etal.StatusandFutureofTopologyOptimizationforAdditiveManufacturing［J］.AeronauticalManufacturingTechnology，2020，63（10）：24-38.

［6］SIGMUNDO，MAUTEK.TopologyOptimizationApproaches［J］.StructuralandMultidisciplinaryOptimization，2013，48（6）：1031-1055.

［7］IGAA，NISHIWAKIS，IZUIK，etal.TopologyOptimizationforThermalConductorsConsideringDesign-dependentEffects，IncludingHeatConductionandConvection［J］.InternationalJournalofHeatandMassTransfer，2009，52（11/12）：2721-2732.

［8］SUBRAMANIAMV，DBOUKT，HARIONJL.TopologyOptimizationofConductiveHeatTransferDevices：anExperimentalInvestigation［J］.AppliedThermalEngineering，2018，131：390-411.

［9］BORRVALLT，PETERSSONJ.TopologyOptimizationofFluidsinStokesFlow［J］.InternationalJournalforNumericalMethodsinFluids，2003，41（1）：77-107.

［10］DHRINGMB，JENSENJS，SIGMUNDO.AcousticDesignbyTopologyOptimization［J］.JournalofSoundandVibration，2008，317（3/4/5）：557-575.

［11］RUPPCJ，EVGRAFOVA，MAUTEK，etal.DesignofPhononicMaterials/StructuresforSurfaceWaveDevicesUsingTopologyOptimization［J］.StructuralandMultidisciplinaryOptimization，2007，34（2）：111-121.

［12］SUTRADHARA，PAULINOGH，MILLERMJ，etal.TopologicalOptimizationforDesigningPatient-specificLargeCraniofacialSegmentalBoneReplacements［J］.ProceedingsoftheNationalAcademyofSciencesoftheUnitedStatesofAmerica，2010，107（30）：13222-13227.

［13］PARKJ，LEEDamian，SUTRADHARA.TopologyOptimizationofFixedCompleteDentureFramework［J］.InternationalJournalforNumericalMethodsinBiomedicalEngineering，2019，35（6）：e3193.

［14］FERRARIF，SIGMUNDO.ANewGeneration99LineMatlabCodeforComplianceTopologyOptimizationandItsExtensionto3D［J］.StructuralandMultidisciplinaryOptimization，2020，62（4）：2211-2228.

［15］WUJun，CLAUSENA，SIGMUNDO.MinimumComplianceTopologyOptimizationofShell–infillCompositesforAdditiveManufacturing［J］.ComputerMethodsinAppliedMechanicsandEngineering，2017，326：358-375.

［16］WANGHui，CHENGWenming，DURun，etal.ImprovedProportionalTopologyOptimizationAlgorithmforSolvingMinimumComplianceProblem［J］.StructuralandMultidisciplinaryOptimization，2020，62（2）：475-493.

［17］KRIEGESMANNB，LDEKERJK.RobustComplianceTopologyOptimizationUsingtheFirst-orderSecond-momentMethod［J］.StructuralandMultidisciplinaryOptimization，2019，60（1）：269-286.

［18］RODRIGUESH，F(xiàn)ERNANDESP.AMaterialBasedModelforTopologyOptimizationofThermoelasticStructures［J］.InternationalJournalforNumericalMethodsinEngineering，1995，38（12）：1951-1965.

［19］YANGXiongwei，LIYueming.TopologyOptimizationtoMinimizetheDynamicComplianceofaBi-materialPlateinaThermalEnvironment［J］.StructuralandMultidisciplinaryOptimization，2013，47（3）：399-408.

［20］PARKJ，NGUYENTH，SHAHJJ，etal.ConceptualDesignofEfficientHeatConductorsUsingMulti-materialTopologyOptimization［J］.EngineeringOptimization，2019，51（5）：796-814.

［21］DEATONJD，GRANDHIRV.StiffeningofRestrainedThermalStructuresviaTopologyOptimization［J］.StructuralandMultidisciplinaryOptimization，2013，48（4）：731-745.

［22］ZHUXuefeng，ZHAOChao，WANGXuan，etal.Temperature-constrainedTopologyOptimizationofThermo-mechanicalCoupledProblems［J］.EngineeringOptimization，2019，51（10）：1687-1709.

［23］TAKEZAWAA，YOONGH，JEONGSH，etal.StructuralTopologyOptimizationwithStrengthandHeatConductionConstraints［J］.ComputerMethodsinAppliedMechanicsandEngineering，2014，276：341-361.

［24］DEATONJD，GRANDHIRV.Stress-basedDesignofThermalStructuresviaTopologyOptimization［J］.StructuralandMultidisciplinaryOptimization，2016，53（2）：253-270.

［25］ZHENGJing，CHENHong，JIANGChao.RobustTopologyOptimizationforStructuresunderThermo-mechanicalLoadingsConsideringHybridUncertainties［J］.StructuralandMultidisciplinaryOptimization，2022，65（1）：29.

［26］YOONGH.TopologyOptimizationforStationaryFluid-structureInteractionProblemsUsingaNewMonolithicFormulation［J］.InternationalJournalforNumericalMethodsinEngineering，2010，82（5）：591-616.

［27］PICELLIR，VICENTEWM，PAVANELLOR，etal.EvolutionaryTopologyOptimizationforNaturalFrequencyMaximizationProblemsConsideringAcoustic-structureInteraction［J］.FiniteElementsinAnalysisandDesign，2015，106（C）：56-64.

［28］LIHao，KONDOHT，JOLIVETP，etal.Three-dimensionalTopologyOptimizationofaFluid–StructureSystemUsingBody-fittedMeshAdaptionBasedontheLevel-setMethod［J］.AppliedMathematicalModelling，2022，101：276-308.

［29］YOONGH.TopologicalDesignofHeatDissipatingStructurewithForcedConvectiveHeatTransfer［J］.JournalofMechanicalScienceandTechnology，2010，24（6）：1225-1233.

［30］YUMinghao，RUANShilun，WANGXinyu，etal.TopologyOptimizationofThermal–fluidProblemUsingtheMMC-basedApproach［J］.StructuralandMultidisciplinaryOptimization，2019，60（1）：151-165.

［31］YAJIK，YAMADAT，KUBOS，etal.ATopologyOptimizationMethodforaCoupledThermal–fluidProblemUsingLevelSetBoundaryExpressions［J］.InternationalJournalofHeatandMassTransfer，2015，81：878-888.

［32］ZHAOXi，ZHOUMingdong，LIUYichang，etal.TopologyOptimizationofChannelCoolingStructuresConsideringThermomechanicalBehavior［J］.StructuralandMultidisciplinaryOptimization，2019，59（2）：613-632.

［33］PANTONRL.IncompressibleFlow［M］.Hoboken：Wiley，2013.

［34］于明豪.液冷散熱器流道與熱源分布拓撲優(yōu)化研究［D］.大連：大連理工大學(xué)，2021.

YUMinghao.ResearchonTopologyOptimizationofCoolingChannelsandHeatSourcesDistributionforLiquid-cooledHeatSink［D］.Dalian：DalianUniversityofTechnology，2021.

［35］LAZAROVBS，SIGMUNDO.FiltersinTopologyOptimizationBasedonHelmholtz-typeDifferentialEquations［J］.InternationalJournalforNumericalMethodsinEngineering，2011，86（6）：765-781.

［36］XUShengli，CAIYuanwu，CHENGGengdong.VolumePreservingNonlinearDensityFilterBasedonHeavisideFunctions［J］.StructuralandMultidisciplinaryOptimization，2010，41（4）：495-505.