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

一次強風作用下大跨度橋梁主梁非平穩抖振可靠性分析

2024-04-20 11:27:26葉澤毅阮偉東張新軍楊名冠
振動與沖擊 2024年7期
關鍵詞:風速有限元橋梁

孫 博, 葉澤毅, 阮偉東, 張新軍, 楊名冠

(1. 浙江工業大學 土木工程學院, 杭州 310023; 2. 中交公路規劃設計院有限公司, 北京 100088)

橋梁作為一種跨越山川河谷等天然障礙的結構形式,一直以來是交通系統中的重要節點和關鍵設施。我國服役橋梁眾多,且隨著我國現代化進程的加快,出現了許多大跨度橋梁建設工程[1]。橋梁對抵抗風災作用存在著天然劣勢和敏感性,隨著橋梁跨度的增長,其結構整體非線性問題也越發的突出,風致抖振效應對橋梁的影響也愈發的突出。運營期大跨度橋梁抖振作為一種持續性的結構動力響應,帶來疲勞問題并影響行車的舒適性,在強風作用下還有可能導致結構的強度破壞,產生巨大的安全隱患[2-3]。與此同時,橋梁抖振屬于動力分析范疇,結構的復雜性和風荷載的隨機時程特性都決定了大跨度橋梁抖振分析和評估問題存在大量不確定性。因此,為了保障大型橋梁結構在強風作用下的安全性,在概率分析框架下考慮動力分析問題的基本特征對大跨度橋梁的抖振可靠性進行正確的分析和評估具有十分重要的意義。

橋梁抖振可靠性分析受限于抖振響應分析方法和對風特性的認知程度,早期研究主要局在頻域內基于平穩風速時程評估橋梁的抖振可靠性。Pourzeymali等[4]以有限元法和頻域譜分析為基礎,結合一次二階矩(FOSM)和全分布法開展了橋梁疲勞可靠性分析。劉高等[5]采用有限元法和虛擬激勵法得到橋梁的抖振響應,并分別基于泊松分布和馬爾可夫過程的串聯失效模式,建立了系統抖振力可靠性分析過程。近年來,隨著橋梁抖振精細化分析手段的發展[6],橋梁抖振可靠性分析逐漸從頻域轉變到時域,結合精細化有限元在時域范圍分析結構的抖振可靠性。Jun[7]基于抖振時域分析,結合韋布爾疲勞壽命分布開展了疲勞可靠性評價。Ren等[8]結合通用有限元軟件和MATLAB進行大跨度橋梁抖振時域疲勞可靠性分析,并采用傳統的抖振頻域分析來驗證其結果。上述研究中,絕大多數采用風場為平穩高斯隨機過程進行可靠性分析,對于非平穩風場和非高斯響應下的抖振可靠性的研究相對較少。Hu等[9]在抖振時域分析的基礎上,考慮了結構和風荷載的非線性因素推導了結構的非高斯抖振響應,開展了動力可靠性的討論,但依舊是基于平穩風場。

隨著分析理論的發展和計算能力的提升,越來越多復雜工程的可靠性問題會直接采用更加精細化的計算分析方法來進行分析評估,以得到更加精確的結果。與此同時,由于大跨橋梁結構本身的非線性以及考慮強風情況下的非平穩風荷載效應,其抖振動力特性響應必然是非平穩和非高斯隨機過程。因此,本文在已有研究基礎上,基于非平穩風場模型得到非平穩風荷載效應,結合精細化的抖振時域分析方法,并考慮抖振響應的非平穩性和非高斯特征,構建了大跨度橋梁主梁非平穩抖振可靠性分析的方法流程。最后采用提出的方法對一次強風作用下某大跨度斜拉橋主梁的抖振可靠性進行了分析評估,證明了其有效性。

1 非平穩風場模擬

作用效應分析模擬是結構計算分析與安全分析的第一步,橋梁抖振可靠性分析首先需要進行考慮強風作用的非平穩風場模擬,包括風場模型構建和非平穩脈動風速兩個部分。

1.1 非平穩風場模型

自然風場沿著笛卡爾坐標系分解成順風向、豎向和橫向三個相互獨立的一維多變量非平穩隨機過程[10]。每一個獨立隨機過程由一個時變的平均風速分量與一個滿足非平穩特性的脈動風速分量組成。則三維非平穩風場模型可表示為

(1)

式中:Ux(t)、Uy(t)、Uz(t)分別為自然風U(t)沿著順風向x、豎向y、橫向z三個坐標軸方向的風速分量;U(t)為時變平均風速分量;u(t)、v(t)、w(t)為非平穩脈動風速分量。

(2)

1.2 非平穩脈動風速模擬

基于三維非平穩風場模型,Wang等[13]采用基于Deodatis雙索引頻率的諧波合成法(weighted amplitude wave superposition, WAWS)模擬出平穩脈動風速。諧波合成法主要利用功率譜密度函數矩陣的Cholesky分解和三角級數疊加來模擬平穩隨機過程,并引入快速傅里葉變換來減少模擬時間

(3)

式中:xj表示第j條平穩脈動風速;p=0,1,2,…,M×n-1,q=0,1,2,…,M-1,M=2×N,N表示采樣頻率,n為總模擬風速數量;Δt為時間間隔;ω為圓曲頻率,Δω=ωup/N為頻率間隔,ωup為截止頻率;Re(·)表示對括號的內容取實部;hjm(qΔt)為Bjm(lΔω)的快速傅里葉變換

(4)

(5)

本次樣本函數的空間相關性采用Davenport相關函數來模擬

(6)

(7)

根據進化譜理論,將平穩隨機過程的功率譜密度函數Sxx(ω)通過與經過參數優化的非均勻調制函數A(ω,t)相乘,可變為非平穩隨機過程的功率譜密度函數Syy(ω,t)

(8)

但在實際工程中,大多采用一個慢變均勻調制函數g(t)來代替非均勻調制函數A(ω,t)。根據Priestly建議的非平穩隨機過程滿足如下的R-S(Riemann-Stieltjes)積分和非平穩隨機過程的邊緣特性可得到[14]

(9)

式中:y(t)為非平穩隨機過程;x(t)為諧波合成法生成的平穩隨機過程;廣義變換Z(ω)為一正交增量過程。

2 非平穩風荷載計算

基于非平穩風場模擬結果,結合主梁截面的風特性參數,可進行非平穩風荷載模擬,包括靜風荷載、非平穩抖振力和自激力三個部分。

2.1 靜風力荷載

(10)

式中:FD、FL、FM分別表示靜風阻力、升力及扭矩;ρ為空氣密度;CD、CL、CM為梁截面的阻力、升力和升力矩系數;α0為有效風攻角;H為主梁特征高度;B為主梁的特征寬度。由于時變的平均風速變化相對較為緩慢,而在Davenport準定常理論中,靜力三分力系數是為了考慮流場流經過橋梁橫斷面后的變化,因此Chen[16]建議使用在平穩流場下獲得的三分力系數來計算非平穩靜風荷載的大小。

2.2 非平穩抖振力

非平穩抖振力是由非平穩脈動風速分量u(t)和w(t)引起。本文采用基于準定常氣動理論計的橋梁抖振力公式[17],并引入氣動導納公式來修正準定常氣動力模型計算的誤差,其修正后的沿主梁單位長度抖振力表達式為

(11)

(12)

(13)

2.3 氣動自激力

主梁風致振動改變了梁體周圍的氣動邊界條件,引起風場的變化,而風場的變化又會使主梁產生新的振動,這種新的振動激發力的力即為氣動自激力。為了實現自激力的時域化, 本文采用與自激力等價的12階單元氣動剛度矩陣K0與氣動阻尼矩陣C0[18],K0與C0通過對準定常氣動力模型進行雙變量泰勒展開得到

(14)

(15)

式中,ml=B/4為主梁旋轉平均半徑。

(16)

(17)

式中,L為主梁單元長度。

3 非平穩抖振動力可靠性分析

獲得了作用大跨度橋梁主梁上的風荷載作用效應后,基于精細化的有限元分析方法進行時域分析,可以獲得結構的抖振響應,考慮其非平穩性和非高斯特性對主梁進行動力可靠性分析。

3.1 抖振時域化分析

在大跨度非平穩抖振特性分析中,其動力求解控制方程可表示為

(18)

(19)

對式(18)的求解主要有頻域法和時域法,早期主要在頻域進行,但由于頻域法在進行橋梁抖振分析時,不能全面地考慮多種非線性因素,而時域法能有效的避免這些問題,加上計算機性能的提高以及各類具有強大功能有限元軟件的開發和應用,使得時域法已成為橋梁抖振分析的主要發展方向[19]。結合通用有限元軟件的橋梁抖振時域分析流程主要分為4步:① 依據橋梁設計資料和風洞試驗得到的主梁風特性參數,建立考慮自激力的大跨度橋梁空間有限元模型;② 結合《規范》和橋址處的各類風參數,模擬非平穩脈動風速時程;③ 采用顯示的靜風荷載和抖振力荷載公式,將風速時程轉換為力的時程;④ 將得到力的時程,施加到有限元模型上,并采用New-Raphsan與Newmark-β法進行瞬態分析,得到橋梁的位移響應。

需要注意的是,由于大跨度橋梁結構本身屬于典型的柔性體系,在通用有限元軟件中開展瞬態分析需要通過計入大變形及應力剛化效應來考慮結構的非線性特性。此外,在幾何非線性求解過程中,針對參數不合適而出現非線性求解困難的情況,需打開自動時間步長激活二分法,使得求解繼續。

3.2 非平穩隨機過程首超概率

大跨度橋梁結構在一次強風作用下的抖振響應屬于一個非平穩隨機過程,本文采用基于Poisson假定的首次超越概率來求解主梁的抖振動力可靠性。在大跨度斜拉橋主梁抖振可靠性分析問題中,邊界類型為雙側邊界,其抖振安全邊界,參考現行JTG/T 3365-01—2020《公路斜拉橋設計規范》[20]中的剛度失效邊界為

(20)

式中:R為橋梁抖振可靠性界限;L為橋梁主跨跨徑。

大跨度斜拉橋抖振響應中,靜風荷載會導致結構安全界限出現不對稱的情況,而泊松過程法能夠完美地適應各種對稱或者不對稱的雙側界限。此時抖振可靠性分析的安全邊界為[-b1,b2]

b1=R+S,b2=R-S

(21)

式中:b1為下側邊界;b2為上側邊界;S為斜拉橋在靜風荷載下的響應。

當結構響應為非平穩高斯隨機過程時,對于在非平穩風場模型下的雙側界限來說,基于泊松假定的時變動力可靠概率Pr(-b1,b2,T)可表示為[21]

(22)

3.3 非高斯抖振響應處理

在非平穩風荷載作用下,主梁的非平穩抖振響應為非高斯響應[22],不符合式(22)中的高斯隨機過程假定。針對這一問題,可以利用一個單調平移函數將非高斯響應轉換為更加符合標準高斯分布函數的響應,進而較為方便的求解抖振動力可靠概率。本文采用Winterstein修正模型[23]對非高斯抖振響應進行處理,其采用的多項式為

X=μX+κσX[Z+c3(Z2-1)+c4(Z3-3Z)]

(23)

式中:μX和σX為隨機過程X的均值和標準差;Z為標準高斯過程;c3、c4、κ為多項式系數,可通過X的偏態系數γ1和峰態系數γ2計算得到

(24)

(25)

(26)

則z-x轉換形式為

(27)

式中

(28)

將逼近高斯概率密度函數后的變量組z,進行逆標準化處理

G=zσX+μX

(29)

式中,G為我們所需要進行可靠度計算的新變量組。

4 算例分析

4.1 工程概況

某斜拉橋為鋼索塔鋼箱梁雙索面五跨斜拉橋,跨徑布置為63+255+648+255+63=1 284 m,主跨648 m。采用半漂浮結構體系,縱向設彈性約束,限制活荷載及風載作用下的鋼箱梁縱向漂移。主梁為帶風嘴的閉口箱梁斷面,梁高為3.2 m,主梁寬度為37.16 m。橋梁布置示意圖如圖1所示,主梁截面如圖2所示。

圖1 案例斜拉橋布置圖Fig.1 Cable-stayed bridge layout

圖2 案例主梁截面Fig.2 Girder section

根據橋梁建設相關資料[24],采用通用有限元軟件ANSYS建立大跨度斜拉橋空間有限元模型。鋼箱梁主梁和橋塔使用BEAM44梁單元,共594個梁單元,拉索使用LINK10桿單元,共168個桿單元,二期恒載采用MASS21模擬。圖3為ANSYS全橋有限元模型。

圖3 案例橋梁ANSYS有限元模型Fig.3 ANSYS finite element model of bridge

采用ANSYS有限元軟件自帶的模態分析模塊,進行動力特性分析,得到了成橋狀態的主要自振振型頻率計算數據,并與風洞試驗數據和實測自振頻率[25]進行對比。表1給出了三者的自振頻率及其偏差。圖4給出了相應自振頻率的振型圖。

表1 橋梁自振頻率檢驗Tab.1 Natural frequency verification of bridge

(a) 一階振型(正對稱豎彎)

(b) 三階振型(正對稱豎彎)

(d) 四階振型(反對稱豎彎)

(c) 五階振型(正對稱豎彎)

(f) 六階振型(反對稱豎彎)

(d) 七階振型(反對稱豎彎)

(h) 八階振型(反對稱豎彎)圖4 橋梁基頻振型圖Fig.4 Fundamental frequency mode under service state

4.2 荷載效應模擬結果

《規范》中采取600 s為風速的標準時距來計算風荷載效應,本文同樣采用時間長度為600 s作為一次強風歷經時長,并分別模擬在主梁處平均風速為40 m/s、50 m/s、60 m/s的非平穩抖振效應。

針對非平穩風場中的非平穩脈動風模擬,采用《規范》推薦的使用的Kaimal順風向譜Su(n)與Panofsky豎向譜Sw(n)

(30)

(31)

式中:n為脈動風頻率;f為莫寧坐標;u*為風的摩阻速度。沿主梁縱向共設置44個模擬風速點,時間間隔0.125 s,截止頻率3 πrad/s,頻率采樣點個數211。圖5給出了平均風速為40 m/s時,主梁跨中處的順風向與豎向平穩脈動風速的模擬及功率譜密度函數檢驗結果。

(b) 主跨跨中平穩豎向脈風速時程

(c) 平穩跨中順風向脈動風速模擬譜與目標譜對比

圖5 平穩脈動風速模擬示意圖Fig.5 Simulation of stationary turbulence

從圖5(a)、(b)可以看出來,考慮了sears函數的平穩脈動風速均值基本為零,其順風向模擬的結果主要在[-10,15]m/s范圍內波動,而橫向模擬結果在[-5,5]m/s范圍內波動,二者的平穩性均較好。由圖5(c)、(d)可知順風向和豎向平穩功率譜的模擬值和理論值在頻域內高度吻合。由圖5(e)、(f)可知主梁跨中順風向和豎向的相關函的模擬值和理論值有較好的吻合效果。

在平穩風速模擬基礎上引入慢變均勻調制函數g(t)將其轉換為非平穩脈動風速,采用的g(t)形式如下

(32)

式中,α、β為非平穩特性的調制參數。經過參數優化篩選,分別取值為α=300、β=40 000,則調制函數如圖6所示。調制函數的最大值為1,最小值為0.1,其一個周期剛好與本次模擬強風作用的時間相同,且與臺風經過時的能量變化規律類似,呈現出一個兩端小中間大的慢變均勻變化趨勢。

圖6 慢變均勻調制函數Fig.6 Uniform nodulating function

由圖7給出了利用式(9)最終得到的非平穩脈動風速時程模擬結果。通過均勻調制函數得到的非平穩順風向與豎向脈動風速時程整體均值為零,具有很明顯的時變趨勢,同時還有一定的時變方差性。與慢變均勻調制函數的變化趨勢相符,說明模擬的結果符合非平穩脈動風速的特性。

圖7 非平穩脈動風速時程圖Fig.7 Time history of simulation of non-stationary turbulence

基于風場模擬結果和風洞試驗得到的氣動力系數,可以進行風載模擬效應。圖8給出了部分模擬點的非平穩抖振力效應代表性時程圖(40 m/s),圖9為模擬采用的主梁氣動力系數圖。

(a) 非平穩抖振扭轉代表性模擬點

(b) 非平穩抖振阻力代表性模擬點

圖9 主梁氣動力系數Fig.9 Drag, lift, and moment coefficient curves

4.3 關鍵點動力位移響應結果

將得到的風荷載結果施加到有限元模型上,其中靜風荷載與非平穩抖振力采用外部文件的方式導入,而氣動自激力則在ANSYS中采用MATRIX27矩陣的形式輸入。利用ANSYS的瞬態動力學分析功能并考慮幾何非線性與剛度硬化效應,計算得到案例橋梁主梁時變抖振響應結果。

圖10給出了主梁跨中抖振位移時程,主梁的豎向位移變化幅度呈現出中間大兩邊小的趨勢,這與非平穩脈動風速時程圖及力的時程圖變化趨勢相近,且具有明顯的方差時變特性。圖11給出了抖振位移對應的標準化概率密度函數擬合與標準高斯分布概率密度函數的對比結果,可以看出兩者相差較大,說明直接采用原數據進行非平穩抖振動力可靠度計算會產生較大的誤差,故需要對數據進行高斯轉換。

圖10 非平穩抖振主梁跨中位移時程Fig.10 Time history of non-stationary buffeting displacement responses in the vertical mid-span girder

圖11 跨中豎向位移概率密度函數擬合圖Fig.11 The standard probability density function in the vertical mid-span girder

根據得到的抖振位移響應,圖12給出了采用Winterstein修正模型的高斯轉換結果。高斯轉換之后的主梁跨中位移時程相比于之前的位移時程圖其幅度有所降低。圖13給出了高斯轉換前后標準化概率密度函數對比結果,變量組在標準化高斯轉換之后,其偏態以及峰度相比之前更加貼近標準高斯分布的概率密度函數。由圖14(a)圖可知,響應數組在高斯轉換后,其變量組的均值未發生明顯的變化,從圖14(b)圖可知,在高斯轉換后的數據在跨中處更具有穩定性。

圖12 高斯轉換后的非平穩抖振主梁跨中位移時程Fig.12 Time history of non-stationary buffeting displacement responses in the vertical mid-span girder after Gaussian translation

圖13 高斯轉換前后標準化概率密度函數對比圖Fig.13 The standard probability density function in the middle of mid-span girder before and after Gaussian translation

(a) 響應標準差

(b) 響應導數標準差圖14 高斯轉換前后響應統計特性對比圖Fig.14 The Statistical characteristic of responses before and after Gaussian translation

在本次非平穩可靠性分析中,主要采用變量組的標準差與其導數的標準差計算可靠概率,因此由圖14(a)、(b)可知本次高斯轉換前后數據在標準差與其導數的標準差并未發生明顯的變化。

4.4 非平穩抖振時變動力可靠度結果

由于得到的高斯轉換后的動力特性響應具有很強的非平穩特性,故需要計算其時變方差。對主梁每個計算位置,采用10 s時間區間,每個區間80個抖振響應樣本進行計算,共得到60個響應方差樣本,圖15給出了主梁跨中的抖振響應時變方差的計算結果。

圖15 主梁跨中時變方差(40 m/s)Fig.15 Time-varying variance in the vertical of mid-span girder

此外,還需要對可靠度邊界進行處理。圖16為主梁在靜風荷載下的豎向位移響應,利用式(15),對主梁豎向剛度的動力可靠性界限進行處理,結果如圖17所示。

圖16 主梁豎向空氣靜力位移響應沿橋跨的分布Fig.16 Distribution of aerostatic force displacement responses on the vertical of mid-span girder

圖17 主梁豎向非平穩抖振可靠度界限沿橋跨分布Fig.17 Distribution of vertical non-stationary buffeting reliability limit on the girder

經過上述結果,通過首次穿越概率的原理計算其時變動力可靠度,圖18給出了不同平均風速時抖振動力可靠性結果沿主梁的分布。表2給出了南塔邊跨跨中(模擬點12)、主跨四分之一跨(模擬點34)、主跨跨中(模擬點45)、主跨四分之三跨(模擬點56)、北塔邊跨跨中(模擬點78)的具體的動力可靠性數值。結果表明斜拉橋主梁非平穩抖振動力可靠性隨一次強風平均風速的上升而下降。主梁跨中在風速為50 m/s的情況下,最先開始出現變化,并隨著風速的增大,其可靠性逐漸減小。對于主跨處四分之一跨(模擬點34)與四分之三跨(模擬點56),也在60 m/s風速下,出現了可靠性的變化。主梁跨中(模擬點45)的可靠性最低,為非平穩抖振的最不利位置。

圖18 不同平均風速時可靠性沿主梁分布圖(40 m/s)Fig.18 Distribution of different mean wind velocity reliability on the girder

表2 主梁節點模擬的非平穩抖振力可靠性Tab.2 Simulated girder points of non-stationary buffeting reliability

5 結 論

科學有效的分析強風作用下大跨度橋梁的抖振可靠性,對保障橋梁安全和道路交通體系的順利運營具有十分重要的意義。本文針對強風荷載作用的非平穩特性,基于精細化的抖振時域顯式分析過程,并考慮抖振相應的非平穩和非高斯特性,提出了一次強風作用下大跨度橋梁主梁非平穩抖振可靠性分析的方法流程,主要工作總結如下:

(1) 非平穩強風荷載效應的時域化模擬。采用諧波合成法模擬得到主梁的兩個獨立平穩脈動風速,結合進化譜理論將平穩脈動風速轉化為非平穩脈動風速。根據現有顯示靜風與非平穩抖振力荷載公式求解時變風荷載,并通過引入單元氣動剛度矩陣與氣動阻尼矩陣實現氣動自激力。

(2) 提出了考慮非平穩非高斯抖振相應的動力可靠性分析方法。對于精細化有限元分析得到的抖振響應,考慮響應方差的時變特性,采用基于Poisson假定的首次超越概率來求解非平穩抖振動力可靠性。對于響應的非高斯特性,則采用Winterstein修正模型對其進行高斯轉換。

(3) 應用所提出的方法流程對某大跨度斜拉橋在一次強風作用下的抖振可靠性進行了分析評估,結果表明:① 非平穩脈動風速模擬結果與均勻調制函數的變化趨勢相符,符合強風歷經時的能量變化規律;② 抖振響應具有明顯的非平穩和非高斯特性,采用高斯轉換處理響應結果并考慮方差的時變特性是十分必要的;③ 案例斜拉橋主梁非平穩抖振動力可靠性隨一次強風平均風速的上升而下降,主梁跨中可靠性最低,為非平穩抖振的最不利位置。

猜你喜歡
風速有限元橋梁
基于Kmeans-VMD-LSTM的短期風速預測
基于最優TS評分和頻率匹配的江蘇近海風速訂正
海洋通報(2020年5期)2021-01-14 09:26:54
手拉手 共搭愛的橋梁
句子也需要橋梁
高性能砼在橋梁中的應用
基于GARCH的短時風速預測方法
考慮風速分布與日非平穩性的風速數據預處理方法研究
磨削淬硬殘余應力的有限元分析
基于SolidWorks的吸嘴支撐臂有限元分析
箱形孔軋制的有限元模擬
上海金屬(2013年4期)2013-12-20 07:57:18
主站蜘蛛池模板: 国产成人精品一区二区三在线观看| 日本成人不卡视频| 久久国产精品影院| 国产成人精品午夜视频'| 综合亚洲网| 久久国产高清视频| 午夜国产大片免费观看| 欧美全免费aaaaaa特黄在线| 色综合中文字幕| 日韩中文字幕免费在线观看| 内射人妻无码色AV天堂| 国产精品专区第1页| 欧美成人h精品网站| jizz国产视频| 欧美在线中文字幕| 欧美综合中文字幕久久| 日韩人妻无码制服丝袜视频| 无码 在线 在线| 久久综合伊人 六十路| 日本伊人色综合网| 亚洲日韩在线满18点击进入| 人妻精品全国免费视频| 国国产a国产片免费麻豆| 国产精品播放| 2020精品极品国产色在线观看 | 2020亚洲精品无码| 国产一区在线视频观看| 亚洲动漫h| 国产精品主播| 九色综合视频网| 91www在线观看| 欧美成人午夜视频| 久久99蜜桃精品久久久久小说| 一级不卡毛片| 亚洲精品在线影院| 三级毛片在线播放| 亚洲第一区欧美国产综合| 亚洲人成网18禁| 人妻无码中文字幕第一区| 亚洲国产精品国自产拍A| 免费可以看的无遮挡av无码| 成人精品免费视频| 久久99精品久久久久久不卡| 国产精品自在线拍国产电影| 久久综合九九亚洲一区| 亚洲一区二区三区国产精华液| 免费一级成人毛片| 亚洲三级电影在线播放| 国内精品视频| 亚洲精品在线91| 国产尤物在线播放| 免费一级毛片不卡在线播放| 成年网址网站在线观看| 国产精品夜夜嗨视频免费视频| 国产中文在线亚洲精品官网| 日韩欧美一区在线观看| 97成人在线视频| 日韩精品一区二区深田咏美| a级免费视频| 精品无码国产自产野外拍在线| 国产96在线 | 日韩欧美中文字幕一本| 日本伊人色综合网| 欧美高清国产| 亚洲天堂久久| 亚洲IV视频免费在线光看| 日韩精品久久无码中文字幕色欲| 美女免费精品高清毛片在线视| 欧美一级爱操视频| 国产精品刺激对白在线| 波多野结衣在线一区二区| 国产成人精品高清不卡在线| 亚洲国产午夜精华无码福利| 国产精品第三页在线看| 日本亚洲最大的色成网站www| 亚洲第一成年免费网站| 美美女高清毛片视频免费观看| 亚洲精品无码久久久久苍井空| 国产成人高清亚洲一区久久| 1769国产精品免费视频| 国产在线精彩视频二区| 婷婷亚洲视频|