李一術，朱 波，劉家磊，房軍齡

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

0 引 言

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

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

為減少計算量，化簡數值模擬中的設置，嘗試采用單流體方法處理流固耦合傳熱問題。結合流固耦合傳熱計算中解析法與數值模擬各自的特點，評價得出:流固耦合數值模擬能夠作為管道傳熱計算的基準解。以此基準，驗證單流體模擬計算結果。通過對模擬參數的改進，得到了單流體計算方式下較為理想的計算結果。說明單流體傳熱計算方法可行，可以應用到復雜的模型計算中。

1 流固耦合數值模擬的驗證

1.1 物理模型

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

圖1 物理模型Fig.1 Physical model

1.2 流固耦合數值模擬

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

連續(xù)方程:

動量方程(Navier-Stokes 方程):

其他變量的輸運方程:

式中:ρ 為流體密度;t 為時間;φ 為通用變量;Γ 為廣義擴散系數;S 為廣義源項。

1.2.2 網格劃分

根據模型自身特點，利用Ansys－Icem 網格劃分軟件，對管道的流體域及固體域都采用結構網格劃分。通過計算結果對網格敏感性分析，最終確定流體域網格數為3.1×105，固體域網格數為2.9 ×105。其網格劃分如圖2 所示。

圖2 管道網格劃分Fig.2 Mesh division of the tube

1.2.3 邊界條件

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

1.2.4 模擬結果

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

圖3 數值模擬計算結果Fig.3 Numerical simulation results

1.3 解析法對流固耦合數值模擬結果的驗證

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

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

對于管內流體:

對于管壁:

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

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

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

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

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

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

q(0，Z)=0。

當Z →∞時，θi有限。

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

式中:

這里:

氣流在管內的溫度變化幅度在20℃～50℃之間，通過查看空氣的熱物理性質表，發(fā)現其物性參數變化很小，基本可以看作是常物性。ρ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)。在解析法的計算中，氣體流動的未知參數為ho與hi。經計算，物理模型中的氣流Re=8.8 ×104，屬于過度流狀態(tài)。因此，對于加熱型的流體，使用最廣泛的為關聯(lián)式迪圖斯－貝爾特公式:

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

圖4 數值模擬結果與解析解對比Fig.4 The results of numerical simulation compared with analytical solution

2 單流體處理方法

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

式中:q 為熱流密度;λ 為管道材料導熱系數;r1為圓柱體內半徑;r2為圓柱體外半徑;tw1為圓柱體內壁溫;tw2為圓柱體外壁溫。

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

式中:To為指定的外部邊界溫度;Tw為流體域邊界溫度;qw為從外界流入流體域的熱流密度;hc為流體域邊界的表面?zhèn)鳠嵯禂担洳煌诒砻鎿Q熱系數，因為它不涉及流體的溫度。

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

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

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

3 單流體模擬計算結果及改進

流固耦合的數值模擬與解析法都能夠得到相對準確的結果，但解析法應用的條件更為苛刻，這是因為只有在流體溫度變化范圍小時，才可以利用其特點取常物性參數。對于數值模擬來說，應用范圍更廣，而且利用計算機強大計算能力，減輕了人員的工作量。

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

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

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

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

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

圖7 x 與y 的關系Fig.7 Relationship of x and y

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

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

圖8 改進后兩種數值模擬所得流經壁面熱流量Fig.8 The heat flux through the wall with two kinds of improved numerical simulation method

4 結 語

對于管道內流固耦合傳熱這類實際問題，解析法與流固耦合數值模擬法都能夠得到較為準確的計算結果。解析法通常只能針對簡單模型抽象處理，當實際問題較為復雜時，如求解艦船艙室內溫度場分布問題，流固耦合數值模擬更有其實用價值。

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

［1］張強，劉巨保，趙曉榮.換熱管流固耦合傳熱的分區(qū)求解數值算法［J］. 科學技術與工程，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］.南京航空航天大學學報，2004，36(1).

［4］王福軍. 計算流體動力學分析—CFD 軟件原理與應用［M］.北京:清華大學出版社，2004.

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

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