(1.大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室，大連 116024；2.鄭州大學(xué) 橡塑模具國(guó)家工程研究中心，鄭州 450002)

1 引 言

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

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

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

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

本文以偽密度作為設(shè)計(jì)變量，對(duì)流固材料進(jìn)行插值描述，采用移動(dòng)漸近線法求解拓?fù)鋬?yōu)化問(wèn)題，結(jié)合連續(xù)伴隨法進(jìn)行靈敏度分析，對(duì)注塑成型模具冷卻管道的拓?fù)鋬?yōu)化設(shè)計(jì)展開研究。本文還對(duì)比了是否考慮模具結(jié)構(gòu)承載力的拓?fù)鋬?yōu)化結(jié)果，發(fā)現(xiàn)不考慮結(jié)構(gòu)承載力的冷卻管道優(yōu)化設(shè)計(jì)不具備工程應(yīng)用意義，而將結(jié)構(gòu)柔度作為約束條件之一，可以較好地解決實(shí)際加工中考慮模具承載能力影響的冷卻管道設(shè)計(jì)問(wèn)題。

2 控制方程

2.1 流動(dòng)控制方程

在注塑成型模具的冷卻管道設(shè)計(jì)中，冷卻劑常視為粘性不可壓縮定常流。一般情況下湍流狀態(tài)冷卻介質(zhì)的散熱能力優(yōu)于層流，但基于密度法的湍流拓?fù)鋬?yōu)化由于無(wú)法識(shí)別流固邊界，目前仍存在諸多問(wèn)題，如無(wú)法設(shè)置壁面函數(shù)、無(wú)法對(duì)邊界層網(wǎng)格加密和難以準(zhǔn)確計(jì)算近壁距離等。故本文選擇定常層流N-S方程作為流動(dòng)控制方程。

連續(xù)性方程

·u=0

(1)

動(dòng)量守恒方程

(2)

在大多數(shù)基于變密度法的流體拓?fù)鋬?yōu)化工作中，通常以達(dá)西滲流理論作為控制方程出發(fā)點(diǎn)，即增加多孔介質(zhì)流阻力項(xiàng)f=-αfu來(lái)表征流固材料分布情況。因此動(dòng)量守恒方程可表示為

(3)

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

(3)

2.2 傳熱控制方程

注塑成型是一個(gè)周期性循環(huán)過(guò)程，這決定了型腔面溫度也隨時(shí)間進(jìn)行變化。注塑過(guò)程進(jìn)入周期穩(wěn)態(tài)后，模具表面瞬時(shí)溫度雖然會(huì)變化，但其平均溫度不變。在實(shí)際的模具冷卻分析應(yīng)用中，一般利用周期平均溫度代替瞬時(shí)溫度，將原瞬態(tài)問(wèn)題轉(zhuǎn)化為穩(wěn)態(tài)問(wèn)題[1]。故注塑成型問(wèn)題的傳熱控制方程可表示為

(5)

(6)

(7)

熱傳導(dǎo)系數(shù)插值

(8)

式中qk為正常數(shù)，本文取0.01；kf為流體熱傳導(dǎo)系數(shù)；ks為固體熱傳導(dǎo)系數(shù)。

2.3 結(jié)構(gòu)平衡方程

當(dāng)塑料熔體以高壓注入模具時(shí)，注塑機(jī)的鎖模元件要對(duì)模具施加足夠的鎖模力克服型腔內(nèi)熔體的脹模力。由于冷卻管道會(huì)影響模具的剛度，為保證模具在巨大的鎖模力下不受損壞，則應(yīng)在散熱系統(tǒng)優(yōu)化設(shè)計(jì)時(shí)考慮結(jié)構(gòu)承載力。本文基于模具結(jié)構(gòu)承載力因素進(jìn)行算例對(duì)比研究，線彈性力學(xué)平衡方程可表示為

(9)

(10)

(11)

(12,13)

(14)

式中qE為正常數(shù)，本文取0.01；Emin為孔洞的楊氏模量，Emax為固體材料的楊氏模量，為避免計(jì)算過(guò)程中出現(xiàn)矩陣奇異，取Emin=10-9Emax。

2.4 拓?fù)鋬?yōu)化列式

散熱性能是冷卻系統(tǒng)的主要關(guān)心指標(biāo)，本文以平均溫度為目標(biāo)函數(shù)衡量其散熱效果,

(15)

式中Ψ表示目標(biāo)函數(shù)，|Ω|為計(jì)算域總體積。

冷卻系統(tǒng)的散熱能力與流場(chǎng)的能量耗散有關(guān)，一般情況下，流場(chǎng)壓力降越大則散熱效果越好。但在工程應(yīng)用中其受到冷水機(jī)功率限制，故本文對(duì)流場(chǎng)的能量耗散進(jìn)行約束。

(16)

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

根據(jù)2.3節(jié)對(duì)結(jié)構(gòu)平衡方程的介紹可知，注塑過(guò)程中要保證模具能承受住鎖模力，因此本文以結(jié)構(gòu)柔度作為約束函數(shù)之一,

(17)

綜上所述，基于變密度法的流-熱-力耦合拓?fù)鋬?yōu)化問(wèn)題的數(shù)學(xué)列式可表示為

MinimizeΨ

(18)

2.5 過(guò)濾方法

基于變密度法的拓?fù)鋬?yōu)化會(huì)出現(xiàn)棋盤格效應(yīng)以及產(chǎn)生灰度單元等問(wèn)題。為了解決上述問(wèn)題，本文聯(lián)合使用Helmholtz偏微分方程[15]和雙曲正切投影[16]進(jìn)行密度過(guò)濾。

Helmholtz偏微分方程表達(dá)式為

(19)

雙曲正切投影表達(dá)式為

(20)

3 拓?fù)鋬?yōu)化數(shù)值實(shí)現(xiàn)

3.1 優(yōu)化方法

拓?fù)鋬?yōu)化問(wèn)題的目標(biāo)函數(shù)及約束函數(shù)一般為設(shè)計(jì)變量的隱式表達(dá)，且拓?fù)鋬?yōu)化問(wèn)題通常具有非凸性和高非線性，計(jì)算量大。當(dāng)前求解拓?fù)鋬?yōu)化問(wèn)題的主要方法分為優(yōu)化準(zhǔn)則法和序列規(guī)劃法兩類。在序列規(guī)劃法中，移動(dòng)漸近線法MMA(method of movingasymptotes)[17]的變量獨(dú)立且只需求解一階導(dǎo)數(shù)就能達(dá)到很好的收斂效果，故本文采用此方法進(jìn)行求解。

3.2 靈敏度分析

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

拉格朗日函數(shù)

(21)

(22)

式中Ψ為目標(biāo)函數(shù)，此處引入如式(22)所示一般化積分形式進(jìn)行表述。θ為拉格朗日乘子，R為控制方程。對(duì)于本文來(lái)說(shuō)，R即為N-S方程、能量守恒方程及結(jié)構(gòu)平衡方程。θ=(ua,pa,Ta)為伴隨變量。

(23)

根據(jù)偏微分方程約束最優(yōu)化問(wèn)題的KKT條件有

(24)

則設(shè)計(jì)變量靈敏度可求解得

(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分別為伴隨速度場(chǎng)和伴隨溫度場(chǎng)，式(24～30)的左端項(xiàng)列入表1。

表1 伴隨方程左端項(xiàng)

由于最小化結(jié)構(gòu)柔度的線彈性力學(xué)問(wèn)題是自伴隨的，即伴隨變量等于結(jié)構(gòu)位移，此處省略該問(wèn)題的伴隨敏度分析。

3.3 基于OpenFOAM進(jìn)行優(yōu)化問(wèn)題求解

本文基于OpenFOAM[19]搭建拓?fù)鋬?yōu)化求解器。流-熱-力耦合拓?fù)鋬?yōu)化流程如圖1所示，首先讀取材料參數(shù)，再對(duì)Navier-Stokes方程以及結(jié)構(gòu)平衡方程進(jìn)行計(jì)算，得到速度場(chǎng)和位移場(chǎng)之后，將速度場(chǎng)帶入能量守恒方程求得溫度場(chǎng)。完成以上初始場(chǎng)計(jì)算后，對(duì)伴隨控制方程進(jìn)行求解，得到各伴隨變量，然后計(jì)算目標(biāo)函數(shù)以及約束函數(shù)對(duì)設(shè)計(jì)變量的靈敏度。最后調(diào)用MMA更新設(shè)計(jì)變量。

圖1 基于OpenFOAM拓?fù)鋬?yōu)化求解流程

4 算 例

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

圖2 對(duì)稱結(jié)構(gòu)尺寸

4.1 不考慮結(jié)構(gòu)承載的流-熱耦合拓?fù)鋬?yōu)化

4.1.1 算例參數(shù)

本文算例的模具冷卻管道不可設(shè)計(jì)域如圖3所示，分別為模具冷卻管道的出入口部分以及模具型腔面。為減少拓?fù)鋬?yōu)化過(guò)程中的計(jì)算量，對(duì)模型施加對(duì)稱邊界條件。冷卻管道入口處法向速度為 0.01 m/s，溫度為0，冷卻管道出口處為絕熱邊界。熱源面施加均布?jí)毫?0 MPa和熱流量40000 W/m2，其他邊界面為絕熱無(wú)滑移面。本算例選用的材料參數(shù)列入表2。

圖3 模具不可設(shè)計(jì)域

圖4 模具工況

表2 材料參數(shù)

4.1.2 傳統(tǒng)冷卻管道設(shè)計(jì)

在注塑成型模具中，傳統(tǒng)的冷卻管道通常為直管設(shè)計(jì)，本文將傳統(tǒng)直管設(shè)計(jì)作為拓?fù)鋬?yōu)化初始構(gòu)型進(jìn)行模具冷卻仿真模擬。初始直管設(shè)計(jì)的流體材料用量百分比為1.4%，模具平均溫度為83.5 ℃，能量耗散為6.861×10-10kg·m2/s3，結(jié)構(gòu)柔度為 0.086 J。因此，在本文優(yōu)化算例中，給定流體材料百分比上限為3%，能量耗散為直管模型的2倍，結(jié)構(gòu)柔度為直管模型的1.5倍。

圖5 傳統(tǒng)直管冷卻溫度場(chǎng)

4.1.3 優(yōu)化結(jié)果

本算例不考慮模具結(jié)構(gòu)承載能力進(jìn)行拓?fù)鋬?yōu)化，其優(yōu)化列式為

MinimizeΨ

(33)

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

圖6 不考慮模具承載能力優(yōu)化結(jié)果

4.2 考慮結(jié)構(gòu)承載的流-熱-力耦合拓?fù)鋬?yōu)化

在4.1節(jié)的優(yōu)化模型基礎(chǔ)上增加考慮結(jié)構(gòu)承載力的結(jié)構(gòu)柔度約束進(jìn)行優(yōu)化設(shè)計(jì)，優(yōu)化列式如式(18)所示。

本算例冷卻管道拓?fù)鋬?yōu)化構(gòu)型如圖7所示，冷卻管道沿著模具型腔面均勻覆蓋。結(jié)構(gòu)溫度場(chǎng)如圖8所示，結(jié)構(gòu)最高溫度出現(xiàn)在模具型腔面臨近出口邊緣處，優(yōu)化結(jié)果顯示此處形成了流道以帶走模具熱量。拓?fù)鋬?yōu)化后的模具平均溫度為76.0 ℃，與傳統(tǒng)直管設(shè)計(jì)相比，拓?fù)鋬?yōu)化后的冷卻管道設(shè)計(jì)使模具的平均溫度降低了9%。

圖7 考慮結(jié)構(gòu)承載能力的拓?fù)鋬?yōu)化管道構(gòu)型

圖8 考慮結(jié)構(gòu)承載能力優(yōu)化后溫度場(chǎng)

5 結(jié) 論

本文基于變密度法進(jìn)行注塑模具冷卻管道拓?fù)鋬?yōu)化設(shè)計(jì)，采用插值函數(shù)描述流固材料分布，以最小結(jié)構(gòu)平均溫度為設(shè)計(jì)目標(biāo)。聯(lián)合使用Helmholtz偏微分方程和雙曲正切投影方程進(jìn)行密度過(guò)濾，得到了清晰的冷卻管道拓?fù)湫螒B(tài)。通過(guò)對(duì)比是否考慮模具承載能力的兩個(gè)算例，得出注塑成型中應(yīng)該將模具承載性能作為評(píng)價(jià)指標(biāo)的結(jié)論。最后基于OpenFOAM平臺(tái)對(duì)注塑成型冷卻系統(tǒng)設(shè)計(jì)問(wèn)題進(jìn)行三維流-熱-力耦合拓?fù)鋬?yōu)化，得到了冷卻效果良好且具有實(shí)際加工意義的流體通道設(shè)計(jì)。與傳統(tǒng)的直管冷卻系統(tǒng)相比，拓?fù)鋬?yōu)化后的冷卻系統(tǒng)能夠使模具平均溫度降低9%。

本文僅是對(duì)于注塑成型冷卻系統(tǒng)設(shè)計(jì)的初步研究，目前還存在諸多不足,如應(yīng)采用周期平均法計(jì)算模具型腔面處的平均壓力和熱源值當(dāng)做優(yōu)化的邊界條件；應(yīng)該采用湍流模型模擬冷卻介質(zhì)的流動(dòng)傳熱性能。未來(lái)將從以上角度對(duì)優(yōu)化進(jìn)行改進(jìn)。

參考文獻(xiàn)(References):

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

[2] 申長(zhǎng)雨,李海梅,高 峰.注射成型技術(shù)發(fā)展概況[J].工程塑料應(yīng)用,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] 王 勇.基于快速隨形熱冷的高光無(wú)痕注塑成型[D].湖北工業(yè)大學(xué),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.