基于ISIGHT圓柱繞流的換熱強化與減阻優(yōu)化設(shè)計

賴海清，張東輝，陳 寧，石 珊

(江蘇科技大學 能源與動力工程學院，江蘇 鎮(zhèn)江212003)

0 引 言

在能源開發(fā)利用中，熱量傳遞現(xiàn)象十分普遍，強化換熱技術(shù)得到發(fā)展。對于對流換熱強化問題，需尋求換熱性能和流動阻力的最佳權(quán)衡，它與具體結(jié)構(gòu)配置緊密相關(guān)。研究表明:利用圓柱強化對流換熱能帶來更好的強化換熱效果，比如縱向渦發(fā)生器、換熱器、渦輪葉片冷卻結(jié)構(gòu)中的擾流柱、微肋管等［1］。在過去20年中，研究者們對雙圓柱繞流的氣動特性進行了大量的實驗研究和數(shù)值計算。Stephen Wornom［2］研究了表明光滑單圓柱體的阻力系數(shù)Cd 與雷諾數(shù)Re 有很大關(guān)系。當Re 很小時，無漩渦發(fā)生，Cd 與Re 成反比;Re 增至104時，壓強阻力是主要的，故繞流阻力幾乎與Re 無關(guān)。Zukauskas［3］發(fā)現(xiàn)在叉排管束中，流體流動非常類似在交替收縮與擴張的彎曲通道中流動，產(chǎn)生強烈的擾動，可抵消陰影效應(yīng)的影響，使叉排情況換熱相應(yīng)變得強烈，但阻力也會同時增加。Zhang 等［4］利用二維層流模型研究了雷諾數(shù)Re=200 時，擾動體對下游圓柱渦的形成和脫落的影響。臺灣學者許立杰［5］研究了在低雷諾數(shù)下，兩圓柱間隔大小不同和兩圓柱比例不同對流動形態(tài)和對主圓柱阻力的影響。南京航空航天大學的張喜東［6］研究了雷諾數(shù)Re=1 000時鈍體直徑比，間隙率，對雙鈍體系統(tǒng)傳熱性能和系統(tǒng)阻力的影響規(guī)律，結(jié)果表明，小直徑擾鈍體對系統(tǒng)的對流換熱更為有利，且系統(tǒng)的對流換熱隨間隙率的增大而遞增。上述研究主要針對圓柱繞流的尾跡結(jié)構(gòu)和流動阻力開展，關(guān)于傳熱的研究僅限于圓柱共線排布及所有圓柱鈍體均發(fā)熱時的對流換熱。且實驗研究只在可壓流體(如空氣)方面開展，并沒有研究控制圓柱與主圓柱之間夾角與間隙的耦合影響。

在探討優(yōu)化結(jié)構(gòu)過程中，由于參數(shù)，需要考慮特征值繁多，即使采用CAM/CAE 技術(shù)，仍然需要大量的手工操作，耗時耗力。而隨著計算機技術(shù)及數(shù)值計算技術(shù)的發(fā)展，采用參數(shù)優(yōu)化設(shè)計成為了一種必須。ISIGHT 是目前國際上優(yōu)秀的綜合性計算機輔助工程軟件之一。ISIGHT 軟件將大量需要人工完成的工作由軟件實現(xiàn)自動化處理，從而替代工程設(shè)計者進行重復(fù)性、易出錯的數(shù)字和設(shè)計處理工作，因此ISIGHT被稱為“軟件機器人”［7］。ISIGHT 強大的集成功能，可以在設(shè)計過程中充分利用各學科先進的分析工具，并將其集成起來以實現(xiàn)設(shè)計流程的自動化。

因此，本文采用ISIGHT 軟件集成GAMBIT，F(xiàn)luent 及自身所帶的數(shù)據(jù)處理模塊，利用有限體積法系統(tǒng)的研究當Re=200 時，雙圓柱繞流系統(tǒng)的對流換熱和阻力性能。通過數(shù)值模擬獲得了主圓柱的時均努謝爾數(shù)和阻力系數(shù)變化規(guī)律，揭示了在合理配置圓柱間夾角和間隙率的情況下，即可得到在減小主柱體阻力同時，還能增強換熱性能的優(yōu)化結(jié)構(gòu)體。

1 控制方程及參數(shù)說明

對層流粘性不可壓縮流體，在直角坐標系下，其運動規(guī)律可用N-S 方程來描述。

連續(xù)性方程:

動量方程:

能量方程:

當橫流流經(jīng)圓柱，其尾流場的流場現(xiàn)象主要與流體特性、圓柱的直徑、流體的流速有關(guān)，即主要與無因次參數(shù)雷諾數(shù)有關(guān)系。

式中:ρ 為流體密度;U 為自由來流速度;D 為圓柱直徑;μ 為流體動力粘度。

柱體表面與周圍流體的對流換熱為:

柱體壁面上的局部、時間平均、面平均及時間和面平均努塞爾數(shù)可分別表示為:

常將流體作用于圓柱表面的作用力分為與流動方向平行的阻力及與流動方向垂直的擾動升力，以無因次阻力系數(shù)Cd 及升力系數(shù)Cl 來表示這2個方向的作用力。

式中:FD為阻力;FL為升力。

系統(tǒng)阻力減少率計算方法如下式:

2 計算模型

物理計算模型如圖1所示，主圓柱直徑為D，取控制圓柱直徑d=0.5 D 為研究對象，流域邊界高為20D，長度方向劃分為:入流場10 D，尾流場20 D，以減小邊界帶來的誤差。為研究間隙與夾角對流場及換熱的影響，取控制小圓柱與主圓柱間隙L/D和圓柱間夾角degree 為優(yōu)化變量，夾角取入流方向與兩圓柱中心連線形成的角度，方向為自上游沿順時針向下游。

圖1 物理計算模型Fig.1 Physical calculationmodel

3 計算方法及算例驗證

為保證本研究的可行有效，首先探討單圓柱繞流的特性，本文網(wǎng)格劃分均采用非結(jié)構(gòu)化網(wǎng)格劃分法，在近場、流動比較復(fù)雜、流動變化大的區(qū)域布置較為密集的網(wǎng)格，在圓柱表面采用結(jié)構(gòu)化邊界層網(wǎng)格進行局部加密。最后整體網(wǎng)格網(wǎng)格數(shù)為88651，網(wǎng)格檢查equisize skew 為0.12。如圖1所示，入口選為速度入口，遠場邊界壁面采用動網(wǎng)格，移動速度等于入口來流速度，這是為了最大程度減小上下壁面對流動的干擾。壓力速度耦合采取PISO 方法，動量、能量方程用二階迎風格式離散，連續(xù)性方程和動量方程收斂殘差標準均為10-5，以提高計算精度，從速度入口為初始點計算。取Re=200 計算，監(jiān)視并輸出圓柱表面的Cd和Cl。當流動穩(wěn)定(尾流震蕩已達穩(wěn)定狀態(tài))后，得到圓柱的時均阻力系數(shù)Cd，升力系數(shù)Cl。采用NU 數(shù)來表征對流換熱的強弱，依照式(7)求其平均數(shù)，得到NU 為12.41。各目標量與權(quán)威文獻對比結(jié)果如表1所示。

表1 本文計算結(jié)果與文獻對比Tab.1 The results compared with the literature in this paper

其中邱吉爾關(guān)聯(lián)式

茹卡烏斯卡斯關(guān)聯(lián)式:

式中:Prf為參考溫度時的Pr;Prw為圓柱表面溫度時流體的Pr。

將Fluent 計算結(jié)果導(dǎo)入后處理軟件Tecplot 中，由渦量計算公式vor=ddx(v)－ ddy(u)，得到圓柱渦量圖(見圖2)，可發(fā)現(xiàn)與蘇炯彰［9］的模擬結(jié)果相符。由上述可知，本研究的網(wǎng)格質(zhì)量、求解條件等設(shè)置滿足計算精度要求。

圖2 渦量比較圖(蘇炯彰)Fig.2 Tthe vorticity comparison chart (SU Jiong-zhang)

4 基于ISIGHT的雙圓柱計算分析

4.1 過程集成

過程集成就是實現(xiàn)設(shè)計過程的自動化，實現(xiàn)自動由變量參數(shù)完成物理模型的建模，然后網(wǎng)格化計算流域，由Fluent 計算輸出計算結(jié)果，最后經(jīng)數(shù)據(jù)處理得到所需的目標變量。整個自動化主要流程如圖3所示。

Mesh.jou 模塊記錄著圓柱系統(tǒng)建模及網(wǎng)格畫法命令流文件;Gambit.bat 實現(xiàn)由Gambit 讀取mesh.jou 命令文件劃分網(wǎng)格并導(dǎo)出網(wǎng)格文件;Solve.jou 模塊包括Fluent 一系列前處理，如邊界條件、初始條件、材料屬性、求解方法等操作的命令流;Fluent.bat 為批處理文件，其作用是讀取Slove.jou 進行求解，并在同一目錄按流動時間生成包含NU，Cd，dat 等結(jié)果文件;Calculator 為數(shù)據(jù)后處理模塊，計算出NU和Cd的時均值;DOE (實驗設(shè)計)選用全因子法進行全面的探索設(shè)計空間，得到設(shè)計空間內(nèi)各個設(shè)計變量對優(yōu)化問題的影響情況，得到一個最優(yōu)點，這樣即可避免出現(xiàn)局部最優(yōu)解的情況。整個流程由ISIGHT 集成，避免了繁瑣且極易出錯的重復(fù)人工操作。數(shù)學描述:約束:1D≤L/D≤3D，0°≤degree≤180°，因每個完整方案計算下來大約需要4 h，為節(jié)省計算成本，且并具有代表性，設(shè)置L/D 步長為0.5D，degree 增長因子為30°，且變量均為整型，取NU 權(quán)重大于Cd，求使主圓柱的時均努謝爾數(shù)NU 最大化，阻力系數(shù)Cd 最小化的最優(yōu)解。

圖3 過程集成示意圖Fig.3 Process integration schemes

4.2 設(shè)計結(jié)果分析

經(jīng)計算，最終得到最優(yōu)設(shè)計方案為degree=30°，L/D=2，其換熱強度NU=13.2，阻力系數(shù)Cd=1.21。相對單圓柱NU 提高了6.37%，同時阻力系數(shù)減少了7.63%。圖4 為全方案下各變量得出的目標值擬合圖，橫坐標為圓柱間夾角大小。

從圖4 可得:一定范圍內(nèi)，夾角不變時，間隙率越大，主圓柱NU 越大，即系統(tǒng)換熱能力越強;間隙率不變時，隨著夾角的逐漸增大，主圓柱的換熱強度和阻力均表現(xiàn)為先增大，再平緩，最后減小的趨勢。當雙圓柱夾角為0°或180°時，系統(tǒng)對主圓柱的影響是換熱性能較低，但減阻性能得到了很大的提升，究其原因，取L/D=2 時的溫度流線圖分析，如圖5 當小直徑擾動體與來流方向共線(夾角為0°)時，經(jīng)過初步熱交換后，由于雙圓柱之間存在回流，導(dǎo)致擾動體之后，主圓柱之前這部分區(qū)域流體的溫度高于來流靜溫，使主圓柱前半周迎風面與流體溫差減小，換熱系數(shù)降低，從而導(dǎo)致?lián)Q熱不佳，而在其他夾角范圍內(nèi)(30°～120°)，則不存在此類問題，故換熱較好，同理可分析夾角為180°時換熱不良的原因;減阻方面，由于控制圓柱與主圓柱組成的系統(tǒng)為空穴模式，張［10］等研究表明，空穴模式的減阻性能優(yōu)良，當處于這種流動模式時，上游圓柱對下游主圓柱起保護作用，即主圓柱處于流動死區(qū)，故其阻力系數(shù)變得很小。

圖4 各組合參數(shù)下NU，Cd 擬合圖Fig.4 Tthe combination of NU，Cd fitting figure

圖5 L/D=2 溫度流線圖Fig.5 L/D=2，temperature flow diagram

為進一步分析在夾角為30°～150°之間的雙圓柱系統(tǒng)與單圓柱相比換熱增強的原因，圖6 給出了兩系統(tǒng)的溫度等值線圖。可以看出，由于控制圓柱的加入，其后產(chǎn)生震蕩渦街，將“新鮮”的流體卷入主圓柱尾流中，由于這種交互作用，直接導(dǎo)致主圓柱溫度場的尾跡區(qū)變寬，從而提高了主圓柱背風面的換熱強度。

圖6 溫度等值線圖Fig.6 Temperature contour map

圓柱繞流阻力可分為壓差阻力和摩擦阻力，在Re≤1 時，阻力以摩擦阻力為主，但隨著雷諾數(shù)的增大，壓差阻力占總阻力份額逐漸增大，實際流體繞圓柱等凸壁流動時，在背風面容易發(fā)生邊界層分離開成旋渦，后部壓強不能恢復(fù)到與前部相對稱的程度，因此圓柱前后產(chǎn)生壓強差，這就是壓差阻力。在尾流分離區(qū)內(nèi)，壓強大致是均勻分布，對于具有熱對流交換的圓柱繞流，表現(xiàn)為:發(fā)生邊界層分離之前，換熱強度沿著圓柱表面由最大值逐漸減小，在分離點因邊界層最厚使換熱性能達到最低值，隨后由于旋渦的產(chǎn)生將低溫的尾流卷入近壁面，導(dǎo)致?lián)Q熱再增強。并在后駐點達到極大值。圖7 為3 種情況下主圓柱表面當?shù)豊U 數(shù)分布圖。從圖中可看出，單圓柱繞流時，分離點在120°左右，此結(jié)果與理論研究文獻及實驗結(jié)果［5］能夠很好的吻合。在受控制圓柱影響后，主圓柱表面分離點發(fā)生改變，在夾角為90°時，分離點有所提前;當夾角變?yōu)?0°，分離點滯后，大概在125°左右。正是由于分離滯后(流線型物體減阻成因)，使圓柱表面壓差阻力減小，從而導(dǎo)致阻力系數(shù)降低，這也解釋了為什么夾角為30°時的圓柱阻力小于90°。

圖7 圓柱表面的NU 分布圖Fig.7 NU distribution on the surface of the cylinder

5 結(jié) 語

1)利用ISIGHT 強大的集成能力及優(yōu)化設(shè)計能力，通過集成Gambit，F(xiàn)luent 軟件并結(jié)合自身數(shù)據(jù)處理模塊，針對雙圓柱繞流模型進行了參數(shù)化設(shè)計優(yōu)化模擬，使整個設(shè)計流程自動化，極大地提高了設(shè)計效率。

2)對直徑比為0.5，利用全因子法得出圓柱間隙和夾角與主圓柱的對流換熱性能及阻力的關(guān)系，發(fā)現(xiàn)在主圓柱周圍不同方位布置控制圓柱，對系統(tǒng)的換熱性能和阻力影響不同，與單圓柱系統(tǒng)相比，在合適的組合參數(shù)(L/D=2，degree=30°)下，雙圓柱系統(tǒng)的利用可以使主圓柱減阻7.63%的同時提高換熱性能6.73%。

3)本文只對Re=200，d/D=0.5的雙圓柱進行探討，今后研究方向應(yīng)該結(jié)合更多的參數(shù)，在不同雷諾數(shù)，多尺度比情況下，結(jié)合參數(shù)化優(yōu)化設(shè)計，討論最優(yōu)參數(shù)組合。

［1］陶文銓，多尺度換熱數(shù)值模擬［M］.北京:科學出版社，2009.

［2］WORNOM S，OUVRARD H，et al.Variational multi-scale largeeddy simulations of the flow past a circular cylinder:Reynolds number effects［J］.Computers ＆Fluids，2011，47:44-50.

［3］ZUKAUSKAS A A.MA Wen-chang，et al.The convective heat transfer of heat exchanger［M］.Beijing:Science Press，1986.

［4］ZHANG P F，WANG J J，HUANG L X.Numerical simulation of flow around cylinder with an upstream rod in tandem at low reynolds numbers［J］.Applied Ocean Research，2006，28:183-192.

［5］許立杰，田文政.以頻譜元素法模擬圓柱流場及其減阻方法研究探討［D］.國立云林科技大學，2009.

［6］ZHANG Xi-dong，HUANG Hu-lin.Cylindrical heat transfer enhancement disturbance in the flow around the body and the drag reduction［J］.Journal of chemical industry，http://www.cnki.net/kcms/detail/11.1946.TQ.20130619.1516.002.html.

［7］Engineering Ltd.ISIGHT——Multidisciplinary design optimization platform［J］.Journal of CAD/CAM and Manu-Facturing Information，2003(12):56-58.

［8］CHUCHILL S W，BENSTEInm.A correlating eauation for forced convection from gases and liquids to a circular cylinder in cross flow［J］.ASME J Heat Tansfer，1997，99(1):300-306.

［9］SU Jiong-akira，ALE broadband element method in the development and application of incompressible［D］.Mechanical Engineering，National Taiwan University，Master′s Thesis，2006.

［10］張攀峰，王晉軍，魯素芬，等.采用上游擾動體被動控制減阻技術(shù)的實驗研究［J］.空氣動力學學報，2006，24(4):114-514.