亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        一種高精度軌道單元及其在車輛-軌道耦合系統(tǒng)的應(yīng)用

        2024-01-06 01:52:22劉林芽崔巍濤秦佳良左志遠(yuǎn)
        鐵道學(xué)報(bào) 2023年12期
        關(guān)鍵詞:膠墊輪軌扣件

        劉林芽,崔巍濤,秦佳良,左志遠(yuǎn)

        (華東交通大學(xué) 軌道交通基礎(chǔ)設(shè)施性能監(jiān)測(cè)與保障國(guó)家重點(diǎn)實(shí)驗(yàn)室,江西 南昌 330013)

        隨著列車運(yùn)行速度的提升、行車密度的提高以及新型車輛和軌道結(jié)構(gòu)形式大量投入工程應(yīng)用,導(dǎo)致輪軌相互作用更加復(fù)雜,甚至?xí)绊懥熊嚨陌踩\(yùn)行[1-3]。進(jìn)行車輛-軌道系統(tǒng)動(dòng)力學(xué)分析是研究復(fù)雜輪軌關(guān)系和相互作用機(jī)制的基礎(chǔ)[4],是優(yōu)化車輛、軌道結(jié)構(gòu)設(shè)計(jì)必不可少的內(nèi)容。

        針對(duì)軌道結(jié)構(gòu)的建模方面,國(guó)內(nèi)外學(xué)者開(kāi)展了大量的研究[5-7]。Clark等[8]建立考慮鋼軌彈性點(diǎn)支承和軌枕的振動(dòng)的有砟軌道結(jié)構(gòu)模型,分析了鋼軌波磨對(duì)輪對(duì)與軌道結(jié)構(gòu)振動(dòng)的影響。Nielsen等[9]建立包含轉(zhuǎn)向架、輪對(duì)、下部軌道基礎(chǔ)為一體的有限元模型,采用復(fù)模態(tài)疊加法分析了車輪扁平、鋼軌磨耗等因素對(duì)軌道結(jié)構(gòu)振動(dòng)影響的影響。雷曉燕[10]基于有限元理論建立了車輪-有砟軌道動(dòng)力計(jì)算模型,分析高速列車對(duì)不同道砟厚度的動(dòng)力響應(yīng)。隨著高速鐵路的發(fā)展,無(wú)砟軌道由于整體性好和養(yǎng)護(hù)維修工作量少等優(yōu)點(diǎn)得到了廣泛應(yīng)用。Guigou-Carter等[11]提出一種二維的板式軌道預(yù)測(cè)模型,扣件的阻尼建模為滯后阻尼且采用復(fù)剛度進(jìn)行考慮,分析了扣件的動(dòng)剛度變化對(duì)下部軌道結(jié)構(gòu)減振性能的影響規(guī)律。Khajehdezfuly等[12]建立二維的車輛-板式無(wú)砟軌道有限元分析模型,鋼軌采用Timoshenko梁模擬,扣件采用點(diǎn)支承的彈簧阻尼單元模擬,研究了扣件剛度變化對(duì)輪軌相互作用的影響規(guī)律。翟婉明等[13]基于模態(tài)疊加法建立一種車輛-板式軌道相互作用模型,鋼軌視為連續(xù)彈性點(diǎn)支承基礎(chǔ)上的歐拉梁,研究了砂漿層阻尼或剛度改變對(duì)系統(tǒng)振動(dòng)的影響規(guī)律。蔡成標(biāo)[14]建立列車-軌道-基礎(chǔ)耦合分析模型,將鋼軌視為彈性點(diǎn)支承的歐拉梁,軌道板和底座板視為彈性薄板,研究了板式軌道及過(guò)渡段的動(dòng)力學(xué)性能。文獻(xiàn)[15-16]基于彈性系統(tǒng)動(dòng)力學(xué)總勢(shì)能不變值原理,提出一種橫向有限條與無(wú)砟軌道板段單元,比較了鋼軌與軌道板的靜、動(dòng)態(tài)位移變化。Lei等[17]提出一種新型板式軌道單元,建立車輛-板式軌道-路基耦合動(dòng)力學(xué)模型,評(píng)估了列車速度和軌道剛度對(duì)軌道振動(dòng)的影響。上述研究對(duì)于軌道結(jié)構(gòu)中扣件的模擬多采用點(diǎn)支承的彈簧阻尼單元模擬而忽略扣件墊板的支承長(zhǎng)度。

        另有一些學(xué)者采用仿真軟件建立精細(xì)化軌道模型進(jìn)行研究。Ganesh等[18]建立有砟軌道結(jié)構(gòu)有限元模型,扣件墊板采用solid92單元模擬,分析諧響應(yīng)激勵(lì)下扣件、道砟和路基參數(shù)變化對(duì)軌道結(jié)構(gòu)的影響。Oregui等[19]提出一種軌道有限元模型,鋼軌和軌枕采用實(shí)體單元,扣件采用多個(gè)彈簧-阻尼單元組成的面區(qū)域模擬,研究了錘擊試驗(yàn)下軌道結(jié)構(gòu)的振動(dòng)特性以及扣件膠墊建模形式對(duì)軌道動(dòng)力學(xué)的影響。Zhao等[20]采用有限元法建立科隆蛋扣件系統(tǒng)的軌道真實(shí)有限元頻域分析模型,利用該模型進(jìn)行了數(shù)值模擬和輪軌振動(dòng)特性分析。上述模型建立往往需借助仿真軟件,單元數(shù)量多,計(jì)算效率不高。

        本文建立一種更精確的軌道單元模型,模型中考慮扣件膠墊的支承長(zhǎng)度和鋼軌的剪切變形,可以得到更準(zhǔn)確的高頻振動(dòng)響應(yīng)。為便于計(jì)算程序的設(shè)計(jì),減小總剛度矩陣的帶寬,取兩扣件之間的一跨鋼軌、軌道板、底座板作為一個(gè)單元,鋼軌視為扣件分布支承的Timoshenko梁,軌道板和底座板視為連續(xù)支承的Euler梁,軌道結(jié)構(gòu)采用黏彈性彈簧-阻尼單元連接。車輛和軌道通過(guò)非線性輪軌接觸關(guān)系耦合在一起,采用顯-隱式積分法求解耦合系統(tǒng)的動(dòng)力響應(yīng),并驗(yàn)證算法的正確性,然后比較扣件集中支承和分布支承下車輛、軌道結(jié)構(gòu)的動(dòng)力響應(yīng)及輪軌相互作用,最后研究不同扣件膠墊支承長(zhǎng)度對(duì)車輛、軌道結(jié)構(gòu)振動(dòng)的影響。

        1 車輛-軌道耦合系統(tǒng)動(dòng)力學(xué)模型

        根據(jù)車輛-軌道耦合動(dòng)力學(xué)理論,建立高速列車-板式軌道耦合系統(tǒng)模型,車輛為考慮一系和二系懸掛的整車模型,軌道為三層板式無(wú)砟軌道模型,其中采用分布式彈簧阻尼元件模擬扣件墊板的支承長(zhǎng)度,建立更加精確的軌道單元模型,兩個(gè)子系統(tǒng)通過(guò)非線性接觸關(guān)系耦合,見(jiàn)圖1。

        圖1 車輛-軌道耦合系統(tǒng)模型

        Mc和Jc分別為車體的質(zhì)量和轉(zhuǎn)動(dòng)慣量;Mt和Jt分別為轉(zhuǎn)向架的質(zhì)量和轉(zhuǎn)動(dòng)慣量;Mw為輪對(duì)的質(zhì)量;Ks1和Cs1分別為一系懸掛的剛度和阻尼;Ks2和Cs2分別為二系懸掛的剛度和阻尼;lc為車體上兩轉(zhuǎn)向架中心之間的長(zhǎng)度;lt為轉(zhuǎn)向架兩軸之間的距離。Z0i(t)、Pwi(t)(i=1,2,3,4)分別為第i個(gè)車輪接觸處的不平順?lè)岛洼嗆壌瓜蛄?ρr、Er、Ir分別為單位長(zhǎng)度鋼軌的質(zhì)量、鋼軌的彈性模量和繞水平軸的截面慣性矩;ρs、Es、Is分別為單位長(zhǎng)度軌道板的質(zhì)量、軌道板的彈性模量和繞水平軸的截面慣性矩;ρh、Eh、Ih分別為單位長(zhǎng)度底座板的質(zhì)量、底座板的彈性模量和繞水平軸的截面慣性矩。kp、cp分別為扣件的剛度和阻尼;km、cm分別為軌道板下瀝青水泥砂漿層的剛度和阻尼;kl、cl分別為路基層的剛度和阻尼。

        1.1 高精度的軌道單元模型

        軌道子系統(tǒng)包括鋼軌、扣件、預(yù)制軌道板、CA砂漿、底座板和路基,采用有限元法建立三層板式軌道單元模型[21]。當(dāng)扣件考慮為集中支承時(shí),板式軌道單元的長(zhǎng)度l取鋼軌兩個(gè)相鄰離散支承點(diǎn)之間的距離,軌道板和底座板的長(zhǎng)度與鋼軌一致,將鋼軌、軌道板和底座板視為組合單元,見(jiàn)圖2(a);當(dāng)扣件考慮為連續(xù)支承時(shí),仍取該長(zhǎng)度作為板式軌道單元的長(zhǎng)度,此時(shí)扣件考慮為連續(xù)支承的彈簧阻尼單元,長(zhǎng)度lr為二分之一的扣件膠墊寬度,見(jiàn)圖2(b)。

        圖2 三層板式軌道單元模型

        兩種板式軌道單元模型為12個(gè)自由度,v1、v4分別為鋼軌的豎向位移;v2、v5分別為軌道板的豎向位移;v3、v6分別為底座板的豎向位移;θ1、θ4分別為鋼軌的轉(zhuǎn)角;θ2、θ5分別為軌道板的轉(zhuǎn)角;θ3、θ6分別為底座板的轉(zhuǎn)角。只需形成單元的質(zhì)量、剛度和阻尼矩陣,便可方便地按有限元集成規(guī)律組集,得到任意長(zhǎng)度的軌道系統(tǒng)的總質(zhì)量、總剛度和總阻尼矩陣。單元內(nèi)任何一點(diǎn)的位移可以通過(guò)插值函數(shù)N和單元節(jié)點(diǎn)位移ae表示,即

        X=Nae

        (1)

        式中:

        ae=

        1.1.1 板式軌道單元質(zhì)量矩陣

        板式無(wú)砟軌道單元的質(zhì)量矩陣為

        (2)

        (3)

        式中:

        (4)

        式中:As為軌道板的截面面積。

        (5)

        式中:Ah為底座板的截面面積。

        1.1.2 板式軌道單元?jiǎng)偠染仃?/p>

        板式軌道單元的剛度矩陣為

        (6)

        (7)

        (8)

        (9)

        單元中任意一點(diǎn)鋼軌與軌道板之間的相對(duì)位移vrs為

        vrs=vr-vs=N1v1+N2θ1+N3v4+N4θ4-N5v2+

        N6θ2+N7v5+N8θ5=

        (10)

        故扣件膠墊層的彈性勢(shì)能為

        (11)

        扣件分布支承引起的剛度為

        (12)

        (13)

        (14)

        1.1.3 板式軌道單元阻尼矩陣

        (15)

        最后,按有限元集成規(guī)則組集單元質(zhì)量、剛度和阻尼矩陣(式(2)、式(6)和式(15)),可得到三層板式軌道結(jié)構(gòu)的總質(zhì)量、總剛度和總阻尼矩陣分別為

        (16)

        運(yùn)用Lagrange方程,可以得到板式軌道結(jié)構(gòu)的有限元方程,寫(xiě)成統(tǒng)一的形式,即

        (17)

        1.2 車輛單元模型

        車輛子系統(tǒng)考慮一個(gè)車體、兩個(gè)轉(zhuǎn)向架的點(diǎn)頭和沉浮及四個(gè)輪對(duì)的沉浮共10個(gè)自由度。定義車輛單元的位移向量為

        (18)

        式中:vc和θc分別為車體的豎向位移和角位移;vti和θti(i=1,2)分別為前、后轉(zhuǎn)向架的豎向位移和角位移;vwi(i=1,2,3,4)為第i個(gè)輪對(duì)的豎向位移。

        運(yùn)用Lagrange方程,可得到車輛子系統(tǒng)的動(dòng)力學(xué)方程為

        (19)

        1.3 輪軌接觸關(guān)系

        采用Hertz非線性彈性接觸理論,輪軌力為

        (20)

        式中:G為輪軌接觸常數(shù),m/N2/3;δZj(t)(j=1~4)為輪軌間的彈性壓縮量,m。

        輪軌間的彈性壓縮量為

        δZj(t)=vwj(t)-Zrj(t)-Z0j(t)

        (21)

        式中:vwj(t)、Zrj(t)分別為輪對(duì)和鋼軌接觸位置的垂向位移;Z0j(t)為軌道不平順?lè)怠?/p>

        2 數(shù)值求解及模型驗(yàn)證

        2.1 顯-隱式數(shù)值積分法

        車輛和軌道兩個(gè)子系統(tǒng)通過(guò)非線性輪軌接觸力實(shí)現(xiàn)耦合,故需要采用數(shù)值積分進(jìn)行求解,常用的數(shù)值積分方法有隱式法和顯式法。顯式法中,新型快速顯式積分法[23]在質(zhì)量矩陣為對(duì)角陣或經(jīng)對(duì)角化,不需要聯(lián)立求解高階線性代數(shù)方程組,求解效率高且更經(jīng)濟(jì)。隱式法中,當(dāng)考慮高頻時(shí)要求極小的時(shí)間積分步長(zhǎng),隱式法每一步長(zhǎng)都需要求解一次動(dòng)力學(xué)方程,且對(duì)于非線性問(wèn)題,需要采用交叉迭代控制精度范圍內(nèi)系統(tǒng)時(shí)程響應(yīng)始終是收斂的,計(jì)算工作量大。

        本文車輛系統(tǒng)因采用集中質(zhì)量法建模,質(zhì)量矩陣為對(duì)角陣,采用新型顯式積分法求解;軌道模型采用有限元法建模,質(zhì)量矩陣為協(xié)調(diào)質(zhì)量矩陣非對(duì)角陣,故采用隱式Newmark數(shù)值積分法求解,即采用顯式-隱式數(shù)值積分法求解車輛-軌道非線性耦合系統(tǒng)動(dòng)力分析模型,避免了交叉迭代,提高了計(jì)算效率,同時(shí)具有較好的結(jié)果。采用Matlab編制相應(yīng)的程序,主要計(jì)算步驟如下:

        Step3根據(jù)軌道結(jié)構(gòu)位移Xl(t)可得到輪軌接觸處的鋼軌位移Zrj(t),由式(20)計(jì)算輪軌力Pwj(t)。

        需注意的是,在起步時(shí),先假設(shè)軌道結(jié)構(gòu)初始位移Xl(0)(一般取Xl(0)=0)。同時(shí),顯式積分法令φ=ψ=0,具有積分“自開(kāi)始”的特性。

        2.2 模型驗(yàn)證

        采用文獻(xiàn)[24]中相同的條件對(duì)模型和求解方法進(jìn)行驗(yàn)證。車輛模型為日本新干線300系列高速列車,軌道為三層板式無(wú)砟軌道,軌道不平順選用鋼軌接頭焊縫不平順。鋼軌焊縫不平順滿足

        (22)

        式中:a1和L1分別為長(zhǎng)波的波深和波長(zhǎng);a2和L2分別為短波的波深和波長(zhǎng)。

        當(dāng)車輛以250 km/h的速度通過(guò)焊接區(qū)時(shí),長(zhǎng)波的波深和波長(zhǎng)分別為0.3 mm和1 m,短波的波深和波長(zhǎng)為0.2 mm和0.1 m。采用顯-隱式積分法求解板式軌道的動(dòng)態(tài)響應(yīng),積分步長(zhǎng)為0.000 1 s。顯隱式積分法計(jì)算結(jié)果見(jiàn)圖3。并與文獻(xiàn)[24]中的結(jié)果進(jìn)行對(duì)比,見(jiàn)圖4。由圖4可知,本文所提出的模型計(jì)算結(jié)果與文獻(xiàn)的結(jié)果吻合較好,表明本文建立模型的正確性。文獻(xiàn)[24]中積分步長(zhǎng)為0.000 02 s,本文采用較小的積分步長(zhǎng)仍得到較好的計(jì)算結(jié)果,提高了計(jì)算效率。

        圖3 顯隱式積分法計(jì)算結(jié)果

        圖4 文獻(xiàn)[24]的計(jì)算結(jié)果

        3 扣件分布支承式軌道模型下車輛-軌道耦合系統(tǒng)動(dòng)力響應(yīng)分析

        3.1 計(jì)算參數(shù)

        在仿真分析中,車輛選用中國(guó)CRH3高速客車,運(yùn)行速度為300 km/h,軌道采用CRTS Ⅱ型無(wú)砟軌道結(jié)構(gòu),具體參數(shù)見(jiàn)表1。

        表1 車輛與軌道具體參數(shù)

        3.2 扣件支承形式對(duì)耦合系統(tǒng)動(dòng)力響應(yīng)影響

        軌道不平順選擇我國(guó)高速鐵路無(wú)砟軌道不平順譜,波長(zhǎng)范圍為2~200 m,同時(shí)將其延伸至波長(zhǎng)0.1 m作為短波不平順激勵(lì)[25]。應(yīng)用上述的模型和計(jì)算參數(shù),研究采用不同扣件支承形式下車輛-軌道耦合系統(tǒng)的動(dòng)力響應(yīng),其中FCS為扣件集中支承,FDS為扣件分布支承,見(jiàn)圖5。圖5(e)、圖5(h)中鋼軌的動(dòng)力響應(yīng)取在軌道結(jié)構(gòu)跨中位置。

        如圖5(a)~圖5(e)所示,在時(shí)域內(nèi),與FCS模型相比,FDS模型的輪軌力、輪對(duì)加速度和鋼軌加速度峰值均減小,波動(dòng)范圍變窄,FCS模型的輪軌力、輪對(duì)和鋼軌加速度最大值分別為98.291 kN、21.419 m/s2、247.184 m/s2,FDS模型的輪軌力、輪對(duì)和鋼軌加速度最大值分別為97.705 kN、20.681 m/s2、239.874 m/s2,差異約為0.6%、4.231%和3.047%。在頻域內(nèi),兩種模型的輪軌力在70~140 Hz以及1 000 Hz附近存在較小差異,FDS模型的振幅通常小于DCS模型振幅。輪對(duì)和鋼軌加速度采用1/3倍頻加速度振級(jí)可以明顯看出,兩種模型在125、1 000 Hz頻段存在區(qū)別,FCS模型大于FDS模型的計(jì)算結(jié)果。第一主頻在63 Hz,由車輛與軌道共同彈性變形的耦合共振引起。第二主頻為1 000 Hz附近,由鋼軌的垂向“pinned-pinned”振動(dòng)引起??紤]膠墊支承長(zhǎng)度的FDS模型在1 000 Hz頻段加速度振級(jí)更小,與文獻(xiàn)[26]規(guī)律一致,更精確地反映了兩個(gè)模型之間的差異。從圖5(g)~圖5(h)可以看出,對(duì)于輪對(duì)位移和跨中鋼軌位移,與FCS模型相比,FDS模型的輪對(duì)和鋼軌的位移均較小,但差異均小于0.5%。

        3.3 扣件膠墊支承長(zhǎng)度對(duì)車輛-軌道耦合系統(tǒng)的影響分析

        在FDS模型的基礎(chǔ)上,進(jìn)一步研究扣件系統(tǒng)中膠墊支承長(zhǎng)度對(duì)車輛-軌道耦合系統(tǒng)垂向動(dòng)態(tài)響應(yīng)的影響,橡膠墊板的支承長(zhǎng)度設(shè)置為50、150、250 mm,計(jì)算結(jié)果見(jiàn)圖6。為了更直觀地反映膠墊不同支承長(zhǎng)度的區(qū)別,不考慮外部激勵(lì)作用,并假設(shè)鋼軌表面平順。另外,圖6(a)、圖6(b)響應(yīng)參照系取在車輛上,圖6(c)~圖6(h)鋼軌、軌道板和底座板的響應(yīng)參照系取在地面上。

        圖6 扣件膠墊支承長(zhǎng)度對(duì)軌道結(jié)構(gòu)的影響

        如圖6(a)、圖6(b)所示,由于彈性膠墊的周期性支承,輪軌力和輪軌接觸處的輪對(duì)位移呈現(xiàn)周期性變化,變化幅值隨著膠墊支承長(zhǎng)度的增加而減小。當(dāng)扣件膠墊長(zhǎng)度分別取50、150、250 mm時(shí),最大輪軌力分別為71.278、71.205、71.167 kN;最大輪對(duì)位移分別為1.005 38、1.004 89、1.004 51 mm。隨著膠墊支承長(zhǎng)度從50 mm增加到250 mm,最大輪軌力和輪對(duì)位移減小了0.15%和0.08%。圖6(c)~圖6(e)反映了隨膠墊支承長(zhǎng)度變化,鋼軌、軌道板和底座板加速度的變化,隨著膠墊長(zhǎng)度的增加,鋼軌、軌道板和底座板的加速度均減小,其中膠墊支承長(zhǎng)度的變化對(duì)鋼軌影響最大。當(dāng)扣件膠墊長(zhǎng)度取50 mm,鋼軌、軌道板和底座板加速度最大值分別為15.151、0.736、0.598 m/s2;扣件膠墊長(zhǎng)度取250 mm,鋼軌、軌道板和底座板加速度最大值分別為14.627、0.718、0.577 m/s2,分別減小3.582%、2.506%和3.639%。圖6(f)~圖6(h)反映了膠墊支承長(zhǎng)度變化鋼軌、軌道板和底座板位移的變化,隨著膠墊長(zhǎng)度的不斷增加,鋼軌、軌道板和底座板的位移也均略有減小,但影響較小。隨著支撐長(zhǎng)度從50 mm增加到250 mm,鋼軌、軌道板和底座板的最大位移分別減小了0.6%、0.2%和0.08%。

        4 結(jié)論

        本文基于車輛-軌道耦合動(dòng)力學(xué)理論,采用有限元法建立了更精確的考慮扣件膠墊的支承長(zhǎng)度和鋼軌的剪切變形的軌道單元模型,將其應(yīng)用于車輛-軌道耦合系統(tǒng)動(dòng)力學(xué),采用顯-隱式積分法求解耦合系統(tǒng)動(dòng)力學(xué)響應(yīng)。研究了扣件模型支承形式和長(zhǎng)度對(duì)車-軌耦合動(dòng)力學(xué)的影響,得到以下結(jié)論:

        1)該軌道單元模型考慮了鋼軌的剪切變形及扣件的分布支承,具有精度高、編程簡(jiǎn)捷清晰的優(yōu)點(diǎn)。

        2)扣件考慮集中支承時(shí)會(huì)高估輪軌相互作用以及輪對(duì)和鋼軌的振動(dòng),采用扣件分布支承式軌道單元模型可以更精確地預(yù)測(cè)“pined-pined”頻率范圍附近的輪軌相互作用。

        3)扣件膠墊支承長(zhǎng)度的增加可以降低輪軌相互作用,實(shí)際工程中可以適當(dāng)增加扣件膠墊的支承長(zhǎng)度來(lái)減小軌道結(jié)構(gòu)的振動(dòng)。

        4)采用本文建立的車輛-軌道耦合模型及顯-隱式積分算法可以在降低計(jì)算成本的情況下給出更為準(zhǔn)確的軌道結(jié)構(gòu)振動(dòng)響應(yīng)。

        猜你喜歡
        膠墊輪軌扣件
        科隆蛋扣件力學(xué)特性分析
        中低速磁浮道岔與輪軌道岔的差異
        一種改進(jìn)的LBP特征實(shí)現(xiàn)鐵路扣件識(shí)別
        一種基于工業(yè)機(jī)器人的膠墊自動(dòng)粘貼系統(tǒng)的研究與開(kāi)發(fā)
        換膠墊250
        非均布荷載下扣件式腳手架的隨機(jī)缺陷穩(wěn)定分析
        中低速磁浮與輪軌交通信號(hào)系統(tǒng)的差異
        WJ-8型小阻力扣件軌下橡膠墊板滑出動(dòng)力學(xué)研究
        WJ -8型小阻力扣件軌下膠墊滑出縱向阻力試驗(yàn)研究
        非線性穩(wěn)態(tài)曲線通過(guò)時(shí)輪軌滾動(dòng)接觸的數(shù)值求解方法
        小黄片免费在线播放观看| 粗大的内捧猛烈进出在线视频| 国产最新AV在线播放不卡| 亚洲av色香蕉一区二区三区蜜桃| 大香蕉国产av一区二区三区| 18禁裸男晨勃露j毛网站| 女人夜夜春高潮爽a∨片传媒| 亚洲www视频| 亚洲av男人的天堂在线| 高h喷水荡肉爽文np肉色学校| 高中生粉嫩无套第一次| 男人天堂av在线成人av| av一区二区三区观看| 久久天天躁夜夜躁狠狠| 国产精品va在线播放我和闺蜜 | 午夜无码熟熟妇丰满人妻| 亚洲一区二区自偷自拍另类| 色欲人妻aaaaaaa无码| 成人免费毛片内射美女-百度| 人妻少妇精品无码系列| 成人做爰黄片视频蘑菇视频| 色偷偷av一区二区三区| 亚洲一区爱区精品无码| 亚洲av午夜福利精品一区二区 | 护士的小嫩嫩好紧好爽| 亚洲av久久无码精品九九| 久久aⅴ无码av高潮AV喷| 亚洲av老熟女一区二区三区| 香港台湾经典三级a视频| 91福利国产在线观看一区二区| 午夜av福利亚洲写真集| 51国产偷自视频区视频| 中文国产日韩欧美二视频| 亚洲欧洲一区二区三区波多野| 日本一区二区在线播放视频| а√天堂资源官网在线资源 | 日韩精人妻无码一区二区三区| 久久亚洲春色中文字幕久久| 亚洲精品久久区二区三区蜜桃臀| 国产精品久免费的黄网站| 精品在线视频免费在线观看视频|