林 玲，梁瑞峰，李克鋒，李 嘉

（四川大學 水力學與山區(qū)河流開發(fā)保護國家重點實驗室，四川 成都 610065）

k-ε雙方程水庫水溫模型浮力模擬分析

林 玲，梁瑞峰，李克鋒，李 嘉

（四川大學 水力學與山區(qū)河流開發(fā)保護國家重點實驗室，四川 成都 610065）

在比較國內外常用水庫水溫紊流模型的基礎上，針對采用k-ε雙方程紊流模型的寬度平均立面二維水庫水溫模型，從理論分析和數值模擬兩方面討論了k-ε紊流模型在模擬溫差異重流時紊動動能生成項和浮力項的作用機理。分析了紊動動能生成項和浮力項對紊動動能和耗散率各自的影響程度，對水溫垂向梯度與紊動動能的相互制約關系進行了探討，研究了紊動動能浮力項和生成項在水庫垂向的分布規(guī)律。W2模型根據實測數據率定的渦黏系數或混合長度在通用性上存在一定不足，但在使用者經驗豐富的條件下仍能得到足夠精確的預測結果。k-ε模型紊動動能生成項和浮力項引起的耗散率改變對紊動動能的影響作為小量可予忽略；k-ε模型可通過紊動動能生成項和浮力項計算與垂向流速和溫度分布相適應的渦黏系數，將溫差浮力影響反映在垂向動量方程中；浮力項主要受垂向溫度梯度控制，紊動動能與垂向溫度梯度的相互制約關系是k-ε紊流模型模擬溫差異重流的關鍵因素。

紊流模型；紊動動能；浮力項；溫差異重流

水庫水溫分層對水庫水質、水生動植物及下游水生環(huán)境都具有重大影響。水溫分層環(huán)境下的水庫一般表層溫度與氣溫變化同步，而具有較穩(wěn)定的庫底水溫、一定的表層與底層溫差和溫躍層，溫躍層范圍內的水溫在垂直方向上存在劇烈變化。這種環(huán)境下的水體流動屬于湍浮力流，表現在垂向上由于溫度差產生浮力，出現分層流動等異于等密度流的特殊流動形式。伴隨這種流動，水庫垂向出現水溫分層形成、發(fā)展和消失的過程。水庫的流場分布與水庫的溫度分布具有強烈的相互作用，因此準確模擬水庫溫度場的關鍵在于精確模擬浮力流流場。

k-ε紊流雙方程模型是在單方程模型的基礎上，增加確定紊動特征長度的微分輸運方程而得到的，因而這種模型仍然采用了各向同性的紊動黏性系數的概念。在雙方程模型中，標準k-ε模型在工程實際中得到了最廣泛的應用和驗證，它考慮了紊動速度比尺和紊動長度比尺的輸運，對于大多數水流問題，標準k-ε雙方程模型或其改進模型一般能得到較好的結果［1-3］。

標準k-ε模型也有其不足，如模型中經驗常數的通用性尚不十分令人滿意，另外，對于有些各向異性明顯的水流或流動區(qū)域，有必要精確地描述應力各分量的輸運，各向同性渦黏性系數的概念和據此建立起來的標準k-ε模型的適用性受到懷疑。有學者嘗試采用各向異性的紊動黏性系數，如Young［4］在其二維水庫模型中對垂向上的動量和熱量的紊動擴散采用了不同經驗公式。

也有學者直接采用各向同性的k-ε模型計算各向異性紊流流動。Farrell等［5］將k-ε模型應用于水庫密度流的模擬，試驗性地計算了一個100 m長的水庫的下潛流過程，結果認為k-ε模型能夠模擬出水庫密度流的下潛、垂向旋渦和溫度分層的特征現象。余真真［6-7］采用k-ε模型計算了溫度分層環(huán)境對潮成內波和懸浮顆粒垂向輸運的影響。鄧云等［8］、劉蘭芬等［9］也對水庫溫差異重流進行了模擬并得到了實測資料的驗證［8-9］。

綜上可知，k-ε模型能夠模擬某些各向異性的紊流流動，如溫差異重流，并得到了實驗資料的良好驗證，但上述文獻并未從模型機理方面深入分析此類模擬成功的理論基礎，故有必要在對比常用紊流模型的基礎上從模型細節(jié)入手找出k-ε模型浮力模擬的主要影響因素。

1 模型對比

雖然有一些紊流模型不采用紊動黏性和紊動擴散的概念而采用紊動動量和紊動熱（質量）通量的微分輸運方程，但鮑辛涅斯克于1877年提出的紊動黏性概念至今仍是模擬雷諾應力的重要基礎。根據紊流模型采用的微分輸運方程的個數，以下從零方程、單方程和雙方程三類模型分別對紊動黏性進行討論。

1.1 零方程模型零方程模型或者據實驗數據建立經驗公式，或者將紊動黏性系數與時均速度建立聯(lián)系，其中最為著名的當為普朗特提出的混合長模型：

式中：?u?y為流速梯度，1/s。

式（1）最難確定的參數即為混合長lm，對紊動量擴散輸運和對流輸運的忽略也使其在模擬復雜流動時易出現謬誤。

CE-QUAL-W2模型是美國陸軍工程兵團于1986年開發(fā)的寬度平均立面二維水質模型。對于湖泊和水庫的溫差異重流模擬，該模型推薦采用W2公式：

式中：vt為垂向渦流系數；κ為卡門常數；lm為混合長度；u為縱向流速；τwy為因風而產生的橫向剪應力；k為波數；τytributary為因支流入流而產生的橫向剪應力； ρ為水體密度；C為常數；Δzmax為垂向網格間距的最大值；Ri為理查森數：

式（2）通過Ri考慮了垂向密度梯度帶來的紊動黏性變化。與普朗特的模型相比，式（2）另一主要改變在于將混合長直接定為計算網格的垂向最大間距，這就使紊動黏性系數vt的確定依賴于網格劃分，實際上成為一個由使用者經驗確定的系數。盡管可以采用實測數據對vt或者lm進行率定，但率定難免將影響水溫的其它因素，如地形概化，轉由vt或者lm承擔，其適用性也存在一定疑問。

此外，CE-QUAL-W2的垂向動量方程式（4）采用了靜水壓力假定，實際是假定了在水庫中垂向對流輸移與擴散處于一種處處平衡的狀態(tài)：

式中：p為壓力；g為重力加速度；α為河底與水平線夾角。

由于密度中考慮了溫度因素，靜水壓力假定并不排除因溫度導致的浮力流速。只考慮兩個方向時，據連續(xù)方程可得垂向流速w：

因而該模型在一定程度上是可以模擬溫差浮力流的。

但在真實的水庫環(huán)境中，處處平衡的對流輸移與擴散畢竟不存在，簡化過程中忽略的對流項將使模型難以模擬有垂向加速度的流場，如表層水體在降溫期或夜晚遇冷下沉、冬季的低溫水入庫下潛等。忽略的擴散項導致缺乏對應的對流項對其進行平衡，通過壓力計算和式（5）、縱向動量方程最終影響到流速的準確計算。

在求解包含完整的垂向動量方程的SIMPLE算法中，需要通過迭代解決動量方程確定的流場與連續(xù)方程確定的流場之間的不匹配問題，是較為耗時的。而CE-QUAL-W2通過靜水壓力假定和連續(xù)方程直接求解垂向流速，不需要解決流場在兩個方程中的差異，因而其計算速度是比SIMPLE算法有較大優(yōu)勢的，且在使用者經驗豐富的條件下仍能得出足夠精確的預測結果。

1.2 單方程模型單方程模型的速度比尺 k由微分方程確定，柯莫哥洛夫和普朗特分別獨立得到由速度和長度比尺確定的紊動黏性系數公式：

式中：C′μ為經驗常數；k為紊動動能；L為長度比尺。

紊動動能的精確輸運方程為：

式中：ve為分子黏性系數v與紊動渦黏系數vt之和；σk為紊動普朗特數；Gk和Gb為紊動動能生成項和浮力項，β為等壓膨脹系數；σT為紊動普朗特數。

封閉方程組的黏性耗散項通過因次分析得到：

式中：ε為紊動動能耗散率，CD亦為經驗常數。

k方程考慮了浮力項Gb的影響，但長度比尺L如同混合長模型中的混合長度一樣難以采用經驗方法確定，限制了單方程模型的使用，其所能準確模擬的剪力層流動采用混合長模型也可得到滿意的結果且計算更為簡單。

1.3 雙方程模型對表征漩渦尺寸的長度比尺紊動動能耗散率ε采用微分方程計算，單方程模型就轉為雙方程模型：

式中：σε為紊動動能耗散普朗特數；Cε1、Cε2、Cε3為模型常數。

雙方程模型解決了單方程模型長度比尺的確定難題，同時，消去式（6）和式（8）中的長度比尺可得到：

式中：Cμ=C′μCD。Cμ避免了對和的各自單獨確定，Cμ和Cε1、Cε2均可通過剪力層、近壁區(qū)等特殊區(qū)域的試驗測定得到較為準確的數值。

但雙方程在使用中如果采用完全的不可壓縮假定，也將導致浮力模擬的偏差。如美國ANSYS公司的FLUENT軟件，所采用的渦黏模型包括標準k-ε模型、RNGk-ε模型和Realizablek-ε模型，由于不可壓流體假定使等壓膨脹系數為零，垂向溫差引起的密度流不能通過浮力項G來影響b渦黏系數，這些模型在計算溫差異重流時是不能反映垂向密度差引起的熱對流和紊動擴散的。

2 采用k-ε雙方程紊流模型的寬度平均立面二維水庫水溫模型

雖然紊流問題均是三維問題，但在目前的技術經濟條件下對大體積水體進行三維模擬在時間和必要性上尚存在較大改進空間，本文擬采用寬度平均的k-ε紊流模型和立面二維水庫水溫控制方程對影響浮力生成的紊動動能浮力生成項進行分析，闡明其模擬溫差異重流的關鍵所在。

出于模型概化的考慮，在密度變化不大的浮力流問題中，Boussinesq認為可以近似地只在重力項中考慮浮力的影響，而在控制方程的其它項中忽略浮力的作用。當質量力僅為重力時，Navier-Stokes方程中的體積力可表達為：

其中：ΔT=T-Ta，Ta為與參考密度 ρa相對應的參考溫度；β為等壓膨脹系數。將Boussinesq近似式（11）引入瞬時Navier-Stokes控制方程并通過k-ε雙方程?；?，在河寬方向作積分，可得到湍浮力流的寬度平均方程。式（12）—（16）分別為縱向動量、垂向動量、紊動動能和紊動動能耗散率以及溫度方程：

由于水溫在4℃時密度最大，在4℃左右的浮力作用規(guī)律是不一致的，本文只討論水溫大于4℃的情形，通過對k-ε紊流模型浮力項影響的分析確定其對溫差異重流浮力模擬的具體影響因素。

3 原型水庫

在對模型理論分析基礎上，需要對水庫進行模擬以直觀說明各因素的實際作用和影響量級。

選取岷江上游的紫坪鋪水庫作為模擬水庫。水庫總庫容為11.12×108m3，調節(jié)庫容7.74×108m3，具有不完全年調節(jié)能力，正常蓄水位877 m，死水位817 m?，F場觀測表明該水庫水溫在升溫期某些月份存在垂向分層現象，可用于驗證模型的浮力計算效果。

收集到2009年11月6日、2010年7月6日紫坪鋪水庫全庫區(qū)的水溫觀測以及同步的入庫水溫、氣象、調度等資料，模型計算以2009年11月6日測量的全庫區(qū)水溫分布作為計算初始水溫；流場計算據庫尾入流水溫進行多年循環(huán)計算，待庫區(qū)水溫穩(wěn)定收斂后以11月6日的流場分布作為初始流場。圖1為據實測庫區(qū)水溫作為初始條件計算得到的2010年5月15日的庫區(qū)水溫，溫躍層在高程790 m左右。計算結果得到了紫坪鋪水庫實測數據的驗證，成果將另文述評。

4 浮力模擬影響因素分析

熱量在水體垂向上的傳遞主要依賴于隨流輸移和紊動擴散。在采用k-ε雙方程的水庫水溫模型中，紊動動能k和耗散率通過渦黏系數影響熱擴散系數和垂向流速w。

圖1 紫坪鋪水庫2010年5月15日的庫區(qū)水溫分布（計算值）

4.1Gb和對Gk的影響紊動動能浮力項Gb使模型考慮了溫度垂向分布對熱擴散和紊動擴散的作用。為考察Gb對紊動動能k的影響，可據式（14）寫出紊動動能的簡化離散方程：

其中Lk表示本地輸移和擴散對k的影響。

水溫分層是通過浮力生成項Gb對紊動動能產生作用的。從式（17）可知，當T＞4℃，溫度在垂向上存在正梯度，即穩(wěn)定分層時，Gb為負值，直接抑制下一時刻紊動動能k的增長。

圖2顯示了水庫垂向的紊動動能變化與溫度梯度變化對比關系。水庫表層溫度梯度主要受太陽輻射和長波輻射等熱通量控制，與同高程處的紊動動能無明顯的相關關系。在發(fā)電引水孔口之上（高程820 m左右）的主流區(qū)域邊緣，流速的垂向梯度因處于主流區(qū)域邊緣而變化劇烈（見圖3），使紊動生成項Gk增加從而帶來較強的紊動，紊動擴散使區(qū)域內趨于同溫，對應該高程處較低的溫度梯度。而在靠近庫底的780 m高程附近，各向流速和流速梯度均較小，缺乏紊動動能生成來源，垂向熱量難以通過隨流輸移和紊動擴散傳遞，形成了局部較高的溫度梯度。

圖2 紫坪鋪水庫2010年5月15日壩前200 m的?T?z和k值

4.2 Gb、Gk對ε的影響紊動動能耗散率的簡化離散方程為：

其中Lε表示本地輸移和擴散對ε的影響。

紊動動能耗散也直接影響下一時刻的紊動動能（式（17），右側-ε項），而溫度垂向正梯度對kt、εt的作用是相反的。但由于Cε1和Cε2均小于1，通過和Cε3對Gb、Gk的“稀釋”（式（18）），Gb、 Gk對ε的影響極為有限。圖4給出了岷江紫坪鋪水庫2010年5月15日壩前200 m處的垂向值，可見溫度分布對εt的作用與對kt的作用相比為小量，可以忽略。

圖3 紫坪鋪水庫2010年5月15日壩前區(qū)域的流場

圖4 紫坪鋪水庫2010年5月15日壩前200m的值

4.3 k、ε對渦黏系數vt的影響k、ε對隨流輸移和紊動擴散的影響通過渦黏系數vt來傳遞。

k-ε雙方程計算得到的vt是各向同性的，但時均流通常具有傾向性的主流向，這種傾向性對大比尺紊動施加影響，使大比尺紊動具有很強的各向異性，紊動強度和紊動長度比尺因方向而異。如果雷諾數足夠大，大比尺運動和小比尺運動在譜域上的間距足夠大，方向靈敏性便完全消失，使得耗散能量的小比尺紊動變?yōu)楦飨蛲裕?0］。在水庫環(huán)境下，雷諾數變化范圍較大，在流速較緩區(qū)域的流動具有各向異性特性，浮力對垂向對流的影響不可忽略，因而需要計算由浮力引起的vt變化。

由于一般情形下縱向流速總是遠大于垂向流速，縱向動量方程中的慣性項主要由縱向壓力梯度項來平衡，但垂向壓力梯度與縱向壓力梯度不在一個量級，因此垂向慣性項便主要由紊動擴散項和重力項平衡。此時垂向隨流輸移中的紊動輸運項和熱量方程中的熱擴散項就對慣性項起著重要的平衡作用。另一方面，溫度控制方程中不含有壓力梯度項和浮力項，除非在源項很大的情形下，紊動輸運項總是與慣性相平衡的重要因素，而除了表面熱源外，水面下除了少量短波輻射進入外，一般不存在外來熱源，因此在重力作用下，主流區(qū)域上邊緣（圖2，820 m附近）由于存在一定的垂向流速，對應地應存在較大的紊動擴散系數與慣性項平衡。

圖5為紫坪鋪壩前分別采用W2模型和k-ε模型計算的渦黏系數和水溫。k-ε模型得到的渦黏系數分布與圖2的流場是相匹配的，主流區(qū)域上邊緣的大擴散系數也帶來820 m高程上下的均溫分布。但W2模型得到的渦黏系數由于主要受流速梯度和混合長度影響呈不規(guī)則分布且其量級遠大于k-ε模型計算值，說明渦黏系數承擔了本應由對流和垂向紊動擴散承擔的輸運部分。計算中，縱向流速u的上邊界采用零梯度邊界，這對W2模型和k-ε模型的影響是一致的，都導致了較小的渦黏系數，但W2模型由于缺乏垂向流速對強紊動的平衡，使表層熱量難以進行對流輸移，因而其表層水溫高于k-ε模型計算值。

圖5 紫坪鋪水庫2010年5月15日壩前200 m的渦黏系數和水溫

5 結論

本文對國內外常用水庫水溫模型采用的紊流模型進行了對比，對其浮力模擬能力進行了分析，并以采用k-ε雙方程紊流模型的寬度平均立面二維水庫水溫模型為基礎，從理論分析和數值模擬兩方面討論了的k-ε紊流模型在模擬溫差異重流時紊動動能浮力項和生成項的作用機理。研究表明：

（1）零方程的W2模型具有浮力模擬的理論基礎，其渦黏系數的計算依賴于經驗確定的混合長度，據實測數據率定的渦黏系數或混合長度將帶入不確定的其它影響，使參數缺乏通用性，但在使用者經驗豐富的條件下仍能得到足夠精確的預測結果。

（2）雙方程k-ε模型的紊動動能浮力項和生成項對紊動動能和耗散率的影響程度具有較大差異，紊動動能浮力項和生成項引起的耗散率改變對紊動動能的影響作為小量可予忽略。

（3）雙方程k-ε模型可以通過紊動動能生成項和浮力項計算與垂向流速和溫度分布相適應的渦黏系數，將溫差浮力影響反映在垂向動量方程中；浮力項Gb主要受垂向溫度梯度控制，紊動動能與垂向溫度梯度的相互制約關系是k-ε紊流模型模擬溫差異重流的關鍵因素。

［1］ 羅麟，趙文謙，李嘉.三維湍浮力回流數學模型在環(huán)境水力學中的應用［J］.水利學報，1992（9）：1-7.

［2］ 胡振紅，沈永明，鄭永紅，等.溫度和鹽度分層流的數值模擬［J］.水科學進展，2001，12（4）：439-444.

［3］ BOURNET P E，DARTUS D，TASSIN B，et al.Numerical investigation of plunging density current［J］.J.Hy?dr.Engng.，1999，125（6）：584-594.

［4］ YOUNG D L.CHAPTER 10：Finite element analysis of stratified lake hydrodynamics［R］//Environmental Fluid mechanics—theories and applications.Edited by Hayley Shen，et al.2002.

［5］ GERARD J FARRELL，HEINZ G STEFAN.Mathematical modeling of plunging reservoir［J］.J.Hydraulic Re?search，1988，26（5）：525-537.

［6］ 余真真，王玲玲，朱海，等.三峽水庫香溪河庫灣潮成內波數值模擬［J］.四川大學學報：工程科學版，2012，44（3）：26-30.

［7］ 余真真，王玲玲，戴氣超，等.水溫分層對水體中懸浮顆粒物垂向輸運影響的研究［J］.四川大學學報：工程科學版，2011，43（1）：65-70.

［8］ 鄧云，李嘉，羅麟，等.河道型深水庫的溫度分層模擬［J］.水動力學研究與進展，2004，19（5）：604-609.

［9］ 劉蘭芬，張士杰，劉暢，等.漫灣水電站水庫水溫分布觀測與數學模型計算研究［J］.中國水利水電科學研究院學報，2007，5（2）：87-94.

［10］ 金忠青.N-S方程的數值解和紊流模型［M］.南京：河海大學出版社，1989.

The buoyance simulation analysis of k-ε double equation reservoir water temperature model

LIN Ling，LIANG Ruifeng，LI Kefeng，LI Jia
（State Key Laboratory of Hydraulics and Mountain River Engineering，Chengdu 610065，China）

Based on a comparative study of the reservoir water temperature turbulence models commonly used，a reservoir water temperature model with averaged width and 2-D elevation adopting k-ε double equa?tion turbulence model is used to simulate thermal density flow.In the view of theoretical analysis and nu?merical simulation，this paper discussed the mechanism of turbulent kinetic energy’s buoyancy term and pro?duction term.It investigated the buoyancy term’s and production term’s impacts on turbulent kinetic energy and energy dissipation，respectively.This paper explored the inter restricted relationship between water tem?perature vertical gradient and turbulent kinetic energy，and studied the vertical distribution law of turbulent kinetic energy’s buoyancy term and production term in reservoir.The eddy viscosity coefficient or mixing length in W2 model was lack of generality，which was rated by measured data，but a sufficient precise pre?diction result can be obtained by a veteran.The effect of the energy dissipation caused by the buoyancy term and production term on turbulent kinetic energy can be ignored as a small factor in k-ε turbulence model.This model can calculate eddy viscosity corresponding to the vertical velocity and temperature distri?bution by turbulent kinetic energy’s buoyancy term and production term.Therefore，the influence of thermal buoyancy was reflected in vertical momentum equation.The buoyancy term was mainly controlled by temper?ature vertical gradient.The inter restricted relationship of water temperature vertical gradient and turbulent kinetic energy was the key factor in a simulation of thermal density flow with the k-ε turbulence model.

turbulent model；turbulent kinetic energy；buoyancy term；thermal density flow

TV131.3

A

10.13244/j.cnki.jiwhr.2017.01.005

1672-3031（2017）01-0037-07

（責任編輯：李福田）

2016-04-08

國家自然科學基金項目（51279114）

林玲（1992-），女，四川瀘州人，碩士生，主要從事環(huán)境水力學研究。E-mail：enidlinling@outlook.com

梁瑞峰（1974-），男，河南南陽人，博士，副研究員，主要從事環(huán)境水力學研究。E-mail：liangruifeng@scu.edu.cn