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

串列雙旋轉圓柱群繞流特性及互擾效應研究

2023-09-22 01:48:32涂佳黃黃林茜肖毅凡吳蘇莉彭世中
船舶力學 2023年9期

涂佳黃,黃林茜,肖毅凡,2,吳蘇莉,彭世中,吳 楷

(1.湘潭大學土木工程學院,湖南 湘潭 411105;2.湖南建工集團有限公司,長沙 410000;3.長沙市無線電監測站,長沙 410013;4.湖南省無線電監測站,長沙 410011)

0 引 言

鈍體繞流現象廣泛存在于實際工程中,如海洋工程、土木工程與機械工程等。當流體流經鈍體結構時,由于流體黏性力的存在,導致鈍體周圍會出現剪切層的分離,并伴隨回流現象的產生與旋渦周期性脫落。圓柱繞流作為流體力學中一個經典而復雜的課題,引起了國內外學者的關注,相關研究取得了一系列的成果[1-3]。同時,對如何優化結構表面所受流體力的相關研究工作也已展開[4-5]。

關于旋轉單圓柱繞流,國內外學者已開展了一系列實驗與數值模擬方面的研究[6]。Kang等[7]對低雷諾數下旋轉單圓柱繞流進行了實驗研究,發現隨著轉速的增加,圓柱沿流體流動方向的受力逐漸減小,當超過臨界轉速后,流場渦脫落被抑制;Kumar等[8]發現當圓柱旋轉速度與流體流動速度的比值大于1.95 時,圓柱渦脫落被抑制,當比值達到4.34~4.7 時,渦脫落現象再次出現;Lu 等[9]對雷諾數的變化對旋轉單圓柱尾流特性與流體力參數的影響進行了研究。另一方面,Rao等[10]對單旋轉圓柱三維繞流問題進行了系統研究,分析了多個參數對流體力與流場特性的影響,揭示了其內在機理。何穎等[11]發現亞臨界雷諾數條件下旋轉運動可以有效抑制其旋渦脫落,且隨著轉速的增加其升阻比大幅度地提高;同時,程友良等[12]發現超臨界雷諾數條件下,圓柱平均阻力系數變化的臨界旋轉速率為2。

相比較于單圓柱工況,多圓柱工況由于結構之間以及結構與流體之間的相互干擾作用使得其流動性能以及尾流場情況變得更加復雜。雙圓柱繞流系統其流動影響的因素主要包括雷諾數、排列方式、間距以及旋轉運動參數等[13-15]。對于串列布置旋轉雙圓柱工況,陶實等[16]對Re=100條件下串列旋轉雙圓柱進行了數值模擬研究,分析了間距比與旋轉速度對流場結構的影響以及雙圓柱的渦脫落形態和升阻力特性;巴悅等[17]分析了不同間距以及旋轉方式對雙圓柱氣動性能以及流場結構的影響,發現上下游圓柱均順時針旋轉且間距為3 時為較優的排列旋轉方式;吳新等[18]對Re=100 條件下不同旋轉方式對并列雙圓柱繞流特性的影響進行了分析,發現逆流向旋轉比順流向旋轉對渦脫落的抑制效果更好;Darvishyadegari 等[19]發現Re=200 條件下串列布置圓柱間豎向流線呈現“V”型,隨著旋轉速度的增大,圓柱上側流線與下側流線發生錯位現象,使得流線呈現“之”型;Behara 等[20]對Re=150 與Re=2000工況下,串列布置旋轉三圓柱的流致振動問題進行了研究,分析了圓柱旋轉速度與折減速度對圓柱群動力響應及其周圍流場的影響。

綜上所述,目前對于雙旋轉圓柱繞流問題的研究主要集中于并列布置、二維工況,關于不同的參數對串列布置雙旋轉結構三維流場特性以及機理的影響研究相對較少,存在一系列的問題需要深入探索。本文運用Fluent軟件,對Re=200工況下,串列布置旋轉雙圓柱三維繞流問題進行數值模擬,探究不同間距比(L/D=1.5~5.0)與旋轉速率(α=0.0~4.0)對上、下游圓柱升阻力系數、尾流特性及流場流動形態的影響,揭示上、下游旋轉圓柱之間的互擾效應,對實際工程初步設計階段具有一定的參考價值。

1 計算理論及模型驗證

1.1 基礎理論

不可壓縮黏性流體的連續性方程與動量方程為

本文采用層流模型進行計算,數值模擬采用壓力與速度的耦合應用SIMPLE 算法,采用二階迎風格式對控制方程進行離散、二階隱式格式對時間項進行離散,連續性方程與動量方程收斂余差均小于10-5。另一方面,圓柱旋轉速率為α=Dω/(2U),其中ω為旋轉角速度。阻力和升力系數計算公式為

式中,Cdp、Cdf與Clp、Clf分別代表阻力與升力的壓力和粘性分量,Fd和Fl分別是施加在圓柱上的總阻力和總升力。

1.2 計算模型與邊界條件

本文計算域尺寸設定為40D×20D×4D,D為圓柱體直徑,如圖1(a)所示。上游圓柱圓心距流體進口與上、下邊界距離均為15D,距流體出口距離為30D。上、下游圓柱直徑D=1 m,n=L/D為間距比,L為兩圓柱截面圓心之間的距離,如圖1(b)所示。流場入口設置為均勻來流,流速U=1 m/s,流體粘度μ=0.02~0.005 Pa·s,流體密度恒定ρ=1 kg/m3。流場入口邊界設置為速度入口,即ux=U=1 m/s,uy=uz=0;流場出口邊界設置為壓力出口,即p=0;壁面邊界設置為無滑移壁面條件,即ux=uy=0。上、下游圓柱表面均設定為無滑移壁面,且均以逆時針方向相同速度旋轉。

圖1 計算模型與旋轉方向示意圖Fig.1 Schematic diagram of calculation model and rotation direction

1.3 計算參數選取與驗證

本文采用結構化網格對計算流域進行網格劃分,靠近圓柱體壁面的網格進行局部加密,能更好地刻畫邊界層處的流場。同時,圓柱體尾流區域進行局部加密處理,較遠處的網格較為稀疏,如圖2 所示。另外,邊界層第一層網格高度小于0.03,時間步長取為0.001 s。

圖2 網格劃分示意圖Fig.2 Schematic diagram of the meshing

為了確保該數值模型的精確性,對串列雙旋轉圓柱繞流數值模型進行了網格無關性驗證。表1為四套不同加密網格計算得到的圓柱阻力系數均值,根據ITTC 推薦規程選取網格加細比為。由于Mesh1、Mesh2網格數量較稀疏,使得值相較于其他網格模型差異較大,從而導致計算結果產生偏差。Mesh3 與Mesh4 的計算結果較為一致,綜合考慮選用Mesh3 網格模型。本文時間步長取為0.001 s,且滿足y+小于1,邊界層第一層高度小于0.03。作為與雙圓柱繞流問題對比的依據,本文對Re=200 條件下的單圓柱繞流問題進行了數值計算,發現本文計算的結果與文獻實驗數據以及模擬結果吻合良好,見表2。

表1 不同網格下,串列雙圓柱阻力系數平均值變化情況Tab.1 Mean drag coefficient of a series of two cylinders at different mesh cases

表2 Re=200工況下,靜止單圓柱繞流驗證Tab.2 Verification of flow around a stationary single cylinder at Re=200

本文進一步對靜止雙圓柱流場分布特性進行了驗證。根據Zdravkovich[25]的實驗,臨界間距一般為3 倍圓柱直徑,因此選取L/D=3.0 與L/D=4.0 兩個間距比工況的流場瞬時分布圖進行對比驗證。由圖3 可知,當小于臨界間距比時(L/D=3.0),此時上游圓柱沒有尾渦脫落,兩圓柱之間由間隙渦填滿。下游圓柱下表面的尾渦發育完整,而上表面的尾渦已形成回流區。當間距比超過臨界間距比(L/D=4.0)時,上游圓柱體結構尾流區的漩渦有了一定的泄渦空間,上下游圓柱都產生了周期性的渦脫落。這與文獻[21]結果基本一致。

圖3 Re=200條件下串列靜止雙圓柱繞流的流場瞬時分布圖Fig.3 Instantaneous distribution diagram of the flow field around a series of stationary two cylinders under Re=200

2 計算結果與分析

本文對不同間距比(1.5≤L/D≤5.0)與旋轉速率(0.0≤α≤4.0)工況下串列布置雙旋轉圓柱群三維繞流問題進行數值模擬研究,主要分析了串列雙旋轉圓柱體結構群繞流的流場分布特性及其互擾效應機理。

2.1 流場結構

結構表面剪切層的發展與相互作用及尾渦結構的分布會影響流場分布特性,導致柱體表面所受的壓力分布發生變化,進一步改變圓柱體表面受到的流體力。本節對α=0~4.0,L/D=2.0 與3.0 工況下串列雙旋轉圓柱流向與展向瞬時三維渦量圖進行分析。其中,流向與展向渦度分量定義為:ωx=?w/?y-?v/?z,ωz=?v/?x-?u/?y,其中u、v和w分別是x、y和z方向的速度分量。圖中紅色表示正等值面,藍色表示負等值面。

圖4 與圖5 分別為L/D=2.0 時串列雙旋轉圓柱在不同轉速下流向渦與展向渦的瞬時三維渦量圖。當α=0 時,由于兩圓柱處于靜止狀態且間距較小,從圖4(a)可觀察到兩圓柱表面沒有渦流的覆蓋,只在下游圓柱尾部形成正負相間具有對稱性的流向渦,由于流向渦的對稱性使得此時兩圓柱上受到的升阻力系數均值為0。同時,展向渦主要由上游圓柱產生而將下游圓柱包裹在內,且由于剪切層受到強烈的拉伸作用,使得兩個剪切層在下游同向卷起,如圖5(a)所示。隨著雙圓柱開始旋轉,尾部渦流對稱性被破壞,分布變得紊亂起來,當圓柱旋轉速率較小(α=1.0)時,兩圓柱表面逐漸有渦流的覆蓋,其尾渦分布逐漸變得紊亂但大致呈對稱型分布,因此兩圓柱受到的升、阻力系數均值仍舊為零。然而展向渦受到圓柱旋轉的影響逐漸分離撕裂,剪切層卷起位置較α=0工況更靠近下游圓柱,如圖5(b)所示。隨著旋轉速率的進一步增大,尾渦的形成受到抑制,會出現破碎重組,最后轉捩為混合交錯排列的小尺度結構,形成正負相間的小渦排列,從上游圓柱脫落的流向渦填充滿整個間隙區域將下游圓柱包裹在內。當α=4.0時,兩圓柱間隙區剪切層沿切向速度方向彎曲,最終纏繞在圓柱體表面,形成包圍整個圓周的渦層,且下游圓柱尾部的展向渦不再卷起,如圖5(d)所示。

圖4 L/D=2.0時,在不同旋轉速率下流向渦|ωx|=1的瞬時渦量等值面圖Fig.4 Instantaneous iso-surfaces of streamwise vorticity|ωx|=1 of cylinder at L/D=2.0 with different rotating rates

圖5 L/D=2.0時,在不同旋轉速率下展向渦|ωz|=1的瞬時渦量等值面圖Fig.5 Instantaneous iso-surfaces of spanwise vorticity|ωz|=1 of cylinder at L/D=2.0 with different rotating rates

圖6 與圖7 分別為L/D=3.0 時串列雙旋轉圓柱在不同轉速下流向渦與展向渦的瞬時三維渦量圖。此時兩圓柱同時受到臨近效應與剪切層干擾效應作用。較之于L/D=2.0工況,上游圓柱尾流對下游圓柱影響加大。從圖6(a)觀察到,此時兩圓柱靜止時其表面便有大量的渦流覆蓋,在下游圓柱尾部形成的流向渦比L/D=2.0工況要更大,且流向渦越往后行寬度越大。當圓柱開始旋轉時(α≤2.0),由于圓柱間距較大,流向渦從上游圓柱表面脫落后,充斥滿整個間隙區,并與下游圓柱流向渦相互耦合,在尾流區后方較規律地卷起,隨著轉速的增大,兩圓柱之間的耦合作用越發增強,流向渦結構分布越發紊亂,導致上、下游圓柱升力系數增大。隨著圓柱的旋轉,流場中上、下兩個剪切層之間的相互作用加強,正負展向渦相互混合,在下游圓柱表面與尾流區均出現了撕裂現象,如圖7(b)與7(c)所示。當α=4.0時,流向渦結構中的渦瓣區進一步收縮變小,尾流區流向渦分布變得比較規則,而由于展向渦卷吸拉伸,導致流場中負向展向渦會在兩側包圍住正向展向渦,使得雙圓柱升力系數均值降為零。

圖6 L/D=3.0時,在不同旋轉速率下流向渦|ωx|=1的瞬時渦量等值面圖Fig.6 Instantaneous iso-surfaces of streamwise vorticity|ωx|=1 of cylinder at L/D=3.0 with different rotating rates

圖7 L/D=3.0時,在不同旋轉速率下展向渦|ωz|=1的瞬時渦量等值面圖Fig.7 Instantaneous iso-surfaces of spanwise vorticity|ωz|=1 of cylinder at L/D=3.0 with different rotating rates

2.2 壓力分布特性

圖8 為不同間距比工況下,串列雙旋轉圓柱體的時均壓力與流線分布情況,選取的截面z/H=0.5。圖中紫色區域表示正壓,灰色區域表示負壓。當圓柱處于靜止狀態時,時均壓力系數沿上游圓柱表面呈對稱分布,圓柱體前壁面受到來流施加的正壓,其后壁面受到回流產生的負壓,而下游圓柱由于受到上游圓柱尾流的影響導致其表面時均壓力系數分布不對稱。當圓柱開始旋轉時,圓柱前壁面的正壓區往圓柱左上方移動,負壓區面積增大,且旋轉速度越大,負壓面積越大。隨著間距比的增大,圓柱周圍流動特征會發生相應的變化,諸多學者對間距比變化的影響進行了研究,本文得到的流場隨間距比的變化規律與文獻[21,26]結果相符合。

圖8 在不同間距比與旋轉速率下,串列雙旋轉圓柱的時均壓力分布與流線圖Fig.8 Time-average pressure distribution and streamline of double-rotating cylinders at different space ratios and rotation rates

當間距比較小時(L/D=1.5 與2.0),兩圓柱之間的干擾效應以臨近效應為主。當α=0 時,上游圓柱迎流面受來流沖擊作用,在其迎流面形成正壓區,由于圓柱不旋轉使得來流在上游圓柱頂端分離,并在其尾部形成中間小、兩頭寬的兩個對稱梯形渦。從上游圓柱分離的剪切層再附著在下游圓柱上,兩者間隙區域形成穩定的回流區,此時兩圓柱體形成一個擴大的“單圓柱體”,渦脫落現象只在下游圓柱出現。由于間隙渦的存在以及下游圓柱脫落的尾渦使得上游圓柱背流面以及下游圓柱整個圓周表面分布的時均壓力為負值。當圓柱開始旋轉時,流體粘性力的存在使得圓周表面流體會隨著圓柱的旋轉而加速,由于兩圓柱均為逆時針旋轉,使得圓柱下表面來流速度與圓柱旋轉速度方向相同而流速疊加,圓柱上表面來流速度與圓柱旋轉速度方向相反而總流速下降,造成兩圓柱下表面負壓區面積增大,從而抑制了圓柱下表面渦流的形成。且由于圓柱的旋轉使得來流邊界層與上游圓柱分離點向圓柱左上方轉移,上游圓柱正壓區也向左上方偏移,隨著圓柱旋轉速度的增大,正壓區面積越來越小直至消失。當α=4.0時,此時由于圓柱旋轉速度過快,兩圓柱表面完全由負壓主導,抑制了尾部渦流的形成,且尾流角度(即尾流與來流方向的角度)變化呈現較大差異。

當間距比較大時(L/D=3.0與4.0),兩圓柱之間的干擾效應以臨近效應和剪切層干擾為主。當α=0時,由于兩圓柱之間的間距增大,使得上游柱的尾渦拉長,填充滿整個間隙區域。在L/D=3.0 工況下,下游圓柱背流面的尾渦急劇收縮,上側形成尾渦,下側形成回流區。然而,L/D=4.0 工況剛好相反,其下游圓柱尾部上側形成回流區,下側形成尾渦,其時均壓力分布情況與小間距比工況相似。當α=1.0時,在L/D=3.0工況下,上游圓柱僅在背流面上側形成尾渦,而下游圓柱在其迎流面與背流面上方各形成一個小渦。當間距比進一步增大到4.0 時,兩圓柱之間的臨近干擾效應微弱,使得下游圓柱迎流面的渦消失,兩圓柱均只在其尾部形成一個上側渦。下游圓柱表面負壓區面積隨著間距比的增大而逐漸減小,此時兩圓柱主要受到升力作用。隨著旋轉速率的進一步增大,上游圓柱的尾渦向上偏移,下游圓柱尾流中的尾渦向上偏移同時收縮變小。隨著間距比的增大尾渦結構逐漸拉長,最終形成回流區。當α=4.0時,兩圓柱尾部均沒有旋渦形成,其表面主要承受負壓作用。

2.3 流體力參數特性

本節對不同間距比工況下,串列雙旋轉圓柱流體力系數隨旋轉速率的變化規律進行分析。由圖可知,上游圓柱與下游圓柱的阻力系數平均值隨旋轉速率α的增加其變化規律相反。隨著α的增大,上游圓柱Cˉd值減小,其平均阻力系數均為負值,且間距比越小其下降幅度越大,最大由0.62 下降到-7.96(L/D=1.5),如圖9(a)所示。當旋轉速率較小時,對不同間距比工況Cˉd值影響較小,從α>2.0開始,隨著旋轉速率的增大,不同間距比工況受到的影響加大。由圖10可知,當α=1時,壓差產生的阻力Cdp為負值,粘性阻力Cdf為正值,兩種成分所占比例均衡,使得平均阻力較小且隨間距比的變化影響不大。當α≥3.0時,上游圓柱的Cdp值均為正值,Cdf值均為負值,隨著圓柱旋轉速率的增加,阻力中的組成部分均會增大,但粘性阻力占主導地位,如圖10所示,從而導致上游圓柱阻力系數平均值為負值。下游圓柱Cˉd值則隨α的增大而增大,其平均阻力系數基本為正值,這是由于雙圓柱逆流向旋轉時,上游圓柱的升力方向為垂直向上,下游圓柱的升力方向為垂直向下,導致雙圓柱之間存在斥力,如圖9(b)所示。相較于上游圓柱,下游圓柱的Cdp值均為負值,Cdf值均為正值,且粘性阻力占主導地位,這也是下游圓柱阻力系數平均值為正值的原因,如圖11所示。

圖9 不同間距比工況下,串列雙旋轉圓柱流體力系數隨旋轉速率的變化Fig.9 Variations of fluid force coefficients with rotation rate of tandem double-rotating cylinder under different space ratios

圖10 不同旋轉速率工況下,上游圓柱平均阻力系數隨間距比的變化Fig.10 Variation of mean drag coefficient of the upstream cylinder with the space ratios under different rotation rates

圖11 不同旋轉速率工況下,下游圓柱平均阻力系數隨間距比的變化Fig.11 Variation of mean drag coefficient of the downstream cylinder with the space ratios under different rotation rates

與阻力系數平均值相比,上游圓柱與下游圓柱在各個間距比工況下的升力系數均方根值隨旋轉速率α變化的整體規律大致相同,如圖9(c)與9(d)所示。小間距比工況(L/D=1.5 與2.0)下,隨著α的增大,Cl'值變化不大,在0.0~0.2之間微弱波動。然而,L/D=3.0與4.0工況下,Cl'值隨著α的增大呈現先增后減的明顯趨勢。不同的是,當α=0 時,上游圓柱Cl'=0,而下游圓柱由于上游圓柱尾流的影響Cl'=0.2。當L/D=5.0時,兩圓柱的升力系數均方根值在α=0.0工況達到最大,隨著α的增大,逐漸減小到0。整體上看,由于上游圓柱尾流的影響,使得下游圓柱升力系數均方根值大于上游圓柱。由圖12 與圖13 觀察到,相比較于阻力系數平均值,升力系數平均值的變化規律一致。隨著L/D增大,壓差升力為正,Clp逐漸增大,然而,粘性升力為負,Clf逐漸減小,兩者變化幅度基本相同,使得升力系數平均值受間距比的影響較小。另一方面,隨著α的增大,升力系數平均值的影響顯著,其值會大幅度減小,其中Clp值逐漸增大,Clf值逐漸減小。

圖13 不同旋轉速率工況下,下游圓柱平均升力系數隨間距比的變化Fig.13 Variation of mean lift coefficient of the downstream cylinder with the space ratios under different rotation rates

2.4 泄渦頻率特性

圖14 給出了不同工況下串列雙旋轉圓柱的斯托羅哈數(St)的變化情況。由圖可知,隨著α的增加,上、下游圓柱的St變化規律趨于一致,但兩者的渦脫落頻率并非相同,與文獻[21]的結論一致。當旋轉速率較小(α<2.0)時,間距比對雙圓柱泄渦頻率有較大影響,St值在0.05~0.10 之間浮動。當α≥2時,間距比的影響逐漸減弱,雙圓柱的泄渦頻率均會趨于0。同時,隨著間距比的增大,St變化的速度會加快。由于圓柱旋轉速率增加,使得圓柱周圍流速加快,從而導致圓柱表面主要受負壓作用,產生較大吸力,并抑制尾流渦脫落,如圖7所示。

圖14 不同間距比工況下,串列雙旋轉圓柱斯托羅哈數隨旋轉速率的變化Fig.14 Variations of Strouhal number with rotation rates of tandem double-rotating cylinder under different space ratios

3 結 論

本文對串列雙旋轉圓柱三維繞流問題進行了數值模擬研究,探究了不同間距比和旋轉速率對旋轉圓柱流體力系數、表面壓力系數、流場特性以及渦脫落的影響,并揭示了旋轉雙圓柱體的互擾效應及其內在機理。研究得出以下結論:

(1)間距比較小工況下,隨著α的增大,流向渦會從大渦結構轉變為混合交錯排列的小尺度結構。展向渦則主要由上游圓柱產生而將下游圓柱包裹在內,并在遠尾流處同向卷起。當L/D≥3.0 時,由于上游圓柱尾流影響的增強,使得流向渦比L/D=2.0 工況時更大更寬,而展現渦則會出現撕裂與融合現象。

(2)由于圓柱旋轉,使得圓柱尾流呈現不對稱形態,尾渦向下側發生偏移。隨著L/D增大,圓柱對轉速的敏感度增加。上、下游圓柱時均壓力分布受間距比的影響較小,受轉速的影響較大。隨著轉速的增加,圓柱表面負壓區面積增大,形成較大吸力,從而抑制尾流渦的產生。

(3)不同間距比工況下,隨著α的增大,上、下游圓柱阻力的變化規律相反,而升力變化規律則相同。另外,當α≥2.0 時,上、下游圓柱的升力系數均方根值均會逐漸趨于0,且間距比增大會加快變化的速度。

(4)隨著α的變化,上游圓柱的阻力系數成分的受力方向及大小均會發生改變,導致阻力系數平均值發生異號。然而,對于下游圓柱,會影響其阻力系數成分的受力大小,在小間距比工況下阻力系數平均值會發生異號。同時,上下游圓柱的升力系數成分的受力大小也會發生改變,導致其平均值的大小改變。

(5)上、下游圓柱斯托羅哈爾數(St)隨α變化的規律類似,當旋轉速率較大(α≥2)時其均會趨于0。同時,隨L/D增加,St趨于0的變化速度會加快。

主站蜘蛛池模板: 毛片在线播放a| 最新国产你懂的在线网址| 丰满人妻久久中文字幕| 40岁成熟女人牲交片免费| 免费国产福利| 精品自窥自偷在线看| 国产sm重味一区二区三区| 精品无码人妻一区二区| 在线人成精品免费视频| 综合色在线| 久久99这里精品8国产| 色首页AV在线| 日韩精品无码免费专网站| 日韩AV无码免费一二三区| 91久久国产热精品免费| 五月婷婷导航| 91精品国产无线乱码在线| 亚洲福利视频一区二区| 亚洲无码精彩视频在线观看| 色婷婷狠狠干| 亚洲国产综合精品一区| 国产一区二区三区夜色| 一级看片免费视频| 伊人大杳蕉中文无码| 99久久精品免费看国产电影| 亚洲自拍另类| 国产成人精品2021欧美日韩| 天天综合天天综合| 久久久久无码精品| 亚洲性日韩精品一区二区| 色香蕉影院| 免费va国产在线观看| 亚洲日韩AV无码精品| 91亚洲影院| 亚洲一区精品视频在线| 国产视频你懂得| 亚洲精品无码在线播放网站| 成人在线亚洲| 欧美翘臀一区二区三区| 国产导航在线| 四虎永久在线精品影院| 国产亚洲欧美在线专区| 五月天天天色| 国产 在线视频无码| 一区二区三区四区精品视频 | 亚洲中文字幕久久精品无码一区| 国产精品毛片一区| 国产迷奸在线看| 色亚洲激情综合精品无码视频| 中文字幕一区二区人妻电影| 青青久久91| 亚洲日韩久久综合中文字幕| 久久网欧美| 亚洲第一天堂无码专区| 午夜福利网址| 精品自拍视频在线观看| 日韩不卡高清视频| 日韩AV手机在线观看蜜芽| 99一级毛片| 国产成人乱码一区二区三区在线| 国产爽爽视频| 青青青国产视频| 国产va欧美va在线观看| 久久情精品国产品免费| 亚洲成a人在线播放www| 亚洲高清无码久久久| 日本在线免费网站| 欧美精品综合视频一区二区| 欧美激情第一区| 国产成人91精品| 亚洲第一成年人网站| 正在播放久久| 蜜桃臀无码内射一区二区三区| 亚洲精品爱草草视频在线| 久久这里只精品国产99热8| 亚洲天堂网在线播放| 亚洲国产精品成人久久综合影院| 国产成年女人特黄特色大片免费| 中文毛片无遮挡播放免费| 精品五夜婷香蕉国产线看观看| 怡红院美国分院一区二区| 亚洲人成网线在线播放va|