湍流普朗特數識別的隨機抽樣算法

朱 嵩，劉國華，程偉平，黃躍飛

湍流普朗特數識別的隨機抽樣算法

朱 嵩1，劉國華2，程偉平2，黃躍飛3

（1.廣東省電力設計研究院水務部，廣州 510663；2.浙江大學建筑工程學院，杭州 310058；3.清華大學水沙科學與水利水電工程國家重點實驗室，北京 100084）

在溫排水等涉及熱交換的環境水力學研究中，湍流普朗特（Prandtl，簡稱Pr）數是控制溫度的主要參數。對于一個特定的問題，傳統湍流Pr數的確定方法主要采用經驗法或試錯法，因而具有一定盲目性和低效性。為了提高湍流Pr數確定的可靠性，采用馬爾科夫鏈蒙特卡羅（Markov Chain Monte Carlo，簡稱MCMC）隨機抽樣的方法（Metropolis-Hastings算法）來對湍流Pr數進行識別，其中湍流場計算采用了穩態標準k-ε模型，溫度場計算采用非穩態熱傳導方程。算例計算結果表明，MCMC方法對湍流Pr數的識別具有良好的適用性和較高的識別精度。

湍流Prandtl數；參數識別；湍流傳熱；Metropolis-Hastings算法；MCMC隨機抽樣

1 概 述

在電廠的水工工藝設計過程中，需要解決大量的湍流傳熱問題。如濕冷塔內的水氣兩相換熱、空冷平臺熱交換等。此類問題的主要特征是換熱與湍流流動的耦合，屬于湍流擴散的研究范疇。湍流擴散遵循著質量守恒、動量守恒和能量守恒方程，可以采用廣義對流擴散方程描述。要準確獲得湍流熱擴散中的溫度場，關鍵是要確定湍流熱擴散系數。研究已經表明，湍流渦黏性系數與湍流熱擴散系數之比為一常數，稱為Prandtl（Pr）數。不同流動下的湍流Pr數是不同的，一般根據流動的類型經驗加以確定，或通過手工調節的方法和觀測值對比從而找到當前流動的湍流Pr數。前一種方法，由于湍流Pr數的范圍仍然較大，選擇起來仍然具有主觀性；后一種雖然通過試驗結果加以率定，但手工調節的方式費時費力。

目前，優化算法、人工智能等方法已廣泛應用于相關的環境水力學反問題研究中［1］，如遺傳算法［2，3］、模擬退火算法［4］以及貝葉斯方法［5－9］等。由于貝葉斯方法的先驗、后驗概念以及馬爾科夫鏈蒙特卡羅（Markov Chain Monte Carlo，MCMC）的強大后驗分布抽樣能力，使得該方法近年來獲得了飛速的發展和廣泛的應用。相比較優化算法而言，MCMC方法能獲得整個參數的后驗分布信息。相比較人工神經網絡等黑箱模型而言，貝葉斯MCMC方法具有明確的數學理論指導。

由于傳熱傳質問題的類比性，在湍流傳熱傳質研究過程中，Schmidt數和Pr數在一些場合下并不區分，實指同一含義［9］。Yoshihide Tominaga和Ted Stathopoulos調查總結了射流、濁流、邊界層羽流彌散和建筑物附近的流動中的湍流Schmidt數的經驗取值，指出湍流Schmidt數對預測結果具有較大的影響，應該根據不同類型的流動具體地加以研究［10］。Yanhu Guo等人采用遺傳算法（Genetic Algorithm，GA）估計了橫流中射流變系數Schmidt數，取得了較好的結果［11］。

采用貝葉斯方法以及MCMC抽樣解決傳熱學方面的文獻可參見文獻［12－14］，但目前尚未見貝葉斯MCMC方法在湍流傳熱反問題研究上的應用。本文采用Metropolis-Hastings算法對湍流Pr數進行了識別研究。湍流傳熱物理場耦合場的求解采用有限單元法計算，其中湍流場采用穩態標準k-ε模型求解，溫度場采用非穩態熱傳導方程求解。采用的二維湍流熱擴散的算例計算表明貝葉斯MCMC方法能夠較好地識別湍流傳熱過程中的Pr數。

2 湍流傳熱數學模型

本文采用非穩態熱傳導方程描述熱擴散現象，其中湍流流速場采用k-ε模型求解，控制方程如下：

式中：ρ為密度；Cp為熱容；T為溫度；u為流速矢量；kt為湍流熱擴散系數；p為壓強；η為流體動力黏度；ηt為湍流運動黏度；c1，c2，cμ，σk和 σε為標準k-ε模型的5個參數；Pr為湍流Prandtl數。

3 貝葉斯推理與MCMC抽樣方法

3.1 貝葉斯定理和反問題

貝葉斯推理的基礎是貝葉斯定理，它可以表述如下

式中：θ，Y分別為模型參數和觀測數據；ρM（θ），L（Y｜θ），σM（θ｜Y）分別為模型的先驗概率密度函數，似然函數和后驗概率密度函數；ρM（θ）表示在未收集到數據之前關于參數的認識，它主要來源于以往的觀測資料、經驗和主觀判斷等；L（Y｜θ）表示模型參數擬合實測數據的程度，值越大表明擬合程度越好，反之則差；σM（θ｜Y）即為反問題的解，它表示在獲得觀測數據之后模型參數的分布規律。

3.2 抽樣方法

一旦獲得了參數的后驗分布，就可以獲得參數的具體特征，如單個參數的邊緣分布或均值、方差等。但對于很簡單的情況，后驗概率分布都不能以解析解的形式給出，所以必須采用數值積分方法。馬爾可夫鏈蒙特卡羅（MCMC）就是一種對復雜問題在高維空間上的數值積分方案。它的主要思路就是創造一個隨機游動（馬爾可夫過程）使得后驗概率密度函數σM（θ｜Y）作為它的靜態分布，只要這個計算過程愈長，那么采樣結果就愈服從后驗概率分布。

構造馬爾可夫鏈的方法很多，主要包括Gibbs采樣，Metropolis（Metropolis-Hasting）算法等。本文采用Metropolis算法對模型參數后驗分布進行采樣，算法如下：

（1）給定參數的初始點 θ（0）；

（2）在當前模型參數θ（t）的狀態下給定一個隨機的無偏擾動，從而生成一個新的狀態θ′，計算出似然函數 L（θ（t））和 L（θ′）；

（3）生成［0，1］間均勻分布的隨機數U；

（4）如果 L（θ′）≤L（θ（t））或則接受θ′，否則拒絕，轉到（2）直到預定的次數。

4 算例計算

4.1 算例一

假設在0.9 m×0.25 m的矩形域如圖1內存在一水體的湍流傳熱過程，流場和溫度場的計算采用穩態標準k-ε模型和非穩態熱傳導方程。采用有限元法離散計算，計算域共剖分了2 800個四邊形網格，采用二次形函數。邊界條件設置如表1。進口流速為1 m／s、進口湍流強度為5%，進口溫度為343.15 K，初始溫度為293.15 K，水的熱容為 4 182J×kg×K－1。湍流模型參數取值如表2。

圖1 計算域Fig.1 Computational domain

表1 邊界條件設置Table 1 Setting of the boundary condition

表2 標準k-ε湍流模型參數取值Table 2 Values of standard k-εturbulence model

為了考查MCMC算法對該問題的適用性，預設湍流Pr數的真值為0.8，通過后驗分析來驗證貝葉斯MCMC參數識別的精度。在距離進口0.5 m處的對稱軸上設置一測點 O。取 O點上 0，0.1，0.2，0.3，0.4，0.5，0.6，0.7，0.8，0.9以及 1 s時 11個時刻上的溫度值作為觀測數據，如圖2所示。

考慮到湍流Pr數幾乎對湍流場沒有影響，因而識別過程中只計算一次湍流場，而反復調用溫度場求解模型。設計的湍流Pr數識別的流程如下：

（1）根據幾何生成有限元網格；

（2）根據流場的參數及邊界條件計算得到湍流場，為標量場的計算提供背景；

圖2 測點O上的溫度觀測值Fig.2 Measured temperature at point O

（3）在先驗范圍內產生湍流Pr數的初值，計算似然函數，采用Metropolis-Hastings算法產生湍流Pr數的后驗樣本；

（4）對后驗樣本進行統計來識別未知湍流Pr數。

MCMC參數設置如下：Pr數的先驗范圍為［0.2，1.0］，先驗分布為均勻分布。似然函數構造基于溫度噪聲為白噪聲N［0，σ2］，其中 σ假設為0.1。隨機游走的步長為先驗范圍的5%，迭代次數為100次。

湍流Pr數的迭代曲線如圖3所示。可見，MCMC迭代在經歷過約20步初始化階段后，進入了統計收斂區域，主要表現為圍繞真值0.8上下小幅震蕩。圖4為湍流Pr數的后驗統計結果。從圖中可以看出，湍流Pr數的后驗概率在真值0.8附近出現的概率最高，符合預設的判斷。表3為后驗統計結果，從表3中可以看出，后驗均值和后驗中位數均具有較高的識別精度。

圖3 湍流Pr數迭代曲線Fig.3 Iteration curve of turbulent Pr number

圖4 湍流Pr數的后驗分布Fig.4 Posterior distribution of turbulent Pr number

表3 湍流Pr數后驗統計結果Table 3 Posterior results of turbulent Pr number

4.2 算例二

在一個0.889 m×0.304 8 m的二維計算域中有一橫流中射流（Jet in Crossflow），射流小孔長度為0.025 4 m，射 流 流 速 為 2.3 m／s，射 流 溫 度 為293.15 K；橫 流 流 速 為 1 m／s，橫 流 溫 度 為283.15 K，流動介質為水。假設湍流Pr數為1。有限元計算考慮溫度對水的密度及熱容的影響，流場方程和溫度場方程進行耦合求解。計算得到的溫度場如圖5所示。

圖5 射流溫度場分布（單位：K）Fig.5 Temperature field of jet flow（Unit in K）

測 點 布 置 在 （0.5，0.05），（0.5，0.1），（0.5，0.15），（0.5，0.2）和（0.5，0.25）5個位置，通過有限元計算得到測點上的溫度值分別為291.148 9 K，290.506 1 K，283.860 3 K，283.150 4 K和283.150 0 K。

MCMC參數設置如下：湍流Pr數的先驗范圍為［0.2，1.0］，先驗分布為均勻分布。似然函數為拉普拉斯函數，標準差為0.1。隨機游走的步長為先驗范圍的5%，迭代次數為1 000次。

湍流Pr數的迭代曲線如圖6所示。可見，MCMC迭代在經歷過約100步初始化階段后，進入了統計收斂區域，主要表現為圍繞真值1.0上下小幅震蕩。圖7為湍流Pr數的后驗統計結果。從圖中可以看出，湍流Pr數的后驗概率在真值1.0附近出現的概率最高，符合預設值。

圖6 湍流Pr數迭代曲線Fig.6 Iteration curve of turbulent Pr number

圖7 湍流后驗Pr數的后驗分布Fig.7 Posterior distribution of turbulent Pr number

5 結論與展望

湍流普朗特數是控制湍流擴散中溫度分布的關鍵參數，本文采用一種馬爾科夫鏈蒙特卡羅方法Metropolis-Hastings算法對湍流穩態傳熱和非穩態傳熱兩種典型模型進行了參數識別。其中，湍流求解采用了標準k-ε模型，離散方法采用了通用性較強的有限元法。計算結果表明，基于Metropolis-Hastings算法的湍流普朗特數識別能夠達到較高的精度。

計算中的數據采用湍流時均值，考慮脈動測量值以及高級湍流模型的關鍵參數識別是需要進一步研究的課題。

［1］ 劉曉東，姚 琪，薛紅琴，等.環境水力學反問題研究進展［J］.水科學研究進展，2009，20（6）：885－893.（LIU Xiao-dong，YAO Qi，XUE Hong-qin，et al.Advance in Inverse Problems of Environmental Hydraulics［J］.Advances in Water Science，2009，20（6）：855－893.（in Chinese））

［2］ NG A W M，PERERA B JC.Selection of Genetic Algorithm Operators for River Water Quality Model Calibration［J］.Engineering Applications of Artificial Intelligence，2003，16（5－6）：529－541.

［3］ 朱 嵩，毛根海，劉國華.基于FVM-HGA的河流水質模型多參數識別［J］.水力發電學報，2007，26（6）：91－95.（ZHU Song，MAOGen-hai，LIU Guo-hua.Parameters Identification of River Water Quality Model Based on Finite Volume Method-Hybrid Genetic Algorithm［J］.Journal of Hydroelectric Engineering，2007，36（6）：91－95.（in Chinese））

［4］ 王 薇，曾光明，何 理.用模擬退火算法估計水質模型參數［J］.水利學報，2004，（6）：61－67.（WANG Wei，ZENG Guang-ming，HE Li.Estimation of Water Quality Model Parameters with Simulated Annealing Algorithm［J］.Journal of Hydraulic Engineering，2004，（6）：61－67.（in Chinese））

［5］ 朱 嵩，毛根海，程偉平，等.利用貝葉斯推理估計二維含源對流擴散方程參數［J］.四川大學學報（工程科學版），2008，40（2）：38－43.（ZHU Song，MAO Genhai，CHENG Wei-ping，et al.Application of Bayesian Inference to Estimate the Parameters in 2D Convection-Diffusion Equation with Source［J］.Journal of Sichuan University（Engineering Science Edition），2008，40（2）：38－43.（in Chinese））

［6］ 朱 嵩，毛根海，劉國華，等.改進的MCMC方法及其應用［J］.水利學報，2009，40（8）：1019－1023.（ZHU Song，MAO Gen-hai，LIU Guo-hua，et al.Improved MCMC Method and Its Application［J］.Journal of Hydraulic Engineering，2009，40（8）：1019－1023.（in Chinese））

［7］ 朱 嵩，劉國華，王立忠.水動力-水質耦合模型污染源識別的貝葉斯方法［J］.四川大學學報（工程科學版），2009，41（5）：30－35.（ZHU Song，LIU Guo-hua，WANG Li-zhong.A Bayesian Approach for Identification of the Pollution Source in Water Quality Model Coupled with Hydrodynamics［J］.Journal of Sichuan University（Engineering Science Edition），2009，41（5）：30－35.（in Chinese））

［8］ 朱 嵩.基于貝葉斯推理的環境水力學反問題研究［D］.杭州：浙江大學，2008.（ZHU Song.Research on Inverse Problems of Environmental Hydraulics by Bayesian Inference［D］.Hangzhou：Zhejiang University，2008.（in Chinese））

［9］ 朱 嵩.湍流標量輸運過程中關鍵參數的識別研究［R］.杭州：浙江大學，2010.（ZHU Song.Study on the Identification for Key Parameters of the Scalar Transport Process in the Turbulence［R］.Hangzhou：Zhejiang University，2008.（in Chinese））

［10］TOMINAGA Y，STATHOPOULOST.Turbulent Schmidt Numbers for CFD Analysis with Various Types of Flowfield［J］.Atmospheric Environment，2007，41（37）：8091－8099.

［11］GUO Yan-hu，HE Guang-bin，ANDREW T H.Application of Genetic Algorithms to the Development of a Variable Schmidt Number Model for Jet-in-crossflows［J］.International Journal of Numerical Methods for Heat＆Fluid Flow，2007，11（8）：744－761.

［12］WANG Jing-bo，ZABARASN.Using Bayesian Statistics in the Estimation of Heat Source in Radiation［J］.International Journal of Heat and Mass Transfer，2005，48：15－29.

［13］WANG Jing-bo，ZABARASN.A Markov Random Field Model of Contamination Source Identification in Porous Media Flow［J］.International Journal of Heat and Mass Transfer，2006，49：939－950.

［14］LING L，YAMAMOTO M，HON Y C，et al.Identification of Source Locations in Two-Dimensional Heat Equations［J］.Inverse Problems，2006，22：1289－1305.

A Random Sampling Algorithm for Identification of Turbulent Prandtl Number

ZHU Song1，LIU Guo-hua2，CHENG Wei-ping2，HUANG Yue-fei3
（1.Water Engineering Department，Guangdong Electric Power Design Institute，Guangzhou 510663，China；2.College of Civil Engineering and Architecture，Zhejiang University，Hangzhou 310058，China；3.State Key Laboratory of Hydroscience and Engineering，Tsinghua University，Beijing 100084，China）

Turbulent Prandtl number（Pr）is a key parameter for controlling the temperature distribution in the research of thermal discharge water and other environmental hydraulics involving heat transfer.For a given problem，turbulent Pr number generally comes from previous experiences or trial-and-error method，which is blindfold and inefficient.To increase the reliability of turbulent Pr number for a given problem，Metropolis-Hastings algorithm，a Markov Chain Monte Carlo（MCMC）random sampling method was employed in this paper to identify turbulent Pr number.In the numerical simulation，steady standard k-εmodel was used for turbulence flow field computation，while unsteady heat transfer equation was adopted for computing the temperature field.The computation results manifested that MCMCmethod is suitable and can offer precise results for the identification of turbulent Pr number.

turbulent Prandtl number；parameter identification；turbulent heat transfer；Metropolis-Hastings algorithm；MCMC random sampling

O242.28，O357.5＋3

A

1001－5485（2011）06－0016－04

2010-07-01

國家重點基礎研究計劃資助項目（2005CB724202）；國家自然科學基金資助項目（50879075）

朱 嵩（1981-），男，安徽潁上人，博士后，工程師，從事計算流體力學和電廠水工工藝研究，（電話）020-32117746（電子信箱）zhusong＠gedi.com.cn。

（編輯：王 慰）