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

高速鐵路復線隧道內壓力波特性的數值模擬研究

2015-03-05 02:35:12劉希賢,郭安寧,梅元貴
鐵道科學與工程學報 2015年1期
關鍵詞:高速列車

?

高速鐵路復線隧道內壓力波特性的數值模擬研究

劉希賢,郭安寧,梅元貴,許建林,李士安,周朝暉

(蘭州交通大學 機電工程學院,甘肅 蘭州 730070)

摘要:采用有限體積方法和任意滑移界面動網格技術的CFD方法,對我國高速列車和復線隧道條件下的壓力波基本特性進行數值模擬研究。為避免數值計算產生不合理的物理現象,應用光滑啟動技術。在驗證本文數值方法正確性的基礎上,采用三維非定常可壓縮湍流流動模型方法,研究高速列車駛入單洞復線隧道時的壓力波基本特征,揭示隧道內壓力波的形成過程,得到不同斷面上不同測點的壓力變化及分布規律;研究結果表明:在初始壓縮波三維效應的影響下,近隧道入口斷面不同測點壓力差異較大,且隨著列車速度的增加,正壓最大值也越大;在距離出口較近的斷面上,當列車速度增加到一定范圍,正壓最大值反而減小;隧道襯砌上的正壓最大值多出現在距離隧道入口較近處,而負壓最大值則出現在隧道的中央斷面上。

關鍵詞:高速列車;復線隧道;壓力波;數值模擬

隧道壓力波是高速鐵路隧道空氣動力學中的一個重要現象,會產生乘客舒適性、列車車體及部件和隧道內固定設備等疲勞強度設計問題[1-4]。隨著我國高速鐵路的發展,國內開展了實車和動模型試驗,揭示了壓力波的形成和變化特性[5-10]。隧道壓力波的數值模擬方法有一維流動模型和三維流動模型方法。一維流動模型方法能合理地提供隧道斷面上的平均壓力波特性,并取得了試驗驗證[1-2,11]。三維流動模型能給出較詳細的洞內壓力和流速特性,為洞口和列車外形設計等提供依據[5,10,12-14],并用于列車運行、洞內輔助結構等方面的空氣動力特性分析[15-17]。同時,還出現了一維和三維流動模型耦合方法,用于分析長大隧道及隧道群的空氣動力效應[18]。高速鐵路隧道大多采用單洞復線的形式。列車和隧道形成的幾何空間為非對稱空間,引起的流動也是非對稱的,作用于不同隧道襯砌表面的壓力特性也不同,這就提出了合理確定其氣動載荷的問題。對于此類問題,目前公開的研究文獻尚沒有見到報道。因此,本文基于我國高速列車CRH380A,采用基于CFD軟件的三維非定常可壓縮湍流流動模型,在與發展成熟可靠的一維可壓縮非定常模型的特征線求解結果比較基礎上,進行隧道壓力波的形成機理、隧道斷面上壓力分布特性等方面的研究,期望為今后的研究提供一定的基礎。

1問題描述和數值計算方法

1.1幾何模型和計算區域

本文將單洞復線隧道簡化為光滑半圓管模型,忽略隧道內的建構造物,同時不考慮隧道外的地形環境,隧道的橫斷面直徑D=13.3 m,橫截面面積FTU=100 m2,長LTU=58.6D,線間距0.38D;列車模型采用光滑表面的CRH380A高速列車簡化氣動模型,忽略轉向架和受電弓等部件,軌面以上車高ht=0.28D,長LTR=7.74D;車/隧模型的阻塞比為0.108。本文采用1∶10的縮尺模型。計算區域由隧道空間及其端口外的開闊空間構成,兩者橫截面為同心圓弧,開闊空間半徑15D。以隧道入口平面內隧道中心線與軌面的交點為坐標原點,列車前進方向為x正方向,鉛垂向上為z正方向,按右手準則建立直角坐標系。圖1為計算區域和簡化的列車模型。圖2為隧道內含測點的斷面位置和測點布置示意圖。在結果分析中,定義高速列車到達隧道入口端的時刻為t=0。

1.2流動模型及定解條件

高速列車通過隧道引起的空氣流動是三維非定常可壓縮湍流流動。流動控制方程有連續性方程、動量方程和能量方程。本文采用的湍流模型和壁面處理模型分別為SST k-ω高雷諾數湍流模型和非平衡壁面函數法[19]。高速列車通過隧道時的邊界區域如圖3所示,其中隧道表面和地面為無滑移靜止壁面,伸入計算區域流體的隧道部分為無滑移擋板邊界,兩側均為無滑移靜止壁面,高速列車表面為無滑移運動壁面。這些固體壁面均為絕熱邊界。隧道外開闊空間的各表面為常溫常壓的黎曼邊界。為避免列車模型突然啟動產生的不合理物理現象,應用了光滑啟動技術。在求解初始時刻,高速列車處于靜止狀態,其周圍流場也靜止,按海平面國際標準大氣條件確定參考壓力和溫度狀態,湍流物理量處處為0。

圖1 計算區域和高速列車幾何模型Fig.1 Sketch of computational domain and train model

(a)縱斷面測點位置;(b)橫斷面測點位置圖2 隧道內典型斷面測點位置示意圖Fig.2 Sketch of pressure monitors in tunnel

1.3網格劃分策略

應用STAR-CD軟件的ASI技術求解高速列車通過隧道的問題,將整個計算區域分為靜止和運動兩大部分,兩者之間的交界面構成滑移界面,動靜網格通過“耦合—求解—解耦—運動—再耦合,實現滑移網格的計算。在網格劃分中,將動、靜2部分作為2個單獨的“體(Body)”處理,同時將列車周圍區域的網格從運動網格中單獨分出來作為一個“體(Body)”,以進行適度加密。這樣形成3個體的多塊網格劃分方法,其中列車周圍區域完全包含于運動部分之中,兩者之間的交界面用網格耦合處理。網格劃分采用軟件ICEM-CFD,網格總數約為420萬。其中,高速列車表面法向第1層網格的厚度按照y+=50取定,流向的網格尺度為Dx+=9~152y+,展向的網格尺度為Dz+=5~52y+。計算區域體網格的尺度依照光順過渡原則向四周延伸,拉伸比為1.2~2。圖4為列車在隧道內運行的面網格劃分結果示意圖。

圖3 網格分布規劃與多體網格劃分策略Fig.3 Layout of grid distribution and the strategy of multi-body grid generation

圖4 列車和隧道表面網格示意圖Fig.4 Surface mesh of the train and the tunnel

1.4數值計算方法

采用壓力修正算法中的PISO算法求解列車通過隧道引起的流動控制方程。各方程時間項的離散格式為歐拉隱格式,擴散項的離散格式采用中心差分格式,對流項的離散格式采用MARS格式。代數方程組的求解方法為代數多重網格法(AMG)。

在非定常問題的求解中,盡管隱式時間離散格式是無條件穩定的,但實際計算中為減小誤差累積,應滿足穩定性條件,即

式中:Co為Courant數;|U|=ucell+a,ucell為某點速度,a為當地聲速;l為網格特征尺度。用列車表面最小網格尺度和固體壁面法向第1層網格厚度估算時間步長Dt,其典型值為0.018LTR/Vtr。

2計算結果和討論

2.1數值方法驗證

這里采用已發展成熟可靠的一維可壓縮非定常不等熵廣義黎曼變量特征線方法的計算結果[12],驗證本文三維流動模型方法的正確性和合理性。隧道和列車模型與圖1保持一致,一維流動模型計算參數采用原型尺度。列車長度為114 m,速度為350 km/ h,隧道長度為1 140 m。

(a)距隧道入口570 m處的測點壓力時間歷程曲線;(b)頭車車身測點的壓力時間歷程曲線圖5 列車過隧道產生的壓力波Fig.5 Pressure wave generated by train passing through the tunnel

圖5為本文三維流動模型與一維流動模型計算結果的對比。由圖5可見,2種流動模型的計算結果吻合較好,說明本文基于軟件的三維流動模型和數值模擬方法能夠較好地模擬隧道壓力波問題。

2.2壓力波的形成

根據圖2選取的測點斷面獲取壓力波數據,研究列車過復線隧道的壓力波特性。圖6~圖8分別為距隧道入口2D斷面、距隧道入口12.5D斷面和距出口2D斷面的中央頂部測點“6”的壓力和相應斷面上平均壓力的時間歷程曲線。其中,對于圖7,t=0.157 2 s為列車車頭鼻尖到達x=12.5D斷面時刻。由圖可見:這3個斷面的平均壓力和測點“6”的時間歷程曲線基本重合,單點測得的壓力和斷面上平均壓力基本相同。同時,從圖7得出:當壓縮波達到該斷面時,測點上壓力及斷面平均壓力均升高,膨脹波經過時,測點壓力及斷面平均壓力均降低。

圖6 隧道內距離入口2D處測點及斷面平均壓力的比較Fig.6 Comparison of presssure change at section 2D from the entry portal and their average pressure

圖7 隧道內距離入口12.5D處測點及斷面平均壓力的比較Fig.7 Comparison of presssure change at section 12.5D from the entry portal and their average pressure

圖9為列車通過隧道全過程的壓力云圖。從圖中可以清晰看出:壓縮波和膨脹波的產生及傳播過程,以及在隧道端口的反射,波系強度的衰減變化的特性。

圖8 隧道內距離出口2D處測點及斷面平均壓力的比較Fig.8 Comparison of presssure change at section 2D from the exit portal and their average pressure

圖9 列車通過隧道全過程的壓力云圖Fig.9 Pressure contour of the whole process the train passing through tunnel

2.3不同斷面上測點壓力差異的比較

圖10為隧道內3個不同斷面上不同測點的壓力時間歷程曲線。表1~表3分別統計了3個斷面上的壓力峰值,表中以遠離運行列車的隧道地面測點“1”的壓力峰值為基準,得出其他測點與其差值的百分比。其中,列車車頭在t=0.025 2 s到達距入口x=2D斷面,t=0.369 s到達中央斷面,t=0.712 8 s到達距隧道出口x=2D斷面。由圖表可見:當列車車頭、車尾經過測點斷面時,斷面上不同位置測點的壓力不同。

特別是距入口x=2D斷面,由于其距隧道入口較近,當t=0.025 2 s,列車頭部鼻尖通過時,斷面上測點的壓力僅受到列車進入隧道產生的初始壓縮波影響,在此處這一期間壓縮波具有三維效應,導致不同測點壓力差異較大;當t=0.034 9 s,列車頭部肩部通過;t=0.113 s,列車尾部肩部過;t=0.122 7 s,列車尾部鼻尖通過,由于端部曲線幾何變化,也導致了不同測點的壓力差異。而在其他時刻,不同測點的壓力值幾乎一樣,并且同一斷面不同測點的壓力變化趨勢基本一致。

(a)距離隧道入口2D斷面上不同測點的壓力時間歷程曲線;(b)隧道中央斷面上不同測點的壓力時間歷程曲線;(c)距離隧道出口2D斷面上不同測點的壓力時間歷程曲線圖10 隧道內不同斷面上測點的壓力時間歷程曲線Fig.10 Pressure history on different cross-sections of the tunnel measuring points

對于其他2個斷面上不同測點的壓力差異,也主要體現在列車頭部鼻尖與肩部、列車尾部肩部與鼻尖通過時帶來的差異。列車端部通過測點斷面帶來的壓力差異主要原因是:列車和復線隧道形成的幾何空間為非對稱空間,導致列車與隧道之間的空氣流動不對稱,加之列車端部曲線幾何形狀的變化導致流動更為復雜;由于測點“3”和“4”位于靠近車體的隧道一側,所以這2個測點的正壓最大值和負壓最大值均要大于其他位置的測點,測點“4”更是受到地面、隧道壁面和列車效應的多重影響,所以使得其正壓最大值和負壓最大值均要大于其他測點。

2.4隧道內壓力分布特征

由圖10及表1~表3中還可看出:對于距入口端較近的斷面上,不同測點的壓力峰值差異較大,三維效應明顯。而對于距隧道入口較遠的斷面和距出口較近的斷面,各斷面不同測點的壓力峰值基本相同。圖11和圖12分別為列車車頭、車尾通過距隧道入口2D斷面和中央斷面處列車表面和隧道表面的壓力分布特征。由圖可見:列車車頭通過距隧道入口2D的斷面時,隧道壁面的壓力差異較大,靠近車體一側的壓力明顯高于遠離車體一側的壓力,三維效應明顯,車尾通過隧道入口2D的斷面時,隧道壁面的壓力差異較小,而列車通過隧道中央斷面時,中央斷面處的隧道壁面壓力相差也較小。

表1隧道入口2D斷面上不同測點的壓力對比

Table 1 Comparison of the pressure change at different monitors 2D from the entry portal

Pa

表2中央斷面上不同測點的壓力對比

Table 2 Comparison of the pressure change at different monitors at the central cross-section of tunnel

Pa

表3隧道出口2D斷面上不同測點的壓力對比

Table 3 Compression of the pressure change at different monitors 2D from the exit portal

Pa

(a)車頭通過;(b)車尾通過圖11 隧道入口附近的壓力分布Fig.11 Pressure distribution near to the tunnel entrance

(a)車頭通過;(b)車尾通過圖12 隧道中央處的壓力分布Fig.12 Pressure distribution near to the central cross-section of tunnel

2.5列車速度的影響特性

圖13為列車速度對距隧道入口2D斷面上不同測點正壓最大值和負壓最大值的影響特性。由圖可見:隨著車速的增加,各個測點的正壓最大值和負壓最大值也越來越大。在不同列車速度時,由于上述所分析的原因,距隧道入口2D斷面處的正壓最大值和負壓最大值均出現在近車體隧道一側的近地測點“4”上。

(a)正壓最大值;(b)負壓最大值圖13 列車速度對正壓最大值和負壓最大值的影響特性Fig.13 Effect of the train velocity on the maximum positive pressure value and the negative pressure value in tunnel

圖14為列車速度對不同斷面上近車體隧道一側的近地測點“4”的壓力峰值影響。由圖可見:隨著車速的增加,距離隧道入口2D斷面和中央斷面的正壓最大值也隨之增加,而在距隧道出口2D處的斷面上,當速度增加到一定程度,正壓最大值反而下降;負壓最大值隨著車速的增加而增加。當列車速度為300 km/h時,正壓最大值出現在隧道中央斷面處,而在速度為380,450及500 km/h時,正壓最大值均出現在隧道內距離隧道入口2D處,負壓最大值則均出現在隧道的中央斷面處。

(a)正壓最大值;(b)負壓最大值圖14 列車速度對正壓最大值和負壓最大值的影響特性Fig.14 Effect of the train velocity on the maximum positive pressure and the negative pressure in tunnel of different cross-sections

3結論

1)揭示了隧道內壓力波形成的過程。在復線隧道內同一斷面上,隧道中央頂部的壓力與其斷面平均壓力的變化基本相同,壓力值也基本相同。

2)對于距隧道入口較近的斷面上,當列車通過時,不同測點間的壓力差異較大,三維流動效應明顯;而在距離隧道入口較遠的斷面和距離出口較近的斷面上,不同測點間的壓力差異不大,差異主要是由列車頭部和尾部通過帶來的。并且根據每個斷面上測點的壓力值對比發現,由于受到列車周圍流場和隧道壁面、地面效應的多重影響,其正壓最大值和負壓最大值均出現在靠近車體的地面附近。

3)車速越高,距離隧道入口較近斷面上的正壓最大值也越大,而在距離出口較近的斷面上,當列車速度增加到一定程度上,正壓最大值反而減小。

4)隨著車速的增加,正壓最大值多出現在距離隧道入口較近處,而負壓最大值則出現在隧道的中央處。

參考文獻:

[1] Ahmed S R, Gawthorpe R G. Aerodynamics of road and rail vehicles[J]. Vehicle System Dynamics, 1985(14):319-392.

[2] Raghunathan R S, Kim H D, Setoguchi T. Aerodynamics of high -speed railway train [J]. Progress in Aerospace Sciences, 2002(38):469-514.

[3] 王建宇. 列車通過隧道時誘發的空氣動力學問題和高速鐵路隧道設計參數[J]. 世界隧道, 1995(1): 3-13.

WANG Jianyu. Train induced by tunnel aerodynamic problems and high-speed railway tunnel design parameters[J]. World Tunnel, 1995(1): 3-13.

[4] Gawthorpe R G. Pressure effects in railway tunnels[J]. Rail International, 2000(4):10-17.

[5] 田紅旗. 列車空氣動力學[M]. 北京: 中國鐵道出版社, 2007.

TIAN Hongqi. Train Aerodynamics[M]. Beijing: China Railway Press, 2007.

[6] 劉堂紅, 田紅旗. 不同外形列車過隧道實車試驗的比較分析[J]. 中國鐵道科學, 2008, 29(1):51-55.

LIU Tanghong, TIAN Hongqi. Comparison analysis of the fullscale train tests for trains with different shapes passing tunnel[J]. China Railway Science, 2008,29(1):51-55.

[7] 周亞宇. 高速列車通過合武鐵路湖北段隧道空氣動力性能測試[J]. 鐵道建筑,2011(4):73-76.

ZHOU Yayu. Aerodynamic test of high-speed railway tunnels in Hubei section of Hefei Wuhan railway[J]. Railway Engineering, 2011(4):73-76.

[8] 趙有明, 馬偉斌, 程愛君,等. 高速鐵路隧道氣動效應[M]. 北京:中國鐵道出版社, 2012.

ZHAO Youming, MA Weibin, CHENG Aijun, et al. High-speed railway tunnel aerodynamic effect[M]. Beijing: China Railway Press, 2012.

[9] 余南陽. 高速鐵路隧道壓力波數值模擬和模型試驗研究[D]. 成都: 西南交通大學, 2004.

YU Nanyang. Study on numerical simulation of high-speed railway tunnel pressure wave and the model test[D]. Chengdu: Southwest Jiaotong University, 2004.

[10] 王英學, 高波, 朱丹,等. 高速鐵路隧道空氣動力效應控制技術[M].北京:科學出版社,2012.

WANG Yingxue, GAO Bo, ZHU Dan, et al. The control technology on aerodynamics effect of high-speed railway tunnel[M]. Beijing: Science Press,2012.

[11] 梅元貴, 周朝暉, 許建林. 高速鐵路隧道空氣動力學[M]. 北京: 科學出版社, 2009.

MEI Yuangui, ZHOU Chaohui, XU Jianlin. Aerodynamics of high-speed railway tunnel[M]. Beijing: Science Press, 2009.

[12] Ogawa T, Fuji K. Numerical investigation of three-dimensional compressible flows induced by a train moving into a tunnel[J]. Computers and Fluids, 1997, 26(6): 565-585.

[13] Chang-Hoon Shin, Warn-Gyn Park. Numerical study of flow characteristics of the high speed train entering into a tunnel[J]. Mechanics Research Communications,2003(30): 287-296.

[14] 張雷, 楊明智, 張輝,等.高速鐵路隧道洞門對隧道空氣動力效應的影響[J]. 鐵道學報,2013,35(11):92-97.

ZHANG Lei, YANG Mingzhi, ZHANG Hui, et al. Influence of tunnel portals on tunnel aerodynamic effects in operation of high-speed railways[J]. Journal of the China Railway Society, 2013,35(11):92-97.

[15] 王一偉, 楊國偉, 黃晨光. 高速列車通過隧道時氣動阻力特性的CFD仿真分析[J]. 中國鐵道科學,2012,33(增刊):33-38.

WANG Yiwei, YANG Guowei, HUANG Chenguang. CFD simulation analysis on the aerodynamic drag characteristics of high-speed train running through tunnel[J]. China Railway Science, 2012,33(Supplement):33-38.

[16] 趙晶, 李人憲, 劉杰.高速列車隧道內等速會車時氣動作用力的數值模擬[J].鐵道學報,2010,32(4):27-32.

ZHAO Jing, LI Renxian, LIU Jie. Numerical simulation of aerodynamic forces of high-speed trains passing each other at the same speed through a tunnel[J]. Journal of the China Railway Society, 2010,32(4):27-32.

[17] 施成華, 楊偉超, 彭立敏,等.高速鐵路隧道空氣動力效應對水溝蓋板穩定性的影響研究[J]. 鐵道學報,2012,34(1):103-108.

SHI Chenghua, YANG Weichao, PENG Limin, et al. Study on aerodynamic influence on stability of ditch covers in high-speed railway tunnels[J]. Journal of the China Railway Society, 2012,34(1):103-108.

[18] 周丹. 長大隧道、隧道群空氣動力效應算法研究及應用[D]. 長沙: 中南大學, 2007.

ZHOU Dan. Research on the long tunnel and tunnel groups aerodynamic alogorithm and its application[D]. Changsha: Central South University, 2007.

[19] CD-adapco Group. STAR-CD version 4.08 methodology[M]. Computational Dynamics Limited, 2008.

Numerical simulation of pressure wave characteristics generated

by a high-speed train passing through a double-track tunnel

LIU Xixian, GUO Anning, MEI Yuangui, XU Jianlin, LI Shian, ZHOU Chaohui

(College of Mechanical Engineering, Lanzhou Jiaotong University, Lanzhou 730070, China)

Abstract:The CFD method with finite volume method and arbitrary sliding interface moving mesh technology employed was utilized to conduct numerical investigation of the tunnel pressure waves induced by a high-speed train passing through a tunnel. In order to avoid unreasonable numerical physical phenomena, smooth start technology was used. After the validation, three-dimensional unsteady compressible turbulent flow model was adopted to study the formation of the pressure wave in tunnels. The formation process of the pressure wave was revealed, and its characteristics in different cross sections and measuring points were obtained. It shows that the pressure excursion is different in different cross sections near the portals of the tunnel yield by the three-dimensional effect of the initial compression wave. And the maximum positive pressure value also increases with the increase of train velocity. But the maximum positive pressure value decreases near the exit portal of the tunnels when the train velocity increases to a certain range. The maximum positive pressure value always appears near the entry portal of the tunnel, while the maximum negative pressure value occurs near the central section along with the tunnel.

Key words:high-speed railway; double-track tunnel; pressure wave; numerical simulation

中圖分類號:U298.1

文獻標志碼:A

文章編號:1672-7029(2015)01-0020-08

通訊作者:梅元貴(1964-),男,河南滎陽人,教授,博士,從事列車空氣動力學研究;E-mail:meiyuangui@163.com

基金項目:國家自然科學基金資助項目(51065013);國家重點基礎研究發展規劃(“973計劃”)項目(2011CB711101)

*收稿日期:2014-07-24

猜你喜歡
高速列車
基于彩色圖像處理的高速列車快速檢測方法
高速列車新型蜂窩壁板材料的研究現狀
輪軌動力學與安全特性階段研究報告
氣動作用下高速列車響應特性研究
科技資訊(2016年29期)2017-02-28 14:36:58
新型動車組牽引集成單元
物聯網技術(2016年9期)2016-11-09 19:09:44
高速列車復合材料地板振動性能分析
高速列車系統集成試驗的工藝淺談
科技傳播(2016年7期)2016-04-28 00:00:02
動車風道系統的合理化設計
試論焊后退火工藝對高速列車轉向架焊接接頭組織和性能的影響
底部導流板形式對高速列車氣動阻力的影響
主站蜘蛛池模板: 色婷婷亚洲综合五月| 亚洲性日韩精品一区二区| 亚洲黄色片免费看| 亚洲美女一级毛片| 亚洲床戏一区| 国产不卡网| 亚洲精品国产日韩无码AV永久免费网 | 亚洲av无码久久无遮挡| 激情亚洲天堂| 天堂中文在线资源| 国产玖玖玖精品视频| 国产丝袜无码精品| 女人18毛片水真多国产| 亚洲av无码人妻| 亚洲娇小与黑人巨大交| 亚洲第一区精品日韩在线播放| 国产91视频免费| 欧美亚洲激情| 99草精品视频| 久久久久人妻一区精品色奶水 | 色悠久久综合| 欧美伦理一区| 99热最新在线| 欧美日韩中文国产va另类| 国产人碰人摸人爱免费视频| 2022国产91精品久久久久久| 国产精品无码制服丝袜| 国产中文一区二区苍井空| 国产色爱av资源综合区| 美女国内精品自产拍在线播放| 人妻21p大胆| 精品福利视频导航| 精品三级在线| 97久久精品人人做人人爽| 亚洲日韩图片专区第1页| 欧美精品影院| 成人蜜桃网| 九月婷婷亚洲综合在线| 免费精品一区二区h| 国产成人91精品| 在线播放精品一区二区啪视频| 日本午夜三级| 国产亚洲欧美在线中文bt天堂| 午夜高清国产拍精品| 亚洲丝袜中文字幕| 小说 亚洲 无码 精品| 国产区91| 国产性生大片免费观看性欧美| 国产三级国产精品国产普男人 | 国产永久在线视频| 亚洲综合经典在线一区二区| 免费无码AV片在线观看中文| 97超碰精品成人国产| 久久国产亚洲欧美日韩精品| a级免费视频| 国内精自视频品线一二区| 欧洲日本亚洲中文字幕| 亚洲av无码成人专区| 91在线一9|永久视频在线| 91青青草视频| 亚洲丝袜第一页| 亚洲成肉网| 在线国产毛片手机小视频| 美女视频黄又黄又免费高清| 黄色网站不卡无码| 国产福利影院在线观看| a毛片免费在线观看| 国产美女91呻吟求| 91网红精品在线观看| 亚洲高清在线播放| 成年午夜精品久久精品| 国产精品第一区在线观看| 国产在线观看一区精品| 天天做天天爱夜夜爽毛片毛片| 国产精品区视频中文字幕| 色偷偷av男人的天堂不卡| 精品国产福利在线| 一级片一区| 国产正在播放| 久久久噜噜噜久久中文字幕色伊伊| 91精品专区国产盗摄| 99久久国产综合精品2020|