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

新型雙足模型下結構動力響應分析

2016-05-20 02:26:35王益鶴
振動與沖擊 2016年7期
關鍵詞:振動

王益鶴, 楊 娜

(北京交通大學 土木建筑工程學院,北京 100044)

?

新型雙足模型下結構動力響應分析

王益鶴, 楊娜

(北京交通大學 土木建筑工程學院,北京100044)

摘要:基于考慮人體動力特性的新型雙足模型,采用Lagrange方程建立了步行過程中人與結構相互作用系統運動方程。研究了人與結構各自的動力響應及其之間相互作用,并通過對比說明了新型雙足模型反映的人體動力特性對結構振動的影響。結果表明:人與結構相互作用增大結構動力響應;行人經過結構各階模態振型波峰時,結構各階模態頻率和阻尼比分別達到最小值和最大值。

關鍵詞:振動;雙足模型;人與結構相互作用;動力特性

在土木工程領域,人行激勵下的結構振動問題尤其是人與結構相互作用問題,越來越受到人們的關注[1-2]。如設計不當,結構會在人群行走、跑步等運動下出現較大的振動響應,導致其舒適性和安全性下降。

眾多學者建立了各種數學模型描述人行激勵。其中最常見的是用傅里葉級數表達為靜載荷與幾個簡諧動載荷之和的確定性模型[3-4]。然而這種方式僅僅把人看作移動載荷,忽略了人作為獨立的動力體系這一特性,未考慮人與結構相互作用,預測具有輕質、低頻等特點的柔性結構振動時,不僅結果有所偏差,且無法研究人體對結構振動的反應和舒適情況。

考慮人體動力特性的生物模型開始應用于人與結構系統。人體被簡化為不同自由度的剛度、質量和阻尼體系(Spring Mass Damper,SMD)[5-7]。進一步研究表明[8-9],雙自由度(2 Degree of Freedom,2-DOF)的SMD模型較之單自由度(Single Degree of Freedom,SDOF)SMD模型,能更好的描述人體振動特性。而自由度更高的模型則用來研究振動在人體不同部位的精細傳遞過程,不適用于人與結構相互作用[10]。Kim等[11]也采用ISO 5982[12]的2-DOF人體模型分析人行橋的豎向振動問題。然而SMD模型僅能表示人體的豎向運動,無法反映人體步態尤其是行人步行中雙足支撐階段時的特性。故其不能很好的理解結構振動對人體動載荷的改變。Geyer等[13]提出用雙足模型模擬人體的步行運動,后來一些研究人員[14-15]也對雙足模型進行了分析,但均將人身體表示為質心,忽略了2-DOF人體模型所精確描述的人體動力特性,不能全面考慮人體對結構振動的影響。

本文綜合考慮ISO 5982的2-DOF模型描述的人體動力特性及雙足模型的雙足步行特點,建立了新型雙足步行模型。分別對人與結構相互作用系統中的人和結構進行了動力分析,并研究了新型雙足步行模型中人體動力特性對結構動力響應的影響。

1人與結構相互作用系統運動方程

新型雙足模型采用ISO 5982模型表示身體,人體質量等效為質量m1和m2,身體剛度為k1,阻尼為c1,長度為lb;腿部由兩個無質量的等長彈簧和阻尼表示的體系構成,其彈簧剛度為k1eg,阻尼為cleg,長度為l0(見圖1)。

圖1 新型雙足模型步行原理圖Fig.1 Schematic of the new bipedal model

采用Lagrange方程建立人與結構系統運動方程,人行橋由Euler-Bernoulli簡支梁表示。系統模型(見圖2)。在雙足支撐階段,人與結構系統的動能和勢能分別表示為:

圖2 人與結構相互作用系統模型Fig.2 Modal of human-structure interaction system

(1)

(2)

式中:z1、z2和u分別為人體m1和m2的垂直位移和水平位移;Ll和Lt分別為支撐腿和隨動腿彈簧長度;y為梁的垂直位移,y″(x,t)為曲率;m為單位長度梁質量;E為彈性模量;g為重力加速度。

梁的垂直位移由振型疊加法可表示為:

(3)

式中:{φi(x),i=1,2,…,n}為梁的振型;{Yi(t),i=1,2,…,n}為廣義坐標。

人與結構系統虛功的變分δW可表示為:

F1δz1+F2δz2+F3δu

(4)

式中:cl、ct和cs分別為支撐腿、隨動腿和結構的阻尼;vl和vt為兩腿的軸向速度;δ(ΔLl)和δ(ΔLt) 為兩腿的虛變形;{δYi,i=1,2,…,n}為梁的虛位移;δy″為梁曲率的變分;Q1,Q2,…,Qn,F1,F2,F3為對應廣義坐標{Y1,Y2,…,Yn,z1,z2,u}的廣義力。

支撐腿長度和隨動腿長度分別由其相對位置,可表示為:

(5)

支撐腿和隨動腿的軸向速度,可由m2與前后落腳點位置的相對豎向速度和及其水平速度合成得出,

(6)

式中:θl,θt分別為支撐腿和隨動腿與地面的夾角。

同理,支撐腿和隨動腿的位移變分,也可根據相對位移關系合成得出,

將式(3)和式(5)分別代入式(1)和式(2),可得:

(8)

(9)

將式(3)、式(6)和式(7)分別代入到式(4)中,虛功變分δW可表示為,

(10)

對比式 (4)和式 (10),廣義力表示為:

(11)

系統Lagrange方程可以表示為:

(12)

把動能、勢能和廣義力的表達式代入Lagrange方程,人與結構系統的運動方程表示為:

(13)

(14a)

(14b)

(14c)

(14d)

式中:

ci,i=2ξiwiMi+clsin2θlφi,i(N)+ctsin2θtφi,i(N-1),(i=1,2,…,n)

ci,j=clsin2θlφi,j(N)+ctsin2θtφi,j(N-1),(i≠j≤n)

ci,n+2=-clsin2θlφi(N)-ctsin2θtφi(N-1),(i=1,2,…,n)

ci,n+3=clsinθlcosθlφi(N)-ctsinθtcosθtφi(N-1),(i=1,2,…,n)

cn+2,i=-clsin2θlφi(N)-ctsin2θtφi(N-1),(i=1,2,…,n)

cn+3,i=clsinθlcosθlφi(N)-ctsinθtcosθtφi(N-1),(i=1,2,…,n)

φi,j(N)=φi(N)φj(N),φi(N)=

(15)

上述分析過程針對雙足支撐階段。而單足支撐階段時,隨動腿不與地面接觸,將所有與隨動腿相關的量設置為零即可。

由于雙足模型中存在阻尼耗能,故需提供外在機制進行控制。否則系統運動方程無法保持平衡,雙足模型經若干步后將失去穩定[13-15]。本文采用水平外力做功的方式維持系統能量恒定。水平控制力表示為如下形式:

(16)

式中:E0為人體初始能量;E(t)為時刻t時人體能量;Δu(t)為時刻t時人體水平位移增量。

人體總能量E(t)由動能、彈性勢能和重力勢能構成,可以表示為:

(17)

則式(14)中載荷向量變為:

F(t)=[0,0,0,…,k1lb-m1g,

(18)

由于系統中阻尼矩陣和剛度矩陣時變且含有未知量,故需用迭代的方法確定位移增量及速度增量。本文選取Newton-Raphson非線性方法進行迭代求解。

2人與結構相互作用分析

為了研究新型雙足模型下人與結構的相互作用,本文對上述推導過程進行Matlab編程,詳細分析整個系統的響應及動力特性。并將基于新型雙足模型的計算結果與基于時域力模型的計算結果進行比較。

2.1人與結構系統振動分析

對于雙足模型,腿剛度kleg、沖擊角θ0和步行初始速度v,是影響人體運動的關鍵參數。Geyer[13]首先給出腿剛度kleg和沖擊角θ0的參數區域圖,并得出雙足模型步行初始速度v的適用范圍為0.8~1.6 m·s-1。Whittington[14]列出了腿剛度kleg<50 kN·m-1,沖擊角θ0在50°~80°內的參數區域圖。Kim[16]將腿剛度kleg取值為14~28 kN·m-1,且認為阻尼比ξ在3%~8%范圍內即可描述人體步行。

新型雙足模型選取以下人體特性:身體部位參考ISO 5982模型[12],人體質量mh=75 kg(其中m1=62 kg,m2=13 kg),身體剛度k1=62 kN·m-1,身體阻尼c1=14.6 kN·s·m-1,身體長度取lb=0.4 m;腿部參數根據已有雙足模型[13-16],取腿剛度kleg=20 kN·m-1,阻尼比ξ=8%,腿長取為l0=1 m。步行初始速度為1.0 m/s,在運動過程中假定沖擊角θ0=69°為常數。

為闡述人與結構相互作用,假定兩個簡支人行橋。①算例一:剛度較大人行橋,11.0 m×1.25 m×0.35 m。材料特性為:單位長度質量m=1.364×103kg·m-1,剛度EI=1.64×108N·m2,模態阻尼比為ξ=0.3%。空載橋的前兩階頻率為:f1=4.50 Hz,f2=18.01 Hz。②算例二:柔度較大人行橋,9.0 m×0.8 m×0.1 2m。材料特性為:密度為2 400 kg·m-3,彈性模量E為3×1010N·m-2,模態阻尼比為ξ=0.3%。空載橋的前兩階頻率為:f1=2.375 Hz ,f2=9.50 Hz。

眾多研究者根據人體表觀質量幅頻特性研究人體模型[17]。圖3所示為新型雙足模型的人體表觀質量幅頻特性曲線。由圖3可知,新型雙足模型的人體表觀質量幅頻特性曲線與ISO 5982模型相近。由于ISO5982模型是由人體振動響應幅頻特性所得出,故新型雙足模型在模擬人體步行運動時,也可反映人體自身動力特性。

圖3 人體表觀質量幅頻特性Fig.3 Normalised apparent masses response functions

圖4所示為根據選取參數所得的行人單步步行載荷。新型雙足模型能夠捕捉人體重心變化產生的加速度及模擬兩腿交替行進,從而可得步行載荷的雙峰值特性。模擬得到的步行載荷特征與眾多的研究者試驗結果一致,可知新型雙足模型可較好的模擬人體的步行特性。

圖4 單步豎向地面反力Fig.4 Vertical ground reaction force with one foot

2.1.1算例一

首先對剛度較大的人行橋進行分析。圖5和圖6為行人經過剛度較大人行橋時跨中的位移和加速度時程。虛線為采用圖4中地面反力時程得到的人行橋跨中動力響應。由圖4可知,兩種方式所得結果差異較小。即結構剛度較大,振動響應較小,人與結構的相互作用較小。

圖5 人行橋的跨中位移Fig.5 Displacement at midspan of bridge

圖6 人行橋的跨中加速度Fig.6 Acceleration at midspan of bridge

圖7為行人步行通過結構時的豎向地面反力。由圖7可知,人與結構的相互作用較小,豎向地面反力時程平穩,結構振動對步行載荷的影響較小。

圖7 步行中產生的豎向地面反力Fig.7 Vertical ground reaction force in walking

由以上分析可知,結構剛度較大時,人與結構的相互作用較小。

2.1.2算例二

為分析人與結構相互作用,以下對柔度較大人行橋進行詳細分析。

圖8和圖9為行人經過柔度較大人行橋時跨中的位移和加速度響應時程。虛線為采用圖4中地面反力時程得到的人行橋跨中動力響應。由圖可知,兩種方式所得計算結果差異明顯。新型雙足模型所得最大位移值和最大加速度值分別比時域力模型大24%和37%。差異表明人與結構相互作用較大。即結構較大的振動引起新型雙足模型中腿部彈簧較大的變形,而腿部離開接觸面時彈簧中儲藏的彈性勢能得以釋放,使得行人又明顯影響結構響應,從而人與結構相互作用顯著。

圖8 人行橋的跨中位移Fig.8 Displacement at midspan of bridge

圖9 人行橋的跨中加速度Fig.9 Acceleration at midspan of bridge

圖10和圖11給出了人體上下兩部分m1和m2的豎向位移和加速度時程。如圖所示,每一步中行人的動力響應均不同,且隨著行人向跨中移動而增大。即隨著行人向人行橋跨中移動,結構振動對人體的動力響應影響增加,人與結構相互作用加大。同時人體上下兩部分m1和m2的動力響應也不相同。m1和m2的位移幅度接近,但是m1的位移滯后于m2。而由m1和m2的加速度局部放大圖可知,m1的加速度值要大于m2。即表明人體上下兩部分m1和m2的運動狀態并不一致。運動中的新型雙足模型即可描述人與結構相互作用的影響,又可反映人體不同部位的相對振動情況。

圖10 m1和m2的豎向位移Fig.10 The displacement of m1 and m2

圖11 m1和m2的加速度Fig.11 The acceleration of m1 and m2

在人與結構相互作用下,為研究結構振動對步行特性的影響,得出了行人步行通過結構時的豎向地面反力和行人步頻(見圖12和圖13)。由圖12可知,行人每一步產生的豎向地面反力均不同,隨著行人接近跨中,結構振動響應明顯,豎向地面反力將越來越大。結構振動對步行載荷產生較大影響。圖13為行人步行經過人行橋與在剛性地面上行走時的步頻對比。由圖13可知,隨著行人向跨中移動,步頻變化較大,結構振動對行人步頻產生影響。人的步行特性顯著發生改變,人與結構之間的相互作用較明顯。

圖12 步行中產生的豎向地面反力Fig.12 Vertical ground reaction force in walking

圖13 步行中的步頻Fig.13 Step frequencies in walking

圖14為人行橋跨中的加速度譜。虛線為采用圖4中豎向地面反力時程得到的結構跨中加速度譜。由于存在人與結構相互作用,新型雙足模型所得的加速度譜的峰值要大于時域力模型所得,即人與結構相互作用對結構振動響應影響顯著。從圖14可知,加速度譜的峰值主要位于行人步頻及其整數倍附近。由于人與結構相互作用下,結構振動使行人步頻產生較大變化(見圖13),故新型雙足模型模擬所得各峰值對應頻率均大于時域力模型所得。

圖14 人行橋跨中的加速度功率譜 Fig.14 Spectrum of acceleration at midspan of bridge

2.2結構動力特性

為研究行人對結構動力特性的影響,采用狀態空間法對人與結構系統進行瞬時模態的求解,得出算例二中行人作用下柔度較大人行橋前三階模態頻率和阻尼比(見圖15~圖17)。人與結構系統中人行橋動力特性同時受行人步態和行人位置影響,表現為各階模態頻率降低,阻尼比增加。而其整體變化趨勢主要與人的位置密切相關。對比圖18中人行橋前三階模態振型可知,當行人步行經過人行橋各模態振型的波峰時,各階模態頻率和阻尼比的變化也同時到達峰值。

以上研究表明,新型雙足模型下的人與結構系統更適合研究人與結構相互作用。而時域力模型則無法考慮人受結構振動影響產生的自身行為調整及對結構振動的反饋作用,即無法描述人與結構相互作用。

圖15 第一階模態Fig.15Thefirstmodal圖16 第二階模態Fig.16Thesecondmodal圖17 第三階模態Fig.17Thethirdmodal

圖18 人行橋的振型Fig.18 The mode shape of bridge

3人體動力特性對結構動力響應的影響

新型雙足模型基于ISO 5982模型和身體由質心表示的雙足模型(SDOF雙足模型),可反映人體不同部位的振動情況及其傳遞。為進一步分析全身振動特性對結構振動的影響,以下計算單人步行通過振動較大的梁。并將基于新型雙足模型所得的分析結果與SDOF雙足模型計算結果進行比較。

新型雙足模型選取人體特性與前相同,SDOF雙足模型選取和其一致的參數。考慮人行橋的基頻與人步頻接近,人行橋長為10.0 m,其余參數與算例二相同。空載橋的前兩階頻率為:f1=2.169 Hz,f2=8.677 Hz。

模擬得到行人在剛性地面上單步行走時的豎向地面反力時程(見圖19)。由圖19可知,在相同人體參數取值情況下,新型雙足模型和SDOF雙足模型所得的步行載荷特性接近,所得豎向地面反力略大于SDOF雙足模型。

圖19 單步豎向地面反力對比Fig.19 Comparison of vertical ground reaction force

圖20 人行橋的跨中位移對比Fig.20 Comparison of displacement at midspan of bridge

圖21 人行橋的跨中加速度對比Fig.21 Comparison of acceleration at midspan of bridge

圖20和圖21為兩種模型步行經過結構時跨中的位移和加速度。由圖可知,采用SDOF雙足模型所得結構振動響應要大于采用新型雙足模型所得振動響應。注意到圖19中,在相同人體參數取值情況下,新型雙足模型所得的豎向地面反力反而大于SDOF雙足模型。由此可得,結構為柔度較大結構時,振動響應較大,人體全身振動特性對結構振動響應將產生影響。

圖22為兩種雙足模型步行時橋跨中的加速度功率譜。由圖可知,由于人的步頻和結構基頻2.169 Hz接近,故將產生較大的振動響應,而SDOF雙足模型所得結構響應顯然大于考慮人體全身振動特性的新型雙足模型。即人體振動特性對結構振動影響不可忽略。考慮人體全身振動特性的新型雙足模型能更好的反映結構的振動響應。

圖22 人行橋跨中的加速度功率譜對比Fig.22 Comparison of spectrum of acceleration

4結論

分析了新型雙足模型下的人與結構相互作用系統,得到以下主要結論:

(1) 人與結構相互作用增大結構動力響應。新型雙足模型既可描述人體在振動中的動力特性,又可闡述結構振動對人的影響及人根據結構振動改變步態和動載荷,進一步增大結構動力響應的情況。

(2) 人步行經過人行橋時結構動力特性主要受人的步態和位置影響。人行橋各階頻率和阻尼比在行人位于各階模態振型波峰時,分別達到最小值和最大值。

參 考 文 獻

[ 1 ] ?ivanovi? S, Pavic A, Reynolds P. Vibration serviceability of footbridges under human-induced excitation: a literature review[J].Journal of Sound and Vibration, 2005, 279(1/2): 1-74.

[ 2 ] Sachse R, Pavic A, Reynolds P. Human-structure dynamic interaction in civil engineering dynamics: A literature review[J].Shock and Vibration Digest, 2003, 35(1): 3-18.

[ 3 ] Bachmann H, Ammann W. Vibration in structures-induced by man and machines [M]. Structural Engineering Documents No.3a. Zurich: International Association for Bridge and Structural Engineering (IABSE), 1987.

[ 4 ] Rainer J H, Pernica G, Allen D E. Dynamic loading and response of footbridges[J].Canadian Journal of Civil Engineering, 1988, 15(1): 66-71.

[ 5 ] Alexander N A. Theoretical treatment of crowd-structure interaction dynamics[J].Proceedings of the ICE-Structures and Buildings, 2006, 159(6): 329-338.

[ 6 ] 陳建英,方之楚. 人—結構相互作用動力學建模研究[J].振動與沖擊,2007, 26(6): 10-13.

CHEN Jian-ying, FANG Zhi-chu. Study on modeling of human occupants-structure dynamic interaction[J].Journal of Vibration and Shock, 2007, 26(6): 10-13.

[ 7 ] 秦敬偉,楊慶山,楊娜. 人體-結構系統靜態耦合的模態參數[J].振動與沖擊,2012, 31(15): 150-157.

QIN Jing-wei, YANG Qing-shan, YANG Na. Modal parameters of human body-structure static interaction systems[J].Journal of vibration and shock, 2012, 31(15): 150-157.

[ 8 ] Mansfield N J, Griffin M J. Non-linearities in apparent mass and transmissibility during exposure to whole-body vertical vibration[J]. Journal of Biomechanics,2000,33(8):933-941.

[ 9 ] Matsumoto Y, Griffin M J. Mathematical models for the apparent masses of standing subjects exposed to vertical whole-body vibration[J].Journal of Sound and Vibration,2003, 260(3): 431-451.

[10] Sachse R, Pavic A, Reynolds P. Human-structure dynamic interaction in civil engineering dynamics: a literature review[J].The Shock and Vibration Digest,2003,35(1):3-18.

[11] Kim S H, Cho K I, Choi M S, et al. Development of human body model for the dynamic analysis of footbridges under pedestrian induced excitation[J].Steel Structures, 2008, 8(4):333-345.

[12] ISO 5982. Vibration and shock-mechanical driving point impedance of the human body [S]. Geneva, Switzerland: ISO, 1981.

[13] Geyer H, Seyfarth A, Blickhan R. Compliant leg behaviour explains basic dynamics of walking and running[J].Proceedings of the Royal Society B: Biological Sciences, 2006, 273(1603): 2861-2867.

[14] Whittington B R, Thelen D G. A simple mass-spring model with roller feet can induce the ground reactions observed in human walking[J].Journal of Biomechanical Engineering, 2009, 131(1): 011013—1-8.

[15] Hong H, Kim S, Kim C, et al. Spring-like gait mechanics observed during walking in both young and older adults[J]. Journal of Biomechanics, 2013, 46(1): 77-82.

[16] Kim S, Park S. Leg stiffness increases with speed to modulate gait frequency and propulsion energy[J]. Journal of Biomechanics, 2011, 44(7): 1253-1258.

[17] Wei L, Griffin M J. Mathematical models for the apparent mass of the seated human body exposed to vertical vibration[J].Journal of Sound and Vibration, 1998, 212(5): 855-874.

Structural dynamic response analysis based on a new bipedal model

WANGYi-he,YANGNa(College of Civil Engineering, Beijing Jiaotong University, Beijing 100044, China)

Abstract:Based on a new bipedal model considering human dynamic characteristics, dynamic equations of human-structure interaction were established by means of Lagrange equation. Dynamic responses and interactions of human and structure were studied. The effects of human dynamic characteristics on structure vibration were also analyzed. The results showed that the human-structure interaction enhances structure vibration level the modal frequencies and damping radios of structure reach the minimum and the maximum, respectively when human walking passes the peak of each structural modal shape.

Key words:vibration; bipedal model; human-structure interaction; dynamic characteristics

中圖分類號:TU311.3

文獻標志碼:A

DOI:10.13465/j.cnki.jvs.2016.07.005

通信作者楊娜 女,博士,教授,博士生導師,1974年生

收稿日期:2015-04-27修改稿收到日期:2015-10-08

基金項目:國家自然科學基金資助項目(51422801;50938008;51278049);新世紀優秀人才獎勵計劃(NCET-11-0571);中央高校基本科研業務費專項資金資助(2012JBM007;2012YJS078)

第一作者 王益鶴 男,博士生,1984年生

猜你喜歡
振動
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
某調相機振動異常診斷分析與處理
大電機技術(2022年5期)2022-11-17 08:12:48
振動與頻率
天天愛科學(2020年6期)2020-09-10 07:22:44
This “Singing Highway”plays music
具非線性中立項的廣義Emden-Fowler微分方程的振動性
中立型Emden-Fowler微分方程的振動性
基于ANSYS的高速艇艉軸架軸系振動響應分析
船海工程(2015年4期)2016-01-05 15:53:26
主回路泵致聲振動分析
UF6振動激發態分子的振動-振動馳豫
計算物理(2014年2期)2014-03-11 17:01:44
帶有強迫項的高階差分方程解的振動性
主站蜘蛛池模板: 草逼视频国产| 久久久噜噜噜久久中文字幕色伊伊| 亚洲最大情网站在线观看| 久久综合五月婷婷| 国产精品专区第1页| 国产成在线观看免费视频| 午夜精品区| 伊人激情综合| 91精品国产一区| 亚洲国产综合自在线另类| 国产区在线观看视频| 九九热免费在线视频| 她的性爱视频| 91丝袜乱伦| 91蝌蚪视频在线观看| 伊人久综合| 久久亚洲精少妇毛片午夜无码| 欧美啪啪精品| 亚洲天堂在线免费| 国产精品黑色丝袜的老师| 一级毛片免费不卡在线视频| 国产精品国产三级国产专业不| 国产精品综合久久久| 在线看片中文字幕| 欧美成人亚洲综合精品欧美激情| 一本大道AV人久久综合| 婷婷午夜天| 色欲不卡无码一区二区| 狠狠色丁香婷婷综合| 亚洲天堂在线视频| 国产日韩丝袜一二三区| 中文字幕av无码不卡免费| 国产97视频在线观看| 无码中文字幕乱码免费2| 精品无码人妻一区二区| 在线观看亚洲人成网站| 久久a毛片| 色综合国产| 亚洲人精品亚洲人成在线| 区国产精品搜索视频| 啪啪啪亚洲无码| 亚洲欧洲日韩久久狠狠爱| 久久国产精品国产自线拍| 婷婷综合色| 国产亚洲欧美日本一二三本道| 久久精品一卡日本电影 | 本亚洲精品网站| 国产凹凸视频在线观看| 亚洲三级色| 亚洲精品成人片在线观看| 天堂av综合网| 国产欧美日韩在线一区| 毛片久久久| 亚洲第一网站男人都懂| 成人日韩精品| 99视频在线免费看| 91九色视频网| 不卡视频国产| 婷五月综合| 精品国产Av电影无码久久久| 一边摸一边做爽的视频17国产| 日本在线欧美在线| 色婷婷亚洲综合五月| 国产无人区一区二区三区| 成人在线亚洲| 久久免费看片| 中文字幕亚洲第一| a亚洲视频| 精品丝袜美腿国产一区| 国产成人永久免费视频| 91精品啪在线观看国产60岁| 喷潮白浆直流在线播放| 97国产成人无码精品久久久| 亚洲国产天堂久久九九九| 久久毛片网| 在线观看视频99| 毛片免费网址| 久久精品人人做人人爽电影蜜月 | 狼友视频一区二区三区| 亚洲狼网站狼狼鲁亚洲下载| 动漫精品啪啪一区二区三区| 国产精品三级专区|