流動與傳熱數值計算通用控制方程的研究

李敬法，石國赟，陳宇杰，宇 波

(1.北京石油化工學院機械工程學院，深水油氣管線關鍵技術與裝備北京市重點實驗室，北京 102617；2.中國石油大學(北京)機械工程學院，北京 102249；3.西安交通大學能源與動力工程學院，陜西 西安 710049)

隨著計算機技術的進步和數值計算方法的發展，流動與傳熱數值模擬在工程領域和科學研究中得到越來越廣泛的應用。為進行模塊化編程和方便結果分析，通常將不同物理問題的控制方程寫成由非穩態項、對流項、擴散項和廣義源項組成的通用控制方程形式[1-2]。通用控制方程簡化了原方程的復雜形式，使不同類型的控制方程具有相似的表達形式，大大提高了相似程序模塊的通用性和整體編程效率，有利于探索不同問題的共同特征，深化理論研究。因此，通用控制方程深受流動與傳熱數值計算研究者的青睞[3-5]。

為了得到真實的數值模擬結果，首先需要保證物理問題控制方程的準確性。但目前通用控制方程多用于常物性參數并忽略摩擦效應的場合，針對變物性參數和考慮摩擦效應情形下的通用控制方程研究較少。在多數情況下，為了進一步簡化計算、降低編程實施的復雜度，研究者們通常忽略物性變化和摩擦效應對計算結果的影響。對于絕大多數常物性、黏度較小、流場變化不是很劇烈的物理問題，這種處理對最終計算結果的影響甚微。但對某些黏度變化大、流場變化劇烈、物性參數在空間存在突變或劇烈變化的物理問題，這種處理會對計算結果產生一定影響，嚴重時甚至會導致計算結果失真。如原油在大型浮頂油罐中儲存時，受太陽輻射、環境因素等影響，不同油層的溫度會存在一定差異。由于原油黏度受溫度的影響較大，導致油罐內不同位置處的原油黏度不同。若在大型儲罐溫度場的數值模擬研究中將油罐內原油黏度當作常數處理，計算結果會與實際情況產生一定偏差，最終對原油的周轉調度和加熱方案等產生影響[4，6]。

基于此，筆者對流動與傳熱通用控制方程進行了進一步深入研究，探討了通用控制方程是否考慮物性參數變化和摩擦效應對數值模擬結果的影響。重點分析了動量方程中由黏度變化引起的附加項和能量方程中由摩擦作用產生的附加熱源項對數值計算結果的影響規律，并通過2個經典算例對上述分析進行驗證和說明。

1 通用控制方程分析

流動與傳熱數值計算常用的通用控制方程形式為[7-8]：

(1)

式中：φ為通用變量，可表示速度、溫度、組分質量濃度等；ρφ為廣義密度；Γφ為廣義擴散系數；U表示速度矢量；Sφ為廣義源項，通常情況下不能歸入前3項的項均可歸入廣義源項中。當φ表示不同的變量時，ρφ、Γφ和Sφ的具體表達式也會發生相應變化。

以層流流動與換熱為例，連續性方程、動量方程和能量方程的表達式分別為：

(2)

μ(U)T]+

(3)

μ(U)T]·U+[μU+μ(U)T]·U-

(4)

將方程(2)～(4)改寫成通用控制方程(1)的形式：

(5)

(6)

(7)

其中：動量方程的廣義源項SU=·[μ(U)T]+·U)+ρf-p；能量方程的廣義源項ST=·[μU+μ(U)T]·U+[μU+μ(U)T]·U-·U-p)·U+Sh。

1.1 物性參數的影響

方程(5)～(7)中的物性參數主要有密度、黏度、比熱容、導熱系數等。文獻[3，7-8]中已專門討論過非常數比熱容的影響，這里不再贅述。由于工程實際中常見的流體多為不可壓縮流體，這里也不再討論密度的影響。因此，對于不可壓縮流體流動，式(6)中的廣義源項可進一步簡化為：

SU=·[μ(U)T]+ρf-p

(8)

為簡化計算和方便編程，研究者們通常忽略式(8)中由黏度變化產生的附加項·[μ(U)T]。分析式(8)可知，當黏度為常數時，·[μ(U)T]為0。因此，對于絕大多數流體黏度為常數或黏度變化非常小的物理問題，·[μ(U)T]對最終計算結果的影響甚微，基本可忽略。但當流體黏度變化劇烈時，這種近似處理會對計算結果產生一定影響，此時一般不能忽略·[μ(U)T]。

對于導熱系數λ，無論λ是否為常數，均對擴散項·(λT)和數值計算結果無影響。

1.2 摩擦效應的影響

對于不可壓縮流體流動，式(7)中的廣義源項可進一步簡化為：

ST=[μU+μ(U)T]·U+Sh

(9)

由式(9)可知，能量方程的廣義源項包含摩擦做功產生的附加熱源項和內熱源項。為簡化計算和便于編程，摩擦做功產生的附加熱源項[μU+μ(U)T]·U常被忽略。分析該項可知，對于流體黏度較小、流場變化不劇烈的物理問題，[μU+μ(U)T]·U的數值較小，因此對計算結果的影響較小。但當流體的黏度很大、流場變化劇烈時，[μU+μ(U)T]·U的值較大，會對計算結果產生影響，此時該項在數值計算中不應該被忽略。

2 算例驗證

2.1 物性參數的影響

以流動與傳熱數值計算的經典算例—封閉方腔自然對流為例，分析動量方程中由黏度變化產生的附加項對模擬計算結果的影響。這里需要指出的是，為了更好地說明黏度變化對計算結果的影響，本算例不考慮能量方程中摩擦效應的影響。計算區域示意圖如圖1所示，方腔邊長l=0.01 m，上、下邊界為絕熱邊界條件，左、右邊界為恒溫邊界條件，其中Th=60 ℃、Tc=0 ℃，所有邊界均為無滑移邊界。采用2種不同性質的流體作為方腔內的流動工質，物性參數如表1所示。采用64×64的均分網格劃分計算區域，采用有限容積內節點法離散控制方程，采用SIMPLE算法進行壓力-速度耦合計算，詳細數值求解步驟及算法介紹參見文獻[9-11]，這里不再贅述。

表1 流體物性參數表

對于自然對流問題而言，工程實際中最關注的是溫度場的分布。方腔溫度場等值線對比圖如圖2所示，其中紅色實線表示考慮黏度變化產生附加項時的計算結果，藍色虛線表示不考慮黏度變化產生附加項時的計算結果。從圖2中可明顯看出，對于本算例的2種流體工質，當流體黏度隨溫度發生變化時，考慮和不考慮黏度變化產生的附加項時得到的溫度等值線不完全重合，尤其是在中間位置處兩者偏差較明顯。說明黏度變化產生的附加項對計算結果有一定影響。因此，當對計算精度要求較高時，為得到準確的數值解，需考慮黏度變化產生的附加項對最終數值計算結果的影響。

方腔垂直中心線處的溫度分布對比圖如圖3所示。與圖2類似，2種情況的溫度分布不能完全重合，說明黏度變化產生的附加項對計算結果會產生一定影響。此外，還可觀察到相同位置處考慮黏度變化產生附加項時的溫度值要高于不考慮黏度變化產生的附加項時的溫度值，原因是黏度變化產生的附加項相當于一個內熱源，在計算區域內產生了熱量。所以，考慮黏度變化產生的附加項時，計算區域流體的溫度會升高。

為定量說明考慮黏度變化產生的附加項影響，計算了圖3方腔垂直中心線溫度分布曲線的平均相對誤差和最大相對誤差，結果如表2所示。其中，平均相對誤差和最大相對誤差的計算式分別為:

(10)

(11)

式中：Ngrid表示方腔垂直中心線方向的節點數；Tv表示考慮黏度變化附加項時的溫度；T0表示不考慮黏度變化附加項時的溫度。

由表2中的數據可知，黏度變化產生的附加項對數值計算結果會產生一定影響，但這種影響不大。本算例方腔垂直中心線上的溫度最大平均相對誤差和最大相對誤差分別為1.26%和1.78%。而且還可發現，當黏度變化劇烈時，平均相對誤差更大。

表2 方腔垂直中心線處溫度誤差

2.2 摩擦作用的影響

以流動與傳熱數值計算的另一個經典算例—封閉方腔混合自然對流為例，分析能量方程中摩擦作用產生的附加熱源項對模擬計算結果的影響。這里需要指出的是，為了更好地說明摩擦作用對計算結果的影響，本算例不考慮動量方程中物性參數非常數的影響。計算區域示意圖如圖4所示。方腔邊長l=0.01 m，上、下邊界為絕熱邊界條件，左、右邊界為恒溫邊界條件，其中Th=50 ℃、Tc=0 ℃，頂蓋拖動速度為ulid，其余邊界均為無滑移邊界。采用某流體作為方腔內的流動工質，物性參數如表3所示。采用與算例2.1完全相同的網格數、離散方法以及數值算法進行求解，這里不再贅述。

表3 流體物性參數表

方腔溫度場等值線對比如圖5所示，其中紅色實線表示考慮摩擦作用產生附加項時的計算結果，藍色虛線表示不考慮摩擦產生附加項時的計算結果。由圖5可知，在一定的方腔頂蓋拖動速度和較大的流體黏度下，是否考慮摩擦作用產生的附加項對計算區域溫度分布規律有較大影響，尤其是方腔中心位置處存在較大偏差。因此，當流場變化較劇烈且流體黏度較大時，需考慮摩擦作用產生的附加項對最終數值計算結果的影響。從圖5還可看出，不同的頂蓋拖動速度、流體黏度對數值計算結果的影響不同。

方腔垂直中心線處的溫度分布對比圖如圖6所示。對比圖6(a)和圖6(c)可發現，頂蓋拖動速度對兩者溫度等值線分布差異具有重要影響。在流體黏度不變的情況下，當頂蓋拖動速度從0.1 m/s提高到0.5 m/s時，流場的變化變得更加劇烈，考慮摩擦作用影響和不考慮摩擦作用影響的溫度分布偏差顯著增大。對比圖6(b)和圖6(c)可發現，在頂蓋拖動速度不變的情況下，流體黏度對兩者溫度等值線分布也具有重要影響。當流體黏度從5 Pa·s提高到15 Pa·s時，考慮摩擦作用影響和不考慮摩擦作用影響的溫度分布偏差增大，但增大的幅度比頂蓋拖動速度弱一些，這是因為此時流體黏度已經較大，黏度對計算結果的影響比頂蓋拖動速度弱一些。此外，考慮摩擦作用時的溫度值要高于不考慮摩擦作用時的溫度值，原因是摩擦作用產生的附加項相當于一個內熱源項，所產生的熱量提升了計算區域的整體溫度。因此，考慮摩擦作用產生的附加項時，計算區域流體的溫度會升高。

為定量說明考慮摩擦作用產生的附加項影響，計算了圖6方腔垂直中心線處溫度分布的平均相對誤差和最大相對誤差，如表4所示。由表4中可以看出，摩擦作用產生的附加項對數值計算結果產生重要影響，當頂蓋拖動速度為0.5 m/s、流體黏度為15 Pa·s時，方腔垂直中心線上的溫度最大相對誤差高達70.35%。而且當流場變化劇烈程度增大或流體黏度增大時，相對誤差會變得更大，在本算例條件下，流場變化對相對誤差的影響比流體黏度對相對誤差的影響更顯著。因此，當流體黏度較大且流場變化較劇烈時，數值計算應考慮摩擦作用產生的附加項，否則會使計算結果失真。

表4 方腔垂直中心線處溫度誤差

3 結論

分析和推導了通用控制方程中考慮物性參數變化和摩擦效應時的附加源項，探討了物性參數變化和摩擦效應對通過控制方程計算結果的影響，得出如下結論：

(1)目前常用通用控制方程一般只適用于黏度為常數的情況。當流體工質的黏度不是常數且流場變化較劇烈時，動量方程中由于黏度變化產生的附加項相當于一個內熱源項，會對計算結果產生一定影響。在進行精確數值模擬時，可考慮此項的影響。但相比摩擦作用的影響，黏度變化的影響相對較弱。

(2)目前常用通用控制方程一般未考慮摩擦作用的影響。當流體工質的黏度較大且流場變化劇烈時，能量方程中由于摩擦作用產生的附加項不可忽略。該項相當于一個內熱源項，會對計算結果產生較大影響，嚴重時甚至導致計算結果失真。