999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

改進Schnerr-Sauer模型在水翼空化模擬中的評估分析

2016-10-11 07:58:57洪鋒袁建平周幫倫石洋
哈爾濱工程大學學報 2016年7期
關鍵詞:模型

洪鋒,袁建平,周幫倫,石洋

(1.江蘇大學流體機體機械工程技術研究中心,江蘇鎮江212013;2.埃因霍分理工大學機械工程系,埃因霍分5612AZ荷蘭)

改進Schnerr-Sauer模型在水翼空化模擬中的評估分析

洪鋒1,2,袁建平1,周幫倫1,石洋1

(1.江蘇大學流體機體機械工程技術研究中心,江蘇鎮江212013;2.埃因霍分理工大學機械工程系,埃因霍分5612AZ荷蘭)

為了提高Schnerr-Sauer模型模擬空化流動的能力,基于Rayleigh-Plesset方程和均相流假設建立了一種修正的Schnerr-Sauer空化模型,修正后的模型考慮了湍動能引起的脈動壓力和不可凝結氣體的影響。聯立渦粘度系數修正的SST k-ω湍流模型,對繞二維Clark-Y翼型空化流進行了非定常數值計算,得到了時均化升阻力系數隨空化數的變化曲線、云狀空化時空泡周期性演變及其水平速度分布規律。通過與已有文獻的實驗結果對比分析得到,采用修正后的空化模型計算得到的時均化升阻力系數與實驗值更為接近,并能較準確地捕捉升阻力系數隨空化數的特殊變化規律;修正后的Schnerr-Sauer模型能準確地模擬翼型表面空化云的初生、成長、脫落和潰滅的全過程。本文的研究結果驗證了修正后的Schnerr-Sauer模型在水翼空化流數值計算中的可靠性和準確性。

空化流;水翼;空化模型;湍流模型;數值計算;阻力系數;Schnerr-Sauer模型

空化是水動力學中特有的一種物理現象,研究空化現象對改善水力機械和船艦推進器性能具有很現實的工程意義,發展空泡流動的數值模擬技術是目前空化研究領域中的熱點[1]。

當前,在工程實際應用中主要使用以下三種類型的空化模型來模擬空化流場:1)基于簡化Rayleigh-plesset方程的空化模型(如Singhal等[2]模型、Zwart-Gerber-Belamri[3]模型及 Schnerr-Sauer[4]模型);2)基于經驗得出的比例化空化模型(如Merk-le[5]模型、Kunz[6]模型);3)通過空泡受力分析得出的空化模型(如Tamura[7]模型)。上述空化模型中,基于簡化Rayleigh-Plesset方程的空化模型是在混合多相流模型的不可壓縮連續性方程中導出了空化反應速率,由于工程上解決空化問題的需要,ANSYS FLUENT v6.0開始引入Zwart-Gerber-Belamri模型和Schnerr-Sauer模型,并成功模擬了水翼的空化問題[8-12]。然而,文獻[13-14]中指出,基于簡化Rayleigh-plesset方程的Schnerr-Sauer空化模型在模擬云狀空化時存在不足,因此為了改善該類空化模型模擬空化流動的能力,本文基于Schnerr-Sauer模型的思想及均相流假設,推導出一種考慮不可凝結氣體及湍流脈動壓力影響的修正Schnerr-Sauer模型;然后應用UDF(User-defined Function),把修正后的Schnerr-Sauer模型加載到ANSYS FLUENT v14.5平臺并聯立渦黏度系數修正的SST k-ω湍流模型對二維Clark-Y翼型的空化流場進行數值計算。通過與原始的Schnerr-Sauer模型模擬結果及已有文獻中的實驗結果作對比,驗證修正后的Schnerr-Sauer模型在非定常空化流動計算中的可行性及準確性。

1 Schnerr-Sauer空化模型改進

參照文獻[2]將空化流還原成液態水,空泡及不可凝結氣體(non-condensable gas,NCG)的三相混合物,并且重新定義混合流體密度為

式中:αv和ρv分別為空泡體積分數和密度,αg和fg分別為NCG體積分數和質量分數,ρl為液態水的密度。NCG的體積分數與質量分數對應關系為

原始的Schnerr-Sauer空化模型蒸發與凝結項分別為

為了考慮不可凝結氣體的影響,將式(3)、(4)修正為

式中:αg為NCG體積分數。為了體現湍流脈動對空話初生的影響,將空化壓力進行修正為

式中:pturb為湍動能引起的脈動壓力;psat為飽和蒸汽壓,取25℃時水的飽和汽化壓力psat=3 169 Pa;ρm為混合介質的密度,k為湍動能。

式(5)、(6)就是修正后的Schnerr-Sauer空化模型空化源項,與文獻[15]中修正后Schnerr-Sauer模型相比,本文的修正方法還體現了NCG對蒸發及凝結過程的影響。確定不可凝結氣體的質量分數和體積分數分別為fg=1×10-6和αg=7.8×10-4。fg直接影響混合密度的改變,進而決定修正粘度μt和pv的改變,從而影響空化初生。αg直接決定空化速率S的改變,進而影響翼型的空化范圍[16]。

2 湍流模型及局部可壓縮性修正

由于最初的RANS(或URANS)湍流模型是在單相、完全不可壓縮介質的基礎上發展起來的,因而不能直接應用在可壓縮的兩相流動數值計算中。在使用RANS模型模擬空化流動時會過渡預測空穴尾部的湍流黏度,導致在空穴尾部區域產生了比實際大很多的粘滯力并阻礙流場中回射流結構因能量不足而無法向上游運動。因此,在使用RANS模型模擬這類空化流場時就無法準確模擬物理現象的本質特征,從而使得模擬失真。為了提高對水翼表面空化的預測精度,根據Reboud等[17]的思想,考慮到氣液混合區局部可壓縮性的影響,本文采用SST k-ω并對湍流黏度μT作如下相應的修正:

式中:n取值為3,ρm為混合流體密度,k為湍動能,ω為耗散率,其余變量與原始SST k-ω湍流模型保持一致。

3 計算網格及邊界條件

計算采用弦長c=70 mm的Clark-Y型水翼,圖1給出了計算區域及其邊界條件。計算區域的入口距翼型前緣為4c,出口距翼型尾緣的距離為5c。為了保證壁面函數對無量綱y+值的要求,對翼型周圍近壁面區域進行了網格加密,如圖2所示。

數值計算邊界條件為速度進口、壓力出口,上下邊界均為自由無滑移壁面條件,翼型表面采用絕熱、無滑移固壁條件。為了保證與實驗條件一致,設定進口速度為Uin=10 m/s,出口壓力大小根據空化數取得。翼型攻角α=8°,由弦長、來流速度及介質的黏度系數,得到對應的雷諾數Re=7×105。

圖1 計算域及其邊界條件Fig.1 Computational domain and boundary conditions

圖2 計算網格Fig.2 Computational meshes

本文對計算域劃分了4套不同尺寸的網格并作網格無關性驗證,網格節點數分別為N1=69 226,N2=83 208,N3=105 588及N4=123 188。網格無關性驗證結果如表1所示,從表1中可以看出,當網格數大于N3時,兩種不同空化數條件下水翼的升阻力系數相對變化均小于1%,表明此時網格數量對計算結果的影響可以忽略。同時,網格數為N3時,翼型的升阻力系數與試驗值誤差均在5.5%以內。圖3、4分別為當空化數σ=3.0時通過計算得到的翼型表面壓力系數(-Cp)及無量綱y+值,從圖中可以看出,壓力系數計算值與實驗值吻合較好且翼型表面網格的y+大部分在30~90,滿足壁面函數對近壁面區域網格的要求。因此,最終選擇網格3作為本文的計算網格。

表1 網格無關性驗證Table 1 Mesh independence check

圖3 水翼表面壓力系數分布Fig.3 Surface pressure coefficient of hydrofoil

圖4 水翼表面y+分布Fig.4 Wall y+distribution of hydrofoil

4 計算結果及分析

4.1時均化升阻力系數

為研究在不同程度空化條件下,采用修正后的Schnerr-Sauer模型及渦黏度系數修正的SST k-ω湍流模型的聯合求解能力。本文通過數值計算的方法得到了原始Schnerr-Sauer模型及修正后的Schnerr-Sauer模型在多種空化數條件下水翼的升阻力系數,并與已有文獻中的實驗結果作對比分析,結果如圖5所示。由于空化是一種包含相變的瞬時物理現象,本文所有計算均為非定常,時間步長設置為Δt= 1×10-4s,總迭代步數設置為2 000,取最后8個周期內升阻力系數的平均值作為時均化后的計算結果。定義空化數σ,升力系數Cl及阻力系數Cd分別為

式中:pout為出口壓力,pv為空化壓力,ρl為水的密度,Uin為來流速度,c為翼型弦長。

分析升力系數可以看出,在無空化條件下,由于介質中空化核的存在以及不可凝結氣體的作用,導致二者計算得到的升力系數并不相同。采用修正后的Schnerr-Sauer模型得到的升力系數與實驗值相比略微偏高,而原始Schnerr-Sauer模型模擬得到的計算值比實驗值低,如圖5(a)所示。根據實驗描述[18],當空化數降低至σ=1.6時,初生空化形成,此時翼型頭部的低壓區流體因發生空化而形成了典型的游離型泡狀結構,并附著于水翼表面,從而導致升力系數出現了小幅上升。從圖中可以看到,修正后的Schnerr-Sauer模型能夠捕捉到升力系數這種特殊的行為特征,而原始的Schnerr-Sauer模型不具有這種能力。隨著空化數降低至σ=1.4時,游離型空泡逐漸過渡到片狀空化,此時水翼吸力面附著有一層呈片狀的空穴薄層,并最終導致升力急劇下降。

分析阻力系數可以看出,兩種空化模型預測得到的阻力系數均低于實驗值,在空化初生時,阻力系數才開始迅速增大,但是修正后的Schnerr-Sauer模型模擬得到的阻力系數與實驗值變化趨勢一致,即在空化數σ=1.6時才開始增大,而原始的Schnerr-Sauer模型在空化數為σ=1.8處時就開始增大;當空化數降低至σ=1.0左右時,云狀空化開始逐漸形成,這時阻力系數達到最大值,從圖中可以看出,修正后的Schnerr-Sauer模型模擬得到的阻力系數也達到最大值,而原始的Schnerr-Sauer模型計算得到的阻力系數保持繼續上升趨勢,表明原始的Schnerr-Sauer模型無法捕捉阻力系數在云狀空化條件下的這種特殊變化規律。綜上,無論是在非空化條件還是空化條件下,修正后的Schnerr-Sauer模型與原始的Schnerr-Sauer模型相比,更能準確地預測翼型所受到的升阻力。

圖5 時均化升阻力系數預測值與實驗值[18]對比Fig.5 Comparisons of time-averaged lift and drag coefficients between simulation and experiment[18]

4.2準周期性的云狀空化

4.2.1云狀空化形態演變

當空化數降低至σ=0.8時,附著于翼型表面上片狀空穴薄層的長度和厚度都出現了明顯增大,靠近尾緣處的空穴形態逐漸變得不穩定并有大量空泡脫落,云狀空化形成[16]。采用兩種空化模型計算得到的云狀空化及實驗結果如圖6所示,每個分圖中左欄為實驗結果[19],右欄上為采用修正后的Schnerr-Sauer模型的模擬結果,右欄下為采用原始的Schnerr-Sauer模型的模擬結果。需要注意的是,根據實驗測量,云狀空化演化周期為40 ms,采用修正的Schnerr-Sauer模型計算得到的云狀空化演化周期為35.9 ms,而原始的Schnerr-Sauer模型計算得到的云狀空化演化周期為30.2 ms,表明修正后的Schnerr-Sauer模型預測得到的云空化周期性演變過程與實驗結果基本一致。

云狀空化一個周期內演化主要分為以下幾個階段:1)在t=t0時刻,翼型前緣處片狀空穴幾乎消失,尾緣處空泡開始逐漸脫離翼型壁面,從圖中可以看到,原始的Schnerr-Sauer模型模擬得到的空穴形態與實驗結果相差較大,而修正后的Schnerr-Sauer模型的缺陷是無法模擬出此時翼型前緣處的少許片狀空穴,但能較好地模擬尾緣處逐漸脫離翼型表面的空泡團;2)在t=t0+38%T時刻,附著于翼型表面的片狀空穴發展至最大長度且尾緣處的空泡團消失,修正后的Schnerr-Sauer模型模擬的空穴長度與實驗結果吻合較好,而原始的Schnerr-Sauer模型無法準確模擬此時的片狀空穴形態;3)在t=t0+70%T時刻,翼型尾緣處的空泡團逐漸脫落并向下游移動,對比發現,修正后的Schnerr-Sauer模型能夠捕捉到尾緣處脫落的空泡,而原始的Schnerr-Sauer模型捕獲這一細節的能力欠佳;4)在t=t0+84%T時刻,翼型前緣的片狀空穴又開始逐漸回縮,尾緣處的大體積空泡團逐漸形成,開始重復t=t0時刻空泡形態特征。

4.2.2云空化水平速度分布特性

為了揭示云空化流場的運動特性,研究了流場內四個不同位置(x/c=0.4、x/c=0.6、x/c=0.8和x/c=1.2,其中x為水平方向到水翼前緣的距離,y為垂直方向到水翼表面的距離)上的水平速度分布,結果如圖7所示。整體上看,與原始的Schnerr-Sauer模型相比,采用修正后的Schnerr-Sauer模型計算得到的時均水平分速度與實驗值吻合程度更好。

結合圖6、7可以看出,采用修正后的Schnerr-Sauer模型更能準確反映水翼發生云空化時的流場情況。在x/c=0.4位置處,靠近翼型表面處流體具有較大的速度梯度,但遠離壁面處流體速度趨于一致,與主流速度較接近。在x/c= 0.8位置處,靠近翼型表面處流體的時均水平速度出現負值,表明該區域的速度方向與主流方向相反,說明該區域內出現了反向射流,隨著離壁面距離的增大,時均水平分速度也逐漸與主流速度達到一致;在x/c=1.2位置處,距離翼型表面不同位置處流場速度發生了急劇變化,出現了較大范圍的逆向水平速度,表明該區域內形成了較大區域的漩渦結構。在x/c=0.8和x/c=1.2位置處模擬得到的速度場與實驗值有較大差異,原因是所采用的SST k-ω模型本質上是標準雙方程湍流模型,該類湍流模型預測能力具有一定的局限性。由于該類湍流模型對速度場進行時均化處理因而無法對空化區域內的微小漩渦結構進行準確預測。

圖6 一個典型周期內實驗和模擬的云空化空泡演變對比Fig.6 Comparisons of cavity shape evolution between simulation and experiment

圖7 云空化條件下不同位置處時均水平分速度分布比較Fig.7 Comparisions of time-averaged u-velocity at different locations at cloud cavitation

5 結論

本文基于Schnerr-Sauer空化模型的思想和均相流假設建立了一種考慮不可凝結氣體及湍動能脈動影響的空化模型(修正后的Schnerr-Sauer模型)。在SST k-ω湍流模型中引入與混合密度相關的修正函數對湍流粘度系數進行了修正,分別采用原始Schnerr-Sauer模型及修正后的Schnerr-Sauer模型對二維Clark-Y水翼繞流空化流場進行數值模擬,并與文獻中的實驗結果作對比,得到的主要結論有:

1)與原始的Schnerr-Sauer模型相比,采用修正后的Schnerr-Sauer模型計算得到的升阻力系數與實驗值吻合程度更高,特別是捕捉空化初生時的升力系數及云狀空化形成時的阻力系數所表現出的特征。

2)修正后的Schnerr-Sauer模型與原始模型計算出的云空化周期分別為35.9 ms和30.2 ms。修正后的Schnerr-Sauer模型能較清晰地模擬云狀空泡的初生、發展、潰滅和脫落的全過程。

3)云狀空化階段,下游處(x/c=1.2)距離翼型表面不同位置處流場速度發生了急劇變化,出現了較大范圍的逆向水平速度,但本文所采用的湍流模型在預測該處流動速度時出現了較大偏差。

[1]季斌,洪方文,彭曉星.二維水翼局部空泡脫落特性數值分析[J].水動力學研究與進展,2008,23(4):412-418. JI Bin,HONG Fangwen,PENG Xiaoxing.The numerical analysis of shedding characteristics on partial cavitation over 2D hydrofoils[J].Chinese journal of hydrodynamics,2008,23(4):412-418.

[2]SINGHAL A K,ATHAVALE M M,LI Huiying,et al. Mathematical basis and validation of the full cavitation model[J].Journal of fluids engineering,2002,124(3):617-624.

[3]ZWART P J,GERBER A G,BELAMRI T.A two-phase flow model for predicting cavitation dynamics[C]//Proceedings of the 5th international conference on multiphase flow.Yokohama,Japan,2004.

[4]SAUER J,SCHNERR G H.Development of a new cavitation model based on bubble dynamics[J].Journal of applied mathematics and mechanics,2001,81(S3):561-562.

[5]MERKLE C L,FENG J,BUELOW P E O.Computational modeling of the dynamics of sheet cavitation[C]//Proceedings of the 3rd International Symposium on Cavitation. Grenoble,France,1998.

[6]KUNZ R F,BOGER A D,STINEBRING D R,et al.A preconditioned Navier-Stokes method for two-phase flows with application to cavitation prediction[J].Computers&fluids,2000,29(8):849-875.

[7]TAMURA Y,MATSUMOTO Y.Improvement of bubble model for cavitating flow simulations[J].Journal of hydrodynamics,ser.B,2009,21(1):41-46.

[8]DULAR M,BACHERT R,STOFFEL B,et al.Experimental evaluation of numerical simulation of cavitating flow around hydrofoil[J].European journal of mechanics-b/fluids,2005,24(4):522-538.

[9]LI D,GREKULA M,LINDELL P.A modified SST k-ω turbulence model to predict the steady and unsteady sheet cavitation on 2D and 3D hydrofoils[C]//Proceedings of the 7th International Symposium on Cavitation.Ann Arbor,Michigan,USA:University of Michigan,2009,16-20.

[10]FROBENIUS M,SCHILLING R,BACHERT R,et al. Three-dimensional unsteady cavitating effects on a single hydrofoil and in a radial pump measurement and numerical simulation[C]//Proceedings of the 5th International Symposium on Cavitation.Osaka,Japan,2003.

[11]LI Daqing,GREKULA M,LINDELL P.Towards numerical prediction of unsteady sheet cavitation on hydrofoils[J]. Journal of hydrodynamics,ser.B,2010,22(5):741-746.

[12]LOHRBERG H,STOFFEL B,FORTES-PATELLA R,et al.Numerical and experimental investigations on the cavitating flow in a cascade of hydrofoils[J].Experiments in fluids,2002,33(4):578-586.

[13]ZNIDARCIC A,METTIN R,DULAR M.Modeling cavitation in a rapidly changing pressure field-Application to a small ultrasonic horn[J].Ultrasonics sonochemistry,2015,22:482-492.

[14]ROOHI E,ZAHIRI A P,PASSANDIDEH-FARD M.Numerical simulation of cavitation around a two-dimensional hydrofoil using VOF method and LES turbulence model[J]. Applied mathematical modelling,2013,37(9):6469-6488.

[15]楊瓊方,王永生,張志宏.改進空化模型和修正湍流模型在螺旋槳空化模擬中的評估分析[J].機械工程學報,2012,48(9):178-185. YANG Qiongfang,WANG Yongsheng,ZHANG Zhihong. Assessment of the improved cavitation model and modified turbulence model for ship propeller cavitation simulation [J].Journal of mechanical engineering,2012,48(9): 178-185.

[16]洪鋒,袁建平,周幫倫.空泡半徑非線性變化的空化模型及應用[J].華中科技大學學報:自然科學版,2015,43(10):15-20. HONG Feng,YUAN Jianping,ZHOU Banglun.Cavitation model considering non-linear variation of bubble radius and its application[J].Journal of Huazhong university of science and technology:nature science edition,2015,43 (10):15-20.

[17]REBOUD J L,STUTZ B,COUTIER O.Two-phase flow structure of cavitation:experiment and modeling of unsteady effects[C]//Proceedings of the Third Symposium on Cavitation.Grenoble,France,1998.

[18]WANG Guoyu,SENOCAK I,SHYY W,et al.Dynamics of attached turbulent cavitating flows[J].Progress in aerospace sciences,2001,37(6):551-581.

[19]HU Changli,WANG Guoyu,CHEN Guanghao,et al.A modified PANS model for computations of unsteady turbulence cavitating flows[J].Science China physics,mechanics&astronomy,2014,57(10):1967-1976.

[20]OCHIAI N,IGA Y,NOHMI M,et al.Numerical Prediction of Cavitation Erosion Intensity in Cavitating Flows around a Clark Y 11.7%Hydrofoil[J].Journal of Fluid Science&Technology,2010,5(3):416-431.

本文引用格式:

洪鋒,袁建平,周幫倫,等.改進Schnerr-Sauer模型在水翼空化模擬中的評估分析[J].哈爾濱工程大學學報,2016,37(7):885-890.

HONG Feng,YUAN Jianping,ZHOU Banglun,et al.Assessment of improved Schnerr-Sauer model in cavitation simulation around a hydrofoil[J].Journal of Harbin Engineering University,2016,37(7):885-890.

Assessment of improved Schnerr-Sauer model in cavitation simulation around a hydrofoil

HONG Feng1,2,YUAN Jianping1,ZHOU Banglun1,SHI Yang1
(1.Research Center of Fluid Machinery Engineering and Technology,Jiangsu University,Zhenjiang 212013,China;2.Department of Mechanical Engineering,Eindhoven University of Technology,Eindhoven 5612AZ,Netherlands)

To enhance the capability of Schnerr-Sauer model in simulating cavitating flows,in this study,we developed an improved Schnerr-Sauer model based on the Rayleigh-Plesset equation and homogeneous flow assumption. The model considers the effects of turbulent fluctuation and noncondensable gas.The unsteady cavitating flow over a two-dimensional Clark-Y hydrofoil was numerically investigated combined with the SST k-ω turbulence model,in which the turbulence eddy viscosity is corrected.We obtain the time-averaged lift and drag coefficients,the cavity shape evolution at cloud cavitation,and its x-velocity profiles from these calculations.Compared with available experimental data in the literature,the time-averaged lift and drag coefficients predicted from the improved Schnerr-Sauer model agree better with the experimental results than with those of the original.In addition,the modified model can accurately predict the characteristics of the time-averaged lift and drag coefficients during different cavitation conditions.Moreover,the improved model can better simulate the inception,growth,shedding,and collapse of cloud cavitation than the original model.The overall results prove the reliability and accuracy of the improved Schnerr-Sauer model in cavitating flow simulations over a hydrofoil.

cavitating flow;hydrofoil;cavitation model;turbulence model;numerical simulation;drag coefficient;Schnerr-Sauer model

10.11990/jheu.201505010

S277.9;TH311

A

1006-7043(2016)07-885-06

2015-05-06.網絡出版日期:2016-05-04.

國家自然科學基金重點項目(51239005);江蘇省普通高校研究生科研創新計劃項目(KYLX_1040);江蘇高校優勢學科建設工程項目(PAPD).

洪鋒(1988-),男,博士研究生;袁建平(1970-),男,研究員,博士生導師.

洪鋒,E-mail:zjjaihf@163.com.

網絡出版地址:http://www.cnki.net/kcms/detail/23.1390.u.20160504.1435.002.html

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 亚洲网综合| 亚洲青涩在线| 亚洲人成网站日本片| 538精品在线观看| 亚洲第一成年网| 四虎精品免费久久| 久久精品国产999大香线焦| 国产全黄a一级毛片| 成人在线综合| 久草视频中文| 激情无码字幕综合| 日韩中文无码av超清| 鲁鲁鲁爽爽爽在线视频观看 | 国产日韩欧美精品区性色| 一区二区日韩国产精久久| 在线免费不卡视频| 亚洲女同一区二区| 波多野结衣视频一区二区 | 中文字幕av无码不卡免费 | 成人在线亚洲| 国产一区二区免费播放| 免费a级毛片视频| 91无码人妻精品一区| 99久久精彩视频| 亚洲免费成人网| 成人夜夜嗨| 国产呦视频免费视频在线观看| 午夜不卡视频| 高清无码不卡视频| 成人在线第一页| 在线中文字幕日韩| 国产精品女主播| 中文无码日韩精品| 亚洲视屏在线观看| 成年人视频一区二区| 人妻无码AⅤ中文字| 欧美爱爱网| 夜夜操狠狠操| 99精品伊人久久久大香线蕉| 免费啪啪网址| 99热最新在线| 亚洲色欲色欲www网| 国产无码精品在线播放| 欧美国产视频| 精品无码一区二区三区电影| 欧美另类图片视频无弹跳第一页| 天天综合网色| 国产丝袜一区二区三区视频免下载| 超碰免费91| 精品夜恋影院亚洲欧洲| 精品国产中文一级毛片在线看| 久久一色本道亚洲| 91色在线观看| 国产在线日本| 国产精品网拍在线| 国产成人综合亚洲欧美在| 久青草国产高清在线视频| 欧美一级一级做性视频| 国产精品免费p区| 国产一在线观看| 欧洲欧美人成免费全部视频 | 高清无码不卡视频| 免费人成又黄又爽的视频网站| 99久久精品国产麻豆婷婷| 国产黄色视频综合| 精品久久久久成人码免费动漫| 精品中文字幕一区在线| 91无码国产视频| 国产精品久久久久久久久久久久| 九九热精品在线视频| 国产三级视频网站| 免费看久久精品99| 日本不卡视频在线| 国产精品福利社| 国产成人无码久久久久毛片| 国产不卡在线看| 99ri国产在线| 国产农村妇女精品一二区| 在线国产你懂的| 亚洲视频免| 91人妻在线视频| 免费精品一区二区h|