陳金月，王 坤，閔春華，楊旭光

（河北工業(yè)大學(xué) 能源與環(huán)境工程學(xué)院，天津 300401）

0 引言

通道表面增設(shè)擾流橫肋是無(wú)源強(qiáng)化換熱的一種常用手段。因其具有操作簡(jiǎn)單、強(qiáng)化效果較好等優(yōu)點(diǎn)，在電子器件散熱[1]、冷卻設(shè)備[2]及太陽(yáng)能空氣加熱器[3]等領(lǐng)域均得到了廣泛應(yīng)用。擾流橫肋的形狀和排布方式是影響強(qiáng)化換熱的重要因素。已有很多學(xué)者對(duì)此開(kāi)展了相關(guān)研究工作。表1給出了文獻(xiàn)中出現(xiàn)的幾種不同形狀的擾流橫肋結(jié)構(gòu)。

表1 文獻(xiàn)中典型橫肋擾流橫肋匯總Tab.1 Summary of typical structures of transverse ribs

Vanaki 等[9]對(duì)矩形、三角形、楔形等不同擾流橫肋形狀進(jìn)行了優(yōu)化和性能對(duì)比，結(jié)果表明三角形肋能夠?qū)崿F(xiàn)最大綜合強(qiáng)化換熱效果；Yadav等[10]對(duì)太陽(yáng)能空氣加熱器內(nèi)的流動(dòng)與換熱進(jìn)行了數(shù)值分析，研究了湍流狀態(tài)下的12種不同尺寸的等邊三角形剖肋構(gòu)型，結(jié)果表明肋片存在最優(yōu)的螺距與高度。除肋片形狀外，排布方式對(duì)加肋通道的性能影響也很大，以下是研究者對(duì)不同排布方式的肋片的研究。Vanaki 等[9]對(duì)湍流狀態(tài)下，上底面加熱的2 種布置方式（直排式和交錯(cuò)式）進(jìn)行數(shù)值模擬對(duì)比，結(jié)果表明直排式的Nu值最高，綜合強(qiáng)化換熱效果最好；Liu 等[11]重點(diǎn)研究了不同截?cái)喾绞嚼咂鞯赖膫鳠峒白枇μ匦?，?duì)比研究發(fā)現(xiàn)，截?cái)嘈屠咂诓唤档蛽Q熱效果的基礎(chǔ)上，可以大大減小流道的壓降；Abraham等[12]對(duì)湍流狀態(tài)下V型和W型肋的中心頂點(diǎn)分別指向上游和下游的4種布置方式進(jìn)行了實(shí)驗(yàn)和模擬研究，最終得出建議使用W型排布肋。

以上研究針對(duì)不同形狀和排布的擾流橫肋的性能進(jìn)行了對(duì)比分析，未對(duì)橫肋的形狀進(jìn)行優(yōu)化。近年來(lái)反演優(yōu)化的思想越來(lái)越受到大家的重視。早期的研究主要關(guān)注于導(dǎo)熱反問(wèn)題[13-15]或是反演未知的熱邊界條件[16-18]，而將其應(yīng)用于對(duì)流換熱幾何參數(shù)反演優(yōu)化的還較少。在本研究團(tuán)隊(duì)之前的研究工作中[19-20]，對(duì)于布置方形肋片的二維流道進(jìn)行了對(duì)流換熱反問(wèn)題的研究，分別以流道出口溫度（Tout定壁溫邊界條件）和火積耗散（Eh定熱流邊界條件）作為目標(biāo)函數(shù)，反演優(yōu)化了流道內(nèi)肋片的最佳節(jié)距值。然而，這些研究工作僅對(duì)傳熱性能進(jìn)行了優(yōu)化，并未考慮阻力特性。本文基于共軛梯度法，以綜合考慮傳熱與阻力特性的綜合換熱性能作為目標(biāo)函數(shù)，對(duì)不同排布方式下矩形通道內(nèi)擾流橫肋的形狀進(jìn)行反演優(yōu)化。

1 模型及計(jì)算方法

1.1 計(jì)算模型

圖1 a）是三維矩形通道，圖1 b）是簡(jiǎn)化的二維模型。主流道幾何尺寸為：高H=20 mm，長(zhǎng)L=12.5H，寬W=2H。計(jì)算過(guò)程中為保證進(jìn)口流體速度穩(wěn)定及出口無(wú)回流，流道的實(shí)際計(jì)算長(zhǎng)度分別向兩端延伸了l，其中l(wèi)=5H。初始擾流橫肋結(jié)構(gòu)為等腰直角三角形橫肋h(yuǎn)=1/4H，w=1/2H，第一個(gè)三角形橫肋距流道入口端的距離為11.25H，S為擾流橫肋頂點(diǎn)之間的距離，擾流橫肋個(gè)數(shù)用N表示，圖中N=4。

圖1 計(jì)算模型示意圖Fig.1 Schematic diagram of model

流道中的流動(dòng)介質(zhì)為空氣，假設(shè)空氣的物性參數(shù)為常數(shù)，通道左側(cè)為入口，流體以恒定流速流入，給定雷諾數(shù)為8 000，進(jìn)口溫度Tin=300 K，右側(cè)為出口，出口截面為壓力出口。主流道底面為第二類(lèi)邊界條件，q=5 000 W/m2，其他邊界采用絕熱邊界條件，忽略該系統(tǒng)的輻射換熱。

1.2 計(jì)算方法

壓力和速度的耦合項(xiàng)采用SIMPLE 算法，無(wú)滑移邊界條件，對(duì)有限體積法中控制體進(jìn)行離散化時(shí)，對(duì)流項(xiàng)采用二階迎風(fēng)格式離散，擴(kuò)散項(xiàng)采用中心差分格式離散，湍流模型選擇standardk-ε。正問(wèn)題的收斂準(zhǔn)則是當(dāng)連續(xù)方程、動(dòng)量方程的殘差均小于10-6，能量方程的殘差小于10-7時(shí)。模型的處理、網(wǎng)格的劃分均采用ANSYS-FLUENT軟件實(shí)現(xiàn)。

流動(dòng)與換熱過(guò)程的控制方程如下。

連續(xù)性方程為

式中，ρ為流體的密度。

動(dòng)量守恒方程為

式中，

湍流動(dòng)能方程為

耗散率方程為

式中：S為應(yīng)變張量的平均速率之模，其定義式為S=(2SijSij)1/2，Sij=1/2(?uj/?xi+?uj/?xi)；αk和αε分別為k和ε的有效普朗特?cái)?shù)（Prandtl number）的倒數(shù)，其值αk=1.0，αε=1.3；C1ε=1.42，C2ε=1.68；應(yīng)變率，η=Sk/ε，η0=4.38，β=0.012。

能量守恒方程為

式中，λ是導(dǎo)熱系數(shù)。

1.3 模型與網(wǎng)格驗(yàn)證

對(duì)于計(jì)算區(qū)域進(jìn)行網(wǎng)格獨(dú)立性驗(yàn)證時(shí)，矩形通道采用結(jié)構(gòu)化網(wǎng)格，擾流元左右20 mm附近采用非結(jié)構(gòu)性網(wǎng)格加密，網(wǎng)格數(shù)大于109 942后，流道的Nu及f均不發(fā)生大的變化，其結(jié)果如圖2所示。因此，選擇此時(shí)的網(wǎng)格數(shù)進(jìn)行數(shù)值模擬計(jì)算。

圖2 網(wǎng)格無(wú)關(guān)性考核Fig.2 Grid independence test

在進(jìn)行優(yōu)化之前，為了驗(yàn)證模型的準(zhǔn)確性，將擾流橫肋周?chē)木植颗麪枖?shù)數(shù)值結(jié)果與Liou等[21]的實(shí)驗(yàn)結(jié)果進(jìn)行了比較，其結(jié)果如圖3所示?？梢钥闯?，模擬結(jié)果與實(shí)驗(yàn)結(jié)果具有很好的吻合度，驗(yàn)證了數(shù)值模型的可靠性。

圖3 模型驗(yàn)證Fig.3 Model validation

2 反演優(yōu)化方法

2.1 優(yōu)化問(wèn)題的目標(biāo)函數(shù)

上述的擾流橫肋幾何邊界反問(wèn)題的目的是找到擾流元幾何邊界r（θ），使目標(biāo)函數(shù)值最大，目標(biāo)函數(shù)可以定義為

式中：a為系統(tǒng)傳熱能力與系統(tǒng)阻力的權(quán)重系數(shù)，其大小反映了優(yōu)化過(guò)程中系統(tǒng)換熱能力及減阻效果的側(cè)重程度；Nu0和f0分別代表光滑通道的努塞爾數(shù)和摩擦系數(shù)，由Gnielinski和Filonenko公式定義：

2.2 擾流橫肋邊界幾何形狀的擬合

本文采用三次Spline函數(shù)將離散的數(shù)據(jù)點(diǎn)連成光滑曲線，與其他方法相比，此方法操作性強(qiáng)，光滑性好。設(shè)有m個(gè)邊界點(diǎn)(θi,ri)(i=1,2,…)，ri(i=1,2,3,…)為極半徑，ri值可變；θi為極角，(0＜i≤m)，下角標(biāo)i表示沿順時(shí)針?lè)较虻碾x散點(diǎn)序號(hào)。由插值函數(shù)的邊界條件可得

則可以唯一確定一個(gè)φΔ(θ)表示成如下的分段形式：

式中，Sk(θ)(1 ≤k≤m)是1個(gè)三次多項(xiàng)式，且滿足插值條件

由式（13）可得

其系數(shù)的確定過(guò)程為：

對(duì)式（14）兩邊同時(shí)求導(dǎo)得，

利用αk(θ)和βk(θ)的性質(zhì)，式（15）等價(jià)于

式中，hk=θk+1-θk。

由邊界條件公式（11）可導(dǎo)出2個(gè)方程：

由式（17）和式（18）得矩陣形式的線性方程組為

由此可得到m0,m1,···,mm-1,mm。

2.3 簡(jiǎn)化共軛梯度法求解反問(wèn)題的計(jì)算流程

本文所用共軛梯度法構(gòu)造一組搜索方向，求出目標(biāo)函數(shù)的極值點(diǎn)。其求解的計(jì)算流程如下。

每個(gè)元素ri的更新方式為

式中：βi為搜索步長(zhǎng)，取搜索步長(zhǎng)為定值，其值為0.0002mm[19]；di為搜索方向，其計(jì)算式（21）為

反問(wèn)題求解的計(jì)算過(guò)程如下：

1）采用數(shù)值模擬的方法，求解整個(gè)計(jì)算區(qū)域的速度和溫度分布，然后求得目標(biāo)函數(shù)值；

2）根據(jù)正問(wèn)題的結(jié)果和收斂準(zhǔn)則ε，（ε是很小的正數(shù)），此時(shí)令k=1；

3）給研究變量一個(gè)小的擾動(dòng)Δr，根據(jù)數(shù)值計(jì)算結(jié)果，按照式（8）和式（23）計(jì)算對(duì)應(yīng)的目標(biāo)函數(shù)J[r(θ)]及其梯度，判斷梯度值是否滿足收斂準(zhǔn)則，若≤ε，則迭代終止，否則繼續(xù)；

4）按照式（22）和式（21）分別計(jì)算共軛系數(shù)γi和搜索方向di；

5）根據(jù)式（20）計(jì)算更新ri的值，并令n=n+1，返回步驟3。

3 結(jié)果與討論

3.1 不同目標(biāo)函數(shù)的對(duì)比

以具有單個(gè)擾流橫肋的通道為例，初始擾流橫肋形狀為圖1所示結(jié)構(gòu)，此次迭代為未擬合邊界所得結(jié)果，離散邊界點(diǎn)為11，其他邊界條件不變。Fan等[22]曾提出一種傳熱綜合效果評(píng)價(jià)圖，給出了等流量約束條件下、等壓降約束條件下、等泵功約束條件下分別為a=1、1/2、1/3。根據(jù)式（8），目標(biāo)函數(shù)分別為。

圖4和圖5分別為不同a值反演優(yōu)化后所對(duì)應(yīng)的優(yōu)化結(jié)果。由圖可知，隨著a值逐漸增大，即隨著系統(tǒng)阻力權(quán)重的逐漸增大，優(yōu)化后擾流橫肋結(jié)構(gòu)的尖角逐漸平緩甚至消失。這是由于尖角的存在可以提高傳熱性能，但同時(shí)也會(huì)增大流動(dòng)阻力。當(dāng)a=1/3時(shí)，目標(biāo)函數(shù)值達(dá)到1.869，但增量最小，僅為3.32%。可根據(jù)不同的工程要求來(lái)選取不同的a值。

圖4 不同a 值下的優(yōu)化形狀Fig.4 Optimal shape with different a

圖5 不同a 值下優(yōu)化結(jié)果Fig.5 Optimization results under different a

3.2 不同排布方式的對(duì)比

以具有2 個(gè)擾流元的通道為例（N=2），初始擾流元形狀為圖1 所示結(jié)構(gòu)，此時(shí)目標(biāo)函數(shù)為，其他邊界條件不變。圖6為擾流元不同排布方式的結(jié)構(gòu)示意圖，之后分別對(duì)4 種布置方式下擾流元的間距S進(jìn)行討論，得出底排隨著S增大J值一直在減小，最終取S=10 mm；叉排后置S增大的過(guò)程中，出現(xiàn)J值先增大后變小，當(dāng)S=65 mm 時(shí)，此時(shí)J值最大；叉排前置隨著S增大J值一直增大，當(dāng)S=120 mm 時(shí)，此時(shí)擾流元已經(jīng)位于加熱壁面的邊緣，因此取此時(shí)的位置對(duì)擾流元進(jìn)行反演優(yōu)化，以下均采用四種間距下的初始結(jié)構(gòu)進(jìn)行優(yōu)化。

圖6 不同排布方式的結(jié)構(gòu)（N=2）Fig.6 Schematic diagram of different arrangement（N=2）

圖7中4種擾流元布置方式，可以看出布置方式為底排的擾流元優(yōu)化后所得的目標(biāo)函數(shù)最大；順排布置的擾流元優(yōu)化后所得的目標(biāo)函數(shù)與底排布置的擾流元優(yōu)化結(jié)果相比稍??；叉排布置的2 種擾流元優(yōu)化后所得的目標(biāo)函數(shù)遠(yuǎn)小于底排布置的結(jié)構(gòu)。

圖7 不同排布方式優(yōu)化性能隨迭代步數(shù)的變化（N=2）Fig.7 Optimization of the performance-changing process with different different arrangement（N=2）

3.3 不同擾流橫肋個(gè)數(shù)的對(duì)比

3.3.1 邊界離散點(diǎn)數(shù)目的影響

以具有單個(gè)擾流元的通道為例，初始擾流元形狀為圖1 所示結(jié)構(gòu)，此時(shí)目標(biāo)函數(shù)為，其他邊界條件不變。分別對(duì)離散邊界點(diǎn)個(gè)數(shù)為m=3、7、11 的3 種情況進(jìn)行反演。表2 給出了擬合前后取不同離散邊界點(diǎn)優(yōu)化結(jié)果的速度與壓力分布?？梢钥闯觯瑪M合前后大致形狀基本類(lèi)似。

表2 擬合前后云圖對(duì)比Tab.2 Contour plot before and after fitting

圖8 給出了擬合后不同離散邊界點(diǎn)所得的目標(biāo)函數(shù)值，m=3時(shí)，J值達(dá)到1.866，此時(shí)的目標(biāo)函數(shù)值遠(yuǎn)高于多個(gè)點(diǎn)所取得的值，分別比m=7 時(shí)增加0.865%，比m=11 時(shí)增加0.92%。因此，選用m=3的反演結(jié)果。對(duì)于擬合，并非所選點(diǎn)越多，擬合階次越高越能達(dá)到我們的要求，傳統(tǒng)的高階曲線擬合方法，如拉格朗日法會(huì)在擬合曲線的兩端出現(xiàn)大幅震蕩的“龍格”現(xiàn)象，造成求解區(qū)域的嚴(yán)重失真，反演過(guò)程無(wú)法實(shí)現(xiàn)。

圖8 不同點(diǎn)數(shù)下的優(yōu)化結(jié)果Fig.8 Optimization results under different m

圖9 為不同初始擾流橫肋形狀示意圖，分別為：等腰直角三角形橫肋、正五邊形、半圓形結(jié)構(gòu)，主要尺寸w，h與等腰直角三角形相同，3個(gè)初始形狀在離散邊界點(diǎn)為3時(shí)進(jìn)行反演優(yōu)化。圖10為離散邊界點(diǎn)為3時(shí)不同初始形狀下的優(yōu)化結(jié)果，擬合優(yōu)化后形狀基本類(lèi)似，如半個(gè)水滴形，最高點(diǎn)像左偏移，文獻(xiàn)[23]中曾模擬了擾流元柱的流動(dòng)和換熱，優(yōu)化得到插排布置的擾流元柱為水滴形，與本文離散邊界點(diǎn)為3時(shí)反演優(yōu)化得到的結(jié)果類(lèi)似。

圖9 不同初始擾流元形狀Fig.9 Different initial rib shapes

圖10 不同初始形狀下的優(yōu)化結(jié)果Fig.10 Optimization results under different initial shapes

3.3.2 擾流橫肋個(gè)數(shù)的影響

取圖1所示的三角形為初始結(jié)構(gòu)，邊界離散點(diǎn)為3，目標(biāo)函數(shù)為，其他邊界條件不變。圖11為不同擾流橫肋個(gè)數(shù)擬合后的目標(biāo)函數(shù)值的比較，由圖11可知，不同擾流元個(gè)數(shù)擬合所得結(jié)果，N=1與N=3 所得目標(biāo)函數(shù)結(jié)果相差不大，能達(dá)到1.866，N=1 時(shí)，減少2.84%的換熱，降低了18.47%的阻力，得到2.89%的目標(biāo)函數(shù)的提升；N=3時(shí)，增加了2.6%的換熱，降低了2.12%的阻力，得到3.33%的目標(biāo)函數(shù)的提升。

圖11 不同擾流元個(gè)數(shù)下的優(yōu)化結(jié)果Fig.11 Optimization results under different number of ribs

圖12 為不同個(gè)數(shù)擾流橫肋的速度等值線圖，可以看出，N=1 時(shí)擾流元的形狀為頂點(diǎn)向左偏移的半個(gè)水滴形；N=2 時(shí)，為2 個(gè)向內(nèi)凹的三角形結(jié)構(gòu)；N=3 時(shí)，前2 個(gè)為向內(nèi)凹的三角形結(jié)構(gòu)，第3 個(gè)結(jié)構(gòu)為頂點(diǎn)較低且向左偏移的半個(gè)水滴形；N=4時(shí)，前2個(gè)同為向內(nèi)凹的三角形結(jié)構(gòu)，第3個(gè)為較低的三角形結(jié)構(gòu)，第4 個(gè)為頂點(diǎn)較低且向左偏移的半個(gè)水滴形，N=5 時(shí)，此時(shí)目標(biāo)函數(shù)值較小，不與考慮。這是因?yàn)镹=3 時(shí)的結(jié)構(gòu)是由N=1 和N=2 結(jié)構(gòu)的組合，初始形狀需對(duì)來(lái)流的擾動(dòng)盡可能大，所以開(kāi)始為內(nèi)凹三角形，最后1個(gè)需要減小由于三角形的強(qiáng)擾動(dòng)引起的阻力，所以結(jié)構(gòu)會(huì)較平緩，最終會(huì)出現(xiàn)這種形狀的多個(gè)擾流元的組合結(jié)構(gòu)。

圖12 不同擾流元個(gè)數(shù)下的速度云圖Fig.12 Velocity contour with different numbers of rib transverse

由于N=1與N=3所得目標(biāo)函數(shù)結(jié)果相差不大，所以我們觀察在流動(dòng)過(guò)程中，如圖13所示，不同擾流元個(gè)數(shù)沿x軸局部Nux和fx變化，當(dāng)氣體流到擾流元之前，矩形通道內(nèi)氣體的局部阻力損失相同，Re一定時(shí)，流體剛進(jìn)入通道時(shí)，由于入口段的邊界層較薄，局部Nux數(shù)較高，隨著流動(dòng)的進(jìn)行，入口效應(yīng)逐漸減弱，Nux數(shù)快速降低，當(dāng)流體流經(jīng)擾流元結(jié)構(gòu)時(shí)，在擾流元尾部，在壓差的作用下，產(chǎn)生了二次流，加強(qiáng)了擾動(dòng)和混合，Nux數(shù)有所增加，之后Nux會(huì)趨于穩(wěn)定，N=1與N=3的Nux數(shù)基本相似，但沿著通道到達(dá)擾流元尾部時(shí)，N=3時(shí)的局部阻力值遠(yuǎn)小于N=1時(shí)的局部阻力，在二次流完全結(jié)束之后，兩者的局部阻力值基本相同。因此，最終選用N=3時(shí)的優(yōu)化結(jié)構(gòu)。

圖13 不同擾流橫肋個(gè)數(shù)沿x 軸局部Nux 和fx 變化Fig.13 The number of different rib transverse varies locally along the X-axis in Nux and fx

4 結(jié)論

本文以綜合換熱性能最佳為目標(biāo)，采用共軛梯度法對(duì)矩形通道內(nèi)擾流橫肋的最佳形狀進(jìn)行了反演優(yōu)化，對(duì)優(yōu)化目標(biāo)函數(shù)的選取、擾流橫肋布置方式及數(shù)目的影響進(jìn)行了討論分析。得到如下結(jié)論：

1）采用不同的綜合換熱性能評(píng)價(jià)因子作為優(yōu)化目標(biāo)，反演優(yōu)化得到不同的最優(yōu)肋片形狀，在橫肋的肋基面積給定的情況下，隨著優(yōu)化目標(biāo)中阻力特性權(quán)重的增大，最優(yōu)橫肋形狀的尖角逐漸平緩甚至消失；

2）擾流橫肋底部順列的布置方式優(yōu)于上下錯(cuò)列的布置方式，能夠獲得更高的綜合換熱性能；

3）在本文研究工況下，擾流橫肋數(shù)目為3時(shí)，能夠取得最佳的強(qiáng)化換熱效果，且沿著流動(dòng)方向3個(gè)肋片的最優(yōu)形狀逐漸從向內(nèi)微凹的三角形變成向外微凸的水滴狀。