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

受電弓下沉對其氣動和聲學行為的影響

2022-11-01 01:34:52戴志遠
中國機械工程 2022年20期

秦 登 戴志遠 周 寧 李 田

西南交通大學牽引動力國家重點實驗室,成都,610031

0 引言

隨著高速列車運行速度的提高,空氣動力學問題成為制約高速鐵路發展的重要因素[1-3]。受電弓周圍復雜的氣流流動會導致滑板和接觸線之間發生振蕩和產生較強的氣動噪聲。已有研究表明,當列車運行速度達到300 km/h以上時,空氣阻力占列車總阻力的85%以上,受電弓氣動阻力占總阻力的15%以上[4]。此外,氣動噪聲在列車時速高于300 km時成為主導噪聲,轉向架、車間連接、頭車鼻尖、受電弓等均為氣動噪聲的主要來源[5-8]。高速列車減阻降噪這個多目標優化問題隨著速度的增大變得更加重要。

為了減小高速列車的運行阻力和氣動噪聲,目前采用的優化措施主要有:增大頭車流線型長度和改變頭車流線型外形[9-10];設置半裙板或全裙板優化轉向架[11];采用半包或全包形式的外風擋優化車端連接處[12]。已有眾多學者對高速受電弓的氣動特性和聲學特性進行了試驗和數值研究[13-14]。受電弓主要由桿件組成,優化桿件的橫截面形狀可以有效減小空氣阻力和氣動噪聲[15]。現有受電弓滑板橫截面大多近似為矩形設計,研究表明,設計成流線型橫截面時,受電弓氣動阻力減小約40%左右,氣動噪聲減小2~5 dBA[16]。盡管受電弓的結構復雜,但通過其他合理設計也可顯著減小其氣動阻力和噪聲。通過增設流線型導流罩、受電弓桿件表面覆蓋多孔介質或進行表面微結構設計均可減阻降噪[17-19]。

隨著列車運行速度的進一步提高,近年來有了新的優化方案。在列車車頂設計一個下沉式的平臺安裝受電弓,與受電弓直接安放于車頂上相比,列車總氣動阻力減小1%~4%,同時研究表明矩形和五邊形的下沉平臺比其他平臺更好[20]。許多研究人員對轉向架和轉向架腔的流動行為和聲學特性進行了較為深入的研究[21-22],但是對下沉平臺所形成的受電弓腔的流動行為和氣動噪聲特性研究較少。并且,形成的受電弓艙很難通過某種手段對其產生的噪聲輻射進行屏蔽,安裝平臺下沉對受電弓氣動性能和遠場氣動噪聲性能的影響規律還有待研究。

受電弓下沉安裝方式目前已較為普遍,但是對受電弓艙產生的氣動噪聲關注卻很少,腔體與受電弓部件之間的相互作用仍不清楚,因此,本文基于數值計算的方法,研究安裝平臺的下沉高度和傾角對高速受電弓流動行為和聲學行為的影響。

1 數值計算方法

1.1 改進延遲分離渦模擬(IDDES)

受電弓由多個部件組成且大部分為桿件,其尾渦流場具有非常復雜的三維流動結構,包含了各種尺度的漩渦。直接數值模擬是對受電弓周圍湍流描述最準確的方法,由于模型的復雜性,目前的計算條件還難以實現。非定常雷諾時均模擬(unsteady Reynolds-averaged Navier-Stokes, URANS)、大渦模擬(large eddy simulation, LES)以及分離渦模擬(detach eddy simulation, DES)是目前廣泛采用的湍流數值求解方法[23]。URANS方法具有計算量小的優點,但它很難準確模擬出流體的瞬時脈動。LES方法可以捕捉到URANS方法無法捕捉的許多非穩態、非平衡過程中出現的大尺度效應和擬序結構,但是完全捕捉到附面層內擬序結構的演化需要非常精細的網格。

DES方法屬于雷諾時均(RANS)/LES混合類方法,結合了RANS和LES方法的各自特點,減少了網格數量和計算量。SPALART等[24]在Spalart-Allmaras(SA)模型的基礎上提出了SA-DES方法,這種方法對網格間距有很大的依賴性,可能會導致邊界層內部從RANS提前轉換為LES,并導致模化應力損耗現象。針對這個問題,SPALART等[25]又提出了延遲分離渦模擬(delay detach eddy simulation, DDES),通過修改長度尺度的定義以確保RANS模式求解的邊界層與網格分辨率無關。MENTER等[26-27]也提出了基于SST(shear stress transfer)模型的DES方法來解決模化應力損耗問題。之后又提出了改進延遲分離渦模擬(improve delay detach eddy simulation, IDDES)方法,該方法包含較多的流場上游和湍流信息,在尺度構造中包含更多的近壁信息以防止RANS和LES交界處出現雷諾應力過度下降和產生對數層不匹配現象[28]。

部分研究發現SSTk-ω和IDDES方法在模擬壓力系數分布、描述流場精細結構等方面均能取得較好的結果[29],因此,本文采用SSTk-ω湍流模型求解高速列車周圍的定常流場,之后采用基于SST湍流模型的IDDES方法求解非定常流場。

1.2 FW-H方程

Ffowcs Willams和Hawkings運用廣義函數理論推導出了在靜止流體中做任意運動的控制面的發聲方程,簡稱FW-H方程[30],其積分形式為

(1)

式中,t為時間;H(f)為Heaviside函數,f=0表示有限固定表面邊界區域,f>0表示無限空間邊界區域,H(f)=0表示在積分曲面內,H(f)=1表示積分曲面在外;δ(f)為Dirac delta函數;2即Δ表示拉普拉斯算子;ρ′為流體密度的脈動量,ρ′=ρ-ρ0,ρ0為流場中未受擾動的密度,ρ為流場中的密度;xi、xj為直角坐標分量(i,j=1,2,3);c0為流場聲速;Q為流體與控制面相互作用導致的單極子噪聲源項,Q=ρ0vn+ρ(un-vn),un為垂直于積分表面的流體速度分量,vn為表面本身的法向速度;Fi為流體與控制面相互作用導致的偶極子噪聲源項,Fi=pijnj+ρui(un-vn),pij為應力張量,nj為積分曲面的法向單位向量;Tij為流場中湍流以及漩渦產生的四極子噪聲源項,分別為流體速度在xi、xj方向的速度分量,p′為流體脈動壓力,δij為克羅內克符號(當i=j時,δij=1;當i≠j時,δij=0),τij為黏性應力張量。

式(1)等號右邊三部分聲源項分別對應于厚度聲源、載荷聲源以及四極子聲源。高速列車的運行速度一般小于367.5 km/h(0.3 Ma),對應受電弓的氣動噪聲主要為偶極子噪聲[31],因此,本文只考慮偶極子噪聲項對高速受電弓氣動噪聲的貢獻。

2 數值計算模型

2.1 幾何模型、計算區域和邊界條件

受電弓主要分為弓頭(滑板、弓頭支撐、弓角等)、上框架、下臂桿、拉桿、平衡桿、底架和絕緣子等幾部分,對受電弓部件附屬的細小結構進行了一定的簡化,如螺栓等,簡化后的幾何模型如圖1所示。

圖1 高速受電弓模型Fig.1 High-speed pantograph model

為了減少網格計算量且反映受電弓下沉對列車氣動阻力的影響,建立了三編組的真實列車車體模型,但忽略了轉向架和轉向架艙。計算區域如圖2所示,計算區域長L=320 m,寬W=60 m,高H=20 m。計算域的設置滿足列車空氣動力學規范中的最低要求。列車位于計算區域的中心位置,計算域前表面距列車頭車鼻尖的距離約為80 m;后表面距尾車鼻尖的距離約為160 m。兩個側面距離列車縱向中心截面的距離均為30 m。車底距離地面0.376 m。受電弓安裝平臺如圖3所示,長l=3.424 m,寬w=2.354 m,高h分別為0、0.1、0.2、0.3、0.4、0.5 m;車頂與受電弓安裝平臺的過渡傾角α分別為0°、30°、45°、60°、90°。

圖2 計算域(m)Fig.2 Calculation domain(m)

(a)三視圖

(b)俯視圖

(c)正視圖圖3 受電弓安裝平臺(mm)Fig.3 Pantograph installation platform(mm)

邊界條件設置如下:入口設定為速度入口,速度大小為350 km/h,出口設定為壓力出口,兩個側面和頂部設為對稱邊界,列車和受電弓表面設置為無滑移固定壁面,地面設置為滑移地面,滑移速度大小與列車運行速度相等,方向相反。當來流方向為正y時,受電弓為開口運行。若來流方向為負y方向,受電弓為閉口運行。

2.2 網格劃分及獨立性檢驗

劃分3套不同尺度的網格開展網格獨立性檢驗,計算結果如表1所示。由表1可得,網格1和網格2計算得到的整車時均氣動阻力之間的誤差小于0.5%。圖4進一步給出了3套計算網格下,離受電弓后滑板不同距離處的流場時均速度幅值。可見,在時均氣動阻力和流場物理量上網格2已經滿足獨立性要求。后續計算選擇網格2。

表1 網格獨立性檢驗

圖4 時均速度幅值對比Fig.4 Comparison of time average speed amplitude

網格劃分的基礎尺寸為512 mm,列車表面最大網格尺寸為64 mm,受電弓表面最大網格尺寸為8 mm。為保證網格過渡良好,在列車和受電弓周圍進行了多次局部加密,加密區位置如圖2所示,外流場網格劃分結果如圖5a所示。為捕捉附面層流動,在受電弓和列車表面劃分邊界層網格。一般劃分邊界層網格時與所采用的湍流模型息息相關,不同的湍流模型對邊界層網格的第一層高度有所限制。理應保證邊界層過渡均勻,每一層網格之間的跨度不宜過大,邊界層網格層數應保證最后一層網格位于邊界層之外。選取邊界層的增長率為1.1,劃分15層,滑板上表面邊界層劃分結果如圖5b所示。計算網格數量約為5675萬。

(a)外流場網格劃分

(b)滑板上表面邊界層圖5 網格劃分Fig.5 Meshing

2.3 求解設置和噪聲測點布置

利用SSTk-ω湍流模型和IDDES方法分別求解高速列車周圍的定常和非定常流場。利用FW-H聲學類比方程求解遠場噪聲。為了保證流場計算的準確性,壓力-速度耦合采用SIMPLEC算法,時間項采用二階隱式方案,對流項采用二階迎風格式離散,擴散項采用中心差分格式離散。非定常流場和聲學計算的時間步長均選取為5×10-5s,當非定常流場達到穩定后,提取受電弓以及受電弓艙表面的聲源數據用于遠場氣動噪聲輸入,計算總時長為0.3 s,最大分析頻率為10 kHz。

氣動噪聲性能評估以遠場噪聲測點的聲壓級大小和頻譜特性來衡量,根據高速列車噪聲測試國際標準ISO 3095-2013相關要求,測點布置如圖6所示。距離軌道中心線25 m、距離軌道3.5 m高,沿列車縱向(y向)布置31個測點(yp1~yp31),每個測點之間相距1 m。受電弓幾何中心在地面的投影坐標為(0, 0, 0)m。橫向布置8個噪聲測點(xp1~xp8),距受電弓幾何中心的距離分別為7.5 m、12 m、15 m、25 m、30 m、40 m、50和60 m,其中xp4和yp16噪聲測點重合。

圖6 遠場噪聲測點示意圖Fig.6 Schematic diagram of far-field noisemonitoring points

3 數值計算方法驗證

在中國空氣動力研究與發展中心低速空氣動力研究所的風洞群中,分別開展了某型高速受電弓的氣動力試驗[32]和高速列車的氣動噪聲試驗[33]。圖7和圖8所示分別為現場試驗模型。氣動力測試選用1/1的受電弓模型,氣動力通過三臺常溫五分量天平測得。高速列車噪聲試驗選用1/8的縮比模型,在模型側面5.8 m遠處布置了30個遠場傳聲器用于測量試驗模型的遠場噪聲。

圖7 氣動試驗模型Fig.7 Aerodynamic test model

圖8 聲學試驗模型Fig.8 Acoustic test model

圖9給出了受電弓氣動力試驗與數值計算的對比結果。開口、閉口狀態下,數值模擬得到的整弓氣動阻力Fd與風洞試驗數據基本吻合,誤差小于5%。圖10給出了風速200 km/h時受電弓區域內氣動噪聲在數值模擬和試驗之間的頻譜比較。為了確保聲壓級不受干擾,數值模擬與風洞試驗采用相同的頻率帶寬。由圖9、圖10可得,數值模擬和試驗展現出較好的一致性。表2進一步給出了遠場30個噪聲測點的算術平均聲壓級對比結果,試驗與仿真之間的誤差小于1.5 dBA。考慮到數值仿真存在一定的誤差且試驗存在隨機誤差等,本文所采用的數值計算方法能保證一定的準確性和可靠性。

圖9 受電弓氣動阻力對比Fig.9 Comparison of aerodynamic resistanceof pantograph

圖10 受電弓頻譜對比Fig.10 Comparison of the spectrum of the pantograph

表2 受電弓氣動噪聲對比

4 計算結果

4.1 氣動行為

4.1.1下沉高度的影響

圖11給出了受電弓的下沉形式,考慮到列車實際運行時接觸網高度一定,平臺下沉高度通過受電弓升弓來補償。受電弓工作高度(軌道上表面到滑板頂部的距離,5300 mm)不變,受電弓下沉高度通過增大升弓高度來補償,受電弓初始升弓高度為1290 mm。

圖11 受電弓下沉方式Fig.11 Subsidence methods of pantograph

圖12給出了受電弓下沉300 mm和500 mm情況下的表面時均壓力分布云圖。受電弓氣動阻力的主要來源是前后表面的壓力差。受電弓安裝平臺下沉后改變了受電弓的下層流場結構,與安裝平臺無下沉相比,隨著下沉高度的增大,列車前方高速氣流作用在底架和絕緣子前部的正壓明顯減小。

表3給出了安裝平臺不同的下沉高度下,受電弓開口運行時各部件時均氣動阻力的變化情況。安裝平臺的下沉高度對弓頭和上框架氣動阻力的影響較小。平衡桿和拉桿的氣動阻力隨安裝平臺下沉高度的增大略有增大,主要原因是以升弓方式補償受電弓的下沉高度會導致這些部件打表4進一步給出了受電弓安裝平臺不同下沉高度下,受電弓開口運行時高速列車各部件的時均氣動阻力。由表4可得,頭車的氣動阻力基本不受受電弓安裝平臺下沉的影響。在受電弓安裝平臺下沉后,尾車以及不含受電弓艙的中間車體部分的氣動阻力受下沉高度變化的影響較小。受電弓安裝平臺下沉高度增大,所形成的受電弓艙產生的氣動阻力先增大后減小。總之,受電弓安裝平臺下沉后,雖然會顯著減小受電弓的氣動阻力,但形成的受電弓艙也會產生額外的氣動阻力。安裝平臺下沉高度h在100~500 mm范圍內,高速列車整車的氣動阻力先增大后減小,但均為增阻狀態,阻力增大1.11%~4.44%。當h=300 mm左右時,整車阻力增加最大。

(a)升弓高度1290 mm,h=0 (b)升弓高度1590 mm,h=300 mm (c)升弓高度1790 mm,h=500 mm圖12 受電弓表面的時均壓力分布Fig.12 Time-averaged pressure distribution on the surface of the pantograph

表3 受電弓開口運行時各部件的氣動阻力

表4 受電弓開口運行時高速列車的氣動阻力

表5給出了受電弓以閉口方式運行時,高速列車各部件的時均氣動阻力分布。由表5可得,各部件氣動阻力的變化規律基本同受電弓開口運行時一致。安裝平臺下沉高度h在100~500 mm范圍內,受電弓以閉口方式運行時,高速列車整車阻力增阻1.05%~3.85%。

圖13給出了受電弓安裝平臺不同下沉高度下,縱向中心截面的壓力和流線分布情況。由圖13可得,受電弓安裝平臺下沉改變了絕緣子和底架周圍的流場結構。隨著下沉高度的增大,腔體前端面對受電弓下層結構的遮擋效果增強。受電弓下層結構的阻力會顯著減小。此外,隨著下沉高度的增大,腔體內的流動也發生顯著的改變。當下沉高度h在100~300 mm時,形成的腔體較淺,腔體前端面后方和后端面前方分別形成負壓和正壓區域,這導致腔體產生了較大的氣動阻力。

表5 受電弓閉口運行時高速列車的氣動阻力

(a)開口下沉h=0 (b)開口下沉h=100 mm (c)開口下沉h=200 mm

(d)開口下沉h=300 mm (e)開口下沉h=400 mm (f)開口下沉h=500 mm圖13 下沉高度對受電弓周圍流場的影響Fig.13 The effect of subsidence height on the flow field around the pantograph

當下沉高度進一步增大到500 mm時,腔體較深,氣流在腔體前后端面附近形成較大尺度的漩渦,后端面受到的正壓作用顯著減小。因此,隨著下沉高度的增大,腔體所導致的阻力先增大后減小。當受電弓以閉口運行時,氣流發展與開口運行時類似。

4.1.2傾角的影響

由上文分析可知,h=300 mm時整車氣動阻力增加最大。圖14給出了安裝平臺下沉300 mm時,不同腔體傾角下受電弓周圍的流場。由圖14可得,隨著傾角的減小,受電弓下層結構受到氣流的沖擊作用增強,但腔體前端形成的剪切氣流強度減弱,腔體前端面后方的漩渦消失且負壓區域明顯減小。同時腔體后端面受到氣流沖擊形成的正壓也減小,因此,減小腔體的過渡傾角可減小受電弓艙的氣動阻力。

圖15給出了高速列車各部件氣動阻力隨腔體傾角變化的結果。由圖15可得,頭車、尾車以及未加受電弓艙的中間車氣動阻力受傾角變化的影響較小。隨傾角的減小,受電弓氣動阻力逐漸增大,受電弓艙和整車的氣動阻力逐漸減小。當傾角由90°減小為30°時,受電弓開口、閉口運行整車氣動阻力分別減小4.8%和3.9%。與安裝平臺不下沉相比,受電弓開口、閉口運行時其氣動阻力分別減小2.0%、1.8%,整車氣動阻力分別減小1.4%、1.1%。因此,選取合適的下沉高度和腔體過渡傾角可使整車氣動阻力達到最好的減阻效果。

(a)開口傾角30° (b)開口傾角45° (c)開口傾角60°

(d)閉口傾角30° (e)閉口傾角45° (f)閉口傾角60°圖14 傾角對受電弓周圍流場的影響Fig.14 The effect of inclination on the flow field around the pantograph

(a)開口運行

(b)閉口運行圖15 傾角對高速列車氣動阻力的影響Fig.15 The effect of inclination on the aerodynamicresistance of high-speed trains

(b)h=300 mm, α=30°圖16 受電弓周圍渦量分布(Q=1800)Fig.16 Vorticity distribution around the pantograph(Q=1800)

4.2 聲學行為

4.2.1噪聲源特性

(a)h=0, α=0°

圖16為基于Q準則的受電弓區域渦量分布對比圖,受電弓周圍區域存在尺度和旋向不同的漩渦。弓頭、鉸接、底架以及絕緣子區域是漩渦脫落和重組的主要部位,是受電弓的主要氣動噪聲源,受電弓安裝平臺下沉后,絕緣子和底架周圍氣流速度顯著降低,其后部漩渦脫落強度和尺度顯著減小。由于腔體的存在,受電弓下層渦旋的發展也受到一定的限制。

圖17給出了受電弓表面的聲功率級分布。由圖17可得,受電弓安裝平臺下沉后,由于流場的改善,絕緣子和底架表面的聲功率顯著減小。絕緣子和底架區域的表面聲功率級最大減小可達50 dB。可見,受電弓下沉后可顯著降低其下層結構的噪聲輻射強度。

(a)h=0 mm, α=0°

(b)h=300 mm, α=30°圖17 受電弓表面聲功率級Fig.17 Surface acoustic power of pantograph

4.2.2遠場噪聲特性

圖18給出了兩種工況下橫向噪聲測點的聲壓級對比結果。由圖18可得,受電弓艙產生的氣動噪聲較小,與受電弓相比相差10 dBA以上。受電弓下沉后所產生的噪聲輻射強度顯著降低。與安裝平臺無下沉相比,安裝平臺下沉后受電弓下層結構(底架和絕緣子)的遠場聲壓級最大減小4.48 dBA,整弓遠場聲壓級最大減小2.61 dBA。

圖18 橫向測點聲壓級對比Fig.18 Comparison of sound pressure levels at lateralmonitoring points

圖19 縱向噪聲測點聲壓級對比Fig.19 Comparison of sound pressure levels oflongitudinal noise monitoring points

圖19給出了兩種工況下,受電弓縱向噪聲測點的聲壓級對比。與受電弓安裝平臺無下沉相比,下沉后受電弓下層結構的遠場聲壓級最大減小4.25 dBA,整弓遠場聲壓級最大減小2.02 dBA。

為了綜合表示受電弓遠場氣動噪聲的性能,以平均聲壓級對受電弓的氣動噪聲進行評估,根據能量疊加原理,平均聲壓級(sound pressure level)Lp的計算公式為

(2)

式中,Lpi為第i個噪聲評估點的聲壓級,i=1,2,…,m,m為噪聲測點的總數。

表6給出了兩種工況下受電弓縱向噪聲測點的平均聲壓級結果。將受電弓分為三大部分,分別為上層結構(弓頭等),中層結構(上框架、下臂桿、拉桿和平衡桿等)和下層結構(底架和絕緣子等)。安裝平臺下沉對受電弓上層結構的噪聲影響較小,對下層結構的影響較大。下層結構減小4.06 dBA,整弓遠場平均聲壓級減小1.31 dBA。

表6 遠場平均聲壓級對比

4.2.3頻譜特性

圖20給出了平臺下沉高度h=0、傾角α=0°下,yp16噪聲測點處的功率譜密度。由圖20可知:受電弓具有明顯的主頻特性,主要頻率為330 Hz左右,受電弓上層結構(弓頭)對主頻的貢獻量最大。受電弓下層和中層結構對低頻能量貢獻較大,上層結構對中頻能量貢獻較大。

圖20 yp16噪聲測點的功率譜密度Fig.20 Power spectral density at noise monitoring point yp16

圖21給出了不同工況下,yp16噪聲測點處功率譜密度的對比。由圖21可知,兩種工況下yp16噪聲測點處功率譜密度的變化規律一致,安裝平臺改善了絕緣子和底架周圍的流場結構。與受電弓安裝平臺無下沉相比,當受電弓下沉300 mm時,遠場噪聲測點處的噪聲能量更低。

圖21 yp16噪聲測點功率譜密度對比Fig.21 Comparison of power spectral density at noisemonitoring point yp16

圖22 yp16噪聲測點處的1/3倍頻程對比Fig.22 Comparison of 1/3 octave frequency at noisemonitoring point yp16

圖22給出了兩種工況下,yp16噪聲測點處的A計權1/3倍頻程。受電弓遠場氣動噪聲的能量集中在400~2500 Hz范圍內。與受電弓安裝平臺無下沉相比,受電弓安裝平臺下沉后降低了受電弓低中頻范圍內的噪聲輻射,減小了受電弓在遠場噪聲測點處的聲壓級,改善了受電弓的氣動噪聲性能。

5 結論

(1)受電弓前后表面的壓力差是氣動阻力的主要來源。安裝平臺下沉后改變了受電弓的下層流場結構,絕緣子和底架受高速氣流沖擊的迎風面積減小致使氣動阻力減小。隨著安裝平臺下沉高度的增大,受電弓氣動阻力逐漸減小。受電弓艙會產生額外的氣動阻力,需通過減小腔體傾角來緩和氣流使附加的阻力減小。當受電弓安裝平臺下沉300 mm且腔體傾角為30°時,受電弓開口、閉口運行整車分別減阻1.4%、1.1%。

(2)受電弓區域氣動噪聲的主要來源有弓頭、框架鉸接、底座和絕緣子等。安裝平臺下沉后,受電弓下層結構周圍流體流速顯著降低,絕緣子和底座表面聲功率級大幅度減小。

(3)受電弓向外輻射噪聲具有明顯的主頻特性,主要頻率約為330 Hz,主要能量集中在400~2500 Hz范圍。受電弓艙所附加的氣動噪聲相比于受電弓氣動噪聲更小。受電弓艙通過改善流場提高了受電弓的聲學性能,降低了其噪聲輻射能量。與安裝平臺無下沉相比,下沉300 mm且腔體傾角為30°的情況下,受電弓向外輻射的聲壓級最大減小2.02 dBA,平均聲壓級減小1.31 dBA。

主站蜘蛛池模板: h视频在线播放| 一区二区三区四区精品视频| 国产国模一区二区三区四区| 亚洲an第二区国产精品| 国产精品免费福利久久播放| 67194在线午夜亚洲| 国产欧美自拍视频| 亚洲精品麻豆| 中文字幕自拍偷拍| 亚洲天堂精品视频| 天天综合亚洲| 国产91精品久久| 精品国产电影久久九九| 亚洲欧美综合精品久久成人网| 日韩小视频在线播放| 成年人午夜免费视频| 免费一级毛片完整版在线看| 性视频久久| 国产小视频免费观看| 国产美女无遮挡免费视频| 国产资源免费观看| 亚洲国产中文精品va在线播放| 久久永久视频| 婷婷综合在线观看丁香| 欧美日韩一区二区三区四区在线观看 | 一本二本三本不卡无码| 国产av剧情无码精品色午夜| 国产性生大片免费观看性欧美| 色屁屁一区二区三区视频国产| 777午夜精品电影免费看| 茄子视频毛片免费观看| 日本福利视频网站| 国产极品美女在线播放| 精品超清无码视频在线观看| 无码又爽又刺激的高潮视频| 国产香蕉在线视频| 国产小视频免费| 日韩在线第三页| 丰满人妻被猛烈进入无码| 黄色国产在线| 日本高清在线看免费观看| 国产成人综合在线视频| 有专无码视频| 日本欧美午夜| 456亚洲人成高清在线| 欧美一级黄色影院| 国产成人免费视频精品一区二区| 国产男人的天堂| 亚洲天堂啪啪| 国产激情在线视频| 国产欧美在线观看精品一区污| 亚洲欧美极品| 国产高潮视频在线观看| 国产99视频免费精品是看6| 成人精品视频一区二区在线| 亚洲欧美极品| 亚洲免费福利视频| 欧美精品v日韩精品v国产精品| 一本大道无码高清| 国产精品999在线| 日韩在线欧美在线| 伊人91视频| 美女啪啪无遮挡| 老司机aⅴ在线精品导航| 免费国产无遮挡又黄又爽| 成人午夜久久| 亚洲伊人天堂| 亚洲熟女中文字幕男人总站| 日韩欧美国产中文| 久久99国产视频| 国产91小视频| 亚洲精品无码AⅤ片青青在线观看| 亚洲无码高清视频在线观看| 不卡视频国产| 成人一级免费视频| 中文成人在线视频| 亚洲三级片在线看| 国产超薄肉色丝袜网站| 亚洲人网站| 女同久久精品国产99国| 操操操综合网| 不卡的在线视频免费观看|