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

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

0 引言

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

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

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

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

1 模型及計算方法

1.1 計算模型

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

圖1 計算模型示意圖Fig.1 Schematic diagram of model

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

1.2 計算方法

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

流動與換熱過程的控制方程如下。

連續(xù)性方程為

式中，ρ為流體的密度。

動量守恒方程為

式中，

湍流動能方程為

耗散率方程為

式中：S為應(yīng)變張量的平均速率之模，其定義式為S=(2SijSij)1/2，Sij=1/2(?uj/?xi+?uj/?xi)；αk和αε分別為k和ε的有效普朗特數(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)格驗證

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

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

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

圖3 模型驗證Fig.3 Model validation

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

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

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

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

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

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

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

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

由式（13）可得

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

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

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

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

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

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

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

2.3 簡化共軛梯度法求解反問題的計算流程

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

每個元素ri的更新方式為

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

反問題求解的計算過程如下：

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

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

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

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

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

3 結(jié)果與討論

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖13 不同擾流橫肋個數(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)，采用共軛梯度法對矩形通道內(nèi)擾流橫肋的最佳形狀進(jìn)行了反演優(yōu)化，對優(yōu)化目標(biāo)函數(shù)的選取、擾流橫肋布置方式及數(shù)目的影響進(jìn)行了討論分析。得到如下結(jié)論：

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

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

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