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

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

（四川大學(xué) 水力學(xué)與山區(qū)河流開發(fā)保護(hù)國家重點(diǎn)實(shí)驗(yàn)室，四川 成都 610065）

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

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

（四川大學(xué) 水力學(xué)與山區(qū)河流開發(fā)保護(hù)國家重點(diǎn)實(shí)驗(yàn)室，四川 成都 610065）

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

紊流模型；紊動(dòng)動(dòng)能；浮力項(xiàng)；溫差異重流

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

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

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

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

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

1 模型對(duì)比

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

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

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

式（1）最難確定的參數(shù)即為混合長lm，對(duì)紊動(dòng)量擴(kuò)散輸運(yùn)和對(duì)流輸運(yùn)的忽略也使其在模擬復(fù)雜流動(dòng)時(shí)易出現(xiàn)謬誤。

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

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

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

此外，CE-QUAL-W2的垂向動(dòng)量方程式（4）采用了靜水壓力假定，實(shí)際是假定了在水庫中垂向?qū)α鬏斠婆c擴(kuò)散處于一種處處平衡的狀態(tài)：

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

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

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

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

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

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

式中：C′μ為經(jīng)驗(yàn)常數(shù)；k為紊動(dòng)動(dòng)能；L為長度比尺。

紊動(dòng)動(dòng)能的精確輸運(yùn)方程為：

式中：ve為分子黏性系數(shù)v與紊動(dòng)渦黏系數(shù)vt之和；σk為紊動(dòng)普朗特?cái)?shù)；Gk和Gb為紊動(dòng)動(dòng)能生成項(xiàng)和浮力項(xiàng)，β為等壓膨脹系數(shù)；σT為紊動(dòng)普朗特?cái)?shù)。

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

式中：ε為紊動(dòng)動(dòng)能耗散率，CD亦為經(jīng)驗(yàn)常數(shù)。

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

1.3 雙方程模型對(duì)表征漩渦尺寸的長度比尺紊動(dòng)動(dòng)能耗散率ε采用微分方程計(jì)算，單方程模型就轉(zhuǎn)為雙方程模型：

式中：σε為紊動(dòng)動(dòng)能耗散普朗特?cái)?shù)；Cε1、Cε2、Cε3為模型常數(shù)。

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

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

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

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

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

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

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

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

3 原型水庫

在對(duì)模型理論分析基礎(chǔ)上，需要對(duì)水庫進(jìn)行模擬以直觀說明各因素的實(shí)際作用和影響量級(jí)。

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

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

4 浮力模擬影響因素分析

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

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

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

其中Lk表示本地輸移和擴(kuò)散對(duì)k的影響。

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

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

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

4.2 Gb、Gk對(duì)ε的影響紊動(dòng)動(dòng)能耗散率的簡化離散方程為：

其中Lε表示本地輸移和擴(kuò)散對(duì)ε的影響。

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

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

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

4.3 k、ε對(duì)渦黏系數(shù)vt的影響k、ε對(duì)隨流輸移和紊動(dòng)擴(kuò)散的影響通過渦黏系數(shù)vt來傳遞。

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

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

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

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

5 結(jié)論

本文對(duì)國內(nèi)外常用水庫水溫模型采用的紊流模型進(jìn)行了對(duì)比，對(duì)其浮力模擬能力進(jìn)行了分析，并以采用k-ε雙方程紊流模型的寬度平均立面二維水庫水溫模型為基礎(chǔ)，從理論分析和數(shù)值模擬兩方面討論了的k-ε紊流模型在模擬溫差異重流時(shí)紊動(dòng)動(dòng)能浮力項(xiàng)和生成項(xiàng)的作用機(jī)理。研究表明：

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

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

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

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

［2］ 胡振紅，沈永明，鄭永紅，等.溫度和鹽度分層流的數(shù)值模擬［J］.水科學(xué)進(jìn)展，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］ 余真真，王玲玲，朱海，等.三峽水庫香溪河庫灣潮成內(nèi)波數(shù)值模擬［J］.四川大學(xué)學(xué)報(bào)：工程科學(xué)版，2012，44（3）：26-30.

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

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

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

［10］ 金忠青.N-S方程的數(shù)值解和紊流模型［M］.南京：河海大學(xué)出版社，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

（責(zé)任編輯：李福田）

2016-04-08

國家自然科學(xué)基金項(xiàng)目（51279114）

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

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