基于工程反問題的圓管湍流流體熱物性參數(shù)反演

梅 寧，焦 思

（中國海洋大學(xué)工程學(xué)院，山東青島266100）

基于工程反問題的圓管湍流流體熱物性參數(shù)反演

梅 寧，焦 思

（中國海洋大學(xué)工程學(xué)院，山東青島266100）

給出一種利用測定圓管對流換熱溫度場求解湍流流體熱物性參數(shù)的方法。根據(jù)邊界層理論，對管內(nèi)湍流強(qiáng)制對流換熱控制方程組進(jìn)行簡化，并通過考慮湍流運(yùn)輸中的密度脈動(dòng)，對傳統(tǒng)的混合長度湍流模型進(jìn)行改進(jìn)，基于該計(jì)算模型求解圓管湍流模型的溫度場。采用Levenberg-Marquardt方法通過測量溫度場反演了湍流對流換熱控制方程中的黏度系數(shù)。結(jié)果表明，利用反問題方法通過測定溫度場反演湍流流動(dòng)流體的熱物性參數(shù)具有可行性，精確度較高。

Levenberg-Marquardt方法；混合長度理論；湍流對流換熱；黏度系數(shù)

作為CFCs的替代物之一，由幾種單質(zhì)制冷劑按比例混合而成的非共沸混合制冷劑具有較為廣闊的市場空間［1］。非共沸制冷劑的溫度滑移特性使其在循環(huán)過程中組分逐漸偏移從而對制冷機(jī)組的性能造成影響，嚴(yán)重時(shí)甚至損壞機(jī)組。目前唯一的解決辦法是將系統(tǒng)內(nèi)的變質(zhì)制冷劑全部抽出，再重新注入標(biāo)準(zhǔn)配比的制冷劑，顯然投資巨大并且浪費(fèi)資源。因此廢棄制冷劑的回收再生問題顯得十分迫切。目前市場上有較多的制冷劑回收再生裝置，但是都不具有組分檢測的功能。本文嘗試尋找一種熱力學(xué)方法，實(shí)現(xiàn)混合制冷劑的組分檢測。用熱力學(xué)方法實(shí)現(xiàn)混合制冷劑的組分檢測主要包括混合流體的熱物性反演和根據(jù)熱物性與組分的關(guān)系求解組分2個(gè)方面內(nèi)容，本文主要針對流體熱物性參數(shù)反演展開研究。

湍流流動(dòng)是工程技術(shù)領(lǐng)域與自然界中常見的流動(dòng)現(xiàn)象，流體作湍流流動(dòng)時(shí)的對流換熱也是工程傳熱過程中最常見的一種熱交換方式。關(guān)于湍流運(yùn)動(dòng)與換熱的數(shù)值計(jì)算，是目前計(jì)算流體力學(xué)與計(jì)算傳熱學(xué)中困難最多因而研究最活躍的領(lǐng)域之一。一方面為了簡化計(jì)算，湍流對流換熱計(jì)算模型通常會(huì)被定為常流動(dòng)模型；另一方面在旺盛湍流區(qū)內(nèi)，由分子擴(kuò)散引起的熱物性變化相對于有湍流脈動(dòng)引起的物性變化較小可以忽略不計(jì)，因此很少有人關(guān)注湍流模型中的熱物性參數(shù)。本文嘗試采用Levenberg-Marquardt方法通過采集的溫度場對換熱流體的熱物性參數(shù)進(jìn)行反演。

基于可觀測量如溫度等運(yùn)用反問題方法推算混合物熱參數(shù)已經(jīng)有很多成熟的反問題方法：基于正則方法的Levenberg-Marquardt方法［2-3］、共扼梯度法［4-5］；對反問題經(jīng)典方法正則化方法進(jìn)行改進(jìn)的迭代正則化方法［6］；從系統(tǒng)學(xué)發(fā)展形成的遺傳算法、神經(jīng)網(wǎng)絡(luò)算法［7］、MCMC方法［8］以及基于廣義概率的貝葉斯方法［9-11］等。傳統(tǒng)的正則方法以其相對較小的計(jì)算量得到廣泛的應(yīng)用；系統(tǒng)學(xué)方法的很強(qiáng)的適應(yīng)性也使其得到廣泛關(guān)注；對于重復(fù)性反問題具有先天優(yōu)勢的貝葉斯方法，在重復(fù)假設(shè)算例中得到了進(jìn)一步發(fā)展。

本文根據(jù)一維穩(wěn)態(tài)熱傳導(dǎo)模型設(shè)計(jì)出可以實(shí)時(shí)測量溫度數(shù)據(jù)的試驗(yàn)裝置。運(yùn)用Levenberg-Marquardt方法以及熱物性參數(shù)的經(jīng)驗(yàn)公式［12］，結(jié)合試驗(yàn)結(jié)果，求得實(shí)驗(yàn)流體的熱物性參數(shù)。反演結(jié)果表明，利用反問題求解湍流流動(dòng)流體的熱物性參數(shù)具有可行性，為下一步分析混合制冷劑組分反演做了鋪墊。

1 正問題的數(shù)學(xué)模型

光滑圓管中的充分發(fā)展湍流的平均流動(dòng)有如下特點(diǎn)：（1）平均速度場是平行流。（2）除壓強(qiáng)外，一切平均量只和圓管中徑向坐標(biāo)有關(guān)。因此本文取充分發(fā)展區(qū)的圓管湍流對流換熱模型為研究對象，并做了進(jìn)一步簡化。

（1）忽略徑向流動(dòng)，將管內(nèi)湍流對流換熱作為一維導(dǎo)熱型問題。

（2）管外壁通過均勻高溫來流加熱，銅質(zhì)管壁的熱阻可忽略。

（3）通過加熱段之后的熱電偶組測定徑向節(jié)點(diǎn)溫度，即測定管內(nèi)有限節(jié)點(diǎn)換熱穩(wěn)定時(shí)無內(nèi)熱源溫度場。

根據(jù)邊界層理論，通過數(shù)量級分析，有管內(nèi)湍流對流換熱控制方程組：

圖1 圓管湍流充分發(fā)展對流換熱示意圖Fig.1 Diagram of heat flow model

試驗(yàn)段內(nèi)，壓力變化較小，認(rèn)為熱物性參數(shù)僅受溫度影響，上式中隨溫度發(fā)生變化的熱物性參數(shù)有密度ρ、定壓比熱容cp（本文中將密度與定壓比熱容的乘積定義為體積熱容，用C表示）、黏度系數(shù)μ、導(dǎo)熱系數(shù)λ。實(shí)驗(yàn)溫度范圍內(nèi)，黏度系數(shù)隨溫度的變化為Δμ／μ=18.8%，而ΔC／C=0.22%，Δλ／λ=2%，因此選擇黏度系數(shù)作為熱物性參數(shù)的反演目標(biāo)［13］。為簡化方程將隨溫度變化較小的體積熱容和導(dǎo)熱系數(shù)設(shè)定為常系數(shù)。

定義下列無量綱量。

管內(nèi)湍流對流換熱控制方程組變?yōu)椋?/p>

從（1）、（2）式可以看出，湍流黏性系數(shù)μt及湍流導(dǎo)熱系數(shù)λt均為未知量，需要引入湍流模型來確定。

2 混合長度理論的變物性修正

在對流換熱過程中，流體的物性隨著溫度的改變而改變，而溫度場的變化將對流體的流動(dòng)產(chǎn)生影響，此時(shí)速度場和溫度場是相互耦合的，在實(shí)際研究過程中，若同時(shí)對流體的每個(gè)物性變化產(chǎn)生的影響進(jìn)行分析，將使問題變的非常復(fù)雜，通常采取的辦法是重點(diǎn)分析其中的某個(gè)物性變化，該物性的變化會(huì)引起溫度場的變化，而溫度場的變化又會(huì)引起整個(gè)物性場的變化。本文采取的方法是考慮密度的變化，進(jìn)而得到適用于變物性計(jì)算的混合長度模型［14］。

對于直管內(nèi)的二維邊界層型湍流流動(dòng)，考慮密度脈動(dòng)，湍流應(yīng)力的表達(dá)式將變成：

對于邊界層型流動(dòng)，υ的值近于零，忽略三階脈動(dòng)值乘積項(xiàng)，則（6）式變?yōu)椋?/p>

根據(jù)混合長度理論有：

將（9）代入（8）整理得到：

變物性流動(dòng)時(shí)湍流黏性系數(shù)表示如下：

定義Fm為密度脈動(dòng)時(shí)的變物性修正因子，μt，o表示原始混合長度模型計(jì)算的湍流黏性系數(shù)，

在本文數(shù)值計(jì)算中，原始混合長度湍流計(jì)算模型用Nikurades公式計(jì)算，并按vanDriest提出的阻尼函數(shù)進(jìn)行修正［15］：

3 液體黏度的計(jì)算

液體黏度隨溫度升高而減小，在低于正常沸點(diǎn)下，液體黏度的對數(shù)值與溫度的倒數(shù)近似呈線性關(guān)系：

當(dāng)溫度高于沸點(diǎn)時(shí)，lnμ與1／T關(guān)系圖中出現(xiàn)了彎曲，上式計(jì)算誤差增大。因而有以下2個(gè)公式：

對于不同的流體上述3個(gè)公式的擬合精度不同，因此對于不同的流體，選取的經(jīng)驗(yàn)公式不同，本文選取水作為實(shí)驗(yàn)介質(zhì)，公式（17）的擬合精度最高。記P=［A B C D］，即為待反演參數(shù)。上述公式中的A，B，C，D為系數(shù)值，可以查表得到。

應(yīng)用反問題通過采集的溫度場求解流體的熱物性參數(shù)即為對上述系數(shù)值A(chǔ)，B，C，D的反演。

4 Levenberg-Marquardt方法

4.1反問題模型

Levenberg-Marquardt方法是特殊化正則方法Gauss-Newton迭代方法，它可以看作對非線性問題作先線性化后正則化的過程［16］。用來求解下述關(guān)于參數(shù)P的函數(shù)的最小值：該方法通過輸入的參數(shù)P求解正問題得到Tnum并通過敏感系數(shù)矩陣對參數(shù)P進(jìn)行修正。對于給定初始猜測值P0，在第k步，有如下計(jì)算式：

上式中的J即為敏感系數(shù)矩陣，敏感系數(shù)矩陣可采用多種方法計(jì)算，對于有確定表達(dá)形式的，可以采用微分鏈?zhǔn)椒▌t求解。本文的敏感系數(shù)矩陣為：

μLM為Levenberg-Marquardt參數(shù)。μLM減少了迭代過程中的震蕩和反問題的不適定性，從而使得迭代過程快速的向最佳的參數(shù)P收斂。μLM的引入使得參數(shù)估計(jì)存在偏差，但隨著迭代向最優(yōu)參數(shù)值收斂這個(gè)誤差逐漸減小并可忽略。

4.2 L-M方法迭代實(shí)施

采用L-M方法求解圓管內(nèi)流體熱物性參數(shù)的過程：首先測量圓管內(nèi)流體的溫度值；然后假設(shè)一個(gè)未知的熱物性參數(shù)，在該假設(shè)下通過求解正問題得到管內(nèi)流體的溫度場，根據(jù)溫度的測量值和計(jì)算值，計(jì)算目標(biāo)函數(shù)，判斷是否符合精度要求，如果不滿足，則對假設(shè)的未知的熱物性參數(shù)進(jìn)行修正，修正后繼續(xù)通過正問題求出測量點(diǎn)的溫度值，計(jì)算目標(biāo)函數(shù)，直至得到滿足條件的收斂值。

假設(shè)給定初始位置的測量溫度數(shù)據(jù)為Tcal，未知的熱物性參數(shù)的假設(shè)值為P0，μLM的初始值為，則Levenberg-Marquardt算法的具體過程為：

圖2 Levenberg-Marquardt算法的迭代程序Fig.2 Flow chart of Levenberg-Marquardt method

5 試驗(yàn)驗(yàn)證及結(jié)果分析

5.1試驗(yàn)臺(tái)的描述

如圖3所示實(shí)驗(yàn)裝置主要有3部分組成：實(shí)驗(yàn)流體的循環(huán)系統(tǒng)，加熱流體循環(huán)系統(tǒng)和信號采集分析系統(tǒng)。裝置選取內(nèi)徑為8 mm壁厚1 mm的銅管作為實(shí)驗(yàn)流體的循環(huán)管道。流體速度為1.5 m／s，初始溫度為291.5 K，恒溫水浴加熱溫度為333 K。均勻的流體流經(jīng)充分長的入口段進(jìn)入長20 cm的加熱段。管內(nèi)溫度由布置在加熱段末端沿徑向分布的4個(gè)K型熱電偶采集。信號通過A／D轉(zhuǎn)換器輸入MATLAB2009a編寫的程序中進(jìn)行處理。

圖3 實(shí)驗(yàn)裝置示意圖Fig.3 Convective heat-transfer experimental apparatus

5.2實(shí)驗(yàn)結(jié)果及分析

本文以水做為實(shí)驗(yàn)流體進(jìn)行了實(shí)驗(yàn)。由文獻(xiàn)［12］可知，在0～100℃的溫度范圍內(nèi)，式（17）中的系數(shù)分別為A=-24.71，B=4 209，C=4.527×10-2，D= -3.376×10-5。表1給出了5組不同的參數(shù)初始猜測值。

表1 參數(shù)的初始猜測值Table 1 The initial value of parameters

圖4、5中的曲線表示水的黏度系數(shù)在計(jì)算溫度區(qū)域內(nèi)（273.15～373.15K）隨溫度變化的關(guān)系曲線、實(shí)驗(yàn)擬合曲線和誤差曲線。從圖中曲線可以看出，反問題方法對于黏度系數(shù)能進(jìn)行有效的反演，當(dāng)輸入不同的初始猜測值時(shí)，誤差可以控制在11%內(nèi)。從誤差曲線可以看出該方法對于待測流體溫度較高接近沸點(diǎn)的區(qū)域內(nèi)有較高的精度。

圖4 不同初值條件下的擬合曲線Fig.4 Fitting parameter curves in different initial values

圖5 不同初值條件下的相對誤差曲線Fig.5 relative error curves in different initial values

表2給出了不同初始猜測值對應(yīng)的擬合參數(shù)值及絕對誤差。將第一組初值作為標(biāo)準(zhǔn)，后4組數(shù)據(jù)只變動(dòng)1個(gè)系數(shù)的初值進(jìn)行反演得到結(jié)果見表2。對比幾組結(jié)果，可以看出當(dāng)初始值偏離真實(shí)值越大時(shí)，反演得到結(jié)果的誤差越大。A、B的反演誤差普遍比C、D的誤差小說明了反問題方法對于系數(shù)A、B的敏感性較強(qiáng)而對C、D的敏感性相對較弱。由1、2組結(jié)果可以看出當(dāng)A的初值偏離實(shí)際值較大時(shí)，A參數(shù)反演的誤差變化并不大，但是D的誤差增倍；對比1、3組結(jié)果，當(dāng)B的初值變大時(shí)，反演得到的其他參數(shù)的誤差全部增大：A值誤差由1.21%上升到3.12%，B的誤差上升了1.47%，C的誤差由1.13%上升到2.10%，D的誤差上升了1.69%；對比1、4組結(jié)果，當(dāng)C的初值偏離實(shí)際值較小時(shí)，對A的反演有一定正影響，對于B、D影響不大，但是對于自身反演誤差并沒有優(yōu)勢；對比1、5組數(shù)據(jù)可以發(fā)現(xiàn)，D初值的變化對于A、B的影響是相反的，但是對于自身和C的影響與C的情況是相同的。

由上可以看出反演初值對于反演結(jié)果有一定的影響，對于A、B 2個(gè)系數(shù)，初值越接近真實(shí)值，反演的精度越高，初始猜測值與真實(shí)值的偏差較大時(shí)，反演結(jié)果的誤差也較大，所以在選擇初值時(shí)，應(yīng)盡量選擇合理的初始猜測值；但是對于C、D 2個(gè)系數(shù)，初始值偏離真實(shí)值的大小對于自身的影響不成正比關(guān)系，初始猜測值與真實(shí)值偏差較小時(shí)，對A、B的反演精度有正影響。綜上，合理的初始猜測值對于降低誤差提高精度有一定影響，并且4個(gè)系數(shù)之間的誤差關(guān)系不是獨(dú)立的而是相關(guān)的，有一定的疊加作用。

表2 不同初值條件下的擬合參數(shù)值及絕對誤差Table 2 Fitting parameters and absolute error in different initial values

6 結(jié)語

本文從湍流流體的熱物性參數(shù)的不易確定性出發(fā)，給出了一種可以反演得到湍流流動(dòng)流體熱物性參數(shù)的一般方法。

根據(jù)參數(shù)的初始猜測值，利用改進(jìn)的混合長度模型求解正問題得到了圓管內(nèi)流體的溫度場，并應(yīng)用LM方法中的敏感系數(shù)矩陣對系數(shù)的猜測值進(jìn)行修正直至得到滿足誤差的溫度場和系數(shù)猜測值，根據(jù)熱物性參數(shù)的半經(jīng)驗(yàn)公式得到待反演熱物性參數(shù)值。

本文的研究結(jié)果表明，利用L-M方法及局部的溫度測量信息，能夠?qū)A管充分發(fā)展區(qū)對流換熱流體的熱物性參數(shù)（黏質(zhì)系數(shù)）進(jìn)行有效的反演；反問題方法對于4個(gè)系數(shù)的敏感性不同，對于A、B的敏感性較好但是對于C、D的敏感性較差；4個(gè)系數(shù)并不是獨(dú)立存在的，而是相關(guān)的，并且彼此之間的反演誤差會(huì)對其他系數(shù)造成影響；合理的初始猜測值有利于降低誤差提高精度。

［1］ 俞炳豐，彭伯彥.CFCs制冷劑的回收與再利用［M］.北京：機(jī)械工業(yè)出版社，2007：47-50.

［2］ Mejias M，Orlande H R B，M?zisik.A Comparison of different parameter estimate on techniques for the identification of thermal conductivity components of orthotropic solids［C］.Washington：3rd Int Conf on Inverse Problems in Eng，1999.

［3］ SAWAF B，OZISIK N，Determining the constant thermal conductivities of orthotropic materials by inverse analysis［J］.International Communications in Heat and Mass Transfer，1995，22（2）：201-211.

［4］ Mzali F，Sassi L，Ben Nasrallah S.Optimal experiment design and simulations identification of thermo-physical properties of orthotropic solids 4th［C］.Rio de Janeiro：Int Conf on Inverse Problems in Eng，2002.

［5］ Huang C H，Wang SP.A three-dimensional inverse heat conduction problem in estimating surface heat flux by conjugate gradient method［J］.Int J of Heat and Mass Transfer，1999，42：3387-3403.

［6］ Yao J，Chia C.Comparison of Newton-Gauss with Levenberg-Marquardt Algonthm for Space Resection［C］.Singapore：The 22nd Asian Conference on Remote Sensing，2001.

［7］ 肖庭延，于慎根，王彥飛.反問題的數(shù)值解法［M］.北京：科學(xué)出版社，2003.

［8］ AndrieuC，De Freitas N，Doucet A，et al.An introduction to MCMC for machine learning［J］.Machine Learning，2003，50：5-43.

［9］ Wang J，Zabaras N.A Bayesian inference approach to the inverse heat conduction problem［J］.Int J of Heat and Mass Transfer，2004，47：3927-3941.

［10］ Wang J，Zabaras N.Hierarchical Baysesian models for inverse problems in heat conduction Issue［J］.Inverse Problem，2005，21（1）：183-206.

［11］ Wang J，N Zabaras.Using Bayesian statistics in the estimation of heat source in radiation Int［J］.Heat Mass transfer，2005，48：1529-1542.

［12］ 童景山.流體熱物性學(xué)-基本理論與計(jì)算［M］.北京：中國石化出版社，2008.

［13］ Lee J，Jung S Y，Sung H J，et al.Turbulent thermal boundary layers with temperature-dependent viscosity［J］.Int Jof Heat and Fluid Flow，2014.

［14］153 毛宇飛，白博峰，郭烈錦.混合長度理論的改進(jìn)及其在變物性傳熱中的應(yīng)用［J］.工程熱物理學(xué)報(bào)，2011，32（1）：58-62.

［15］ 陶文銓.數(shù)值傳熱學(xué)［M］.西安：西安交通大學(xué)出版社，2009.

［16］ Tarantola A.Inverse Problem Theory and Methods for Model Parameter Estimation［M］.Pennsylvania：SIAM，2005.

A Research on ldentifying Thermophysical Parameter of Pipe Turbulence Fluid by lnverse Heat Transfer Problem

MEI Ning，JIAO Si
（College of Engineering，Ocean University of China，Qingdao 266100，China）

A method to identify the thermophysical parameters by measuring the pipe convection heat transfer temperature field was presented in this paper.The governing equations for turbulent forced convection heat transfer in pipe were simplified according to the boundary-layer theory.and an improved mixing length turbulence model was developed by introducing the effect of density fluctuations into the equations of turbulent transport，which can be applied to calculate the temperature field and velocity field inside the pipe.Using Levenberg-Marquardt method got the viscosity coefficient in the control equations by measured temperature field.The results show that IHTPis feasible for identifying the thermophysical parameters and has a higher precision.

Levenberg-Marquardt method；mixing length turbulent model；turbulent convective heat transfer；turbulent viscosity coefficient

TK124

A

1672-5174（2014）10-142-06

責(zé)任編輯 陳呈超

國家自然科學(xué)基金項(xiàng)目（51276174）；青島市科技發(fā)展計(jì)劃項(xiàng)目（12-1-3-77）資助

2014-06-20；

2014-08-20

梅 寧（1961-），男，教授，博導(dǎo)。E-mail：nmei＠ouc.edu.cn