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

    異重流潛入現(xiàn)象探討Ⅰ:水槽實(shí)驗(yàn)與理論分析成果回顧

    2018-05-11 07:42:12范家驊
    水利學(xué)報 2018年4期
    關(guān)鍵詞:異重流交界面渾水

    范家驊,祁 偉,戴 清

    (中國水利水電科學(xué)研究院,北京 100048)

    1 異重流潛入及其工程問題

    異重流的形成是兩種不同密度的流體相遇時,在一定條件下較重或較輕的流體從有壓水流或無壓水流向另一種不同密度的流體開始過渡到分層流(異重流)的相對運(yùn)動過程。需要研究的問題是:從非分層流過渡到分層流需要遵循哪些條件。這種過渡現(xiàn)象即在潛入點(diǎn)下游形成底部異重流流動,或稱之為異重流渾水楔。

    (1)水庫入口處。洪峰期間或非汛期挾沙水流進(jìn)入水庫末端回水區(qū),流至一定的水深處潛入庫底形成底部異重流,即異重流渾水楔。潛入點(diǎn)下游渾水楔長度隨洪峰歷時而變,洪峰歷時短時,渾水楔流不到壩址,全部淤在水庫內(nèi);洪峰歷時長時,可通過壩底孔排出部分泥沙;而潛入點(diǎn)上游壅水區(qū)會形成三角洲泥沙淤積并導(dǎo)致河段水位抬高[1-2]。

    (2)盲腸河段與河道的交匯處(船閘和引航道)。河道挾沙水流與盲腸河段內(nèi)靜止水體相遇因密度不同產(chǎn)生相對運(yùn)動(即所謂交換水流),挾沙水流潛入底部并流進(jìn)盲腸河段,同時上層清水則自盲腸河段流出進(jìn)入河道。由于河道含沙水流這樣持續(xù)不斷地進(jìn)入盲腸段而造成盲腸河段的累積淤積。

    (3)河道交匯處。河道主流含沙量大于支流含沙量時,主流渾水潛入支流形成底部異重流并向支流上游方向運(yùn)動,且上溯至一定距離,形成一段具有一定長度的渾水楔;而當(dāng)支流含沙量大于主流含沙量時,支流渾水則潛入主流,也形成底部異重流和具有一定長度的渾水楔。工程設(shè)計(jì)或研究的內(nèi)容主要為楔內(nèi)泥沙沿程淤積和淤積處理方法。

    (4)河口與海域的交匯處。海區(qū)鹽水與河道水流在河口處相遇,在一定條件下會形成向上游的鹽水楔運(yùn)動;而鹽水中含沙量較大時(受海灘風(fēng)浪掀沙影響),上溯水流與河道水流相遇,在一定條件下會形成含沙分層運(yùn)動的渾水楔,渾水楔沿程泥沙淤積會抬高航道底高,需挖除淤積礙航泥沙,才能保持一定航深的航道通航。另外河道挾沙水流流出河口,因其含沙量較小,渾水密度小于海水密度,渾水會向海面擴(kuò)散,形成上層異重流。

    2 潛入現(xiàn)象研究回顧

    異重流潛入現(xiàn)象的水槽實(shí)驗(yàn)、水力理論分析和模型數(shù)值計(jì)算工作,國內(nèi)國外不少學(xué)者進(jìn)行過研究,以下將分述和回顧各家研究成果;同時本文還將著重于渾水異重流在等寬水槽內(nèi)的現(xiàn)象,在筆者先前工作的基礎(chǔ)上做進(jìn)一步的分析,與前人工作進(jìn)行對比,以獲得對異重流潛入現(xiàn)象進(jìn)一步的認(rèn)識,并對潛入現(xiàn)象的水力分析中若干問題進(jìn)行討論,為工程設(shè)計(jì)中水力計(jì)算提供參考[1-2]。

    2.1 水庫壅水區(qū)渾水潛入點(diǎn)實(shí)驗(yàn) 渾水異重流潛入點(diǎn)水槽實(shí)驗(yàn)有范家驊[3]、蘆田和男[4]、曹如軒[5]、曹如軒與錢善琪[6]、俞維升[7]、姚鵬與王興奎[8-9]和焦恩澤[10]等。筆者分別在寬15 cm水槽和寬50 cm水槽內(nèi)進(jìn)行過渾水潛入實(shí)驗(yàn),異重流潛入過程如圖1所示。

    水槽實(shí)驗(yàn)中一般觀測記錄潛入點(diǎn)的水深、含沙量、流速以及潛入點(diǎn)潛入后交界面線的沿程變化,并進(jìn)行部分測次若干斷面的含沙量垂線分布和流速分布的測量。對潛入點(diǎn)測定數(shù)據(jù)進(jìn)行量綱分析,可獲得潛入點(diǎn)斷面密度Froude數(shù)Fp的值。根據(jù)Schijf和Schonfeld[11]的異重流雙層漸變流方程,潛入點(diǎn)水深過渡到潛入點(diǎn)異重流水深之間會存在一個交界面線的拐點(diǎn),該處密度Froude數(shù)等于1;而潛入點(diǎn)水深大于該處臨界密度Froude數(shù)中的水深,故可見潛入點(diǎn)處密度Froude數(shù)Fp小于1。資料分析得潛入點(diǎn)處密度Froude數(shù):

    圖1 異重流潛入圖

    式中:up為異重流潛入點(diǎn)流速;hp為潛入點(diǎn)水深;g′=(?ρ/ρ)g。1961年筆者又在長50 m、寬0.5 m、高2 m水槽內(nèi)進(jìn)行異重流潛入點(diǎn)實(shí)驗(yàn),進(jìn)沙粒徑d90=0.0145~0.05mm,d50=0.0055~0.008 mm,觀測沿程渾水水深、含沙量和泥沙粒徑的變化,4次實(shí)驗(yàn)后得Fp=0.58~0.66。

    蘆田和男[4]進(jìn)行渾水實(shí)驗(yàn),分析得:其中hp為潛入點(diǎn)水深;S為底部比降;q為單寬流量。曹如軒等曾對上式利用他們的實(shí)驗(yàn)資料進(jìn)行比較,其系數(shù)大于蘆田和男的0.365,一般在0.4~0.57之間,平均值為0.44。

    曹如軒[5]利用3個水槽進(jìn)行含沙量6.5~715 kg/m3的渾水潛入實(shí)驗(yàn),其中根據(jù)含沙量在30 kg/m3以下的資料得:而渾水含沙量大于100 kg/m3時Fp值較小,且隨著含沙量的增加Fp值減小。含沙量在100~360 kg/m3時,F(xiàn)p=0.4~0.2。曹如軒等對含沙量高且流體黏滯系數(shù)大的高含沙水流,利用Bingham體有關(guān)參數(shù)和流體阻力系數(shù)聯(lián)系來分析潛入點(diǎn)密度Froude數(shù)與含沙濃度的某種關(guān)系,分析得到符合實(shí)驗(yàn)數(shù)據(jù)變化趨勢的結(jié)果。曹如軒與錢善琪[6]進(jìn)行含粗沙高含沙量潛入點(diǎn)水槽實(shí)驗(yàn),得到與低含沙量潛入類似的Fp值結(jié)果。俞維升[7]用高岺土和石英沙挾沙水流以及鹽水進(jìn)行潛入點(diǎn)水槽實(shí)驗(yàn),觀察到潛入點(diǎn)向下游移動至一定處趨于穩(wěn)定的現(xiàn)象,且平均Fp≈0.71;同時通過沿程流速測量與分析,計(jì)算得潛入后異重流流量沿程增加的現(xiàn)象。姚鵬與王興奎[8-9]在長而寬的水槽內(nèi)進(jìn)行4種底坡的渾水異重流潛入點(diǎn)試驗(yàn),槽寬1.2 m、長63.8 m,觀察到潛入點(diǎn)Fp值因底坡加大而降低,變化范圍在0.67~0.57之間。焦恩澤[10]進(jìn)行不同含沙量(12~479 kg/m3)的潛入點(diǎn)水槽實(shí)驗(yàn),其中大量試驗(yàn)為高含沙量異重流實(shí)驗(yàn)。同曹如軒實(shí)驗(yàn)類似,F(xiàn)p值隨含沙量增加而變小,并利用Bingham體參數(shù)與潛入點(diǎn)建立一定關(guān)系。范家驊[12]考慮小浪底水庫的含沙量范圍可用曹如軒、焦恩澤等人的高含沙異重流潛入點(diǎn)數(shù)據(jù),并采用量綱分析方法,得到含沙量在400 kg/m3以下范圍內(nèi)的Fp平均關(guān)系為Fp=0.9c0.1,c為含沙量。

    近年來國內(nèi)關(guān)于水庫潛入點(diǎn)Fp的研究公開發(fā)表的論文近20余篇,主要為小浪底水庫高含沙量異重流的控制運(yùn)用中進(jìn)行異重流潛入點(diǎn)參數(shù)的研究,并對包括小浪底水庫及其它水庫的觀測資料、水槽實(shí)驗(yàn)和模型試驗(yàn)資料進(jìn)行有關(guān)水沙因素與Fp的影響分析。如李濤和張俊華等[13]、李書霞等[14]、李濤和夏軍強(qiáng)等[15]的研究,可見在高含沙量異重流潛入點(diǎn)處,F(xiàn)p參數(shù)與含沙量值有一定影響。

    各家實(shí)驗(yàn)結(jié)果以及渾水模型、數(shù)值計(jì)算結(jié)果列于表1。

    表1 渾水異重流潛入點(diǎn)水槽實(shí)驗(yàn)、數(shù)值計(jì)算與水力分析

    2.2 鹽水與冷熱水潛入點(diǎn)實(shí)驗(yàn) 除了渾水實(shí)驗(yàn),尚有利用鹽水和冷熱水模擬渾水進(jìn)入水庫壅水區(qū)潛入的水槽實(shí)驗(yàn),有Keulegan[17]、Singh與Shah[18]、Itakura與Kishi[19]、Fukuoka等[20]、Kan與Tamai[21]等人的鹽水實(shí)驗(yàn),以及Farrell與Stefan[22]、Akiyama與Stefan[23-24]等人的冷熱水實(shí)驗(yàn)。各家實(shí)驗(yàn)結(jié)果列于表2。

    Keulegan[17]進(jìn)行無潮河口鹽水水槽實(shí)驗(yàn),觀測進(jìn)槽的鹽水異重流前鋒加速運(yùn)動規(guī)律。令進(jìn)口水深即為潛入點(diǎn)水深,以hp表示;異重流前鋒速度以ud表示。用資料分析可得:hp=2hd且0.57,故有

    表2 鹽水和冷熱水潛入點(diǎn)水槽實(shí)驗(yàn)

    Singh與Shah[18]用鹽水在水槽內(nèi)進(jìn)行了異重流潛入實(shí)驗(yàn)。水槽長14 m,寬0.4 m,進(jìn)槽單寬流量0.5~135 cm2/s,密度差為0.0005~0.013 g/cm3,異重流水深為1.5~17 cm,潛入點(diǎn)水深3~22.5 cm,水槽底坡S為0.005~0.02。他們使用量綱分析潛入點(diǎn)水深Re(Reynolds數(shù))和底坡S諸因素的函數(shù),并點(diǎn)繪潛入點(diǎn)水深與Re和S底坡的關(guān)系,表明在Re=600~11000以及S=0.0056~0.0215之間,此兩因素對潛入點(diǎn)水深無多大影響;在低流量和高密度的某些測次,水面有明顯的分界線,但未觀測到有向下游流動的異重流,這些測次的值小于2 cm,均未包括在分析之中。最后分析對較大值有:

    以密度Froude數(shù)Fp可表示為:

    Fukuoka等[20]在長6.9 m、高0.9 m、寬0.20 m水槽內(nèi)進(jìn)行鹽水潛入點(diǎn)實(shí)驗(yàn)。槽底比降為1/10的實(shí)驗(yàn)7次、比降1/60的實(shí)驗(yàn)7次,觀測潛入點(diǎn)交界面形狀、流速分布、鹽分分布等;并分析點(diǎn)繪潛入點(diǎn)水深與的關(guān)系;還加入石橋和巖崎在水庫中數(shù)次實(shí)驗(yàn)數(shù)據(jù),得平均值:Fu?kuoka等認(rèn)為他們的實(shí)驗(yàn)結(jié)果與他們在文中介紹Benjamin的理論分析結(jié)果接近,并計(jì)算出摻混系數(shù)E在0.046×10-2~1.8×10-2之間,還點(diǎn)繪實(shí)測和計(jì)算的沿程交界面高度的變化。

    Kan與Tamai[21]在槽長6 m、高0.9 m、寬0.3 m、底坡為0.1~0.25的水槽內(nèi)進(jìn)行鹽水異重流潛入實(shí)驗(yàn),進(jìn)槽鹽水密度為1.00006~1.056 g/cm3,q=2.4~2.5 cm3/cm,得潛入點(diǎn)Fp平均值0.63,此值系26次數(shù)據(jù)0.45~0.92的平均值。

    Farrell與Stefan[22]應(yīng)用k-ε模型計(jì)算了潛入點(diǎn)水流現(xiàn)象,并進(jìn)行了水槽冷熱水實(shí)驗(yàn)。槽長40呎、寬16.5呎、底坡S=0.047,觀測潛入點(diǎn)處水深、流速和溫度差以及潛入?yún)^(qū)下游異重流水深、流速及溫度等。實(shí)驗(yàn)共進(jìn)行7次,實(shí)驗(yàn)資料分析得:Fp=0.69。實(shí)驗(yàn)還觀測了1次摻混系數(shù):式中Qm為異重流流量,Q0為進(jìn)槽流量。

    Akiyama與 Stefan[23-24]在槽寬61 cm、深30.5 cm且具有一邊擴(kuò)展角1°、3°和7°的水槽內(nèi)進(jìn)行冷熱水異重流潛入點(diǎn)實(shí)驗(yàn),研究槽寬擴(kuò)寬對潛入的影響,實(shí)測獲得Fp=0.56~0.89之間,平均值為0.68;實(shí)驗(yàn)中還測得潛入點(diǎn)下游異重流水深hd和潛入點(diǎn)水深hp的比值在0.65~0.9之間;并實(shí)測摻混系數(shù)γ=0~0.31,摻混系數(shù)的定義為Q0為進(jìn)槽流量,Qd、Qp分別為異重流流量和潛入點(diǎn)流量。該文中點(diǎn)繪Fp和F0的關(guān)系,可看出在1°時,F(xiàn)p在0.6~0.7之間;3°時Fp=0.65~0.8之間;7°時Fp=0.6~0.8之間;總的范圍在0.6~0.8之間,可見F0對Fp的影響不大。關(guān)于此關(guān)系,Stefan與Johnson[25]在他們論文中補(bǔ)充了擴(kuò)展角3°共計(jì)5個新的測點(diǎn),點(diǎn)繪在一起,顯示Fp與F0和δ兩者無多大關(guān)系。

    2.3 “交換水流”渾水及鹽水潛入實(shí)驗(yàn)與水力分析 引航道口門渾水潛入形成渾水楔、船閘鹽水交換以及無潮河口入海處含鹽水流潛入河口段形成鹽水楔示意圖見圖2。

    圖2 引航道河口門異重流潛入圖

    這類船閘“交換水流”潛入實(shí)驗(yàn)的研究重點(diǎn)在于:求取船閘閘門打開后鹽水進(jìn)入儲有清水的船閘時潛入形成底層水流的鹽水和前鋒速度兩者的關(guān)系,可用密度Froude數(shù)的值表示,式中符號:前鋒速度uf(即異重流流速ud),河口水深H(即潛入點(diǎn)水深hp),清水與鹽水的密度差?ρ/ρ。

    進(jìn)行相關(guān)理論分析和鹽水實(shí)驗(yàn)工作的有O'Brien 與 Cherno[26]、 Yih C S(易 家 訓(xùn))[27]、Keulegan[28]、Middleton[29]、Rozovsky[30]、Kersey和 Hsu[31]、 Lam[32]、Barr[33]、Rottman 與 Simp?son[34]。利用渾水進(jìn)行引航道異重流潛入實(shí)驗(yàn)的有范家驊[1]。此外陳國謙和李行偉[35]利用k-ε(BN6)模型進(jìn)行了數(shù)值計(jì)算。實(shí)驗(yàn)結(jié)果與水力分析列于表3。

    至于各家水槽實(shí)驗(yàn),由于上下層水深不同以及交界的阻力系數(shù)的影響,F(xiàn)p值略大于或略小于0.25。

    水庫壅水區(qū)和引航道口門這兩類異重流潛入,因?yàn)闈撊胩幍膸缀魏退鳁l件不同,主要由于上層水流流速ua值在水庫中很小,分析時可忽略不計(jì);而在引航道中ua接近ud。故兩者的潛入點(diǎn)密度Froude數(shù)Fp有差異。綜上所述:水庫中Fp在0.5~0.8之間,引航道中Fp在0.25左右。

    表3 船閘交換水流鹽水、渾水實(shí)驗(yàn)與水力分析

    3 異重流潛入現(xiàn)象的理論分析成果

    3.1 von Karman的理論分析成果 von Karman[36]應(yīng)軍事部門的求詢,討論一種理想流體模型。圖3為較重液體ρ2在較輕液體ρ1中運(yùn)動,考慮較輕液體在無限深度假設(shè)前提下,推導(dǎo)下層較重液體恒定運(yùn)動的流動速度c1。

    von Karman應(yīng)用A點(diǎn)和B點(diǎn)的Bernoulli方程:

    在以一定速度運(yùn)動的框架之中水流為恒定,因?yàn)橄鄬τ谶@框架,水流為靜止:

    因此有:

    3.2 Benjamin的理論分析成果 Benjamin[37]認(rèn)為von Karman利用能量方程并忽略滯流點(diǎn)前端頭部高密度水流受到與上層水流混合的能量損失是不能接受的,雖然Benjamin的推理得到同樣的c1結(jié)果。易家訓(xùn)[38]的看法是von Karman雖然在這一點(diǎn)上是一個疏忽,但使用Bernoulli方程所取上下游兩點(diǎn),如果在上游的點(diǎn)取的位置更上游一些,下游的點(diǎn)更下游一些,使用Bernoulli方程則無問題。

    Benjamin分析兩平板之間的空穴流的空穴移動速度,分析中忽略黏性和表面張力,上游無窮遠(yuǎn)處水深為d,流速為c1,下游無限遠(yuǎn)處水深為h,流速為c2。概化圖見圖4。

    沿表面的Bernoulli方程,O點(diǎn)表面的壓力為零:

    無限遠(yuǎn)處的壓力p0:

    上下游遠(yuǎn)處的斷面壓力分布為靜水壓力分布,動量守恒方程為:

    圖3 von Karman異重流概化圖

    圖4 Benjamin空穴流概化圖

    連續(xù)條件為c1d=c2h,解得:

    其次Benjamin考慮水流的能量損失:假定水流有水頭損失,流速c2在很遠(yuǎn)的下游會變成均勻流,定義Δ為水頭損失,則有:

    此式與力平衡所得的下式:

    相等,可得:當(dāng)h<0.5d時,Δ為負(fù)值,故在實(shí)際中不能發(fā)生。利用以上各式,可得當(dāng)h>0.5d時,有:

    3.3 Singh和Shah的理論分析成果 Singh和Shah[18]除了進(jìn)行鹽水潛入點(diǎn)水槽實(shí)驗(yàn)外,還用動量方程分析潛入現(xiàn)象,動量方程為:

    式中:左邊代表動量的改變率,β1和β2為動量系數(shù);v1和v2為兩斷面的平均流速;P1和P2為兩斷面的壓力;F1和F2為底部和交界面的剪力阻力;Ws為兩斷面間液體重量的分力。

    Singh在分析中作出如下假定:水壓力為靜水壓力分布;每一點(diǎn)的水流方向平行于槽底;潛入長度L和底坡S均為小值,故第2斷面總水深等于第1斷面的水深,即潛入點(diǎn)水深;交界面上的阻力F2假定平行于槽底。概化圖見圖5。

    為簡化分析,忽略F1、F2和Ws小值,可得:

    對上式的分析采用三種簡化:第一種是忽略潛入點(diǎn)斷面與下游異重流斷面之間底部阻力和交界面阻力,并忽略動量流速分布系數(shù),得:第二種考慮動量流速分布系數(shù),計(jì)算得該系數(shù)β1和β2等于1.18,則:第三種考慮潛入點(diǎn)潛入?yún)^(qū)距離之內(nèi)的底部和交界面阻力的影響,Singh從均勻異重流中計(jì)算得平均總阻力系數(shù)f為底部阻力系數(shù)與交界面系數(shù)之和,其值為0.13,最后得:

    Singh按三種關(guān)系式計(jì)算的潛入點(diǎn)水深值同實(shí)驗(yàn)測得的潛入點(diǎn)水深值點(diǎn)繪在同一圖上對比,顯示三條平均線在水深較小部分計(jì)算值偏大,在水深較大時則偏差較小。

    3.4 Savage和Brimberg的理論分析成果 Savage與Brimberg[39]采用類似Benjamin的方法分析潛入點(diǎn)的水力特性。應(yīng)用Bernoulli原理,分析恒定二維兩種密度ρ1和ρ2的流體,斷面1處于潛入點(diǎn)上游遠(yuǎn)處,流速為u1,水深為h1;斷面2處于潛入點(diǎn)的下游遠(yuǎn)處,流速為u2,水深為h2。概化圖見圖6。設(shè)

    O點(diǎn)為滯流點(diǎn),沿自由水面的O點(diǎn)的上游應(yīng)用Bernoulli原理;其次應(yīng)用Bernoulli原理自O(shè)點(diǎn)沿交界面下游;并設(shè)在潛入點(diǎn)上游相當(dāng)遠(yuǎn)處和下游相當(dāng)遠(yuǎn)處水壓力為靜水壓力,運(yùn)用動量平衡方程可導(dǎo)得:

    圖5 Singh潛入點(diǎn)概化圖

    圖6 Savage潛入點(diǎn)概化圖

    令潛入點(diǎn)水深h1即為hp,可表示為:

    Savage用以上推導(dǎo)的關(guān)系式同Singh的潛入點(diǎn)實(shí)驗(yàn)水深比較,兩者較接近,但測點(diǎn)分散。實(shí)驗(yàn)數(shù)據(jù)顯示除以外尚有其它參數(shù)的影響,如底坡和交界面剪應(yīng)力。Savage考慮到第一部分的水力分析工作忽略了底部比降和剪應(yīng)力,故利用Schijf與Schonfeld[11]的雙層漸變流一維運(yùn)動方程來分析潛入點(diǎn)現(xiàn)象,分析中令水流為恒定及上層流速為零,對動力方程進(jìn)行數(shù)值計(jì)算。

    Savage為了定性描述潛入點(diǎn)下游的交界面形狀,利用底層漸變流方程,將潛入點(diǎn)處的底部作為坐標(biāo)零點(diǎn),垂向坐標(biāo)z代表水深,縱向坐標(biāo)為距離x,推求漸變流方程的解。這樣就類似于明渠漸變水流,可分為數(shù)種水面線類型。

    從潛入點(diǎn)起向下游的交界面線在緩坡和急坡情況下可能出現(xiàn)的是接近M2和S2的曲線。數(shù)值計(jì)算的結(jié)果可建立Fp為f,α和S(底坡)的函數(shù)。對于緩坡有:

    上式是在f=0.01~0.09,α=0.2~0.8,α=fi/f以及緩坡范圍內(nèi)的計(jì)算結(jié)果。

    而在急流條件下,計(jì)算表明其性質(zhì)與緩坡情況類似。但有一點(diǎn)不同,當(dāng)潛入點(diǎn)密度Froude數(shù)為Fp=0.8至0.845時,交界面曲線出現(xiàn)劇烈變化。Savage分析Singh的實(shí)驗(yàn)資料,其潛入點(diǎn)Fp值小于0.8,在0.3~0.8之間,點(diǎn)繪的關(guān)系線,并以參數(shù)Fp=0.3~0.8作直線,測點(diǎn)絕大部分在Fp=0.3~0.8線的范圍之內(nèi)。

    3.5 Jain的理論分析成果 Jain[40]著文檢驗(yàn)Savage的論文,提出一些不同意見。他認(rèn)為用能量方程于潛入點(diǎn)附近的水流在物理上是不可能的。潛入點(diǎn)附近的能量耗損不可忽視,因此他利用動量方程和能量方程(考慮能量損失Δ)進(jìn)行分析,概化圖見圖7。并得到以下關(guān)系式:

    圖7 Jain潛入點(diǎn)概化圖

    Jain指出Savage和Brimberg的分析,F(xiàn)1=0.5為以上關(guān)系式當(dāng)HL=0和且β=0.5時的特解。Jain對Δ值作出的假定,表示β與的相應(yīng)關(guān)系。Jain計(jì)算出在不同β值時動量方程的解,并且指出對于某一已知β值存在一個F1的極大值(F1max),大于此值時不能發(fā)生潛入現(xiàn)象。

    Jain的論文第二部分分析交界面曲線的計(jì)算,認(rèn)為積分的方向是敏感的,計(jì)算應(yīng)從控制斷面算起。他利用Schijf與Schonfeld[11]的雙層漸變流的一維方程進(jìn)行計(jì)算,并假定其上層水流的平均流速為零。上層水面為水平,水深為a1,下層水深為a2,在潛入點(diǎn)處a1=0,如圖5所示圖形。可得:

    上式僅包含兩個無量綱參數(shù)λ和ηc,可用不同的λ和ηc值和不同下游條件求得數(shù)值解,其積分是從下游控制斷面向上游方向計(jì)算。

    對緩坡的計(jì)算,設(shè)一定的阻力系數(shù)α=0.5、λ=1 3和ηc<1(即臨界水深小于正常水深,如ηc=0.5)以及若干下游條件向上游進(jìn)行積分。所求得的對于大范圍的下游條件的交界面曲線在潛入點(diǎn)附近連接,且潛入點(diǎn)水深值與下游邊界條件無關(guān)。但當(dāng)積分是從上游向下游進(jìn)行計(jì)算,即使很小的差別,卻會造成不同的交界面曲線。

    對陡坡(ηc>1)計(jì)算兩種典型的交界面曲線ηc=1.3以及λ=1 3(α=0.5),交界面曲線S1和S3均為緩流。

    對于潛入水深,無量綱潛入點(diǎn)水深ηp僅為ηc和λ(或α)的函數(shù),有下列近似關(guān)系:

    3.6 朱鵬程的理論分析成果 朱鵬程[41]首先對異重流潛入條件按以下概化圖形(見圖8)進(jìn)行了理論分析,得:

    得:

    曹如軒等[5]用同樣方法,得上述相同公式。朱鵬程注意到上列動量方程導(dǎo)得的方程類似水躍前后的共軛水深,因此其間須經(jīng)過臨界水深。他認(rèn)為要產(chǎn)生共軛水深,F(xiàn)p必須小于1且Fd必須大于1。

    蘆田和男[4]對異重流潛入現(xiàn)象進(jìn)行分析和渾水實(shí)驗(yàn),求解潛入點(diǎn)水深。解得:

    3.7 Akiyama和Stefan的理論分析成果 Akiyama與Stefan[42-43]分析潛入現(xiàn)象時,首先利用動量公式并引進(jìn)摻混系數(shù)γ,推導(dǎo)潛入點(diǎn)密度Froude數(shù)與潛入點(diǎn)下游異重流厚度等因子的關(guān)系。

    圖9概化圖形中的CVⅠ動量方程假定可忽略混合區(qū)的底部阻力以及忽略比降的影響,得下式:

    他們定義摻混系數(shù)γ(水庫清水混入底部異重流的水量):

    CVⅡ的動量方程:

    圖8 朱鵬程異重流潛入點(diǎn)概化圖

    圖9 Akiyama異重流潛入點(diǎn)概化圖

    以上兩動量方程合并,得:

    式(19)中Fd為潛入點(diǎn)下游的異重流密度Froude數(shù),此式表明Fd和和γ的關(guān)系。

    為了分析潛入點(diǎn)水沙因子與底部阻力、交界面阻力以及下游異重流的影響,Akiyama作出假定:潛入點(diǎn)下游異重流在急坡時Fd等于臨界值;在緩坡時Fd等于正常值Fn。他們利用Ellison與Turner[44]的具有摻混系數(shù)的異重流緩變方程,求解緩坡和急坡時的與摻混系數(shù)、阻力系數(shù)、底部比降以及異重流流速分布系數(shù)的關(guān)系式。由于所求得的關(guān)系式中阻力系數(shù)與摻混系數(shù)的各項(xiàng)關(guān)系式過長,為了同Singh等人實(shí)驗(yàn)資料做比較,Akiyama假定在進(jìn)口處無摻混現(xiàn)象,即γ=0,故得下式:

    緩坡:

    陡坡:

    式中的異重流流速分布系數(shù)采用Ellison與Turner的實(shí)驗(yàn)值。(Ellison與Turner用鹽水做實(shí)驗(yàn)所得流速分布系數(shù)S1、S2值遠(yuǎn)較Parker等[45]渾水異重流實(shí)驗(yàn)測得值為?。?。Akiyama將此兩式與Singh的鹽水潛入實(shí)驗(yàn)資料進(jìn)行對比,實(shí)驗(yàn)值與計(jì)算值符合。

    3.8 Parker和Toniolo的理論分析成果 Parker和Toniolo[46]專門針對Akiyama與Stefan的論文[43]提出意見。他們指出:Akiyama文中在分析阻力系數(shù)和底部比降等影響時所假定潛入點(diǎn)下游的異重流運(yùn)動在緩坡時的Fd值為正常值、在急坡時的Fd值為臨界值是錯誤的,也是沒有必要的。他們認(rèn)為Akiya?ma的分析,用兩種CVⅠ和CVⅡ的動量守恒方程就足以闡明hd和hp的關(guān)系,而且能得到Fp和Fd兩者僅為摻混系數(shù)γ的關(guān)系式。

    Parker分析CVⅠ的動量方程,導(dǎo)得以下兩式:

    當(dāng)γ=0時,式(22)可得與朱鵬程推導(dǎo)相同的式(16)。

    從式(22)可計(jì)算k為Fp和γ的函數(shù),如果hp、up、εγ(來水密度以及γ為已知,則hd、ud和εd(異重流密度可從式(22)、連續(xù)方程求得。

    Parker根據(jù)CVⅡ的動量方程,導(dǎo)得:

    將式(24)代入式(22)

    Parker認(rèn)為利用CVⅡ的動量方程,有了式(24)就可直接計(jì)算Fp,如僅有CVⅠ的動量方程,所求得的式(22)還須先知道Fp和γ(或Fd和γ)用于計(jì)算k值,而用CVⅡ動量方程求得的式(24)和式(25),F(xiàn)p和k兩者可從已知γ時求得。

    4 異重流潛入現(xiàn)象的數(shù)學(xué)模型計(jì)算成果

    方春明等[16]對水庫異重流潛入條件采用動量方程分析得:

    根據(jù)兩種概化圖形,利用能量方程和動量方程進(jìn)行理論分析,探討了異重流在水庫中潛入點(diǎn)密度Froude數(shù)小于1的水沙條件。方春明等還建立了二維水流和懸沙數(shù)學(xué)模型進(jìn)行數(shù)值計(jì)算。為檢驗(yàn)?zāi)P徒Y(jié)果,選取曹如軒9組水槽潛入點(diǎn)實(shí)驗(yàn)所測資料進(jìn)行對比計(jì)算,計(jì)算得Fp=0.48~0.77,與實(shí)驗(yàn)Fp=0.45~0.75比較,兩者接近。而且他們對影響水流和泥沙模型計(jì)算結(jié)果的幾個因素進(jìn)行了潛入點(diǎn)流態(tài)和潛入點(diǎn)水深的對比計(jì)算;還計(jì)算了異重流潛入點(diǎn)的流場、清渾水交界面水深的縱向變化以及水面線的縱向變化,發(fā)現(xiàn)在壅水區(qū)水面線在潛入點(diǎn)附近出現(xiàn)倒比降。此計(jì)算結(jié)果與Savage提出的異重流潛入現(xiàn)象分析的概化圖形中下游斷面水位隆起高度Δ在定性上相符。

    Farrell與Stafan[22,47]對水庫中異重流潛入現(xiàn)象用k~ε模型進(jìn)行計(jì)算,模型包括水流動量方程和溫度方程,共進(jìn)行12次紊流模型方程計(jì)算,得,即Fp=0.5。潛入點(diǎn)水深與Singh和Shah、Farrell和Stefan的實(shí)驗(yàn)數(shù)據(jù)系數(shù)1.3相比大20%,F(xiàn)arrell和Stafan尚難解釋其差別的原因;其次這12次計(jì)算中,計(jì)算初期混合系數(shù)(用符號γ表示)在0.028~0.19之間,γ的變化可用下式表示,Q0為進(jìn)槽流量,Qd為異重流流量。

    Bournet等[48]的異重流潛入數(shù)值模型計(jì)算,應(yīng)用水流動量方程、k-ε模型以及溫度輸移方程,計(jì)算流速場和溫度在等寬槽和水槽其中一邊擴(kuò)展一個小的角度(1°~7°)兩種情況下的分布。在等寬槽的計(jì)算中引進(jìn)摻混系數(shù),其值在0.005~0.012之間,計(jì)算得潛入點(diǎn)水深此式hp計(jì)算值較Singh實(shí)驗(yàn)值大些。

    表4 水庫異重流潛入點(diǎn)理論分析與數(shù)值計(jì)算

    表5 引航道口門異重流潛入點(diǎn)理論分析

    以上各家理論分析和數(shù)值模型計(jì)算成果,匯總列于表4和表5。有關(guān)渾水異重流潛入實(shí)驗(yàn)結(jié)果分析,將在本文第Ⅱ篇中給予介紹。

    參 考 文 獻(xiàn):

    [1]范家驊.異重流與泥沙工程:實(shí)驗(yàn)與設(shè)計(jì)[M].北京:中國水利水電出版社,2011.

    [2]范家驊,等.異重流運(yùn)動的實(shí)驗(yàn)研究[J].水利學(xué)報,1959(5):30-48.

    [3]范家驊,等.異重流的研究和應(yīng)用[M].北京:水利電力出版社,1959.

    [4]EGASHIRA S,ASHIDA K.Condition that suspended load plunges to form a gravity current[C]//19th Proc.of the Conference on Natural Disaster Science.Japan,1978.

    [5]曹如軒,任曉楓,盧文新.高含沙異重流的形成與持續(xù)條件分析[J].泥沙研究,1984(2):1-10.

    [6]曹如軒,錢善琪,郭崇,等.粗沙高含沙異重流的運(yùn)動特性[J].泥沙研究,1995(2):64-73.

    [7]俞維升.水庫沉滓運(yùn)動特性之研究[D].臺北:臺灣大學(xué),1991.

    [8]姚鵬,王興奎.異重流潛入規(guī)律研究[J].水利學(xué)報,1996(8):77-83.

    [9]姚鵬.異重流運(yùn)動的試驗(yàn)研究[D].北京:清華大學(xué),1994.

    [10]焦恩澤.黃河水庫泥沙[M].鄭州:黃河水利出版社,2004.

    [11]SCHIJF J B,SCHONFELD J C.Theoretical Considerations on the Motion of Salt and Fresh Water[C]//Proceed?ings of Minnesota International Hydraulics Convention.Minnesota:Minneapolis,1953:321-333.

    [12]范家驊.關(guān)于水庫渾水潛入點(diǎn)判別數(shù)的確定方法[J].泥沙研究,2008(1):74-81.

    [13]李濤,張俊華,馬懷寶,等.異重流潛入重力修正系數(shù)研究[J].人民黃河,2012,34(7):28-29.

    [14]李書霞,夏軍強(qiáng),張俊華,等.水庫渾水異重流潛入點(diǎn)判別條件[J].水科學(xué)進(jìn)展,2012,23(3):363-368.

    [15]李濤,夏軍強(qiáng),張俊華,等.水庫異重流潛入點(diǎn)流速分布及其判別式改進(jìn)[J].工程科學(xué)與技術(shù),2017,49(2):62-68.

    [16]方春明,韓其為,何明民.異重流潛入條件分析及立面二維數(shù)值模擬[J].泥沙研究,1997(4):68-75.

    [17]KEULEGAN G H.Twelfth progress report on model laws for density currents:the motion of saline fronts in still water[R].National Bureau of Standards:U.S.Dept.of Commerce,1958.

    [18]SINGH B,SHAH C R.Plunging phenomenon of density currents in reservoirs[J].La Houille Blanche,1971,26(1):59-64.

    [19]ITAKURA T,RISHI T.Study on the turbidity density current in a reservoir[C]//16th Symp.on Science of Natural Disaster.Japan,1979:233-234.

    [20]FUKUOKA S,F(xiàn)UKUSHIMA Y,NAKAMURA K.Study on the plunge depth and interface form of density cur?rents in a two-dimensional reservoir[C]//土木學(xué)會論文報告集第302號 .1980:55-65.

    [21]KAN K,TAMA N.On the plunging point and initial mixing of the inflow into reservoirs[C]//Proc.25th Japanese Conf.on Hydraulics.Japan,1981:631-636.

    [22]FARRELL C J,STEFAN H G.Buoyancy induced plunging glow into reservoirs and coastal regions[R].St.An?thony Falls Hydraulic Laboratory,Minnesota:University of Minnesota,1986.

    [23]AKIYAMA J,STEFAN H.Gravity currents in lakes,reservoirs and coastal regions:Two-layer stratified flow analysis[R].St.Anthony Falls Hydraulic Laboratory,Minnesota:University of Minnesota,1987.

    [24]AKIYAMA J,STEFAN H.Onset of underflow in slightly diverging channels[J].Journal of Hydraulic Engineer?ing-Asce,1987,113(7):825-844.

    [25]STEFAN H,JOHNSON T.Negative buoyant flow in diverging channel.Ⅲ.Onset of underflow[J].Journal of Hydraulic Engineering,1989,115(4):423-436.

    [26]O'BRIEN M P,CHERNO J.Model law for motion of salt water through fresh[J].Trans.ASCE,1934(99):576-594.

    [27] YIH C S.A study of the characteristics of gravity waves at a liquid surface[D].Iowa:State Univ.of Iowa,1947.

    [28]KEULEGAN G H.An experimental study of the motion of saline water from locks into fresh water channels[C]//13th Progress Report on Model Laws of Density Current,National Bureau of Standards,1957.

    [29]MEDDLETON G V.Experiments on density and turbidity currents.I.motion of the head[J].Canadian J.of Earth Sciences,1966(3):523-546.

    [30]ROGOVSKY.International Symposium of Stratified Flows[C].Novosibirsk,1972:440-458.

    [31]KERSEY D G,HSü K.Energy relations of density current flows:an experimental investigation[J].Sedimentolo?gy,1976(23):761-789.

    [32]LAM S T.Experimental investigation of lock-release gravity current[D].Hong Kong:Univ.of Hong Kong,1995.

    [33]BARR D I H.Densimetric exchange flow in rectangular channels.II:Some observations of the structure of lock exchange flow[J].La Houille Blanche,1963(7):757-766.

    [34]ROTTMAN J W,SIMPSON J E.The initial development of gravity currents from fixed-volume release of heavy fluids[C]//IUTAM Symp.Delft,1983:347-359.

    [35]陳國謙,李行偉,李植.開閘式湍動異重流[J].中國科學(xué)(E輯),2002,32(6):754-764.

    [36]von KARMAN T.The engineer grapples with non-linear problems[J].Bull.Am.Math.Soc.,1940(46):615-683.

    [37]BENJAMIN T B.Gravity current and related phenoma[J].J.Fluid Mech.,1968(31):209-248.

    [38]YIH C S.Stratified Flows[J].Annual Review of Fluid Mechanics,1969.

    [39]SAVAGE S B,BRIMBERG J.Analysis of plunging phenomena in water reservoirs[J].J.Hydraul.Res.,1975,13(2):187-205.

    [40]JAIN S C.Plunging phenomena in reservoirs[C]//Proc.Symp.On Surface Water Impoundments,Minnesota:Minneapolis,1981.

    [41]朱鵬程.異重流的形成與衰減[J].水利學(xué)報,1981(5):52-59.

    [42]AKIYAMA J,STEFAN H.Theory of plunging flow into a reservoir:Internal memoratum[R].St.Anthony Falls Hydraulic Lab.,Univ.of Minnesota,1981.

    [43]AKIYAMA J,STEFAN H.Plunging flow into a reservoir:theory[J].J.Hydraul.Eng.,1984,110(4):484 499.

    [44]ELLISON T H,TURNER J S.Turbulent entrainment in stratified flows[J].J.Fluid Mech.,1959(6):423-428.

    [45]PARKER G,GARCIA M,F(xiàn)UKUSHIMA Y,et al.Experiments on turbidity currents over an erodible bed[J].J.Hydraul.Res.,1987,25(1):123-147.

    [46]PARKER G,TONIOLO H.Note on the analysis of plunging of density flows[J].Journal of Hydraulic Engineer?ing,2007,133(6):690-694.

    [47]FARRELL C J,Stefan H.Mathematical modelling of plunging reservoir flows[J].Journal of Hydraulic Re?search,1988,26(5):525-537.

    [48]BOURNET P,DARTUS D,TASSIN B,et al.Numerical investigation of plunging density current[J].J.Hy?draul.Eng.,1999,125(6):584-594.

    猜你喜歡
    異重流交界面渾水
    Promoting the International Dissemination of Chinese Culture Through International Chinese Language Education: A Case Study of Chinese-English Idiomatic Equivalence
    小浪底水庫異重流排沙效率分析
    鋼-混凝土交界面法向粘結(jié)性能研究
    泡沫(外一首)
    高速公路機(jī)電工程相關(guān)交界面管理組織建設(shè)探討
    水生植被影響異重流動力特性的試驗(yàn)分析
    雙塊式無砟軌道軌枕與道床交界面損傷特性分析
    中國鐵路(2019年1期)2019-03-23 01:11:58
    渾水變清
    幼兒畫刊(2018年4期)2018-04-11 03:38:39
    改進(jìn)的徑向基神經(jīng)網(wǎng)絡(luò)模型在水庫異重流泥沙淤積量模擬中的應(yīng)用
    異重流沉積過程和沉積特征研究
    化工管理(2017年9期)2017-03-05 12:05:20
    亚洲国产欧美日韩在线播放| 国产精品久久久久久人妻精品电影| 男女那种视频在线观看| 亚洲精品美女久久久久99蜜臀| 老司机深夜福利视频在线观看| 免费人成视频x8x8入口观看| 亚洲久久久国产精品| 搡老岳熟女国产| 久久精品亚洲精品国产色婷小说| 亚洲男人的天堂狠狠| a级毛片在线看网站| 黄片小视频在线播放| 亚洲精品在线观看二区| 成人18禁高潮啪啪吃奶动态图| 亚洲欧美激情综合另类| 久久香蕉激情| 欧美黑人巨大hd| 51午夜福利影视在线观看| 成人特级黄色片久久久久久久| 老司机福利观看| 男女视频在线观看网站免费 | 亚洲最大成人中文| 成人欧美大片| 亚洲九九香蕉| 国产亚洲精品第一综合不卡| 久久婷婷成人综合色麻豆| 亚洲中文字幕日韩| 少妇 在线观看| 女人爽到高潮嗷嗷叫在线视频| 国产1区2区3区精品| 久久国产精品人妻蜜桃| 亚洲人成伊人成综合网2020| 国产精品久久久久久人妻精品电影| 国产亚洲精品第一综合不卡| 国产精品免费视频内射| 免费看日本二区| 久久精品国产亚洲av高清一级| 黄色女人牲交| 丝袜在线中文字幕| 国产成+人综合+亚洲专区| 国产片内射在线| 免费在线观看成人毛片| 欧美性猛交╳xxx乱大交人| 午夜免费观看网址| 国产一区二区三区视频了| 怎么达到女性高潮| 国产麻豆成人av免费视频| 色精品久久人妻99蜜桃| 在线观看午夜福利视频| 性色av乱码一区二区三区2| 嫩草影院精品99| 日本精品一区二区三区蜜桃| 他把我摸到了高潮在线观看| 久久久国产成人精品二区| 精品乱码久久久久久99久播| 亚洲五月色婷婷综合| 欧美成人免费av一区二区三区| 视频区欧美日本亚洲| 日本免费一区二区三区高清不卡| 岛国在线观看网站| 亚洲国产欧美网| 亚洲欧美日韩高清在线视频| 欧美另类亚洲清纯唯美| 国产成人av教育| 18禁裸乳无遮挡免费网站照片 | 给我免费播放毛片高清在线观看| 日本a在线网址| 国产三级在线视频| 精品卡一卡二卡四卡免费| 免费在线观看成人毛片| 免费人成视频x8x8入口观看| 香蕉av资源在线| 国产又色又爽无遮挡免费看| 精品卡一卡二卡四卡免费| 日韩欧美在线二视频| 法律面前人人平等表现在哪些方面| 亚洲国产中文字幕在线视频| 午夜久久久久精精品| 国产精品 欧美亚洲| 午夜久久久久精精品| 级片在线观看| 国产日本99.免费观看| 欧美激情极品国产一区二区三区| 18禁观看日本| 欧美 亚洲 国产 日韩一| 色尼玛亚洲综合影院| 国内精品久久久久久久电影| 十八禁人妻一区二区| 啪啪无遮挡十八禁网站| 亚洲成人精品中文字幕电影| 又紧又爽又黄一区二区| 欧美日本视频| 精品久久久久久久久久免费视频| 中出人妻视频一区二区| 亚洲一区高清亚洲精品| 国产成人欧美在线观看| 亚洲成a人片在线一区二区| 成人国语在线视频| 人妻久久中文字幕网| 久久伊人香网站| 亚洲欧美激情综合另类| 亚洲成人久久爱视频| 午夜激情福利司机影院| 精品一区二区三区四区五区乱码| 国产亚洲欧美98| 老司机午夜十八禁免费视频| 亚洲五月婷婷丁香| 男人操女人黄网站| 日本撒尿小便嘘嘘汇集6| 国产成人欧美在线观看| 午夜精品久久久久久毛片777| 这个男人来自地球电影免费观看| 中文字幕另类日韩欧美亚洲嫩草| 天天添夜夜摸| 51午夜福利影视在线观看| 黄色视频不卡| 国产av一区二区精品久久| 日韩欧美国产一区二区入口| 麻豆成人av在线观看| 欧美不卡视频在线免费观看 | 这个男人来自地球电影免费观看| 国产成人影院久久av| 两性夫妻黄色片| 老司机午夜十八禁免费视频| 久久草成人影院| 好男人在线观看高清免费视频 | 国产激情久久老熟女| 99精品久久久久人妻精品| 一本一本综合久久| 久久久久久久精品吃奶| 搡老岳熟女国产| 国产精品,欧美在线| 长腿黑丝高跟| or卡值多少钱| 亚洲七黄色美女视频| 制服丝袜大香蕉在线| 制服丝袜大香蕉在线| 日韩av在线大香蕉| 好男人在线观看高清免费视频 | 熟女电影av网| 成人18禁高潮啪啪吃奶动态图| 久久久久国产一级毛片高清牌| 最近在线观看免费完整版| 国产精品av久久久久免费| 国产精品久久久人人做人人爽| 美女免费视频网站| 国产又爽黄色视频| av免费在线观看网站| 精品少妇一区二区三区视频日本电影| 亚洲中文字幕一区二区三区有码在线看 | 亚洲五月色婷婷综合| 一区二区三区高清视频在线| 久久久国产欧美日韩av| 久久精品夜夜夜夜夜久久蜜豆 | 十分钟在线观看高清视频www| 99精品久久久久人妻精品| 视频在线观看一区二区三区| 国产一区二区三区在线臀色熟女| 亚洲人成77777在线视频| 国产亚洲欧美在线一区二区| 久久久国产精品麻豆| 琪琪午夜伦伦电影理论片6080| 淫秽高清视频在线观看| 露出奶头的视频| 丝袜美腿诱惑在线| 久久久国产精品麻豆| av有码第一页| 88av欧美| 午夜a级毛片| 亚洲精品久久成人aⅴ小说| 国产av一区二区精品久久| 欧美性猛交黑人性爽| 免费人成视频x8x8入口观看| 久久久精品国产亚洲av高清涩受| 国产精品久久电影中文字幕| 啦啦啦韩国在线观看视频| 午夜激情av网站| 超碰成人久久| 香蕉久久夜色| 国产精品野战在线观看| 高潮久久久久久久久久久不卡| 国产av一区在线观看免费| 丰满人妻熟妇乱又伦精品不卡| 最新美女视频免费是黄的| 国产亚洲精品一区二区www| 亚洲色图av天堂| 国产精品国产高清国产av| 身体一侧抽搐| 亚洲国产精品999在线| 久久婷婷人人爽人人干人人爱| 女性生殖器流出的白浆| 亚洲精华国产精华精| 国产欧美日韩一区二区精品| 精品久久久久久成人av| 成人欧美大片| 少妇熟女aⅴ在线视频| 18禁观看日本| 巨乳人妻的诱惑在线观看| 国产精品久久久av美女十八| 亚洲全国av大片| 亚洲精品国产区一区二| 久久久久久人人人人人| 精品少妇一区二区三区视频日本电影| 亚洲最大成人中文| 亚洲国产精品合色在线| 伊人久久大香线蕉亚洲五| 色综合站精品国产| 亚洲精品色激情综合| 国产精品爽爽va在线观看网站 | av福利片在线| 特大巨黑吊av在线直播 | 久久中文看片网| 日韩视频一区二区在线观看| 国产成人精品久久二区二区91| 色综合亚洲欧美另类图片| 美女午夜性视频免费| 99精品欧美一区二区三区四区| 国产成年人精品一区二区| 亚洲男人天堂网一区| 国产精品九九99| 国产伦人伦偷精品视频| 99久久无色码亚洲精品果冻| 午夜免费激情av| 国产精品美女特级片免费视频播放器 | xxx96com| 国产精品亚洲美女久久久| 岛国在线观看网站| 熟女电影av网| 女人高潮潮喷娇喘18禁视频| 丁香欧美五月| avwww免费| 宅男免费午夜| 99国产综合亚洲精品| 十八禁人妻一区二区| 色综合欧美亚洲国产小说| 精品国产国语对白av| 亚洲 欧美 日韩 在线 免费| 亚洲国产看品久久| 听说在线观看完整版免费高清| 久久精品91无色码中文字幕| 国产av又大| 人人妻人人澡人人看| www国产在线视频色| 久久久久精品国产欧美久久久| 超碰成人久久| 国产三级在线视频| 好男人电影高清在线观看| 岛国视频午夜一区免费看| 国产精品久久久久久人妻精品电影| 亚洲欧美精品综合久久99| 色综合婷婷激情| 国产高清视频在线播放一区| 亚洲精品色激情综合| 免费女性裸体啪啪无遮挡网站| 亚洲人成电影免费在线| 亚洲欧美精品综合久久99| 韩国av一区二区三区四区| 久久久久久九九精品二区国产 | 中文亚洲av片在线观看爽| 国产精品98久久久久久宅男小说| 国产精品1区2区在线观看.| 欧美一区二区精品小视频在线| 欧美一区二区精品小视频在线| 国产成年人精品一区二区| 久久精品影院6| 女生性感内裤真人,穿戴方法视频| 久久精品成人免费网站| 1024视频免费在线观看| 亚洲欧美精品综合久久99| 91九色精品人成在线观看| 亚洲专区中文字幕在线| 精品乱码久久久久久99久播| 日本 欧美在线| 在线av久久热| 午夜日韩欧美国产| 国产精品影院久久| 欧美激情极品国产一区二区三区| 亚洲七黄色美女视频| 中文字幕精品亚洲无线码一区 | 91大片在线观看| 国产免费男女视频| 欧美成人午夜精品| 人人妻,人人澡人人爽秒播| 男女做爰动态图高潮gif福利片| 欧美zozozo另类| 午夜久久久久精精品| 黄网站色视频无遮挡免费观看| 亚洲国产毛片av蜜桃av| 久久中文看片网| 国产伦人伦偷精品视频| 国产真实乱freesex| 国产精品 欧美亚洲| 男男h啪啪无遮挡| 国产精品综合久久久久久久免费| 精品国产美女av久久久久小说| 国产午夜精品久久久久久| 久久亚洲真实| 老司机午夜十八禁免费视频| 可以在线观看毛片的网站| 亚洲av片天天在线观看| 久热爱精品视频在线9| 午夜福利在线在线| 在线免费观看的www视频| 十八禁人妻一区二区| 日本在线视频免费播放| www日本在线高清视频| 人人妻,人人澡人人爽秒播| 露出奶头的视频| 久久久精品国产亚洲av高清涩受| 日日爽夜夜爽网站| 久久香蕉精品热| 国产真人三级小视频在线观看| 亚洲久久久国产精品| 国产亚洲av嫩草精品影院| 成年女人毛片免费观看观看9| 久久狼人影院| 丝袜美腿诱惑在线| 在线观看日韩欧美| 夜夜躁狠狠躁天天躁| 国产一级毛片七仙女欲春2 | 欧美日韩一级在线毛片| 精品国内亚洲2022精品成人| 国产高清视频在线播放一区| 在线免费观看的www视频| 大香蕉久久成人网| 级片在线观看| 亚洲熟女毛片儿| 日本精品一区二区三区蜜桃| 精品熟女少妇八av免费久了| 波多野结衣巨乳人妻| 日韩欧美 国产精品| 99re在线观看精品视频| 国内毛片毛片毛片毛片毛片| 又大又爽又粗| 久久人妻av系列| 国产亚洲精品av在线| av中文乱码字幕在线| 他把我摸到了高潮在线观看| 一级片免费观看大全| 欧美丝袜亚洲另类 | 一卡2卡三卡四卡精品乱码亚洲| 在线十欧美十亚洲十日本专区| avwww免费| 啪啪无遮挡十八禁网站| 美女扒开内裤让男人捅视频| 亚洲 欧美一区二区三区| 老司机深夜福利视频在线观看| 国产精品一区二区精品视频观看| 国产视频一区二区在线看| 麻豆成人av在线观看| 一二三四在线观看免费中文在| 国产午夜福利久久久久久| 老熟妇仑乱视频hdxx| 天天躁夜夜躁狠狠躁躁| 一区二区三区激情视频| 亚洲电影在线观看av| 久久久久久免费高清国产稀缺| 亚洲国产高清在线一区二区三 | 精品福利观看| 一二三四在线观看免费中文在| 国产精品亚洲一级av第二区| 在线观看www视频免费| 亚洲一区二区三区色噜噜| 亚洲一区二区三区色噜噜| 一进一出抽搐gif免费好疼| 欧美一区二区精品小视频在线| 国产精品影院久久| 国产精品久久久av美女十八| 国产又色又爽无遮挡免费看| 免费av毛片视频| 天堂动漫精品| 久久性视频一级片| 丰满人妻熟妇乱又伦精品不卡| 中文亚洲av片在线观看爽| 亚洲av美国av| 久久午夜亚洲精品久久| 夜夜躁狠狠躁天天躁| 国内少妇人妻偷人精品xxx网站 | 亚洲第一欧美日韩一区二区三区| 啦啦啦韩国在线观看视频| 大型黄色视频在线免费观看| 色综合站精品国产| 午夜激情福利司机影院| 欧美性长视频在线观看| 久久精品aⅴ一区二区三区四区| 黄片大片在线免费观看| 中文字幕久久专区| 欧美又色又爽又黄视频| 日日干狠狠操夜夜爽| 亚洲欧美精品综合久久99| 又黄又爽又免费观看的视频| avwww免费| a在线观看视频网站| 中文在线观看免费www的网站 | 不卡av一区二区三区| 久久99热这里只有精品18| 久99久视频精品免费| 亚洲午夜理论影院| 久久精品夜夜夜夜夜久久蜜豆 | 亚洲av成人不卡在线观看播放网| 亚洲性夜色夜夜综合| 观看免费一级毛片| 精品电影一区二区在线| 两个人视频免费观看高清| 日本三级黄在线观看| 动漫黄色视频在线观看| 欧美性猛交黑人性爽| 1024手机看黄色片| 韩国av一区二区三区四区| 男女做爰动态图高潮gif福利片| 午夜免费激情av| 色综合婷婷激情| 亚洲国产看品久久| 麻豆成人午夜福利视频| 成人永久免费在线观看视频| 久久天堂一区二区三区四区| 日本a在线网址| 国产成人欧美在线观看| 国产精品亚洲美女久久久| 久久久久亚洲av毛片大全| 亚洲无线在线观看| 在线观看免费午夜福利视频| 可以在线观看毛片的网站| 国产一区二区三区视频了| 99精品欧美一区二区三区四区| 黄色片一级片一级黄色片| 国产成人精品久久二区二区免费| 啦啦啦免费观看视频1| 午夜福利视频1000在线观看| 欧美日本视频| 午夜a级毛片| 制服丝袜大香蕉在线| 最近最新中文字幕大全电影3 | 午夜免费观看网址| 热99re8久久精品国产| 国产亚洲精品久久久久久毛片| АⅤ资源中文在线天堂| 精品日产1卡2卡| 妹子高潮喷水视频| 久久久久久国产a免费观看| 中文字幕精品亚洲无线码一区 | 欧美日韩黄片免| 少妇被粗大的猛进出69影院| 这个男人来自地球电影免费观看| www.999成人在线观看| 一级a爱片免费观看的视频| 国产欧美日韩精品亚洲av| 日韩高清综合在线| 午夜久久久在线观看| 成人永久免费在线观看视频| 精品电影一区二区在线| 国产精品一区二区三区四区久久 | 天堂动漫精品| 成人手机av| 制服人妻中文乱码| 99国产极品粉嫩在线观看| 日韩av在线大香蕉| 亚洲欧美日韩高清在线视频| 久久欧美精品欧美久久欧美| 日韩av在线大香蕉| 国产成人av激情在线播放| 色在线成人网| 大香蕉久久成人网| 国产亚洲精品一区二区www| 男女做爰动态图高潮gif福利片| 在线观看免费午夜福利视频| 搞女人的毛片| 香蕉av资源在线| 亚洲成a人片在线一区二区| 久热爱精品视频在线9| 亚洲无线在线观看| 国产片内射在线| 精品人妻1区二区| 黄片大片在线免费观看| 久久精品亚洲精品国产色婷小说| 国产91精品成人一区二区三区| 国内精品久久久久久久电影| 欧美成人午夜精品| 免费观看人在逋| 我的亚洲天堂| 精品国产一区二区三区四区第35| 久9热在线精品视频| 在线观看舔阴道视频| 亚洲国产精品999在线| 麻豆成人午夜福利视频| 波多野结衣巨乳人妻| 色在线成人网| 老鸭窝网址在线观看| 中文在线观看免费www的网站 | 怎么达到女性高潮| www.自偷自拍.com| 中国美女看黄片| 亚洲中文日韩欧美视频| 1024香蕉在线观看| 久久草成人影院| 成人av一区二区三区在线看| cao死你这个sao货| 成人一区二区视频在线观看| 亚洲第一av免费看| 中文字幕人妻熟女乱码| 老熟妇乱子伦视频在线观看| 色综合站精品国产| √禁漫天堂资源中文www| 神马国产精品三级电影在线观看 | a级毛片a级免费在线| 日韩欧美国产在线观看| 亚洲国产精品sss在线观看| 午夜免费激情av| 久久精品亚洲精品国产色婷小说| 波多野结衣高清作品| 成熟少妇高潮喷水视频| 成人av一区二区三区在线看| 欧美成人性av电影在线观看| 最新美女视频免费是黄的| 亚洲美女黄片视频| 成年人黄色毛片网站| 亚洲自偷自拍图片 自拍| 亚洲av电影在线进入| 日韩大尺度精品在线看网址| 国产野战对白在线观看| 久久九九热精品免费| 国产成+人综合+亚洲专区| 国产av不卡久久| 国产高清videossex| 亚洲国产精品999在线| 久99久视频精品免费| 国产午夜福利久久久久久| 国产精品久久视频播放| 国产精品亚洲美女久久久| 久久久久久久午夜电影| 最近最新中文字幕大全电影3 | 黑丝袜美女国产一区| 成人一区二区视频在线观看| 99热只有精品国产| 在线播放国产精品三级| 男人舔女人下体高潮全视频| 在线av久久热| 国产爱豆传媒在线观看 | 欧美黑人巨大hd| 岛国在线观看网站| 中文字幕另类日韩欧美亚洲嫩草| 色综合婷婷激情| 久久久久久久午夜电影| 免费观看精品视频网站| 欧美绝顶高潮抽搐喷水| 国产一区二区激情短视频| 日韩大尺度精品在线看网址| 亚洲成av片中文字幕在线观看| 观看免费一级毛片| 一区福利在线观看| 久久这里只有精品19| 午夜福利高清视频| 午夜精品久久久久久毛片777| 欧美日韩黄片免| 午夜免费激情av| 精品少妇一区二区三区视频日本电影| 欧美zozozo另类| 亚洲国产精品成人综合色| 国产亚洲精品一区二区www| 国产精品日韩av在线免费观看| 久久青草综合色| 日韩三级视频一区二区三区| 我的亚洲天堂| 亚洲精品色激情综合| 伦理电影免费视频| 午夜a级毛片| 久久久久国产一级毛片高清牌| 亚洲精品在线观看二区| bbb黄色大片| 国产一区二区三区视频了| 午夜福利视频1000在线观看| 精品少妇一区二区三区视频日本电影| 亚洲aⅴ乱码一区二区在线播放 | 人人澡人人妻人| 男女床上黄色一级片免费看| 国产又黄又爽又无遮挡在线| 在线永久观看黄色视频| svipshipincom国产片| 国产成人系列免费观看| 国产国语露脸激情在线看| 妹子高潮喷水视频| 亚洲五月天丁香| 欧美 亚洲 国产 日韩一| 亚洲狠狠婷婷综合久久图片| 亚洲欧美日韩高清在线视频| 欧美午夜高清在线| 国产高清videossex| av在线播放免费不卡| 久久草成人影院| 神马国产精品三级电影在线观看 | 亚洲片人在线观看| 免费在线观看黄色视频的| 黄色成人免费大全| 在线观看午夜福利视频| 黑人巨大精品欧美一区二区mp4| 18美女黄网站色大片免费观看| 国产黄片美女视频| av超薄肉色丝袜交足视频| 精品少妇一区二区三区视频日本电影| 丝袜人妻中文字幕| 欧美三级亚洲精品| 亚洲五月婷婷丁香| 身体一侧抽搐| 中文字幕人妻熟女乱码| 亚洲专区中文字幕在线| 最近在线观看免费完整版| 国产精品久久久久久亚洲av鲁大| 亚洲成a人片在线一区二区| 国产亚洲精品第一综合不卡| 法律面前人人平等表现在哪些方面| 久热这里只有精品99| 91大片在线观看| 亚洲精品色激情综合| 日日爽夜夜爽网站| 免费高清在线观看日韩|