基于對流傳熱反問題的矩形通道內擾流橫肋的形狀優化

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

（河北工業大學 能源與環境工程學院，天津 300401）

0 引言

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

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

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

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

1 模型及計算方法

1.1 計算模型

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

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

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

1.2 計算方法

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

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

連續性方程為

式中，ρ為流體的密度。

動量守恒方程為

式中，

湍流動能方程為

耗散率方程為

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

能量守恒方程為

式中，λ是導熱系數。

1.3 模型與網格驗證

對于計算區域進行網格獨立性驗證時，矩形通道采用結構化網格，擾流元左右20 mm附近采用非結構性網格加密，網格數大于109 942后，流道的Nu及f均不發生大的變化，其結果如圖2所示。因此，選擇此時的網格數進行數值模擬計算。

圖2 網格無關性考核Fig.2 Grid independence test

在進行優化之前，為了驗證模型的準確性，將擾流橫肋周圍的局部努塞爾數數值結果與Liou等[21]的實驗結果進行了比較，其結果如圖3所示。可以看出，模擬結果與實驗結果具有很好的吻合度，驗證了數值模型的可靠性。

圖3 模型驗證Fig.3 Model validation

2 反演優化方法

2.1 優化問題的目標函數

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

式中：a為系統傳熱能力與系統阻力的權重系數，其大小反映了優化過程中系統換熱能力及減阻效果的側重程度；Nu0和f0分別代表光滑通道的努塞爾數和摩擦系數，由Gnielinski和Filonenko公式定義：

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

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

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

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

由式（13）可得

其系數的確定過程為：

對式（14）兩邊同時求導得，

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

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

由邊界條件公式（11）可導出2個方程：

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

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

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

本文所用共軛梯度法構造一組搜索方向，求出目標函數的極值點。其求解的計算流程如下。

每個元素ri的更新方式為

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

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

1）采用數值模擬的方法，求解整個計算區域的速度和溫度分布，然后求得目標函數值；

2）根據正問題的結果和收斂準則ε，（ε是很小的正數），此時令k=1；

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

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

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

3 結果與討論

3.1 不同目標函數的對比

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

圖4和圖5分別為不同a值反演優化后所對應的優化結果。由圖可知，隨著a值逐漸增大，即隨著系統阻力權重的逐漸增大，優化后擾流橫肋結構的尖角逐漸平緩甚至消失。這是由于尖角的存在可以提高傳熱性能，但同時也會增大流動阻力。當a=1/3時，目標函數值達到1.869，但增量最小，僅為3.32%。可根據不同的工程要求來選取不同的a值。

圖4 不同a 值下的優化形狀Fig.4 Optimal shape with different a

圖5 不同a 值下優化結果Fig.5 Optimization results under different a

3.2 不同排布方式的對比

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

圖6 不同排布方式的結構（N=2）Fig.6 Schematic diagram of different arrangement（N=2）

圖7中4種擾流元布置方式，可以看出布置方式為底排的擾流元優化后所得的目標函數最大；順排布置的擾流元優化后所得的目標函數與底排布置的擾流元優化結果相比稍小；叉排布置的2 種擾流元優化后所得的目標函數遠小于底排布置的結構。

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

3.3 不同擾流橫肋個數的對比

3.3.1 邊界離散點數目的影響

以具有單個擾流元的通道為例，初始擾流元形狀為圖1 所示結構，此時目標函數為，其他邊界條件不變。分別對離散邊界點個數為m=3、7、11 的3 種情況進行反演。表2 給出了擬合前后取不同離散邊界點優化結果的速度與壓力分布。可以看出，擬合前后大致形狀基本類似。

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

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

圖8 不同點數下的優化結果Fig.8 Optimization results under different m

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

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

圖10 不同初始形狀下的優化結果Fig.10 Optimization results under different initial shapes

3.3.2 擾流橫肋個數的影響

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

圖11 不同擾流元個數下的優化結果Fig.11 Optimization results under different number of ribs

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

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

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

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

4 結論

本文以綜合換熱性能最佳為目標，采用共軛梯度法對矩形通道內擾流橫肋的最佳形狀進行了反演優化，對優化目標函數的選取、擾流橫肋布置方式及數目的影響進行了討論分析。得到如下結論：

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

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

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