浮力回流

方法數值模擬^>[1]

一、引言

火力發電廠的冷卻水取水水溫，是影響電廠出力的重要因素；受納水域的水溫升高，會影響到水的環境。核電站冷卻水有時還帶有效射性物質，受納水域除熱污染外，還存在放射性物質污染問題。因此，對水電、核電建設來說，預報水域水體的溫度分布與水體異物的濃度分布是十分必要的。

隨著計算技術的發展，數學模擬方法日益被人們所重視，它與物理模擬一樣，是環境預測技術必不可少的工具。數學模型一般分為整區模型與區域模型。目前普遍采用后者，即將受納水域按流動的特性劃分為近區與遠區分別加以模擬，在每一區內按各自的特點簡化控制方程。針對這兩個不同區域，已有了不少數學模型>^［1］。由于近區流動的復雜性，近區模型不及遠區模型成熟，不少問題尚有待解決。近區不但是污染最重的地區，直接反映著熱水及其他污源排放的物理性質，是研究工程及環境的主要關注所在，而且遠區模型也以近區計算成果作為其邊界條件。k-ε模型用于近區相當成功>^［2］，本文的突擴部分計算結果，也證實了這一點。

區域模型雖然有利于數學處理，經濟易行，但區域劃分還不明確，成果銜接亦存在不少問題。因此，整區模型仍應該是發展方向。k-ε方程為建立整區模型創造了條件。

在瀕臨江河地區，一般是直接引取河水作為冷卻水水源。圖1是這種直流取水或混合取水的工程布置示意圖。

這是一個典型的整區問題。20世紀60年代初，陳惠泉、許玉林對此作了勢流分析、電模擬及水槽試驗驗證>^［3］；在數學模擬方面，麥國克作過嘗試>^［2］，但未見正式發表。本文針對這一問題，用k-ε模型進行模擬，并與水槽試驗成果進行對比。結果基本合理，總的趨勢與試驗結果定性符合，證實了這類精細模型的可用性與合理性。在計算排取水問題前，為了驗證數學模型及數值方法的正確性，參照文獻［4］，先用此法計算單邊突擴問題，所得結果與文獻［4］的試驗結果符合良好，比用φ-ω方法的計算結果稍優^>[2]，詳見成果分析部分。

二、數學模型

考察淺水定常出流，棱柱體渠道，選用水深平均形式的方程。流體不可壓縮，不計流散、浮力及表面風應力的影響，自由面作剛蓋近似，數學模型具有下述形式：

其中

以上各式中：x、y為坐標變量；U、V為相應于x、y方向的水深平均時均速度；T為水深平均時均溫升（實際水溫與平衡水溫之差）；k、ε為水深平均形式的湍流動能及湍流動能耗散率，；c_f為經驗摩擦系數，取為0.003分別為渦黏性系數與渦擴散系數；c₁、c₂、c_μ、σ_k、σ_ε、σ_T為通用常數，它們皆由基本試驗確定，不依賴于特定的問題。參見文獻［2］，分別取值見表1。

本模型中的k、ε兩個輸運方程，是從納維—斯托克斯方程出發推導而來。解這兩個方程以確定全場各點的渦黏性系數與渦擴散系數，從而充分反映各點湍流紊動的差異對動量擴散與熱量（或質量）擴散的不同影響。這就是k-ε模型區別于一般常系數擴散模型的本質所在。

三、數值計算方法及程序說明

（一）統一形式的微分方程

控制方程中的每個方程都具有類似的項，為便于差分，可以表示成如下統一的形式：

或

式中：φ為任一因變量，如U、V、T、k、ε等；Γ_φ為φ變量對應的交換系數。

式（7）中〈1〉項為時間變率項，雖然考慮定常問題，但仍保留此項，為便于今后應用時間相關法；〈2〉項為對流項；〈3〉項為擴散項；〈4〉項為源項，但并非真實的物理源，而是將不能納入前三項的一切項都歸并于源項。

為清楚起見，表2列出了對應不同的φ方程式（7）中各項的具體內容。

表2中，μ_eff=ρυ_eff=ρ（ν+ν_t）=μ+μ_t，稱為有效黏性系數，它是分子黏性系數與湍流黏性系數之和。

基本方程組是一多元非線性方程組，各方程間皆具有耦合性，因此只能用數值方法求解，且必須使用迭代方法。

（二）差分方程及其邊界條件

1.網格系

建立一個差分方程，首先應將求解區域離散化，形成差分網格。本文采用交錯網格如圖2所示。

為了提高效率，網格在x、y方向皆取成非均勻，節點量隨L、D、B₁、B₂有所改變，一般為10×60或12×54個。排水及取水口處應密分，一般為3～5個節點。圖2中，X方向實線為U網格線，Y方向實線為V網格線，虛線及邊界實線為G變量（K、ε、T、h）網格線。

2.差分方程

采用五點差分格式，把變量φ在P點上的值φ_P表示成鄰節點上φ值的函數：

式中，E、W、N、S分別代表P點的東、西、北、南四點，見圖3。

使用控制容積法，將統一形式的微分方程在某一單元P上積分，可得

式中：為單元體體積、單元表面積矢量；為對全部單元邊界面求和（i=E，W，N，S）。

式（9）的〈1〉項為時間項，可取向后差分：

式中：P為括號內的量取t時刻的值；P^-為括號內的量取t-Δt時刻的值。

當網格不隨時間而變時，不變：

式（9）的〈3〉項為源項，要線性化之后才能求解，一般將S_φ線性表示成：

式（9）的〈2〉項為總通量項，根據φ分布函數的不同取法，一般有中心格式、迎風格式、混合格式、乘方定律格式、指數格式等5種常用差分格式。本文的突擴部分計算采用迎風格式，排取水問題計算采用乘方定律格式。為簡單起見，以迎風格式為例，推導式（8）中的A系數形式。

參考圖3，式（9）的第〈2〉項應為：

對于迎風格式：

w面項為：

［［］］表示取兩者之大。

同理，e面項為：

s面項為：

n面項為：

由上述各式，〈2〉項變為：

式（16）最后一項中，-F_w+F_e-F_s+F_n正是連續方程左邊項的差分應為零，將式（10）、式（11）及式（16）代入式（9）得：

其中

式（17）即為所求的差分方程。

3.邊界條件

AG：V=0或=0，T=0或=0，U、k、ε給定

FH：V=0或=0，=0（φ=U、k、ε、T）

BC：V、k、ε給定，U=0=0

DE：V、k、ε、T給定，U=0

全場水深在平均水深h₀上、下波動，但變幅Δh極小，一般｜Δh｜/h₀&lt；6%。因此，邊界處的水深可以通過附近的計算值外推而得。在排、取水口及上游來流處，根據流量、平均水深及給定的速度分布型式，即可確定U₀或V₀。k₀取為×（0.05U₀）²或×（0.05V₀）²，ε₀取為，式中，c_μ如前所述，取為0.09，L為湍流長度尺度，可取各自的特征長度，如水深、水力半徑等的0.03～0.05倍。

AB，CD，EF，GH：U=V=0，=0，對k、ε應作特殊處理，因這里碰到的是固體壁面問題，不能簡單給定k=0，ε=0，否則將會使全場的k、ε值算小。本文采用常用的“壁函數”處理方法。設距壁面δ處一點P，平行于壁面的速度分量可用下述對數分布公式計算

且假定在壁面附近，湍流處于局部平衡狀態，即認為湍流的對流與擴散可以忽略，產生與耗散相平衡。于是，p點的湍能及湍能耗散率可用下列式子計算：

ε_p=U^*3/（kδ）

式中，k=0.41，E=9.0，U^*=（τ_w/ρ）^1/2。

（三）差分方程的求解方法

如前所述，基本方程組是一相互耦合的非線性方程組。因此，求解工作應包括：①系數線性化，即在求解之前，用各變量的當前值計算出所有系數；②解除各方程之間的耦合，使每個方程都能獨立求解。但是以原始變量U、V、P描述的方程式，U、V與P的耦合不易解除，本文采用h—V（或P—V）修正方法來處理這一困難。即先假定全場的水深分布，使動量方程能夠獨立求解。這時所得的值不是真實的，應通過連續方程來修正，從而給出新的水深分布。如此往復迭代，即能期望得到真實的解。

針對本問題的復雜性，采用線迭代四向掃描的方式進行計算，以提高穩定性，改善收斂性。

（四）程序說明

本文計算程序是在帕坦卡（Patankar）的標準層流程序的基礎上修改而成。修改內容包括：①幾何域、邊界條件、網格劃分；②改壓力p為水深h；③增加湍流量k、ε，溫度T的計算以及它們的壁函數處理；④增加若干子程序。

為方便起見，本程序定名為四向掃描程序。

四、計算成果及其分析

1.突擴通道的計算

為了驗證所用數學模型的正確性及數值方法、計算程序的合理性與可用性，在處理河道排取水問題之前，先對如圖5所示的二維單邊突擴問題進行了計算。

對圖5所示的通道形式，本文僅給出一種工況下的斷面流速無量綱分布曲線（圖6）。由圖6可見，在沿程10個斷面上，除進口附近回流區高度比試驗值小外；其他地方都符合良好，并且比255頁資料的計算結果更符合實際，特別是在上壁面附近流速最大值處更為明顯。

2.排取水問題的計算

本文就寬排取水口問題作了多種工況的計算，這里僅給出圖7作為示例。其中所用符號參照圖4，只是x、y坐標分別應改為無量綱值x/D、y/D。Δh為相對于G點的水深，單位為mm；ΔT₀為排水口DE處的溫升；Re指DE處的雷諾數；i_b為河道縱向底坡；α=Q_上/Q為河道來水量與電廠熱水排放量之比，Q的單位為m³/s，ζ=L/D，ζ₁=B₁/D，ζ₂=B₂/D。

當上游來水比電廠取水為小時（即α&lt；1），排水口處的水流勢必上溯到取水口，回水在全斷面發生（圖7）。因排取水口之間無丁壩，即使α&gt；1，在一定范圍內，仍有回水存在，但回水緊貼排、取水邊岸。受出口處強烈擾動摻混影響，在各種工況下，下游回流始終存在。水深的變化甚微，完全符合剛蓋假定。當有極緩的正縱底坡時，水深沿程略有增加；平坡時，略有減小。當上游回水存在時，取水口的溫度上升；否則，溫度幾乎維持平衡水溫不變。這說明在近區，對流作用是主要的，而擴散影響相對較小。

圖8給出了流線及無量綱k、ε等值線。這是窄排取水口問題，更符合冷卻水工程實際。由于問題本身的復雜性，計算模擬困難更大。現階段本文的計算結果定性與文獻［3］的試驗結果相符：排、取水口上、下游各存在一個回流區，但回流區的大小及位置等，還不能定量符合；對此種特定情形，計算的穩定性與收斂性也有待進一步改善。圖8（b）、（c）分別給出了k、ε的分布。由圖可見，k、ε在排、取水口處最大，說明這里的湍流最為強烈，從而湍能耗散率亦最大，符合物理規律。

五、結束語

k-ε這類精細模型用于近區小尺度問題比較成熟，現已有不少成功的算例>^［2］，本文的突擴通道，亦是典型的近區問題，計算結果與試驗資料符合良好，是對k-ε方法可以成功用于近區的又一次證實；k-ε模型用于全場以建立整區模型，本文亦作過嘗試（見排取水問題的計算），結果表明，這種嘗試是可行的。然而，由于對遠區大尺度問題，對水深較大、水流分層的問題，流散的影響與湍流的影響具有同等意義，有時它們的影響還相互抵償>^［6］。因此，如何處理好流散效應，如何綜合考慮流散與湍流的相互作用，是建立整區模型的關鍵所在，這也正是今后努力的方向。

k-ε方法目前仍處于發展階段，不少問題尚不成熟。對于實際工程中復雜的邊界條件，k、ε的邊界處理這有待進一步改善。另外，k-ε模型增加了k、ε兩個輸運方程，顯然要比積分模型或擴散模型更費機時。因此，目前離冷卻水問題的實際應用尚有一定距離。但作為基本研究，方向無疑是正確的。對于湍流結構，k、ε的變化規律等，應進行試驗研究，才能更好地了解湍流的物理本質，從而推動數學模型的發展。

Numerical Moddelling of the Intake-discharge Problem in Channels with a k-ε Turbulence Model^>[4]

Abstract

參考文獻

>［1］Gerhard H.Jirka，Gerrit Abraham，Donald R.F.Harleman.An Assessment of Techniqes for Hydrothermal Prediction.Technical Report of MIT，1976.

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

>［3］陳惠泉，許玉林.河道冷卻水平面問題的研究∥水利水電科學研究院論文集（第7集）.北京：水利電力出版社，1982.

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

>［5］Rastog，A.K.and Rodi，W..Prediction of heat and mass transfer in open channels.Journal of the Hydraulics Division，ASCE，No.HY3，1978.

>［6］Aspects of Modelling Horizonal Momentum Transfer in Shallow.Report on Investigations，R1150 Delft Hydraulics Laboratory，December，1981.

---