湍流普朗特?cái)?shù)識(shí)別的隨機(jī)抽樣算法

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

湍流普朗特?cái)?shù)識(shí)別的隨機(jī)抽樣算法

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

（1.廣東省電力設(shè)計(jì)研究院水務(wù)部，廣州 510663；2.浙江大學(xué)建筑工程學(xué)院，杭州 310058；3.清華大學(xué)水沙科學(xué)與水利水電工程國家重點(diǎn)實(shí)驗(yàn)室，北京 100084）

在溫排水等涉及熱交換的環(huán)境水力學(xué)研究中，湍流普朗特（Prandtl，簡稱Pr）數(shù)是控制溫度的主要參數(shù)。對(duì)于一個(gè)特定的問題，傳統(tǒng)湍流Pr數(shù)的確定方法主要采用經(jīng)驗(yàn)法或試錯(cuò)法，因而具有一定盲目性和低效性。為了提高湍流Pr數(shù)確定的可靠性，采用馬爾科夫鏈蒙特卡羅（Markov Chain Monte Carlo，簡稱MCMC）隨機(jī)抽樣的方法（Metropolis-Hastings算法）來對(duì)湍流Pr數(shù)進(jìn)行識(shí)別，其中湍流場計(jì)算采用了穩(wěn)態(tài)標(biāo)準(zhǔn)k-ε模型，溫度場計(jì)算采用非穩(wěn)態(tài)熱傳導(dǎo)方程。算例計(jì)算結(jié)果表明，MCMC方法對(duì)湍流Pr數(shù)的識(shí)別具有良好的適用性和較高的識(shí)別精度。

湍流Prandtl數(shù)；參數(shù)識(shí)別；湍流傳熱；Metropolis-Hastings算法；MCMC隨機(jī)抽樣

1 概 述

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

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

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

采用貝葉斯方法以及MCMC抽樣解決傳熱學(xué)方面的文獻(xiàn)可參見文獻(xiàn)［12－14］，但目前尚未見貝葉斯MCMC方法在湍流傳熱反問題研究上的應(yīng)用。本文采用Metropolis-Hastings算法對(duì)湍流Pr數(shù)進(jìn)行了識(shí)別研究。湍流傳熱物理場耦合場的求解采用有限單元法計(jì)算，其中湍流場采用穩(wěn)態(tài)標(biāo)準(zhǔn)k-ε模型求解，溫度場采用非穩(wěn)態(tài)熱傳導(dǎo)方程求解。采用的二維湍流熱擴(kuò)散的算例計(jì)算表明貝葉斯MCMC方法能夠較好地識(shí)別湍流傳熱過程中的Pr數(shù)。

2 湍流傳熱數(shù)學(xué)模型

本文采用非穩(wěn)態(tài)熱傳導(dǎo)方程描述熱擴(kuò)散現(xiàn)象，其中湍流流速場采用k-ε模型求解，控制方程如下：

式中：ρ為密度；Cp為熱容；T為溫度；u為流速矢量；kt為湍流熱擴(kuò)散系數(shù)；p為壓強(qiáng)；η為流體動(dòng)力黏度；ηt為湍流運(yùn)動(dòng)黏度；c1，c2，cμ，σk和 σε為標(biāo)準(zhǔn)k-ε模型的5個(gè)參數(shù)；Pr為湍流Prandtl數(shù)。

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

3.1 貝葉斯定理和反問題

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

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

3.2 抽樣方法

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

構(gòu)造馬爾可夫鏈的方法很多，主要包括Gibbs采樣，Metropolis（Metropolis-Hasting）算法等。本文采用Metropolis算法對(duì)模型參數(shù)后驗(yàn)分布進(jìn)行采樣，算法如下：

（1）給定參數(shù)的初始點(diǎn) θ（0）；

（2）在當(dāng)前模型參數(shù)θ（t）的狀態(tài)下給定一個(gè)隨機(jī)的無偏擾動(dòng)，從而生成一個(gè)新的狀態(tài)θ′，計(jì)算出似然函數(shù) L（θ（t））和 L（θ′）；

（3）生成［0，1］間均勻分布的隨機(jī)數(shù)U；

（4）如果 L（θ′）≤L（θ（t））或則接受θ′，否則拒絕，轉(zhuǎn)到（2）直到預(yù)定的次數(shù)。

4 算例計(jì)算

4.1 算例一

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

圖1 計(jì)算域Fig.1 Computational domain

表1 邊界條件設(shè)置Table 1 Setting of the boundary condition

表2 標(biāo)準(zhǔn)k-ε湍流模型參數(shù)取值Table 2 Values of standard k-εturbulence model

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

考慮到湍流Pr數(shù)幾乎對(duì)湍流場沒有影響，因而識(shí)別過程中只計(jì)算一次湍流場，而反復(fù)調(diào)用溫度場求解模型。設(shè)計(jì)的湍流Pr數(shù)識(shí)別的流程如下：

（1）根據(jù)幾何生成有限元網(wǎng)格；

（2）根據(jù)流場的參數(shù)及邊界條件計(jì)算得到湍流場，為標(biāo)量場的計(jì)算提供背景；

圖2 測點(diǎn)O上的溫度觀測值Fig.2 Measured temperature at point O

（3）在先驗(yàn)范圍內(nèi)產(chǎn)生湍流Pr數(shù)的初值，計(jì)算似然函數(shù)，采用Metropolis-Hastings算法產(chǎn)生湍流Pr數(shù)的后驗(yàn)樣本；

（4）對(duì)后驗(yàn)樣本進(jìn)行統(tǒng)計(jì)來識(shí)別未知湍流Pr數(shù)。

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

湍流Pr數(shù)的迭代曲線如圖3所示?？梢?，MCMC迭代在經(jīng)歷過約20步初始化階段后，進(jìn)入了統(tǒng)計(jì)收斂區(qū)域，主要表現(xiàn)為圍繞真值0.8上下小幅震蕩。圖4為湍流Pr數(shù)的后驗(yàn)統(tǒng)計(jì)結(jié)果。從圖中可以看出，湍流Pr數(shù)的后驗(yàn)概率在真值0.8附近出現(xiàn)的概率最高，符合預(yù)設(shè)的判斷。表3為后驗(yàn)統(tǒng)計(jì)結(jié)果，從表3中可以看出，后驗(yàn)均值和后驗(yàn)中位數(shù)均具有較高的識(shí)別精度。

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

圖4 湍流Pr數(shù)的后驗(yàn)分布Fig.4 Posterior distribution of turbulent Pr number

表3 湍流Pr數(shù)后驗(yàn)統(tǒng)計(jì)結(jié)果Table 3 Posterior results of turbulent Pr number

4.2 算例二

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

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

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

MCMC參數(shù)設(shè)置如下：湍流Pr數(shù)的先驗(yàn)范圍為［0.2，1.0］，先驗(yàn)分布為均勻分布。似然函數(shù)為拉普拉斯函數(shù)，標(biāo)準(zhǔn)差為0.1。隨機(jī)游走的步長為先驗(yàn)范圍的5%，迭代次數(shù)為1 000次。

湍流Pr數(shù)的迭代曲線如圖6所示?？梢?，MCMC迭代在經(jīng)歷過約100步初始化階段后，進(jìn)入了統(tǒng)計(jì)收斂區(qū)域，主要表現(xiàn)為圍繞真值1.0上下小幅震蕩。圖7為湍流Pr數(shù)的后驗(yàn)統(tǒng)計(jì)結(jié)果。從圖中可以看出，湍流Pr數(shù)的后驗(yàn)概率在真值1.0附近出現(xiàn)的概率最高，符合預(yù)設(shè)值。

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

圖7 湍流后驗(yàn)Pr數(shù)的后驗(yàn)分布Fig.7 Posterior distribution of turbulent Pr number

5 結(jié)論與展望

湍流普朗特?cái)?shù)是控制湍流擴(kuò)散中溫度分布的關(guān)鍵參數(shù)，本文采用一種馬爾科夫鏈蒙特卡羅方法Metropolis-Hastings算法對(duì)湍流穩(wěn)態(tài)傳熱和非穩(wěn)態(tài)傳熱兩種典型模型進(jìn)行了參數(shù)識(shí)別。其中，湍流求解采用了標(biāo)準(zhǔn)k-ε模型，離散方法采用了通用性較強(qiáng)的有限元法。計(jì)算結(jié)果表明，基于Metropolis-Hastings算法的湍流普朗特?cái)?shù)識(shí)別能夠達(dá)到較高的精度。

計(jì)算中的數(shù)據(jù)采用湍流時(shí)均值，考慮脈動(dòng)測量值以及高級(jí)湍流模型的關(guān)鍵參數(shù)識(shí)別是需要進(jìn)一步研究的課題。

［1］ 劉曉東，姚 琪，薛紅琴，等.環(huán)境水力學(xué)反問題研究進(jìn)展［J］.水科學(xué)研究進(jìn)展，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的河流水質(zhì)模型多參數(shù)識(shí)別［J］.水力發(fā)電學(xué)報(bào)，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ì)水質(zhì)模型參數(shù)［J］.水利學(xué)報(bào)，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ì)二維含源對(duì)流擴(kuò)散方程參數(shù)［J］.四川大學(xué)學(xué)報(bào)（工程科學(xué)版），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］ 朱 嵩，毛根海，劉國華，等.改進(jìn)的MCMC方法及其應(yīng)用［J］.水利學(xué)報(bào)，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］ 朱 嵩，劉國華，王立忠.水動(dòng)力-水質(zhì)耦合模型污染源識(shí)別的貝葉斯方法［J］.四川大學(xué)學(xué)報(bào)（工程科學(xué)版），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］ 朱 嵩.基于貝葉斯推理的環(huán)境水力學(xué)反問題研究［D］.杭州：浙江大學(xué)，2008.（ZHU Song.Research on Inverse Problems of Environmental Hydraulics by Bayesian Inference［D］.Hangzhou：Zhejiang University，2008.（in Chinese））

［9］ 朱 嵩.湍流標(biāo)量輸運(yùn)過程中關(guān)鍵參數(shù)的識(shí)別研究［R］.杭州：浙江大學(xué)，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

國家重點(diǎn)基礎(chǔ)研究計(jì)劃資助項(xiàng)目（2005CB724202）；國家自然科學(xué)基金資助項(xiàng)目（50879075）

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

（編輯：王 慰）