胡和敏，杜小澤，楊立軍，楊勇平

（華北電力大學 電站設備狀態(tài)監(jiān)測與控制教育部重點實驗室，北京102206）

近年來，直接空冷機組以其顯著的節(jié)水優(yōu)勢，在我國北方“富煤缺水”地區(qū)取得了很大發(fā)展.然而，火電機組空冷島的熱流系統(tǒng)具有明顯的多尺度特性，其影響因素包括：10-3m 數(shù)量級的翅片結(jié)構(gòu)；101m數(shù)量級的Λ 型空冷凝汽器單元；空冷島的換熱性能還易受到電廠建筑，甚至103m 數(shù)量級的地形地貌的限制.不同特征尺度上空氣流場的湍流渦旋及其級聯(lián)破裂過程共同對空冷系統(tǒng)散熱性能產(chǎn)生影響.

在空冷島的設計和運行階段，計算流體力學（CFD）方法經(jīng)常被用于研究空冷島的對流換熱特性［1-4］.由于時間和計算資源有限，很難實現(xiàn)對跨尺度區(qū)域的同步模擬.且這些研究都是在某種假設的前提下，在某1個或2個尺度上進行的，而忽略其他尺度的影響.隨著對空冷島熱流系統(tǒng)研究的逐步深入，需要尋求一種能夠解決空冷島跨尺度熱流問題的模擬方法.

多尺度問題廣泛存在于自然科學的各個領域［5-6］，如環(huán)境科學、生物科學、材料科學、地殼運動以及流動傳熱等.近年來，多尺度科學引起了學者們的廣泛關注，其研究手段也得到了長足發(fā)展，其中一個典型的模擬方法是對不同尺度計算區(qū)域建立不同的物理模型，其數(shù)值解在大小尺度區(qū)域界面上耦合［7］.Nie等［8］采用這種分區(qū)建模界面耦合的方法來研究頂蓋驅(qū)動流動的多尺度問題，其中采用分子動力學（MD）方法來描述頂蓋附近的流動，用N-S方程模擬其他區(qū)域的流動.Wu等［9］將直接模擬技術和N-S方程相結(jié)合來模擬分析包括連續(xù)和稀薄形態(tài)的高速氣體流動.欒輝寶等［10］構(gòu)造了格子-Boltzmann方法（LBM）和有限容積法（FVM）的耦合模型，給出了重構(gòu)算子，并在典型流動傳熱問題中得到了驗證.Tang 等［11］用一種多層網(wǎng)格嵌入的多尺度方法來分析印刷板的流動和換熱特征，這種方法用粗網(wǎng)格來離散整體模型，而后逐漸向所關心的區(qū)域細化網(wǎng)格.

Rambo［12］基于正交分解方法（POD）提出了一種新的多尺度方法來預測數(shù)據(jù)中心冷卻系統(tǒng)的溫度場和速度場.其中，對風室和服務器構(gòu)件分別建立關于溫度場和速度場的POD 低維模型.在得到上游模塊POD 解的基礎上，利用兩個模塊之間界面上的通量守恒法（FMP）來得到下游模塊的POD 解.

為了更好地解決火電廠直接空冷系統(tǒng)的跨尺度熱流問題，筆者在前人工作的基礎上，提出一種新的多尺度湍流模擬方法.在小尺度區(qū)域建立POD 模型并求解，再將大小尺度界面上POD 解的參數(shù)信息提取出來，作為邊界條件賦值給大尺度區(qū)域.在此基礎上，用Fluent商業(yè)軟件來對大尺度區(qū)域進行CFD求解，預測了二維空冷單元的跨尺度流動特征.

1 物理數(shù)學模型

圖1給出了Λ 型空冷單元的示意圖.由圖1可知，Λ 型布置的空冷單元由汽輪機排汽管、翅片管束、軸流風機以及風筒組成.由軸流風機引入的空氣流經(jīng)翅片管束并與之發(fā)生對流換熱，將汽輪機排汽的凝結(jié)潛熱帶走.

圖1 Λ 型空冷單元的示意圖Fig.1 Schematic diagram of theΛ-frame air-cooled condenser

根據(jù)Oliveira等［13］的研究結(jié)果，為了得到與計算域大小無關的數(shù)值解，為POD 的低維模型提供精確的快照數(shù)據(jù)，計算域x方向和y方向的坐標取值范圍分別為（0，472.1）和（0，138），模型的細節(jié)尺寸見圖2.

圖2 空冷單元的物理模型及其計算域Fig.2 Physical model of the ACC unit and its computational domain

1.1 控制方程

假設空氣為不可壓縮流體，模型計算域的流動特性滿足如下控制方程

其中

式中：ρ為空氣密度；ui和uj分別表示x（或y）方向的速度分量；p為壓力；μ和μt 分別表示動力黏度和湍流黏度；δij為Kroneker符號.

數(shù)值模擬采用標準k-ε兩方程模型，湍動能k和湍流耗散率ε的控制方程如下

式中：Gk和Gb分別表示由平均速度梯度和浮升力所產(chǎn) 生 的 湍 動 能；C1，ε、C2，ε和C3，ε為 常 數(shù)，其 中C1，ε＝1.44，C2，ε＝1.92.

1.2 邊界條件

計算域迎風面（x＝0）采用速度入口邊界條件，大氣層內(nèi)空氣平均速度大小隨海拔高度按指數(shù)規(guī)律變化

式中：y和分別為任意海拔高度及其平均速度大小；y0和分別為參考高度及其平均速度大小，本文取y0＝10m；α為當?shù)氐牡孛娲植诙龋鶕?jù)日本荷載規(guī)范［14］，取α＝0.2.

翅片管束采用集總參數(shù)的散熱器模型，流經(jīng)管束的無量綱壓降為垂直于翅片表面速度分量的多項式.

軸流風機簡化為風機邊界條件.根據(jù)風機性能曲線，將流經(jīng)風機平面的壓力躍變Δp擬合為軸向速度分量w的4次多項式

計算域出口面（x＝472.1）采用出流邊界條件.

1.3 數(shù)值解法

根據(jù)文獻［15］中對模型進行的網(wǎng)格無關性檢驗結(jié)果，模型整體網(wǎng)格數(shù)為7.64×104，其中大尺度區(qū)域采用特征尺度為1.00大小的網(wǎng)格進行劃分.而空冷單元及其周圍的小尺度區(qū)域所包含的1.05×103個網(wǎng)格的特征尺度為0.35.整個模型計算域均采用四邊形網(wǎng)格進行劃分.

模型的控制方程采用有限差分法進行離散；除了壓力項采用標準差分格式外，其他方程均采用一階迎風差分格式；壓力與速度的耦合方式采用Simple算法；用Fluent商業(yè)軟件對模型進行求解.當所有變量的殘差均小于10-10時，認為計算結(jié)果收斂.

2 正交分解方法

POD 模型的計算域如圖1所示，由CFD 計算得到的關于變量Θ的快照被正交分解為一組POD基φ，即建立了POD 低維模型.在研究范圍內(nèi)的任意工況下Θc可以用這組POD 基線性表示，從而得到Θ的POD 解.

式中：源項為N個CFD 模擬結(jié)果的平均值.快照空間｛θ1，θ2，…，θi，…，θN｝被正交分解為N個POD基φi，其中θi＝Θi-.

對速度u、湍動能k和湍流耗散率ε3個變量建立POD 低維模型，用3次樣條插值得到POD 基的權系數(shù)ai，c.

更多關于POD 方法的理論基礎、特點及其實際操作過程（snapshot法）可參考文獻［15］.

3 多尺度湍流模擬方法

傳統(tǒng)CFD 方法中，基于計算資源和計算精度的綜合考慮，模型中流動結(jié)構(gòu)復雜的空冷單元附近區(qū)域需采用細網(wǎng)格進行劃分，而其他區(qū)域則可采用相對粗網(wǎng)格劃分.筆者認為圖1所示的細網(wǎng)格劃分的空冷單元區(qū)域為小尺度區(qū)域，除此之外的圖3所示的區(qū)域則為大尺度區(qū)域；其中交疊區(qū)共同屬于大小尺度區(qū)域.采用POD 方法快速準確地預測細網(wǎng)格劃分的小尺度區(qū)域的流動形態(tài)；而粗網(wǎng)格劃分的大尺度區(qū)域的流動形態(tài)則用N-S 方法描述.將從POD解中提取出的大小尺度界面上的變量信息賦值于大尺度區(qū)域，以實現(xiàn)大小尺度的界面耦合.多尺度湍流模擬方法的具體計算步驟如下：

（1）在風向為（1，0），取4m／s、5m／s、6m／s和7m／s的情況下，對圖2所示的物理模型進行CFD 計算，得到數(shù)值解.

（2）從CFD 數(shù)值解中取出圖1的小尺度區(qū)域中關于u、k和ε3個變量的快照數(shù)據(jù)，采用snapshot法對其建立POD 低維模型.

（3）采用3 次樣條插值法得到試驗工況為4.7m／s、5.7m／s和6.6m／s時的POD 解.

（4）從3個試驗工況的POD 解中提取出3個變量在CFD 模型邊界面上的參數(shù)分布.

（5）將3個參數(shù)分布以速度入口邊界條件的形式賦值于圖3所示的大尺度區(qū)域，而后利用Fluent商業(yè)軟件來求解包含交疊區(qū)域在內(nèi)的N-S方程.

假定傳統(tǒng)CFD 方法同時對整個區(qū)域（包含大小尺度在內(nèi)）的N-S方程進行迭代求解的計算結(jié)果為真值，多尺度湍流模擬方法所得的解中，某變量Θ相對于傳統(tǒng)CFD 解的誤差可表示為

式中：‖‖2表示2-范數(shù)；Θm和Θt分別為多尺度湍流模擬方法和傳統(tǒng)CFD 方法得到的解.

圖3 CFD模型及其計算域Fig.3 CFD model and its computational domain

4 結(jié)果與分析

圖4給出了3個試驗工況下POD 解的截斷誤差隨POD 基數(shù)目的變化規(guī)律.由于最后一個POD基的失真，將其舍去.從圖4可以看出，隨著POD 基數(shù)目的增大，誤差逐漸減小.在只有4個快照的情況下，所得POD 解中3個變量的誤差均小于4.35%，而速度的誤差可以控制在0.50%以內(nèi).對比3個試驗工況的誤差曲線可以看出，工況＝5.7m／s下3個變量所能達到的計算精度最高，其次是＝6.6 m／s，工況＝4.7m／s的精度最差.這是因為工況＝5.7m／s所對應的湍流狀態(tài)與源項所對應的湍流狀態(tài)最接近，因此可以得到更高精度的POD解.當試驗工況所對應的湍流狀態(tài)與源項所對應的湍流狀態(tài)偏差較大時，用POD 方法預測得到的結(jié)果精度會有所下降.這是POD 方法所固有的弊端，可以通過POD 補集（PODc）方法［10］得到改進.

表1給出了大尺度區(qū)域CFD 計算所得的誤差.由表1可知，與傳統(tǒng)CFD 計算結(jié)果相比，所有變量的誤差均小于3.55%，其中速度的誤差不大于0.40%，湍動能和湍流耗散率的誤差有所增大，壓力的誤差可以控制在1.20%左右.大尺度區(qū)域的CFD解與小尺度區(qū)域的POD 解的誤差隨u的變化規(guī)律不盡相同.這是因為小尺度區(qū)域的湍流狀態(tài)只是通過大小尺度之間界面上的變量信息對大尺度區(qū)域的CFD 計算起作用.而大尺度區(qū)域的CFD 迭代過程還會受其他邊界條件、計算模型和離散方式等的綜合影響.

圖4 POD解的截斷誤差隨POD基數(shù)目的變化Fig.4 Change of truncation error of POD solution with POD mode number

表1 大尺度區(qū)域的誤差Tab.1 Errors for the large-scale region

圖5為試驗工況＝5.7m／s時，多尺度湍流模擬方法與傳統(tǒng)CFD 方法計算所得的速度流線圖.圖5（a）為大尺度區(qū)域的流線圖，圖5（b）則是在圖5（a）基礎上增加了小尺度區(qū)域的POD 計算結(jié)果.與圖5（c）所示的傳統(tǒng)CFD 計算結(jié)果進行對比發(fā)現(xiàn)，2種方法的計算結(jié)果具有很高的吻合度.因此，筆者所提出的多尺度湍流模擬方法能夠很好地描述空氣流經(jīng)空冷單元所產(chǎn)生的尾渦.在大小尺度區(qū)域分開求解的情況下，揭示出空氣在風機、翅片管束、擋風墻以及環(huán)境風共同影響下的流動特性，能夠?qū)⒋笮〕叨鹊耐牧鳡顟B(tài)關聯(lián)起來，從而很好地揭示出空冷單元的多尺度流動特征.

圖5 試驗工況＝5.7m／s時的速度流線圖Fig.5 Velocity streamlines in the test case of＝5.7m／s

圖6給出了試驗工況＝5.7m／s時的壓力云圖.由圖6可知，與傳統(tǒng)CFD 方法計算結(jié)果相比，二者所得的壓力云圖十分吻合.另外，從圖4所示的小尺度區(qū)域各變量的誤差和表1列出的大尺度區(qū)域的誤差中，也很好地說明多尺度湍流模擬方法在處理二維空冷單元多尺度流動問題上具有足夠高的精度.

圖6 試驗工況u＝5.7m／s時的壓力云圖Fig.6 Pressure contours in the test case of ＝5.7m／s

由于POD 計算方法的準確性和高效性，幾乎可以忽略POD 模型求解所用的時間，因而只需考察多尺度方法中大尺度區(qū)域的計算時間.以工況＝5.7 m／s為例，用所提出的多尺度方法迭代300次需要102s，而用傳統(tǒng)CFD 方法則需要110s.由于模型中小尺度區(qū)域的網(wǎng)格數(shù)所占的比例較小，而用多尺度方法對大尺度區(qū)域進行CFD 計算時，在處理邊界信息上需要花費一定的時間.因而在迭代次數(shù)相等的情況下，多尺度方法的計算時間節(jié)省不多.若將本文所述的多尺度湍流模擬技術推廣至直接空冷系統(tǒng)中三維小尺度空冷單元與大尺度的電廠建筑乃至地形地貌的尺度耦合中去，效果會非常明顯.由于三維空冷單元結(jié)構(gòu)的復雜性，需要用極細網(wǎng)格進行劃分才能保證計算的準確性，其離散網(wǎng)格將會占模型網(wǎng)格總數(shù)的半數(shù)以上.在大尺度區(qū)域的CFD 計算過程中，細網(wǎng)格的小尺度區(qū)域處于反激活狀態(tài)，所用的計算資源和時間將會明顯節(jié)省.

表2為用2種方法對3個試驗工況進行求解，達到收斂時所需要的迭代次數(shù).從表2可以看出，由于給定了邊界信息，使多尺度方法比傳統(tǒng)CFD 方法更容易收斂，從而提高了計算效率.

表2 2種方法計算所需的迭代次數(shù)Tab.2 Iteration number required by both the methods

5 結(jié) 論

（1）考慮橫向環(huán)境風的作用，利用POD 方法對直接空冷凝汽器單元中空氣流動的變量分布進行了精確高效的預測.

（2）利用從空冷單元小尺度區(qū)域的POD 解中提取跨尺度界面上的變量信息，再賦值于大尺度區(qū)域的方法，可以準確地描述環(huán)境風流經(jīng)空冷單元的湍流特征，實現(xiàn)跨尺度界面耦合.

（3）分開求解小尺度區(qū)域的POD 模型和大尺度區(qū)域的CFD 模型，可大大節(jié)省計算資源.在保證一定計算精度的基礎上提高計算效率.

［1］周蘭欣，張情，李衛(wèi)華.直接空冷凝汽器噴霧增濕系統(tǒng)的結(jié)構(gòu)優(yōu)化［J］.動力工程學報，2011，31（2）：148-152. ZHOU Lanxin，ZHANG Qing，LI Weihua.Structural optimization for spray humidification system of a direct air-cooled condenser［J］.Journal of Chinese Society of Power Engineering，2011，31（2）：148-152.

［2］石維柱，安連鎖，張學鐳，等.直接空冷機組噴淋冷卻系統(tǒng)的數(shù)值模擬和性能分析［J］.動力工程學報，2010，30（7）：523-527. SHI Weizhu，AN Liansuo，ZHANG Xuelei，etal，Numerical simulation and performance analysis of spray cooling system in direct air cooling units［J］.Journal of Chinese Society of Power Engineering，2010，30（7）：523-527.

［3］周蘭欣，崔皓程，仲博學.主導風向?qū)χ苯涌绽淠鲹Q熱效率的影響［J］.動力工程學報，2010，30（2）：138-142. ZHOU Lanxin，CUI Haocheng，ZHONG Boxue.Influence of prevailing wind direction on heat transfer efficiency of direct air-cooled condensers［J］.Journal of Chinese Society of Power Engineering，2010，30（2）：138-142.

［4］劉達，楊建國，張兆營，等.水平擋板改善直接空冷凝汽器流場特性的數(shù)值模擬［J］.動力工程學報，2011，31（12）：949-954. LIU Da，YANG Jianguo，ZHANG Zhaoying，etal.Numerical simulation on flow field of direct air-cooled condenser improved by horizontal baffles［J］.Journal of Chinese Society of Power Engineering，2011，31（12）：949-954.

［5］何國威，夏蒙棼，柯孚久，等.多尺度耦合現(xiàn)象：挑戰(zhàn)和機遇［J］.自然科學進展，2004，14（2）：121-124. HE Guowei，XIA Mengfen，KE Fujiu，etal.Phenomena of coupling multiscale：challenges and opportunities［J］.Progress in Natural Science，2004，14（2）：121-124.

［6］柴立和.多尺度科學的研究進展［J］.化學進展，2005，17（2）：186-191. CHAI Lihe.Recent progress of multiscale science［J］.Progress in Chemistry，2005，17（2）：186-191.

［7］陶文銓.傳熱與流動問題的多尺度數(shù)值模擬：方法與應用［M］.北京：科學出版社，2008.

［8］NIE Xiaobo，CHEN Shiyi，ROBBINS M O.Hybrid continuum-atomistic simulation of singular corner flow ［J］.Physics of Fluids，2004，16（10）：3579-3591.

［9］WU J S，LIAN Y Y，CHENG G，etal.Development and verification of a coupled DSMC-NS scheme using unstructured mesh［J］.Journal of Computational Physics，2006，219（2）：579-607.

［10］欒輝寶，徐輝，陳黎，等.有限容積法與格子Boltzmann方法耦合模擬傳熱流動問題［J］.科學通報，2010，55（32）：3128-3140. LUAN Huibao，XU Hui，CHEN Li，etal.Coupling between FVM and LBM for heat transfer and fluid flow problems［J］.Chinese Science Bulletin，2010，55（32）：3128-3140.

［11］TANG L，JOSHI Y K.A multi-grid based multi-scale thermal analysis approach for combined approach for combined mixed convection，conduction，and radiation due to discrete heating［J］.Journal of Heat Transfer，2005，127（1）：18-26.

［12］RAMBO J D.Reduced-order modeling of multiscale turbulent convection：application of data center thermal management［D］.Atlanta，USA：Georgia Institute of Technology，2006.

［13］OLIVEIRA P J，YOUNIS B A.On the prediction of turbulent flows around full-scale buildings［J］.Journal of Wind Engineering and Industrial Aerodynamics，2000，86（2／3）：203-220.

［14］JUN Kanda.AIJ-RLB-2004AIJ recommendations for loads on buildings［S］.Japan：Architectural Institute of Japan，2004.

［15］DU Xiaoze，HU Hemin，SHEN Yinqi，etal.，Reduced order analysis of flow and heat transfer for aircooled condenser of power generating unit［J］.Applied Thermal Engineering，2013，51（1／2）：383-392.