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

        ?

        孔板對(duì)水槽內(nèi)流場(chǎng)影響規(guī)律

        2018-01-06 00:32:56林日億李曉辰徐偉棟
        關(guān)鍵詞:孔板水槽射流

        林日億, 郭 彬, 李曉辰, 徐偉棟, 唐 振

        (1.中國(guó)石油大學(xué)儲(chǔ)運(yùn)與建筑工程學(xué)院,山東青島 266580; 2.中國(guó)石油大學(xué)山東省油氣儲(chǔ)運(yùn)安全省級(jí)重點(diǎn)實(shí)驗(yàn)室,山東青島 266580)

        孔板對(duì)水槽內(nèi)流場(chǎng)影響規(guī)律

        林日億1,2, 郭 彬1, 李曉辰1, 徐偉棟1, 唐 振1

        (1.中國(guó)石油大學(xué)儲(chǔ)運(yùn)與建筑工程學(xué)院,山東青島 266580; 2.中國(guó)石油大學(xué)山東省油氣儲(chǔ)運(yùn)安全省級(jí)重點(diǎn)實(shí)驗(yàn)室,山東青島 266580)

        為了在室內(nèi)試驗(yàn)臺(tái)上模擬海洋海流橫向沖刷海洋立管,采用水槽內(nèi)加孔板結(jié)構(gòu)模擬海洋平流。為獲取穩(wěn)定、均勻的平流場(chǎng),利用FLUENT軟件進(jìn)行水槽內(nèi)孔板均流流場(chǎng)數(shù)值模擬,選用標(biāo)準(zhǔn)κ-ε紊流模型分析孔板孔徑、厚度、位置及布置層數(shù)等因素對(duì)孔板后流場(chǎng)的影響,并優(yōu)化孔板結(jié)構(gòu)。結(jié)果表明:未安裝孔板時(shí),在整個(gè)水槽內(nèi)出現(xiàn)射流且有擴(kuò)散的趨勢(shì),但射流并不能帶動(dòng)水槽內(nèi)全部流體,僅有部分流體有速度,在水槽上、下兩端會(huì)有漩渦及回流產(chǎn)生,速度分布不均勻。布置3塊厚度為0.2 m的孔板且板與板之間的間隔為0.3 m,按照第一塊板距離進(jìn)口截面0.2 m布置且中間孔孔徑設(shè)計(jì)為0.03 m時(shí)能夠獲得較為合理的均勻流場(chǎng)。

        孔板; 數(shù)值模擬; 計(jì)算流體力學(xué); 均勻流場(chǎng)

        近年來(lái),對(duì)深海石油和天然氣的開采與日俱增,海洋立管作為深海石油開采中不可或缺的一個(gè)部件,因海流沖刷而設(shè)計(jì)和建設(shè)技術(shù)難度大[1]。為研究海洋立管流動(dòng)規(guī)律,試驗(yàn)室內(nèi)用均流孔板模擬海洋流場(chǎng),通過(guò)水槽研究孔板對(duì)流場(chǎng)的影響規(guī)律。對(duì)于流場(chǎng)的研究,數(shù)值模擬相比于實(shí)驗(yàn)研究來(lái)說(shuō)具有經(jīng)濟(jì)、高效的特點(diǎn)[2],尤其利用計(jì)算流體力學(xué)(CFD)通用軟件后,避免了復(fù)雜的編程工作,已成為數(shù)值模擬研究流場(chǎng)的有力工具[3]??装逅魇堑湫偷姆蛛x流動(dòng),流場(chǎng)高度不均勻,且具有很高的紊動(dòng)性,分區(qū)特性也很明顯[4]。在孔口和靠近壁面位置,會(huì)伴隨有漩渦和回流的出現(xiàn)[5]。朱孟府等[6]和過(guò)小玲等[7]選用標(biāo)準(zhǔn)κ-ε湍流模型對(duì)孔板的水力特性進(jìn)行了研究,Yakhot等[8]從理論上導(dǎo)出了κ-ε模型,并用理論分析的方法確定了κ-ε模型中的參數(shù)。進(jìn)行數(shù)值模擬時(shí),對(duì)標(biāo)準(zhǔn)κ-ε、RNG、RSM三種湍流模型進(jìn)行比較,表明標(biāo)準(zhǔn)κ-ε湍流模型是更合適的湍流模型[9],Singh等[10]以及Werth等[11]則分別用數(shù)值模擬和實(shí)驗(yàn)的方法對(duì)孔板流量計(jì)進(jìn)行了研究。綜上所述,國(guó)內(nèi)外的學(xué)者對(duì)孔板的研究主要集中在湍流模型的選擇、孔板流量計(jì)以及孔板的水力特性方面,對(duì)孔板的均流作用研究甚少。筆者在室內(nèi)試驗(yàn)臺(tái)基礎(chǔ)上,通過(guò)FLUENT軟件開展水槽內(nèi)孔板均流流場(chǎng)的數(shù)值模擬,選用標(biāo)準(zhǔn)κ-ε紊流模型,使用SIMPLEC算法進(jìn)行計(jì)算,分別對(duì)孔板孔徑、厚度、位置及布置層數(shù)等影響因素進(jìn)行模擬,分析這幾種因素對(duì)流場(chǎng)的影響,優(yōu)化孔板結(jié)構(gòu)。

        1 物理數(shù)學(xué)模型

        1.1 物理模型

        模擬海水平流水槽內(nèi)充滿水,水自左側(cè)口流入,經(jīng)水槽從右側(cè)口流出(圖1)。綜合孔板位置、孔徑、板數(shù)等因素的影響,在水槽內(nèi)進(jìn)行孔板布置,最終使流經(jīng)最后一塊孔板后的流場(chǎng)達(dá)到均勻。

        圖1 水槽模型Fig.1 Flume model

        計(jì)算水槽模型如圖1所示,其主要尺寸:長(zhǎng)L=5 m,高H=1.3 m,厚δ=0.3 m,進(jìn)出口直徑d=0.2 m,進(jìn)出口截面為a=0.3 m的正方形截面,頂面長(zhǎng)L′=4 m,模型上下左右前后均對(duì)稱。

        1.2 數(shù)學(xué)模型

        1.2.1 控制方程和湍流模型

        對(duì)于所有的方程,均用FLUENT解質(zhì)量守恒和動(dòng)量守恒方程[12],計(jì)算控制方程選用連續(xù)性方程和動(dòng)量方程,湍流模型選擇標(biāo)準(zhǔn)k-ε雙變量湍流模型,因其適用性強(qiáng),求解穩(wěn)定,精度也能滿足要求。壁面函數(shù)的選擇保持為標(biāo)準(zhǔn)k-ε湍流模型默認(rèn)的壁面函數(shù),即標(biāo)準(zhǔn)壁面函數(shù),流動(dòng)為定常流動(dòng),考慮重力。

        (1) 質(zhì)量守恒方程。質(zhì)量守恒方程又稱為連續(xù)性方程,其形式為

        (1)

        式中,ρ為流體的密度,kg/m3;ui為流體的速度,m/s。

        該方程是質(zhì)量守恒方程的一般形式,適用于可壓縮流動(dòng)和不可壓縮流動(dòng)。源項(xiàng)Sm是從分散相中加入到連續(xù)相的質(zhì)量, 也可以是任何的自定義源項(xiàng)[13]。

        (2) 動(dòng)量守恒方程。在慣性坐標(biāo)系下,i方向的動(dòng)量守恒方程為

        (2)

        式中,p為靜壓力,Pa;τij為應(yīng)力張量,kg/(m·s2);ρgi和Fi分別為i方向上的重力體積力和外部體積力,kg/(m2·s2)。

        Fi包含了其他的模型相關(guān)源項(xiàng),如多孔介質(zhì)和自定義源項(xiàng)[13]。

        1.2.2 控制方程離散

        選用SIMPLEC算法進(jìn)行計(jì)算,SIMPLEC算法收斂速度快且計(jì)算結(jié)果準(zhǔn)確,使用Green-Gauss Node Based方法進(jìn)行迭代,相比于Green-Gauss Cell Based求解更精確且能夠最小化偽擴(kuò)散,壓力的離散采用Standard格式,適用性廣,其余均采用二階迎風(fēng)格式,求解精度高。

        2 模 擬

        2.1 計(jì)算域和計(jì)算網(wǎng)格

        采用的計(jì)算域?yàn)槿鐖D1所示的水槽模型,計(jì)算域網(wǎng)格圖以及局部加密圖如圖2所示。

        2.2 邊界條件

        入口邊界為速度入口,大小設(shè)置為1.5 m/s,給定湍流強(qiáng)度和水力直徑,分別為3.3%和0.2 m,出口邊界設(shè)為壓力出口,出口壓力設(shè)置為大氣壓。

        2.3 網(wǎng)格獨(dú)立性檢驗(yàn)

        以進(jìn)出口壓差為考核依據(jù)對(duì)網(wǎng)格進(jìn)行獨(dú)立性考核,如圖3所示。

        圖2 模型網(wǎng)格圖Fig.2 Model grid graph

        圖3 網(wǎng)格獨(dú)立性考核Fig.3 Grid independence assessment

        由圖3可以看出,隨網(wǎng)格數(shù)增加,進(jìn)出口壓差減小;當(dāng)網(wǎng)格數(shù)為240萬(wàn)時(shí),曲線趨于平穩(wěn),之后隨網(wǎng)格數(shù)增加,進(jìn)出口壓差穩(wěn)定,因此選取240萬(wàn)網(wǎng)格進(jìn)行計(jì)算。

        3 結(jié)果分析與討論

        3.1 單層孔板時(shí)流場(chǎng)的影響因素

        3.1.1 位置布置

        針對(duì)第一塊孔板選取3種不同的布置位置,距離入口分別為20、30、40 cm,對(duì)應(yīng)的速度沿高度分布如圖4所示。

        由圖4可以看出,在高度為-1.26 ~-1.0 m內(nèi),3種情況均產(chǎn)生回流,且集中在底部,板后回流區(qū)產(chǎn)生的位置與文獻(xiàn)[14]中結(jié)果一致。對(duì)比3條曲線,距離進(jìn)口20 cm時(shí),進(jìn)口射流擴(kuò)散范圍廣,回流速度較其余兩種情況小,因此20 cm的情況更優(yōu)。

        3.1.2 孔 徑

        針對(duì)第一塊板中間位置孔分別選取3種不同孔徑:2、3和4 cm。3種情況下速度沿高度方向分布如圖5所示。

        由圖5可以看出,孔徑為2 cm時(shí),水槽兩端射流流速大,中間流速小,說(shuō)明2 cm的孔徑較小,射流主要從水槽兩端流過(guò);孔徑為3 cm時(shí),中間流速大,水槽底部有回流產(chǎn)生,說(shuō)明在回流區(qū)對(duì)應(yīng)位置產(chǎn)生了低壓區(qū)[15];孔徑為4 cm時(shí),速度峰值增大,但峰值較窄,上下兩端有回流產(chǎn)生,因此,孔徑為3 cm時(shí)最優(yōu)。

        圖4 不同布置位置下高度方向速度分布Fig.4 Velocity distribution of height direction under different arrangement positions

        圖5 不同孔徑下高度方向速度分布Fig.5 Velocity distribution of height direction under different aperture

        3.1.3 厚 度

        不同厚度時(shí)高度方向速度分布如圖6所示。

        圖6 不同厚度時(shí)高度方向速度分布Fig.6 Velocity distribution of height direction under different

        由圖6可以看出,不同厚度間差別較小,底部均產(chǎn)生回流,速度峰值及高度稍有差別,相比來(lái)說(shuō),厚度為20 cm時(shí)回流量小,速度峰值的位置也靠近中間,因此選擇厚度為20 cm的孔板。

        3.2 多層孔板時(shí)流場(chǎng)的影響因素

        3.2.1 布置層數(shù)

        不加孔板和加孔板后的流場(chǎng)分布如圖7所示。

        由圖7可以看出,流場(chǎng)均勻性方面,加入孔板相比于未加入孔板,進(jìn)口射流得到擴(kuò)散,隨布置層數(shù)的增多,進(jìn)口射流進(jìn)一步擴(kuò)散至整個(gè)水槽,加人3塊孔板后,水槽內(nèi)流場(chǎng)均勻化。圖7(a)中,上下兩端流體速度小,接近于零,進(jìn)口射流在黏性力作用下速度會(huì)逐漸減小,在靠近出口處,因出口截面收縮,速度會(huì)急劇增大;進(jìn)口射流速度大,兩端速度小,速度不均導(dǎo)致水槽內(nèi)產(chǎn)生漩渦,進(jìn)而在上下兩端靠近壁面處有回流產(chǎn)生;圖7(b)中,加入第一塊孔板,核心區(qū)自孔口開始,逐漸收縮,主流區(qū)內(nèi)徑向流速梯度接近為零,剪切應(yīng)力也很小[14];加入兩層孔板后,射流擴(kuò)散至整個(gè)水槽,僅有個(gè)別位置速度較小,回流區(qū)消失,符合前述文獻(xiàn)研究:隨孔板數(shù)增加,回流區(qū)減少[14]。加入3塊孔板后,速度分布均勻化,滿足設(shè)計(jì)要求。

        圖7 XOY截面速度分布云圖Fig.7 Velocity distribution cloud map of XOY section

        沿長(zhǎng)度方向速度變化曲線如圖8所示。由圖8可以看出,在接近出口位置,速度均有一個(gè)急劇回升。不加孔板時(shí),軸線位置速度大,變化劇烈;加入孔板后,速度變化平穩(wěn)。加入兩塊孔板,相較于一塊孔板速度更加穩(wěn)定,但在前期速度會(huì)出現(xiàn)波動(dòng),綜上所述,加入3塊孔板是最為理想的布置方式。

        圖8 沿長(zhǎng)度方向速度變化曲線Fig.8 Curve of velocity along length direction

        3.2.2 布置方式

        針對(duì)第3塊孔板選取兩種不同的分布方式:順排和叉排。兩種情況下速度分布如圖9所示。

        圖9 不同分布方式時(shí)高度方向速度分布Fig.9 Height direction velocity distribution in different distribution modes

        由圖9可以看出,兩種布置方式下速度沿高度方向均有波動(dòng),從波動(dòng)幅度上來(lái)說(shuō),順排布置時(shí),最大速度約為0.14 m/s,最小速度為0.1 m/s;叉排布置時(shí),最大速度約為0.176 m/s,最小速度約為0.098 m/s,因此順排的布置方式波動(dòng)幅度相對(duì)較小,且從整個(gè)高度范圍內(nèi)的速度分布來(lái)看,順排布置下流場(chǎng)更加穩(wěn)定。

        3.2.3 孔眼結(jié)構(gòu)

        3塊孔板孔眼結(jié)構(gòu)布置如圖10所示。由圖10可以看出,第一塊孔板孔眼分布有兩個(gè)特點(diǎn):一是中間分布稀,兩端分布密;二是中間孔徑小,兩端孔徑大,目的是將進(jìn)口射流打散。射流經(jīng)過(guò)第一塊孔板后,由圖7(b)所示,速度仍集中于中間位置,加入第二塊孔板,孔眼結(jié)構(gòu)見圖10(b),目的將射流進(jìn)一步擴(kuò)散至整個(gè)水槽,流向兩端。經(jīng)過(guò)兩塊孔板后的流場(chǎng)如圖7(c)所示,僅有個(gè)別位置不均勻;為進(jìn)一步使流場(chǎng)均勻,加入如圖10(c)所示的第三塊孔板,加入3塊孔板后,流場(chǎng)均勻。

        圖10 三塊孔板孔眼結(jié)構(gòu)布置Fig.10 Three-plate hole structure layout

        4 結(jié) 論

        (1)未安裝孔板時(shí),在整個(gè)水槽內(nèi)出現(xiàn)射流且有擴(kuò)散的趨勢(shì),但射流并不能帶動(dòng)水槽內(nèi)全部流體,僅有部分流體有速度,在水槽上、下兩端會(huì)有漩渦及回流產(chǎn)生,速度分布不均勻。

        (2)相比于位置、厚度等因素,孔徑對(duì)流場(chǎng)的影響最明顯。

        (3)孔板結(jié)構(gòu)優(yōu)化后可實(shí)現(xiàn)最佳均勻平流,孔板結(jié)構(gòu)參數(shù)為:布置3塊間隔0.3 m、厚度為0.2 m的孔板,第一塊板距離進(jìn)口截面0.2 m,前兩塊板中間孔孔徑設(shè)置為0.03 m。

        [1] DENG Y, HUANG W,ZHAO J. Combined action of uniform flow and oscillating flow around marine riser at low Keulegan-Carpenter number[J]. Ocean University, 2014,13(3):390-396.

        [2] SALAH-ADDIN B, AL-OMARI K K, et al. Soot processes in a methane-fueled furnace and their impact on radiation heat transfer to furnace walls [J].International Journal of Heat and Mass Transfer, 2001,44(13):2567-2581.

        [3] 于娟.揮發(fā)分、CO 火焰與碳粒燃燒的相互作用及其?;痆D]. 上海:上海交通大學(xué),2002.

        YU Juan. Study and modelling on the interaction of volatile flame,CO flame and char particle combustion[D]. Shanghai:Shanghai Jiaotong University,2002.

        [4] 才君梅,馬俊,張子冀.孔板流場(chǎng)的二維激光測(cè)速試驗(yàn)研究[J]. 水力發(fā)電學(xué)報(bào), 1999,4:52-58.

        CAI Junmei, MA Jun, ZHANG Ziji. An experimental research on the flow field of orifice plate by using the 2-dimension LDV system[J]. Journal of Hydroelectric Engineering,1999,4:52-58.

        [5] 楊永剛.方形孔口多孔板水力空化反應(yīng)器的試驗(yàn)研究與數(shù)值模擬[D]. 杭州:浙江工業(yè)大學(xué),2012.

        YANG Yonggang. Experimental study and numerical simulation of hydrodynamic cavitations reactor with multi-square-hole orifice plates[D]. Hangzhou: Zhejiang University of Technology,2012.

        [6] 朱孟府,苗秀娟,鄧橙,等.不同孔分布孔板的水力空化效果的數(shù)值模擬[J].輕工機(jī)械,2012,30(4):8-12.

        ZHU Mengfu, MIAO Xiujuan, DENG Cheng, et al. Numerical simulation of hydrodynamic cavitations effect of orifice plate with different layout[J]. Light Industry Machinery,2012,30(4):8-12.

        [7] 過(guò)小玲,金保升,沈丹,等.裝有多孔板的脫硫噴淋塔流場(chǎng)數(shù)值模擬研究[J]. 鍋爐技術(shù),2007,38(6):5-9.

        GUO Xiaoling, JIN Baosheng, SHEN Dan, et al. Flow simulation for FGD spray scrubber with porous plate[J]. Boiler Technology,2007,38(6):5-9.

        [8] YAKHOT V, ORZAG S. Renormalization group analysis of turbulence: basic theory[J].Journal of Scientific Computing, 1986,1(1):3-51.

        [9] AI W, ZHOU Q J. Hydraulic characteristics of multi-stage orifice plate[J]. Journal of Shanghai Jiaotong University(Science), 2014,19(3):361-366.

        [10] SINGH V K, THARAKAN T J. Numerical simulations for multi-hole orifice flow meter [J]. Flow Measurement and Instrumentation, 2015(45):375-383.

        [11] WERTH D E, KHAN A A, GREGG W B. Experimental study of wall curvature and bypass flow effects on orifice discharge coefficients [J]. Experiments in Fluids, 2005,39(3):483-489.

        [12] 韓占忠,王敬,蘭小平.FLUENT流體工程仿真計(jì)算實(shí)例與應(yīng)用[M]. 北京:北京理工大學(xué)出版社,2004.

        [13] 楊肖曦,李松巖,林日億,等.泡沫流體攜砂能力的數(shù)值模擬[J]. 中國(guó)石油大學(xué)學(xué)報(bào)(自然科學(xué)版),2006,30(3):89-96.

        YANG Xiaoxi, LI Songyan, LIN Riyi, et al. Numerical simulation for prop-carrying capacity of foam fluid[J]. Journal of China University of Petroleum (Edition of Natural Science),2006,30(3):89-96.

        [14] 鐘偉.流量測(cè)量節(jié)流裝置的流場(chǎng)數(shù)值模擬與實(shí)驗(yàn)研究[D]. 南京:南京航空航天大學(xué), 2007.

        ZHONG Wei. Numerical simulation and experimental study on measuring the flow field of throttling device[D]. Nanjing:Nanjing University of Aeronautics & Astronautics, 2007.

        [15] 郭婷婷,李雪梅,格日勒,等.電站擋風(fēng)抑塵板后流場(chǎng)特性的數(shù)值研究[J].電站系統(tǒng)工程,2009,25(2):7-10.

        GUO Tingting, LI Xuemei, GE Rile, et al. Numerical simulation on flow characteristics of porous wind fences[J].Power System Engineering,2009,25(2):7-10.

        Influencelawsofporeplatesonbasininternalflowfield

        LIN Riyi1,2, GUO Bin1, LI Xiaochen1, XU Weidong1, TANG Zhen1

        (1.CollegeofPipeline&CivilEngineeringinChinaUniversityofPetroleum,Qingdao266580,China;2.ShandongProvincialKeyLaboratoryofOil&GasStorageandTransportationSafety,ChinaUniversityofPetroleum,Qingdao266580,China)

        In order to simulate the situation how the vertical tube is washed by the seawater inside the tube, this study is aimed to simulate the ocean advection by using the pore plate structure in a flume.To obtain the stable and uniform velocity field, the FLUENT software is used to simulate the stable velocity field by using the pore plate structure in a flume. The standardκ-εturbulent model is adopted to analyze the pore diameter, thickness, position and arrangement layers of the pore plate, then the structure of the pore plate is optimized. The simulation results show that, when the orifice plate is not installed, there is a trend of jet and diffusion in the whole tank, but the jet does not drive all the fluid in the tank. Only part of the fluid has the velocity, in the water tank, the lower end has whirlpool and reflux generation, and the velocity distribution is not uniform. It is more favorable to obtain the uniform flow field with three pore plates whose thickness is 0.2 m, where the interval between them is 0.3 m,the first pore plate is set to 0.2 m away from the inlet, and the diameter of the hole on the first pore plate is set to 0.03 m.

        pore plate; numerical simulation; computational fluid mechanics; uniform flow field

        2016-09-30

        “十三五”國(guó)家科技重大專項(xiàng)(2016ZX05012-002-005);中央高校基本科研業(yè)務(wù)費(fèi)專項(xiàng)(15CX05002A)

        林日億(1973-),男,教授,博士,研究方向?yàn)闊崃Σ捎秃蜔崮芾谩-mail:linry@upc.edu.cn。

        1673-5005(2017)06-0154-06

        10.3969/j.issn.1673-5005.2017.06.019

        TE 952

        A

        林日億,郭彬,李曉辰,等. 孔板對(duì)水槽內(nèi)流場(chǎng)影響規(guī)律[J]. 中國(guó)石油大學(xué)學(xué)報(bào)(自然科學(xué)版), 2017,41(6):154-159.

        LIN Riyi, GUO Bin, LI Xiaochen, et al. Influence laws of pore plates on basin internal flow field [J]. Journal of China University of Petroleum (Edition of Natural Science), 2017,41(6):154-159.

        (編輯 沈玉英)

        猜你喜歡
        孔板水槽射流
        核電廠高壓安注系統(tǒng)再循環(huán)管線節(jié)流孔板的分析與改進(jìn)
        深海逃逸艙射流注水均壓過(guò)程仿真分析
        低壓天然氣泄漏射流擴(kuò)散特性研究
        煤氣與熱力(2022年4期)2022-05-23 12:45:00
        可升降折疊的飲水機(jī)水槽
        限流孔板的計(jì)算與應(yīng)用
        廣州化工(2020年6期)2020-04-18 03:30:20
        為什么水槽管要做成彎曲狀
        長(zhǎng)距離礦漿管道系統(tǒng)中消能孔板的運(yùn)行優(yōu)化
        要挑好水槽,就看這里了!
        幸福(2016年6期)2016-12-01 03:08:13
        廚房水槽設(shè)計(jì)
        射流齒形噴嘴射流流場(chǎng)與氣動(dòng)聲學(xué)分析
        亚洲国产美女高潮久久久| 人妻少妇不满足中文字幕| 欧美精品aaa久久久影院| 国产精品性一区二区三区| 日本av一级视频在线观看| 一本色道久久88加勒比一| 国精品人妻无码一区免费视频电影| 台湾佬自拍偷区亚洲综合| 538任你爽精品视频国产| 中文字幕成人乱码亚洲| 国产精品女主播在线播放| 中文字幕 亚洲精品 第1页| 东北女人毛多水多牲交视频| 亚洲精品夜夜夜| 精品久久久亚洲中文字幕| 日韩av中文字幕波多野九色| 精品福利一区二区三区蜜桃| 99久久精品免费观看国产| 欧美性性性性性色大片免费的| 国产美女裸身网站免费观看视频| 少妇人妻偷人精品无码视频| 亚洲AV秘 无套一区二区三区| 97中文乱码字幕在线| 白白白在线视频免费播放| 欧美又大又色又爽aaaa片| 精品久久久久久久中文字幕| 中国人妻沙发上喷白将av| 亚洲av偷拍一区二区三区| 不卡免费在线亚洲av| 亚洲av无码日韩av无码网站冲| 人妻少妇久久中文字幕一区二区 | 刺激一区仑乱| 少妇被猛男粗大的猛进出| 国产成人精品曰本亚洲| 亚洲精品成人久久av| 日本熟妇另类一区二区三区| 国产综合无码一区二区辣椒| 黑人巨大videos极度另类| 国产激情视频在线观看首页| 国产一区亚洲一区二区| 少妇被粗大进猛进出处故事|