管內(nèi)充分發(fā)展流動與傳熱數(shù)值模擬的教學方法探討1)

毛宇飛 上官燕琴 肖 洪

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

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

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

1 邊界層積分控制方程組

1.1 物理過程

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

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

1.2 動量方程

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

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

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

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

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

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

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

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

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

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

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

1.3 能量方程

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

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

1.3.1 均勻熱流工況

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

式中，Um為半徑r流動截面上流體的平均無量綱速度，計算公式為

定義無量綱溫度Θ

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

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

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

1.3.2 均勻壁溫工況

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

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

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

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

2 數(shù)值計算方法

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

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

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

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

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

2.2 速度場和溫度場的求解

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

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

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

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

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

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

3.1 層流流動與傳熱

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

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

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

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

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

3.2 湍流流動與傳熱

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

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

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

4 結(jié)論

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