單邊突擴通道中kε方法數值模擬^>[1]

一、引言

如圖1所示的單邊突擴通道形式，是在燃燒室設計、電廠冷卻水研究等工程問題上經常遇到的流動現象。由于流動具有分離，湍流結構復雜，不宜采用傳統的近區積分模型。因為這種積分模型需要假定射流邊緣上的摻入率及流速、溫度的相似剖面，經驗性較大；在確切判斷半經驗公式的系數上，遇到不少困難。近年來，人們開始應用湍流模型來處理這類問題，它可以算得所論區域的湍流場、摻入率及速度、溫度剖面而無需經驗系數的事先假定。這種方法已由單方程模型發展到雙方程模型、多方程模型，是國際上迅速發展的一種新的理論和方法。在上述模型中，k-ε雙方程湍流模型應用最為廣泛。對浮射流、淹沒射流、尾流、回流、次生流等問題，無論定常或非定常、無論二維或三維，都已有若干成功的算例。結果令人鼓舞，預示了這類模型強大的生命力。

二、數學模型

不可壓縮流體作平面二維定常運動，當忽略質量力的影響時，其控制方程可以表示為：

稱為有效黏性系數，由分子黏性系數及湍流黏性系數所組成。其中k、ε分別稱為湍流動能及湍流動能耗散率，它們由以下兩輸運方程確定：

為湍流動能的產生項。

不考慮散熱損失，能量方程可以表示成一無源的對流擴散方程：

方程（1）～（8）組成一封閉系統。若要考慮密度ρ隨時間T的變化，例如溫差異重流問題，則還應加入狀態方程ρ=ρ（T）聯立求解。

數學模型中的六個通用常數，系由基本實驗確定，根據文獻［2］，它們的值見表1。

三、數值計算

差分方程的推導及“壁函數”的應用在文獻［3］中已有詳細敘述，本文擬就P-V迭代計算中出現的速度修正方程及壓力修正方程作一推導，并對CHAMPION程序的初值設定問題提出新的方法。

1.P-V迭代

本文采用斯波爾丁學派創立的SIMPLE算法，以壓力P和速度V為基本變量進行P-V迭代計算。即先假定全場壓力分布，使運動方程能獨立求解。這時解出的速度值自然不滿足連續方程，應通過速度修正方程來修正，將修正后的速度值代入壓力修正方程又給出新的壓力值。如此往復迭代，即能期望得到真實的解。

2.速度修正方程

參見圖2，U分量差分方程為：

式中：b為除壓力梯度外的源項。

V分量差分方程為：

式（9）與式（10）的A_i是分別對U、V單元而言的。

以初設壓力P^*代入上兩式得：

令

式（9）-式（9′）及式（10）-式（10′）得：

為了簡化計算，通常略去（參見文獻［4］），故得：

式中

3.壓力修正方程

連續方程（1）的差分形式為：

參見圖2，知：

將式（11）、式（12）、式（14）代入式（13），整理即得關于壓力修正P′的方程：

其中

式（15）即為所求的壓力修正方程，它具有統一形式，可以與其他方程一樣求解。

式（16）的B是壓力修正方程的源項，它具有連續方程的形式。若B=0，說明U^*、V^*滿足連續方程，顯然它們又滿足運動方程，因此必然就是精確解。通常稱B為質量源，以全場各點B值的絕對值之和趨于零作為迭代收斂的判據。

4.邊界條件

參見圖1，各變量的邊界條件可表述為：

AB面：

k=k₀=1.5×（0.05×U₀）²

DE面：

其他各面：

U=V=0

T=T_W

式中：T_W為壁面溫度。至于k、ε兩個變量，則采用“壁函數”方法處理>^［3］。

5.迭代方法及初值設定

迭代方法一般有點迭代與線迭代兩種。點迭代計算較穩定，但太費機時；本文采用線迭代，在每一條線上使方程的系數矩陣成為三對角矩陣形式，可用通用追趕法程序求解。線迭代不如點迭代穩定，為此，應用低松弛法提高穩定性和收斂性。

變量全場初值是否合理，直接影響到計算收斂問題。因為完全不合理的初值，可能一開始就使問題超出收斂半徑的范圍之外，計算無法進行。對初值設定問題，通常的辦法是用粗網格或令交換系數為常數進行粗略計算，將所得結果作為初值。這種做法，不僅同樣存在不收斂的危險，同時還使數據的輸入輸出變得很復雜，增加許多工作量。本文在迭代計算的前若干步，控制變量不超過其可能最大值，保證計算可以進行。經過這若干步的計算，變量初值很快趨于合理。此時取消其控制，任其自由收斂。此法對網格過分不均勻或算法上有本質錯誤的問題不起任何作用，即使每一步都加以控制，其余源及質量源也永遠不會變小，不可能得到收斂解。

收斂判據以相對質量源的絕對值小于10^-2為準，在M-160機上，需CPU時間約6min。

四、計算結果及其分析

本文對標準CHAMPION程序進行改進，將原來的單向掃描改為來回雙向掃描。這對回流問題，因為反向掃描可以迅速地將下游影響傳遞到區域各點，大大提高了計算的穩定性與收斂性。另外，去掉原來的逐線總體連續修正，僅保留出口處一條線不變，計算反而容易收斂。

定義進口雷諾數Re為：

式中：U_m為進口平均流速；ν為水的運動黏性系數；H為入口高度。

為了便于與計算結果^>[2]>[3]及實驗資料>^［5］對比，分別取Re=5322、3441、1803進行計算。圖3～圖5是各個工況下本文的計算結果、文獻［1］的計算結果及文獻［5］的實驗結果。圖中各量系下列無量綱數：

1.流線圖及時均速度分布

由圖3（a）、圖4（a）、圖5（a）可知，本文計算所得回流區長度在臺階下游7h～9h之間，詳見表2。

由表2可以看出，本文計算結果稍偏大。這一方面是由于數學模型的常數c₁對L／h值有一定影響，另一方面是實驗資料的精度亦不太高。文獻［5］的L／h值是由10個斷面的X方向時均速度分布內插而得；文獻［2］的結果是由雷諾數Re為5900和8340的L／h值線性外推得來。但事實上，當Re高達一定數值時，L/h值與Re無關，線性外推不合理。因此，文獻［2］中，Re=5322時的L/h值顯然應比7.0為大。

在圖3（b）、圖4（b）、圖5（b）中，分別給出了所取三個雷諾數情況下，十個橫斷面上的時均速度分布及文獻［1］、［5］的結果。由圖可知，除進口附近由于存在三維分離區而使回流區高度較實驗值小外，其他的地方均符合良好，并且比文獻［1］的計算結果更符合實驗資料，特別是在上壁面附近流速最大值處更為明顯。

2.湍流動能k、湍流動能耗散率ε及渦黏性系數μ_t分布

圖3（c）、圖4（c）、圖5（c）分別給出了三個雷諾數下的k剖面曲線，圖3（d）、圖4（d）、圖5（d）分別給出了ε剖面曲線，圖3（e）給出了Re=5322時的μ_t剖面曲線。剖面圖中同時還給出了文獻［1］的計算結果。由圖可見，吻合情況良好，只是ε在壁面附近符合稍差，這顯然是對ε的壁函數的不同處理而引起的。

3.相對壓力及溫度分布

圖5（e）給出了Re=1803時的相對壓力剖面曲線，圖4（e）給出了Re=3441時的等溫線。這里所指的壓力，是相對于進口處的壓力。由圖可見，相對壓力是在回流區最小。雷諾數越大（亦即進口速度越大），相對壓力越易出現負值，即越易引起通道內部出現低壓區，這與高速水流容易引起氣蝕的現象是吻合的。對溫度的計算，由于本文所選通道尺寸為長×高=20×2cm²，區域范圍很小，若按一般絕熱要求取作為壁面條件，則計算的溫度分布幾乎維護進口溫度不變。作為數學模型及程序的驗證，本文取T_W=0（即壁面處溫度等于自然水溫）作為邊界條件來計算溫度。由圖4（e）可見，溫度亦是在回流區最小。

Numerical Simulation of a Single Side Sudden Enlargement Channel with a k-ε Turbulence Model

Abstract

參考文獻

>［1］John F.Ripken，Norio Hayakawa.Cavitation in High-head Conduit Control Dissipators.ASCE，HY₁，1972：239-256.

［2］Rodi，W..Turbulent Models and Their Application in Hydraulics.A state of the Art Review，SFB 80/T/127，May，1978.

>［3］王能家.河道排取水問題的k-ε方法數值模擬∥水利水電科學研究院論文集（第26集）.北京：水利電力出版社，1986.

［4］Suhas V.Patankar.Nemerical Heat Transfer and Fluid Flow，McGraw Hill，New York，1980.

>［5］沈熊，廖湖生.應用頻移激光測速系統測量高湍流度回流流動.力學學報，1982（5）：505-511.

---