李新磊 吳 坤 趙林英 范學(xué)軍,

* (中國(guó)科學(xué)院力學(xué)研究所高溫氣體動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室,北京 100190)

? (中國(guó)科學(xué)院大學(xué)工程科學(xué)學(xué)院,北京 100049)

** (西北工業(yè)大學(xué)機(jī)電學(xué)院,西安 710072)

?? (合肥中科重明科技有限公司,合肥 230601)

引言

超燃沖壓發(fā)動(dòng)機(jī)作為一種高超聲速巡航飛行器,受氣動(dòng)加熱和燃料的燃燒釋熱影響因而承受著巨大的熱載荷[1-3].在發(fā)動(dòng)機(jī)的燃燒室壁中,通常嵌入矩形直通道構(gòu)成再生冷卻系統(tǒng),以保證發(fā)動(dòng)機(jī)能夠長(zhǎng)時(shí)間穩(wěn)定運(yùn)行.然而,隨著飛行馬赫數(shù)的提升,傳統(tǒng)的再生冷卻結(jié)構(gòu)面臨著冷卻能力不足、結(jié)構(gòu)超溫等問(wèn)題[4-6].因此,亟需開(kāi)展新型熱防護(hù)方案的設(shè)計(jì)及優(yōu)化.

國(guó)內(nèi)外學(xué)者對(duì)此進(jìn)行了大量研究,其中最為常見(jiàn)的一類做法是添置粗糙元件.如Sunden 等[7]通過(guò)在圓管內(nèi)設(shè)置非對(duì)稱的波浪形翅片,使得管道的換熱系數(shù)顯著提高,壁面溫度顯著降低.Li 等[8]對(duì)帶有截?cái)嗬叩木匦瓮ǖ肋M(jìn)行了研究,分析了截?cái)嗬咧車膹?qiáng)化傳熱機(jī)制,并探究了肋高和肋間距等參數(shù)對(duì)通道流動(dòng)換熱性能的影響規(guī)律.楊澤南等[9]對(duì)由Kagome 型網(wǎng)格陣列填充的通道進(jìn)行了研究,發(fā)現(xiàn)在相同工況下,錯(cuò)排網(wǎng)格填充的通道,不僅減重效果明顯,網(wǎng)格后復(fù)雜的渦系結(jié)構(gòu)也使得流體努塞爾數(shù)明顯提高,壁面溫度顯著下降.此外,部分研究還通過(guò)改變?cè)偕鋮s通道的流動(dòng)模式來(lái)達(dá)到增強(qiáng)換熱的目的.如Li 等[10]基于變分法生成了能夠自適應(yīng)幾何和熱邊界條件的變管徑冷卻通道,在犧牲部分壓降的條件下使得壁面平均溫度和溫度不均勻度均有所下降.Li 等[11]對(duì)嵌套式冷卻通道進(jìn)行了研究,發(fā)現(xiàn)同向平行通道有著更好的傳熱性能,且內(nèi)通道的嵌入有助于克服熱分層現(xiàn)象.Zhang 等[12]提出一種雙向流動(dòng)模式,該模式使得壁面溫度分布更加均勻,正癸烷的轉(zhuǎn)化率提升以及較大的化學(xué)熱沉,但壁面最高溫度和壓降均有不同程度的增加.

除了以上研究,近年來(lái)以拓?fù)鋬?yōu)化為代表的新興優(yōu)化方法,逐漸成為散熱系統(tǒng)設(shè)計(jì)的新途徑[13-16].該方法是以微結(jié)構(gòu)的材料性質(zhì)作為設(shè)計(jì)變量,在固定的設(shè)計(jì)空間內(nèi)自動(dòng)尋找最優(yōu)的材料分布方案[17],因此在理論上擁有最高的設(shè)計(jì)自由度,被認(rèn)為是最具前途的優(yōu)化方法之一[18].Tang 等[14]曾對(duì)考慮材料變物性的導(dǎo)熱熱沉結(jié)構(gòu)進(jìn)行拓?fù)鋬?yōu)化,得到了如同樹(shù)枝枝干般的熱沉結(jié)構(gòu),同時(shí)發(fā)現(xiàn)固體材料的變物性顯著改變了拓?fù)鋬?yōu)化構(gòu)型,進(jìn)而影響了其在真實(shí)熱環(huán)境中的性能表現(xiàn).Zhang 等[19]基于偽密度法對(duì)沖壓發(fā)動(dòng)機(jī)的再生冷卻結(jié)構(gòu)進(jìn)行了拓?fù)鋬?yōu)化,發(fā)現(xiàn)其設(shè)計(jì)域內(nèi)生成了眾多垂直于主流方向的微通道,顯著改善了整體換熱,但通道內(nèi)局部的逆流和滯留區(qū)導(dǎo)致出現(xiàn)了局部高溫.本文認(rèn)為,針對(duì)發(fā)動(dòng)機(jī)再生冷卻通道的拓?fù)鋬?yōu)化是結(jié)合了湍流流動(dòng)、燃料物性劇烈變化等特征的復(fù)雜流熱耦合問(wèn)題[19],其計(jì)算量大且極易出現(xiàn)網(wǎng)格依賴性布局、棋盤(pán)格單元等構(gòu)型缺陷,進(jìn)而影響結(jié)構(gòu)性能.針對(duì)上述問(wèn)題,本文將開(kāi)發(fā)一套考慮變熱物性的流熱耦合拓?fù)鋬?yōu)化求解器,采用建表插值法求解熱物性,用以緩解熱物性求解計(jì)算量大的問(wèn)題.此外,將完整再生冷卻設(shè)計(jì)域簡(jiǎn)化為一個(gè)換熱設(shè)計(jì)單元,通過(guò)適當(dāng)?shù)臑V波技術(shù)和松弛求解策略,以盡量避免出現(xiàn)構(gòu)型缺陷.

本文主要內(nèi)容如下,首先建立拓?fù)鋬?yōu)化模型,基于連續(xù)伴隨法對(duì)考慮變物性的伴隨方程和靈敏度進(jìn)行推導(dǎo),并基于開(kāi)源CFD 平臺(tái)OpenFOAM[20]構(gòu)建拓?fù)鋬?yōu)化求解器.針對(duì)流熱耦合問(wèn)題開(kāi)展拓?fù)鋬?yōu)化設(shè)計(jì),重點(diǎn)考察能量耗散約束對(duì)優(yōu)化構(gòu)型的影響.為了評(píng)估拓?fù)鋬?yōu)化構(gòu)型的流動(dòng)換熱性能,本文還將對(duì)拓?fù)鋬?yōu)化構(gòu)型輪廓進(jìn)行提取,開(kāi)展三維共軛換熱數(shù)值模擬分析.

1 拓?fù)鋬?yōu)化模型

拓?fù)鋬?yōu)化方法種類眾多,包括偽密度法、水平集法、相場(chǎng)法以及進(jìn)化結(jié)構(gòu)優(yōu)化法等[18].在各方法中,以偽密度法的應(yīng)用最為廣泛和成熟.因此,本節(jié)將基于偽密度法建立拓?fù)鋬?yōu)化模型,并利用連續(xù)伴隨法對(duì)考慮變物性的伴隨方程及靈敏度進(jìn)行推導(dǎo).

1.1 控制方程

針對(duì)流熱耦合的拓?fù)鋬?yōu)化問(wèn)題,其物理模型示意圖如圖1 所示,設(shè)計(jì)域承受來(lái)自壁面的熱載荷q或由內(nèi)部產(chǎn)生的體積熱Q,流體在流經(jīng)該設(shè)計(jì)域時(shí)與固體實(shí)現(xiàn)熱交換,進(jìn)而達(dá)到熱平衡.該優(yōu)化問(wèn)題的核心是在滿足壓降或材料體積比等約束下,找到最優(yōu)的流固材料單元的空間分布,從而使區(qū)域的平均溫度或最高溫度等目標(biāo)降到最低.

圖1 拓?fù)鋬?yōu)化設(shè)計(jì)示意圖Fig.1 Schematic of the topology optimization design

首先建立流熱耦合拓?fù)鋬?yōu)化模型,該模型的控制方程包含如下的質(zhì)量方程、動(dòng)量方程和能量方程

其中,R(w,γ) 表示等式約束向量,γ 為設(shè)計(jì)變量,取值范圍為0～1,w=(p,u,T) 為隱式依賴于設(shè)計(jì)變量的狀態(tài)向量,在該問(wèn)題中分別對(duì)應(yīng)壓力、速度和溫度;ρ,Cp,μ,λ 和h分別為密度、定壓比熱容、動(dòng)力黏度、熱導(dǎo)率和比焓,均考慮其溫度的相關(guān)性.τ為黏性應(yīng)力張量,其定義式為τ=2μeffS(u)=(μ+μt)·為有效動(dòng)力黏度,表示分子動(dòng)力黏度μ與湍流黏度μt之和,Prt為湍流普朗特?cái)?shù),默認(rèn)值為0.85,Q為體積熱源.

與常規(guī)形式N-S 方程不同的是,拓?fù)鋬?yōu)化模型的動(dòng)量方程中引入了體積力項(xiàng) -α(γ)u,α (γ) 是與多孔介質(zhì)滲透率的倒數(shù)相關(guān)的插值函數(shù).引入設(shè)計(jì)變量 γ 對(duì) α (γ) 進(jìn)行插值,用以區(qū)分流體與固體材料單元,從而實(shí)現(xiàn)相似多孔介質(zhì)流動(dòng)的描述,因此該變量又稱為“偽密度”.插值函數(shù)采用如下形式[21]

式中,θ 為形狀參數(shù),取值范圍為0.005～0.1,該函數(shù)形狀如圖2 所示.αmax的取值則與達(dá)西數(shù)Da相關(guān)

圖2 不同形狀參數(shù)下的插值函數(shù)Fig.2 Interpolation function with different shape parameterθ

式中，l為幾何特征尺寸,Da表征了多孔介質(zhì)的滲透能力,其取值小于 10-5較為合適[15].

相應(yīng)地,在能量方程中對(duì)熱導(dǎo)率進(jìn)行插值處理,λmin和 λmax分別為流體和固體材料的熱導(dǎo)率.當(dāng)設(shè)計(jì)變量 γ=0 時(shí),α 取值為一個(gè)較大數(shù) αmax,表明該處單元為滲透率較低的固體單元;當(dāng)設(shè)計(jì)變量 γ=1 時(shí),α 取值為0,表明該處為流體單元.

本文采用了比焓形式的能量方程,與定壓比熱容存在如下熱力學(xué)關(guān)系

因此,在求解溫度時(shí),可利用如下形式的牛頓迭代法計(jì)算得到

式中,pold,Told,Cp(pold,Told) 和h(pold,Told) 分別為上個(gè)迭代步中的壓力、溫度、比熱和比焓,hnew為利用式(1)計(jì)算得到的更新后的比焓,該迭代在相對(duì)誤差小于10-4時(shí)終止.

除了上述控制方程,本文定義了如下形式的Pnorm 函數(shù)[22]作為該優(yōu)化問(wèn)題的目標(biāo)函數(shù)

該目標(biāo)有利于避免設(shè)計(jì)域內(nèi)局部區(qū)域溫度過(guò)高,式中n取值為30.此外,引入如下的不等式約束函數(shù)

其中,G(w,γ) 為不等式約束向量,分別用來(lái)約束設(shè)計(jì)域的材料體積比例以及進(jìn)出口能量耗散.Vmax為材料體積比例的期望值,φmax為進(jìn)出口能量耗散的期望值.

綜上所述,一個(gè)完整的基于密度法的拓?fù)鋬?yōu)化問(wèn)題可以由以下數(shù)學(xué)表達(dá)式來(lái)概述.

1.2 伴隨方程與靈敏度推導(dǎo)

本節(jié)將基于連續(xù)伴隨法對(duì)流熱耦合問(wèn)題的伴隨關(guān)系式及靈敏度公式進(jìn)行推導(dǎo),首先構(gòu)造如下形式的拉格朗日函數(shù)

式中,Ψ 為目標(biāo)函數(shù),λ=(pa,ua,Ta) 為拉格朗日乘子,又稱伴隨向量.由于R(w,γ)=0,因此原目標(biāo)函數(shù)關(guān)于設(shè)計(jì)變量的靈敏度等價(jià)于如下拉格朗日函數(shù)的全導(dǎo)數(shù)

依據(jù)偏微分方程約束的KKT (Karush–Kuhn–Tucker)條件[23],構(gòu)造如下的伴隨變量方程組

利用積分變換對(duì)上述方程組進(jìn)行化簡(jiǎn),便可得到求解伴隨變量所需的控制方程,其具體形式及推導(dǎo)過(guò)程參見(jiàn)附錄A.能夠發(fā)現(xiàn),由于考慮了溫度相關(guān)的變熱物性和輸運(yùn)性質(zhì),伴隨變量的控制方程形式相比于常物性模型[15]復(fù)雜很多.而伴隨變量控制方程與原始狀態(tài)變量的控制方程存在形式上的相似性,因此可采用相似的數(shù)值算法對(duì)其進(jìn)行求解.

此時(shí)拉格朗日函數(shù)的全導(dǎo)數(shù),即目標(biāo)函數(shù)的靈敏度,經(jīng)推導(dǎo)后轉(zhuǎn)化為如下形式

在求解得到原始狀態(tài)變量與伴隨變量解后,便可計(jì)算得到該靈敏度數(shù)值.

2 數(shù)值方法

本節(jié)將對(duì)拓?fù)鋬?yōu)化過(guò)程中遇到的關(guān)鍵數(shù)值問(wèn)題以及所采取的數(shù)值方法進(jìn)行介紹,包括濾波與投影函數(shù)的使用、熱物性的求解、拓?fù)鋬?yōu)化求解器的構(gòu)建以及共軛換熱數(shù)值求解器的驗(yàn)證等.

2.1 濾波與投影

在拓?fù)鋬?yōu)化模型中,由于對(duì)材料物性采用了特殊的懲罰策略,導(dǎo)致優(yōu)化過(guò)程中容易出現(xiàn)一些不切實(shí)際的布局,如網(wǎng)格依賴性布局或棋盤(pán)格單元等[24-26],這嚴(yán)重影響了材料布局的重構(gòu)和收斂.研究表明,一些正則化技術(shù),如密度濾波器[27]可以有效地解決上述問(wèn)題.本文采用亥姆霍茲型偏微分方程作為濾波器,其方程形式如下

式中 γ0為濾波前的偽密度,γf為濾波后的偽密度,Rmin為濾波半徑.當(dāng)網(wǎng)格單元尺寸為 ?x時(shí),濾波半徑Rmin越大,則 γf的形狀特征越容易被抹平,因此,本文控制各算例的最大濾波半徑不超過(guò) 2 ?x.

此外,采用濾波器通常會(huì)產(chǎn)生灰色單元,即處于流固材料單元之間的模糊單元.因此,可以配合使用投影函數(shù)[28]來(lái)獲得清晰的流固邊界.本文采用如下的雙曲正切函數(shù)來(lái)對(duì) γf進(jìn)行投影

式中 γp為投影后的偽密度,β 和 η 為用來(lái)控制投影函數(shù)形狀的形狀因子.投影函數(shù)的 η 取值為Vmax,用來(lái)保證材料的體積比例,而控制斜率的 β 因子隨優(yōu)化過(guò)程進(jìn)行而逐漸變大,本文中 β 初始值為1,最大值為32,投影后的優(yōu)化變量 γp近乎收斂到接近于{0,1}的集合,可以得到較為清晰的流固邊界.

在對(duì)設(shè)計(jì)變量 γ 進(jìn)行濾波和投影操作后,此時(shí)的靈敏度形式為

為了獲得原始靈敏度,需要進(jìn)行如下的鏈?zhǔn)角髮?dǎo)[15]

2.2 熱物性求解

在發(fā)動(dòng)機(jī)冷卻通道內(nèi),碳?xì)淙剂系臒嵛镄噪S溫度與壓力在很大的范圍內(nèi)變化[29].以正癸烷為例,作為一種典型的碳?xì)淙剂?在再生冷卻的流動(dòng)吸熱過(guò)程中,其自身溫升高達(dá)數(shù)百攝氏度,熱力學(xué)性質(zhì)和輸運(yùn)性質(zhì)會(huì)在臨界點(diǎn)附近發(fā)生陡峭變化,容易引發(fā)方程求解過(guò)程中的數(shù)值積分剛性問(wèn)題.此外,在利用連續(xù)伴隨方法推導(dǎo)得到的伴隨方程中,存在熱力學(xué)性質(zhì)及輸運(yùn)性質(zhì)關(guān)于溫度的偏導(dǎo)項(xiàng),包括 ?h/?T,?ρ/?T,?Cp/?T,? μ/?T和 ? λ/?T,當(dāng)采用解析形式的物性求解方程時(shí),偏導(dǎo)項(xiàng)的求解將變得十分困難.

因此,本文對(duì)物性參數(shù)及各偏導(dǎo)項(xiàng)的計(jì)算采取了建表-插值[30]的方法,首先基于CoolProp 軟件[31]建立數(shù)據(jù)表,在不考慮燃料高溫裂解的前提下,該物性數(shù)據(jù)表以壓力和溫度作為基礎(chǔ)變量,變量間隔分別為1 kPa 和2 K,模擬過(guò)程中采用線性插值計(jì)算得到相關(guān)的參數(shù).圖3 所示為正癸烷在3 MPa 壓力下(Tpc=648.4 K)的密度、比熱以及二者關(guān)于溫度的熱力學(xué)偏導(dǎo)性質(zhì),其中利用建表-插值得到的數(shù)據(jù)與利用CoolProp 直接計(jì)算得到的數(shù)據(jù)十分接近,最大相對(duì)誤差均不超過(guò)0.2%,表明本文所采用的物性計(jì)算方法是準(zhǔn)確的.

圖3 熱物理性質(zhì)及熱力學(xué)偏導(dǎo)Fig.3 Thermo-physical properties and derivatives

2.3 求解器構(gòu)建

基于上述的模型與數(shù)值方法,本文利用開(kāi)源平臺(tái)OpenFOAM 搭建了拓?fù)鋬?yōu)化求解器.圖4 所示為求解器的主要求解框架,其主程序包含以下4 部分:(1)對(duì)原始狀態(tài)變量方程求解;(2)對(duì)伴隨變量方程求解;(3)預(yù)估目標(biāo)函數(shù)與約束值;(4)對(duì)靈敏度和偽密度進(jìn)行濾波和投影.其中對(duì)原始狀態(tài)變量和伴隨變量控制方程的求解模塊均基于內(nèi)置的標(biāo)準(zhǔn)可壓縮流動(dòng)求解器rhoSimpleFoam 進(jìn)行改造.對(duì)熱物性及熱力學(xué)偏導(dǎo)的求解模塊則按上節(jié)所述的建表-插值法編譯為動(dòng)態(tài)鏈接庫(kù).耦合了MMA (method of moving asymptotes)[32]算法作為優(yōu)化工具,同樣將其編譯成庫(kù)供主程序調(diào)用,對(duì)設(shè)計(jì)變量進(jìn)行更新.求解器迭代優(yōu)化的終止條件以目標(biāo)函數(shù)的相對(duì)變化小于0.1%或固定迭代步數(shù)為準(zhǔn)則.

圖4 考慮變物性的拓?fù)鋬?yōu)化求解流程圖Fig.4 Flowchart of the topology optimization approach considering variable physical properties

2.4 共軛換熱數(shù)值驗(yàn)證

本文所開(kāi)展的拓?fù)鋬?yōu)化設(shè)計(jì)均基于二維模型,為了對(duì)其實(shí)際三維構(gòu)型的流動(dòng)換熱性能進(jìn)行評(píng)估,將對(duì)拓?fù)鋬?yōu)化構(gòu)型進(jìn)行提取,并利用OpenFOAM 內(nèi)置的標(biāo)準(zhǔn)共軛換熱求解器chtMultiRegionFoam 進(jìn)行數(shù)值模擬.為了驗(yàn)證該共軛換熱求解器的準(zhǔn)確性,首先與Zhu 等[33]進(jìn)行的流動(dòng)換熱實(shí)驗(yàn)進(jìn)行對(duì)比分析.該實(shí)驗(yàn)中固體材料為不銹鋼,燃料為正癸烷.按照對(duì)應(yīng)的實(shí)驗(yàn)條件,入口給定質(zhì)量流量0.6083 g/s,入口溫度為625.93 K,出口給定背壓4.19 MPa,測(cè)試段體積熱源大小為1.758×108W/m3,連接段熱流邊界的熱流密度大小為11550 W/m2.采用k-? 湍流模型進(jìn)行驗(yàn)證,其結(jié)果如圖5 所示,顯然,該數(shù)值計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果較為吻合,外壁溫及冷卻劑溫度的最大相對(duì)誤差均小于5%,證明該共軛換熱求解器可以較好地預(yù)測(cè)流動(dòng)換熱性能.

圖5 冷卻劑及外壁面溫度的數(shù)值與實(shí)驗(yàn)結(jié)果對(duì)比Fig.5 Comparison of predicted coolant and out wall temperatures against the experimental data

3 結(jié)果與分析

3.1 流熱耦合拓?fù)鋬?yōu)化設(shè)計(jì)

本節(jié)將針對(duì)發(fā)動(dòng)機(jī)再生冷卻通道開(kāi)展流熱耦合拓?fù)鋬?yōu)化設(shè)計(jì),為了節(jié)省計(jì)算成本,以圖6 所示的二維換熱單元作為設(shè)計(jì)域,即不考慮構(gòu)型在z方向的優(yōu)化.其幾何及邊界條件設(shè)置見(jiàn)圖6,計(jì)算域進(jìn)出口尺寸d=5 mm,厚度為2 mm,進(jìn)、出口段為長(zhǎng)度10d的非設(shè)計(jì)域,中間設(shè)計(jì)域長(zhǎng)度為36d,流體材料為正癸烷,固體材料為GH3128,其熱導(dǎo)率隨溫度的變化如下

圖6 流-熱耦合結(jié)構(gòu)計(jì)算域示意圖Fig.6 Schematic of the computational domain of the coupled thermal–fluid structure

設(shè)定入口速度為1 m/s,初始溫度為300 K,對(duì)應(yīng)的入口雷諾數(shù)為2500,計(jì)算域出口壓力設(shè)為3 MPa,其余邊界為對(duì)稱邊界,對(duì)設(shè)計(jì)域施加體積熱源Q=5×108W/m3.利用均勻結(jié)構(gòu)化網(wǎng)格對(duì)區(qū)域進(jìn)行離散化,考慮到要盡可能保留小尺度特征,同時(shí)還要保證加工可制造性,本文將網(wǎng)格尺寸設(shè)定為0.1 mm×0.1 mm,網(wǎng)格總數(shù)為140000,采用k-ε湍流模型進(jìn)行模擬.

該優(yōu)化模型以P-norm 函數(shù) Ψ 作為目標(biāo),施加材料體積約束gvol,Vmax取值為0.68,αmax取初始值為103,每步增長(zhǎng)1.05 倍,最大值設(shè)為108.為了考察不同能量耗散約束值對(duì)拓?fù)鋬?yōu)化構(gòu)型的影響,本文首先對(duì) γ0=0.5,αmax=103的初始分布構(gòu)型進(jìn)行計(jì)算,得到其能量耗散值約為0.015 k g·m/s2,隨后將該數(shù)值作為參考值,對(duì) φmax的取值范圍在0.02～0.25 內(nèi)的算例進(jìn)行了優(yōu)化,并從中提取出 φmax為0.03,0.05,0.14,0.19 和0.21 的5 個(gè)算例,分別命名為Case 1～Case 5.

圖7 為優(yōu)化得到的偽密度場(chǎng)分布及對(duì)應(yīng)的速度分布與溫度分布.為了方便觀察,這里將計(jì)算結(jié)果按關(guān)于x軸對(duì)稱形式顯示.能夠發(fā)現(xiàn),在材料比例基本一致的條件下,Case 1～Case 5 的固體域均發(fā)生了不同程度的分裂,且各分裂的固體胞元塊呈現(xiàn)出不同的排列方式.在Case 1 中,大部分固體域分布在通道兩側(cè),流體流通面積較大,流速較低.而隨著能量耗散約束值增大,固體域分裂成的塊數(shù)逐漸增多,流體域內(nèi)的微細(xì)通道數(shù)量增多,流體流通面積減小,流速增大,與固體域的接觸面積明顯增加.從Case 1～Case 5,流體的流動(dòng)分離與再混合行為逐漸增多,多處區(qū)域存在局部的流動(dòng)加速,同樣促進(jìn)了流體與固體接觸面的熱量交換.

圖7 拓?fù)鋬?yōu)化構(gòu)型及狀態(tài)變量云圖Fig.7 Contours of the optimized layouts and state variables

圖7 拓?fù)鋬?yōu)化構(gòu)型及狀態(tài)變量云圖 (續(xù))Fig.7 Contours of the optimized layouts and state variables (continued)

通過(guò)觀察各偽密度法場(chǎng)對(duì)應(yīng)的溫度云圖能夠發(fā)現(xiàn),盡管采用了投影等數(shù)值策略,在優(yōu)化構(gòu)型中仍存在部分灰色單元區(qū)域,如在Case 5 的后段,由于灰色單元的存在,導(dǎo)致該區(qū)域中出現(xiàn)流動(dòng)死區(qū),其溫度呈現(xiàn)出非連續(xù)變化,這在一定程度上干擾了對(duì)流動(dòng)換熱性能的評(píng)估.因此,為了評(píng)估真實(shí)的拓?fù)鋬?yōu)化通道的流動(dòng)換熱性能,本文將進(jìn)一步對(duì)拓?fù)鋬?yōu)化構(gòu)型進(jìn)行提取,將其對(duì)應(yīng)的三維結(jié)構(gòu)進(jìn)行共軛換熱數(shù)值模擬.

3.2 拓?fù)淅鋮s通道的流動(dòng)換熱性能

本節(jié)將對(duì)三維拓?fù)鋬?yōu)化通道開(kāi)展共軛換熱數(shù)值模擬,并對(duì)其流動(dòng)換熱性能進(jìn)行分析.首先將5 種拓?fù)鋬?yōu)化構(gòu)型中 γ <0.5 的區(qū)域進(jìn)行提取,利用樣條曲線對(duì)拓?fù)浒M(jìn)行建模,圖8 所示為5 種拓?fù)渫ǖ?Case 1～Case 5)的三維模型,此外,加入了傳統(tǒng)直通道構(gòu)型Case 0 作為對(duì)比算例,其固體壁肋寬約為1.70 mm,高度為2 mm,冷熱壁厚度均為0.75 mm.

圖8 拓?fù)鋬?yōu)化通道的三維幾何視圖Fig.8 Three-dimensional layouts of topology optimization channels

隨后對(duì)上述模型進(jìn)行計(jì)算,其入口速度仍為1 m/s,初始溫度為300 K,出口壓力為3 MPa.按能量守恒原則,將體積熱源Q換算為在固體域底部面施加的面熱流q,其大小約為1 MW/m2,上下層固體域左右面為對(duì)稱邊界,其余壁面為絕熱邊界,采用kε湍流模型并結(jié)合壁面函數(shù),對(duì)通道內(nèi)各固體壁面添加邊界層網(wǎng)格,第一層網(wǎng)格高度為0.004 mm,以保證無(wú)量綱距離y+<10.

圖9 所示為計(jì)算得到的各通道中心截面和底部被加熱面的溫度分布云圖.能夠發(fā)現(xiàn),由于固體胞元的不規(guī)則排列,5 個(gè)通道的熱壁面溫度呈現(xiàn)非均勻分布,在某些較大的胞元間隙位置,如Case 3 的190 mm 和210 mm 位置處,由于流體流通面積較大,流速較低,導(dǎo)致底面加熱位置的壁溫較高,但從Case 1到Case 5,各通道內(nèi)的高溫區(qū)域顯著減少.

圖9 xy 截面的溫度分布云圖Fig.9 Contours of temperature on xy plane

圖10 為中心xy截面上的速度云圖、壓力云圖以及濕周底部面的努塞爾數(shù)云圖.在Case 1～Case 5中,冷卻劑平均流速逐漸增大,壓力損失也逐漸增大.在Case 3～Case 5 中,通道前端位置的固體胞元首先起到了節(jié)流的作用,使冷卻劑流速加快,隨后在流動(dòng)過(guò)程中發(fā)生了大量的流動(dòng)分離和再混合,帶來(lái)了較大的能量耗散和壓力損失.但由于冷卻劑流動(dòng)邊界層和熱邊界層的分離和再附著,使得傳熱熱阻降低,增強(qiáng)了流體與固體間的換熱,局部努塞爾數(shù)增大.此外,在多個(gè)固體胞元的鈍前緣位置,局部努塞爾數(shù)也有明顯提高,這表明流體對(duì)固體的沖擊效應(yīng),同樣使得胞元前緣位置的傳熱增強(qiáng).

表1 所示為各拓?fù)淅鋮s通道的流動(dòng)和換熱性能數(shù)據(jù),圖11 為熱壁面平均溫度和最高溫度與通道壓降的性能關(guān)系曲線.能夠發(fā)現(xiàn),從Case 1～Case 5,各通道的壓降逐漸增加,且熱壁面平均溫度和最高溫度均隨壓降升高而顯著降低,其中Case 5 熱壁面平均溫度和最高溫度分別降低了91.1 K 和55.1 K.各通道的平均努塞爾數(shù)逐漸增加,表明各通道整體的換熱能力逐漸增強(qiáng).但相比于Case 0,僅Case 3～Case 5 起到了強(qiáng)化換熱的效果,其平均努塞爾數(shù)增益百分比分別為12.6%,16.0%和23.4%.

表1 6 種通道的流動(dòng)換熱性能Table 1 Flow and heat transfer performances of six channels

圖11 各通道熱壁面平均溫度和最高溫度隨壓降的性能關(guān)系曲線Fig.11 Relationship between the average and maximum temperature in terms of the pressure drop of the cooling channels

3.3 拓?fù)淅鋮s通道的強(qiáng)化換熱機(jī)理

圖12 所示為各通道在流向不同位置處yz截面上的渦量云圖.能夠觀察到,相比于Case 0～Case 2,Case 3～Case 5 通道內(nèi)冷卻劑在多個(gè)位置呈現(xiàn)二次渦形態(tài),如在Case 4 和Case 5 的x=70～85 mm、Case 5 的105～115 mm 等位置.該渦系結(jié)構(gòu)的形成是由于流體在流經(jīng)固體胞元間隙時(shí),受其他分支的流體摻混與擾動(dòng)而形成,其渦系形態(tài)一直保持到胞元后緣,但沿流向渦強(qiáng)度逐漸減弱,直到下一個(gè)固體域分裂的區(qū)域,流體再次形成明顯的二次渦結(jié)構(gòu).該渦系結(jié)構(gòu)促使近壁面流體與主流發(fā)生摻混,促進(jìn)熱交換.

圖12 不同流向位置yz 截面上的渦量云圖Fig.12 Contours of vorticity on yz planes at different streamwise locations

圖13 所示為各拓?fù)淅鋮s通道在xy中心截面上冷卻劑的湍動(dòng)能分布云圖.在Case 3～Case 5 通道內(nèi),多個(gè)固體胞元的側(cè)壁以及尾部回流區(qū)等位置,由于局部的流動(dòng)加速或流動(dòng)混合,使流體的湍動(dòng)能得以激發(fā),換熱能力增強(qiáng).此外,流體湍動(dòng)能增強(qiáng)區(qū)域與二次渦形成區(qū)域基本吻合,表明二次渦結(jié)構(gòu)有助于激發(fā)湍動(dòng)能,進(jìn)而增強(qiáng)換熱.

圖13 xy 截面湍動(dòng)能云圖Fig.13 Contours of turbulence energy on xy plane

4 結(jié)論

本文針對(duì)發(fā)動(dòng)機(jī)的再生冷卻結(jié)構(gòu),開(kāi)展了考慮變物性的流熱耦合拓?fù)鋬?yōu)化設(shè)計(jì),并對(duì)三維拓?fù)浣Y(jié)構(gòu)的流動(dòng)換熱性能進(jìn)行了數(shù)值模擬與分析,得到主要結(jié)論如下.

(1)在基于連續(xù)伴隨法的求解框架中,考慮工質(zhì)的變物性會(huì)使伴隨方程的形式更加復(fù)雜,求解更加困難.利用建表-插值法對(duì)方程中熱物性、輸運(yùn)性質(zhì)及各偏導(dǎo)項(xiàng)進(jìn)行求解,是一種有效且準(zhǔn)確的方法.

(2)通過(guò)對(duì)流熱耦合結(jié)構(gòu)進(jìn)行拓?fù)鋬?yōu)化設(shè)計(jì),發(fā)現(xiàn)其優(yōu)化構(gòu)型隨著能量耗散約束值的增大而愈加復(fù)雜,且伴有多重的流動(dòng)分離和再混合現(xiàn)象.

(3)通過(guò)對(duì)三維拓?fù)鋬?yōu)化構(gòu)型進(jìn)行共軛換熱數(shù)值模擬,結(jié)果表明,相比于直通道構(gòu)型Case 0,拓?fù)鋬?yōu)化冷卻通道內(nèi)的冷卻劑流態(tài)復(fù)雜,存在二次渦等渦系結(jié)構(gòu),在帶來(lái)更大壓力損失的同時(shí),各通道局部的換熱效果得到增強(qiáng).在5 個(gè)拓?fù)渫ǖ乐?Case 1～Case 5),Case 3～Case 5 整體實(shí)現(xiàn)了強(qiáng)化換熱的效果,其平均努塞爾數(shù)增益百分比分別為12.6%,16.0%和23.4%.

綜上所述,面向發(fā)動(dòng)機(jī)再生冷卻的拓?fù)鋬?yōu)化設(shè)計(jì)具有良好的應(yīng)用前景,后續(xù)將繼續(xù)開(kāi)展更多三維流熱耦合結(jié)構(gòu)的設(shè)計(jì)工作,并逐步完善求解器的其他功能.

附錄A

如1.2 節(jié)所述,伴隨變量方程組滿足正文關(guān)系式(12),在對(duì)其進(jìn)行化簡(jiǎn)前,首先將目標(biāo)函數(shù)的定義域分為邊界和內(nèi)部域兩部分

則目標(biāo)函數(shù)關(guān)于狀態(tài)變量的偏導(dǎo)數(shù)為

式中wd表示目標(biāo)關(guān)于狀態(tài)變量的導(dǎo)數(shù)的方向,變量p,u和T對(duì)應(yīng)的導(dǎo)數(shù)方向分別為pd,ud和Td. 下面對(duì)方程組(12)中的各積分項(xiàng)進(jìn)行積分變換,結(jié)果如下

此時(shí),拉格朗日函數(shù)關(guān)于狀態(tài)變量的偏導(dǎo)數(shù)可以推導(dǎo)為如下形式

若要滿足正文式(12),可得到伴隨變量控制方程

以及對(duì)應(yīng)的邊界條件方程

對(duì)于方程(A12)～式(A14),本文將其分別作為伴隨變量ua,pa和Ta的顯式邊界條件,其數(shù)值隨優(yōu)化迭代而更新.方程中所涉及的目標(biāo)函數(shù)關(guān)于狀態(tài)變量的各項(xiàng)偏導(dǎo)數(shù)表達(dá)式如表A1所示.此時(shí)化簡(jiǎn)后的靈敏度表達(dá)式如下所示

對(duì)于能量耗散約束函數(shù)的靈敏度,其推導(dǎo)過(guò)程與換熱目標(biāo)函數(shù)靈敏度相似,本文不予以贅述.

附表1 目標(biāo)及約束函數(shù)關(guān)于狀態(tài)變量的偏導(dǎo)數(shù)Table 1 Derivatives of the objectives and constraints w.r.t state variables