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

        ?

        膜態(tài)沸騰球體水下運動減阻特性

        2019-03-13 07:11:20李佳川魏英杰王聰夏維學
        兵工學報 2019年2期
        關(guān)鍵詞:水蒸汽球體壁面

        李佳川, 魏英杰, 王聰, 夏維學

        (哈爾濱工業(yè)大學 航天學院, 黑龍江 哈爾濱 150001)

        0 引言

        減小液體中運動物體的阻力是流體動力學中的一個重要問題。較低的阻力意味著更快的速度和更少的能源消耗。減小阻力的方法有很多,例如應用聚合物添加劑[1]、超空泡[2-3]、微氣泡注射[4-5]、疏水涂層[6-8]等。這些方法主要是通過改變物體的表面特性或通過改變物體附近流體的密度來實現(xiàn)減阻功能。

        當物體的表面溫度遠高于液體沸點時,周圍液體汽化后的蒸汽不斷匯聚,形成連續(xù)的蒸汽膜,稱為膜態(tài)沸騰過程。由于蒸汽膜的絕熱性會使熱量交換變得很差,物體壁面會達到很高的溫度,甚至有燒毀的危險,在熱工程領(lǐng)域中,一般要避免膜態(tài)沸騰的發(fā)生。最近的研究表明,利用膜態(tài)沸騰形成的蒸汽膜能夠?qū)崿F(xiàn)減阻功能[9-13],因此膜態(tài)沸騰減阻成為一種較新的減阻方法。

        膜態(tài)沸騰產(chǎn)生的蒸汽膜具有滑移效應,類似于嵌入超疏水性壁面的空氣墊。然而超疏水性壁面的減阻效果有限,且在水壓及水流沖刷作用下空氣層很快會消失。持續(xù)加熱物體可以在壁面處形成連續(xù)穩(wěn)定的蒸汽膜,減小摩擦阻力。Saranadhi等[9]研究表明,當雷諾數(shù)在26 100~52 000范圍內(nèi)時,摩擦阻力可以減小80%~90%,減阻效果非常明顯。

        膜態(tài)沸騰生成的蒸汽膜,不僅可以減小摩擦阻力,還可以減小壓差阻力。Vakarelski等[10]采用一種沸點和汽化潛熱均低于水的氟化液體,對膜態(tài)沸騰階段球體在該液體中的自由下落過程進行了研究,結(jié)果表明蒸汽膜的存在可以降低球體受到的壓差阻力,使總阻力減小85%. 隨后,Vakarelski等[11-12]采用不同液體對膜態(tài)沸騰球體繞流過程進行了詳細的實驗研究,分析了液體溫度、球體溫度和雷諾數(shù)等對減阻特性的影響。

        目前,對膜態(tài)沸騰球體繞流特性的研究還處于初步階段,相關(guān)減阻機理還需要進一步研究。由于膜態(tài)沸騰球體繞流過程涉及相變、湍動、流動分離和多相流動等復雜物理問題,至今還沒有學者對該過程進行數(shù)值仿真研究。Gruncell等[14]對超疏水性球體的繞流過程進行了數(shù)值仿真研究,通過固定一層厚度相等且不變形的空氣層等效疏水特性。對于膜態(tài)沸騰球體繞流過程,不同表面位置蒸汽膜的厚度不同,且存在蒸汽堆積現(xiàn)象,因此用Gruncell等[14]的方法進行等效模擬的準確性較低。

        本文采用Mixture多相流模型、剪切壓力傳輸(SST)k-ω湍流模型,耦合蒸發(fā)- 冷凝相變模型,對膜態(tài)沸騰球體繞流過程進行數(shù)值仿真研究。研究膜態(tài)沸騰球體繞流特性,分析雷諾數(shù)對阻力特性的影響,并對減阻機理進行研究。

        1 數(shù)值計算方法

        1.1 控制方程

        為了識別球體表面蒸汽膜,球體壁面附近需要較密的網(wǎng)格,使得整體網(wǎng)格量較大,因此本文將球體繞流問題簡化為二維軸對稱流動問題,總體網(wǎng)格量約為44萬。

        采用Mixture多相流模型描述各相流體體積分數(shù)。水和水蒸汽的混合密度為

        ρm=ρl(1-αv)+ρvαv,

        (1)

        式中:ρl、ρv分別為水的密度和水蒸汽的密度;αv為水蒸汽的體積分數(shù)。

        質(zhì)量守恒方程為

        (2)

        式中:t為時間;ui(i= 1, 2)為平面笛卡爾坐標系下的速度分量。

        動量守恒方程為

        (3)

        式中:uj(j=1, 2)為平面笛卡爾坐標系下的速度分量;p為壓力;g為重力加速度;F為表面張力;τij為黏性剪切應力,

        (4)

        δij為Kronecker符號,μm為混合黏度;uk(k=1,2)為平面笛卡爾坐標系下的速度分量。

        能量守恒方程為

        (5)

        式中:cp為定壓比熱;T為溫度;k為導熱系數(shù);Φ為相變潛熱產(chǎn)生的能量源項,

        (6)

        Lv為汽化潛熱,fv=ρvαv/ρm為汽相質(zhì)量分數(shù)。

        水汽化成水蒸汽的相間傳遞過程采用蒸發(fā)- 冷凝模型進行模擬,其輸運方程表達式為

        (7)

        (8)

        (9)

        cflv與cfvl分別為蒸發(fā)和冷凝系數(shù),Ts為水的飽和溫度。

        球體繞流問題,需要較精確地預測流動分離位置,因此本文湍流模型選用可以較好地模擬分離問題的SSTk-ω模型。

        1.2 計算域及網(wǎng)格劃分

        計算域、邊界條件和局部網(wǎng)格劃分如圖1所示。本文計算域采用二維軸對稱模型,球體直徑為D,計算域徑向直徑為10D,高度為15D. 球體壁面采用無滑移壁面邊界條件,側(cè)邊壁采用滑移壁面邊界條件,下邊界為速度入口,上邊界為壓力出口。球體位置固定不變,重力方向向下,來流方向向上,入口速度u為定值,可以等效模擬球體在水中勻速自由下落的過程。本文網(wǎng)格采用C-Block劃分策略,為了識別蒸汽膜,球體表面附近網(wǎng)格進行了加密,球體壁面附近第1層網(wǎng)格高度設置為1 μm,滿足SSTk-ω湍流模型所要求的y+<1,壁面附近網(wǎng)格高度增長因子設置為1.05,C-Block外層網(wǎng)格高度增長因子設置為1.1.

        圖1 計算域及局部網(wǎng)格劃分Fig.1 Computational domain and local mesh generation

        1.3 數(shù)值方法驗證

        目前僅有少數(shù)學者針對膜態(tài)沸騰球體水下運動過程進行了實驗研究,本文將數(shù)值仿真得到的膜態(tài)沸騰球體繞流運動阻力系數(shù)與Vakarelski等[11]的膜態(tài)沸騰球體水下自由下落實驗得到的阻力系數(shù)進行對比,以驗證本文數(shù)值方法的準確性。球體水中自由下落過程中,受到水的阻力與運動速度的平方呈正比,球體下落的初始階段速度較小,球體做向下加速運動,當球體速度達到一定值時,水的阻力與球體重力相等,球體將做勻速運動。

        球體水中運動的阻力系數(shù)為

        Cd=2Fd/(πR2ρlu2),

        (10)

        式中:Fd為球體受到的阻力;R為球體半徑;u為球體水中運動速度。

        當球體勻速運動時阻力與重力相等時,阻力系數(shù)可以表達為

        Cd=[8g(ρ-ρl)R]/(3ρlu2),

        (11)

        式中:ρ為球體密度。

        根據(jù)Vakarelski等[11]得到的球體勻速運動速度,利用(11)式計算出球體的阻力系數(shù)。

        如圖1所示,本文數(shù)值仿真球體位置固定,水相對于球體向上運動,水流速度為u,重力方向向下,可以等效模擬Vakarelski等[11]的球體在水中勻速下落過程。表1為本文數(shù)值仿真與Vakarelski等[11]實驗得到的膜態(tài)沸騰球體水中勻速運動阻力系數(shù)對比,表中雷諾數(shù)Re=ρluD/μl,μl為水的黏度。

        表1 數(shù)值與實驗阻力系數(shù)對比

        從表1中可以看到,本文數(shù)值仿真選用的球體溫度與Vakarelski等[11]實驗選用的球體溫度不同。Vakarelski等[11]采用超疏水性球體,相對于親水性球體,超疏水性球體可以在較低溫度下形成膜態(tài)沸騰[11,15]。由于同時考慮流動、傳熱、相變和表面潤濕性比較復雜,不易于模擬,本文數(shù)值仿真沒有考慮球體表面潤濕性,為了形成膜態(tài)沸騰,球體溫度設置為450 ℃的較高溫度。Vakarelski等[10-12]的研究表明,球體發(fā)生膜態(tài)沸騰后會繼續(xù)升高球體溫度,球體達到勻速時的速度u變化較小,可以忽略,由(11)式可知,膜態(tài)沸騰球體的阻力系數(shù)隨球體溫度變化較小,可以忽略。因此,可以采用表1所示的不同球體溫度與Vakarelski等[11]的實驗結(jié)果進行對比。在Vakarelski等[11]的研究中,球體水中下落時間小于1.5 s,且由于蒸汽具有較低的導熱系數(shù),球體溫度降低較小。由于沸騰球體繞流物理過程復雜,本文數(shù)值仿真進行了簡化,僅對定常狀態(tài)下的膜態(tài)沸騰球體繞流進行研究,球體溫度恒為450 ℃.

        從表1中可以看到,3種不同Re下的數(shù)值仿真Cd值與Vakarelski等[11]實驗的Cd值最大相差-12.5%. 由于本文數(shù)值計算簡化為二維軸對稱問題,忽略了尾流不對稱性對阻力的影響,會與實際流動過程略有偏差。但阻力系數(shù)范圍與實驗范圍一致,并隨Re的增加而大幅度減小。因此可以采用本文數(shù)值仿真方法對膜態(tài)沸騰球體在水中的運動過程進行研究,并分析相應的減阻機理。

        2 計算結(jié)果及分析

        2.1 膜態(tài)沸騰球體繞流特性

        本文對比分析Re=1.28×105條件下膜態(tài)沸騰球體與普通球體的繞流特性。其中水溫均為95 ℃,膜態(tài)沸騰球體球溫為450 ℃,普通球體球溫為95 ℃. 兩種球體的球體直徑D為20 mm,水流速度u=2 m/s. 計算所得的普通球體阻力系數(shù)為0.48,計算所得的膜態(tài)沸騰球體阻力系數(shù)為0.22(見表1)。

        當普通球體流動處于亞臨界Re范圍內(nèi),即5×103≤Re≤3×105時,阻力系數(shù)約為常數(shù),Cd≈0.44. 本文計算所得的普通球體阻力系數(shù)在0.44附近,相差約為9.1%. 膜態(tài)沸騰球體的阻力系數(shù)約為普通球體的50%,由此可知膜態(tài)沸騰可以大幅度減小球體受到的阻力。

        為了便于分析,對球體表面位置進行規(guī)定,如圖2所示,水流方向向上,球面底端位置P對應角度θ=0°,球面中端位置M對應θ=90°,球面頂端位置N對應θ=180°.

        圖2 球體表面位置示意Fig.2 Schematic diagram of surface position of sphere

        圖3(a)所示為膜態(tài)沸騰球體繞流的水蒸汽云圖,圖3(b)所示為球體周圍溫度云圖。當球體表面溫度遠高于水的沸點時,汽化成的水蒸汽匯聚到一起,使球體被一層水蒸汽膜包裹,如圖3(a)中M點處的水蒸汽云圖放大部分所示。球體尾部由于旋渦與回流的影響,沿流動方向生成汽泡云。由于球體周圍溫度大部分較低,且溫度范圍為95~450 ℃,跨度較大,圖3(b)僅顯示95~130 ℃部分的溫度云圖,其中溫度大于130 ℃部分顯示的顏色與130 ℃溫度一致。由于水蒸汽的導熱系數(shù)遠小于水,蒸汽膜的存在使高溫球體與周圍水隔開,降低了球體的傳熱能力,溫度場中大于130 ℃分布在球體表面附近的蒸汽膜內(nèi),球體尾部的汽泡云溫度在130 ℃以下。

        圖3 膜態(tài)沸騰球體繞流水蒸汽與溫度云圖Fig.3 Water vapor phase and temperature distribution

        圖4所示為普通球體與膜態(tài)沸騰球體的繞流速度v矢量圖。由圖4可見:球體側(cè)面位置附近流速較高,尾部旋渦區(qū)速度較低;由于蒸汽膜的潤滑作用,膜態(tài)沸騰球體繞流速度最大值要大于普通球體;膜態(tài)沸騰球體尾部低速區(qū)比普通球體窄。

        圖4 普通球體與膜態(tài)沸騰球體繞流速度矢量圖Fig.4 Velocity fields of ordinary sphere and film boiling sphere

        圖5所示為普通球體與膜態(tài)沸騰球體的繞流流線圖。由圖5可見,球體在水中受到的阻力分為摩擦阻力與壓差阻力,摩擦阻力由水的粘滯性引起,壓差阻力由球體上表面和下表面的壓力差產(chǎn)生。對于普通球體,當5×103≤Re≤3×105時,摩擦阻力占總阻力的比例較小,球體主要受壓差阻力的作用。針對球體繞流問題,壓差阻力主要由球體表面發(fā)生流動分離的位置決定,流動分離發(fā)生的越遲緩、越靠近球體尾部,球體繞流受到的壓差阻力就越小。對于普通球體,球體表面流動分離發(fā)生在M點稍前的位置,即分離角度略小于90°[16],與圖5(a)所示的結(jié)果相吻合。

        圖5 普通球體與膜態(tài)沸騰球體繞流流線圖Fig.5 Streamline distributions of flows around ordinary sphere and film boiling sphere

        Vakarelski等[10]對膜態(tài)沸騰球體在氟化液體中的自由下落過程進行了研究,發(fā)現(xiàn)膜態(tài)沸騰球體達到勻速時的速度要大于普通球體。結(jié)果表明,膜態(tài)沸騰階段產(chǎn)生的蒸汽膜使流動分離位置向尾部移動,減小了球體受到的壓差阻力,從而減小了球體受到的總阻力。與Vakarelski等[10]的研究相吻合,圖5(b)所示膜態(tài)沸騰球體繞流,分離位置位于球體表面的尾部區(qū)域。相比于普通球體,膜態(tài)沸騰球體尾部旋渦區(qū)較窄,尾部形狀趨于流線型。流線形狀與圖3所示水蒸汽相與溫度相云圖尾部形狀一致。本文的數(shù)值方法可以較好地模擬膜態(tài)沸騰球體繞流問題。

        2.2 不同Re膜態(tài)沸騰球體繞流特性

        對于普通球體,阻力系數(shù)大小主要與Re有關(guān)。當Re很小(Re<1)時,黏性力起主要作用,阻力系數(shù)可以表達為Cd=24/Re. 隨著Re增大,Cd不存在解析解,Cd隨著Re增大而減小。當5×103≤Re≤3×105時,流動分離點幾乎不發(fā)生改變,Cd幾乎與Re無關(guān)。當Re達到臨界值,即Re≈3.5×105時,邊界層轉(zhuǎn)化為湍流邊界層,使流動分離位置向后移動,阻力系數(shù)減小為Cd≈0.1. 不同于普通球體,Vakarelski等[10-12]的實驗表明,在亞臨界流動范圍內(nèi),即5×103≤Re≤3×105時,膜態(tài)沸騰球體阻力系數(shù)隨Re增大而減小。

        本文采用數(shù)值仿真方法研究亞臨界流動范圍內(nèi)Re對膜態(tài)沸騰球體阻力系數(shù)Cd的影響。Re分別取值5.06×104、1.28×105和2.24×105. 如表1所示,得到的阻力系數(shù)分別為0.38、0.22和0.14,與Vakarelski等[11]的實驗結(jié)果相吻合。這一結(jié)果表明,隨著Re增大,膜態(tài)沸騰球體繞流運動的阻力系數(shù)逐漸減小。在2.1節(jié)已經(jīng)對比分析了Re為1.28×105時普通球體與膜態(tài)沸騰球體的繞流特性,該節(jié)主要分析Re分別取值5.06×104和2.24×105時的膜態(tài)沸騰球體繞流特性,其中球體直徑均為20 mm,水流速度分別為0.79 m/s和3.50 m/s.

        圖6~圖8所示分別為兩種Re條件下膜態(tài)沸騰球體繞流的水蒸汽相圖、速度矢量圖和流線圖。由圖6~圖8可以看出,Re越大,球體尾部水蒸汽云越趨于流線型,尾部低速旋渦區(qū)越窄,流動分離位置越靠近球體尾部。

        圖6 不同Re膜態(tài)沸騰球體繞流水蒸汽相圖Fig.6 Water vapor phase distributions of flows around film boiling spheres with different Re

        圖7 不同Re膜態(tài)沸騰球體繞流速度矢量圖Fig.7 Velocity fields of flows around film boiling spheres with different Re

        圖8 不同Re膜態(tài)沸騰球體繞流流線圖Fig.8 Streamline distributions of flows around film boiling spheres with different Re

        為了清晰對比不同雷諾數(shù)膜態(tài)沸騰球體繞流的流動分離特性,圖9給出了不同雷諾數(shù)膜態(tài)沸騰球體及常溫球體繞流近壁面流線對比圖。由于在亞臨界流動范圍內(nèi),常溫球體繞流的流動分離位置幾乎保持一致,圖9僅呈現(xiàn)Re=1.28×105的普通球體近壁面流線。分析圖9可知,普通球體繞流的流動分離位置為θ≈85°,略大于Achenbach[16]實驗所得的θ≈82.5°.Re取值分別為5.06×104、1.28×105和2.24×105時膜態(tài)沸騰球體繞流的流動分離位置分別為θ≈99°、θ≈108°和θ≈120°.膜態(tài)沸騰球體繞流與普通球體相比,流動分離位置偏向尾部,且隨著Re增大,分離角度逐漸增大。在亞臨界流動范圍內(nèi),阻力大小主要由壓差阻力決定,即由流動分離位置決定。因此,膜態(tài)沸騰繞流受到的阻力要小于普通球體,且隨著Re增大,減阻效果越明顯。

        圖9 近壁面流線對比Fig.9 Comparison of streamlines near the sphere surface

        2.3 膜態(tài)沸騰球體繞流運動減阻機理

        流動分離的條件為逆壓梯度,流動分離發(fā)生的根本原因為流體黏性的存在。在分離區(qū)附近,流體的動能不僅逐漸轉(zhuǎn)變?yōu)閴簭娔?,同時邊界層內(nèi)的粘性流體(特別在臨近壁面處)受到壁面的滯止作用(表現(xiàn)為切應力),不斷消耗動能。當靠近壁面處流體的動能全部耗盡,流體微團將在分離點處停滯不前,并發(fā)生堆積。分離點后在反向壓差作用下,壁面附近產(chǎn)生回流,排擠上游流體使流動在分離點處發(fā)生分離。Vakarelski等[10-11]的實驗研究指出,膜態(tài)沸騰球體被一層蒸汽膜包裹,減小了切應力對流體的滯止作用,使流動分離點向后移動。

        圖10所示為Re=1.28×105時普通球體與膜態(tài)沸騰球體繞流在圖2所示M點附近的速度矢量對比圖。由圖10可見,普通球體壁面附近速度梯度較小,距離壁面較遠位置速度才達到外部主流速度。由于水的黏度遠大于水蒸汽,膜態(tài)沸騰球體在蒸汽膜內(nèi),速度梯度較大,在距離壁面較近的蒸汽膜外部具有較大的流動速度。且膜態(tài)沸騰球體外層速度與來流速度平行,普通球體速度稍偏離來流方向偏向遠離球面的方向。圖11給出了M點處沿y軸方向流體的x軸方向無量綱速度v(v=vx/u,vx為x軸方向的流動速度)對比,其中來流方向為x軸負方向。球體直徑為20 mm,M點處y軸坐標值為10 mm(位置與方向規(guī)定見圖2)。由圖11可知:普通球體M點位于流動分離點之后近壁面位置發(fā)生回流,因此普通球體近壁面處x軸方向速度為正值(與來流方向相反);蒸汽膜的存在,使膜態(tài)沸騰球體沿y軸方向流體速度增長較快,且達到的速度最大值略大于普通球體。

        圖10 普通球體與膜態(tài)沸騰球體在M點附近速度矢量Fig.10 Velocity vectors of ordinary sphere and film boiling sphere near point M

        圖11 M點處沿y軸方向流體的x軸方向無量綱速度對比Fig.11 Comparison of fluid velocities in the direction of x at point M

        將蒸汽膜與球體看作整體,膜態(tài)沸騰球體的繞流運動可以等效為流體繞蒸汽膜的運動,則壁面處的無滑移邊界條件可以轉(zhuǎn)換為蒸汽膜處的滑移邊界條件。蒸汽膜邊界處黏滯力的減小使膜態(tài)沸騰球體繞流流動分離點向球體尾部移動,減小了球體受到的阻力。

        由2.2節(jié)可知,Re對膜態(tài)沸騰球體繞流的阻力特性有較大影響。Gruncell等[14]通過在球體表面固定一層厚度相等且不變形的氣體層,對超疏水球體的繞流過程進行了數(shù)值仿真研究,結(jié)果表明氣體層的滑移效應使流動分離點向球體尾部移動,極大地減小了球體繞流受到的阻力。對于膜態(tài)沸騰球體的繞流過程,球體表面的蒸汽膜厚度并不處處相等。圖12所示為不同Re膜態(tài)沸騰球體繞流蒸汽相體積分數(shù)為0.5的等值線圖。由圖12可知:蒸汽會在逆壓梯度作用下發(fā)生堆積,使蒸汽膜外部流動發(fā)生分離;類似于普通球體,隨著Re增大,堆積位置會向球體尾部移動;蒸汽膜的減阻作用,會由于蒸汽的堆積而受到抑制。因此膜態(tài)沸騰球體繞流過程中球體的阻力系數(shù)隨著Re增大而減小。

        圖12 不同Re水蒸汽相體積分數(shù)為0.5的等值線圖Fig.12 Contour plot of 0.5 vapor fraction for different Re

        3 結(jié)論

        本文采用數(shù)值仿真方法對亞臨界流動范圍內(nèi)膜態(tài)沸騰球體的繞流特性進行了研究,分析了Re對阻力特性的影響,并對膜態(tài)沸騰減阻機理進行了研究。主要得到以下結(jié)論:

        1)本文數(shù)值仿真方法可以較好地模擬膜態(tài)沸騰球體繞流特性,得到的阻力系數(shù)與已有文獻[11]的一致性較好。

        2)由于水蒸汽的黏度遠小于水,壁面附近蒸汽膜內(nèi)部的速度梯度較大,相比于普通球體,膜態(tài)沸騰球體繞流在距離壁面較近的水蒸汽膜外部具有較大的流動速度,水蒸汽膜的滑移效應減小了壁面對流體的粘滯力作用,使流動分離點向尾部移動,使球體受到的阻力減小。

        3)由于逆壓梯度的存在,膜態(tài)沸騰球體生成的水蒸汽會在球體壁面尾部區(qū)域堆積,隨著Re增大,堆積位置向后移動,延緩流動分離,阻力系數(shù)逐漸減小。

        猜你喜歡
        水蒸汽球體壁面
        二維有限長度柔性壁面上T-S波演化的數(shù)值研究
        計算機生成均值隨機點推理三、四維球體公式和表面積公式
        消費電子(2020年5期)2020-12-28 06:58:27
        汽車擋風玻璃上起霧的原因及消除方法的分析
        物理教師(2018年5期)2018-06-14 04:57:34
        變壓法真空干燥鐵芯生銹和干燥時間應對策略
        卷宗(2018年9期)2018-06-07 08:58:16
        廣告創(chuàng)意新方法——球體思維兩極法
        微波輔助水蒸汽結(jié)合上浮溶劑固化法快速萃取谷物中三嗪類除草劑
        分析化學(2018年12期)2018-01-22 12:31:46
        鑄鐵排氣歧管表面氧化處理的研究
        壁面溫度對微型內(nèi)燃機燃燒特性的影響
        Optimization of rice wine fermentation process based on the simultaneous saccharification and fermentation kinetic model☆
        顆?!诿媾鲎步Ec數(shù)據(jù)處理
        影视av久久久噜噜噜噜噜三级| 亚洲中文字幕精品一区二区| 免费人成黄页网站在线观看国产| 久久精见国产亚洲av高清热| av深夜福利在线| 国产人成在线免费视频| 九七青青草视频在线观看| 亚洲中字幕日产av片在线| 久久aⅴ人妻少妇嫩草影院| 日韩精品无码一区二区中文字幕| 久久天天躁狠狠躁夜夜96流白浆| 亚洲va在线va天堂va四虎| 久久亚洲国产高清av一级 | 国产一品二品精品在线| 国产成人a在线观看视频免费 | 精品国产一级毛片大全| 一本久久a久久精品综合| 日韩精品在线视频一二三| 少妇私密会所按摩到高潮呻吟| 国产精品自在线免费| 精品亚洲视频免费观看网站 | 在教室伦流澡到高潮hgl动漫| 曰韩亚洲av人人夜夜澡人人爽| 无码中文字幕在线DVD| 青青草原亚洲在线视频| 日本高清乱码中文字幕| 中文字幕+乱码+中文字幕一区| 中文亚洲爆乳av无码专区| 亚洲综合新区一区二区| 日本边添边摸边做边爱喷水 | 高清高速无码一区二区| 午夜日本理论片最新片| 日本视频一区二区三区观看| 乱中年女人伦av一区二区 | 日韩精品久久久一区| 国产成人精品aaaa视频一区| 国产亚洲av综合人人澡精品| 人妻仑乱a级毛片免费看| 中文字幕日本最新乱码视频| 激情亚洲的在线观看| 人妻少妇艳情视频中文字幕|