針對注塑成型隨形冷卻的流熱力耦合拓撲優化方法研究

(1.大連理工大學 工業裝備結構分析國家重點實驗室，大連 116024；2.鄭州大學 橡塑模具國家工程研究中心，鄭州 450002)

1 引 言

塑料制品廣泛應用于汽車、建筑、機械和農業等領域，是我國經濟發展支柱性產業之一。隨著材料技術的發展，研發出各種工程塑料和特種塑料，為信息電子和航空航天等高科技產業的發展提供制造可能。與此同時，對塑料制品的性能、制造精度和外觀都有了更高要求。注塑成型作為塑料加工的主要方法之一，因其生產精度高、便于制造形狀復雜的產品以及可批量生產等特點，占據了當前塑料工業生產的重要位置。注塑過程中，高溫熔體傳入模具的熱量95%以上經由散熱管道內的冷卻劑帶出模具，其余部分以對流和輻射方式從模具外表面傳出[1]。模具冷卻作為注塑成型的關鍵環節[2]，其型腔的冷卻速度和均勻性影響著塑料制品的生產效率以及產品質量，因此如何使模具更快和更均勻地散熱成為當前注塑成型工藝的一個重點問題[3]。

傳統的工藝設計方法主要依賴于試模，但開模成本較高，為滿足高質量塑料制品的生產需求，許多學者基于數值模擬技術對注塑成型模具冷卻系統進行參數優化設計。Tang等[4]采用伽遼金有限元法的無矩陣雅可比共軛梯度形式模擬瞬態熱傳導，對冷卻管道的尺寸、位置和冷卻劑流量進行設計。Dang等[5]以汽車塑料擋泥板為例，將銑削加工中U形銑槽冷卻管道的直徑、間距和深度作為設計變量，借助三維CAE仿真工具，結合DOE試驗設計方法建立響應面模型，對冷卻管道進行優化。Li等[6]在注塑模具冷卻系統布局階段，考慮不干擾其他模具部件的情況下，在模具鑲件處建立冷卻系統，提出配置空間(C-space)的概念，結合遺傳算法進行冷卻系統布局可行設計搜索。

還有一些學者基于直觀感受或工程經驗，引入啟發式方法進行散熱管道設計。Peng等[7]以植物葉片的脈絡結構為基礎，設計了帶有滲透壁和多孔介質的分支網絡散熱裝置。Barber等[8]受到動物血管的仿生學原理啟發，對矩形、梯形截面通道微流控系統進行研究。對于三維的復雜結構散熱問題，Au等[9]用多面體形狀近似模具表面，并根據照明原理在設計空間內找出能夠照亮多面體形狀的若干光源點，再通過連接光源點進行冷卻通道設計，使模具表面達到均勻冷卻目的。上述參數優化工作都是針對傳統的平行管道進行設計，雖然較于原始設計有明顯的改進效果，但其優化結構依賴于經驗給出的初始管道結構，具有較強的局限性。與參數優化相比，啟發式方法可獲得優于傳統管道設計方法的結果，且計算效率高，但其并未對模具冷卻過程進行數值模擬分析，無法確定是否獲得最優解。

相比參數優化的局限性以及啟發式方法的精確性不足，拓撲優化結合數值模擬技術進行材料分布優化具有較高的計算可靠性和設計自由性。近年來，許多學者將結構拓撲優化方法引入流熱耦合領域中解決散熱器優化設計問題。Dede等[10]對微流控分析芯片進行多通道拓撲的單相流微通道冷卻板研究。Ya等[11]基于二維和三維穩態N-S(Navier-Stokes)方程及能量方程，將散熱結構設計表述為一個熱交換最大化問題，采用水平集法對管道邊界進行描述。Tawk等[12]建立bi-eta插值函數以及流動連續性函數，將密度法用于包含兩種流體和一種固體的多相熱交換器的結構設計中，實現了二維的雙流體傳熱傳質結構拓撲優化。Zheng等[13]基于周期平均法和邊界元研究了注塑成型隨形冷卻管道拓撲優化設計。該方法首先基于隨形面給定初始管網分布，并將管道節點坐標和半徑作為優化變量，最終得到復雜的三維管道。該工作采用了工程計算模型來分析管道流速，雖然可以大幅度提高計算效率，但是對于曲率很大的管道網絡，計算精度無法保證。

本文以偽密度作為設計變量，對流固材料進行插值描述，采用移動漸近線法求解拓撲優化問題，結合連續伴隨法進行靈敏度分析，對注塑成型模具冷卻管道的拓撲優化設計展開研究。本文還對比了是否考慮模具結構承載力的拓撲優化結果，發現不考慮結構承載力的冷卻管道優化設計不具備工程應用意義，而將結構柔度作為約束條件之一，可以較好地解決實際加工中考慮模具承載能力影響的冷卻管道設計問題。

2 控制方程

2.1 流動控制方程

在注塑成型模具的冷卻管道設計中，冷卻劑常視為粘性不可壓縮定常流。一般情況下湍流狀態冷卻介質的散熱能力優于層流，但基于密度法的湍流拓撲優化由于無法識別流固邊界，目前仍存在諸多問題，如無法設置壁面函數、無法對邊界層網格加密和難以準確計算近壁距離等。故本文選擇定常層流N-S方程作為流動控制方程。

連續性方程

·u=0

(1)

動量守恒方程

(2)

在大多數基于變密度法的流體拓撲優化工作中，通常以達西滲流理論作為控制方程出發點，即增加多孔介質流阻力項f=-αfu來表征流固材料分布情況。因此動量守恒方程可表示為

(3)

式中αf為材料流阻系數，其通過式(3)進行材料密度插值[14]。

(3)

2.2 傳熱控制方程

注塑成型是一個周期性循環過程，這決定了型腔面溫度也隨時間進行變化。注塑過程進入周期穩態后，模具表面瞬時溫度雖然會變化，但其平均溫度不變。在實際的模具冷卻分析應用中，一般利用周期平均溫度代替瞬時溫度，將原瞬態問題轉化為穩態問題[1]。故注塑成型問題的傳熱控制方程可表示為

(5)

(6)

(7)

熱傳導系數插值

(8)

式中qk為正常數，本文取0.01；kf為流體熱傳導系數；ks為固體熱傳導系數。

2.3 結構平衡方程

當塑料熔體以高壓注入模具時，注塑機的鎖模元件要對模具施加足夠的鎖模力克服型腔內熔體的脹模力。由于冷卻管道會影響模具的剛度，為保證模具在巨大的鎖模力下不受損壞，則應在散熱系統優化設計時考慮結構承載力。本文基于模具結構承載力因素進行算例對比研究，線彈性力學平衡方程可表示為

(9)

(10)

(11)

(12,13)

(14)

式中qE為正常數，本文取0.01；Emin為孔洞的楊氏模量，Emax為固體材料的楊氏模量，為避免計算過程中出現矩陣奇異，取Emin=10-9Emax。

2.4 拓撲優化列式

散熱性能是冷卻系統的主要關心指標，本文以平均溫度為目標函數衡量其散熱效果,

(15)

式中Ψ表示目標函數，|Ω|為計算域總體積。

冷卻系統的散熱能力與流場的能量耗散有關，一般情況下，流場壓力降越大則散熱效果越好。但在工程應用中其受到冷水機功率限制，故本文對流場的能量耗散進行約束。

(16)

式中Γin,out為冷卻管道的出入口邊界。

根據2.3節對結構平衡方程的介紹可知，注塑過程中要保證模具能承受住鎖模力，因此本文以結構柔度作為約束函數之一,

(17)

綜上所述，基于變密度法的流-熱-力耦合拓撲優化問題的數學列式可表示為

MinimizeΨ

(18)

2.5 過濾方法

基于變密度法的拓撲優化會出現棋盤格效應以及產生灰度單元等問題。為了解決上述問題，本文聯合使用Helmholtz偏微分方程[15]和雙曲正切投影[16]進行密度過濾。

Helmholtz偏微分方程表達式為

(19)

雙曲正切投影表達式為

(20)

3 拓撲優化數值實現

3.1 優化方法

拓撲優化問題的目標函數及約束函數一般為設計變量的隱式表達，且拓撲優化問題通常具有非凸性和高非線性，計算量大。當前求解拓撲優化問題的主要方法分為優化準則法和序列規劃法兩類。在序列規劃法中，移動漸近線法MMA(method of movingasymptotes)[17]的變量獨立且只需求解一階導數就能達到很好的收斂效果，故本文采用此方法進行求解。

3.2 靈敏度分析

由于拓撲優化中設計變量個數遠大于目標和約束函數個數，因此一般采用伴隨法敏度分析一次計算出目標函數(或約束)對所有設計變量的導數。伴隨法又可分為連續伴隨法和離散伴隨法兩種，其區別在于對控制方程進行離散和求導的先后次序。離散伴隨法需要用到雅克比矩陣以及雅克比矩陣對設計變量的導數。如果控制方程采用分離式求解算法(如求解N-S方程的SIMPLE[18]算法)，則雅克比矩陣在計算過程中并未組集，這給伴隨敏度分析帶來困難。本文采用連續伴隨法計算靈敏度，其出發點為將連續的狀態方程與伴隨變量相乘并增廣到目標函數上，通過對拉格朗日函數求加托導數以及分部積分，最終得到連續形式的伴隨方程與伴隨邊界條件。該方法優點在于伴隨方程和狀態方程相互獨立，可以選擇對伴隨方程使用耦合或分離式求解器。該優點和支持自定義PDE求解器的CAE軟件結合時，可以極大地減少用戶在敏度分析方面的編程操作。

拉格朗日函數

(21)

(22)

式中Ψ為目標函數，此處引入如式(22)所示一般化積分形式進行表述。θ為拉格朗日乘子，R為控制方程。對于本文來說，R即為N-S方程、能量守恒方程及結構平衡方程。θ=(ua,pa,Ta)為伴隨變量。

(23)

根據偏微分方程約束最優化問題的KKT條件有

(24)

則設計變量靈敏度可求解得

(25)

從式(22,23)可推得伴隨Navier-Stokes方程為

?A/?p=·uainΩ

(26)

inΩ

(27)

(28)

-?A/(?uT)·n-?B/?u(u·n)ua+

(29)

伴隨能量守恒方程

inΩ(30)

(31)

-?A/(?T)·n-?B/?T=cTau+kTa]·n

式中ua和Ta分別為伴隨速度場和伴隨溫度場，式(24～30)的左端項列入表1。

表1 伴隨方程左端項

由于最小化結構柔度的線彈性力學問題是自伴隨的，即伴隨變量等于結構位移，此處省略該問題的伴隨敏度分析。

3.3 基于OpenFOAM進行優化問題求解

本文基于OpenFOAM[19]搭建拓撲優化求解器。流-熱-力耦合拓撲優化流程如圖1所示，首先讀取材料參數，再對Navier-Stokes方程以及結構平衡方程進行計算，得到速度場和位移場之后，將速度場帶入能量守恒方程求得溫度場。完成以上初始場計算后，對伴隨控制方程進行求解，得到各伴隨變量，然后計算目標函數以及約束函數對設計變量的靈敏度。最后調用MMA更新設計變量。

圖1 基于OpenFOAM拓撲優化求解流程

4 算 例

本文采用的注塑成型初始模具的幾何形狀與尺寸如圖2所示，尺寸單位為毫米(mm)。

圖2 對稱結構尺寸

4.1 不考慮結構承載的流-熱耦合拓撲優化

4.1.1 算例參數

本文算例的模具冷卻管道不可設計域如圖3所示，分別為模具冷卻管道的出入口部分以及模具型腔面。為減少拓撲優化過程中的計算量，對模型施加對稱邊界條件。冷卻管道入口處法向速度為 0.01 m/s，溫度為0，冷卻管道出口處為絕熱邊界。熱源面施加均布壓力10 MPa和熱流量40000 W/m2，其他邊界面為絕熱無滑移面。本算例選用的材料參數列入表2。

圖3 模具不可設計域

圖4 模具工況

表2 材料參數

4.1.2 傳統冷卻管道設計

在注塑成型模具中，傳統的冷卻管道通常為直管設計，本文將傳統直管設計作為拓撲優化初始構型進行模具冷卻仿真模擬。初始直管設計的流體材料用量百分比為1.4%，模具平均溫度為83.5 ℃，能量耗散為6.861×10-10kg·m2/s3，結構柔度為 0.086 J。因此，在本文優化算例中，給定流體材料百分比上限為3%，能量耗散為直管模型的2倍，結構柔度為直管模型的1.5倍。

圖5 傳統直管冷卻溫度場

4.1.3 優化結果

本算例不考慮模具結構承載能力進行拓撲優化，其優化列式為

MinimizeΨ

(33)

優化后管道設計如圖6所示。優化后的管道構型與模具型腔緊密貼合，結構整體剛度不足，在注塑過程中模具易損壞，且不具備模具實際制造條件，故下文將考慮模具結構承載力重新進行冷卻管道優化設計。

圖6 不考慮模具承載能力優化結果

4.2 考慮結構承載的流-熱-力耦合拓撲優化

在4.1節的優化模型基礎上增加考慮結構承載力的結構柔度約束進行優化設計，優化列式如式(18)所示。

本算例冷卻管道拓撲優化構型如圖7所示，冷卻管道沿著模具型腔面均勻覆蓋。結構溫度場如圖8所示，結構最高溫度出現在模具型腔面臨近出口邊緣處，優化結果顯示此處形成了流道以帶走模具熱量。拓撲優化后的模具平均溫度為76.0 ℃，與傳統直管設計相比，拓撲優化后的冷卻管道設計使模具的平均溫度降低了9%。

圖7 考慮結構承載能力的拓撲優化管道構型

圖8 考慮結構承載能力優化后溫度場

5 結 論

本文基于變密度法進行注塑模具冷卻管道拓撲優化設計，采用插值函數描述流固材料分布，以最小結構平均溫度為設計目標。聯合使用Helmholtz偏微分方程和雙曲正切投影方程進行密度過濾，得到了清晰的冷卻管道拓撲形態。通過對比是否考慮模具承載能力的兩個算例，得出注塑成型中應該將模具承載性能作為評價指標的結論。最后基于OpenFOAM平臺對注塑成型冷卻系統設計問題進行三維流-熱-力耦合拓撲優化，得到了冷卻效果良好且具有實際加工意義的流體通道設計。與傳統的直管冷卻系統相比，拓撲優化后的冷卻系統能夠使模具平均溫度降低9%。

本文僅是對于注塑成型冷卻系統設計的初步研究，目前還存在諸多不足,如應采用周期平均法計算模具型腔面處的平均壓力和熱源值當做優化的邊界條件；應該采用湍流模型模擬冷卻介質的流動傳熱性能。未來將從以上角度對優化進行改進。

參考文獻(References):

[1] 石憲章.注塑冷卻數值分析方法的研究與應用[D].鄭州大學,2005.(SHI Xian-zhang.Study and Application of Numerical Method for Cooling Analysis in Injection Molding[D].Zhengzhou University,2005.(in Chinese))

[2] 申長雨,李海梅,高 峰.注射成型技術發展概況[J].工程塑料應用,2003，31(3):53-57.(SHEN Chang-yu,LI Hai-mei,GAO Feng.Development survey of injection molding[J].EngineeringPlasticsApplication,2003，31(3):53-57.(in Chinese))

[3] 王 勇.基于快速隨形熱冷的高光無痕注塑成型[D].湖北工業大學,2018.(WANG Yong.High Gloss Seamless Injection Molding Based on Rapid Shape Following Thermalcool[D].Hubei University of Technology,2018.(in Chinese))

[4] Tang L Q,Chassapis C,Manoochehri S.Optimal cooling system design for multi-cavity injection molding [J].FiniteElementsinAnalysisandDe-sign,1997,26(3):229-251.

[5] Dang X P,Park H S.Design of U-shape milled groove conformal cooling channels for plastic injection mold[J].InternationalJournalofPrecisionEngineeringandManufacturing,2011,12(1):73-84.

[6] Li C G,Li C L.Plastic injection mould cooling system design by the configuration space method[J].Computer-AidedDesign,2008,40(3):334-349.

[7] Peng Y,Liu W Y,Chen W,et al.A conceptual structure for heat transfer imitating the transporting pri-nciple of plant leaf[J].InternationalJournalofHeatandMassTransfer,2014,71:79-90.

[8] Barber R W,Emerson D R.Optimal design of microfluidic networks using biologically inspired principles[J].MicrofluidicsandNanofluidics,2008,4(3): 179-191.

[9] Au K M,Yu K M,Chiu W K.Visibility-based conformal cooling channel generation for rapid tooling[J].Computer-AidedDesign,2011,43(4):356-373.

[10] Dede E M.Optimization and design of a multipass branching microchannel heat sink for electronics cooling[J].JournalofElectronicPackaging,2012,134(4):041001.

[11] Yaji K,Yamada T,Kubo S,et al.A topology optimization method for a coupled thermal-fluid problem using level set boundary expressions[J].Interna-tionalJournalofHeatandMassTransfer,2015,81:878-888.

[12] Tawk R,Ghannam B,Nemer M.Topology optimization of heat and mass transfer problems in two fluids — one solid domains[J].NumericalHeatTransfer,PartB:Fundamentals,2019,76(3):130-151.

[13] Li Z,Wang X Y,Gu J F,et al.Topology optimization for the design of conformal cooling system in thin-wall injection molding based on BEM[J].TheInternationalJournalofAdvancedManufacturingTechnology,2018,94(1-4):1041-1059.

[14] Borrvall T,Petersson J.Topology optimization of fluids in Stokes flow[J].InternationalJournalforNumericalMethodsinFluids,2003,41(1):77-107.

[15] Lazarov B S,Sigmund O.Filters in topology optimization based on Helmholtz-type differential equations[J].InternationalJournalforNumericalMethodsinEngineering,2011,86(6):765-781.

[16] Wang F W,Lazarov B S,Sigmund O.On projection methods,convergence and robust formulations in topology optimization[J].StructuralandMultidisciplinaryOptimization,2011,43(6):767-784.

[17] Svanberg K.The method of moving asymptotes—a new method for structural optimization[J].InternationalJournalforNumericalMethodsinEnginee-ring,1987,24(2):359-373.

[18] Patankar S V,Spalding D B.A calculation procedure for heat,mass and momentum transfer in three -dimensional parabolic flows[J].InternationalJournalofHeatandMassTransfer,1972,15(10):1787-1806.

[19] Open FOAM.AC++LibraryforComplexPhysicsSimulations[M].2007.