• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

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

    2017-05-12 07:14:26梁瑞峰李克鋒
    關(guān)鍵詞:異重流浮力溫差

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

    (四川大學(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

    猜你喜歡
    異重流浮力溫差
    手表與溫差
    小浪底水庫異重流排沙效率分析
    “浮力”知識(shí)鞏固
    我們一起來“制服”浮力
    浮力大小由誰定
    水生植被影響異重流動(dòng)力特性的試驗(yàn)分析
    北方冬季養(yǎng)豬當(dāng)心“四差”
    溫差“催甜”等
    低級(jí)熱溫差發(fā)電器的開發(fā)應(yīng)用研究
    電子制作(2018年23期)2018-12-26 01:01:20
    改進(jìn)的徑向基神經(jīng)網(wǎng)絡(luò)模型在水庫異重流泥沙淤積量模擬中的應(yīng)用
    中国美女看黄片| 国产成人aa在线观看| 久99久视频精品免费| 一进一出抽搐动态| 啦啦啦观看免费观看视频高清| 精品福利观看| 香蕉av资源在线| 99国产精品一区二区蜜桃av| 国产精品久久久久久久久免| 欧美日韩国产亚洲二区| 村上凉子中文字幕在线| 欧美性感艳星| 欧美色欧美亚洲另类二区| 久久久久久久亚洲中文字幕| 色哟哟哟哟哟哟| 精品一区二区免费观看| 国内毛片毛片毛片毛片毛片| 亚洲精品久久国产高清桃花| aaaaa片日本免费| 中文字幕免费在线视频6| 国产精品一区二区免费欧美| 搡老熟女国产l中国老女人| 亚洲最大成人中文| 欧美日韩综合久久久久久 | 网址你懂的国产日韩在线| 一级黄色大片毛片| 午夜福利在线观看免费完整高清在 | 亚洲乱码一区二区免费版| 久久精品夜夜夜夜夜久久蜜豆| 国产精品人妻久久久影院| 最近视频中文字幕2019在线8| 日韩av在线大香蕉| 别揉我奶头~嗯~啊~动态视频| 99国产极品粉嫩在线观看| 亚洲乱码一区二区免费版| 久久久成人免费电影| 亚洲无线在线观看| 日韩欧美精品免费久久| 麻豆av噜噜一区二区三区| 中出人妻视频一区二区| 精品午夜福利视频在线观看一区| 国产高清三级在线| 免费av观看视频| 精品福利观看| 99久久久亚洲精品蜜臀av| 亚洲国产精品sss在线观看| 丰满人妻一区二区三区视频av| 综合色av麻豆| а√天堂www在线а√下载| 春色校园在线视频观看| 久久精品国产亚洲网站| 亚洲久久久久久中文字幕| 免费在线观看日本一区| 国产精品亚洲美女久久久| 欧美三级亚洲精品| 午夜免费成人在线视频| 日本色播在线视频| 免费电影在线观看免费观看| 我要看日韩黄色一级片| 国产精品久久视频播放| 亚洲人成伊人成综合网2020| 久久久久久久久久久丰满 | 老熟妇仑乱视频hdxx| 成人国产一区最新在线观看| 午夜福利18| 欧美黑人欧美精品刺激| 美女黄网站色视频| 中文字幕熟女人妻在线| 中国美女看黄片| 久久久久性生活片| av黄色大香蕉| 国产高清不卡午夜福利| 男女做爰动态图高潮gif福利片| 亚洲黑人精品在线| 老师上课跳d突然被开到最大视频| 国产视频一区二区在线看| 我的女老师完整版在线观看| 国产成人av教育| 色综合婷婷激情| 久久欧美精品欧美久久欧美| 亚洲国产精品成人综合色| 内地一区二区视频在线| 欧美精品啪啪一区二区三区| 精品一区二区免费观看| 成人av在线播放网站| 亚洲av中文字字幕乱码综合| 日本色播在线视频| 成人无遮挡网站| 中文字幕人妻熟人妻熟丝袜美| 亚洲人成网站在线播放欧美日韩| 亚洲一级一片aⅴ在线观看| 国产亚洲av嫩草精品影院| 啦啦啦观看免费观看视频高清| 国产精品综合久久久久久久免费| 我要搜黄色片| 日韩欧美在线二视频| 波野结衣二区三区在线| 国产成人福利小说| 国产69精品久久久久777片| 熟妇人妻久久中文字幕3abv| 国产精品一及| 午夜精品久久久久久毛片777| 欧美潮喷喷水| 亚洲成人中文字幕在线播放| 在线观看午夜福利视频| 又黄又爽又刺激的免费视频.| 午夜福利在线观看免费完整高清在 | 亚洲精品一卡2卡三卡4卡5卡| 美女黄网站色视频| 日本撒尿小便嘘嘘汇集6| 亚洲,欧美,日韩| 美女免费视频网站| 毛片女人毛片| 国产精品亚洲美女久久久| 婷婷色综合大香蕉| 亚洲va在线va天堂va国产| 联通29元200g的流量卡| 亚洲精品日韩av片在线观看| 亚洲国产精品成人综合色| 又爽又黄a免费视频| 999久久久精品免费观看国产| 国产成人福利小说| 亚洲av电影不卡..在线观看| 人妻夜夜爽99麻豆av| 哪里可以看免费的av片| 欧美成人免费av一区二区三区| 亚洲专区国产一区二区| 亚洲无线在线观看| 高清在线国产一区| 人人妻,人人澡人人爽秒播| 国产黄色小视频在线观看| 亚洲avbb在线观看| 国内少妇人妻偷人精品xxx网站| 最近视频中文字幕2019在线8| 91麻豆av在线| 悠悠久久av| 欧美性感艳星| 午夜久久久久精精品| 身体一侧抽搐| 一级毛片久久久久久久久女| 国产 一区 欧美 日韩| 91久久精品电影网| 少妇猛男粗大的猛烈进出视频 | 亚洲男人的天堂狠狠| 黄色日韩在线| 琪琪午夜伦伦电影理论片6080| 夜夜夜夜夜久久久久| 白带黄色成豆腐渣| 国产精品免费一区二区三区在线| 成人午夜高清在线视频| 老司机午夜福利在线观看视频| 久久久久免费精品人妻一区二区| 欧美三级亚洲精品| 搡老岳熟女国产| 亚洲最大成人av| 别揉我奶头~嗯~啊~动态视频| 老师上课跳d突然被开到最大视频| 国语自产精品视频在线第100页| 国产午夜精品久久久久久一区二区三区 | 赤兔流量卡办理| 麻豆久久精品国产亚洲av| 国产精品免费一区二区三区在线| videossex国产| 欧洲精品卡2卡3卡4卡5卡区| 亚洲成人久久性| 听说在线观看完整版免费高清| 校园春色视频在线观看| 白带黄色成豆腐渣| 亚洲真实伦在线观看| 最近视频中文字幕2019在线8| 国产精品一区www在线观看 | 日本色播在线视频| 久久久久久大精品| 性插视频无遮挡在线免费观看| 免费在线观看成人毛片| 免费观看精品视频网站| 国产精品久久久久久久电影| avwww免费| 校园人妻丝袜中文字幕| 很黄的视频免费| 色综合婷婷激情| 三级国产精品欧美在线观看| 亚洲综合色惰| 国模一区二区三区四区视频| 国产日本99.免费观看| 可以在线观看毛片的网站| 国产免费一级a男人的天堂| 国产成年人精品一区二区| 桃红色精品国产亚洲av| 国产成人a区在线观看| 麻豆久久精品国产亚洲av| 亚洲国产欧美人成| 日韩亚洲欧美综合| 人妻丰满熟妇av一区二区三区| av福利片在线观看| 国产伦精品一区二区三区视频9| 深爱激情五月婷婷| 97超级碰碰碰精品色视频在线观看| 人妻丰满熟妇av一区二区三区| 日韩中文字幕欧美一区二区| 99热网站在线观看| 久久久久九九精品影院| 日韩欧美国产在线观看| 又紧又爽又黄一区二区| 国产精品免费一区二区三区在线| 久久午夜福利片| 欧美日韩亚洲国产一区二区在线观看| 亚洲无线观看免费| 亚洲乱码一区二区免费版| 欧美性猛交╳xxx乱大交人| 很黄的视频免费| 精品国内亚洲2022精品成人| 一区二区三区高清视频在线| 国产一区二区激情短视频| 免费观看的影片在线观看| 真实男女啪啪啪动态图| 久久香蕉精品热| 久久久久久伊人网av| 一区二区三区免费毛片| 99热这里只有是精品50| 两个人的视频大全免费| 国产一区二区在线av高清观看| 久久午夜福利片| 精品一区二区三区av网在线观看| 亚洲国产高清在线一区二区三| 午夜福利欧美成人| 一夜夜www| 亚洲人成网站在线播| av在线亚洲专区| 97碰自拍视频| 深爱激情五月婷婷| 99久国产av精品| 国产高清不卡午夜福利| 国产黄色小视频在线观看| 成人综合一区亚洲| 国产91精品成人一区二区三区| 男人和女人高潮做爰伦理| 亚洲午夜理论影院| 国产精品无大码| av视频在线观看入口| 精品99又大又爽又粗少妇毛片 | 亚洲精品色激情综合| 国产亚洲欧美98| 国产一级毛片七仙女欲春2| 岛国在线免费视频观看| 免费av观看视频| 亚洲av美国av| 亚洲精品乱码久久久v下载方式| 亚洲精品亚洲一区二区| 亚洲av免费高清在线观看| 少妇熟女aⅴ在线视频| 少妇的逼好多水| 直男gayav资源| 人人妻人人澡欧美一区二区| 18禁黄网站禁片午夜丰满| 精品不卡国产一区二区三区| 嫩草影院精品99| 午夜久久久久精精品| 国产爱豆传媒在线观看| 欧美性猛交黑人性爽| 最新在线观看一区二区三区| 尾随美女入室| 久久精品国产清高在天天线| 亚洲美女搞黄在线观看 | 国产一区二区亚洲精品在线观看| 成人美女网站在线观看视频| 我的女老师完整版在线观看| 久久久久久久久久久丰满 | 51国产日韩欧美| 九九在线视频观看精品| 久久午夜福利片| 国产伦精品一区二区三区四那| 亚洲国产精品sss在线观看| 日本免费a在线| 直男gayav资源| 99久久久亚洲精品蜜臀av| 亚洲狠狠婷婷综合久久图片| 国内精品久久久久久久电影| 看免费成人av毛片| 久久久精品欧美日韩精品| 99热这里只有是精品50| 久久中文看片网| 成熟少妇高潮喷水视频| 18禁黄网站禁片免费观看直播| 亚洲人与动物交配视频| 我的老师免费观看完整版| 最近最新免费中文字幕在线| 欧美激情久久久久久爽电影| 嫩草影院新地址| 午夜视频国产福利| 精品久久国产蜜桃| 精品国内亚洲2022精品成人| 97超视频在线观看视频| 观看免费一级毛片| 亚洲精品一区av在线观看| 观看美女的网站| 久久国产精品人妻蜜桃| 此物有八面人人有两片| 国产精品久久久久久久电影| 少妇被粗大猛烈的视频| 禁无遮挡网站| 免费一级毛片在线播放高清视频| 国产精品久久久久久亚洲av鲁大| 狂野欧美激情性xxxx在线观看| 成人二区视频| 久久热精品热| 亚洲av免费高清在线观看| 波野结衣二区三区在线| 日本a在线网址| 中文资源天堂在线| 99九九线精品视频在线观看视频| 人人妻人人看人人澡| 亚洲18禁久久av| 尤物成人国产欧美一区二区三区| 午夜a级毛片| 国产人妻一区二区三区在| 少妇熟女aⅴ在线视频| 国产精品一及| 一进一出抽搐动态| 一个人看的www免费观看视频| 在线免费观看的www视频| 国内精品久久久久久久电影| 又黄又爽又免费观看的视频| 成人国产麻豆网| 亚洲国产色片| 一进一出抽搐gif免费好疼| 男人的好看免费观看在线视频| 免费无遮挡裸体视频| 琪琪午夜伦伦电影理论片6080| 九九爱精品视频在线观看| 嫩草影视91久久| 亚洲aⅴ乱码一区二区在线播放| 午夜精品一区二区三区免费看| 人人妻,人人澡人人爽秒播| 69人妻影院| 中文字幕人妻熟人妻熟丝袜美| 嫩草影视91久久| 国产淫片久久久久久久久| 永久网站在线| 午夜老司机福利剧场| 午夜免费男女啪啪视频观看 | 婷婷色综合大香蕉| 久久午夜亚洲精品久久| av在线亚洲专区| 99久久中文字幕三级久久日本| 亚洲精品456在线播放app | 国产精品嫩草影院av在线观看 | 国产一区二区在线观看日韩| 午夜福利在线观看吧| 欧美高清成人免费视频www| 男人舔女人下体高潮全视频| 亚洲欧美日韩高清专用| 国产女主播在线喷水免费视频网站 | 亚洲av二区三区四区| 草草在线视频免费看| 久久午夜亚洲精品久久| 国产伦在线观看视频一区| 国产伦精品一区二区三区视频9| 三级毛片av免费| 欧美最新免费一区二区三区| 久久久久久大精品| 国产伦精品一区二区三区四那| 哪里可以看免费的av片| 直男gayav资源| 亚洲美女黄片视频| 精品久久久久久久末码| 韩国av一区二区三区四区| 免费不卡的大黄色大毛片视频在线观看 | 99热这里只有是精品50| 久久久久久久久久黄片| 日本一本二区三区精品| 色吧在线观看| 桃色一区二区三区在线观看| 国产爱豆传媒在线观看| 最近最新免费中文字幕在线| 国产亚洲精品综合一区在线观看| 亚洲av免费在线观看| 在线免费观看不下载黄p国产 | 亚洲国产精品成人综合色| 看片在线看免费视频| 国国产精品蜜臀av免费| 欧美潮喷喷水| 最近视频中文字幕2019在线8| 在线观看舔阴道视频| 偷拍熟女少妇极品色| 伦精品一区二区三区| 国产伦人伦偷精品视频| 少妇高潮的动态图| 性色avwww在线观看| 夜夜看夜夜爽夜夜摸| 免费观看的影片在线观看| 免费无遮挡裸体视频| 99九九线精品视频在线观看视频| 真实男女啪啪啪动态图| 最新中文字幕久久久久| 亚洲av不卡在线观看| 少妇被粗大猛烈的视频| 又爽又黄a免费视频| 国产精品国产三级国产av玫瑰| 十八禁国产超污无遮挡网站| АⅤ资源中文在线天堂| 亚洲人成网站高清观看| www.色视频.com| 91精品国产九色| 97碰自拍视频| 免费观看的影片在线观看| 国语自产精品视频在线第100页| 欧美性猛交╳xxx乱大交人| 99热精品在线国产| 一区福利在线观看| 国产精品伦人一区二区| 动漫黄色视频在线观看| 中亚洲国语对白在线视频| 欧美+日韩+精品| 亚洲av熟女| .国产精品久久| 春色校园在线视频观看| 在线观看av片永久免费下载| 干丝袜人妻中文字幕| 少妇被粗大猛烈的视频| 观看免费一级毛片| 国产精品三级大全| 欧美激情国产日韩精品一区| 亚洲五月天丁香| 欧美区成人在线视频| 久久午夜亚洲精品久久| 日韩欧美国产一区二区入口| 99久国产av精品| 黄色女人牲交| 精品国内亚洲2022精品成人| 极品教师在线免费播放| 精品福利观看| 日韩亚洲欧美综合| 亚洲真实伦在线观看| 精品午夜福利视频在线观看一区| 91麻豆精品激情在线观看国产| 精品久久久噜噜| 欧美性猛交黑人性爽| 欧美+日韩+精品| 天天躁日日操中文字幕| 国产一区二区三区视频了| 欧美xxxx性猛交bbbb| 色在线成人网| 国产精品永久免费网站| 少妇的逼好多水| 国产亚洲精品av在线| av视频在线观看入口| 在线观看免费视频日本深夜| 两个人视频免费观看高清| 特级一级黄色大片| 欧美性感艳星| 亚洲va日本ⅴa欧美va伊人久久| 国产精品久久久久久av不卡| av在线蜜桃| 国产日本99.免费观看| 少妇被粗大猛烈的视频| 老司机午夜福利在线观看视频| 国产精品福利在线免费观看| aaaaa片日本免费| 午夜福利高清视频| 午夜福利在线观看吧| 国产一级毛片七仙女欲春2| 美女 人体艺术 gogo| 精华霜和精华液先用哪个| 又爽又黄a免费视频| 久久人人爽人人爽人人片va| 91狼人影院| 日韩av在线大香蕉| av在线蜜桃| 成年免费大片在线观看| 他把我摸到了高潮在线观看| 真人做人爱边吃奶动态| 久久久久久久精品吃奶| 人妻夜夜爽99麻豆av| 日本黄色片子视频| 亚洲综合色惰| 看黄色毛片网站| 国产高清有码在线观看视频| 国产精品亚洲美女久久久| 日韩欧美免费精品| 国产亚洲精品综合一区在线观看| 国产国拍精品亚洲av在线观看| 大又大粗又爽又黄少妇毛片口| 亚洲av成人av| 欧美一级a爱片免费观看看| 亚洲av免费在线观看| av黄色大香蕉| 免费av毛片视频| av在线老鸭窝| 深爱激情五月婷婷| 欧美区成人在线视频| 亚洲人与动物交配视频| 国内揄拍国产精品人妻在线| 成人欧美大片| 69av精品久久久久久| 麻豆成人av在线观看| 俺也久久电影网| 国内精品美女久久久久久| 一本久久中文字幕| 日本-黄色视频高清免费观看| 午夜福利视频1000在线观看| 日韩中文字幕欧美一区二区| 日本一二三区视频观看| 又紧又爽又黄一区二区| 性插视频无遮挡在线免费观看| 小说图片视频综合网站| 精品无人区乱码1区二区| 国产高清激情床上av| 国产欧美日韩精品一区二区| 欧美成人性av电影在线观看| 国产午夜福利久久久久久| 亚洲成a人片在线一区二区| 色综合亚洲欧美另类图片| 中文字幕高清在线视频| 久久天躁狠狠躁夜夜2o2o| 精品人妻偷拍中文字幕| 亚洲久久久久久中文字幕| 国产黄a三级三级三级人| 日韩在线高清观看一区二区三区 | 久久精品国产亚洲网站| 91午夜精品亚洲一区二区三区 | 欧美一区二区精品小视频在线| 别揉我奶头~嗯~啊~动态视频| 精品人妻偷拍中文字幕| 国产精品美女特级片免费视频播放器| 一区二区三区高清视频在线| 可以在线观看的亚洲视频| 一级a爱片免费观看的视频| 午夜福利成人在线免费观看| 久久亚洲精品不卡| 国产在视频线在精品| 美女黄网站色视频| 黄色欧美视频在线观看| 亚洲熟妇中文字幕五十中出| av专区在线播放| 琪琪午夜伦伦电影理论片6080| 欧美不卡视频在线免费观看| 精品一区二区三区视频在线| 少妇丰满av| 国产高清三级在线| 男女边吃奶边做爰视频| 欧美日韩国产亚洲二区| 校园人妻丝袜中文字幕| 成人高潮视频无遮挡免费网站| 毛片一级片免费看久久久久 | 亚洲av美国av| 国产欧美日韩精品一区二区| 亚洲精品国产成人久久av| 老熟妇乱子伦视频在线观看| 国产成人av教育| 免费av毛片视频| 尾随美女入室| 久久午夜亚洲精品久久| 一a级毛片在线观看| 国产激情偷乱视频一区二区| 毛片一级片免费看久久久久 | 好男人在线观看高清免费视频| 欧美日本亚洲视频在线播放| 国产又黄又爽又无遮挡在线| 亚洲av一区综合| 国产老妇女一区| 成年免费大片在线观看| 欧美黑人欧美精品刺激| 少妇熟女aⅴ在线视频| 国产精品一区www在线观看 | 免费一级毛片在线播放高清视频| 亚洲精华国产精华精| 国产av在哪里看| 国产在线精品亚洲第一网站| 91在线精品国自产拍蜜月| 一进一出抽搐gif免费好疼| 日本三级黄在线观看| 舔av片在线| 精品一区二区三区视频在线观看免费| 欧美成人免费av一区二区三区| 一进一出好大好爽视频| 啦啦啦观看免费观看视频高清| 亚洲性夜色夜夜综合| 精品午夜福利视频在线观看一区| 不卡一级毛片| 国产麻豆成人av免费视频| 亚洲精品日韩av片在线观看| 老司机午夜福利在线观看视频| 99久久中文字幕三级久久日本| 少妇的逼水好多| 一区二区三区四区激情视频 | 18禁黄网站禁片午夜丰满| 久久人人精品亚洲av| 精品日产1卡2卡| 亚洲综合色惰| 日韩一区二区视频免费看| 我要搜黄色片| 麻豆成人av在线观看| 午夜激情欧美在线| 精品国产三级普通话版| 国产又黄又爽又无遮挡在线| 舔av片在线| 国产老妇女一区| 春色校园在线视频观看| 亚洲精品一区av在线观看| 日韩欧美国产在线观看| 综合色av麻豆| 婷婷精品国产亚洲av| 18禁黄网站禁片午夜丰满| 午夜a级毛片| 免费观看精品视频网站| 天美传媒精品一区二区| 99视频精品全部免费 在线| 日韩 亚洲 欧美在线| 精品人妻1区二区| 亚洲av不卡在线观看| 久久久久久久久大av| 窝窝影院91人妻| 狠狠狠狠99中文字幕| 国产精品亚洲一级av第二区|