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

基于IVCBC渦方法高雷諾數下串聯雙圓柱的數值計算模型及其尾流特征研究

2017-11-02 06:30:00龐建華王大明
船舶力學 2017年10期

龐建華,宗 智,周 力,王大明

(1.大連理工大學 工業裝備與結構分析國家重點實驗室,遼寧 大連116024;

2.四川探險家航空科技公司,成都 610021)

基于IVCBC渦方法高雷諾數下串聯雙圓柱的數值計算模型及其尾流特征研究

龐建華1,宗 智1,周 力1,王大明2

(1.大連理工大學 工業裝備與結構分析國家重點實驗室,遼寧 大連116024;

2.四川探險家航空科技公司,成都 610021)

針對串聯雙圓柱的結構特點,利用IVCBC渦方法適合高雷諾數下數值計算的特點,建立了高雷諾數下串聯雙圓柱繞流的數值計算模型;采用經典算例驗證了IVCBC渦方法的收斂性;探索了雷諾數為2.5×104、間隙比分別為1.1,1.25,2,2.5,3,3.25和4的串聯雙圓柱繞流的尾流特征;清晰地展示了尾流中較小的漩渦的形成、分裂和融合,詳細地闡述了串聯雙圓柱流體特征發生突變的原因;深刻地揭示了串聯雙圓在高雷諾數下的繞流機理。研究表明:尾流模式與經典的實驗尾流模型吻合較好,兩圓柱中間的渦對是串聯雙圓尾流發生突變的主要原因;間隙之間時,間隙上部小漩渦形成、間隙中間流體的振動與下游圓柱表面上渦的脫落是同步的;該流體模式能清晰地展示尾流中較小漩渦,說明與有網格方法比較,該計算模型具有較大的優越性,為進一步研究高雷諾數下串聯雙圓柱流體力的特征提供了重要的研究工具。

串聯雙圓柱;IVCBC渦方法;尾流特征;高雷諾數;繞流機理

0 引 言

串聯多圓柱繞流特征與單圓柱繞流有本質的區別,探索其特點對實際的工程設計和制造有重要的指導意義。在低雷諾數下圓柱繞流已經被廣泛地研究,在實際的海洋環境中,裝備大多處于高雷諾數的工況下。然而基于高雷諾數的雙圓柱繞流研究甚少。

串聯雙圓柱繞流特性也得到廣泛的研究。Mittal等人[1]研究了雷諾數為100-10 000串聯雙圓柱的非定常流動的特點,他們發現串聯雙圓柱定性的流動特征強烈地依賴于圓柱的間距和雷諾數的大小。Li等人[2]采用伽遼金速度壓力有限元法和粗網格對雷諾數100下的兩個串聯圓柱的繞流進行了研究。Johnson等人[3]等采用基于渦量流函數的彼得洛夫—伽遼金有限元法對雙圓柱低雷諾數進行了研究。Farrant等人[4]采用分子邊界元計算對雷諾數為200、間距比為4的串聯雙圓柱渦的脫落特點進行了數值研究。Hall等人[5]發現有聲駐波對在低雷諾數下和間隙比為1.2~1.75的串聯雙圓柱的尾流有明顯的影響。Meneghini等人[6]用分步的方法研究雷諾數為100~200的渦流脫落特征,以及渦流和串聯圓柱之間的尾流相互干擾的特點。Sharman等人[7]采用同位非結構化的計算流體動力學代碼研究雷諾茲數100的兩個串聯圓柱,他們發現兩圓柱的臨界間距比3.75~4。當間距大于臨界間距時,他們發現在尾流圓柱只有一個觸接點和一個分離點。特別是,考慮到一個或兩個圓柱自由擺動時,兩圓柱體旋渦脫落和下游圓柱的尾跡結構模式會有較大改變。Ding等人[8]采用了基于無網格最小二乘有限差分法,對雷諾數為100~200、間距比為L/D=2.2的串聯雙圓柱的尾跡進行了研究,并對低雷諾數范圍的流場進行了定量的分析。Deng等人[9]對雷諾數為220的串聯雙圓柱的繞流進行了三維數值模擬研究。他們改進了虛擬邊界法來滿足無滑移條件,如同二維的情況,他們發現臨界間距范圍內的三維不穩定發生的臨界范圍是在3.5≤L/D≤4,這意味著間距比小于3.5,尾流保持二維狀態,間距比L/D≥4三維現象將出現。Singha等人[10]采用一階隱式有限體積法研究了雷諾數為40~150的串聯雙圓柱層流的繞流特點。他們觀察了間距比L/D=0.2~4、雷諾數為100的尾流特點:當L/D≤2時流體保持穩定,當L/D≥3層流變得不穩定。當L/D≤2沒有明顯的旋渦從上游圓柱脫落。

雖然雙圓柱的繞流在低雷諾數下已得到廣泛的研究,但在高雷諾數下的研究非常少。盡管繞流的機理在實驗中也能較好地展示,但許多特征未能很清晰地展示,比如較小的渦流在網格方法中很難被展示。據目前的文獻資料,仍然沒有一個數值計算模型能較好地研究高雷諾數下串聯雙圓柱的繞流特征。為此,本文采用IVCBC渦方法[11-13]在高雷諾數下不需要湍流模式的優點,結合串聯雙圓柱的邊界條件,建立了串聯雙圓柱繞流的數值計算模型。對雷諾數為2.5×104,間隙比在1.1<L/D<6范圍內7種間距比的串聯雙圓柱繞流進行了研究。清晰地展示了尾流中較小的漩渦的形成、分裂和融合,詳細地闡述了串聯雙圓柱流體特征發生突變的原因;深刻地揭示了串聯雙圓柱在高雷諾數下的繞流機理。

1 數值模型

根據平面點渦的性質,流函數的形式為:

圖1 漩渦產生的位置Fig.1 The model of vortex generation

根據流函數的可疊加性,在第K個控制點的流函數為:

對于光滑的圓柱表面,其表面具有流線性,流函數的值在二維的圓柱表面處處相等。因此根據兩個相鄰控制點上的流函數,可建立如下方程:

由于2M控制點位于兩圓柱表面上,所以有2M個方程,并且包含2M個未知量,而且第2M+1個方程代表渦量守恒。根據最小二乘法可獲得一個線性方程組:

式中:Γi為渦元的強度,通過高斯消元法獲得方程組的解。

根據所有的渦元的強度和位置,獲得每個圓柱受到的升力和阻力,其計算公式如下:

渦元的對流、擴散、融合和消去的計算公式與IVCBC法[11-12]中描述的一致。

2 IVCBC渦方法的收斂性

本文采用無量綱時間步長,Δt=aD/U,其中a是時間步長的一個修正值,D是立管半徑,U是均勻來流的速度。為了研究時間步長和渦數對IVCBC渦方法的收斂性影響,本文選擇了兩個時間步長,修正值分別設定為a=0.1,0.05,圓柱表面產生新生渦的個數分別為32,64,128,這樣形成六種工況。圖2展示六種組合在雷諾數為Re=1×105的壓力分布與Blevins等人[15]的實驗結果。從這個圖中能觀察到組合為N=128、Δt=0.05的結果最接近實驗結果。這同樣表明,渦數越多、時間步長越小,計算結果越收斂。

圖2 平均壓力系數(Re=1×105時間步長為Δt=0.1和Δt=0.05)Fig.2 Mean surface pressure coefficients at Re=1×105with Δt=0.1 and Δt=0.05 by the present method

3 尾流模式

為了獲得高雷諾數下的串聯雙圓柱的特征,根據IVCBC渦方法[12-13]的收斂性的驗證,本文選擇如下計算參數為:雷諾數為2.5×104,新生渦元的數量為128,時間步長為0.05。其中兩圓柱的中心間距與圓柱的直徑之比 L/D 分別為 1.1,1.25,2,2.5,3,3.25和4。

圖 3~9 給出了當 L/D=1.1,1.25,2,2.5,3,3.25和 4時,渦流分布圖。 為了便于比較,Igarashi等人[16]經典的串聯雙圓柱繞流模型也標注在圖中。圖3展示L/D=1.1的渦元分布圖。從圖中,我們觀察到從上游圓柱體的剪切層分離的渦流包裹下游圓柱,但渦流沒有附觸到下游圓柱表面,而是在下游圓柱的尾流中相互作用形成卡門渦街。與單圓柱繞流的尾流特征相似,只是渦街的形狀變大。這是因為渦流從上游圓柱的剪切層分離后,連接到下游圓柱的剪切層,等價于延長了上游圓柱剪切層的長度。從圖中也可以觀察到,在下游圓柱體后側形成了較小的漩渦,如圖所標示。但在兩圓柱的間隙中沒有形成較小的漩渦。這是因為上游圓柱的剪切層脫落后形成的渦流沒有進入兩圓柱之間的間隙,盡管在間隙的內部圓柱表面也會產生渦元,但其被上游圓柱和下游圓柱之間的間隙包裹,渦元流動性較弱。

圖3 L/D=1.1時瞬時渦元分布輪廓Fig.3 Instantaneous vorticity contours in the case of L/D=1.1

圖4 展示的是L/D=1.25的串聯雙圓柱渦流分布情況。與L/D=1.1的渦流分布略有不同,從上游圓柱體上側的剪切層分離形成的渦流有部分附觸在下游圓柱的上側。這是因為當兩圓柱之間的間距變大,上游圓柱上側剪切層分離的渦流不能完全包裹下游圓柱。同時由于下游圓柱周期性地漩渦脫落,導致上圓柱剪切層脫落的渦流也周期性地上下擺動,導致有小一部分渦流在下游圓柱上側,周期性地阻擋,回流到雙圓的間隙中。在兩圓柱之間的間隙上端形成較小的漩渦如圖4(iii)所示。這些進入間隙中的渦流隨著小漩渦的形成和消失在圓柱間隙中間上下擺動。小漩渦的形成也和間隔中渦流的振蕩與下游圓柱尾流漩渦的脫落是同步的,這是因為這三種狀況都與下游圓柱旋渦的脫落相關。本文的渦元分布圖與經典的Zdravkovich等人[17]的實驗中獲得的流體分布圖一致,如圖4(i)所示。

圖4 L/D=1.25時的瞬時渦元分布輪廓Fig.4 Instantaneous vorticity contours in the case of L/D=1.25

圖5 (ii),(iii)展示了L/D=2時的渦流分布圖。從圖中我們可以觀察到從上游圓柱剪切層的上側脫落的渦流受到下游圓柱的阻擋,大部分渦流進入兩圓柱之間的間隙,形成較大的兩個不對稱的環流。這是因為下游圓柱漩渦脫落的周期性影響了兩圓柱之間的渦流長度和擺動劇烈程度,導致渦流被阻擋,回到間隙中的速度和流量都不一樣,從而導致間隙中的渦流不對稱。渦流的特點與Zdravkovich等人[17-18]及 Igarashi等人[16]的實驗結果相同,如圖 5(i)所示。

圖5 L/D=2時的瞬時渦元分布輪廓Fig.5 Instantaneous vorticity contours in the case of L/D=2

L/D=2.5時的兩圓柱的尾流特點與L/D=2時的尾流特點幾乎相同,如圖6所示。只是從上游圓柱剪切層分離的渦流的運動寬度比L/D=2時窄小一些。這是因為上游圓柱和下游圓柱之間的間距增大,使得回流到間隙中的渦流有更大的回流空間,由此可看出,兩圓柱之間形成的環流,其速度隨間距比的增大而減小。

圖6 L/D=2.5時的瞬時渦元分布輪廓Fig.6 Instantaneous vorticity contours in the case of L/D=2.5

圖7 展示了L/D=3時的三種渦流模型,其中圖(a)展示從上游圓柱剪切層分離的渦流被下游圓柱阻擋,在兩圓柱間隙中形成兩個對稱的環流,而且在下游圓柱的尾流中形成多個較小的漩渦,這種漩渦不再是卡門渦街。這是因為上游圓柱剪切層分離的渦流被對稱性地阻擋,下游圓柱上下表面周期性分離的渦流強度減弱,導致了尾流中沒有較大的渦街形成,而是出現多個較小的漩渦。而在圖(b)中,在下游圓柱的尾流中形成兩個較大的漩渦,并且一個距離下游圓柱較近,另一個距離下游圓柱較遠。這是因為在上游圓柱的上下剪切層的渦流周期性地被后圓柱阻擋,在兩個圓柱之間的間隙內形成不對稱的環流,導致下游圓柱上下表面的流量差變大,增強了從下游圓柱剪切層中脫落渦的強度,因而尾流中產生較大的漩渦。這種輪廓特點與Zdravkovich等人[17-18]的相似,如圖(b)(iii)所示。然而圖(c)下游圓柱體的尾流與經典的流動模型的相似,但是上游圓柱剪切層的渦流是不對稱的接近下游圓柱體的上下側,一側的渦流接近下游圓柱的上側,而另外一側的渦流被全部阻擋。這表明上游圓柱的尾流即將被全部阻擋,在兩圓柱間隙中形成渦對。這個狀態是串聯雙圓柱繞流發生突變的前兆。這種流動特點與 Zdravkovich 等人[17-18]及 Igarashi等人[16]的結論相同,如圖(c)(iii)所示。

圖7 L/D=3時的瞬時渦元分布輪廓Fig.7 Instantaneous vorticity contours in the case of L/D=3

前面所給出的4種間隙比的尾流特征,都有一個共同的特點,即:在兩個圓柱間隙之間存在環流,但沒有形成渦對。這一特征是區別串聯雙圓柱繞流的流體力產生突變的主要依據,也表明在未發生突變前,所有的間隙比的流體力具有相似的特點。

圖8展示L/D=3.25時渦元分布,從圖中可以看出其渦流的流動明顯不同于前面的5種渦流形式,上游圓柱體的剪切層的渦流不再附觸在下游圓柱體上,而是從上游圓柱上下側兩側交替地脫落,在兩圓柱之間形成卡門渦對。隨著時間的推移,這些漩渦撞擊下游圓柱體的前側,被分裂和破碎,與下游圓柱脫落的渦,在下游圓柱后側融合,形成穩定的卡門渦街。在下游圓柱后側清晰地觀察到成對的較小的漩渦。這些較小的漩渦最后被上游圓柱的渦流融合形成新的較大的漩渦。再一次表明:雙圓柱之間形成的渦對,是串聯雙圓柱繞流發生突變重要的標志。也是串聯雙圓柱繞流的流體力發生突變的重要標志。

圖8 L/D=3.25時的瞬時渦元分布輪廓Fig.8 Instantaneous vorticity contours in the case of L/D=3.25

圖9給出了L/D=4的渦流分布。為了便于分析漩渦的形成與破裂,在圖9中分別用紅色標出兩圓柱上側的渦流,用藍色標出兩圓柱下側的渦流。在上游和下游圓柱之間觀察到一對逆向旋轉的渦對,從渦對的顏色可以看出紅色的漩渦來自上游圓柱上側的剪切層,藍色的漩渦來自上游圓柱的下側。隨后,紅色的漩渦撞擊下游圓柱體的前側,紅色漩渦被切分成上下兩個部分,其中上面的那邊部分被下游圓柱上側剪切側的渦流脫落形成的漩渦融合,在尾流中形成較大的漩渦。而下面部分,與下游圓柱的下側剪切層的渦流相融合,在尾流形成新的漩渦。這個漩渦與之前形成的漩渦形成旋轉方向相反的渦對,這兩渦對具有相同的特性。這種融合的過程與Ljungkrona等人[19]的實驗一致,如圖9(i)所示。

圖9 L/D=4時的瞬時渦元分布輪廓Fig.9 Instantaneous vorticity contours in the case of L/D=4

4 結 語

針對串聯雙圓柱的結構特點,利用IVCBC[11-12]渦方法適合高雷諾數下數值計算的特點,建立了高雷諾數下串聯雙圓柱繞流的數值計算模型;采用經典算例驗證了IVCBC[11-12]渦方法的收斂性;探索了雷諾數為2.5×104、間隙比分別為1.1,1.25,2,2.5,3和4的串聯雙圓柱繞流的尾流特征。研究發現:

(1)當L/D<3.25時,我們可以觀察到漩渦脫落現象只發生在下游圓柱體處。當L/D≤1.25時,從上游圓柱體的剪切層分離的渦流連接下游圓柱的剪切層,并包裹著下游圓柱體,渦流沒有附著在下游圓柱上下表面。當1.25<L/D≤2時,間隙中的渦流作周期性地上下振蕩,頻率與下游圓柱體的泄渦頻率相同。這表明,間隙之間的小漩渦的形成,兩圓柱間隙中流體的振動與下游圓柱表面渦的脫落同步。

(2)當2≤L/D<3時,從上游圓柱體剪切層分離的渦流,周期性地附觸在下游圓柱體表面的上下兩側,在兩圓柱之間的間隙中形成環流,環流的速度隨間距比的增大而減小。

(3)當L/D=3時有三種流動模型出現,兩圓柱之間只有渦流。其中在L/D≥3.25時,兩圓柱之間出現渦對,導致了所有的流體力和脈動流體力產生突變。這表明,間隙中形成渦對是渦量特性發生脫變的主要原因。

(4)較好地展示了較小漩渦的形成、發展和變化;清晰地展示了間隙中間較大渦對的形成、分裂和融合的過程。并與經典實驗結果吻合較好。說明與有網格方法比較,該計算模型具有較大的優越性,為進一步研究高雷諾數下串聯雙圓柱流體力的特征提供了重要的研究工具。

[1]Mittal S,Kumar V,Raghuvashi A.Unsteady incompressible flows past two cylinders in tandem and staggered arrangements[J].International Journal for Numerical Methods in Fluids,1997,25(11):1315-1344.

[2]Li J,Chambarel A,Donneaud M.Numerical study of laminar flow past one and two circular cylinders[J].Computers&Fluids,1991,19(2):155-170.

[3]Johnson A A,Tezduyar T E,Liou J.Numerical simulation of flows past periodic arrays of cylinders[J].Computational Mechanics,1993,11(5-6):371-383.

[4]Farrant T,Tan M,Price W G.A cell boundary element method applied to laminar vortex shedding from circular cylinders[J].Computers&Fluids,2001,30(2):211-236.

[5]Hall J W,Ziada S,Weaver D S.Vortex-shedding from single and tandem cylinders in the presence of applied sound[J].Journal of Fluids and Structures,2003,18(6):741-758.

[6]Meneghini J R,Saltara F,Siqueira C L R.Numerical simulation of flow interference between two circular cylinders in tandem and side-by-side arrangements[J].Journal of Fluids and Structures,2001,15(2):327-350.

[7]Sharman B,Lien F S,Davidson L.Numerical predictions of low Reynolds number flows over two tandem circular cylinders[J].International Journal for Numerical Methods in Fluids,2005,47(5):423-447.

[8]Ding H,Shu C,Yeo K S.Numerical simulation of flows around two circular cylinders by mesh-free least square-based finite difference methods[J].International Journal for Numerical Methods in Fluids,2007,53(2):305-332.

[9]Jian D,An L R,Jian F Z.Three-dimensional flow around two tandem circular cylinders with various spacing at Re=220[J].Journal of Hydrodynamics,Ser.B,2006,18(1):48-54.

[10]Singha S,Sinhamahapatra K P.High-resolution numerical simulation of low Reynolds number incompressible flow about two cylinders in tandem[J].Journal of Fluids Engineering,2010,132(1):011101.

[11]Pang Jianhua,Zong Zhi.Improving Discrete vortex method for investigation of the fluctuating forces acting on a circular cylinder at subcritical Reynolds number[C]//Proceedings of 3rd International CFD Conference.Dalian,China,2014.

[12]Pang Jianhua,Zong Zhi,Zou Li.A novel vortex scheme with instantaneous vorticity conserved boundary conditions[J].European Journal of Mechanics/B Fluids,2016.

[13]Pang Jianhua,Zong Zhi,Zou Li.Numerical simulation of the flow around two side-by-side circular cylinders by IVCBC vortex method[J].Ocean Engineering,2016.

[14]Kideok R,Zhu B S,et al.Numerical analysis of unsteady viscous flow through a Weis-Fogh-type ship propulsion mechanism using the advanced vortex method[J].ASME Journal of Fluids Engineering,2006,128:481-487.

[15]Blevins R D.Applied fluid dynamics handbook[M].Van Nostrand Reinhold CO,1984.

[16]Igarashi T.Characteristics of the flow around two circular cylinders arranged in tandem:1st report[J].Bulletin of JSME,1981,24(188):323-331.

[17]Zdravkovich M M.The effects of interference between circular cylinders in cross flow[J].Journal of Fluids and Structures,1987,1(2):239-261.

[18]Zdravkovich M M.Flow around circular cylinders:volume 2:applications[M].Oxford:Oxford University Press,2003.

[19]Ljungkrona L,Norberg C H,Sunden B.Free-stream turbulence and tube spacing effects on surface pressure fluctuations for two tubes in an in-line arrangement[J].Journal of Fluids and Structures,1991,5(6):701-727.

Numerical computation model and wake characteristics of two circular cylinders in tandem based on IVCBC vortex method under high Reynolds number

PANG Jian-hua1,ZONG Zhi1,ZHOU Li1,WANG Da-ming2
(1.State Key Laboratory of Structural Analysis for Industrial Equipment,School of Naval Architecture,Dalian University of Technology,Dalian 116024,China;2.Sichuan Explorer Avitation Technology Co.,Chengdu 610001,China)

Based on the characteristic of two circular cylinders in tandem,the numerical computation model under high Reynolds number was established with the characteristics that IVCBC vortex method suits for computing under high Reynolds number.The convergence of IVCBC vortex method was validated by the classical cases.Wake characteristics of double circular cylinders were investingated under Reynolds number 2.5×104and gap ratios 1.1,1.25,2,2.5,3,3.25 and 4.Generation,split and merger of small vortex in wake were shown clearly.The reason of wake mutation and characteristics of flow around two circular cylinders in tandem was explained detailedly.The study shows that wake models are in good agreement with that of classical experimental model,and the vortex pair between two cylinders is the main reason of wake mutaion.At small gap ratios,formation of the small vortex on top of gap and the vibration of flow inside gap is cynchronization with shedding of vortex on the surface of downstream cylinder.It is large advantages of the simulation model compared with the grid methods because the simulation model can show the small vor-tex in wake.It provides an important tool to study further the characteristics of flow force of double cylinders in tandem.

two side by side circular cylinders;IVCBC(instantaneous vorticity conserved boundary conditions)vortex method;wake model;high Reynolds number;mechanism of flow around

O35

A

10.3969/j.issn.1007-7294.2017.10.001

1007-7294(2017)10-1181-09

2017-04-12

龐建華(1985-),男,博士研究生,E-mail:njpjh@sina.com;宗 智(1964-),男,教授,博士生導師。

主站蜘蛛池模板: 国产综合精品一区二区| 亚洲成人www| 亚洲av日韩av制服丝袜| 日韩精品专区免费无码aⅴ| 在线观看免费国产| 久草网视频在线| 中文字幕人妻无码系列第三区| 欧美成人午夜视频免看| 中文字幕在线欧美| 97国产精品视频自在拍| 99国产精品国产| 国产精品v欧美| 亚洲国产日韩在线成人蜜芽| 欧美国产日韩一区二区三区精品影视 | 国内精品视频在线| 免费女人18毛片a级毛片视频| 欧洲av毛片| 国产国产人在线成免费视频狼人色| 国产欧美性爱网| 中文字幕在线观| 蜜臀AVWWW国产天堂| 国产女人综合久久精品视| 成人午夜视频网站| 在线观看欧美精品二区| 永久免费精品视频| 欧美一区二区丝袜高跟鞋| 美女无遮挡拍拍拍免费视频| 日韩一级毛一欧美一国产| 国产一级视频久久| 又粗又大又爽又紧免费视频| 另类欧美日韩| 国产网站黄| 国产成人麻豆精品| 国产欧美视频在线| 国产午夜人做人免费视频中文| 91无码视频在线观看| a级毛片免费播放| 亚洲婷婷六月| 色成人综合| 性喷潮久久久久久久久| 成人国产免费| 欧美精品H在线播放| 高清不卡一区二区三区香蕉| 亚洲丝袜第一页| 亚洲第一中文字幕| 亚洲三级成人| 曰韩人妻一区二区三区| 91免费国产在线观看尤物| 国产欧美日韩专区发布| 亚洲精品无码AⅤ片青青在线观看| 都市激情亚洲综合久久| 美女一区二区在线观看| 欧美有码在线观看| 精品福利网| 久久青草热| 久久99国产综合精品女同| 亚洲国产成人精品无码区性色| 国产成人欧美| 国产福利一区二区在线观看| 人禽伦免费交视频网页播放| 日韩精品成人网页视频在线| 国产成人夜色91| 欧美天堂在线| 草逼视频国产| 国产91无毒不卡在线观看| 性激烈欧美三级在线播放| 欧美激情一区二区三区成人| 亚洲精品老司机| 日本午夜三级| 亚洲无码四虎黄色网站| 欧美成人精品一级在线观看| 国产精品无码AV片在线观看播放| 亚洲无码37.| 呦系列视频一区二区三区| 911亚洲精品| 免费三A级毛片视频| 玖玖精品在线| 国产av一码二码三码无码| 大香伊人久久| 欧美日本在线播放| 婷婷色在线视频| 91福利免费视频|