單流體方法處理管道流固耦合傳熱

李一術(shù)，朱 波，劉家磊，房軍齡

(1．海軍工程大學(xué) 核能科學(xué)與工程系，湖北 武漢430033;2．中國人民解放軍92730 部隊(duì)，海南 三亞572000)

0 引 言

大型艦船為了滿足動(dòng)力需求，通常以核動(dòng)力作為其動(dòng)力輸出來源。在特定的艙室中一般都包含大量管路，管道的主要作用是輸送流體和熱量，它的連接使各設(shè)備之間形成完整回路。在對(duì)艙室溫度場的計(jì)算分析時(shí)，各管路與艙室的傳熱對(duì)艙室內(nèi)溫度場的分布有重要的影響。計(jì)算管道傳熱，通常需要進(jìn)行流固耦合計(jì)算。目前有2 種流固耦合傳熱方法［1］:一種是強(qiáng)耦合方法，即整場離散整場求解方法;另一種是弱耦合方法，即分區(qū)求解邊界耦合方法。對(duì)于計(jì)算模型較大時(shí)，一般采用弱耦合計(jì)算方法。通常艙室內(nèi)的設(shè)備與管路眾多，弱耦合方式則要求設(shè)定相應(yīng)數(shù)量的流固交界面，提高了人員工作量。同時(shí)相比于普通的單流體流動(dòng)計(jì)算，流固耦合計(jì)算需要更多的時(shí)間，對(duì)計(jì)算機(jī)的性能也提出了挑戰(zhàn)。

針對(duì)管道流固耦合傳熱問題，Al －Nimr 和El －shaarawi［2］的解析方法解決了在2 個(gè)半無限大平行平板間或者在一個(gè)半無限大有限壁厚有限熱容的管道中，不可壓縮的熱充分發(fā)展中的穩(wěn)定一維的瞬時(shí)強(qiáng)制對(duì)流耦合熱交換問題。利用解析法求解穩(wěn)態(tài)條件下管道溫度分布，夏文慶等［3］得出其計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果吻合得很好。張強(qiáng)等［1］以整場離散整場求解為基準(zhǔn)，認(rèn)為在數(shù)值模擬中松弛因子取0.75、收斂容差取0.1 可以提高計(jì)算精度和求解速度。

為減少計(jì)算量，化簡數(shù)值模擬中的設(shè)置，嘗試采用單流體方法處理流固耦合傳熱問題。結(jié)合流固耦合傳熱計(jì)算中解析法與數(shù)值模擬各自的特點(diǎn)，評(píng)價(jià)得出:流固耦合數(shù)值模擬能夠作為管道傳熱計(jì)算的基準(zhǔn)解。以此基準(zhǔn)，驗(yàn)證單流體模擬計(jì)算結(jié)果。通過對(duì)模擬參數(shù)的改進(jìn)，得到了單流體計(jì)算方式下較為理想的計(jì)算結(jié)果。說明單流體傳熱計(jì)算方法可行，可以應(yīng)用到復(fù)雜的模型計(jì)算中。

1 流固耦合數(shù)值模擬的驗(yàn)證

1.1 物理模型

假設(shè)有一段長l = 3 m 的圓管道，管道內(nèi)半徑ri= 1.5 cm，外半徑ro= 2 cm，管道材料的導(dǎo)熱系數(shù)k2= 0.07 W/ (m·K)。管道內(nèi)部流動(dòng)的是空氣，導(dǎo)熱系數(shù)為k1，入口溫度To= 20℃，速度V =5 m/s。環(huán)境溫度Ti= 50℃，環(huán)境與圓管外壁的表面換熱系數(shù)為ho。為了計(jì)算的簡便，假定ho非常大，使得圓管外壁溫度與環(huán)境溫度趨于一致，這樣省略了圓管外壁與環(huán)境之間的流固耦合傳熱計(jì)算。其物理模型如圖1 所示。

圖1 物理模型Fig.1 Physical model

1.2 流固耦合數(shù)值模擬

1.2.1 基本的物理守恒定律建立控制方程［4］

連續(xù)方程:

動(dòng)量方程(Navier-Stokes 方程):

其他變量的輸運(yùn)方程:

式中:ρ 為流體密度;t 為時(shí)間;φ 為通用變量;Γ 為廣義擴(kuò)散系數(shù);S 為廣義源項(xiàng)。

1.2.2 網(wǎng)格劃分

根據(jù)模型自身特點(diǎn)，利用Ansys－Icem 網(wǎng)格劃分軟件，對(duì)管道的流體域及固體域都采用結(jié)構(gòu)網(wǎng)格劃分。通過計(jì)算結(jié)果對(duì)網(wǎng)格敏感性分析，最終確定流體域網(wǎng)格數(shù)為3.1×105，固體域網(wǎng)格數(shù)為2.9 ×105。其網(wǎng)格劃分如圖2 所示。

圖2 管道網(wǎng)格劃分Fig.2 Mesh division of the tube

1.2.3 邊界條件

本文所用的流體力學(xué)計(jì)算軟件為Ansys-CFX，該軟件具有先進(jìn)的全隱式耦合多網(wǎng)格線性求解器，且收斂速度較快，具有友好的人機(jī)交互界面。為了進(jìn)行模擬，需要設(shè)定計(jì)算參數(shù)與計(jì)算條件。入口邊界條件采用速度入口，出口邊界條件采用壓力出口。流入的空氣設(shè)為理想氣體，設(shè)置空氣流動(dòng)模型為k - ε 模型。流體域與固體域網(wǎng)格采用GGI 連接模式對(duì)理想氣體與管壁進(jìn)行流固耦合傳熱計(jì)算。

1.2.4 模擬結(jié)果

利用可視化軟件分析模擬結(jié)果，得到如圖3 所示結(jié)果。圖3 (a)為流體邊界層溫度場分布云圖，圖3 (b)為對(duì)管內(nèi)流體主流溫度隨軸向坐標(biāo)變化提取的數(shù)值分布。

圖3 數(shù)值模擬計(jì)算結(jié)果Fig.3 Numerical simulation results

1.3 解析法對(duì)流固耦合數(shù)值模擬結(jié)果的驗(yàn)證

夏文慶等通過實(shí)驗(yàn)的方法證明，解析法在穩(wěn)態(tài)條件下可以較好的預(yù)測實(shí)際結(jié)果，因此這里用解析法代替實(shí)驗(yàn)來驗(yàn)證流固耦合數(shù)值模擬的準(zhǔn)確性。

解析法中的能量方程為［2］:

對(duì)于管內(nèi)流體:

對(duì)于管壁:

式中:ρ1為管內(nèi)流體密度;ρ2為管壁固體材料的密度;c1為流體定壓比熱容;c2為固體比熱容;ri為管道內(nèi)半徑;ro為管道外半徑;V 為流體軸向平均速度;T1為流體平均溫度;T2為固體平均溫度;Ta為管道外環(huán)境溫度;hi為管道內(nèi)壁對(duì)流換熱系數(shù);ho為管道外壁對(duì)流換熱系數(shù);t 為時(shí)間;qd為熱流密度;z 為軸向坐標(biāo);k2為固體的導(dǎo)熱系數(shù)。

將式(4)～式(6)的方程無量綱化，得:

式(7)～式(9)中各參數(shù)定義:

式(7)～式(9)受到以下初始條件的限制:

θi(0，Z)=0，i = 1，2，

θi(τ，0)=θoi= 1，i = 1，2，

q(0，Z)=0。

當(dāng)Z →∞時(shí)，θi有限。

當(dāng)管道內(nèi)的流體傳熱處于穩(wěn)態(tài)時(shí)，將無量綱方程轉(zhuǎn)化為特征方程，解出其特征值。利用初始條件的限制，對(duì)特征值進(jìn)行根的判別，最終得出流體的溫度分布如下:

式中:

這里:

氣流在管內(nèi)的溫度變化幅度在20℃～50℃之間，通過查看空氣的熱物理性質(zhì)表，發(fā)現(xiàn)其物性參數(shù)變化很小，基本可以看作是常物性。ρ1=1.165 kg/m3，ρ2= 250 kg/m3，c1= 1 000 J/ (kg·K)，c2= 1 000 J/ (kg·K)，k1= 0.028 W/ (m·K)。在解析法的計(jì)算中，氣體流動(dòng)的未知參數(shù)為ho與hi。經(jīng)計(jì)算，物理模型中的氣流Re=8.8 ×104，屬于過度流狀態(tài)。因此，對(duì)于加熱型的流體，使用最廣泛的為關(guān)聯(lián)式迪圖斯－貝爾特公式:

分布函數(shù)做出的二維曲線如圖4 所示，圖4 中同時(shí)還加入了之前數(shù)值模擬的計(jì)算結(jié)果，以便兩者進(jìn)行比較。在管道流體初始發(fā)展階段，2 種方法計(jì)算得出的流體溫度有一定的偏差。沿著管道軸向發(fā)展，偏差逐漸變小，最終幾乎趨于一致。流體從進(jìn)入管口開始，需經(jīng)歷一段距離，管斷面流速分布和流動(dòng)狀態(tài)才能達(dá)到定型，初始偏差較大的原因，可能是由于解析解在這個(gè)問題上考慮有所欠缺，所以導(dǎo)致偏差較大。就解析解作為基準(zhǔn)解而言，數(shù)值模擬計(jì)算的流體整體溫度偏差范圍較小，因此，可以很大程度上信賴其計(jì)算結(jié)果。

圖4 數(shù)值模擬結(jié)果與解析解對(duì)比Fig.4 The results of numerical simulation compared with analytical solution

2 單流體處理方法

分析管道的傳熱，當(dāng)其長度遠(yuǎn)大于管道壁厚時(shí)，壁內(nèi)溫度可近似看成一維穩(wěn)態(tài)溫度場，所以圓管壁面導(dǎo)熱公式為［5］:

式中:q 為熱流密度;λ 為管道材料導(dǎo)熱系數(shù);r1為圓柱體內(nèi)半徑;r2為圓柱體外半徑;tw1為圓柱體內(nèi)壁溫;tw2為圓柱體外壁溫。

在對(duì)流體的傳熱計(jì)算中，CFX 允許流體域外圍有一層假想固體域［6］。如圖5 所示，內(nèi)部物理流體邊界至最外部物理邊界之間有一段厚度的墻壁，但其沒有在實(shí)際模型中進(jìn)行模擬。通過模型邊界的熱流量用下式計(jì)算:

式中:To為指定的外部邊界溫度;Tw為流體域邊界溫度;qw為從外界流入流體域的熱流密度;hc為流體域邊界的表面?zhèn)鳠嵯禂?shù)，其不同于表面換熱系數(shù)，因?yàn)樗簧婕傲黧w的溫度。

圖5 假想邊界模型Fig.5 Imaginary boundary model

在處理管道內(nèi)流固耦合傳熱計(jì)算時(shí)，單流體傳熱計(jì)算方法為:假想有一層特定厚度的管壁，實(shí)際省略對(duì)管壁的建模和流固熱耦合的計(jì)算。在單流體模擬計(jì)算中，To為已知條件，Tw由迭代計(jì)算而得出，而hc是需要初始設(shè)定的未知參數(shù)。根據(jù)一維圓柱筒壁面?zhèn)鳠岬慕?，可得hc為:

利用式(15)圓管道的具體尺寸，換算得到在單流體傳熱計(jì)算中需要設(shè)定的hc值。

3 單流體模擬計(jì)算結(jié)果及改進(jìn)

流固耦合的數(shù)值模擬與解析法都能夠得到相對(duì)準(zhǔn)確的結(jié)果，但解析法應(yīng)用的條件更為苛刻，這是因?yàn)橹挥性诹黧w溫度變化范圍小時(shí)，才可以利用其特點(diǎn)取常物性參數(shù)。對(duì)于數(shù)值模擬來說，應(yīng)用范圍更廣，而且利用計(jì)算機(jī)強(qiáng)大計(jì)算能力，減輕了人員的工作量。

為了驗(yàn)證單流體假設(shè)計(jì)算的合理性，做多組數(shù)值模擬方案進(jìn)行分析。在相同物理?xiàng)l件下，參考流固耦合方式的數(shù)值模擬結(jié)果，驗(yàn)證單流體假設(shè)計(jì)算的可行性。選取多組不同尺寸管道，其尺寸及相應(yīng)尺寸下計(jì)算所得hc如表1 所示。環(huán)境溫度Ti提升為100℃，其他物理?xiàng)l件不改變。

表1 管道尺寸及計(jì)算所得hcTab.1 Tube size and calculation value of hc

對(duì)于加熱型管道，流入流體域的熱流量是影響管內(nèi)流體溫度的關(guān)鍵因素。為了方便多組計(jì)算的比較，直接選用流經(jīng)管壁的熱流量進(jìn)行對(duì)比分析。經(jīng)過數(shù)值模擬，分別使用流固耦合方法和單流體方法計(jì)算出的總熱流量如圖6 所示?？梢钥闯?，使用單流體方法得出的總熱流量都普遍偏小。這可能是由于管道壁面導(dǎo)熱近似處理存在誤差所致。設(shè)管道外半徑與管道內(nèi)半徑的比值為x，流固耦合下的總熱流量與單流體得出的總熱流量比值為y。仔細(xì)分析單流體方式所帶來的管壁熱流量偏差，可以發(fā)現(xiàn):y 與x 之間存在一定的關(guān)系。其關(guān)系如圖7 所示。

從圖7 觀察發(fā)現(xiàn)，兩比值之間近似地服從線性關(guān)系，通過線性擬合的方式，得出2 個(gè)比值之間的數(shù)值關(guān)系如下:

圖6 兩種數(shù)值模擬方法所得流經(jīng)壁面熱流量Fig.6 The heat flux through the wall with two kinds of numerical simulation method

圖7 x 與y 的關(guān)系Fig.7 Relationship of x and y

單流體計(jì)算的熱流量由式(14)所決定，To為定值，要想改變qw，只能對(duì)hc進(jìn)行修改。當(dāng)調(diào)整hc時(shí)，Tw的值也會(huì)隨著邊界條件的改變而改變，那么qw的變化不能夠精確掌控。經(jīng)測試，hc在小范圍變化時(shí)，Tw在現(xiàn)有物理?xiàng)l件下，相比原有溫度變化幅度很小。因此小幅度的改動(dòng)hc值，Tw的變化對(duì)于To與Tw之間的差值影響很小。這樣，可以依照管道半徑比值與熱流量比值之間的關(guān)系來修改hc，以達(dá)到修正qw的目的。修改后的hc變?yōu)閔′c，其表達(dá)式如下:

利用h′c重新計(jì)算各方案通過管壁的熱流量，并再次與之前計(jì)算的流固耦合下熱流量進(jìn)行對(duì)比，得到如圖8 所示的結(jié)果。可以看到，利用這種方法很好的解決了單流體傳熱計(jì)算帶來的管壁熱流量偏差，使其計(jì)算結(jié)果與基準(zhǔn)結(jié)果更加接近。這也說明單流體方法處理管道流固耦合傳熱有效。

圖8 改進(jìn)后兩種數(shù)值模擬所得流經(jīng)壁面熱流量Fig.8 The heat flux through the wall with two kinds of improved numerical simulation method

4 結(jié) 語

對(duì)于管道內(nèi)流固耦合傳熱這類實(shí)際問題，解析法與流固耦合數(shù)值模擬法都能夠得到較為準(zhǔn)確的計(jì)算結(jié)果。解析法通常只能針對(duì)簡單模型抽象處理，當(dāng)實(shí)際問題較為復(fù)雜時(shí)，如求解艦船艙室內(nèi)溫度場分布問題，流固耦合數(shù)值模擬更有其實(shí)用價(jià)值。

在流固耦合數(shù)值模擬的基礎(chǔ)上，改進(jìn)模擬思路，提出了單流體數(shù)值模擬方法。經(jīng)驗(yàn)證，其同樣可行。應(yīng)用單流體數(shù)值模擬方法，不僅令設(shè)置更加簡單，而且減少了計(jì)算量，降低了對(duì)計(jì)算機(jī)硬件的要求。所以，在傳熱模型復(fù)雜時(shí)，單流體數(shù)值模擬方法更為有效。

［1］張強(qiáng)，劉巨保，趙曉榮.換熱管流固耦合傳熱的分區(qū)求解數(shù)值算法［J］. 科學(xué)技術(shù)與工程，2010，10(28):1671-1815.

［2］Al-Nimr M A，El-Shaarawi M. Analytical solutions for the transient conjugated heat transfer in parallel plate an circular ducts［J］.Int Comm Heat Mass Transfer，1992，19:869 -878.

［3］夏文慶，明曉，朱春玲.基于狀態(tài)空間的流固耦合傳熱的分析［J］.南京航空航天大學(xué)學(xué)報(bào)，2004，36(1).

［4］王福軍. 計(jì)算流體動(dòng)力學(xué)分析—CFD 軟件原理與應(yīng)用［M］.北京:清華大學(xué)出版社，2004.

［5］章熙民，任澤霈，梅飛鳴.傳熱學(xué)［M］.北京:中國建筑工業(yè)出版社，2007.

［6］ANSYS CFX - Pre User′ s Guide［M］. Pennsylvania:ANSYSY，2009.