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

車載作用下公路橋梁耦合振動(dòng)精細(xì)化建模及驗(yàn)證分析

2021-10-06 08:39:56殷新鋒晏萬(wàn)里任厚乾劉揚(yáng)
關(guān)鍵詞:有限元橋梁振動(dòng)

殷新鋒,晏萬(wàn)里,任厚乾,劉揚(yáng),2

(1.長(zhǎng)沙理工大學(xué) 土木工程學(xué)院,湖南 長(zhǎng)沙 410114;2.湖南工業(yè)大學(xué) 土木工程學(xué)院,湖南 株洲 412007)

車輛在橋梁上行駛時(shí)產(chǎn)生的激勵(lì)導(dǎo)致車輛和橋梁產(chǎn)生相互振動(dòng),該振動(dòng)稱為車橋耦合振動(dòng)[1].國(guó)內(nèi)外學(xué)者對(duì)車橋耦合振動(dòng)問(wèn)題進(jìn)行了大量研究,并取得非凡成果.夏禾[2]將車輛簡(jiǎn)化為懸掛振動(dòng)模型,研究了車-橋-墩相互作用的動(dòng)力響應(yīng).Cai 等[3-4]采用兩軸車輛模型,基于功率譜密度函數(shù)生成隨機(jī)橋面不平整度,分析了車輛荷載作用下橋梁振動(dòng)響應(yīng).Huang等[5]采用三軸車輛有限元模型研究了簡(jiǎn)支梁的車橋隨機(jī)振動(dòng).韓萬(wàn)水等[6]結(jié)合實(shí)測(cè)數(shù)據(jù)對(duì)橋梁模型進(jìn)行修正以得到橋梁基準(zhǔn)模型,并采用梁格法對(duì)車橋振動(dòng)進(jìn)行分析.李奇等[7]考慮車體柔性的影響,分析了高速列車和簡(jiǎn)支梁橋相互作用的影響.鄧露等[8]運(yùn)用動(dòng)態(tài)稱重技術(shù)(BWIM)實(shí)時(shí)監(jiān)測(cè)車輛總重和軸重,進(jìn)行中小型跨徑混凝土梁橋的振動(dòng)研究.趙越等[9]基于等參映射與改進(jìn)折半法在傳統(tǒng)車橋耦合分析的基礎(chǔ)上進(jìn)一步提升其分析精度及計(jì)算效率,從而進(jìn)行公路車橋耦合分析.但關(guān)于充氣輪胎胎壓對(duì)車-橋耦合振動(dòng)影響的研究較少.主要原因?yàn)楝F(xiàn)有車-橋耦合振動(dòng)分析中車輛模型多為簡(jiǎn)化的質(zhì)量-彈簧-阻尼多自由度振動(dòng)模型[10-13],該模型常將車輪簡(jiǎn)化為點(diǎn)或者等效線面接觸,故不能精確考慮車輛動(dòng)力特性和柔性輪胎對(duì)車橋耦合振動(dòng)響應(yīng)的影響.因此,建立精確車輛模型和橋梁模型是至關(guān)重要的,這樣才能提高數(shù)值模擬精度,為橋梁結(jié)構(gòu)的健康運(yùn)營(yíng)提供有利建議.

本文以一座連續(xù)剛構(gòu)箱梁橋?yàn)楣こ瘫尘?首先,基于LS-DYNA 程序以車輛實(shí)際構(gòu)造及動(dòng)力特性為基準(zhǔn),建立車輛精細(xì)三維有限元模型;然后,結(jié)合響應(yīng)面法并依據(jù)實(shí)橋試驗(yàn)結(jié)果對(duì)橋梁模型進(jìn)行修正,以得高精度的橋梁有限元模型;最后,通過(guò)LS-DYNA 程序?qū)蛄耗P秃蛙囕v模型進(jìn)行耦合,求解車橋的振動(dòng)響應(yīng),并通過(guò)改變車輪氣壓,分析該參數(shù)對(duì)橋梁振動(dòng)響應(yīng)的影響.

1 車輛模型的建立及校驗(yàn)

1.1 車輛模型的建立

在車橋耦合振動(dòng)模型中,影響橋梁振動(dòng)響應(yīng)的主要因素為車輛模型的動(dòng)力特性和荷載分布.車輛模型需著重模擬懸架系統(tǒng)、車輪和軸重.本文參照東風(fēng)牌三軸載重自卸貨車,建立車輛有限元模型,其相關(guān)參數(shù)如表1 所示.前中軸距為3.5 m,中后軸距為1.4 m,后輪距為1.8 m.基于CAD 軟件SolidWorks 分別建立各部位的幾何模型并進(jìn)行網(wǎng)格劃分,最后使用梁、殼、實(shí)體單元及離散單元等賦予不同部位網(wǎng)格不同的屬性,從而構(gòu)建車輛模型,如圖1 所示.

表1 車輛模型參數(shù)Tab.1 Vehicle model parameters

圖1 車輛的三維有限元模型Fig.1 Three dimensional finite element model of vehicle

車輪由輪胎、輪盤和輪轂組成.輪胎采用線彈性橡膠材料,輪盤和輪轂采用線彈性鋼材材料,車輪模型中各部位連接均為剛性連接,邊界為剛性固態(tài)約束,有限元模型如圖2 所示.根據(jù)輪胎內(nèi)實(shí)際氣壓,使 用 LS -DYNA 程 序 關(guān) 鍵 字 *AIRBAG_SIMPLE_AIRBAG_MODEL 定義由輪胎、輪轂組成的封閉體內(nèi)的氣壓,其氣體壓力值為0.6 MPa.

圖2 車輪模型的組成Fig.2 Composition of wheel model

由于車輛懸架系統(tǒng)結(jié)構(gòu)復(fù)雜,本文采用殼單元、剛體、彈簧阻尼單元和多點(diǎn)約束來(lái)模擬懸架系統(tǒng).前后懸架的彈簧和阻尼器的參數(shù)參考文獻(xiàn)[14],前懸架的彈簧剛度取800 N/mm,阻尼系數(shù)取20 Ns/mm;后懸架的彈簧剛度取1 200 N/mm,阻尼系數(shù)取25 Ns/mm.前后懸架實(shí)體構(gòu)造及有限元模型如圖3 所示.

圖3 懸架實(shí)體構(gòu)造及FE 模型Fig.3 Suspension structure and FE model

為了使車輪正常轉(zhuǎn)動(dòng),在輪盤與車軸之間設(shè)置旋轉(zhuǎn)鉸,再通過(guò)定義*INITIAL_VELOCTIY_GENERATION 關(guān)鍵字設(shè)置車輪的轉(zhuǎn)動(dòng)和平動(dòng)速度,實(shí)現(xiàn)車輪滾動(dòng)向前的狀態(tài),如圖4 所示.

圖4 前輪上一點(diǎn)的轉(zhuǎn)動(dòng)軌跡曲線Fig.4 Curve of turning track of a point on the front wheel

1.2 車輛模型的校驗(yàn)

在車輛模型的3 個(gè)軸上選擇6 個(gè)節(jié)點(diǎn),約束豎向位移.然后對(duì)車輛施加重力荷載使車輛產(chǎn)生瞬時(shí)振動(dòng),再通過(guò)定義關(guān)鍵字*DAMPING_GLOBAL 對(duì)車輛模型施加全局阻尼,使其快速達(dá)到穩(wěn)定狀態(tài),計(jì)算出約束反力,并與實(shí)測(cè)車輛軸重進(jìn)行比較,最終得到車輛模型的軸載曲線如圖5 所示.

由表2 知,車體總重偏差為0.87%,說(shuō)明實(shí)測(cè)車輛軸重分布特性能體現(xiàn)于建立的車輛模型.

通過(guò)對(duì)比車輛模型和實(shí)測(cè)車輛的自振頻率,驗(yàn)證車輛模型動(dòng)力特性的有效性.如圖6 所示,該車實(shí)測(cè)自振頻率為1.635 Hz,車輛有限元模型的自振頻率為1.647 Hz,顯然,兩者基頻相差很小.

圖6 車輛自振頻率Fig.6 Natural frequency of vehicle

2 橋梁模型

2.1 橋梁概況

該橋主橋?yàn)槿缱兘孛骖A(yù)應(yīng)力混凝土連續(xù)剛構(gòu)結(jié)構(gòu),跨徑布置為(65+120+65)m,其立面及斷面示意圖分別如圖7 和圖8 所示.

圖7 橋梁立面示意圖(單位:m)Fig.7 Elevation diagram of bridge(unit:m)

圖8 橋梁斷面示意圖(單位:cm)Fig.8 Schematic diagram of bridge section(unit:cm)

2.2 實(shí)橋試驗(yàn)

利用有限元軟件ANSYS 建立初始有限元模型如圖9 所示.主梁和墩身都用實(shí)體單元模擬,墩底采用固結(jié)連接.混凝土密度取為2 500 kg/m3,主梁和墩的混凝土彈性模量分別為34.5 GPa 和32.5 GPa.

圖9 橋梁有限元模型Fig.9 Finite element model of bridge

現(xiàn)場(chǎng)對(duì)橋梁進(jìn)行靜力試驗(yàn),測(cè)試中加載車輛總數(shù)為6 輛,每輛總重為35 t,車輛前軸重7 t,中軸和后軸各重14 t,車輛照片如圖10 所示.

圖10 試驗(yàn)車輛Fig.10 Test vehicle

在正式試驗(yàn)前先進(jìn)行預(yù)加載,以消除非彈性變形,確保試驗(yàn)及設(shè)備處于良好工作狀態(tài).正式試驗(yàn)中將加載車分3 級(jí)加載,每一級(jí)加載持荷2~3 min,待實(shí)測(cè)應(yīng)變及撓度數(shù)據(jù)穩(wěn)定后進(jìn)行數(shù)據(jù)的采集工作,再進(jìn)行下一級(jí)的加載工作.限于篇幅,僅就其中一個(gè)工況做詳要概述.為了使加載截面承受最大正彎矩,在中跨跨中截面中心加載,相應(yīng)的荷載效率系數(shù)為0.967.撓度測(cè)點(diǎn)布置如圖11 所示.

圖11 撓度測(cè)點(diǎn)布置示意圖(單位:m)Fig.11 Layout of deflection measuring points(unit:m)

實(shí)橋試驗(yàn)的撓度采用水準(zhǔn)儀進(jìn)行測(cè)量,現(xiàn)場(chǎng)測(cè)試照片如圖12 所示.脈動(dòng)試驗(yàn)采用多通道數(shù)據(jù)采集分析系統(tǒng)NI 公司的PXI 系統(tǒng)進(jìn)行試驗(yàn),采用8330B3 型超低頻加速度傳感器進(jìn)行數(shù)據(jù)采集,最低采樣頻率從0 Hz 開始,采集主梁橫向和豎向振動(dòng)數(shù)據(jù),再經(jīng)信號(hào)分析得到全橋的各階固有振動(dòng)特性,現(xiàn)場(chǎng)采集照片如圖13 所示.

圖12 現(xiàn)場(chǎng)撓度測(cè)試Fig.12 Field deflection test

圖13 現(xiàn)場(chǎng)動(dòng)態(tài)測(cè)試Fig.13 Field dynamic test

2.3 基于響應(yīng)面法的有限元模型更新

基于響應(yīng)面法的有限元模型修正是用響應(yīng)面函數(shù)來(lái)模擬實(shí)際結(jié)構(gòu)的響應(yīng)函數(shù),將試驗(yàn)設(shè)計(jì)與數(shù)理統(tǒng)計(jì)相結(jié)合,通過(guò)樣本選取、方差分析參數(shù)選取、響應(yīng)面的擬合及采用優(yōu)化算法尋求響應(yīng)面模型中的最優(yōu)解來(lái)進(jìn)行有限元模型修正.

據(jù)文獻(xiàn)[15-16]知,影響有限元模型與實(shí)測(cè)橋梁結(jié)構(gòu)差別的主要因素為結(jié)構(gòu)混凝土密度、主梁和墩混凝土彈性模量,因此取這三個(gè)因素作為變量.為獲得響應(yīng)與所選定的三參數(shù)之間的聯(lián)系,首先需要參數(shù)設(shè)計(jì).根據(jù)參數(shù)取值的變化規(guī)律,假設(shè)三參數(shù)的單位長(zhǎng)度值為10%,則可得該三參數(shù)值的變化范圍見表3.

表3 各次試驗(yàn)參數(shù)的取值Tab.3 Values of parameters in each experiment

以橋梁模態(tài)和靜力變形為目標(biāo)來(lái)更新橋梁模型,選取橋梁結(jié)構(gòu)的第一階自振頻率(R1)、測(cè)點(diǎn)5 撓度值(R2)和測(cè)點(diǎn)6 撓度值(R3)作為目標(biāo)函數(shù).這3 個(gè)目標(biāo)函數(shù)充分利用了現(xiàn)場(chǎng)的實(shí)橋試驗(yàn)條件,且包含了橋梁靜、動(dòng)力性能指標(biāo),可較為全面、準(zhǔn)確地反映橋梁的力學(xué)性能.

根據(jù)上述試驗(yàn)設(shè)計(jì),應(yīng)用回歸分析技術(shù)對(duì)樣本數(shù)據(jù)進(jìn)行響應(yīng)面擬合,并用R2準(zhǔn)則和準(zhǔn)則進(jìn)行響應(yīng)面擬合精度的檢驗(yàn),其表達(dá)式見方程(1)(2).

式中:SST=SSE+SSR表示模型的總方差;dT表示模型的總自由度.如果R2和的值都接近1 且兩者差值很小,則表示響應(yīng)面方程擬合得很好.

采用二次多項(xiàng)式對(duì)上述樣本結(jié)果進(jìn)行響應(yīng)面擬合,并用R2和準(zhǔn)則對(duì)響應(yīng)面的精度進(jìn)行檢驗(yàn),檢驗(yàn)結(jié)果如表4 所示.

表4 響應(yīng)面精度檢驗(yàn)表Tab.4 Accuracy test table of response surface

從表4 知,R2和值均接近1,且兩者差值很小,說(shuō)明模型擬合精度較高,能很好地表示各參數(shù)與目標(biāo)函數(shù)的關(guān)系.

響應(yīng)面函數(shù)采取省略交差項(xiàng)影響的二次多項(xiàng)式,則響應(yīng)面擬合方程如方程(3)~(5)所示:

該橋的第一階實(shí)測(cè)頻率為0.867 Hz;實(shí)橋跨中試驗(yàn)加載的情形下,測(cè)點(diǎn)5、6 的撓度值分別為5.687 mm 和9.679 mm.根據(jù)實(shí)測(cè)值和數(shù)值模擬值可以建立一個(gè)目標(biāo)函數(shù):

式中:R(i)為根據(jù)響應(yīng)面法計(jì)算的響應(yīng);M(i)為實(shí)測(cè)值;coef(i)為權(quán)重系數(shù),此處取為[1.5 1 1].目標(biāo)函數(shù)可以利用Matlab 中的遺傳算法(Genetic Algorithm)工具箱優(yōu)化.更新后的上述三參數(shù)分別如表5 所示.

表5 更新后參數(shù)值與初始值的比較Tab.5 Comparison of updated parameter values and initial values

從表5 可得,主梁和墩的混凝土彈性模量和密度都有所增加,說(shuō)明依托工程的施工質(zhì)量較為可靠,同時(shí)參數(shù)的變動(dòng)范圍仍在工程實(shí)踐的常見范圍之內(nèi),說(shuō)明了修正結(jié)構(gòu)的可靠性.利用所更新的參數(shù)計(jì)算得到大橋的第一階振動(dòng)頻率及上述工況的試驗(yàn)加載測(cè)點(diǎn)撓度如表6 所示.

表6 數(shù)值模型和實(shí)測(cè)值的比較Tab.6 Comparison of numerical model and measured value

3 振動(dòng)響應(yīng)對(duì)比分析

3.1 現(xiàn)場(chǎng)測(cè)點(diǎn)布置及荷載工況

該連續(xù)剛構(gòu)橋主要對(duì)主箱梁在如圖14 所示截面1、2 處(即邊跨跨中、中跨跨中)的動(dòng)位移與應(yīng)變變化情況進(jìn)行了分析,全橋共計(jì)布設(shè)動(dòng)位移測(cè)點(diǎn)5個(gè),應(yīng)變測(cè)點(diǎn)18 個(gè).

圖14 動(dòng)載分析立面示意圖(單位:m)Fig.14 Vertical sketch of dynamic load analysis(unit:m)

現(xiàn)場(chǎng)布置了3 種荷載工況:工況1,一輛35 t 重車置于橋縱向中軸線上;工況2,兩輛35 t 重車并排對(duì)稱置于橋縱向中軸線兩側(cè),兩車橫向間距為1.3 m;工況3,四輛35 t 重車以橋縱向中軸線對(duì)稱放置,車輛橫向間距為1.3 m,縱向間距為4 m.所有工況,車輛均以40 km/h 的速度從橋外駛?cè)氩⑼ㄟ^(guò)橋梁.限于篇幅,僅列出工況3 的車輛布置圖,如圖15 所示.

圖15 工況3 車輛布置示意圖Fig.15 Vehicle layout of condition three

3.2 結(jié)果對(duì)比分析

3.2.1 動(dòng)位移對(duì)比

將上述各個(gè)工況進(jìn)行對(duì)比分析,邊/中跨跨中的位移實(shí)測(cè)值與計(jì)算值對(duì)比如圖16 和圖17 所示.

圖16 各工況邊跨跨中撓度時(shí)程曲線Fig.16 Time history curve of side span mid span deflection under various working conditions

圖17 各工況中跨跨中撓度時(shí)程曲線Fig.17 Time history curve of mid span deflection under various working conditions

由圖16 和圖17 可知,計(jì)算值與實(shí)測(cè)值隨時(shí)間變化的波動(dòng)規(guī)律一致,但存在同一時(shí)刻計(jì)算值與實(shí)測(cè)值吻合不是很好的情況,且計(jì)算值的曲線波峰稍大于實(shí)測(cè)值,兩者靜態(tài)豎向位移差值在1.71%~6.55%之間,振動(dòng)位移(消除靜態(tài)位移值)幅值差值在3.71%~6.35%之間.其誤差主要原因可能是:1)橋梁模型雖經(jīng)過(guò)合理的參數(shù)更新,但與實(shí)測(cè)橋梁結(jié)構(gòu)特性還有一定差別;2)實(shí)橋試驗(yàn)受車輛速度控制、行駛軌跡、儀器誤差等因素的影響,測(cè)量數(shù)據(jù)精度有較大的誤差.

3.2.2 動(dòng)應(yīng)變對(duì)比

將上述各個(gè)工況進(jìn)行對(duì)比分析,邊、中跨跨中的動(dòng)應(yīng)變實(shí)測(cè)值與計(jì)算值對(duì)比如圖18 和圖19 所示.

圖18 各工況邊跨跨中應(yīng)變時(shí)程曲線Fig.18 Time history curve of strain in side span under various working conditions

圖19 各工況中跨跨中應(yīng)變時(shí)程曲線Fig.19 Time history curve of strain in mid span under various working conditions

由圖18 和圖19 可知,計(jì)算值與實(shí)測(cè)值兩者的變化規(guī)律一致,且數(shù)值模擬結(jié)果與實(shí)測(cè)數(shù)據(jù)吻合得較好.但計(jì)算值的曲線波峰稍大于實(shí)測(cè)值,兩者動(dòng)應(yīng)變的差值在3.35%~6.94%之間,振動(dòng)應(yīng)變(消除靜態(tài)應(yīng)變值)幅值差值在2.83%~5.96%.

從上述對(duì)比情況可知,豎向振動(dòng)位移與動(dòng)應(yīng)變的計(jì)算值與實(shí)測(cè)值隨時(shí)間的變化趨勢(shì)一致,且整體吻合較好,峰值接近,說(shuō)明該法建立的橋梁模型和車輛模型是有效的,能應(yīng)用于車橋耦合計(jì)算中.

4 輪胎氣壓對(duì)橋梁振動(dòng)的影響

輪胎氣壓的大小與輪胎變形成幾何關(guān)系,因此,輪胎氣壓對(duì)跨中振動(dòng)響應(yīng)影響反映了輪胎變形對(duì)梁跨中振動(dòng)的影響.以上述工況2 為例,取3 種分別為P1=0.4 MPa、P2=0.6 MPa、P3=0.8 MPa 的不同氣壓,來(lái)分析車輪氣壓變化對(duì)中跨跨中動(dòng)位移和動(dòng)應(yīng)變的影響.

從圖20 和圖21 可得,當(dāng)車輪氣壓由0.6 MPa 減少到0.4 MPa 時(shí),車輪柔性變大,與橋面的接觸面增大,橋梁中跨跨中動(dòng)位移峰值減小3.6%,動(dòng)應(yīng)變減小2.7%;當(dāng)車輪氣壓由0.6 MPa 增加到0.8 MPa 時(shí),車輪剛度變大,與橋面的接觸面減少,橋梁中跨跨中動(dòng)位移峰值增大4.4%,動(dòng)應(yīng)變?cè)龃?.85%.

圖20 中跨跨中動(dòng)位移對(duì)比曲線圖Fig.20 Comparison curve of mid span dynamic displacement

圖21 中跨跨中動(dòng)應(yīng)變對(duì)比曲線圖Fig.21 Comparison curve of dynamic strain in mid span

5 結(jié)論

本文基于LS-DYNA 程序建立了車-橋相互作用模型.設(shè)立多種荷載工況,通過(guò)計(jì)算值與實(shí)測(cè)值的對(duì)比,驗(yàn)證了車輛模型的有效性及該法建立車橋耦合模型的可行性,并進(jìn)一步分析了車輪氣壓變化對(duì)車橋振動(dòng)的影響.結(jié)果表明:1)本文方法可較真實(shí)再現(xiàn)車-橋耦合振動(dòng)情形,無(wú)需自編復(fù)雜的車橋耦合分析程序,只需定義LS-DYNA 軟件自帶的面面接觸關(guān)鍵字,就可以實(shí)現(xiàn)車橋間的耦合作用;2)可方便獲得橋梁任意部位的應(yīng)力、應(yīng)變及其他響應(yīng)數(shù)據(jù);3)可高精度分析車輪氣壓變化及車輛其他參數(shù)對(duì)耦合系統(tǒng)振動(dòng)影響.

猜你喜歡
有限元橋梁振動(dòng)
振動(dòng)的思考
振動(dòng)與頻率
手拉手 共搭愛的橋梁
句子也需要橋梁
中立型Emden-Fowler微分方程的振動(dòng)性
高性能砼在橋梁中的應(yīng)用
磨削淬硬殘余應(yīng)力的有限元分析
UF6振動(dòng)激發(fā)態(tài)分子的振動(dòng)-振動(dòng)馳豫
基于SolidWorks的吸嘴支撐臂有限元分析
箱形孔軋制的有限元模擬
上海金屬(2013年4期)2013-12-20 07:57:18
主站蜘蛛池模板: 国产精品午夜电影| 国产成人精品高清不卡在线| 精品视频在线一区| 国产日韩久久久久无码精品| 亚洲精品在线影院| 99精品欧美一区| 日韩国产精品无码一区二区三区 | 欧美一级夜夜爽| 久热99这里只有精品视频6| 久爱午夜精品免费视频| 成人免费黄色小视频| 亚洲第一成年人网站| 国产jizz| 国产欧美视频综合二区 | 欧美国产在线看| 中文字幕亚洲精品2页| 欧美精品v| 91久久国产热精品免费| 国产在线高清一级毛片| 亚洲日本中文字幕乱码中文| 亚洲国产中文欧美在线人成大黄瓜 | 91福利一区二区三区| 99久久99视频| 亚洲中文无码h在线观看| 国产情侣一区二区三区| 人妻21p大胆| 欧美日韩v| 国产精品亚洲αv天堂无码| 91www在线观看| 国产无遮挡猛进猛出免费软件| 亚洲综合第一区| 91视频日本| 国产成人AV综合久久| 亚洲成人黄色在线观看| 久久亚洲高清国产| 欧美性爱精品一区二区三区| 亚洲中文字幕无码爆乳| 香蕉久人久人青草青草| 欧美色99| 爆操波多野结衣| 国产成人亚洲无码淙合青草| 国产精品女同一区三区五区| 91久久大香线蕉| 米奇精品一区二区三区| 国产精品成人一区二区| 亚洲VA中文字幕| 美女无遮挡免费视频网站| 亚洲成人免费看| 91免费在线看| 日韩欧美中文| 女人18毛片水真多国产| 亚洲三级视频在线观看| 亚洲日本中文字幕天堂网| 久久综合丝袜长腿丝袜| 国产av一码二码三码无码| 国产视频一二三区| 九九这里只有精品视频| 午夜免费视频网站| 91在线播放国产| 国产成人一区免费观看| 国产精品第一区| 日韩欧美综合在线制服| 亚洲综合婷婷激情| 国产欧美精品专区一区二区| 欧美区在线播放| 国产综合日韩另类一区二区| 91在线一9|永久视频在线| 色妺妺在线视频喷水| 亚洲高清资源| 无码AV日韩一二三区| 白浆视频在线观看| 伊人久久久大香线蕉综合直播| 国产美女无遮挡免费视频| 成人91在线| 国产免费a级片| 99热国产这里只有精品无卡顿" | 日韩a级片视频| 69综合网| 性色一区| 亚洲三级a| 久爱午夜精品免费视频| av尤物免费在线观看|