亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        泵站虹吸式出水流道駝峰排氣過程氣液兩相流研究

        2022-10-27 09:42:38劉霞黃可燦趙亮常正柏蔣勁
        中國農(nóng)村水利水電 2022年10期
        關(guān)鍵詞:虹吸式氣團(tuán)駝峰

        陳 奇,劉霞,黃可燦,趙亮,常正柏,蔣勁

        (1.湖北省水利水電科學(xué)研究院,湖北武漢 430072;2.湖北省國際灌排研究培訓(xùn)中心,湖北武漢 430072;3.武漢大學(xué)動(dòng)力與機(jī)械學(xué)院,湖北武漢 430072)

        0 引言

        泵站是指通過泵進(jìn)行抽水、增壓,再通過相關(guān)水工筑物進(jìn)行輸水、配水、排水的工程設(shè)施[1]。泵站在實(shí)際工程中應(yīng)用較為廣泛,在農(nóng)業(yè)灌溉、輸水調(diào)水、城市排澇等方面起著重要作用,為國家經(jīng)濟(jì)建設(shè)和人民的幸福生活做出了重要貢獻(xiàn)。出水流道作為泵站的關(guān)鍵過流部件之一,其對(duì)泵站的效率、出水流態(tài)以及泵站的穩(wěn)定運(yùn)行影響較大。出水流道的形式主要有虹吸式、直管式、貓背式、屈膝式等[2]。對(duì)于虹吸式出水流道而言,其頂部裝有真空破壞閥,在水泵停機(jī)時(shí)可用來破壞虹吸切斷水流,運(yùn)行可靠方便,且不會(huì)影響到防洪堤,工程造價(jià)較低,適用于出口水位變化不太大的泵站。

        目前已經(jīng)有諸多學(xué)者對(duì)虹吸式流道的相關(guān)問題進(jìn)行了研究。朱紅耕[3,4]采用了三維數(shù)值模擬方法分析了虹吸管內(nèi)部流道的流態(tài),發(fā)現(xiàn)在駝峰段和下降段水流速度變化較大時(shí),容易發(fā)生脫流現(xiàn)象。譚淋露[5]采用三維數(shù)值模擬方法分析了虹吸式出水管內(nèi)部水力特性,通過對(duì)比計(jì)算虹吸管不同的上升角、駝峰斷面高寬比、下降角和出口斷面高寬比對(duì)其水力特性的影響,確定了更優(yōu)參數(shù),降低了水力損失。模擬馮建剛[6]通過模型試驗(yàn),研究了虹吸形成的3 個(gè)主要階段,分別為水力驅(qū)氣階段、水力挾氣階段、穩(wěn)定虹吸階段。其中水力挾氣階段對(duì)形成穩(wěn)定虹吸所需時(shí)間影響較大。汪尚紅、李澤等[7]對(duì)虹吸式出水流道的軸流泵過渡過程進(jìn)行了研究,發(fā)現(xiàn)開機(jī)過程出水閘門的開啟時(shí)間不宜過長,以防止機(jī)組過載;關(guān)機(jī)過程出水閘門先快關(guān)后慢關(guān),且對(duì)于不同泵站開閘和關(guān)閘相關(guān)時(shí)間都有最優(yōu)值。王曉升、馮建剛等[8,9]采用RNG k-? 湍流模型以及歐拉多相流模型對(duì)虹吸式出水流道虹吸形成過程進(jìn)行了數(shù)值模擬計(jì)算,發(fā)現(xiàn)其他參數(shù)結(jié)構(gòu)固定的情況下,對(duì)水力挾氣階段影響最大的是流速,流速越大所需時(shí)間越短。M.Cihan 等[10]在數(shù)值模擬時(shí)使用VOF模型,模擬出了排水過程中空氣進(jìn)入虹吸管的現(xiàn)象,在其物理模型試驗(yàn)中也觀察到了與數(shù)值模擬類似的空氣夾帶現(xiàn)象,空氣進(jìn)入虹吸管時(shí)會(huì)降低排水效率并且引起振動(dòng)。Babaeyan-Koopaei等[11]通過模型試驗(yàn),分階段對(duì)特定的虹吸溢流道的進(jìn)氣口設(shè)計(jì)和流道結(jié)構(gòu)設(shè)計(jì)進(jìn)行改善,在對(duì)實(shí)驗(yàn)?zāi)P椭贿M(jìn)行簡單修改的情況下,發(fā)現(xiàn)在虹吸式流道特定高度的上方開一定大小和數(shù)量的水平方形槽孔,可以獲得最佳的空氣調(diào)節(jié)穩(wěn)定性,并且提高了溢洪道的容量。Rahim等[12]采用RNG k-?湍流模型和VOF模型計(jì)算出了虹吸式出水流道的流量系數(shù),并與模型試驗(yàn)相對(duì)比,數(shù)據(jù)表明一致性較高,吻合較好,驗(yàn)證了數(shù)值模擬的可行性。

        虹吸式出水流道形成穩(wěn)定虹吸需要經(jīng)歷4 個(gè)階段:啟動(dòng)階段,水力驅(qū)氣階段、水力挾氣階段、穩(wěn)定虹吸階段[8],其中水力挾氣階段最為關(guān)鍵,水流需要將駝峰內(nèi)的大氣團(tuán)帶走,才能進(jìn)入到穩(wěn)定階段。若將氣團(tuán)整體吹走,則該階段經(jīng)歷的時(shí)間短,但需要的流速較大,對(duì)于大型的虹吸式出水流道,整體將氣團(tuán)帶走所需要的流速更大;若設(shè)計(jì)流速較小,則挾氣階段的時(shí)間將增加,流速過小甚至不能進(jìn)入穩(wěn)定虹吸階段。處于該階段的時(shí)間過長,會(huì)導(dǎo)致駝峰處水力損失偏大,且水流狀態(tài)不穩(wěn)定,使泵的工況不穩(wěn)定,這將引起機(jī)組產(chǎn)生振動(dòng)、噪音等。若此時(shí)水泵工作在其特性曲線的馬鞍區(qū)內(nèi),則將產(chǎn)生更強(qiáng)烈的振動(dòng)和噪聲,處于該工況時(shí)間過長甚至有可能造成事故。因此,研究流速大小以及出水流道的自身結(jié)構(gòu)參數(shù)對(duì)虹吸形成階段的影響,以及形成穩(wěn)定虹吸所需時(shí)間,對(duì)于泵站安全高效穩(wěn)定運(yùn)行具有重大意義。

        1 瞬態(tài)數(shù)值模擬方法

        1.1 基本控制方程

        流體流動(dòng)遵循的三大基本定理為:能量守恒、質(zhì)量守恒、動(dòng)量守恒,這三大定理結(jié)合數(shù)學(xué)方法即可得到流體運(yùn)動(dòng)的基本控制方程。本文所研究的虹吸式出水流道中的虹吸形成過程的可壓縮性不明顯,可忽略能量方程,僅考慮連續(xù)性方程及動(dòng)量方程。

        (1)連續(xù)性方程如下:

        式中:ρ為流體密度;u為流體運(yùn)動(dòng)速度。

        (2)動(dòng)量方程如下:

        式中:fi為單位質(zhì)量流體所受的質(zhì)量力;σij為流體的應(yīng)力張量。

        1.2 VOF多相流模型及RNG k-?湍流模型

        1.2.1 VOF多相流模型

        VOF 模型通過求解流域內(nèi)的各相的體積分?jǐn)?shù)及各相的動(dòng)量方程來模擬兩種及以上的流體,求解出的速度場等屬性由各相共享,所有相在每個(gè)網(wǎng)格單元里的體積分?jǐn)?shù)之和均為1。VOF 模型可應(yīng)用于水中大氣泡運(yùn)動(dòng),含有自由液面的流動(dòng)等,故可用于模擬虹吸式出水流道內(nèi)的挾氣過程。

        在VOF 模型中通過求解各相的體積分?jǐn)?shù)連續(xù)方程來追蹤各相之間的界面,其方程如下所示:

        式中:下標(biāo)i表示模型中流體的相;ρ為密度;v為速度。

        1.2.2 RNG k-?湍流模型

        線性渦黏湍流模型是目前在實(shí)際工程湍流問題中使用最為廣泛的湍流模型,如k-?、k-ω 模型等。而k-?模型經(jīng)過不斷完善后,成為目前使用最廣泛的模型,包括Standardk-?、Realizablek-?模型和RNGk-?模型。

        在之前對(duì)虹吸式出水流道的挾氣過程研究中,已經(jīng)比較了Standardk-?、Realizablek-?和RNGk-?模型模擬結(jié)果與實(shí)驗(yàn)的相似程度,發(fā)現(xiàn)選用RNGk-?模型的仿真結(jié)果與實(shí)驗(yàn)有較高的相似性[13],且RNGk-?模型考慮到了湍流漩渦,對(duì)低雷諾數(shù)也適用,可信度和精度較高,故選用RNGk-?模型。

        RNGk-?模型不可壓縮流體的湍動(dòng)能k方程及耗散率ε方程如下所示:

        k方程:

        ε方程:

        式中:μeff是有效湍流黏度,μeff=μ+μt;常數(shù)項(xiàng)取值:αk=αε=0.139,=1.42,C2ε=1.68。

        1.3 幾何模型及邊界條件

        本文涉及到的模型有3個(gè),改變的參數(shù)為駝峰截面高寬比,分別為0.4和0.35,0.3,模型的駝峰截面尺寸均相等,面積與3 m直徑的圓管相同。為了更快速地帶走氣泡,下降段需要保持較高的流速,故下降段沒有采用擴(kuò)散,快到底部時(shí)采用擴(kuò)散管型以減小流速。

        3 個(gè)模型總體較為相似,正視圖如圖1 所示,入口截面高度為H1,駝峰截面高度為H2,駝峰頂距離水面H0,出水管淹沒在水面下的高度為H4,出口截面高度為H3,模型長度為L0,模型寬度為W0。模型具體尺寸如表1所示。

        圖1 虹吸式出水流道正模型正視圖Fig.1 Front view of the front model of the siphon outlet channel

        表1 模型幾何尺寸數(shù)據(jù)表Tab.1 Model geometry data

        3 個(gè)模型均采用速度入口和壓力出口,入口流速按照工況不同則設(shè)置不同,由駝峰截面平均流速作為衡量標(biāo)準(zhǔn),出口壓力由淹沒深度及出口高度共同確定,高寬比為0.4、0.35、0.3 的模型出口壓力分別為19 300、18 400、17 420 Pa,重力方向設(shè)置為豎直向下。

        1.4 網(wǎng)格劃分

        將高寬比為0.3,下降角為30度的模型分別按照0.03、0.04、0.05 m 的網(wǎng)格單元尺寸劃分網(wǎng)格,并在相同邊界條件與求解器設(shè)置下,計(jì)算至各項(xiàng)殘差均低于10-4,設(shè)立如圖2 所示的ABCD四個(gè)監(jiān)測點(diǎn)檢測流速。計(jì)算結(jié)果如表2所示。

        圖2 速度監(jiān)測點(diǎn)分布圖Fig.2 Distribution of velocity monitoring points

        從表2中的數(shù)據(jù)可以看出,A、B、C三點(diǎn)流速幾乎沒有差別,D 點(diǎn)流速0.03 m 和0.04 m 單元尺寸的網(wǎng)格較為接近,0.05 m 單元尺寸的網(wǎng)格與前兩者差異相對(duì)較大,故選用0.04 m 作為網(wǎng)格劃分的主要尺寸。

        表2 速度監(jiān)測點(diǎn)結(jié)果表Tab.2 Velocity monitoring point results

        采用網(wǎng)格單元尺寸為0.04 m 生成的網(wǎng)格如圖3 圖4 所示。圖3 為高寬比0.4 模型的網(wǎng)格,圖4 是高寬比0.3 模型僅僅縮短入口管道后的網(wǎng)格,縮短入口的目的在于減小網(wǎng)格數(shù)量,進(jìn)而縮短模擬破碎帶走的整個(gè)過程所需時(shí)間。高寬比為0.35 和0.3的模型網(wǎng)格與圖3所示網(wǎng)格類似。

        圖3 高寬比為0.4的模型網(wǎng)格圖Fig.3 Mesh diagram of the model with an aspect ratio of 0.4

        圖4 高寬比為0.3的模型縮短入口管段后的網(wǎng)格圖Fig.4 Mesh diagram of the model with an aspect ratio of 0.3 after shortening the inlet pipe

        1.5 求解設(shè)置

        VOF 模型比較適合用來計(jì)算氣團(tuán)的運(yùn)動(dòng),且在先前的工作中已經(jīng)確定RNG k-? 模型用來計(jì)算挾氣過程比較適合,故選用上述兩種模型來模擬氣團(tuán)帶走的過程。使用瞬態(tài)求解,速度入口的速度大小根據(jù)不同工況設(shè)置,壓力出口根據(jù)出口高度以及淹沒深度共同確定,上文已經(jīng)提及。管壁為無滑移邊界,粗糙度為1 mm,粗糙系數(shù)為0.5,選用SIMPLE算法,計(jì)算時(shí)間步長為0.01 s,殘差控制全部設(shè)置為10-4,各項(xiàng)殘差達(dá)到10-4以下視為收斂。

        2 數(shù)值模擬結(jié)果及分析

        2.1 短時(shí)間內(nèi)氣團(tuán)一次性排出計(jì)算

        該部分計(jì)算了3 個(gè)不同駝峰高寬比的虹吸式出水流道,當(dāng)駝峰頂部存在大氣團(tuán)時(shí),一次性帶走所需要的最小流速,其目的在于快速形成虹吸。氣團(tuán)初始狀態(tài)如圖5 所示,氣團(tuán)底部比駝峰底部略低。

        圖5 氣團(tuán)初始狀態(tài)圖Fig.5 Diagram of the initial state of the air masses

        2.1.1 高寬比為0.4的模型一次性排出氣團(tuán)所需最小流速

        高寬比為0.4 的虹吸式出水流道的計(jì)算結(jié)果如圖6 和圖7所示。設(shè)定駝峰截面流速為2.6 m/s,模擬至10 s 時(shí)的氣液分布及速度云圖如圖6 所示。從圖6(a)中可以看出,水流未能充滿駝峰頂部,氣泡將持續(xù)破碎并被帶向出口,該流速下不能一次性排出氣團(tuán)。

        圖6 駝峰截面平均流速為2.6 m/s時(shí)的結(jié)果圖Fig.6 The result diagram when the average velocity of the hump section is 2.6 m/s

        逐漸增大駝峰截面流速,重復(fù)計(jì)算,發(fā)現(xiàn)當(dāng)流速增大至4.3 m/s,駝峰頂部可以直接充滿水。如圖7 所示,模擬至6 s時(shí),水流可以充滿駝峰頂部,后續(xù)帶走氣團(tuán)會(huì)比較迅速,但此時(shí)駝峰截面流速為4.3 m/s,流速過大,實(shí)際工程較少采用。

        圖7 駝峰截面平均流速為4.3 m/s時(shí)的結(jié)果圖Fig.7 The result diagram when the average velocity of the hump section is 4.3 m/s

        2.1.2 高寬比為0.35的模型一次性排出氣團(tuán)所需最小流速

        為減小駝峰截面平均流速,將駝峰高寬比降低至0.35,計(jì)算結(jié)果如圖8 和圖9 所示。當(dāng)駝峰截面平均流速為3.5 m/s時(shí),模擬到6 s時(shí)發(fā)現(xiàn)水流未能充滿駝峰頂部,如圖8所示。當(dāng)駝峰截面流速增大到4.1 m/s時(shí),模擬至6 s時(shí)水流充滿整個(gè)駝峰,如圖9所示,后續(xù)氣團(tuán)排出會(huì)較為迅速。

        圖8 駝峰截面平均流速為3.5 m/s時(shí)的結(jié)果圖Fig.8 The result graph when the average velocity of the hump section is 3.5 m/s

        圖9 駝峰截面平均流速為4.1 m/s時(shí)的結(jié)果圖Fig.9 The result diagram when the average velocity of the hump section is 4.1 m/s

        2.1.3 高寬比為0.3的模型一次性排出氣團(tuán)所需最小流速

        為了進(jìn)一步減小駝峰截面平均流速,減小駝峰高寬比至0.3,計(jì)算結(jié)果如圖10 和圖11 所示。圖10 是駝峰截面流速為3.5 m/s,模擬至6 s時(shí)的氣液分布及速度云圖,從圖10(a)中可以看出,水流未能充滿駝峰頂部,需增大流速。

        圖10 駝峰截面平均流速為3.5 m/s時(shí)的結(jié)果圖Fig.10 Result diagram when the average velocity of the hump section is 3.5 m/s

        當(dāng)截面流速增大至3.8 m/s時(shí),模擬至6 s 時(shí)水流可以充滿駝峰頂部,如圖11(a)所示。

        綜上,高寬比為0.4、0.35、0.3 的模型一次性排出的流速分別為4.3、4.1、3.8 m/s,降低高寬比后,一次性排出氣團(tuán)所需要的臨界流速有所減小,但此流速仍顯偏大,所以大口徑的虹吸式出水流道想要短時(shí)間一次性整體將氣團(tuán)從駝峰頂部排出需要較大的流速。

        2.2 氣團(tuán)破碎帶走計(jì)算

        從上述結(jié)果來看,駝峰截面高寬比越小,一次性排出氣團(tuán)所需流速越小,故選用高寬比為0.3 的模型來進(jìn)行氣泡破碎帶走全過程模擬,有利于進(jìn)一步降低流速。本節(jié)模擬出了120 s和180 s 排氣時(shí)間內(nèi)將氣團(tuán)全部排出所對(duì)應(yīng)的駝峰平均流速以及排氣全過程。

        2.2.1 排氣時(shí)間為120 s的氣團(tuán)排出過程

        當(dāng)排氣時(shí)間為120 s時(shí),氣團(tuán)破碎排出過程及形成穩(wěn)定虹吸時(shí)的速度云圖分別如圖12 和圖13 所示。圖12 所示的(a)~(f)的6 個(gè)狀態(tài)分別為0.3、8、34、50、80、120 s 時(shí)的氣液分布狀態(tài),圖12(f)所示的狀態(tài)為最終的虹吸形成狀態(tài),氣泡完全排出,整個(gè)過程耗時(shí)120 s,此時(shí)駝峰截面平均流速為3.0 m/s。從圖13 的速度云圖可以看出,形成穩(wěn)定虹吸后,駝峰截面在豎直方向有較大的速度梯度。

        圖12 駝峰截面平均流速為3.0 m/s時(shí)的排氣過程Fig.12 The exhaust process when the average velocity of the hump section is 3.0 m/s

        圖13 駝峰截面平均流速為3.0 m/s時(shí)形成虹吸后的速度云圖Fig.13 Velocity diagram after siphon formation when the average flow velocity of the hump section is 3.0 m/s

        2.2.2 排氣時(shí)間為180 s的氣團(tuán)排出過程

        當(dāng)排氣時(shí)間為180 s時(shí),氣團(tuán)破碎排氣過程及形成穩(wěn)定虹吸時(shí)的速度云圖分別如圖14 和圖15 所示。圖14 所示的(a)~(g)的6個(gè)狀態(tài)分別為0.3、8、34、50、80、120、180 s時(shí)的氣液分布狀態(tài),圖14(g)所示的狀態(tài)為最終的虹吸形成狀態(tài),氣泡幾乎完全排出,整個(gè)過程耗時(shí)180 s,此時(shí)駝峰截面平均流速為2.5 m/s。相較于120 s 的排氣時(shí)間,駝峰截面流速從3.0 m/s 降低到了2.5 m/s,降低了16.67%。從圖14(e)和圖14(f)可以看出,排出駝峰頂部的少量氣體所耗的時(shí)間較長,即駝峰頂部的氣體相較于下降段內(nèi)的氣體更難排出,結(jié)合速度云圖可以看出,由于存在垂直方向上的速度梯度,駝峰頂部流速較小,這是駝峰頂部氣體更難排出的主要原因。

        圖14 駝峰截面平均流速為2.5 m/s時(shí)的排氣過程Fig.14 The exhaust process when the average velocity of the hump section is 2.5 m/s

        圖15 駝峰截面平均流速為2.5 m/s時(shí)形成虹吸后的速度云圖Fig.15 Velocity diagram after siphon formation when the average flow velocity of the hump section is 2.5 m/s

        2.3 排氣過程不同時(shí)間點(diǎn)湍流動(dòng)能分布

        對(duì)駝峰截面流速為2.5 m/s 的排氣過程中t=10 s 和t=80 s 時(shí)刻進(jìn)行湍流動(dòng)能進(jìn)行分析,其分布分別如圖16、17所示。

        圖16 t=10 s時(shí)的氣液分布及湍流動(dòng)能分布Fig.16 Gas-liquid distribution and turbulent kinetic energy distribution at t=10 s

        從上面的結(jié)果可以看出,t=10 s時(shí),下降段中氣液交界面附近的湍流動(dòng)能較大,氣團(tuán)破碎迅速;而當(dāng)t=80 s時(shí),駝峰段中氣液交界面附近的湍流動(dòng)能大大減小,氣團(tuán)破碎速度變慢。這說明湍流動(dòng)能的大小會(huì)影響氣團(tuán)的破碎速度。

        3 結(jié)論

        圖17 t=80 s時(shí)的氣液分布及湍流動(dòng)能Fig.17 Gas-liquid distribution and turbulent kinetic energy distribution at t=80 s

        基于CFD 研究了以兩種方式帶走大口徑矩形斷面駝峰頂部大氣團(tuán),分別是一次性快速排出氣團(tuán)和一定時(shí)間內(nèi)破碎排出氣團(tuán),其中前一種方式計(jì)算了0.4、0.35、0.3 三種不同高寬比的虹吸管模型,而后一種方式采用0.3 高寬比的模型模擬了兩種不同流速下氣團(tuán)運(yùn)動(dòng)、破碎帶走、形成虹吸的整個(gè)過程,結(jié)論如下。

        (1)降低高寬比可以降低一次性排出氣團(tuán)需要流速va。模型駝峰高寬比為0.4、0.35、0.3時(shí),對(duì)應(yīng)的va分別為4.3、4.1、3.8 m/s,0.3 高寬比的模型與0.4 高寬比的模型相比,va降低了11.6%。雖然有所降低,但流速仍然偏大,所以較大口徑的虹吸式出水流道不建議采用一次性排出氣團(tuán)的方式。

        (2)允許排氣時(shí)間越長,對(duì)應(yīng)的駝峰截面流速越小,駝峰頂部的氣體相比于下降段中的更難排出。當(dāng)排氣時(shí)間tb為120 s時(shí),對(duì)應(yīng)的駝峰截面流速vb為3.0 m/s,當(dāng)tb為180 s時(shí),對(duì)應(yīng)vb為2.5 m/s,流速降低了16.67%。由于駝峰頂部存在垂直方向上的速度梯度,頂部流速較小,且排氣過程中駝峰段的湍流動(dòng)能較小,使得駝峰頂部氣體比下降段中的氣體更難排出。

        猜你喜歡
        虹吸式氣團(tuán)駝峰
        虹吸式雨水排水技術(shù)在建筑給排水施工中的應(yīng)用
        泗陽二站虹吸式出水流道水力優(yōu)化設(shè)計(jì)
        輸水管線啟動(dòng)填充過程含滯留氣團(tuán)瞬變流數(shù)值模擬
        噴射虹吸式坐便器水包裝置試驗(yàn)與結(jié)構(gòu)優(yōu)化分析
        基于后向軌跡模式的合肥市臭氧來源分析
        駝峰場編發(fā)線向到發(fā)場交授權(quán)電路探討
        駝峰第三制動(dòng)位減速器夾停鉤車問題的改進(jìn)
        氣團(tuán)來源對(duì)沿海城市PM2.5中二次水溶性無機(jī)離子形成特征的影響研究
        HZ330CL橫折缺陷原因分析與改進(jìn)
        淺析虹吸式雨水排水系統(tǒng)與重力式雨水排水系統(tǒng)
        成人免费丝袜美腿视频| 美腿丝袜网址亚洲av| 午夜理论片yy44880影院| 无码人妻精品一区二区三区免费| 一区二区黄色素人黄色| 亚洲成a人片在线观看无码| 蜜桃传媒网站在线观看| 国模精品二区| 三年片免费观看影视大全视频 | 99国产精品欲av麻豆在线观看| 亚洲av日韩专区在线观看| 日韩精品人妻视频一区二区三区 | 国产精品白浆免费观看| 亚洲精品国产综合一线久久| 日本精品一区二区三区二人码| 国产午夜激无码AV毛片不卡| 在线观看精品视频网站| 开心五月激情五月天天五月五月天| 久久精品久久精品中文字幕| 丰满熟妇乱又伦精品| 国产在线观看不卡网址| 无码少妇一区二区三区 | 大尺度无遮挡激烈床震网站 | 中文字幕有码一区二区三区| 天天爽夜夜爽人人爽一区二区| 手机在线看片国产人妻| 亚州毛色毛片免费观看| 男女啪啪无遮挡免费网站| 亚洲精品国产成人久久av| 国产xxxxx在线观看免费| 大地资源在线观看官网第三页| 中文乱码字幕精品高清国产| 亚洲va中文字幕欧美不卡| 大陆极品少妇内射aaaaa| 国产桃色一区二区三区| 天天草夜夜草| 国产精品欧美福利久久| 久久精品亚洲精品国产区| 野外三级国产在线观看| 特级做a爰片毛片免费看无码| 人人人妻人人澡人人爽欧美一区 |