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

繞水翼非定常空化流動的渦動力學分析

2015-04-26 02:45:40王國玉胡常莉陳廣豪北京理工大學機械與車輛工程學院北京100081
船舶力學 2015年8期
關鍵詞:區域

趙 宇,王國玉,黃 彪,胡常莉,陳廣豪,吳 欽(北京理工大學 機械與車輛工程學院,北京 100081)

繞水翼非定常空化流動的渦動力學分析

趙 宇,王國玉,黃 彪,胡常莉,陳廣豪,吳 欽
(北京理工大學 機械與車輛工程學院,北京 100081)

文章采用數值方法對繞水翼的非定常空化流動進行了流動計算,并采用渦動力學方法對非定??栈鲌鎏匦赃M行了分析。計算中,采用基于當地渦旋修正的湍流模型封閉控制方程,通過實驗對該計算結果進行驗證;采用渦量法分析了空化的發生和發展對渦量輸運過程的影響。結果表明:基于當地渦旋修正的湍流模型可以準確地模擬水翼周圍非定常空化流場結構;水翼空化流場具有強烈的非定常特性,空化的發展和流場中的渦結構有著緊密的聯系;反向射流的作用會導致附著空穴尾緣附近速度梯度的非定常變化,是導致渦量產生和消失的重要因素;空化區域內部水氣兩相的質量交換會導致流場內體積變化率的變化,使得空化區域內部流體的非正壓性顯著增強,上述兩者導致空化區域內部的渦量分布呈現強烈的非定常特性;空化現象引起的水氣兩相之間的相互轉換以及對湍流的作用都是影響渦結構空間分布的重要因素。

空化;湍流模型;非定常;渦動力學分析

0 引 言

作為高速水動力學研究領域的核心問題,空化是指在一定條件下,液相介質中發生的一種相變現象。在水流中,當其中的壓力降低到水的飽和蒸汽壓時,水流中將產生充滿氣體和蒸汽的空穴,這時,水流中發生了空化。這些空穴是一些單個的球形空泡,多個空泡形成的空泡群或外形不規則的空腔??栈前嘧儭⒎嵌ǔ?、多維湍流等多種復雜影響因素的流動現象,由空化非定常特性引起的一定尺度的空泡團的脈動、脫落是引起導彈出水過程水動力載荷、螺旋槳空泡噪音、振動和水力機械空蝕的主要原因[1-3]。渦運動是流體中普遍存在的一種運動形式,是湍流的一種基本結構[4]。渦運動因其固有的非定常特性和非線性機制,又因其對自然界的探索和工程應用都具有巨大的意義,一直是流體力學的研究前沿[5]??栈F象的非定常發展過程與渦運動息息相關[6],空泡的產生和發展會產生復雜的渦結構[7],而渦的形成可強化流場的湍流強度,使得湍流和空泡之間的交互作用變得復雜[8]。

近年來,國內外學者開始關注空化流場渦結構的演變過程與空化發展的相互作用,Gopalan和Katz[9]采用粒子測速系統和高速錄像系統對噴口結構的片狀空化流場進行研究,實驗結果表明渦的產生是空化流動的重要特性。Wang等[10]采用實驗方法研究了水翼周圍的附著空化湍流場的動態特性,結果表明云狀空化下大尺度渦結構和反向射流存在密切的關系。張博等[11]采用數值與實驗的方法研究了水翼周圍的云狀空化流動,發現空穴尾部大尺度的渦運動引起了近壁處的逆壓梯度,從而在壁面和附著空穴之間產生了液相的反向射流,反向射流逐漸向水翼前緣運動,從而導致了附著空穴的斷裂。Huang等[12]采用實驗和數值相結合的方法對水翼非定常云狀空化流動進行分析,發現渦結構的脫落是導致云狀空化流場強烈非定常性的主要原因。由此可見,非定??栈F象與流場中的旋結構緊密相關,了解分析空化發展過程與流場渦的產生、演化的相互作用是更加深入理解非定??栈鲃訖C理的重要途徑。渦動力分析方法是研究空化流場渦結構的重要方法,渦量輸運方程是分析空化和渦結構相互作用的有力工具。Nagendra等[13]采用大渦模擬的方法對文丘里管內部典型空化流場進行計算,并采用渦量輸運方程對文丘里管空化流動進行研究,結果表明空化現象的發生降低了渦動力方程中的拉伸項。黃彪[14]對水翼周圍云狀空化流場進行研究,和無空化的流場結構對比發現,空化現象的產生可以有效地增加渦量輸運方程中的斜壓矩產生項。空化流動的渦動力特性逐漸成為空化研究的熱點方法,但上述研究缺少從渦量輸運的角度分析典型空化流動的非定常特性。

本文采用數值方法對繞CLARK-Y水翼的非定??栈鲌鰷u結構進行研究,對水翼周圍非定常空化現象與渦結構的相互關系進行分析;并引入渦動力方程,分析空化現象對渦旋結構發展的影響。

1 數值方法

1.1 控制方程

計算采用單流體模型,適用于牛頓流體的笛卡爾坐標系下忽略體積力和熱傳導守恒形式的N-S方程如下:

式中:ρm是水汽混合相密度,ρl是液相密度,ρv是汽相密度,αv是汽相體積分數,αl是液相體積分數,u是速度,p是壓力,μm是混合相的粘性系數,μl和μv分別是液相和汽相的粘性系數,μT是湍流粘性系數。下標(i,j,k)表示笛卡爾坐標系下的坐標軸方向。和-分別表示凝結和蒸發速率。本文中,采用Kubota空化模型[15]:

式中:αnuc是空化核子的體積分數,RB是空泡直徑,pv是飽和蒸汽壓,p是當地壓力。Cdest和Cprod是凝結和蒸發系數。本文中,上述參數取值如下:αnuc=5×10-4,RB=1×10-6m,Cdest=50,Cprod=0.01。

1.2 湍流模型

本文采用修正標準k-ε湍流模型來封閉上述方程組,該模型在標準k-ε湍流模型基礎上引入當地渦旋運動系數,考慮了當地渦旋運動對湍流的影響。在標準k-ε湍流模型中,湍動能方程和耗散方程如下:

式中:Pk表示湍動能產生項,湍流粘性μt定義如下:

本文將當地渦旋運動效應修正引入到標準k-ε湍流模型中。在k方程中的湍動能產生項上乘以修正系數Pk→Pk·fr。修正因子的表達式如下[16]:

圖1 計算模型及邊界條件設置Fig.1 Geometry and Boundary Conditions

1.3 幾何及邊界條件設置

計算采用上述計算模型,流動介質參數與25℃的水和水蒸氣保持一致。水相密度和運動粘度分別設置為ρl=999.19 kg/m3和μl=ρlvl=1.139×10-3Pa·s;汽相設置為ρv=0.023 08 kg/m3和μv=9.862 6×10-6Pa·s。汽化壓力設置為pv=3 169 Pa。幾何模型和邊界條件如圖1所示,CLARK-Y水翼弦長0.07 m,最大相對厚度為11.7%,翼展0.07 m。水翼攻角設置為8°,入口處設置為法向速度入口,速度大小為10 m/s,從而保證雷諾數Re=U∞c/vl=7×105。出口處采用壓力邊界條件,設定環境壓力為43 169 Pa,從而保證空化數為(σ=(p-p∞)/(0.5ρv2))0.8。水翼表面和流場頂部和底部均設置為無滑移壁面邊界條件,前后邊界設置為對稱面。在計算中,時間步長Δt=1*10-4s,保證計算平均科朗特數CFL=U∞Δt/Δx=1。

1.4 模型驗證

為了驗證數值方法的可行性,將數值模擬結果和實驗結果進行對比。實驗在一閉式空化水洞[9]內進行。該空化水洞試驗段尺寸為0.70×0.07×0.19 m,流速范圍為2~15 m/s,可達到的最小空化數為0.3。在實驗中,采用高速攝像系統來觀察空泡形態的演變歷程,采用流體動力測量系統測量水翼受到的升力。

從高速攝像系統得到的系列空穴形態圖中可以發現水翼吸力面空泡表現為準周期性的發展、斷裂過程,呈現出顯著的非定常特性。根據空穴形態選取特征時刻,圖2給出了各特征時刻的實驗所得空穴形態圖和數值計算結果的對比情況。從圖中可以看出,數值方法預測的空化發展過程與實驗結果相吻合,能準確地捕捉空泡的初生、發展以及斷裂潰滅過程。

圖2 特征時刻空穴形態對比圖Fig.2 Instantaneous outline of cavities of experiment and simulation

定義升力系數為:

式中:U∞為距實驗段上游入口210 mm處參考斷面上的平均速度,A為水翼的有效升力面積,Fy為水翼所受到的阻力和升力。為了研究數值模型對預測翼型空化水動力特性的影響,圖3給出了計算得到的翼型升力隨時間的波動情況,并與實驗結果進行對比。從圖中可以看出,實驗和數值計算得到的水翼升力系數均呈現準周期性波動情況且波動的幅值基本一致,實驗測得的平均升力系數約為0.74,數值計算得到的平均升力系數為0.71,二者基本吻合。

圖3 翼型升力系數曲線對比Fig.3 Lift coefficient of experiment and simulation

2 渦量分析方法

為了分析水翼周圍的流場結構,定義特征渦量如下:

式中:δ≈0.1c=0.007 m為實驗測得的無空化條件下水翼尾緣的湍流邊界層厚度。

為了更清晰地反映流場內的渦旋結構,引入Q判據[17]:

Q是在伽利略變換下的速度梯度張量第二不變量,式中,Si j表示形變率張量,ωi j表示渦量,分別定義為:,其中i和j代表了不同的方向。對Q>0的區域定義為渦流所在,因為從上式可看出,Q可說是渦度和軸向的形變率的差值,當差值為正,表示旋轉的趨勢大于軸向的形變,也就是此區域受渦流所主宰。

引入渦量輸運方程:

式中第一項(I)表示由于流場的速度梯度引起的渦線的伸縮和彎曲,從而導致渦量的大小和方向都發生變化,在繞水翼的非定常空化流動中,空化現象引起流場內的速度梯度發生明顯的變化,從而導致渦量的絕對值發生變化,并且該速度梯度分布具有強烈的非定常性,從而導致速度梯度產生項具有強烈的非定常性。第二項(II)表示流體微團的體積變化引起的渦量大小變化。對于不可壓縮流體,該項絕對值為0,但由于空化現象的發生,使得流場內部空化區域的體積變化率發生變化,從而引起渦量的絕對值發生變化。第三項(III)表示由于不平行的壓力梯度和密度的梯度導致的斜壓矩對渦量的作用效果。對于正壓流體,其密度僅是壓力的函數,則密度和壓力項具有相同的變化梯度,該項為0,但在空化流場中,壓力和密度的梯度并不總是平行的,特別是在空穴界面及閉合區域內,從而會引起渦量的變化。第四項(Ⅳ)表示渦量的粘性擴散效應。定義特征渦量輸運項如下:

3 結果分析

圖4分別給出了特征時刻下的沿翼展方向的渦量(根據特征渦量進行量綱一化處理)等值面分布情況,以及對應時刻的Q分布情況。綜合圖2和圖4可以看出在t1=0.25T0時刻,水翼前緣位置開始出現少量片狀附著空穴,水翼吸力面尾緣附近仍能觀察到上個周期脫落的空泡團。水翼前緣的附著空穴逐漸發展,脫落空泡團逐漸脫離水翼表面,向下游運動。從渦量分布圖中可以看出,分別從水翼前緣和尾緣開始,產生方向相反的渦量聚集帶狀區域并逐漸向下游延伸。前緣渦(-)和尾緣渦(+)相互作用,在水翼吸力面附近匯合,并逐漸脫落,匯合區域對應于水翼吸力面的空化區域。

到t2=0.5T0時刻,附著空穴長度增加到最大,在該時刻,空穴邊界開始出現不穩定的波動,特別是在附著空穴的末端附近;此時前緣渦帶逐漸向下游發展,并較為穩定地附著在水翼吸力面上,附著空穴位于前緣渦帶中,尾緣渦控制范圍較小。從Q分布中可以看出,附著空穴內部Q絕對值較小,說明在附著空穴內部即前緣渦帶核心位置的流動較為穩定,旋轉效應和變形效應相對平衡,流動沒有較強的渦旋特性。而在附著空穴的末端附近,Q絕對值較大,且正負區域間隔分布,這說明在該區域內部,旋轉效應和變形效應處于非平衡狀態,兩者相互轉化,相互制衡;該區域對應于前緣渦帶和尾緣渦的交匯處,方向相反的渦結構相互作用,導致該區域內部出現強烈的不穩定性。

圖4 特征時刻渦量分布和Q分布圖Fig.4 Instantaneous contours of dimensionless vorticity and Q

隨著空化現象的發展,水翼尾部附近的指向前緣渦帶核心區域的反向壓力梯度逐漸發展增大到足以克服由于流動產生的正方向動量,在近壁面區域產生反向射流,在反向射流前端形成局部高壓,并推動近壁面的液相流體向水翼前緣運動;從圖中可以看出在前緣渦帶(-)和水翼壁面之間存在具有一定厚度的正向渦帶(+),這是由于前緣渦和尾緣渦的相互作用導致尾緣渦結構的斷裂和脫落,一部分帶有正向渦量的流體脫離水翼尾緣,另一部分帶有正向渦量的流體在反向射流的作用下沿壁面向水翼前緣運動。到t3=0.75T0時刻,反向射流前端逐漸向水翼前緣移動,到達附著空穴前端附近,附著空穴在高速水流的作用下發生斷裂,反向射流區域內流體帶有正向的渦量(+)。在斷裂位置處和水翼尾緣附近,存在方向相反的渦結構相互作用,旋轉效應和變形效應相互制約,相互轉換。

隨著空化的進一步發展,由于慣性的作用,反向射流繼續向水翼前緣運動,延緩了水翼前緣附著空穴的迅速增長,另一方面,附著空穴斷裂的另一部分在主流的作用下向下游運動,從而形成了大尺度的脫落空泡團。到t4=0.9T0時刻,附著空穴逐漸發展,開始逐漸覆蓋水翼吸力面,大尺度的脫落空泡團逐漸脫離水翼吸力面。從渦量云圖中可以看出,附著空穴處于前緣渦帶的位置處,脫落的空泡內部渦結構較為復雜,具有相反方向渦量的流體存在強烈的相互作用,這是由于具有正向渦量的反向射流沖擊具有相反方向的前緣渦帶,使之發生斷裂,從而形成脫落的渦旋結構,脫落的渦旋結構內部同時保留著反向射流的正向渦量以及前緣渦帶的反向渦量,所以導致脫落空泡內部存在較為復雜的渦旋結構。

圖5中給出了不同時刻下沿翼展方向的渦量各輸運項等值面云圖分布情況,圖中各輸運項根據特征渦量輸運項進行量綱一化處理。從圖中可以看出,I項發展規律具有強烈的非定常特性,在t1= 0.25T0時刻,絕對值較大的位置對應于在附著空穴尾緣和水翼尾緣附近,在t2=0.5T0時刻,附著空穴逐漸發展到覆蓋整個吸力面,絕對值較大區域出現在水翼尾緣,從t3=0.75T0時刻到t4=0.9T0時刻,絕對值較大位置對應于反向射流的前端,從水翼尾緣逐漸向前緣發展。這是由于在附著空穴尾緣,在反向射流的作用下出現速度旋渦,該區域內部的速度梯度分布較為復雜,因而導致該區域內部的I項的絕對值較大。

圖5 特征時刻渦量輸運項分布圖Fig.5 Instantaneous contours of dimensionless vorticity transport term

絕對值較大的II項分布逐漸覆蓋整個吸力面,并且具有一定的非定常特性:在t1=0.25T0時刻,在脫落空泡團的核心位置的II項絕對值較大,從t2=0.5T0時刻到t4=0.9T0時刻,附著空穴逐漸發展,附著空穴核心位置對應于II項絕對值較大的區域。這是由于II項表征的是流體微團的可壓縮性,在附著空穴和脫落空泡核心位置,水蒸氣含量較大,水氣之間存在強烈的質量交換,從而導致強烈的可壓縮效果。

絕對值較大的III項區域分布規律具有強烈的非定常特性。在t1=0.25T0時刻,絕對值較大的區域主要分布在水翼尾緣附近,水翼前緣出現絕對值較大區域但體積較小,在t2=0.5T0時刻,絕對值較大區域相互交錯,覆蓋整個水翼吸力面,從t2=0.5T0時刻到t4=0.9T0時刻,水翼吸力面的絕對值較大的III項區域出現斷裂,并且斷裂位置逐漸向水翼前緣運動。綜上所述,絕對值較大的III項區域和空化區域相一致,這是由于空化現象的發生,使得空化區域內部的流場結構處于非平衡狀態,水蒸氣含量不斷變化,并且速度旋渦,流體的非正壓性得到顯著加劇,從而導致流場內部的密度梯度與壓力梯度方向不一致。

絕對值較大的Ⅳ項分布逐漸覆蓋整個吸力面,并且隨著空化現象的發展呈現出一定的非定常特性。這是由于空化現象的發生發展導致水翼吸力面處于氣液混合狀態,水氣之間存在的相互轉換導致該區域內部流體的粘度發生顯著的非定常變化,而同時,空化現象的發生也導致水翼吸力面的湍流強度產生強烈的非定常變化,上述兩種因素的作用,導致水翼吸力面附近的Ⅳ項取值較大。

4 結 論

本文采用數值方法對繞水翼的非定??栈瘻u動力流場特性進行研究,分析了水翼周圍非定??栈F象與渦結構的相互關系,并引入渦動力方程,分析空化現象對渦結構發展的影響,得到主要結論如下:

(1)水翼空化的發展和流場中的渦結構有著緊密的聯系,基于當地渦旋修正的湍流模型可以準確地模擬水翼周圍非定??栈鲌鼋Y構。

(2)反向射流的作用會導致附著空穴尾緣附近速度梯度的非定常變化,是導致渦量產生和消失的重要因素。

(3)空化區域內部水氣兩相的質量交換會導致流場內體積變化率的變化,空化區域內部流體的非正壓性顯著增強,上述兩者導致空化區域內部的渦量分布呈現強烈的非定常特性。

(4)空化現象引起的水氣兩相之間的相互轉換以及對湍流的作用都是影響渦結構空間分布的重要因素。

[1]劉筠喬,魯傳敬,李 杰,等.導彈垂直發射出筒過程中通氣空泡流研究[J].水動力學研究與進展A輯,2007,22(5): 549-554. Liu Yunqiao,Lu Chuanjing,Li Jie,Cao Jiayi.An investigation of ventilated cavitating flow in vertical launching of a missile[J].Journal of Hydrodynamics(Ser.A),2007,22(5):549-554.

[2]Cervone A,Bramanti C,Rapposelli E.Experimental characterization of cavitation instabilities in a two-bladed axial indeucer[J].Journal of Propulsion and Power,2006,22(6):1389-1395.

[3]顧 巍,何友聲,胡天群.軸對稱體空泡流的噪聲特性與空泡界面瞬態特征[J].上海交通大學學報,2000,34(8):1026-1030. Gu Wei,He Yousheng,Hu Tianqun.Noise feature of cavitation flows on axisymmetric bodies and transient characteristics of cavity interface[J].Journal of Shanghai Jiaotong University,2000,34(8):1026-1030.

[4]童秉綱,尹協遠,朱克勤,等.渦運動理論[M].合肥:中國科學技術大學出版社,2009. Tong Binggang,Yin Xieyuan,Zhu Keqin.Vortex dynamics theroy[M].Anhui:Press of University of Science and Technology of China,2009.

[5]Wu J Z,Ma H Y and Zhou M D.Vorticity and vortex dynamics[M].Berlin:Springer-Verlag,2006.

[6]Li C Y,Ceccio S L.Interaction of single travelling bubbles with the boundary layer and attached cavitation[J].Journal of Fluid Mechanics,1996,322:329-353.

[7]Belahadji B,Franc P,Michel M.Cavitation in the rotational structures of a turbulent wake[J].J Fluid Mech.,1995,287: 383-403.

[8]Arakeri V H,Acosta A J.Viscous effects in the inception of cavitation on axisymmetric bodies[J].Journal of Fluids Engineering,1973,95(4):519-527.

[9]Gopalan S,Kate J.Flow structure and modeling issues in the closure region of attached cavitation[J].Physics of Fluids, 2000,12(4):3414-3431.

[10]Wang G,Senocak I,Shyy W.Dynamics of attached turbulent cavitating Flows[J].Progress in Aerospace Sciences,2001, 37:551-581.

[11]張 博,王國玉,黃 彪,等.云狀空化非定常脫落機理的數值與實驗研究[J].力學學報,2009,41(5):651-659. Zhang Bo,Wang Guoyu,Huang Biao,et al.Numerical and experimental studies on unsteady shedding mechanisms of cloud cavitation[J].Chinese Journal of Theoretical and Applied Mechanics,2009,41(5):651-659.

[12]Huang B,Wang G.Experimental and numerical investigation of unsteady cavitating flows through a 2D hydrofoil[J].Science China Technological Science,2011,54(7):1801-1812.

[13]Nagendra Dittakavi,Aditya Chunekar,Steven Frankel.Large eddy simulation of turbulent-cavitation interactions in a venturi nozzle[J].Journal of Fluids Engineering,2010,132(12):

[14]黃 彪.非定常空化流動機理及數值計算模型研究[D].北京:北京理工大學,2012. Huang Biao.Physical and numerical investigation of unsteady cavitating flows[D].Beijing:Beijing Institute of Technology, 2012.

[15]Kubota A,Kato H,Yamaguchi H.A new modeling of cavitating flows:A numerical study of unsteady cavitation on a hydrofoil section[J].Journal of Fluid Mechanics,1992,240(1):59-96.

[16]Smirnov P E,Menter F R.Sensitization of the SST turbulence model to rotation and curvature by applying the Spalart-Shur correction term[C]//ASME Paper GT 2008-50480.Berlin,Germany,2008.

[17]Haller G.An objective definition of a vortex[J].J Fluid Mech,2005,525:1-26.

Vortex dynamic analysis of unsteady cavitating flows around a hydrofoil

ZHAO Yu,WANG Guo-yu,HUANG Biao,HU Chang-li,CHEN Guang-hao,WU Qin
(School of Mechanical and Vehicular Engineering,Beijing Institute of Technology,Beijing 100081,China)

Unsteady cavitating flows around a hydrofoil are calculated.Vortex dynamics are used to analyze the flow structure.Local swirling correction turbulence model is adopted to close the governing equations.Validation of the LSC turbulence model can be reached comparing with experimental results.Vorticity transport equations are conducted to analyze the influence of cavitation on vorticity transport process. The results show that:The LSC turbulence model can reflect the influence of the local swirling on cavitating flows well.It can get a good compromise with experimental results.Unsteady cavitation process contains developing of complex vortex structures.Re-entrant jet near trailing of attached cavity leads to distinct changes of velocity gradient,which has great influence on production and dissipation of vorticity.Intensive mass transfer between liquid and vapor phase may induce the volumetric expansion or contraction and baroclinic torque.Two factors mentioned above keep changing and result in unsteady distributions of vorticity.Cavitation may lead to unsteady mass and momentum transfers between liquid and vapor phase,which will change turbulence and vorticity distributions distinctly.

cavitation;turbulence model;unsteady;vortex dynamics analysis

O35 TV131.2

A

10.3969/j.issn.1007-7294.2015.08.002

1007-7294(2015)08-0896-09

2015-04-23

國家自然科學基金重點項目(51239005,51479002)

趙 宇(1989-),男,博士研究生;

王國玉(1961-),男 教授,博士生導師,E-mail:wangguoyu@bit.edu.com。

猜你喜歡
區域
分割區域
探尋區域創新的密碼
科學(2020年5期)2020-11-26 08:19:22
基于BM3D的復雜紋理區域圖像去噪
軟件(2020年3期)2020-04-20 01:45:18
小區域、大發展
商周刊(2018年15期)2018-07-27 01:41:20
論“戎”的活動區域
敦煌學輯刊(2018年1期)2018-07-09 05:46:42
區域發展篇
區域經濟
關于四色猜想
分區域
公司治理與技術創新:分區域比較
主站蜘蛛池模板: www.youjizz.com久久| 亚洲另类国产欧美一区二区| 亚洲国产中文在线二区三区免| 欧美黄网站免费观看| 亚洲成年人网| 伊伊人成亚洲综合人网7777| 99热这里只有精品在线播放| 一级一级特黄女人精品毛片| 伊人精品成人久久综合| 亚洲综合18p| 极品国产一区二区三区| 欧美激情视频二区| 伊人国产无码高清视频| 亚洲黄色片免费看| 亚洲国产精品一区二区高清无码久久| 亚洲综合欧美在线一区在线播放| 91在线免费公开视频| 亚洲欧美日韩成人在线| 中国一级毛片免费观看| 久久亚洲国产视频| 亚洲高清国产拍精品26u| 99视频精品全国免费品| 58av国产精品| 五月激情婷婷综合| 狼友视频一区二区三区| 久久久久国产精品嫩草影院| 大学生久久香蕉国产线观看| 久久大香伊蕉在人线观看热2| 18禁色诱爆乳网站| 青青青视频蜜桃一区二区| 重口调教一区二区视频| 精品欧美视频| 久久一日本道色综合久久| 亚洲网综合| 日韩高清一区 | 99手机在线视频| 女人18一级毛片免费观看| 97成人在线视频| 国内精品免费| 国产久操视频| 在线免费观看AV| 日本三级精品| 亚洲精品制服丝袜二区| 国产成人精品高清不卡在线| 欧美h在线观看| 在线国产你懂的| 国产偷倩视频| 日本亚洲成高清一区二区三区| 亚洲视频色图| 免费人成在线观看视频色| av尤物免费在线观看| 日本午夜影院| 国产成人91精品| 在线a视频免费观看| 乱人伦99久久| 五月婷婷中文字幕| 中文字幕亚洲无线码一区女同| 另类重口100页在线播放| 亚洲无线一二三四区男男| 国产成人综合久久精品尤物| 在线播放国产一区| 九九这里只有精品视频| 试看120秒男女啪啪免费| 国产精品露脸视频| 午夜福利网址| 99伊人精品| 99热这里只有精品5| 亚洲人成高清| 四虎永久免费地址| 少妇精品久久久一区二区三区| 19国产精品麻豆免费观看| 国产呦精品一区二区三区网站| 亚洲三级影院| 久久久国产精品免费视频| 久久精品无码国产一区二区三区 | 女同久久精品国产99国| 最新日韩AV网址在线观看| 亚洲天堂成人在线观看| 国模在线视频一区二区三区| 欧美中文字幕在线视频| 色天天综合久久久久综合片| 国产黄视频网站|