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

利用有限振動響應(yīng)快速反演時變荷載的方法與試驗

2022-06-29 09:57:18項子儒錢王蘋顧鎮(zhèn)媛
振動與沖擊 2022年12期
關(guān)鍵詞:橋梁方法模型

項子儒, 王 軼, 錢王蘋, 顧鎮(zhèn)媛

(1.南通大學(xué) 交通與土木工程學(xué)院,江蘇 南通 226019;2.暨南大學(xué) 重大工程災(zāi)害與控制教育部重點(diǎn)實驗室,廣州 510632)

2019年10月,江蘇省無錫312國道錫港路上跨橋路段高架橋因嚴(yán)重超載發(fā)生橋面?zhèn)确鹿剩斐闪司薮蟮娜藛T與經(jīng)濟(jì)損失[1],車輛超載對橋梁健康的威脅被再一次擺在聚光燈下。我國交通體系經(jīng)過飛速的發(fā)展和壯大,多數(shù)橋梁的使用壽命已進(jìn)入中齡,其承載能力不可避免地降低,超重荷載時刻挑戰(zhàn)著橋梁結(jié)構(gòu)的安全性。自國家提出“十三五”戰(zhàn)略以來,我國的橋梁建設(shè)已從“建設(shè)為主”逐漸向“建養(yǎng)并重”轉(zhuǎn)型[2]。移動荷載是影響橋梁使用壽命的主要因素之一,對其準(zhǔn)確的監(jiān)測與識別是提供“養(yǎng)橋”的依據(jù),更是預(yù)防災(zāi)難性事故的關(guān)鍵。

實際工程中,多數(shù)橋梁采用動態(tài)稱重設(shè)備采集車輛質(zhì)量,然而這種設(shè)備的安裝與養(yǎng)護(hù)費(fèi)用較高,且需要預(yù)安裝于調(diào)查路段上,具有很大的局限性。事實上,移動車輛對橋梁的真實荷載是自身車重與車-橋梁耦合力的疊加,僅將車重作為移動荷載會產(chǎn)生較大的誤差。因而針對移動荷載識別(moving force identification,MFI)的研究更傾向于利用動荷載作用下結(jié)構(gòu)響應(yīng),反演車輛與橋梁接觸面的相互作用力。自19世紀(jì)70年代起,MFI經(jīng)過快速的發(fā)展,已積累的不少的研究成果。目前較成熟的方法主要包括一類解析法、二類解析法、時域法和頻時域法[3]。后兩類的基本思想是通過計算實測的結(jié)構(gòu)響應(yīng)和提前獲取的系統(tǒng)脈沖響應(yīng)函數(shù)的反卷積計算荷載,也就是所謂的反卷積方法。反卷積法以其形式簡單、可靈活采用時域和頻域內(nèi)分析的特點(diǎn)應(yīng)用較廣泛。然該方法歸根結(jié)底是反問題的線性方程求解,無法避免求解的病態(tài)性:若測量噪聲大、初始條件不準(zhǔn)可引起系統(tǒng)矩陣(脈沖響應(yīng)矩陣)奇異性,導(dǎo)致求解失敗;若測量時間長、頻率高,則系統(tǒng)矩陣維度巨大,計算耗時明顯,難以實現(xiàn)動態(tài)監(jiān)控。

在約束反問題病態(tài)性方面,常用方法為正則化方法,即在原識別問題中引入合理的附加約束,從而達(dá)到改良問題病態(tài)性的目的。常用正則化有:Tikhonov正則化[4]、稀疏正則化[5-6]和基函數(shù)展開法[7]等。其中,稀疏正則化識別法是目前相關(guān)研究的前沿?zé)狳c(diǎn),其在應(yīng)用時會充分考慮荷載屬性在某種特定情況下的稀疏性。Bao等[8]考慮荷載的稀疏性分布,利用斜拉橋索力識別移動重載。Pan等[9]考慮荷載時程信號的多域稀疏性,結(jié)合冗余字典識別了簡支梁上移動荷載。應(yīng)用正則化法的優(yōu)勢在于能抑制噪聲的影響,獲取穩(wěn)定的識別結(jié)果;而其劣勢在于需要在求解過程中加入正則化流程,增加識別成本。另一個思路是利用函數(shù)逼近方法,構(gòu)造連續(xù)平滑的函數(shù)如三角函數(shù)、冪函數(shù)、小波函數(shù)等[10]代替單位脈沖函數(shù),從而削減奇異值。其中張青霞等[11]提出梁單元形函數(shù)的概念,利用有限元理論中的形函數(shù)概念逼近動態(tài)荷載曲線,從而把識別未知荷載的時間歷程轉(zhuǎn)變?yōu)樽R別數(shù)目有限的形函數(shù)系數(shù),在削減奇異值的同時降低了系統(tǒng)矩陣維度,提高反演效率;該研究利用數(shù)值仿真探討了形函數(shù)對三角荷載、沖擊荷載等周期性荷載的擬合與識別,均得到較好的結(jié)果。王蕾等[12]通過數(shù)值模擬,使用荷載形函數(shù)(load shape function,LSF)識別了大跨橋梁結(jié)構(gòu)上的移動荷載。作者將該方法結(jié)合于移動荷載與預(yù)應(yīng)力的協(xié)同識別[13]計算中,不僅得到較高的噪聲魯棒性,且成功將計算效率提升數(shù)倍。

然而該方法在建立LSF響應(yīng)矩陣時需計算橋梁與荷載各接觸點(diǎn)的脈沖響應(yīng),運(yùn)算成本巨大。本文在現(xiàn)有研究基礎(chǔ)上,定義新的構(gòu)造矩陣提取傳統(tǒng)LSF響應(yīng)矩陣中的有效元素,進(jìn)一步降低識別成本;其次利用數(shù)值模擬分析該改進(jìn)方法應(yīng)用于簡支梁橋MFI問題的優(yōu)勢;最終通過試驗?zāi)P万炞C本方法的識別準(zhǔn)確性、時效性以及對噪聲的魯棒性。

1 時變荷載識別

1.1 時域反卷積法

本文將橋梁結(jié)構(gòu)考慮為線性系統(tǒng)。對于受連續(xù)力激勵的系統(tǒng),如其初始狀態(tài)為零,其振動響應(yīng)可利用杜哈梅爾積分,由激振力與相應(yīng)的單位脈沖響應(yīng)函數(shù)積分求得

(1)

對式(1)進(jìn)行離散化處理,進(jìn)一步推導(dǎo)出多激勵作用下多測點(diǎn)響應(yīng)矩陣

(2)

(3)

1.2 構(gòu)造LSF響應(yīng)矩陣

引用LSF能有效抑制上述兩種問題的影響。該方法將激勵荷載的時程信號比擬為一跨“有限元梁”,并劃分為若干單元,利用每個單元的形函數(shù)擬合“梁”的位移,即荷載信號的幅值

(4)

將式(4)代入式(2),將未知外荷載的識別轉(zhuǎn)化為擬合系數(shù)的識別問題,

(5)

(6)

(7)

圖1 系統(tǒng)矩陣元素分布進(jìn)化Fig.1 The evolution of system matrix

2 移動荷載識別的數(shù)值仿真

本章采用數(shù)值模擬展示構(gòu)造LSF方法在矩形截面簡支梁上識別單個移動荷載的表現(xiàn)。

2.1 案例設(shè)計

本算例根據(jù)Lu等[14]的研究設(shè)置了一跨20 m長的矩形截面簡支梁,其截面慣性矩為0.372 m4,材料的楊氏模量、密度與泊松比分別為3.45×104MPa,2.3×103kg/m3和0.2。在梁的1/2和3/4跨底面設(shè)置了2個動撓度響應(yīng)輸出點(diǎn),輸出頻率為200 Hz。移動荷載P(t)=2×[30+2 sin(10πt)+4 cos(15πt)+cos(3πt+0.2π)]kN沿縱軸線從左至右駛過橋面(見圖2)。仿真數(shù)據(jù)均由有限元軟件ANSYS19.0計算獲得。

圖2 仿真模型示意圖(m)Fig.2 Simulation test model(m)

以測量噪聲、車速、響應(yīng)數(shù)量為主要考慮因素設(shè)置相應(yīng)工況。其中測量噪聲以白噪聲模擬

(8)

2.2 構(gòu)造LSF擬合過程

設(shè)定LSF的頻率fLSF[13]為

(9)

式中:fs為采樣頻率;l為LSF的單元步長,即“梁”每個單元的長度。該步長與“梁”的單元數(shù)m存在以下關(guān)系

T=l×m

(10)

為能準(zhǔn)確逼近荷載,要求fLSF不小于荷載的最高分析頻率。對于未知荷載而言,由于線性結(jié)構(gòu)的響應(yīng)線性依賴于外荷載,因此fLSF的確定可以通過對結(jié)構(gòu)實測響應(yīng)的頻譜分析決定。

圖3 1/2跨動撓度響應(yīng)的傅里葉變換Fig.3 Fourier transform of mid-span deflection

2.3 識別結(jié)果分析

圖4展示了各工況MFI的結(jié)果。其中圖4(a)~圖4(c)分別顯示當(dāng)移動荷載速度為13.3 m/s時,通過1/2跨、3/4跨動撓度與雙動撓度響應(yīng)識別的結(jié)果,圖4(d)~圖4(f)則顯示了移動荷載速度為22.2 m/s時的識別結(jié)果。為了更清楚的對比識別結(jié)果的精度,定義相對誤差(relevant percentage error, RPE),結(jié)果如表1。

(11)

表1 移動荷載識別結(jié)果的相對誤差Tab.1 RPE of MFI results

整體上,所有工況中的移動荷載都被識別,最高誤差不超過10%,符合實際工程需求。在不對式(7)使用正則化約束的情況下,即使響應(yīng)中存在20%的大噪聲,求解依然沒有出現(xiàn)病態(tài),驗證了構(gòu)造LSF方法在消除反卷積法求逆病態(tài)問題的優(yōu)勢。

(1)響應(yīng)位置對識別結(jié)果的影響

對比圖4(a)、圖4(b)和對應(yīng)RPE,可以發(fā)現(xiàn)不同位置響應(yīng)對于識別精度影響不大,其RPE穩(wěn)定在5.5%~8.5%。證明該方法在識別過程中對測點(diǎn)位置并無要求,于實際應(yīng)用中,可選擇相對簡便的位置進(jìn)行布測,大幅度的提高了方法的靈活性。分析圖4(c)則顯示使用多個響應(yīng)進(jìn)行識別時,在無噪聲工況下對精度影響不明顯;隨著噪聲的引入,其對精度有提升。由理論分析可看出只要保證測點(diǎn)數(shù)量不小于待測未知量數(shù)即可保證識別的穩(wěn)定,在具體工程中,可根據(jù)精度要求和計算成本酌情考慮測點(diǎn)的數(shù)量,需要注意的是各測點(diǎn)獨(dú)立性,保證響應(yīng)數(shù)據(jù)不相關(guān)。

圖4 移動荷載識別結(jié)果Fig.4 MFI results

(2)移動荷載速度對識別結(jié)果的影響

整體上,各工況識別時程的開始階段(0~0.2 s)均出現(xiàn)了較大誤差,究其原因是由車-橋耦合振動不穩(wěn)定引起的[15-16]。分析RPE可發(fā)現(xiàn),高速荷載的識別精度明顯低于低速荷載,其誤差增加了1%~2%,在噪聲較大時尤為明顯。深入分析發(fā)現(xiàn),隨著車速增長,識別時長縮短,而首端不穩(wěn)定振動區(qū)間變化不大,其產(chǎn)生的大誤差在整體誤差中占比提升,從而導(dǎo)致整體識別精度的下降。考慮到實際應(yīng)用中,識別時間相對較長,識別的誤差將控制在可接受的范圍內(nèi)。

(3)噪聲對識別結(jié)果的影響

觀察圖4(a)~圖4(c),在低速工況中,不同的噪聲小幅度地影響識別峰值的準(zhǔn)確度,對整體線型影響可忽略不計。而在高速工況中(圖4(d)~圖4(f)),雖然成功識別線型,但在大噪聲下識別的峰值出現(xiàn)較明顯的誤差。特別地,當(dāng)測量噪聲為20%時,高速工況(見圖4(d)~圖4(f))的末端已經(jīng)出現(xiàn)因誤差累計產(chǎn)生的線型畸變,由此產(chǎn)生了較大的RPE。不僅如此,在高速工況中,識別結(jié)果有小幅度的相位延遲,其原因在于有限元軟件計算構(gòu)造LSF矩陣時有一定延遲,可通過減小求解荷載步在一定程度上消除,進(jìn)一步減小誤差。

3 試驗驗證

本章展示了兩個試驗結(jié)果,分別驗證本方法在集中荷載與移動荷載識別中的應(yīng)用。試驗①建造了一跨混凝土箱梁模型并利用液壓千斤頂在跨中施加一個集中的時變荷載;試驗②利用亞克力板建造了一跨小型箱梁,利用模型車施加了一個移動荷載。兩個試驗均通過模態(tài)試驗建立目標(biāo)橋梁的有限元模型,然后計算構(gòu)造LSF矩陣,最后通過少量測量響應(yīng)完成荷載識別。

3.1 集中荷載識別的試驗驗證

依據(jù)文獻(xiàn)[17-18]的研究,本文設(shè)計了印度Mass Rapid Transit System Bridge的1/6比例模型,其跨度為6 m,截面尺寸如圖5(a)所示,成型結(jié)構(gòu)如圖5(b);使用Moog液壓系統(tǒng)在跨中施加固定時變荷載(見圖5(c));在1/2跨、3/4跨布置了兩個TML CDP-50位移計(見圖5(d)),并在5/16跨、7/16跨位置布置了兩個TML PL-60-11應(yīng)變片(見圖5(e));數(shù)據(jù)由Ni DAQ系統(tǒng)(見圖5(f))與LabVIEW軟件同步采集,采樣頻率為1 000 Hz。

圖5 試驗①模型及設(shè)備Fig.5 Lab test ① model and equipment

如2.2節(jié)所述,構(gòu)造LSF矩陣一般由結(jié)構(gòu)的有限元模型計算,需要先對試驗?zāi)P瓦M(jìn)行有限元模型修正。首先使用柱件壓縮測試與回彈儀測試,依據(jù)AS3600-2009[19]測定結(jié)構(gòu)混凝土的真實楊氏模量如表2所示。

表2 混凝土材料參數(shù)修正Tab.2 Correction of concrete material properties

隨后進(jìn)行環(huán)境振動測試(operational modal analysis, OMA),使用錘擊模擬隨機(jī)環(huán)境激勵,通過ARTeMIS模態(tài)測試軟件提取模型的頻率與模態(tài)[20],使用響應(yīng)面方法進(jìn)行有限元模型修正[21]。表3對比了實測頻率與模型計算頻率,可以看出模型的前3階模態(tài)均能較好的吻合試驗?zāi)P偷膶崪y值。雖然第4階、第6階模態(tài)在模擬中并未識別,但實際工程中,結(jié)構(gòu)響應(yīng)以前3階模態(tài)分析為主,該模型可滿足要求。

表3 試驗①的有限元模型參數(shù)修正Tab.3 FE model updating of test ①

3.2 集中荷載識別結(jié)果與分析

本試驗中,施加于跨中的時變集中荷載信號和對應(yīng)的響應(yīng)信號如圖6所示。荷載分為線性加載,循環(huán)荷載和線性減載3個階段,加載時長共53 s。響應(yīng)信號為1/2,3/4跨的豎向位移和5/16,7/16跨的軸向應(yīng)變。注意響應(yīng)信號經(jīng)過靜力加載標(biāo)定,不計梁初始撓度與應(yīng)變。

步驟1依據(jù)試驗工況,在有限元模型的跨中位置施加單位脈沖激勵,獲取1/2,3/4跨動撓度的脈沖響應(yīng)與5/16,7/16跨軸向應(yīng)變的脈沖響應(yīng)。

從圖7可以看出,各工況識別結(jié)果與真實荷載有較高的吻合度,說明所提方法能有效地識別未知荷載。

在第一階段的線性加載中,通過動撓度和應(yīng)變響應(yīng)均能有效識別荷載。但同數(shù)值仿真案例一樣,在識別初期出現(xiàn)了由于振動不穩(wěn)定產(chǎn)生的大誤差問題,尤其體現(xiàn)在應(yīng)變識別結(jié)果中(見圖7(c))。

在第二階段的循環(huán)荷載中,觀察圖7(a)和圖7(b)可發(fā)現(xiàn),由1/2跨動撓度識別的結(jié)果在波峰處吻合度較高,而在波谷處有一定誤差,將響應(yīng)數(shù)量增加為2個后,波谷處識別準(zhǔn)確度有一定提高。而對比圖7(c)和圖7(d)可以發(fā)現(xiàn),由單個應(yīng)變信號識別結(jié)果會因殘留應(yīng)力累積而導(dǎo)致識別值出現(xiàn)增大的趨勢。而增加響應(yīng)數(shù)量可以有效抑制該趨勢,并提高峰值的吻合度。

在第三階段的線性減載中,使用單個響應(yīng)的工況(見圖7(a)和圖7(c))均因誤差累積而出現(xiàn)較明顯的差距,而使用雙響應(yīng)的工況吻合度則較高,進(jìn)一步說明使用多個響應(yīng)進(jìn)行識別能有效抑制累積誤差。

圖6 激勵與響應(yīng)信號Fig.6 Excitation and responses signals

圖7 集中荷載識別結(jié)果Fig.7 Point load identification results

表4 記錄了不同工況下識別結(jié)果的相對誤差與運(yùn)算效率,所采用的計算環(huán)境為:便攜式計算機(jī),4 CPU:Intel(R) Core(TM) i7-4720HQ@2.60 GHz,RAM:8 G,MATLAB version:R2016b。運(yùn)算效率采用計算耗時與加載時長的比值表述。

表4 各工況識別結(jié)果相對誤差與運(yùn)算效率Tab.4 RPE and computational efficiency

表4可看出,以動撓度響應(yīng)反演的結(jié)果較應(yīng)變響應(yīng)更加準(zhǔn)確。這是由于箱梁截面變形不同步,布于底板應(yīng)變片捕捉的局部應(yīng)變與梁體整體變形差距導(dǎo)致的,而豎向動撓度受該影響較小,因此能得到較高精度的結(jié)果,以雙動撓度響應(yīng)反演工況為例,其精度可達(dá)到10%左右,已經(jīng)滿足相關(guān)行業(yè)標(biāo)準(zhǔn)對車輛荷載監(jiān)測的精度要求[22]。

在運(yùn)算效率方面,主要CPU耗時發(fā)生在求解反卷積方程中,其效率與系統(tǒng)矩陣大小正相關(guān)。以單個響應(yīng)反演工況為例,傳統(tǒng)反卷積法(見式(3))的系統(tǒng)矩陣維度高達(dá)(53 000×53 000);引入了構(gòu)造LSF矩陣(見式(7))后的系統(tǒng)矩陣維度僅為(53 000×1 062),減少了50多倍。

3.3 移動荷載識別的試驗驗證

該試驗?zāi)P褪怯蓙喛肆Π逯瞥傻膯蜗鋯问液喼Я海骺鐬? m,兩端分別配置3 m的引橋與1 m的緩沖板,其整體布置如圖8(a)所示,截面如圖8(b)所示。主梁沿線均勻分布9塊橫隔板,提高了梁的整體性和延性,支座采用橡膠支座。橋梁的單位長度密度為4.383 3 kg/m,抗彎剛度為3 550 N·m2。使用電動越野車模擬移動荷載,其動力由自身馬達(dá)提供。小車的等效荷載為18.13 kN(見圖8(c))。梁跨中部位設(shè)置120-3AA型電阻應(yīng)變片,測量跨中應(yīng)變;在橋頭、中、尾部架設(shè)3個L形框架,其橫梁上各布置一個測速應(yīng)變片,如圖8(d),通過記錄小車觸碰應(yīng)變片的時間間隔來計算小車速度。

車輛在引橋上加速,然后以相對恒定的速度通過橋梁。在梁的上表面安裝了兩個輕質(zhì)泡沫導(dǎo)軌,以防止移動的車輛墜落。信號由INV3060V(2)型智能信號采集處理分析儀收集,采樣頻率為2 048 Hz。

圖8 試驗②模型及設(shè)備Fig.8 Lab test ② model and equipment

模型同樣使用OMA進(jìn)行模態(tài)特征采集與模型修正,前4階的模型修正結(jié)果如表5。

表5 試驗②的模型修正Tab.5 Model updating of test ②

3.4 移動荷載識別結(jié)果與分析

模型小車在通過橋梁時,其速度由3個測速應(yīng)變片的響應(yīng)時差確定。如圖9所示,選擇各應(yīng)變片的最大響應(yīng)作為標(biāo)尺。需要說明的是,①本工況假定車輛以恒定速度穿過測試梁;②該方法測出的時差是前車軸駛?cè)搿Ⅰ偝鰳蛄旱臅r差,而本文將小車作為一個等效移動荷載進(jìn)行識別,并未考慮前后軸區(qū)別,所以存在一定誤差,在以后的研究中將進(jìn)行優(yōu)化。由響應(yīng)時差計算出的車速為2.87 m/s,同時確定車上、下橋的時間點(diǎn)為3.521 s和4.565 s,取識別時間為3.6~4.6 s。相應(yīng)時間內(nèi)的跨中應(yīng)變響應(yīng)如圖10所示,識別結(jié)果如圖11所示。

圖9 測速應(yīng)變片的觸發(fā)時間Fig.9 Strain gauge active time

圖10 跨中應(yīng)變響應(yīng)Fig.10 Strain response of mid-span

圖11 移動荷載識別結(jié)果Fig.11 Moving force identification result

通過圖11可知,本方法識別出的是車輛與橋面的瞬時接觸力,幅值與車輛的等效荷載相當(dāng)。由于本工況在計算中并未考慮車輛的軸重,而是將其等效為一個移動荷載進(jìn)行識別,所以在車輛上、下橋的時段,由于前后軸分別位于主橋與引橋上,其接觸力產(chǎn)生了較大變化,導(dǎo)致識別失敗;但車輛全部位于主橋時,識別的結(jié)果還是較為滿意的。

4 結(jié) 論

針對既有反卷積法識別移動荷載過程中耗時長、奇異性大的問題,本文提出了引入構(gòu)造LSF的MFI改進(jìn)方法。該方法使用有限個數(shù)的構(gòu)造LSF,將反卷積方程中的脈沖響應(yīng)矩陣擬合為較平滑的類對角矩陣,同步減少了矩陣的維度和奇異性,大幅度提升了識別效率和魯棒性。采用數(shù)值仿真與模型試驗共同驗證所提方法的有效性,綜合全文所得主要結(jié)論如下:

(1)利用構(gòu)造LSF可提升反卷積法對噪聲的魯棒性,從而省略正則化流程。在數(shù)值仿真中,即使附加20%的白噪聲,該方法仍能保障識別穩(wěn)定性,并保持較高的精度。

(2)所提方法僅用1~2個響應(yīng)信號即能較為準(zhǔn)確地識別時變荷載。在模型試驗①中,針對加載時間長、激勵類型復(fù)雜的集中荷載,使用單個響應(yīng)信號即可得到較好的結(jié)果,增加響應(yīng)個數(shù)可提升精度,但對比運(yùn)算量的增加,收益不高。基于動撓度響應(yīng)反演的精度可保持10%左右,與壓電薄膜型車輛稱重系統(tǒng)誤差持平,但本方法成本不超過稱重系統(tǒng)的10%,有優(yōu)秀的應(yīng)用前景。

(3)在模型試驗②中,該方法能有效地識別2軸移動小車的等效荷載,但是在車輛上、下橋時由于前后軸與主橋接觸力的劇烈變化導(dǎo)致識別失效,在將來的研究中需要針對此問題進(jìn)行進(jìn)一步地探討。

Vol.41 No.12 2022

Vol.41 No.12 2022

猜你喜歡
橋梁方法模型
一半模型
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計的漸近分布
手拉手 共搭愛的橋梁
句子也需要橋梁
高性能砼在橋梁中的應(yīng)用
3D打印中的模型分割與打包
用對方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
捕魚
主站蜘蛛池模板: 尤物在线观看乱码| 91久久国产热精品免费| 国产爽歪歪免费视频在线观看| 狠狠色综合网| 69av在线| 亚洲一级毛片在线播放| 中文字幕在线看视频一区二区三区| 伊人大杳蕉中文无码| 亚洲熟女偷拍| 在线国产毛片手机小视频| 97国产在线视频| 国产成人精品在线1区| 在线中文字幕网| 国产精品综合色区在线观看| 国产成人艳妇AA视频在线| 亚洲黄色视频在线观看一区| 国产日韩欧美精品区性色| 国产精品自在自线免费观看| 午夜a级毛片| 亚洲欧美极品| 精品一区二区无码av| 色男人的天堂久久综合| 精品伊人久久久香线蕉| 免费一极毛片| 亚洲一区波多野结衣二区三区| 色综合婷婷| 99这里只有精品免费视频| 91久久精品日日躁夜夜躁欧美| 日本一区二区不卡视频| 国产成人一二三| 国产亚洲精品97AA片在线播放| 怡红院美国分院一区二区| 欧美成人精品高清在线下载| 国产精品亚洲va在线观看| 国产91精选在线观看| 99热这里只有精品在线观看| 亚洲熟女中文字幕男人总站| 国产激爽大片高清在线观看| 国产精品福利社| 国产精品亚洲精品爽爽 | 免费99精品国产自在现线| 国产爽爽视频| 国产日韩丝袜一二三区| 亚洲午夜国产精品无卡| 54pao国产成人免费视频| 国产精品视频导航| 国产极品粉嫩小泬免费看| 成人免费网站在线观看| 精品国产成人国产在线| 欧美综合中文字幕久久| 国产色婷婷视频在线观看| 亚洲av无码久久无遮挡| 国产精品手机在线观看你懂的| 欧美黑人欧美精品刺激| 九一九色国产| 国产女同自拍视频| 在线五月婷婷| 欧美成人精品在线| 2022国产无码在线| 日韩一级二级三级| 国产精品美女自慰喷水| 欧美另类第一页| 免费在线a视频| 91福利免费视频| 精品精品国产高清A毛片| 日本妇乱子伦视频| 亚洲婷婷丁香| 999国产精品永久免费视频精品久久 | 亚洲人成影院午夜网站| 国产高清国内精品福利| 亚洲高清无码久久久| 中文国产成人久久精品小说| 午夜国产精品视频黄| 天天摸天天操免费播放小视频| 亚洲综合专区| 亚洲天堂视频网| 四虎精品国产永久在线观看| 国产菊爆视频在线观看| 欧美成人国产| 国产成人精品一区二区| 视频二区中文无码| 乱系列中文字幕在线视频|