直接空冷系統流動特性的多尺度模擬研究

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

（華北電力大學 電站設備狀態監測與控制教育部重點實驗室，北京102206）

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

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

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

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

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

1 物理數學模型

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

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

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

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

1.1 控制方程

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

其中

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

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

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

1.2 邊界條件

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

式中：y和分別為任意海拔高度及其平均速度大小；y0和分別為參考高度及其平均速度大小，本文取y0＝10m；α為當地的地面粗糙度，根據日本荷載規范［14］，取α＝0.2.

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

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

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

1.3 數值解法

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

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

2 正交分解方法

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

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

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

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

3 多尺度湍流模擬方法

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

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

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

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

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

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

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

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

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

4 結果與分析

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

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

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

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

圖5為試驗工況＝5.7m／s時，多尺度湍流模擬方法與傳統CFD 方法計算所得的速度流線圖.圖5（a）為大尺度區域的流線圖，圖5（b）則是在圖5（a）基礎上增加了小尺度區域的POD 計算結果.與圖5（c）所示的傳統CFD 計算結果進行對比發現，2種方法的計算結果具有很高的吻合度.因此，筆者所提出的多尺度湍流模擬方法能夠很好地描述空氣流經空冷單元所產生的尾渦.在大小尺度區域分開求解的情況下，揭示出空氣在風機、翅片管束、擋風墻以及環境風共同影響下的流動特性，能夠將大小尺度的湍流狀態關聯起來，從而很好地揭示出空冷單元的多尺度流動特征.

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

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

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

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

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

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

5 結 論

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

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

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

［1］周蘭欣，張情，李衛華.直接空冷凝汽器噴霧增濕系統的結構優化［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］石維柱，安連鎖，張學鐳，等.直接空冷機組噴淋冷卻系統的數值模擬和性能分析［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］周蘭欣，崔皓程，仲博學.主導風向對直接空冷凝汽器換熱效率的影響［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］劉達，楊建國，張兆營，等.水平擋板改善直接空冷凝汽器流場特性的數值模擬［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］何國威，夏蒙棼，柯孚久，等.多尺度耦合現象：挑戰和機遇［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］陶文銓.傳熱與流動問題的多尺度數值模擬：方法與應用［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.