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

基于2.5維有限元法和虛擬激勵法的地鐵交通場地隨機(jī)振動分析

2020-08-06 02:32:38朱志輝劉禹兵王力東王東旭
中國鐵道科學(xué) 2020年4期
關(guān)鍵詞:振動模型

朱志輝,劉禹兵,王力東,王 凡,王東旭

(1.中南大學(xué) 土木工程學(xué)院,湖南 長沙 410075;2.中南大學(xué) 高速鐵路建造技術(shù)國家工程實(shí)驗(yàn)室,湖南 長沙 410075;3.長沙理工大學(xué) 土木工程學(xué)院,湖南 長沙 410114;4.武九鐵路客運(yùn)專線湖北有限責(zé)任公司,山東 武城 430000)

隨著我國軌道交通的不斷發(fā)展與列車行駛速度的逐步提升,鐵路交通的環(huán)境振動問題對人們的日常生活和工作產(chǎn)生重要影響[1],因此也越來越受人們重視。為了系統(tǒng)地了解軌道交通產(chǎn)生振動的規(guī)律和機(jī)制,就需要建立起相應(yīng)的振動模型來分析研究,以便對軌道沿線做出正確的環(huán)境振動分析與評估。

近些年來,研究人員針對軌道交通引起的環(huán)境振動問題進(jìn)行了大量的工作[2]。Sheng 等[3]采用半解析法研究了移動簡諧荷載作用下層狀土體的振動響應(yīng);之后,Sheng等[4]研究了考慮軌道高低不平順條件下列車荷載引起的地面振動;馮青松等[5-8]通過建立車輛—軌道—路基—地基耦合系統(tǒng)振動解析模型,研究了軌道不平順下耦合系統(tǒng)各部分的振動響應(yīng)。這些解析和半解析模型具有計算工作量小、效率高的特點(diǎn)[7-8],但通常需要滿足很多前提假定,對地基的幾何形態(tài)有嚴(yán)格限制[9],當(dāng)大地橫截面并非平行土層或包含有隧道等異質(zhì)結(jié)構(gòu)物時,該方法便很難適用。

針對地鐵列車在隧道中運(yùn)行時引起的環(huán)境振動問題,數(shù)值計算方法強(qiáng)大的分析能力使得其比解析和半解析法能更好地適應(yīng)這種復(fù)雜的結(jié)構(gòu)狀況[10]。2D 模型不能模擬波在三維空間的傳播[11-13],而3D 模型出現(xiàn)計算模型大,計算步長短的情況時,會導(dǎo)致時間成本非常大[14-15]。假定軌道和地基的幾何形狀以及材料特性沿軌道方向均勻不變,這時2.5D 模型能獲得較好應(yīng)用。由于2.5D 模型僅在平面上進(jìn)行離散,通過沿軌道方向的波數(shù)變換將3 維問題變?yōu)槠矫鎽?yīng)變問題[16],使得計算效率相比3D模型大大提高。邊學(xué)成[16]、高廣運(yùn)[17]等采用2.5D 有限元方法分析了移動軸荷載作用下的地基振動,但未考慮車輛與軌道結(jié)構(gòu)的耦合。Ama?do[18],Jin[19]和Lopes[20]等采用2.5D 有限元法研究地鐵列車引起的地面振動時,將軌道高低不平順激勵考慮為諧波不平順激勵或各諧波不平順激勵分量的疊加,但未考慮場地隨機(jī)振動響應(yīng)的統(tǒng)計特征及其變化規(guī)律。

本文基于2.5 維有限元法和虛擬激勵法,通過建立軌道—隧道—地基土2.5維數(shù)值模型,對地面的隨機(jī)振動特征以及車速和軌道不平順等級對其影響規(guī)律進(jìn)行研究。

1 基于虛擬激勵法求解輪軌力功率譜

為構(gòu)造軌道—隧道—地基土2.5 維有限元模型的外部激勵,需先求出輪軌力功率譜密度函數(shù)。車輛—軌道—隧道—地基土垂向耦合模型如圖1所示。移動車輛考慮為彈簧和阻尼器連接的多剛體系統(tǒng),車輪與鋼軌為線性赫茲彈性接觸,鋼軌為無限長的歐拉梁,對無砟軌道而言,當(dāng)頻率低于200 Hz 時,離散支撐的影響可以忽略不計[22],因此,將軌道板、隧道簡化為沿軌道方向質(zhì)量均勻分布的梁,用均布的線性彈簧和阻尼器模擬扣件、墊層以及彈性土體的剛度和阻尼。地基采用Winkle 地基模型,將其考慮為層狀半空間體,忽略地基土的振動對耦合系統(tǒng)振動的影響。

考慮黏滯阻尼與應(yīng)變阻尼的影響,建立軌道—隧道—地基土系統(tǒng)的振動微分方程[15]

式中:Er,Es和Ed分別為鋼軌、軌道板和隧道的彈性模量;Ir,Is和Id分別為鋼軌、軌道板和隧道的水平慣性矩;mr,ms和md分別為鋼軌、軌道板和隧道的單位長度質(zhì)量;zr,zs和zd分別為鋼軌、軌道板和隧道的豎向撓度;cre,cri,cse,csi,cde和cdi分別為鋼軌、軌道板、隧道的黏滯阻尼和應(yīng)變阻尼系數(shù);kp,km,kc和cp,cm,cc分別為單位長度扣件、彈性墊層、土體的剛度和阻尼;Fl為軌道不平順引起的第l個輪對的動態(tài)輪軌力;v為車速;pl為初始時刻第l個輪載距原點(diǎn)的距離,假設(shè)第1 個輪載初始時刻位于原點(diǎn)處;N為車廂數(shù)量。

黏滯阻尼系數(shù)ce和應(yīng)變阻尼系數(shù)ci與結(jié)構(gòu)前2階豎向自振頻率ω1和ω2以及結(jié)構(gòu)阻尼比ξ的關(guān)系為

式中:m為結(jié)構(gòu)的單位長度質(zhì)量;E為結(jié)構(gòu)的彈性模量。

假定軌道不平順為零均值高斯平穩(wěn)隨機(jī)過程,根據(jù)虛擬激勵法[21]原理,當(dāng)?shù)罔F列車在具有高低不平順的軌道上運(yùn)行時,各激勵之間完全相干,車輛受到多點(diǎn)異相位平穩(wěn)隨機(jī)激勵?z(t),考慮由車輪間相對位置所產(chǎn)生的各輪對隨機(jī)激勵相位差,則有

其中,

式中:tl為第l個激勵點(diǎn)距第1 個激勵點(diǎn)的滯后時間;p1為初始時刻第1個輪載距原點(diǎn)的距離,p1=0。

軌道不平順功率譜為波數(shù)kx表示的單邊功率譜Sv(kx),由于Ω=kxv,則以時間圓頻率Ω表示的單邊功率譜Sv(Ω)為

設(shè)?z(t-tl)相應(yīng)的虛擬激勵的軌道諧波不平順為(t-tl),則得

假設(shè)車輪和鋼軌為線性赫茲彈性接觸,由輪軌接觸點(diǎn)處位移限制條件,得到動態(tài)輪軌力求解方程[15]

其中,

式中:F(Ω)為與(Ω)對應(yīng)的虛擬動態(tài)輪軌力幅值;AW和AR分別為車輛在輪對處和軌道—隧道—地基土在輪軌接觸點(diǎn)的柔度矩陣;A?為輪軌線性赫茲接觸的柔度系數(shù)矩陣;Aw,j為第i節(jié)車廂的柔度矩陣;Kv,j,Cv,j和Mv,j分別為第j節(jié)車廂的剛度、阻尼和質(zhì)量矩陣;kh為輪軌間線性赫茲彈性接觸時的接觸剛度;G為撓度系數(shù);Fh為輪軌接觸彈簧剛度線性化時接觸點(diǎn)的法向力,近似取軸重;zrs,q為單位力作用于第q個輪軌接觸點(diǎn)引起的第s個輪軌接觸點(diǎn)處鋼軌位移。

由此可得動態(tài)輪軌力的功率譜

式中:(Ω)*和(Ω)T分別為(Ω)的共軛和轉(zhuǎn)置向量。

利用式(10)得到的輪軌力功率譜密度,基于虛擬激勵法構(gòu)造虛擬動態(tài)輪軌力激勵(t,Ω)

假定各輪對具有相同的幾何參數(shù)與材料參數(shù),車輛的各輪對所歷經(jīng)的軌道不平順除固定的時間差之外是完全一致的。這樣,時間—空間域中虛擬動態(tài)輪軌力的分布函數(shù)(t,Ω,x)為

在頻率—波數(shù)域中式(12)可表示為

式中:ω為頻率—波數(shù)域內(nèi)的頻率。

2 基于2.5D FEM-PEM 的場地隨機(jī)振動模型

2.1 軌道—隧道—地基土數(shù)值模型

本文基于2.5 維有限元—最佳匹配層法[20]建立軌道—隧道—地基土數(shù)值模型。最佳匹配層(Perfectly Matched Layers,PML)單元具有吸收任意方向的波而不產(chǎn)生虛擬反射的功能,且最佳匹配層區(qū)域內(nèi)的波將在兩側(cè)邊界之間衰減至零,PML單元中波的衰減如圖2所示。

圖2 PML單元中波的衰減示意圖

在2.5 維軌道—隧道—地基土有限元模型關(guān)注區(qū)域外圍施加最佳匹配層,在該特殊層的外邊緣施加Dirichlet 和Neummann 邊界條件,將關(guān)注區(qū)域和最佳匹配層區(qū)域各單元振動方程整裝,得到虛擬動態(tài)輪軌力作用下軌道—隧道—地基土系統(tǒng)頻率—波數(shù)域總振動方程[15]為

式中:KFEM(kx)和KPML(kx,ω)分別是有限元關(guān)注區(qū)域和最佳匹配層區(qū)域的總剛度矩陣;MFEM和MPML(kx,ω)分別為對應(yīng)的總質(zhì)量矩陣;(kx,Ω,ω)為頻率—波數(shù)域中軌道—隧道—地基土系統(tǒng)的虛擬位移向量;(kx,Ω,ω)為頻率—波數(shù)域中的虛擬動態(tài)輪軌力向量。

由式(14)中頻率—波數(shù)域的虛擬位移可得到虛擬速度和加速度分別為

分別將,和經(jīng)雙重傅里葉逆變換,得到相應(yīng)的時空域虛擬響應(yīng)ud,vd和ad,時空域虛擬響應(yīng)向量的共軛與轉(zhuǎn)置相乘即可得到時變功率譜。因此,地面振動加速度的時變功率譜為

式中:ad(Ω,t)*和ad(Ω,t)T分別為ad(Ω,t)的共軛和轉(zhuǎn)置向量。

利用式(17)的加速度時變功率譜對環(huán)境振動進(jìn)行評價,第i個三分之一倍頻帶的加速度有效值ai[23]為

其中,

SA(ω,t)=πSA(Ω,t)

式中:ωs、ωe分別為第i個三分之一倍頻帶的上、下限頻率,頻率范圍為1~80 Hz;SA(ω,t)為以頻率表示的地面隨機(jī)振動加速度雙邊功率譜密度。

按全身振動垂向計權(quán)因子修正后的加速度a′為

式中:Wi為第i個三分之一倍頻帶的計權(quán)系數(shù)。

Z振級計算公式為

式中:a0為基準(zhǔn)加速度,取值為1×10–6m ?s-2。

以列車通過時段內(nèi)Z 振級的最大值VLzmax作為環(huán)境振動預(yù)測評價指標(biāo),VLzmax的計算公式為

采用頻率采樣三角級數(shù)法[24]反演地面振動響應(yīng)時程R(t),得

其中,

式中:SR為地面隨機(jī)振動響應(yīng)功率譜,可為位移、速度和加速度功率譜的任1 個;R(t)和Rk分別為與SR對應(yīng)的隨機(jī)振動響應(yīng)時程和頻譜幅值;φk為0到2π 之間總數(shù)為η的隨機(jī)均勻分布序列的第k個數(shù);Ω+和Ω-分別為SR(Ω,t)的圓頻率采樣上、下限;?ω為頻率采樣間隔。

由地面隨機(jī)振動響應(yīng)功率譜SR(Ω,t)可得到地面隨機(jī)振動響應(yīng)方差σ2(t)為

將確定性軸重向量F(t)代入式(14)中,根據(jù)前面的推導(dǎo)可得出確定性軸重引起的地面振動響應(yīng)向量(即響應(yīng)均值μ),基于3σ準(zhǔn)則可得到地面振動響應(yīng)上限值(μ+3σ)。

2.2 模型驗(yàn)證

選用文獻(xiàn)[11]中的軌道—隧道—地基土參數(shù),軌道結(jié)構(gòu)采用UIC60 鋼軌,軌距為1.5 m,沿軌道方向軌枕按0.6 m 離散支撐,厚和寬為0.4 和2.5 m 的軌道板。混凝土軌道板的彈性模量為28.5 GPa,混凝土軌道板的密度為2 500 kg ?m-3,隧道的中心距地表13.5 m,隧道直徑5.4 m,襯砌厚度0.3 m,具體參數(shù)如圖3所示。建立如圖4所示的軌道—隧道—地基土2.5D 有限元模型,分別采用解析法和2.5D有限元法計算鋼軌的柔度曲線,通過在零頻率處匹配鋼軌響應(yīng)得到等效的軌道—隧道解析模型參數(shù)[15]。用解析法和2.5D 有限元法得到的鋼軌柔度系數(shù)曲線如圖5所示,軌道—隧道解析模型等效參數(shù)見表1。

圖3 文獻(xiàn)[11]中的隧道—地基土模型參數(shù)

圖4 軌道—隧道—土體有限元模型示意圖

圖5 鋼軌柔度曲線

文獻(xiàn)[11]中地鐵列車為6節(jié)相同車廂的編組形式,車輛參數(shù)見表2。圖6為利用式(22)反演的仰拱處(圖3中B點(diǎn))垂向振動速度時程。該結(jié)果與文獻(xiàn)[11]結(jié)果吻合較好,驗(yàn)證了本文模型的可靠性。

3 地面振動影響

3.1 車速的影響

地鐵列車運(yùn)行速度分別選取40,60 和80 km·h-1,每次計算僅改變列車速度,軌道不平順激勵采用美國聯(lián)邦鐵路局(Federal Railroad Administra?tion,F(xiàn)RA)5 級軌道高低不平順功率譜,其他參數(shù)與2.2 節(jié)相同,得到的第1 個輪對處的輪軌力功率譜密度如圖7所示。

表1 軌道—隧道解析模型等效參數(shù)

表2 車輛參數(shù)

由圖7可知:輪軌力功率譜密度隨著速度的增加而顯著增加。車速為60 和80 km·h-1時,主頻f1=0.9 Hz處動態(tài)輪軌力功率譜密度分別比40 km·h-1時增大5.22和7.39倍。

運(yùn)用虛擬激勵法,將圖7的輪軌力功率譜密度作為不同速度條件下軌道—隧道—地基土系統(tǒng)2.5D有限元模型的外部激勵,進(jìn)行數(shù)值計算。

圖6 仰拱處垂向振動速度響應(yīng)時程

圖8為距軌道中心線20 m 處(圖3中A 點(diǎn))振動響應(yīng)的最大功率譜密度。由圖8可知:地面振動位移以低頻振動為主;振動速度和加速度以高頻(40~60 Hz)振動為主且主頻相同,車速為40,60 和80 km·h-1時的主頻分別為45.3,46.1 和50.0 Hz,主頻隨著車速的增加顯著增大。

圖7 不同速度下第1個輪對處的輪軌力功率譜密度

圖8 不同速度下A點(diǎn)振動響應(yīng)最大功率譜密度

圖9 不同速度下A點(diǎn)垂向振動響應(yīng)標(biāo)準(zhǔn)差

不同車速條件下A 點(diǎn)的振動響應(yīng)標(biāo)準(zhǔn)差σ和均值μ分別如圖9和圖10所示。由圖9可知:地面振動響應(yīng)標(biāo)準(zhǔn)差隨著車速增加顯著增大,車速為40和80 km·h-1時的振動加速度標(biāo)準(zhǔn)差最大值分別為0.008 9 和0.025 2 m·s-2,表明車速越大,地面振動響應(yīng)的離散性越大。

將圖9中A 點(diǎn)振動響應(yīng)標(biāo)準(zhǔn)差與圖10的均值結(jié)果對比可知:地面振動響應(yīng)上限值(μ+3σ)受確定性軸重和動態(tài)輪軌力的雙重影響。地面振動位移上限值主要由確定性軸重引起,振動加速度上限值主要由動態(tài)輪軌力引起。

圖11給出了地鐵運(yùn)行引起的場地振動上限值。由圖11可知:車速的增加對場地的位移上限值影響很小,最大位移上限值出現(xiàn)在軌道中心線處;振動速度和加速度上限值隨著車速的增加顯著增大。隨機(jī)動態(tài)輪軌力作用下,地面的振動速度和加速度上限值在垂直于地鐵運(yùn)行方向的衰減出現(xiàn)明顯波動,最大值并不位于軌道中心線處。

圖10 不同速度下A點(diǎn)垂向振動響應(yīng)均值

圖11 場地垂向振動響應(yīng)上限值

依據(jù)GB 10070-1988《城市區(qū)域環(huán)境振動標(biāo)準(zhǔn)》[25]的規(guī)定,要求居民、文教區(qū)晝間的振動容許限值為70 dB。圖12給出了地鐵通過過程中地面振動Z 振級最大值。由圖12可知:地鐵以80 km·h-1運(yùn)行時,除距軌道中心線10~15 m范圍內(nèi)超限,其他區(qū)域符合規(guī)范要求;對于40 和60 km·h-1的運(yùn)行速度,各距離處的振動均符合規(guī)范要求。

圖12 不同速度下的Z振級最大值

3.2 軌道不平順的影響

地鐵列車運(yùn)行速度為60 km·h-1,軌道不平順激勵選取FRA 4 級,F(xiàn)RA 5 級和FRA 6 級軌道高低不平順功率譜,每次計算僅改變軌道不平順等級,其他參數(shù)與2.2 節(jié)相同,得到的第1 個輪對處的輪軌力功率譜密度如圖13所示。

圖13 不同軌道不平順下第1個輪對處的輪軌力功率譜密度

由圖13可知:0.9 Hz 處的動態(tài)輪軌力功率譜密度在軌道不平順條件為FRA 5 級和FRA 4 級時,相比FRA 6 級分別增大了6.18 和15.9 倍,這表明軌道不平順加劇了車速對動態(tài)輪軌力的影響。

結(jié)合虛擬激勵法,將圖13的輪軌力功率譜密度作為不同軌道不平順條件下軌道—隧道—地基土系統(tǒng)2.5D有限元模型的外部激勵,進(jìn)行數(shù)值計算。

地面A 點(diǎn)不同軌道不平順條件下振動響應(yīng)的最大功率譜密度曲線如圖14所示。由圖14可知:軌道不平順等級從FRA 6 級變?yōu)镕RA 4 級,振動位移、速度、加速度功率譜最大值均增大15.9倍。不同軌道不平順下振動響應(yīng)的最大功率譜密度曲線波形近似,幅值不同,說明軌道不平順不改變地面振動響應(yīng)功率譜的頻譜分布特性。

圖14 不同軌道不平順下A點(diǎn)振動響應(yīng)的最大功率譜

圖15為不同軌道不平順條件下A 點(diǎn)的振動響應(yīng)標(biāo)準(zhǔn)差。由圖15可知:隨軌道不平順等級的降低顯著增大,地面隨機(jī)振動響應(yīng)的離散性明顯增大,軌道不平順加劇了車速對振動響應(yīng)離散性的影響。

圖15 不同軌道不平順下A點(diǎn)垂向振動響應(yīng)標(biāo)準(zhǔn)差

圖16給出了不同軌道不平順條件下場地的振動響應(yīng)上限值。由圖16可知,軌道不平順等級的降低對場地的位移上限值幾乎不產(chǎn)生影響,位移最大值出現(xiàn)在軌道中心線處。結(jié)合圖11可知,振動速度和加速度上限值隨著軌道不平順等級的降低顯著增加,兩者最大值分別出現(xiàn)在距軌道中心線約13 和14 m 處。距軌道中心線60 m 外,地面振動速度和加速度的衰減趨勢變緩。

圖16 場地垂向振動響應(yīng)上限值

圖17為不同軌道不平順下的Z 振級最大值。由圖12和圖17可知:由于土體振動加速度在垂直于地鐵運(yùn)行方向的衰減出現(xiàn)明顯波動現(xiàn)象,VLzmax曲線最大值出現(xiàn)在距軌道中心線13.5 m 處。距軌道中心線60 m 外,Z 振級最大值的衰減趨勢變緩。當(dāng)車速為60 km ?h-1,由于FRA 4 級軌道高低不平順的軌道平順性較差,除距軌道中心線8~20 m 范圍超限,其他區(qū)域符合規(guī)范要求;對于FRA 5 級和FRA 6 級軌道高低不平順條件,各距離處的振動均滿足規(guī)范要求。

圖17 不同軌道不平順下的Z振級最大值

4 結(jié) 論

(1)地面的振動位移主要取決于車輛的確定性軸重激勵,以低頻振動為主;地面振動速度和加速度以高頻(40~60 Hz)振動為主。地鐵車速和軌道不平順等級的改變對地面振動位移上限值的影響很小,對速度和加速度上限值影響顯著,其中加速度上限值主要由軌道不平順激勵引起。

(2)軌道不平順加劇了地鐵車速對地面振動響應(yīng)標(biāo)準(zhǔn)差的影響。車速越高,軌道不平順等級越低,地面振動響應(yīng)的離散性越高,從而導(dǎo)致振動響應(yīng)上限值越大。

(3)地鐵車速的增加導(dǎo)致振動速度和加速度主頻增大;軌道不平順等級僅改變振動響應(yīng)頻譜幅值,不改變頻譜分布特性。

(4)軌道不平順隨機(jī)激勵下,地面振動速度和加速度的上限值以及Z振級最大值在垂直于地鐵運(yùn)行方向的衰減出現(xiàn)明顯波動,距軌道中心線60 m外衰減趨勢變緩。

猜你喜歡
振動模型
一半模型
振動的思考
噴水推進(jìn)高速艇尾部振動響應(yīng)分析
重要模型『一線三等角』
This “Singing Highway”plays music
重尾非線性自回歸模型自加權(quán)M-估計的漸近分布
振動攪拌 震動創(chuàng)新
中國公路(2017年18期)2018-01-23 03:00:38
中立型Emden-Fowler微分方程的振動性
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉(zhuǎn)換方法初步研究
主站蜘蛛池模板: 国产va在线观看| 亚洲天堂高清| av在线5g无码天天| 2021国产精品自拍| 亚亚洲乱码一二三四区| 四虎影视8848永久精品| 久久午夜影院| 亚洲第一福利视频导航| 天天综合色网| 九九香蕉视频| 久久精品最新免费国产成人| 在线一级毛片| 久久精品最新免费国产成人| 天堂成人av| 成人无码一区二区三区视频在线观看| 国产免费a级片| 手机精品福利在线观看| 国产爽爽视频| 超碰aⅴ人人做人人爽欧美| 亚洲欧美另类色图| 国产精品无码作爱| 国产成人艳妇AA视频在线| 亚洲无码久久久久| 高清不卡毛片| 国产成人高清精品免费| 一本一道波多野结衣一区二区| 全部无卡免费的毛片在线看| 欧美一区二区人人喊爽| 一区二区欧美日韩高清免费| 日韩福利视频导航| 国产一区二区网站| 日韩欧美在线观看| 依依成人精品无v国产| 欧美成人a∨视频免费观看| 国产精品嫩草影院av| 亚洲小视频网站| 在线观看热码亚洲av每日更新| 成人综合久久综合| 夜夜高潮夜夜爽国产伦精品| 国产不卡一级毛片视频| 免费国产在线精品一区| 欧美日韩免费观看| 亚洲人成在线精品| 日本久久免费| 日韩小视频在线播放| 国产精品观看视频免费完整版| 刘亦菲一区二区在线观看| 天堂va亚洲va欧美va国产| 国产欧美日韩另类| 中文成人在线| 波多野结衣亚洲一区| 992Tv视频国产精品| 国产JIZzJIzz视频全部免费| 婷婷五月在线| 欧美啪啪网| 99r在线精品视频在线播放 | 国产成人精品免费视频大全五级| 国产成人精彩在线视频50| 久久久久免费精品国产| 国产特级毛片aaaaaa| 国产午夜人做人免费视频| 婷婷综合色| 波多野结衣一区二区三区88| 2021国产精品自拍| 久久天天躁狠狠躁夜夜躁| 18禁黄无遮挡免费动漫网站| 在线播放国产一区| swag国产精品| 日本一区二区不卡视频| 99热免费在线| 国产精品国产三级国产专业不| 午夜性刺激在线观看免费| 999精品色在线观看| 日韩在线2020专区| 亚洲精品福利视频| a在线观看免费| 巨熟乳波霸若妻中文观看免费 | 亚洲免费毛片| 日韩高清欧美| 伊人久久久大香线蕉综合直播| 亚洲国产日韩在线观看| 国产日韩欧美视频|