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

基于分離迭代和耦合時變的列車—軌道—橋梁耦合系統高效動力分析混合算法

2018-04-04 07:32:13朱志輝余志武蔡成標
中國鐵道科學 2018年1期
關鍵詞:橋梁

朱志輝, 龔 威,張 磊,余志武,蔡成標

(1.中南大學 土木工程學院,湖南 長沙 410075;2.中南大學 高速鐵路建造技術國家工程實驗室,湖南 長沙 410075;3.西南交通大學 牽引動力國家重點實驗室,四川 成都 610031)

自我國第一條高速鐵路建成通車以來,目前中國高速鐵路總里程已超過2萬km,居世界第一。與普通鐵路相比,評估高速鐵路橋梁動力響應和列車行車安全性與乘坐舒適性指標必不可少[1]。隨著車—橋耦合動力學的發展,多位學者指出,為準確評估上述指標,需要在車—橋耦合動力學模型中考慮軌道結構的柔性變形及參振作用[2]。當考慮軌道結構以后,由于輪軌間較大的接觸剛度,導致車橋耦合動力分析的穩定性和收斂性問題變得更為復雜[3];同時,對于大跨度橋梁而言,大量增加的軌道結構自由度將使耦合系統質量、剛度和阻尼矩陣規模呈幾何倍數增加,導致動力分析計算量大幅增加[4]。因此,需要針對列車—軌道—橋梁耦合系統提出一種高效的計算方法。

很多學者針對車—橋耦合系統動力分析計算效率問題開展了研究工作,主要包括2個方面:模型簡化和數值算法優化。Lou[5]建立了不等長的鋼軌橋梁耦合單元,以減少橋梁結構自由度,并討論了鋼軌單元長度對鋼軌與橋梁動力響應的影響。Yang[6]等假設輪對與鋼軌始終密貼,同時將輪對的自由度壓縮至鋼軌。振型疊加法也常被用于減小橋梁結構自由度數目[7-8],但很難選取合適的振型數量,以同時體現軌道結構的高頻振動和橋梁結構的低頻振動響應[9]。張楠等[10]提出了一種全過程迭代法,將列車和橋梁子系統間每一時間步內的迭代轉換到了整個動力分析過程的迭代。Neves[11]等采用了優化的塊式因式分解算法,提高了車橋耦合動力分析的效率。

以上研究在一定程度上提升了車—橋耦合系統動力響應的計算效率和穩定性,所采用的車—橋耦合系統動力方程求解方法主要分為2類:分離迭代法(Separation and Iteration Method, SIM)和耦合時變法(Coupled and Time-Dependent Method, CTM)[12]。其中SIM將車輛和橋梁視作2個獨立的子系統,2個子系統間通過力的平衡與位移協調條件相互耦合,在每一時間步內需進行迭代求解。整個計算過程中,車輛子系統和橋梁子系統的質量、剛度和阻尼矩陣保持不變。該方法的主要缺點在于耦合系統動力計算穩定性受輪軌接觸點處迭代收斂平衡條件控制,通常需要較小的時間積分步長(Δt≤10-4s),以確保輪軌接觸迭代計算收斂[2]。CTM直接建立車輛和橋梁耦合系統整體時變動力方程,在每一時間步內直接對動力方程進行求解[13-15],避免了每一時間步內的輪軌力收斂迭代計算。但CTM需要在每一時間步根據車輛位置對耦合系統整體矩陣進行更新,尤其當車輛模型較為復雜或存在大量輪軌接觸點時,會增大車—橋耦合模型動力系數矩陣的帶寬,進而增加每一時間步的求解時間[16-17]。

為進一步提高列車—軌道—橋梁相互作用問題的計算效率,本文基于SIM和CTM,提出列車—軌道—橋梁耦合系統高效動力分析混合算法(Hybrid Solution Algorithm, HSA)。

1 列車、軌道、橋梁子系統模型及其動力方程

1.1 列車子系統

一列車由Nv節車輛組成,每節車輛僅考慮在豎向平面內的10個自由度,第j節車的10個自由度為車體沉浮zcj、車體點頭θcj、轉向架沉浮zb1j和zb2j、轉向架點頭θb1j和θb2j以及輪對的沉浮zw(4j-3)—zw4j(所有輪對按照行車方向從前到后按1—4Nv統一編號)。根據多剛體假定和線性懸掛系統假定,每節車的參數可以簡化為:一系懸掛剛度kp和阻尼cp;二系懸掛剛度ks和阻尼cs;車體的質量mc和點頭轉動慣量Jc;構架的質量mb和點頭轉動慣量Jb;輪對的質量mw。輪軌之間采用Hertzian接觸關系模擬。

根據D’Alembert原理,可建立基于平衡位置的列車子系統動力方程

(1)

列車子系統質量矩陣Mvv為

(2)

其中,

Mvj=diag(mcJcmbJbmbJbmw

mwmwmw)

式中:Mvj為第j節車的質量矩陣。

列車子系統剛度矩陣Kvv為

(3)

其中,

式中:Kvj為第j節車的剛度矩陣;Lc為車體定距之半;Lb為構架定距之半。

列車子系統阻尼矩陣可表示為

(4)

式中:Cvj為第j節車的阻尼矩陣,Cvj的形式與剛度矩陣Kvj相同,只需將Kvj的一系懸掛剛度kp與二系懸掛剛度ks替換成一系懸掛阻尼cp與二系懸掛阻尼cs即可。

列車子系統位移向量為

(5)

其中,

Xvj=(zcjθcjzb1jθb1jzb2jθcj

zw(4j-3)zw(4j-2)zw(4j-1)zw4j)

式中:Xvj為第j節車的位移向量。

1.2 軌道子系統

在車—橋耦合振動分析中,軌道結構的彈性變形作用既可以起到過濾輪軌高頻振動的目的,還對輪軌力有顯著影響[18]。以往的研究發現,模態疊加法雖然可以提高結構動力響應的計算效率,但是很難準確體現鋼軌的局部高頻振動特性[9]。因此,本文將采用直接剛度法建立軌道子系統模型。

軌道子系統動力方程為

(6)

1.3 橋梁子系統

通過有限元方法建立精細化橋梁模型,采用直接剛度法建立橋梁子系統動力方程為

(7)

2 列車—軌道—橋梁耦合系統高效動力分析混合算法

HSA將列車—軌道—橋梁耦合系統分解為列車—軌道耦合時變子系統和橋梁子系統2部分,并通過軌道和橋梁在連接節點處的相互作用力將2個子系統耦合起來,實現列車—軌道—橋梁耦合系統的動力仿真分析,如圖1所示。圖中Frb為橋梁對鋼軌的作用力向量,Fbr為鋼軌對橋梁的作用力向量,v為列車車速。

圖1 列車—軌道—橋梁耦合系統示意圖

由于HSA在計算過程中僅需更新列車—軌道時變子系統的動力系數矩陣,而式(7)所示的橋梁子系統動力系數矩陣具有時不變特征,因此,本節主要介紹列車—軌道時變子系統動力方程和HSA求解策略。

2.1 列車—軌道耦合時變子系統動力方程

列車子系統與軌道子系統通過輪軌接觸關系耦合為列車—軌道子系統。根據已有的列車子系統動力方程與軌道子系統動力方程,利用“對號入座法則”[20]將輪軌接觸關系直接寫入列車與軌道的剛度矩陣及其耦合項中,即可得到列車—軌道子系統的動力方程

(8)

其中,

Fv=Fvg+Fvr

Fr=Frv+Frb

在列車—軌道子系統中,輪軌之間采用Hertzian接觸模型[21],由于Hertzian接觸模型僅包含剛度項,因此公式(8)中僅剛度矩陣及荷載向量中存在耦合項。由輪軌接觸引起的車輛剛度附加矩陣Kvh為

(9)

其中,

Kvhj=diag(0000002kh2kh

2kh2kh)

式中:Kvhj為輪軌接觸對第j輛車的影響剛度矩陣;kh為赫茲彈簧剛度。

由輪軌接觸引起的剛度附加矩陣Krr2為

(10)

式中:Nij為(1×Nr)維的向量,表示第i個輪對所在的第j個鋼軌單元的三次Hermit插值形函數[17];Nr為鋼軌自由度數目。

列車運行過程中,列車與鋼軌的剛度矩陣耦合項Kvr為

(11)

其中,

在列車—軌道子系統中,輪軌間的相互作用力為內力,故列車荷載項僅包含列車自重及軌道不平順激勵為

(12)

其中,

Fvgj=(mcg0mbg0mbg0mwg

mwgmwgmwg)

外荷載中車輛對軌道作用僅包含軌道不平順激勵:

(13)

2.2 列車—軌道子系統與橋梁子系統耦合作用

列車—軌道子系統與橋梁子系統間通過梁軌間相互作用力實現耦合,橋梁對鋼軌的作用力Frb為

(14)

鋼軌對橋梁的作用力Fbr為

(15)

2.3 HSA計算流程

圖2 HSA計算流程

與傳統基于CTM的整體列車—軌道—橋梁耦合時變動力學模型相比,HSA將鋼軌與橋梁解耦,大大降低了列車—軌道時變子系統和橋梁子系統動力系數矩陣的帶寬;同時在計算過程中僅需更新列車—軌道時變子系統的動力系數矩陣,橋梁子系統的動力系數矩陣具有時不變特征。相較于CTM,HSA具有更高的計算效率。HSA通過鋼軌與橋梁間作用力的平衡迭代實現列車—軌道時變子系統和橋梁子系統耦合,相對于以往基于輪軌力平衡迭代的SIM而言,由于鋼軌與橋梁間的接觸剛度遠遠小于輪軌間的接觸剛度,因此HSA的迭代穩定性更優。

3 方法驗證

為驗證本文HSA的正確性,以圖3所示的朔黃重載鐵路32 m簡支梁橋為例,分別采用HSA與CTM計算車輛過橋的動力響應,并與現場實測數據對比,驗證HSA的正確性。

圖3 朔黃重載鐵路32 m簡支梁現場圖

3.1 工程概況

朔黃重載鐵路32 m簡支梁橋由4根超低高度的預應力T型梁組成,混凝土強度等級為C60。鋼軌為75 kg·m-1標準軌,軌枕為Ⅲ型混凝土軌枕,道床厚度為22 cm。軌道及橋梁詳細參數見表1。

試驗列車編組為“2節提速機車+116節C80重載列車”,詳細參數見表2。列車編組總長度為1 424.6 m。列車以71.57 km·h-1的速度單線行駛過橋。

3.2 軌道、橋梁有限元模型

采用ANSYS建立軌道、橋梁有限元模型。其中,鋼軌、軌枕、主梁均采用梁單元進行模擬,單元長度均為0.6 m。軌枕與橋面的連接以及鋼軌與軌枕的連接通過彈簧—阻尼器單元模擬。為模擬列車上橋時的初始振動狀態,在橋梁兩側各添加32 m的軌道延長段。由于只考慮列車經過重載線的情況,且重載線和輕載線相互分離,所以忽略輕載線部分。軌道不平順采用美國FRA5級軌道譜。根據以上的建模原則,建立了如圖4所示的軌道、橋梁整體三維有限元模型。

表1 軌道及橋梁參數

表2 提速機車及貨運列車參數表

3.3 結果分析

圖5中給出了HSA與CTM計算所得橋梁跨中豎向位移結果與實測結果。從圖5中可以看出,HSA與CTM計算的橋梁跨中位移響應時程曲線幾乎完全重合,二者誤差最大值僅為0.03%;同時,HSA計算結果與試驗數據誤差為1.79%,表明采用HSA分析車—橋耦合動力響應問題,可以得到滿意的精度。

圖4 橋梁有限元模型(單位:mm)

圖5 橋梁跨中豎向位移時程

4 蒙華鐵路龍門黃河大橋動力分析

4.1 工程概況

為進一步對比HSA和CTM的計算效率,本節以蒙華鐵路龍門黃河大橋為例進行對比分析。該橋為中承式鋼管混凝土拱橋,鋼管混凝土拱跨度202 m,矢跨比1/4,其跨度布置如圖6所示。鐵路等級為國鐵Ⅰ級,重載鐵路標準,采用雙正線,線間距為4.0 m,設計行車速度120 km·h-1。橋梁所有材料均假設為線彈性,鋼材彈性模量取210 GPa,混凝土強度等級為C50;墩身、橋面板采用C40混凝土;吊桿采用整束擠壓式鋼絞線拉索體系,鋼絞線抗拉強度為1 860 MPa。橋面二期恒載為134 kN·m-1(已扣除軌道結構二期恒載16 kN·m-1),包括防水層、保護層、人行道、照明設施等。

圖6 橋梁構造示意圖 (單位:m)

4.2 軌道、橋梁模型

軌道、橋梁有限元模型如圖7所示。其中,橋梁支座、軌枕與橋面的連接以及鋼軌與軌枕的連接通過彈簧—阻尼器單元模擬。根據文獻[23]和設計圖紙資料,彈簧—阻尼器單元參數見表3。橋梁支座布置如圖6所示,共設置16個支座,其中活動支座只設置y和z方向的約束,固定支座設置了x,y和z方向的約束,橋梁拱腳與橋墩墩底均固結。為模擬列車上橋前的初始振動狀態,在橋梁兩側各添加30 m的軌道延長段。

圖7 拱橋有限元模型

參數剛度/MPa阻尼/(kN·s·m-1)x向y向z向x向y向z向鋼軌與軌枕1000306026075軌枕與橋面128128128540540540橋梁支座 100001000010000100010001000

4.3 結果分析

選取與橋長大致相等的列車編組開展耦合系統動力分析,使列車長度足以覆蓋整座橋面。列車編組為2節提速機車+20節C80型貨車的重載列車,列車編組總長度272.6 m。列車以120 km·h-1的設計速度單線行駛過橋,采用美國五級軌道不平順譜,積分步長取1/1 000 s。

圖8和圖9別給出了HSA和CTM計算的列車和橋梁跨中位置的動力響應時程曲線。從圖8和圖9中可以看出,2種方法計算所得第一個輪對豎向輪軌力時程、第一節車車體豎向加速度時程、橋梁跨中豎向位移時程和豎向加速度時程結果吻合度較好,最大誤差均在3%以內。

圖10為跨中位置處鋼軌豎向振動加速度時程。從圖10可以看出:2種方法計算結果吻合程度稍差,最大誤差為8.13%。需要指出的是,經過試算,SIM需采用小于或等于2×10-4s的時間積分步長才可以得到收斂結果,從而表明HSA較SIM具有更好的迭代穩定性。

HSA和CTM的計算耗時見表4。由表4可知,CTM計算耗時為369 min,HSA計算耗時為92 min,HSA較之CTM求解耗時縮短了75%,HSA可有效提升列車—軌道—橋梁耦合系統動力響應計算效率。

圖8 車輛動力響應

圖9 橋梁跨中動力響應

圖10 跨中鋼軌豎向振動加速度時程

表4 2種方法計算耗時對比

5 結 論

(1)HSA計算結果與耦合時變計算結果及現場實測數據吻合度極高,證明了HSA的準確性。

(2)HSA僅需將列車和軌道通過輪軌接觸關系組成耦合時變子系統,從而簡化了列車—軌道—橋梁耦合系統動力方程的理論推導過程,降低了建模難度;同時,列車—軌道子系統耦合時變動力方程不受下部橋梁結構矩陣的影響,因此HSA對于不同橋梁結構具有一定通用性。

(3)HSA在求解過程中,僅上部列車—軌道子系統的質量、剛度、阻尼矩陣是時變的,需要在每一時間步進行更新;下部橋梁子系統各項矩陣保持不變,不需要進行更新。

(4)與CTM相比,HSA每一步更新矩陣所需時間大大減少。以蒙華鐵路龍門黃河大橋算例表明,HSA相較于傳統的CTM,求解耗時降低了75%。

[1]DIMITRAKOPOULOS E G, ZENG Qing. A Three-Dimensional Dynamic Analysis Scheme for the Interaction between Trains and Curved Railway Bridges[J]. Computers & Structures, 2015, 149: 43-60.

[2]吳定俊, 李奇, 陳艾榮. 車橋耦合振動迭代求解數值穩定性問題[J]. 力學季刊, 2007, 28(3): 405-411.

(WU Dingjun, LI Qi, CHEN Airong. Numerical Stability of Iteration Scheme for Solution of Vehicle-Bridge Coupling Vibration [J]. Chinese Quarterly of Mechanics, 2007, 28(3): 405-411.in Chinese)

[3]杜憲亭, 夏禾, 張田. 車橋耦合振動迭代求解穩定性研究[J]. 振動與沖擊, 2012, 31(22): 62-65.

(DU Xianting, XIA He, ZHANG Tian. Numerical Stability of Iterative Scheme in Solving Coupled Vibration of a Train-Bridge System [J]. Journal of Vibration and Shock, 2012, 31(22): 62-65. in Chinese)

[4]LI Yongle, SU Yang, XIA Feilong, et al. Vertical Dynamic Response of the Ballastless Track on Long-Span Plate-Truss Cable-Stayed Bridges [J]. Science China Technological Sciences, 2015, 58(2): 236-247.

[5]LOU Ping, YU Zhiwu, AU F T K. Rail-Bridge Coupling Element of Unequal Lengths for Analysing Train-Track-Bridge Interaction Systems [J]. Applied Mathematical Modelling, 2012, 36(4): 1395-1414.

[6]YANG Yeongbin, WU Yeanseng. A Versatile Element for Analyzing Vehicle-Bridge Interaction Response [J]. Engineering Structures, 2001, 23(5): 452-469.

[7]SALCHER P, ADAM C. Modeling of Dynamic Train-Bridge Interaction in High-Speed Railways [J]. Acta Mechanica, 2015, 226(8): 2473-2495.

[8]JIN Zhibin, PEI Shiling, LI Xiaozhen, et al. Vehicle-Induced Lateral Vibration of Railway Bridges: an Analytical-Solution Approach[J]. Journal of Bridge Engineering, 2015, 21(2):04015038.

[9]ZHANG Nan, XIA He, GUO Weiwei, et al. A Vehicle-Bridge Linear Interaction Model and Its Validation[J]. International Journal of Structural Stability & Dynamics, 2012, 10(2): 335-361.

[10]ZHANG Nan, XIA He. Dynamic Analysis of Coupled Vehicle-Bridge System Based on Inter-System Iteration Method [J]. Computers & Structures, 2013, 115(1): 26-34.

[11]NEVES S G M, MONTENEGRO P A, AZEVEDO A F M, et al. A Direct Method for Analyzing the Nonlinear Vehicle-Structure Interaction [J]. Engineering Structures, 2014, 69(11): 83-89.

[12]CHEN Zhiwei, CHEN Bo. Recent Research and Applications of Numerical Simulation for Dynamic Response of Long-Span Bridges Subjected to Multiple Loads [J]. Scientific World Journal, 2014, 2014(2): 763810.

[13]WANG Wei, YAN Wangchen, DENG Lu, et al. Dynamic Analysis of a Cable-Stayed Concrete-Filled Steel Tube Arch Bridge under Vehicle Loading[J]. Journal of Bridge Engineering, 2015, 20(5): 1-20.

[14]DENG Lu, CAI C S. Development of Dynamic Impact Factor for Performance Evaluation of Existing Multi-Girder Concrete Bridges[J]. Engineering Structures, 2010, 32(1): 21-31.

[15]GUO W H, XU Y L. Fully Computerized Approach to Study Cable-Stayed Bridge-Vehicle Interaction[J]. Journal of Sound & Vibration, 2001, 248(4): 745-761.

[16]GU G. Resonance in Long-Span Railway Bridges Carrying TGV Trains[J]. Computers & Structures, 2015, 152: 185-199.

[17]LOU Ping. Finite Element Analysis for Train-Track-Bridge Interaction System[J]. Archive of Applied Mechanics, 2007, 77(10): 707-728.

[18]張攀, 周昌盛, 王平. 軌下墊板剛度的時變特性及其影響研究[J]. 鐵道標準設計, 2015(9): 49-52.

(ZHANG Pan, ZHOU Changsheng, WANG Ping. Study on Time Variant Characteristics and Effects of Rail Pad Stiffness [J]. Railway Standard Design, 2015(9): 49-52. in Chinese)

[19]朱志輝, 龔威, 王力東, 等. 列車—軌道—橋梁耦合系統動力方程求解方法對計算精度和效率的影響[J]. 中國鐵道科學, 2016, 37(5): 17-26.

(ZHU Zhihui, GONG Wei, WANG Lidong, et al. Influence of Solution Method for Dynamics Equation of Train-Track-Bridge Coupled System on Calculation Precision and Efficiency [J]. China Railway Science, 2016, 37(5): 17-26. in Chinese)

[20]曾慶元, 楊平. 形成矩陣的“對號入座”法則與桁梁空間分析的桁段有限元法[J]. 鐵道學報, 1986,8(2): 50-61.

(ZENG Qingyuan, YANG Ping. The “Set-in-Right-Position” Rule for Forming Structural Matrices and the Finite Truss-Element Method for Space Analysis of Truss Bridges [J]. Journal of the China Railway Society, 1986, 8(2): 50-61. in Chinese)

[21]朱志輝, 王力東, 龔威, 等. 基于改進迭代模型的車橋耦合系統豎向隨機振動研究[J]. 湖南大學學報:自然科學版, 2016, 43(11): 120-130.

(ZHU Zhihui, WANG Lidong, GONG Wei, et al. Study on Vertical Random Vibration of Train-Bridge Coupled System Based on Improved Iteration Model [J]. Journal of Hunan University:Natural Sciences, 2016, 43(11): 120-130. in Chinese)

[22]NGUYEN D V, KIM K D, WARNITCHAI P. Simulation Procedure for Vehicle-Substructure Dynamic Interactions and Wheel Movements Using Linearized Wheel-Rail Interfaces[J]. Finite Elements in Analysis and Design, 2009, 45(5): 341-356.

[23]劉瑋, 曲村. 高速鐵路橋上有砟軌道軌枕選型方案研究[J]. 高速鐵路技術, 2011, 2(3): 38-42.

(LIU Wei, QU Cun. Study on Selection of Sleepers for Ballasted Track on Bridges of High-Speed Railway [J]. High Speed Railway Technology, 2011, 2(3): 38-42. in Chinese)

猜你喜歡
橋梁
一種橋梁伸縮縫防滲水裝置
工程與建設(2019年4期)2019-10-10 01:45:56
手拉手 共搭愛的橋梁
句子也需要橋梁
加固技術創新,為橋梁健康保駕護航
中國公路(2017年11期)2017-07-31 17:56:30
無人機在橋梁檢測中的應用
中國公路(2017年10期)2017-07-21 14:02:37
高性能砼在橋梁中的應用
現代鋼橋制造對橋梁鋼的更高要求
焊接(2016年8期)2016-02-27 13:05:15
城鄉建設一體化要注重橋梁的建筑設計
南昌54座橋梁進行兩個月的夏季體檢
橋梁伸縮縫損壞因素與加固
主站蜘蛛池模板: 九色综合伊人久久富二代| 91精品国产91久无码网站| 无套av在线| 久久精品视频亚洲| 26uuu国产精品视频| 无码电影在线观看| 91无码视频在线观看| 71pao成人国产永久免费视频 | 国产福利微拍精品一区二区| 亚洲αv毛片| 国产av一码二码三码无码 | 亚洲欧美成人综合| 免费国产高清精品一区在线| a毛片免费在线观看| 亚洲欧美不卡视频| 久久精品无码中文字幕| 伊人网址在线| 国产精品片在线观看手机版| 麻豆精品久久久久久久99蜜桃| 免费观看男人免费桶女人视频| 日韩天堂在线观看| 欧美色99| 91视频精品| 欧洲在线免费视频| 国产一区成人| 在线免费亚洲无码视频| 美女被狂躁www在线观看| 一级毛片a女人刺激视频免费| 日韩国产欧美精品在线| 97国产成人无码精品久久久| 精品国产一区91在线| 国产一区二区在线视频观看| 成人欧美日韩| 国产成人麻豆精品| 亚洲无码熟妇人妻AV在线| 日韩国产一区二区三区无码| 高清无码不卡视频| 日韩免费毛片| 激情六月丁香婷婷四房播| 国产成人综合网| 毛片网站观看| 欧美日韩激情在线| 91 九色视频丝袜| 人妻一区二区三区无码精品一区| 久久亚洲黄色视频| 久久一本日韩精品中文字幕屁孩| 无码高潮喷水在线观看| 国产亚洲男人的天堂在线观看| 亚洲欧美一区二区三区图片| av一区二区无码在线| 成人一区专区在线观看| 无码精品国产VA在线观看DVD| 免费毛片在线| 国产亚洲精品91| 久久精品娱乐亚洲领先| 丝袜高跟美脚国产1区| 九色综合视频网| 国产日韩久久久久无码精品| 免费毛片视频| 国产一区二区影院| 色婷婷成人网| 国产精品.com| 亚洲精品中文字幕午夜| 欧美不卡二区| 成人午夜精品一级毛片| 亚洲精品无码在线播放网站| 亚洲AV无码乱码在线观看代蜜桃| 亚洲国产成人无码AV在线影院L| av在线无码浏览| 午夜欧美在线| 91精品aⅴ无码中文字字幕蜜桃| 无码网站免费观看| 中国国产A一级毛片| 亚洲精品第一页不卡| 日韩少妇激情一区二区| 国产精品嫩草影院视频| 3344在线观看无码| 大香伊人久久| 香蕉eeww99国产在线观看| 国产乱人伦AV在线A| 亚洲视频三级| 香蕉伊思人视频|