毛宇飛 上官燕琴 肖 洪

(河海大學(xué)機(jī)電工程學(xué)院，江蘇常州 213022)

筆者多年來從事能源與動(dòng)力工程專業(yè)的本科教學(xué)工作，在指導(dǎo)學(xué)生進(jìn)行本科畢業(yè)設(shè)計(jì)過程中，發(fā)現(xiàn)學(xué)生學(xué)習(xí)了流體力學(xué)、工程熱力學(xué)和傳熱學(xué)等專業(yè)基礎(chǔ)知識(shí)，也學(xué)習(xí)了計(jì)算機(jī)編程語言，但在將程序設(shè)計(jì)應(yīng)用于工程計(jì)算這方面卻掌握得非常薄弱。因此，筆者結(jié)合教學(xué)實(shí)踐和研究方向，給本科學(xué)生開設(shè)了一門專業(yè)提升課程——熱工計(jì)算與程序設(shè)計(jì)。與本科引入計(jì)算流體力學(xué)（computational fluid dynamics，CFD）軟件教學(xué)[1]不同，該課程要求學(xué)生針對流體力學(xué)、熱力學(xué)和傳熱學(xué)領(lǐng)域的一些典型算例，進(jìn)行數(shù)學(xué)建模、程序設(shè)計(jì)及計(jì)算分析。通過該課程的學(xué)習(xí)，訓(xùn)練學(xué)生編寫MATLAB 程序，將流體物性計(jì)算、物理過程數(shù)值計(jì)算、數(shù)據(jù)處理及圖像繪制設(shè)計(jì)成一個(gè)整體[2-3]，課程的教學(xué)思路與同濟(jì)大學(xué)付小莉老師的工作[4]有相似之處。

管內(nèi)不可壓縮流體充分發(fā)展流動(dòng)與傳熱是一個(gè)很常見的工程問題，層流工況時(shí)存在理論解析解，湍流工況時(shí)存在近似解和經(jīng)典公式，是計(jì)算流體力學(xué)和數(shù)值傳熱學(xué)課程中一個(gè)很好的入門算例。目前，學(xué)生圍繞該算例進(jìn)行程序設(shè)計(jì)和數(shù)值模擬，通常是在研究生階段完成的[5]，這是因?yàn)橄嚓P(guān)計(jì)算數(shù)學(xué)理論超出了本科生的知識(shí)儲(chǔ)備。本文將邊界層積分方法應(yīng)用于管內(nèi)充分發(fā)展流動(dòng)與傳熱的數(shù)值模擬，從而簡化了數(shù)值求解過程；并基于數(shù)值積分和有限差分法（或有限容積法），發(fā)展出適用于層流和湍流工況下管內(nèi)對流傳熱的數(shù)值計(jì)算方法。論文對該數(shù)值方法的數(shù)學(xué)原理和程序設(shè)計(jì)方案進(jìn)行了詳細(xì)闡述，對該工程算例的教學(xué)過程進(jìn)行了較完整的教學(xué)設(shè)計(jì)。

1 邊界層積分控制方程組

1.1 物理過程

當(dāng)流體進(jìn)入管道與管壁表面接觸時(shí)，在黏性摩擦的作用下，近壁處的流體速度變慢，形成邊界層。沿流動(dòng)方向邊界層逐漸增厚，而管中心的無黏流動(dòng)區(qū)則逐漸縮小，直至邊界層充滿整個(gè)管道，形成流動(dòng)充分發(fā)展段。層流充分發(fā)展區(qū)速度場遵循Poiseuille 拋物線分布規(guī)律，湍流充分發(fā)展區(qū)速度場則遵循冪律（power-law）分布規(guī)律[6-8]。將邊界層概念從速度邊界層推廣至溫度邊界層，這奠定了對流傳熱分析的理論基礎(chǔ)。

邊界層積分求解方法已被成功應(yīng)用于不可壓縮流體外掠平板（直角坐標(biāo)系）層流對流傳熱，引入合理假設(shè)后可采用數(shù)學(xué)方法直接求出摩擦阻力系數(shù)和對流傳熱系數(shù)的近似解[6-8]。本文將邊界層積分法引入管內(nèi)充分發(fā)展流動(dòng)與傳熱（圓柱坐標(biāo)系），展開理論分析與數(shù)值計(jì)算。建立邊界層積分方程組有兩種方法，一種是對邊界層微分方程沿邊界層厚度積分得到，另一種是對坐標(biāo)系內(nèi)控制體進(jìn)行質(zhì)量、動(dòng)量和能量守恒分析而導(dǎo)出。兩種方法殊途同歸，具體推導(dǎo)過程可參閱相關(guān)文獻(xiàn)中關(guān)于外掠平板的邊界層積分方程組[6-8]。

1.2 動(dòng)量方程

在管內(nèi)流動(dòng)充分發(fā)展區(qū)，對于不可壓縮流體，沿流動(dòng)方向（x坐標(biāo)方向）的速度u分布不再變化，即?u/?x=0，同時(shí)徑向（r坐標(biāo)方向）速度v可以忽略，這意味著穩(wěn)態(tài)軸對稱條件下連續(xù)性方程自動(dòng)滿足。此外，沿流動(dòng)方向的壓力梯度可視為常數(shù)。因此，流動(dòng)方向的動(dòng)量方程可簡化為

式中，τ 為剪切應(yīng)力，τw為壁面剪切應(yīng)力；層流流動(dòng)時(shí)μeff=μ，μ為流體動(dòng)力黏性系數(shù)；湍流流動(dòng)時(shí)μeff=μ+μt，μt為湍流（動(dòng)力）黏性系數(shù)；R為管道半徑；y為離開管壁面的距離,y=R-r。

定義無量綱距離Y，無量綱半徑η 和無量綱速度U為

式中，d為管道直徑，d=2R。將式(2)代入式(1)化簡得到

式(3)即為采用積分方法得到的動(dòng)量（無量綱速度）方程，這是一個(gè)一階偏微分方程，其邊界條件為：Y=0，U=0；Y=1，方程自動(dòng)滿足軸對稱條件。

工程計(jì)算中，可以選擇簡單實(shí)用的混合長度理論來計(jì)算圓管內(nèi)充分發(fā)展區(qū)的湍流黏性系數(shù)[5-8]

式中，ρ 為流體密度，l表示混合長度。將式(4)無量綱化

式中，Re為Reynolds 數(shù)，v為流體運(yùn)動(dòng)黏性系數(shù)（v=μ/ρ），ub和Ub分別為整個(gè)管道流動(dòng)截面上流體的平均速度和平均無量綱速度。不可壓縮流體的Ub計(jì)算公式為

黏性底層外的混合長度l可根據(jù)Nikurades公式計(jì)算[5]；而在近壁處的黏性底層中，需引入Van Driest 阻尼函數(shù)DF來修正分子黏性對湍流脈動(dòng)的影響[5]。因此，一個(gè)適用于圓管充分發(fā)展湍流邊界層的混合長度計(jì)算式為

式中，y+為無因次距離。

根據(jù)壁面剪切應(yīng)力τw與Darcy 摩擦阻力系數(shù)f之間的關(guān)系，導(dǎo)出f的計(jì)算式

1.3 能量方程

基于穩(wěn)定流動(dòng)開口系統(tǒng)能量方程，對管內(nèi)充分發(fā)展對流傳熱溫度邊界層進(jìn)行積分并無量綱化得到

式中，q為（半徑r截面）熱流密度，qw為壁面熱流密度，T為流體的溫度。工程中管內(nèi)對流傳熱通常存在兩種典型工況：均勻熱流邊界和均勻壁溫邊界。下面就這兩種不同工況對式(9)進(jìn)行化簡。

1.3.1 均勻熱流工況

若流體熱物性隨溫度的變化忽略不計(jì)，進(jìn)入熱充分發(fā)展段后，對流傳熱系數(shù)將保持常量。因此，在均勻熱流條件下，根據(jù)流體溫度變化關(guān)系[6-8]可將式(9)轉(zhuǎn)化為

式中，Um為半徑r流動(dòng)截面上流體的平均無量綱速度，計(jì)算公式為

定義無量綱溫度Θ

式中，Tw為壁面溫度，λ 為流體導(dǎo)熱系數(shù)。根據(jù)熱流密度的定義，將式(12)代入式(10)化簡得到均勻熱流條件下的能量（無量綱溫度）方程

式中，Pr為流體的Prandtl 數(shù)；Prt為湍流Prandtl 數(shù)。式(13)的邊界條件：Y=0，Θ=0；Y=1，方程自動(dòng)滿足軸對稱條件。工程計(jì)算中，Prt通常取0.9～1.0 之間的定值[5]。考慮到不同流體Pr數(shù)的變化，本文推薦采用下面經(jīng)驗(yàn)方法[9]計(jì)算Prt

由式(14)，Pr=1 時(shí)，Prt=1；Prt隨Pr的增大而減小。應(yīng)用式(14)，可將本文的數(shù)值方法推廣至低Pr數(shù)工況。

1.3.2 均勻壁溫工況

在均勻壁面溫度條件下，根據(jù)熱充分發(fā)展區(qū)的溫度變化關(guān)系[6-8]，對式(9)進(jìn)行化簡并無量綱化，得到均勻壁溫工況的能量（無量綱溫度）方程

式(15)的邊界條件與式(13)完全一致。式中，Θm為半徑r流動(dòng)截面上流體的平均無量綱溫度，Θb為整個(gè)流動(dòng)截面上流體的平均無量綱溫度，定義為

Θb與表征對流傳熱系數(shù)h的Nusselt 數(shù)Nu之間存在確定的數(shù)學(xué)關(guān)系

綜上所述，通過邊界層積分方法，將表征管內(nèi)充分發(fā)展流動(dòng)與傳熱的二階偏微分方程組簡化成邊界條件完全一致的一階偏微分方程組，從而可以使數(shù)值計(jì)算過程得到明顯簡化。

2 數(shù)值計(jì)算方法

2.1 網(wǎng)格生成及控制方程離散

根據(jù)前面推導(dǎo)的邊界層積分控制方程組可知，管壁處無量綱變量U或Θ均等于零，故數(shù)值計(jì)算可從管壁向管心推進(jìn)。采用內(nèi)節(jié)點(diǎn)法對無量綱坐標(biāo)Y方向計(jì)算區(qū)域進(jìn)行離散，將第1 個(gè)節(jié)點(diǎn)布置在管壁處（Y=0），最后一個(gè)節(jié)點(diǎn)（M1節(jié)點(diǎn)）布置在管心（Y=1）?？紤]到近壁處區(qū)域內(nèi)速度場和溫度場變化劇烈，此處網(wǎng)格應(yīng)布置稠密一些（越靠近管壁，網(wǎng)格越密）。通?？梢圆捎枚ū燃用芑蛑笖?shù)加密兩種方法生成網(wǎng)格[5]，兩種方法均可得到滿意的數(shù)值模擬結(jié)果。

采用邊界層積分法得到的動(dòng)量方程和能量方程均為一階偏微分方程，且壁面處（節(jié)點(diǎn)1）U和Θ均等于0，因此可以采用數(shù)值積分方法對控制方程進(jìn)行離散，得到求解節(jié)點(diǎn)2～M1處U和Θ的代數(shù)方程組?？刂品匠探M對應(yīng)的離散方程可以采用如下統(tǒng)一形式

式中，I為節(jié)點(diǎn)編號(hào)，I=1～M1；YCV表示控制容積（I=2～M1-1）的寬度，YCV(2)=YCV(M1)=0；針對不同的流動(dòng)與傳熱工況，Φ，A，B和C對應(yīng)的表達(dá)式或數(shù)值見表1。

表1 式(18)中的Φ，A，B 和C

2.2 速度場和溫度場的求解

對于不可壓縮（常物性）流體，流動(dòng)過程和傳熱過程相互獨(dú)立，因此可以單獨(dú)求解動(dòng)量方程獲得速度場分布，再將速度場分布代入能量方程求出溫度場分布。

層流流動(dòng)時(shí)，根據(jù)式(18)可依次求解出節(jié)點(diǎn)2 至節(jié)點(diǎn)M1的無量綱速度，無需迭代。湍流流動(dòng)時(shí)，由于湍流黏性系數(shù)取決于速度場分布，因此應(yīng)用式(18)數(shù)值求解各節(jié)點(diǎn)U時(shí)需要迭代計(jì)算；計(jì)算過程中，湍流黏性系數(shù)式(5)中的無量綱速度梯度?U/?Y采用二階精度的差分表達(dá)式；迭代收斂的判據(jù)是節(jié)點(diǎn)相鄰兩次計(jì)算值的最大相對偏差的絕對值小于10-6。求得各節(jié)點(diǎn)處的無量綱速度U，并獲得網(wǎng)格獨(dú)立解后，根據(jù)式(6)通過數(shù)值積分求得Ub，再根據(jù)式(8)求得摩擦阻力系數(shù)f。

求出無量綱速度場后，根據(jù)式(11)采用數(shù)值積分求出各節(jié)點(diǎn)對應(yīng)半徑流動(dòng)截面上流體的平均無量綱速度值Um，進(jìn)而代入能量方程進(jìn)行計(jì)算。對于均勻熱流條件，根據(jù)能量方程式(13)，將Um，Ub和μt/μ（湍流工況）代入式(18)直接求解出各點(diǎn)的無量綱溫度，無需迭代。對于均勻壁溫條件，求解能量方程式(15)時(shí)，不僅需要確定速度場，還需要獲得溫度場信息；因此計(jì)算時(shí)，先設(shè)定Θ迭代初場，根據(jù)式(16)通過數(shù)值積分求得Θb和Θm，再應(yīng)用式(18)計(jì)算求解出各節(jié)點(diǎn)的Θ，反復(fù)迭代，直至收斂，收斂條件與求解U方程一致。獲得網(wǎng)格獨(dú)立解后，再根據(jù)式(17)求得表征對流傳熱能力的Nu數(shù)。

根據(jù)上述分析，均勻壁溫條件下管內(nèi)湍流流動(dòng)與傳熱的數(shù)值模擬最為復(fù)雜，對應(yīng)速度場和溫度場均需迭代求解，其具體的數(shù)值求解過程如圖1 所示。針對不同的流動(dòng)與傳熱工況，采用MATLAB 編寫計(jì)算程序和繪圖程序，對相關(guān)計(jì)算結(jié)果進(jìn)行比較分析。

圖1 均勻壁溫管內(nèi)湍流流動(dòng)與傳熱數(shù)值求解過程

3 計(jì)算結(jié)果與分析

3.1 層流流動(dòng)與傳熱

計(jì)算表明，基于邊界層積分方程的算法具有較快的收斂速度。圖2 給出了層流工況下網(wǎng)格數(shù)對摩擦阻力系數(shù)f和傳熱Nu數(shù)的影響。在計(jì)算速度場時(shí)，即使網(wǎng)格（控制容積）數(shù)為1，此時(shí)內(nèi)部節(jié)點(diǎn)的計(jì)算值與精確值相差較大，但由于整個(gè)邊界層滿足總體的動(dòng)量守恒，Ub的計(jì)算值等于精確值，因此阻力系數(shù)預(yù)測值始終與解析解（f·Re=64）吻合。隨著網(wǎng)格數(shù)的增大，速度場不斷趨近精確值，溫度場也逐漸收斂；當(dāng)網(wǎng)格數(shù)大于30，層流對流傳熱Nu數(shù)已基本與精確值（均勻熱流Nu=4.364，均勻壁溫Nu=3.658）吻合。根據(jù)獲得的網(wǎng)格獨(dú)立解情況，層流工況的網(wǎng)格數(shù)取50，湍流工況（Re＜106）的網(wǎng)格數(shù)取200。

圖2 層流工況下網(wǎng)格數(shù)對計(jì)算結(jié)果的影響

圖3 和圖4 分別給出了管內(nèi)層流充分發(fā)展條件下無量綱速度場和無量綱溫度場的模擬結(jié)果。速度場與解析解、均勻熱流工況溫度場與解析解非常吻合。均勻壁溫條件下的溫度分布比起均勻熱流條件下的更不均勻，因此均勻壁溫工況的管道截面上Θb值較大，相應(yīng)Nu數(shù)較小。

圖3 管內(nèi)層流充分發(fā)展速度場模擬結(jié)果

圖4 管內(nèi)層流充分發(fā)展溫度場模擬結(jié)果

3.2 湍流流動(dòng)與傳熱

已有研究[6-8]表明，管內(nèi)湍流充分發(fā)展工況下，在適中的Re數(shù)范圍內(nèi)，速度場和溫度場近似呈1/7 次律指數(shù)型分布規(guī)律，圖5 給出的模擬結(jié)果很好地呈現(xiàn)出這一特征。

圖5 管內(nèi)湍流充分發(fā)展速度場和溫度場模擬結(jié)果

在Re=104～106范圍內(nèi)，將數(shù)值計(jì)算結(jié)果與有代表性的管內(nèi)充分發(fā)展湍流摩擦阻力和傳熱關(guān)聯(lián)式[8,10]（見表2）進(jìn)行對比。如圖6，摩擦阻力系數(shù)的數(shù)值計(jì)算值與Filonenko 關(guān)聯(lián)式計(jì)算值相當(dāng)一致，兩者間的相對偏差基本在 ±3%以內(nèi)，而Blasius 關(guān)聯(lián)式在超出其適用范圍后（Re＞105）與模擬結(jié)果相差較大。傳熱Nu數(shù)的數(shù)值計(jì)算結(jié)果如圖7 所示。圖7(a)表明：對于常規(guī)流體介質(zhì)（Pr＞0.5），均勻熱流和均勻壁溫條件下的模擬結(jié)果十分接近，均與Gnielinski 關(guān)聯(lián)式吻合較好；這說明當(dāng)流體Pr數(shù)較大時(shí)，傳熱熱阻乃至溫度變化主要集中在近壁處的黏性底層中，這兩種邊界條件對湍流傳熱的影響可以忽略。圖7(b)表明：對于低Pr數(shù)流體（通常為液態(tài)金屬），其較強(qiáng)的熱擴(kuò)散能力使得熱阻分布在整個(gè)流動(dòng)截面上，從而導(dǎo)致這兩種邊界條件下傳熱能力之間的差別不能忽視，均勻熱流條件的傳熱能力高于均勻壁溫條件（與層流工況類似）；Gnielinski 關(guān)聯(lián)式的預(yù)測值明顯偏低，兩種邊界條件下的數(shù)值計(jì)算結(jié)果與相應(yīng)條件的Sleicher 關(guān)聯(lián)式基本吻合。計(jì)算結(jié)果表明，本文的數(shù)值模型在不同工況條件下均能取得比較精確的模擬結(jié)果。

圖7 管內(nèi)湍流充分發(fā)展Nu 數(shù)計(jì)算結(jié)果

綜上所述，基于邊界層積分法，對管內(nèi)充分發(fā)展流動(dòng)與傳熱進(jìn)行數(shù)學(xué)建模、程序設(shè)計(jì)和數(shù)值模擬。在過去的3 年里，筆者先后指導(dǎo)2017，2018 和2019 級(jí)3 名能動(dòng)專業(yè)本科生圍繞該算例完成了畢業(yè)設(shè)計(jì)，基本步驟是：（1）對管內(nèi)不可壓縮流體充分發(fā)展層流對流傳熱的解析解進(jìn)行理論推導(dǎo)，熟悉控制方程組；（2）基于邊界層積分法，引入無量綱變量，導(dǎo)出管內(nèi)流動(dòng)與傳熱的無量綱控制方程組；（3）對計(jì)算區(qū)域和控制方程組進(jìn)行離散；（4）分別對層流和湍流工況、均勻熱流和均勻壁溫條件下的管內(nèi)流動(dòng)與傳熱進(jìn)行程序設(shè)計(jì)和數(shù)值模擬；（5）改進(jìn)數(shù)學(xué)模型，將數(shù)值方法從不可壓縮流體推廣至超臨界壓力流體。經(jīng)過3 輪畢業(yè)設(shè)計(jì)，數(shù)值方法和計(jì)算程序不斷得到完善，計(jì)算內(nèi)容也不斷深入，這3 名同學(xué)的工作均獲得了河海大學(xué)優(yōu)秀畢業(yè)設(shè)計(jì)論文。在2024 年《熱工計(jì)算與程序設(shè)計(jì)》課程教學(xué)中，我們將系統(tǒng)地講授這部分內(nèi)容。同時(shí)，課題組不斷完善課程自編講義，為提升教學(xué)質(zhì)量、編寫出版教材奠定基礎(chǔ)。

4 結(jié)論

本文將邊界層積分方法與數(shù)值計(jì)算方法結(jié)合起來，通過編寫計(jì)算程序，對管內(nèi)充分發(fā)展流動(dòng)與傳熱進(jìn)行了數(shù)值模擬?？刂品匠探M概念清晰、形式簡潔、易于理解。數(shù)值計(jì)算方法主要由數(shù)值積分和迭代法構(gòu)成，易于掌握；數(shù)值計(jì)算結(jié)果與解析解和相關(guān)經(jīng)典公式進(jìn)行比較，證明了該數(shù)值方法的可靠性、準(zhǔn)確性和快速收斂性。通過該算例的實(shí)踐練習(xí)，可以幫助學(xué)生鞏固所學(xué)知識(shí)，鍛煉學(xué)生在熱工流體領(lǐng)域的自主編程能力、數(shù)據(jù)處理能力以及理論聯(lián)系實(shí)際的能力。本文的工作將邊界層理論、對流傳熱理論、湍流工程計(jì)算及MATLAB 程序設(shè)計(jì)、數(shù)值計(jì)算和圖像繪制結(jié)合在一起，為本科階段開展計(jì)算流體力學(xué)或數(shù)值傳熱學(xué)的入門教學(xué)提供了一種新的思路。