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

溝槽微結構尺寸對高速列車橫風特性影響研究

2021-11-13 07:19:50王業騰孫振旭鞠勝軍王夢瑩楊國偉
空氣動力學學報 2021年5期
關鍵詞:模型

王業騰,孫振旭,*,鞠勝軍,王夢瑩,楊國偉

(1. 中國科學院力學研究所 流固耦合系統力學重點實驗室,北京 100190;2. 中國科學院大學 工程科學學院,北京 100084)

0 引言

當高速列車在橫風環境中運行時,其氣動性能會顯著惡化,嚴重時就會引發翻車事故。比如,1981年,印度一輛旅客列車在颶風的影響下摔落橋底,致使800多人死亡。因此,系統地研究橫風作用下列車的氣動特性顯得尤為重要。數值模擬已被廣泛應用于各個工業領域[1-3],是高速列車空氣動力學的主要研究方法之一,具有低成本、高效率等優點。Ben通過數值模擬得到了高速列車在不同角度橫風作用下的氣動升力[4]。Justin 等采用不同的湍流模型對列車受側風作用時周圍的流場進行了數值模擬[5]。劉榮等基于氣動力和力矩數值結果,對高速列車臨界傾覆風速進行了研究[6]。公衍軍等建立了三維高架橋-列車耦合模型,并進行了列車速度和風速耦合作用下的數值模擬[7]。

為改善橫風作用下高速列車的氣動特性,學者們開始探索適用的優化方法。目前,常用于高速列車氣動外形優化設計的方法大多圍繞宏觀尺度展開,但由于實際工程中的諸多限制,這類方法帶來的增益效果已接近極限,需要從新的角度出發探尋提升列車氣動性能的有效途徑。因此,有學者試圖從微觀角度出發,探索提高高速列車氣動性能的新途徑。汪久根等設計并研究了不同類型微結構表面的降噪效果[8-9]。苗秀娟發現在列車頂部加肋可以幫助降低表面壓力[10]。孫朋朋等分析了不同類型微結構的設計參數對高速列車氣動阻力的影響[11]。朱海燕等和唐焜等分別研究了具有球窩和凸殼的非光滑表面的減阻機理[12-13]。基于現有研究,微結構主要用于高速列車的減阻降噪。關于微結構如何影響高速列車側風氣動性能,目前尚無相關文獻報道。Browand等在他們的研究中提到,微結構設計可能對車體的穩定性有很強的影響[14]。因此,評價微結構設計參數對高速列車橫風穩定性的影響是一個重要課題,本文從數值模擬出發研究微結構影響列車橫風特性的機理并分析微結構設計尺寸與列車側向力/傾覆力矩之間的影響關系。

實驗常被用來探索和證實科學的規律。然而,隨著科學技術的快速發展,試驗涉及的因素越來越多,各個因素間的關系也越來越復雜,光憑經驗已無法實現預期目標,故出現了試驗設計這一概念。目前,常用的試驗設計方法包括:全面試驗設計、正交試驗設計、均勻試驗設計、回歸試驗設計及混料試驗設計。其中,正交試驗設計是最早引入我國的試驗方法,并在流體力學實驗領域得到了廣泛應用。例如范澤兵等采用正交試驗設計,以總進氣壓力、功率分流、工作方式和節氣門杠桿位置為試驗因素,對發動機穩態性能進行評價[15]。程效銳以進水口寬度、葉片軸向長度、葉片與出水口軸向距離為試驗因素,基于正交設計方法,優化潛水泵葉片的形狀[16]。

本文以390級Pendolino列車的1∶25縮尺模型為原型,通過在頭車頂部增加矩形條,設計了非光滑表面。第1節介紹本文采用的數值算法、計算模型、計算域和邊界條件,并通過與風洞實驗分析比較,驗證數值模擬方法的正確性。第2節從渦量分布和壁面剪應力分布出發,分析微結構改善列車橫風特性的機理。第3節采用正交試驗設計與統計分析結合的方法,分析矩形條帶設計參數與列車側向力/傾覆力矩之間的敏感性,并給出條帶外形設計的優選方案。最后第4節對研究進行總結,列出主要研究成果,并對今后的工作提出建議。

1 數值模擬

1.1 數值算法

本文采用商業軟件STAR-CCM+研究高速列車在橫風作用下的空氣動力學性能。控制方程為三維不可壓縮N-S方程,采用有限體積法離散。對于時間離散,采用二階隱式格式。同時,為了提高計算精度,在空間離散方面,分別采用二階中心差分格式、二階迎風格式對控制方程的粘性項和非粘性項進行離散。此外,基于SSTk-ω兩方程模型的改進延遲分離渦模擬(IDDES)方法將用于模擬表面光滑和粗糙的高速列車周圍流場。

Spalart等在2014年提出了分離渦模擬(DES)方法[17]。DES結合大渦模擬(LES)方法和雷諾平均Navier-Stokes (RANS)方法,可以在保證計算精度的同時減少計算時間。這種方法多用于解決大規模分離流動問題,但不具有通用性,仍然存在一些問題,比如網格敏感性問題會導致非物理分離現象的出現。為了解決這些問題,Mikhail等提出了延遲分離渦模擬(DDES)方法,即在DES模型中加入延遲函數,使得RANS向LES轉換的速度變慢[18]。

但使用DDES進行模擬時,DDES模型中會出現RANS和LES所解得的對數區斜率不同的情況,從而導致對邊界層效應敏感的分離流受到較大干擾。因此,提出了一種改進的延遲分離渦(IDDES)方法,該方法能夠有效處理對數區內同時由RANS和LES求解的網格,大大降低網格相關性。該方法在高速列車橫向風穩定性研究中得到了廣泛應用。例如,Munoz等采用該方法對受橫風作用的高速列車周圍流場進行了數值模擬[19]。通過分析,他們發現這種方法能夠捕獲剪切層中的小渦。Li等利用該方法得到了高速列車在側風作用下的表面壓力和氣動力,與試驗結果吻合較好,表明該方法能夠準確預測列車周圍的平均流場[20]。

1.2 計算模型和計算域

本文采用計算模型源于伯明翰大學的風洞模型,分為光滑模型和粗糙表面模型兩種。研究的橫風場景為風洞中的橫風場景,即列車不動,來流垂直于列車軸線方向。這樣做的目的是一方面方便與風洞實驗結果對標,確定數值方法和網格的準確性;另一方面,基于粗糙模型采用的條帶結構進行參數化研究,通過靈敏度分析確立關鍵設計參數,并確立最佳參數尺寸。研究結果可以為風洞模型下條帶的最佳設計提供理論指導。本文光滑表面列車模型為390級Pendolino列車的1∶25比例模型,包括一節頭車和一節半掛車。為提高計算效率,在不影響計算結果的前提下,對列車轉向架、空調機、受電弓及風擋等區域進行了簡化,且不考慮天線、車門、車窗、底部懸浮架等零部件對列車氣動性能產生的影響。

數值模擬計算外場區域與風洞實驗段一致,如圖1黑框所示。從圖1可以看到,計算模型除了列車模型外,還包括一條單軌軌道及一塊分隔板。值得注意的是,分隔板固定在計算域側壁,且距地面約0.3 m,這樣做的目的是消除地面效應對列車周圍流場產生的干擾。

圖1 計算域Fig. 1 Computational domain

通過在前車頂部的兩個區域添加矩形條,得到一個表面粗糙的模型,如圖2所示。每個區域有15條帶狀,相鄰兩條帶狀之間的距離相等。圖3顯示了矩形條的幾何形狀。所有條帶截面積相同,但不同區域條帶長度不同。區域1的條帶長度為0.12 m,區域2的條帶長度為0.24 m。

圖2 非光滑計算模型Fig. 2 Computational model with a rough surface

圖3 矩形條帶的幾何尺寸Fig. 3 The geometry of a rectangular strip

1.3 邊界條件

邊界條件如圖4所示。除進出口邊界外,其他邊界均規定無滑移壁面條件。速度入口的氣流速度為7.2 m/s,方向沿x軸負向;壓力出口壓力值為0。在計算過程中,列車模型、軌道及分隔板始終保持靜止,并受到具有90°風向角的相對風作用。

圖4 邊界條件Fig. 4 Boundary conditions

1.4 計算網格

本文使用切割體網格生成器進行空間網格劃分,以正交切體六面體網格為主。由于近壁區具有明顯的邊界層效應,為更好地反映車體周圍的流場結構,可以在列車壁面進行棱柱層網格生成,主要為三棱柱網格和金字塔網格。先以1.1的增長比生成了10層棱柱層網格,由此尺寸確定的無量綱壁面距離y+約為1.13,滿足壁面函數的要求。

為捕捉高速列車周圍流場的更多細節并保證計算精度,在車身附近區域進行了網格加密,設置了三個加密區,即圖中區域A2、A3和A4,網格尺寸分別為0.0024 m、0.0048 m和0.0192 m。圖5分別給出了x= 0截面和z= -0.38截面上的網格分布示意圖。在此基礎上,進一步細化了頭車頂部加設條帶區域的網格,如圖中區域A1所示,該區域內網格尺寸為0.0006 m。此外,圖6給出了條帶附近的網格分布示意圖。從局部放大圖可以看到,條帶表面也生成了棱柱層網格,如區域P所示。同時,條帶的幾何外形在網格生成過程中并未遭到破壞,證明所用網格劃分方案確能保證較高質量的計算網格。

圖5 不同截面上的計算網格分布示意圖Fig. 5 Grids on different sections

圖6 條帶區域計算網格局部放大圖Fig. 6 A partially enlarged view of the grids around the strip area

1.5 方法驗證

Hashmi等[21]在伯明翰大學風洞中完成了圖1和圖2所示模型的實驗。他們得到了列車不同回路上某些點的壓力系數Cp,所有回路的位置如圖7所示。無量綱壓力系數Cp的可定義為:

圖7 所有回路位置Fig. 7 The position of all loops

其中,p為局部靜壓,p0為來流總壓,ρ =1.184 kg/m3為來流密度,Vref=7.2 m/s為來流速度。

圖8給出了環E和G上壓力系數數值計算結果與實驗結果的對比。本文采用笛卡爾網格進行空間網格劃分,生成粗糙網格和精細網格兩組網格。前者網格數量為5.216×107,后者網格數量為8.887×107。從圖中可以看出,采用不同網格方案得到的壓力系數分布是一致的,但與實驗結果對比時,明顯看出采用精細網格的方案計算精度更高。因此,所有后續的計算都以與精細網格一致的方式劃分網格。

圖8 回路E和G上數值計算和實驗結果所得到的壓力系數比較Fig. 8 The numerical and experimental results of the pressure coefficient on the loops E and G

1.6 氣動載荷特性

側向力和傾覆力矩是評價高速列車在側風作用下穩定性的重要指標。為簡化分析,對側向力系數Cs和 傾覆力矩系數Cmz進行定義:

其中,Aref是參考面積,也就是列車在x方向的投影面積,取為0.1232 m2[22];Href是車廂高 度,為0.156 m;Fs是側向力,以沿x軸正向為正;Mz是傾覆力矩,力矩 中 心 在(0,0.07155,0)處;ρ為 來 流 密 度,取 為1.184 kg/m3;Vref為來流速度,即7.2 m/s。

圖9給出了橫風作用下光滑表面列車和非光滑表面列車的Cs和Cmz的計算結果。與光滑表面列車模型的比較,非光滑表面列車模型的側向力系數Cs和傾覆力矩系數Cmz分別減小了3.71%和10.56%,表明局部矩形溝槽表面設計在一定程度上確實利于提升高速列車的橫風穩定性。

圖9 光滑表面模型和粗糙表面模型的Cs和Cmz的計算結果Fig. 9 Cs and Cmz on models with smooth and rough surfaces

2 溝槽微結構流場控制機理

將微結構作用于具有復雜幾何外形的高速運動物體表面,繞體氣流的運動規律較之平板邊界層會發生顯著變化。非光滑表面列車能夠減小的側向力和傾覆力矩的原因是小矩形條帶抑制了流動分離,使得氣流能夠以更加穩定的狀態在列車表面運動。為進一步探究微結構減弱側向力和傾覆力矩的機理,本節比較光滑表面列車模型和粗糙表面列車模型的渦量分布與壁面剪應力分布。

2.1 渦量分布比較

選擇以圖7的G回路截面為研究對象,分析車體周圍的渦流特性,圖10分別給出了光滑表面列車模型和粗糙表面列車模型中所選截面上的渦量分布云圖,圖中不同顏色代表不同的渦量幅值范圍。可以看到,不論是光滑表面列車模型還是粗糙表面列車模型,氣流在迎風面拐角處都會形成強度極高的渦結構,這些渦結構將沿約135°的方向向背風面運動。

圖10 渦量分布云圖Fig. 10 Contours of vorticity magnitude at loop D with (a)smooth and (b) rough surfaces

對比圖10(a、b)可以發現,光滑表面列車模型頂部的高強度渦最初是緊貼壁面運動的,而矩形條帶能使其偏離壁面一段距離運動,起到削弱邊界層內湍流變化的作用。同時,條帶的阻斷作用能夠減少由湍流運動引起的瞬時橫流,并促使大型渦結構分裂為較小的渦結構,從而使渦量幅值顯著降低。具有適當高度的條帶對車體周圍的渦結構起到了阻擋作用,使得湍流過渡出現了延遲。當車體周圍的大尺寸渦結構被條帶分解時,邊界層內部的能量交換加劇,使渦結構原有的能量在一定程度上得到衰減,同時渦結構產生的吸力作用也大大減小,這有助于提升列車的橫向穩定性。與粗糙表面列車模型相比,光滑表面列車模型周圍的渦流強度更強且渦核尺寸更大,這是導致其氣動性能較差的關鍵因素。

此外,從圖10還可以看到,由于轉向架和軌道的影響,列車底部來流形成的渦結構也會具有較大強度。當來自頂部和底部的兩股氣流相遇時,列車背風面將出現復雜的渦量分布特性。與光滑表面列車模型進行對比分析后不難發現,粗糙表面列車模型背風面渦結構的數量有所減少且渦量幅值也得到了一定程度上的衰減。造成這種現象的主要原因是,矩形條帶減弱了邊界層內部的湍流變化,使得流經列車頂部的氣流能夠更加穩定地向背風面運動,如此一來氣流交匯時產生的擾動作用也會隨之減小,從而降低背風面流場的復雜程度。

2.2 壁面剪應力

壁面剪應力,是反映車體表面粘性作用的關鍵物理量。圖11為頭車壁面剪切應力分布云圖,圖中不同顏色代表壁面剪應力幅值 τw的不同取值范圍。可以看到,當來流垂直作用于列車側壁時,在迎風面形成的 τw極 小;隨著氣流向頂部運動, τw的值逐漸增大,并在迎風面與頂部的交界處達到最大值;氣流通過流動分離區后, τw又將逐漸減小并趨于穩定。對比圖11(a)和圖11(b)不難發現,在矩形條帶的影響下,列車迎風面拐角處及頂部區域內的 τw值將顯著減小,意味著車體表面粘性作用受到了抑制,其對氣動載荷(比如側向力)造成的影響也會隨之減小,從而達到提升列車橫風穩定性的目的。

圖11 壁面剪應力分布云圖Fig. 11 Contours of the wall-shear stress at (a)smooth and (b)rough surfaces

3 非光滑表面列車橫風性能靈敏度分析

在列車上表面安裝矩形條帶可以抑制氣流流動分離,進而改善列車的氣動特性。當矩形條帶長度方向與橫風作用方向垂直時,其幾何尺寸會直接影響車體周圍的流場特性。

因此,為研究附加矩形條帶幾何參數對高速列車橫風穩定性的影響,本文首先采用正交試驗設計方法生成原始樣本,然后通過計算流體力學(CFD)仿真軟件進行高精度的氣動力/力矩系數計算,最后結合方差分析和極差分析兩種統計分析方法完成矩形條帶幾何設計參數的靈敏度分析。

3.1 正交試驗設計

本文正交試驗設計以列車模型的側向力系數Cs和 傾覆力矩系數Cmz為試驗指標,選取矩形條帶的四個幾何參數L1、L2、W和H作為設計變量,各參數的意義及取值見表1。

表1 試驗因素及其水平Table 1 Experimental factors and their levels

每個試驗因素均取三個值作為其試驗水平。根據試驗因素及其水平,選擇使用L18(38)正交表進行方案設計。需要注意的是,本文的所有設計列數中,前4列為試驗的設計變量列,后4列為隨機誤差列。對每個設計方案分別進行數值模擬,采用的計算方法、計算域、模型位置、網格劃分方式及邊界條件等均與上一節相同。表2列出了所有試驗方案,并在最后兩列給出了各個方案側向力系數Cs和傾覆力矩系數Cmz的計算結果。

表2 試驗方案及數值計算結果Table 2 Cs and Cmz for rough surfaces with different strip shapes

3.2 方差分析

方差分析技術作為一種統計方法可以用來分析不同觀測條件下實驗結果的差異[23-25],評價各因子對試驗指標的影響程度。其基本思想是:將數據的總波動分解為由因素水平變化引起的波動和由誤差引起的波動兩部分,然后構造出F統計量進行F檢驗,以判斷各個因素對評價指標的影響是否顯著。

表3和表4分別給出了側向力系數和傾覆力矩系數的方差分析結果。df表示自由度,SS表示偏方差和,MS表示方差,F表示各因素的方差和與實驗誤差的方差的比值,通常將空列的偏差平方和視為實驗誤差的方差,P表示F統計量對應的顯著性水平。一般情況下,在方差分析中,“顯著”置信水平常取為0.05,而“極顯著”置信水平常取為0.01。

表3 側向力系數C s的 方差分析結果Table 3 Results of the variance analysis for the side force coefficient Cs

表4 傾覆力矩系數C mz的 方差分析結果Table 4 Results of the variance analysis for the roll moment coefficient Cmz

由表3可以看到,與區域1中矩形條帶的長度L1、區域2中條帶的長度L2、條帶寬度W及條帶高度H的F值對應的顯著性水平p均大于0.05,表明這四個設計變量對側向力系數Cs的影響都不顯著。但相比之下,區域2中條帶的長度L2及條帶高度H對側向力系數Cs的影響效果較之另外兩個設計變量要明顯一些。

由表4給出的傾覆力矩系數Cmz的方差分析結果,不難發現,與區域2中條帶的長度L2及條帶高度H的F值對應的顯著性水平p均小于0.05,意味著這兩個設計變量對傾覆力矩系Cmz有顯著影響,而另外兩個設計變量的影響并不顯著。此外,值得一提的是,與條帶高度H的F值對應的顯著性水平p小于0.01,說明條帶高度H對傾覆力矩系數Cmz影響是極顯著的。

基于上述分析結果,按照F值的大小對四個試驗因素的影響效果進行排序:就Cs而言,L2的影響最明顯,H次之,接著是W,最后是L1;對于Cmz,H的影響效果變為明顯,其次是L2,然后是L1,而W的影響效果最弱。因此,在進行矩形條帶外形設計時,應以區域2中條帶的長度L2及條帶高度H為重點優化對象。

3.3 極差分析

為評價各因子對試驗指標的影響程度,可以對數據進行極差分析,其基本思想是:將數據的總波動分解為由因素水平變化引起的波動和由誤差引起的波動兩部分,然后構造出F統計量進行F檢驗,以判斷各個因素對評價指標的影響是否顯著。

極差分析包括計算和判斷兩個步驟,可按如下流程完成分析:

1)計算Kjm、和Rj,并將結果記錄在表5中。Kjm為第j列因素第m水平所對應的試驗指標和;Kˉjm為第j列因素第m水平所對應的試驗指標平均值;Rj為第j列因素水平下的最大平均值和最小平均值的差值,反映了第j列因素水平波動時試驗指標的變化幅度。

表5 極差分析結果Table 5 Results of the range analysis

2)根據極差Rj的大小,可以判斷試驗因素對試驗指標影響的主次順序。Rj越大,說明該因素水平變化時對試驗指標的影響越大。由表5計算結果可知:對于側向力系數Cs,各因素主次順序為L2>H>W>L1;而對于傾覆力矩系數Cmz,各因素主次順序為H>L2>L1>W。因此,L2和H為主要因素,L1和W為次要因素,這與方差分析結論完全一致。

4)為提升高速列車橫風穩定性,兩個試驗指標的絕對值必然是越小越好,也就是說Kˉj1、Kˉj1、Kˉj1需取最大值。此外,由于L1和W是次要因素,可固定為原值成為約束條件,即L1= 0.09 m,W= 0.003 m,則只需根據圖12趨勢曲線確定L2和H的優選水平組合。由圖12(a)可知,對于側向力系數Cs,L2和H的優選水 平 組 合 為L2= 0.18 m、H= 0.00075 m;而 從圖12(b)可以看出,對于傾覆力矩系數Cmz,L2和H的優選水平組合為L2= 0.12 m、H= 0.0012 m。

圖12 因素與指標關系趨勢曲線Fig. 12 Trend curves of the relationship between factors and indicators

5)由于兩個試驗指標單獨分析得到的優選水平組合不一致,可結合綜合平衡法確定最終優選方案。L2對側向力系數Cs的影響大小排第一位,對傾覆力矩系數Cmz的影響大小排第二位,故選擇L2= 0.18 m。H對側向力系數Cs的影響大小排第二位,對傾覆力矩系數Cmz的影響大小排第一位,故選擇H=0.0012 m。因此,初步預測當區域1中條帶的長度L1取0.09 m、區域2中條帶的長度L2取0.18 m、條帶寬度W取0.003 m、條帶高度H取0.0012 m時,能對高速列車側向力系數Cs和傾覆力矩系數Cmz起到較好的優化作用。

所得優選水平組合對應正交試驗設計的方案5,表6給出了原值方案與優選方案的側向力系數Cs和傾覆力矩系數Cmz的計算結果。優選方案的條帶高度較之原值方案增加了60%。可以看到,與原值方案相比,優選方案的側向力系數Cs的值增加了0.80%,而傾覆力矩系數Cmz的值減少了6.14%。由于傾覆力矩系數Cmz是造成列車傾覆的最主要原因,從綜合性能來看,該優選方案確實具有提升列車橫風穩定性的作用。

表6 原值方案與優選方案計算結果對比Table 6 Comparison of the original and optimized schemes

4 結論

本文采用基于SSTk-ω兩方程模型的IDDES方法,論證了矩形條對受側風作用的高速列車氣動性能的影響,并結合正交試驗設計,采用方差分析法和極差分析了矩形條幾何參數的影響。得到如下結論

1)與光滑模型相比,粗糙模型可使列車側力系數和傾覆力矩系數分別降低3.71 %和10.56 %,從而提高列車的橫風穩定性。

2) 矩形條帶可以抑制氣流流動分離,促使大尺寸渦結構分解為眾多小尺寸渦結構。在此過程中,邊界層內部的能量交換加劇,使渦結構中的能量得到衰減,減小了近壁區內的湍流擾動和車體表面粘性作用,保證低速氣流能夠更加穩定地發展,以達到提升列車橫向穩定性的目的。

3)采用方差分析方法和極差分析可以準確分析矩形條的幾何參數對橫風作用下高速列車氣動性能的影響。在利用矩形條帶參數優化進行高速列車橫風氣動外形設計時,需重點關注兩個主要參數,即區域2中條帶的長度L2和條帶高度H,而區域1中條帶的長度L1和條帶寬度W為次要因素,對列車橫風性能的影響很小,在非必要時可以忽略。

4)通過極差分析得到了矩形條帶的優選參數設計方案,即區域1中條帶的長度L1取0.09 m、區域2中條帶的長度L2取0.18 m、條帶寬度W取0.003 m、條帶高度H取0.0012 m。數值計算結果顯示,與原值方案相比,優選方案的側向力系數Cs增加了0.80%、傾覆力矩系數Cmz減少了6.14%,表明該方案確能有效改善高速列車氣動性能。

上述結論可為后續高速列車側風氣動優化設計提供新的思路。然而,目前的研究主要集中在風洞實驗中的縮比模型上,由于風洞中列車模型靜止,因而本文研究橫風作用與實際運行場景下的橫風作用存在一定差異。微結構表面在實際運行中如何影響高速列車的氣動性能尚不確定,這將是我們今后研究的主要目標。

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 四虎成人免费毛片| 日韩高清中文字幕| 思思99热精品在线| 在线观看视频一区二区| 亚洲人成人无码www| 欧美精品亚洲精品日韩专区| 无码中字出轨中文人妻中文中| 人人91人人澡人人妻人人爽| 国产午夜一级毛片| 国产精品hd在线播放| 欧美激情伊人| a级免费视频| 国产精品久久久久鬼色| 亚洲欧美自拍中文| 国产精品任我爽爆在线播放6080| 精品1区2区3区| 中文字幕在线观看日本| 午夜毛片免费看| 欧美国产在线精品17p| 国产综合另类小说色区色噜噜| 强奷白丝美女在线观看| 99热最新网址| 亚洲Av综合日韩精品久久久| 精品视频福利| 在线观看精品自拍视频| 亚洲第一黄色网| 午夜精品一区二区蜜桃| 精品丝袜美腿国产一区| 国产在线98福利播放视频免费| 伊人久久婷婷五月综合97色 | 成AV人片一区二区三区久久| 亚洲AV电影不卡在线观看| 91免费观看视频| 欧美人人干| 国产综合欧美| 国产女同自拍视频| 69国产精品视频免费| 国产欧美成人不卡视频| 欧美人与牲动交a欧美精品| 久久久久国产精品嫩草影院| 99精品在线视频观看| 亚洲人成色在线观看| 欧美日本视频在线观看| 色综合成人| 国产91高跟丝袜| a级毛片免费看| 国产成人精品一区二区免费看京| 亚洲视频二| 亚洲综合极品香蕉久久网| 色爽网免费视频| 91小视频在线| 91青草视频| 亚洲美女一级毛片| 天堂av高清一区二区三区| 日本亚洲国产一区二区三区| 国产乱子伦手机在线| 四虎在线高清无码| 久久五月天综合| 四虎亚洲国产成人久久精品| 亚洲伊人电影| 亚洲全网成人资源在线观看| 亚洲综合色婷婷| 国产在线观看第二页| 欧洲成人在线观看| 18禁黄无遮挡免费动漫网站| www亚洲天堂| 日本高清免费不卡视频| 日韩一区精品视频一区二区| 亚洲成AV人手机在线观看网站| 日韩一区精品视频一区二区| 国产亚洲高清在线精品99| 香蕉eeww99国产精选播放| 欧美激情视频二区三区| 久99久热只有精品国产15| 国产精品无码久久久久久| 亚洲国产在一区二区三区| 国产最新无码专区在线| 国产男女XX00免费观看| 成人91在线| 国产香蕉在线| 国产手机在线ΑⅤ片无码观看| 欧美成人午夜视频|