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

重型車輛與橋梁耦合振動及橋梁沖擊系數的數值模擬研究

2014-06-07 07:15:48
結構工程師 2014年5期
關鍵詞:有限元橋梁模型

王 娟 錢 江

(同濟大學結構工程與防災研究所,上海200092)

重型車輛與橋梁耦合振動及橋梁沖擊系數的數值模擬研究

王 娟*錢 江

(同濟大學結構工程與防災研究所,上海200092)

為了更合理地評估橋梁在重型車輛作用下的耦合動力響應,在有限元軟件LS-DYNA的平臺上,根據車橋振動實驗數據,建立了具有11個自由度的3軸重型車輛和橋梁上部結構的車橋耦合系統有限元模型。對橋梁模型進行了模態分析,并對車橋耦合振動的三個工況進行了數值模擬,所得結果與實測數據一致,驗證了有限元模型的合理性、可靠性。進一步利用該有限元模型進行車橋耦合體系參數分析,選擇了車輛運行速度,橋梁橋臺搭板的沉降以及橋梁跨度進行橋梁沖擊系數的參數敏感性研究。參數分析結果表明:沖擊系數整體上隨車輛速度的增大而增大,同時有局部極大值的存在;沖擊系數隨橋臺搭板沉降量的增大而增大;沖擊系數呈現隨橋梁跨度的增大而增大的趨勢。

車-橋耦合振動,沖擊系數,LS-DYNA,重型車輛

1 引 言

近年來,隨著公路貨運重型大型車輛的高速發展,車橋耦合振動的相關研究成為研究的熱點,尤其交通運輸中橋梁損毀事件屢見報端,因而正確評估重型車輛對橋梁的動態性能影響顯得尤為重要。橋梁動態響應的確定是根據車輛靜力作用乘以沖擊系數得出,按照我國《公路橋涵設計通用規范》[1]的規定,沖擊系數僅僅考慮了橋梁的自振頻率的影響。實際上,橋梁和車輛是相互作用的耦合系統,影響橋梁動態響應的因素基本分為三類:第一類是橋梁的固有特性,比如自振頻率、跨長[2]等;第二類為車輛的特性,包括振動頻率、阻尼和運行速度等;第三類包括其他影響橋梁響應的因素,比如橋面的不平整度、橋臺搭板沉降等。由于橋梁節點處兩側橋面板不均勻沉降或橋臺搭板的沉降會造成車輛激勵的顯著增大,其中針對橋臺搭板的沉降對橋梁的影響,美國橋梁規范[3]允許橋梁節點的設計采用0.75作為沖擊系數,而我國規范并未做出相應規定。

在理論研究方面,很多學者通過建立簡單的車橋模型求解車橋耦合振動方程的解析解或者數值解析解。比如橋梁模擬成簡支梁,車輛以單軸質量彈簧系統代表[4,5]。盡管這種簡單模型提供了理論參考價值,但由于過于簡化無法應用于實際工程。有些學者通過編程解決車橋耦合運動方程,可以解決相對復雜的問題,但由于編程的復雜性使其不易推廣。比如侯永姣等[6]采用擬力法處理車輛模型求解車輛耦合振動方程;楊建榮等[7]采用單軸車輛模型研究了車橋耦合系統的迭代解法。隨著計算機硬件和軟件的快速發展,基于有限元方法的數值模擬軟件已被廣泛采用。值得一提的是,超強計算機的應用使精細化的車輛和橋梁模型的建立成為可能,提供了精確模擬車橋耦合系統的方法,但同時也增加了計算成本和時間。實驗手段因其可靠性并能夠驗證有限元模型而成為一種重要研究方法,但由于實驗代價較高,不可能進行詳盡的參數研究。

基于上述原因,本文對車橋耦合振動系統以有限元軟件LS-DYNA為平臺建立了11個自由度的3軸重型車輛和橋梁上部結構模型。車橋耦合模型規避了模型的復雜性并兼顧了計算的可靠性和計算效率。通過與實驗數據對比驗證了有限元模型的合理性、可靠性。選擇了車輛運行速度,橋臺搭板沉降量和橋面不平整度為參數,研究了橋梁沖擊系數的敏感性。

2 車橋耦合振動方程

基于有限元方法的車橋耦合振動基本方程可以表示成矩陣形式[8],橋梁的控制方程為[Mb]{d··

b}+[Cb]{d·b}+[Kb]{db}={Fb}(1)式中,[Mb],[Cb],[Kb]分別為橋梁的質量、阻尼和剛度矩陣;{db}為橋梁節點位移向量;{Fb}為作用在橋梁上的輪胎與橋梁的接觸力向量。

車輛的控制方程為

式中,[Mv],[Cv],[Kv]分別為車輛的質量、阻尼和剛度矩陣;{dv}為車輛節點位移向量;{FG}為車輛節點重力向量;{Fv}為作用在輪胎上的車輛與橋梁之間的接觸力向量。

利用位移連續條件和接觸力的關系,方程式(1)、式(2)可以寫成:

式中,Cbb,Cvb,Cbv,Kbb,Kbv,Kvb,Fbr是由于車橋接觸造成的增加項。

方程式(3)中的各項隨車輛的運行而變化,因而屬于時變系統。

3 車橋耦合系統有限元模型

重型車輛和橋梁模型根據文獻[9-11]的實驗數據建立。文獻中的實驗是以重型車輛的運行作為激勵測量橋梁動力特性和動態響應的實驗。車輛和橋梁的耦合振動通過有限元軟件LS-DYNA中的關鍵字RAIL[12]實現,其中,RAIL_TRAIN關鍵字定義了輪胎與軌道接觸的部分,接觸算法是通過罰函數實現。RAIL_TRACK定義了車輛運行的軌道,軌道由一系列梁單元組成嵌入橋梁單元,采用空材料(Null Material)以避免增加系統剛度。RAIL關鍵字同時允許施加路面的不平整度和輪胎的粗糙度。這種算法由英國Arup公司開發,為研究包括行人車輛等移動荷載造成建筑結構動態響應提供了方便。

3.1 車輛模型

場地實驗中采用的車輛是可以簡化成三軸的五軸重型卡車(圖1(a)),總重量為319 kN(32.5 t),其中前軸重量50 kN,中軸重量100 kN,后軸重量169 kN。車輛模型根據美國規范[3]中車輛荷載HS 20-44建立。這種簡化車輛模型兼顧了計算效率的同時體現了重型車輛多軸動力特性。如圖1(b)所示,車輛由拖車、掛車、懸掛系統和輪胎組成,其中拖車和掛車用球型鉸連接,并用剛體材料模擬,輪胎用單點質量單元模擬,懸掛系統和輪胎的剛度用彈簧阻尼系統模擬。考慮到拖車和掛車之間的連接關系并對車輛模型施加約束,使其具有11個自由度,包括6個輪胎豎直方向的位移y11,y21,y31,y12,y22,y32,拖車和掛車豎直方向的位移y1和y2,拖車和掛車的仰俯轉角θ1和θ2,車體傾覆轉角θ3。車輛的參數如表1所示,懸掛系統和輪胎的剛度滿足車輛的三個軸重,幾何尺寸及其他參數參考文獻[11]并做部分調整。通過模態分析得出車輛模型的整體豎向彈跳模態對應的頻率為3.42 Hz。

圖1 重型車輛模型Fig.1 Analyticalmodel of a heavy truck

表1 卡車模型參數Table 1 Parameters of the truck model

3.2 橋梁模型

實驗中的橋梁位于美國佛羅里達州,由護欄、板和大梁構成的混凝土跨線橋[11]。整個橋梁有三跨,每跨長度21.7 m,寬度14.15 m。橋面板厚度0.2 m,由6根間距2.4 m的I字形大梁支撐(圖2)。由于橋梁每跨單獨制作,并非連續性橋梁,因而只建立與公路相連的第一跨橋梁的上部結構模型。

圖2 橋梁上部結構橫截面Fig.2 Cross section of a bridge

圖3 橋梁的有限元計算模型Fig.3 Analytical finite elementmodel of the bridge

為了規避三維實體模型的計算代價,并確保數值模擬的可靠性,橋梁的護欄,板和大梁均采用分層殼單元進行模擬。分層殼單元的截面簡化計算方法參考文獻[13],護欄和大梁截面簡化過程遵循的原則是面積守恒和慣性矩守恒,并保持截面高度不變。分層殼單元的設置通過LS-DYNA軟件中的關鍵字PART_COMPOSITE實現,這個關鍵字允許對殼單元中的每一層分別設置不同材料參數,比如彈性模量和密度等。計算模型如圖3所示,沿橋梁橫向殼單元采用三種截面類型,分別代表護欄、板和大梁,其中大梁和板的單元類型交替布置。板對應的殼單元分成均勻6層,其中包含兩層鋼筋層,護欄位置的殼單元分成20層,其中包含位于下部的6層板單元,大梁位置的殼單元分成26層,其中包含位于上部的6層板單元。不同殼單元的每一層按照對應構件(大梁、護欄和板)的不同分別定義材料參數。由于實驗觀測到的動態變形較小,材料類型采用線彈性材料,材料常數按照混凝土核心試驗測得的數據采用[10]。

正確模擬橋梁模型的關鍵因素是橡膠支座模型的建立。本文采用LS-DYNA軟件中的離散梁單元(discrete beam element)模擬支座。這種單元類型可以方便的定義三個坐標軸方向的拉伸剛度和旋轉剛度,剛度參數的選擇參照文獻[14]。橋梁模型的邊界約束條件施加在支座上,并且沒有考慮支座與大梁和墩臺的摩擦和接觸。

3.3 有限元模型的驗證

橋梁模態分析是模型驗證的重要步驟。從模態分析結果中識別出三種模態和頻率與實驗結果接近,如圖4所示。由表2可以看出數值模擬結果與實驗結果吻合良好,第一階頻率的誤差為2.17%,表明了橋梁模型的可靠性。

圖4 橋梁模態和頻率的有限元結果和實驗結果[11]的對比Fig.4 Comparison of vibration modes between FE and experimental data

表2 有限元結果和實驗結果的頻率對比Table 2Comparison of natural frequency between finite element and testing results

與實驗動態時程結果的對此可以驗證車橋耦合模型的可靠性。本文選取了三個實驗工況進行模擬,第一個工況是車輛以48 km/h的速度運行在橋梁橫向的中間位置上;第二個工況是車輛以80 km/h的速度運行在橋梁橫向的中間位置上,并且在橋梁跨中橫向布置了截面為40 mm厚、400 mm寬的木板(木板的設置是為了增加車輛的激勵作用并模擬惡劣的路面粗糙度[9]);第三個工況是車輛以48 km/h的速度運行在橋梁橫向的遠離測量點的對面車道位置上。位移時程測量點在3號大梁的跨中位置(圖2),加速度時程的測量點在跨中護欄內側的橋面板上。

數值模擬的結果沒有進行濾波處理。第一個工況的模擬結果與實驗結果的對比如圖5和圖6所示。從圖5可以看出位移時程的曲線趨勢與實驗結果接近,實驗結果的最大位移為2.22 mm,數值模擬結果為2.26 mm,從圖6的加速度時程可以看出,盡管兩種結果的頻率不一致,但曲線的幅值和整體的趨勢吻合良好。圖7和圖8分別為第二個和第三個工況的位移時程結果對比,這兩個結果對比均顯示了良好的吻合性,從而驗證了車橋耦合模型的準確性。

圖5 在48 km/h速度下豎向動態位移的有限元結果和實驗結果的對比Fig.5 Vertical dynamic displacement comparison between experimental data and FE results at velocity of 48 km/h

圖6 在48 km/h速度下加速度的有限元結果和實驗結果的對比Fig.6 Acceleration comparison between experimental data and FE results at velocity of48 km/h

圖7 在80 km/h速度下橋梁豎向動態位移的有限元結果和實驗結果的對比Fig.7 Vertical dynamic displacement comparison between experimental data and FE results at velocity of80 km/h

圖8 運行在對面車道的車輛在48 km/h速度下橋梁豎向動態位移的有限元結果和實驗結果的對比Fig.8 Vertical dynamic displacement comparison between experimental data and FE results at velocity of 48 km/h with truck on the other lane

4 橋梁在車輛作用下的動態響應

沖擊系數是評估橋梁承受車輛作用的指標,公式如下:

式中,μ為沖擊系數;RD為橋梁在車輛作用下的最大動態響應;RS為橋梁在車輛作用下的最大靜態響應。

可以選擇位移、彎矩和應變作為響應進行計算,本文采用豎向位移響應計算沖擊系數。

車輛速度是影響橋梁響應的重要參數,本文按增量5 km/h選取車輛速度從30 km/h到120 km/h計算沖擊系數,如圖9所示。從圖9曲線可以看出沖擊系數的整體趨勢隨速度的增加而增加,同時存在一些局部極值,說明沖擊系數不是隨速度單調增加。這種由速度引起的局部極值點的現象發生在單質量、單軸及雙周車輛模型的車橋耦振動中,已在Fryba[5]的書中有所論述。但對三軸車輛模型引起的局部極值點的現象并未見其他文獻報道。

圖9 不同速度下的沖擊系數Fig.9 Impact factors at various velocities

由于橋梁使用時間的增加會導致橋臺搭板的沉降,或者橋梁節點處兩邊橋面的不均勻沉降造,從而造成成橋面凹陷,導致車輛的激勵增大,并增加了橋梁或節點處的動態響應。

參照文獻[15]中的橋臺搭板沉降計算模型(圖10),由于Δ2,Δ3的變化對橋梁的動態響應的影響比Δ1小[15],本文設定Δ2=0,Δ3=0,只考慮Δ1的變化。本文選擇了沉降量Δ1=5mm,10mm,20 mm,30 mm,并在車輛速度為60 km/h的條件下進行計算。如圖11所示,不同的沉降量對應的橋梁跨中節點處的位移時程表明:隨Δ1的增加,最大位移顯著增大。不同沉降量對應的最大位移分別為-2.23 mm,-2.36 mm,-2.72 mm和-2.89 mm,相應的沖擊系數分別為0.21,0.28,0.48和0.57。車輛第一個軸的軸力時程如圖12所示,表明沉降量的增大導致車輛軸力的增大,從而造成沖擊系數的增大。

圖10 橋臺搭板沉降計算模型[15]Fig.10 Model of approach slab deformation[15]

文獻[15]表明短跨橋梁對沉降量的變化比較敏感,為了考察沉降量對不同跨度的橋梁影響,通過改變現有橋梁模型的跨度[11,16]和大梁的截面高度[16],分別建立了跨度為16 m,26 m,30 m和34 m的橋梁模型。包括跨度為21.7 m的實驗橋梁模型的5個跨度的橋梁對應的第一個豎向振動模態的頻率分別為8.60 Hz,5.87 Hz,4.81 Hz,3.40 Hz,2.88 Hz,車輛的整體豎向振動頻率為3.42 Hz,設r為車輛頻率和橋梁頻率的比,則r=0.4,0.58,0.71,1.01和1.18。

圖11 在60 km/h速度下不同橋臺搭板沉降量的橋梁跨中豎向動態位移Fig.11 Dynamic vertical displacements atmidpoint of bridge at various approach slab settlements under 60 km/h traveling speed

圖12 在60 km/h速度下不同橋臺搭板沉降量的車輛第一軸的軸力時程Fig.12 Time histories of the first axle forces at various approach slab settlements under 60 km/h traveling speed

圖13 在120 km/h速度下,橋臺搭板沉降量為20 mm時不同跨度橋梁的跨中動態位移Fig.13 Dynamic vertical displacements atmidpoint of bridge for various span lengths with 20 mm approach slab settlement under 120 km/h traveling speed

為了便于比較,本文采用了無量綱時間參數,即

式中,τ為無量綱時間參數;t為車輛運行時間;L為橋梁的跨長;v為車輛運行速度。

在車輛運行速度為120 km/h的作用下,不同跨度橋梁在沉降量Δ1=20 mm的跨中位移時程如圖13所示,不同的沉降量對應的沖擊系數如表3所示。圖13中曲線表明,隨著橋梁跨度的增大,即橋梁自振頻率的減小,橋梁的動態位移隨之增大。從表3可以明顯看出,相同的Δ1值對應的沖擊系數在r=1.01時最大,此時橋梁和車輛耦合共振。除去共振的因素,相同的Δ1值對應的沖擊系數基本隨頻率的減小而增大,表明橋梁柔度的增加會導致動態響應的增加。我國橋梁規范[1]中規定跨度5 m≤L<20 m屬于小橋,20 m<L≤40 m屬于中橋,本文的算例中16 m的橋梁對應的沖擊系數沒有呈現出對沉降量敏感的趨勢。

表3 橋臺搭板沉降量對不同跨度橋梁的沖擊系數的影響Table 3 The effect of settlem ent of app roach slab on im pact factor of bridge of various span length

5 結 語

本文基于車橋耦合振動的場地實驗在LSDYNA的平臺上建立了車輛和橋梁的有限元模型,通過與實驗的模態頻率和時程數據的比較,驗證了有限元模型的合理性、可靠性,并對車輛運行速度、橋臺搭板沉降量和車橋頻率比對沖擊系數的影響進行了研究。

參數研究結果表明,車輛運行速度的增加整體上會導致橋梁沖擊系數的增大,同時在某些速度下沖擊系數存在局部極大值,表明橋梁的沖擊系數并非隨速度單調增加;橋臺搭板沉降量對橋梁的沖擊系數影響較大,是由于橋臺搭板的增加導致車輛軸力的增加造成;相同沉降量的條件下,沖擊系數呈現隨橋梁跨度增加而增大的趨勢,在車橋頻率比為1附近的共振點沖擊系數顯著增加。

[1] 中華人民共和國交通部.JTG D60—2004公路橋涵設計通用規范[S].北京:人民交通出版社,2004.Ministry of communications of the people′s Republic of China.JTG D60—2004 General code for design of highway bridges and culverts[S].Beijing:China Communications Press,2004.(in Chinese)

[2] 呂佳,吳定俊,李奇.簡支梁橋位移與內力動力系數差異研究[J].結構工程師,2012,28(4):78-83.Lu Jia,Wu Dingjun,Li Qi.Study on difference between displacement and internal force dynamic impact factor of simply-supported bridges[J].Structural Engineers,2012,28(4):78-83.(in Chinese)

[3] American Association of State Highway and Transportation Officials.LRFD bridge design specification[S].Washington,D.C.,2004.

[4] Olsson M.Finite element,modal co-ordinate analysis of structures subjected to moving loads[J].Journal of Sound and Vibration,1985,99(1):1-12.

[5] Fryba L.Vibration of Solids and Structures under Moving Loads[M].3rd Edition.Thomas Telford,London,1999.

[6] 侯永姣,吳定奇,李奇.車橋耦合振動分析中基于虛擬力的車輛建模方法[J].結構工程師,2012,28(3):55-59.Hou Yongjiao,Wu Dingjun,LiQi.Dynamic analysis method for vehicles based on force analogymethod in coupled vehicle-bridge systems[J].Structural Engineers,2012,28(3):55-59.(in Chinese)

[7] 楊建榮,劉章軍.車-橋耦合系統迭代解法研究[J].結構工程師,2007,23(4):41-44.Yang Jianrong,Lu Zhangjun.Iterative procedure for dynamic analysis of vehicle-bridge interaction[J].Structural Engineers,2007,23(4):41-44.(in Chinese)

[8] Deng L,Cai C S.Development of dynamic impact factor for performance evaluation of existing multigirder concrete bridges[J].Engineering Structures,2010,32:21-31.

[9] Kwasniewski L,Wekezer J,Roufa G,et al.Experimental evaluation of dynamic effects for a selected highway bridge[J].Journal of Performance of Constructed Facilities,2006,20(3):253-260.

[10] Kwasniewski L,Li H,Wekezer J,et al.Finite-element analysis of vehicle-bridge interaction[J].Finite Elements in Analysis and Design,2006,42(11):950-959.

[11] Li H,Wekezer J,Kwasniewski L.Dynamic response of highway bridge subjected tomoving vehicles[J].Journal of Bridge Engineering,2008,13(5):439-448.

[12] Hallquist JO.LS-DYNA Keyword User′s Manual,Livermore software technology corporation,Livermore,California,2013.

[13] Kwasniewski L.Nonlinear dynamic simulations of progressive collapse for amultistory building[J].Engineering Structures.2010,32:1223-1235.

[14] Cook R A,Allen D T,Ansley M H,et al.Stiffness evaluation of neoprene bearing pads under long term loads[R].University of Florida,2009,55-74.

[15] Shi X,Cai C S,Chen S,et al.Vehicle induced dynamic behavior of short-span slab bridges considering effect of approach slab condition[J].Journal of Bridge Engineering,2008,13(1):83-92.

[16] Liu C,Huang D,Wang Ton-Lo,et al.Analytical dynamic study based on correlated road roughness[J].Computers and Structures,2002,80:1639-1650.

Numerical Simulation of Bridge-and-heavy-vehicle Coupling Vibration for the Im pact Factor Study

WANG Juan*QIAN Jiang
(Research Institute of Structural Engineering and Disaster Reduction,Tongji University,Shanghai200092,China)

In order to efficiently evaluate the dynamic response of bridge under passage of vehicles,the superstructure of a bridge coupled with an heavy vehiclemodel of11 degrees of freedom with three axleswas developed by using LS-DYNA software based on available field tests.To validate reliability and efficiency of the finite element(FE)models and the corresponding algorithm,analyseswere carried out and three testing cases were selected for comparison.Good correlation between FE and experimental results confirms that the FEmodel is reliable for parametric studies.Three parameters includingmoving velocity of vehicle,the settlement of approach slab and the span length of the bridge,which has influence on the impact factor,were investigated.Numerical results show that impact factor increaseswith velocity whilemay reach its localmaximum somehow.It goes up with the increase of approach slab settlement.It has a tendency to increase with the increase of bridge span length.

dynamic response,LS-DYNA,finite elementmodel,impact factor,parameter study

2013-08-07

科技部國家重點實驗室基金項目(SLDRCE10-B-07)*聯系作者,Email:09 jwang@tongji.edu.cn

猜你喜歡
有限元橋梁模型
一半模型
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
手拉手 共搭愛的橋梁
句子也需要橋梁
高性能砼在橋梁中的應用
3D打印中的模型分割與打包
磨削淬硬殘余應力的有限元分析
基于SolidWorks的吸嘴支撐臂有限元分析
箱形孔軋制的有限元模擬
上海金屬(2013年4期)2013-12-20 07:57:18
主站蜘蛛池模板: 91在线日韩在线播放| AV不卡无码免费一区二区三区| 久久精品电影| 华人在线亚洲欧美精品| 伊人五月丁香综合AⅤ| 亚洲欧洲自拍拍偷午夜色| 亚洲精品中文字幕午夜| 精品丝袜美腿国产一区| 51国产偷自视频区视频手机观看 | www.狠狠| 国产精品理论片| 国产精品视频免费网站| 自偷自拍三级全三级视频| 在线不卡免费视频| 欧美日韩国产综合视频在线观看| 一本视频精品中文字幕| 久综合日韩| 欧美69视频在线| 人妻精品久久无码区| 久久精品无码一区二区日韩免费| 在线观看精品国产入口| 国产精品无码久久久久AV| 国产一区二区三区在线观看免费| 亚洲综合九九| 尤物国产在线| 欧美日本不卡| 久久青草精品一区二区三区| 999精品色在线观看| 日韩精品欧美国产在线| aa级毛片毛片免费观看久| 国产成人综合网在线观看| 黄色三级毛片网站| 六月婷婷综合| 美美女高清毛片视频免费观看| 国产麻豆aⅴ精品无码| 午夜三级在线| 久久亚洲综合伊人| 狠狠做深爱婷婷久久一区| 欧美精品二区| 精品久久人人爽人人玩人人妻| 天堂网亚洲系列亚洲系列| 国产在线97| 亚洲综合色区在线播放2019| 又爽又大又黄a级毛片在线视频| 99青青青精品视频在线| 久久五月天综合| 日韩毛片免费观看| 欧美翘臀一区二区三区| 国产日韩精品欧美一区喷| 亚洲男人在线天堂| 毛片免费网址| 男女猛烈无遮挡午夜视频| 国产在线观看成人91 | 亚洲天堂精品视频| 中文字幕一区二区人妻电影| 国产第八页| 午夜视频www| 亚洲区视频在线观看| 国产精品尤物在线| 综合色88| 亚洲自偷自拍另类小说| 亚洲日韩国产精品无码专区| 欧美黄网站免费观看| 日韩高清在线观看不卡一区二区| 国产精品福利尤物youwu | 中文字幕丝袜一区二区| 国产三区二区| 国产视频自拍一区| 中文纯内无码H| 国产精鲁鲁网在线视频| 在线观看国产小视频| 国产91av在线| 天堂久久久久久中文字幕| 日本国产在线| 国产精品网拍在线| 四虎精品国产AV二区| 少妇精品网站| 亚洲aaa视频| 国产区成人精品视频| 日韩无码黄色网站| 欧美不卡二区| 内射人妻无码色AV天堂|