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

桅桿結構有限元模型修正與參數敏感性研究

2024-12-31 00:00:00劉慕廣喬磊王雷余先鋒張春生謝壯寧張麗
振動工程學報 2024年7期

摘要: 基于臺風“圓規”激勵下的實測加速度響應識別了深圳356 m高氣象梯度觀測塔的模態參數,結合帶精英策略的快速非支配排序遺傳算法(NSGA?Ⅱ)對氣象塔有限元模型進行修正。結果表明:氣象塔模態非常密集,且纖繩模態的參與程度較為顯著。氣象塔X向和Y向的基頻分別為0.614 Hz和0.603 Hz,其前3階彎曲模態阻尼比在1%~2%之間。塔身密度、纖繩彈模對塔身模態頻率和振型有顯著影響,高層纖繩的線質量和塔身彈模對其也有一定影響,但纖繩張力的影響較低。氣象塔有限元模型修正后的風致響應高于修正前,并更接近實測結果,驗證了修正模型的準確性。

關鍵詞: 模型修正; 風振響應; 桅桿; 模態識別; 實測

中圖分類號: TU311.3; TU323.4""" 文獻標志碼: A""" 文章編號: 1004-4523(2024)07-1107-08

DOI:10.16385/j.cnki.issn.1004-4523.2024.07.003

收稿日期: 2022-11-09; 修訂日期: 2023-01-19

基金項目:"國家自然科學基金資助項目(51978285,52378514);廣東省現代土木工程技術重點實驗室基金資助項目(2021B1212040003);亞熱帶建筑與城市科學全國重點實驗室開放基金資助項目(2021ZB11);廣東省基礎與應用基礎研究基金資助項目(2024A1515011828,2024A1515011525)。

引 言

精確的有限元模型對于研究土木結構的靜動力響應和開展結構安全二級評估等至關重要。但是,由于土木結構體型龐大、參數繁多,且存在各類不確定性和非線性因素[1],按照設計圖紙建立的有限元模型不可避免地存在各種誤差[2],最終會影響到有限元模型分析結果的準確性。因此,極有必要以現場實測的數據為基準優化有限元模型,使有限元與實際結構的模態特征盡量吻合,這對于高聳或大跨等土木結構尤為重要。

有限元模型修正方法可分為矩陣優化法[3]和基于敏感性分析的參數型法[4]。矩陣優化法需要事先精確估計結構的質量、剛度矩陣,對于復雜工程,實際應用中通常難以滿足[5]。相比來說,參數型方法以結構物理參數(如材料特性、幾何特性、邊界條件等)作為修正對象,以有限元響應特征值與實測值的殘差為目標函數,運用優化算法迭代調整物理參數,使目標函數最優。該類方法具有明確的物理意義,在有限元修正中被較多采用。張建等[6]以2019年北京世界園藝博覽會國際竹藤組織館為實測對象,基于實測模態參數對該結構的有限元模型進行了修正。Saudi[7]基于一座90 m桅桿通信塔現場實測的固有頻率修正了有限元模型,并進一步評估了桅桿的結構安全。Ni等[8]基于實測數據對廣州塔有限元模型進行了修正。葉錫鈞等[9]也以廣州塔為實測對象,利用遺傳算法對廣州塔的初始參數進行了修正。Foti等[10]對一座意大利塔樓進行多次環境振動實測,基于實測模態采用敏感性參數分析法修正了有限元模型。Ren等[11]基于實測的固有頻率構造二次多項式響應曲面,通過修正彈性模量、節點區域面積來優化有限元模型,并通過振動試驗驗證有限元修正的有效性。Brownjohn等[12]利用實測頻率和振型,采用基于靈敏度的參數型修正方法修正橋梁有限元模型,且以修正模型為基準評估了實際結構的損傷狀況。

本文以深圳市氣象梯度觀測塔這一高聳桅桿結構為研究對象,采用2021年臺風“圓規”激勵下現場實測的加速度信號識別塔身的模態參數。結合參數型模型修正方法,對影響桅桿結構模態特征的參數進行了敏感性研究,通過帶精英策略的快速非支配排序遺傳算法(NSGA?Ⅱ)[13]對有限元模型進行修正,并對修正前后的風致響應特征進行了初步對比分析。

1 桅桿結構有限元模型

深圳市氣象梯度觀測塔位于深圳市寶安區,是亞洲第一、世界第二高的桅桿格構塔。該塔高356 m,由5層纖繩固定,在東南西北各方向分別設置3個錨固點,如圖1所示。塔身內部設置一臺載重量為500 kg的齒條式升降機。塔身為桁架結構,采用Q345B熱軋無縫鋼管。塔身橫斷面為正方形,在地面0標高處為5 m5 m,至15 m標高處漸變為2.5 m2.5 m并保持不變。纖繩采用半平行鍍鋅鋼絲拉索,與地面傾角在44°~56°范圍內,在65,195和330 m各方向單根布置,在130,260 m各方向雙根布置。

采用ANSYS建立氣象塔有限元模型,原型中的觀測平臺、爬梯、電梯、螺栓、法蘭板等附屬設備未在模型中直接建模,而是通過對塔身密度適當放大來簡化考慮。塔身各桿件采用BEAM188單元,彈性模量為206 GPa,考慮附屬結構后的等效密度取7890 kg/m3。纖繩采用LINK180單元,彈性模量為185 GPa,每層纖繩線質量和張力如表1所示,有限元模型中通過施加初應變的方式來模擬纖繩的預張力。塔頂避雷針等設備通過質量單元MASS21施加在模型上。有限元模型共包含9555個單元、4190個節點、25116個自由度,如圖2所示。

通過ANSYS的Block Lanczos大變形預應力模態分析法對有限元模型進行模態計算。表2和3分別給出了頻率低于1.0 Hz的纖繩(南北向)和塔身模態分析結果。有限元結果表明,氣象塔塔身的X向、Y向前3階模態主要集中在0.7~0.9 Hz,纖繩在1 Hz以內出現了10多階模態,且部分頻率與塔身較為接近。雖然南北向和東西向纖繩地面標高不同(見圖1)會導致兩個方向的纖繩的模態頻率略不同,但考慮到東西向纖繩頻率與南北向差異并不大,表2中不再列出。

2 基于實測的模態參數識別

2.1 實測數據

華南理工大學在氣象塔塔身50,160,250,300和350 m 5個高度布置了X,Y雙向加速度儀(LAC?Ⅱ型)。設備采樣頻率為25 Hz,其X,Y正向分別對應氣象塔正東向和正北向。

2021年第18號臺風“圓規”于10月8日下午在菲律賓以東洋面生成,10月13日5時,“圓規”加強為臺風級,中心附近最大風力12級,10月13日15時40分前后,“圓規”在海南省瓊海市沿海登陸。臺風“圓規”的中心距離深圳氣象梯度塔最近時不足400 km,且2021年10月12—13日,氣象塔基本處于臺風“圓規”的7級風圈內。圖3為2021年10月12日0時至14日0時氣象塔320 m高度處的風速和風向樣本,圖4為對應氣象塔的加速度響應。由圖3可見,自12日0時起,氣象塔處風速穩步增大,在13日12時前后達到最大,對應的10 min平均風速為16 m/s,隨后風速逐漸降低;風向則由12日0時的北風逐漸在13日下午變為東風。由圖4可見,隨風速增加,氣象塔的加速度響應也逐漸增大,且響應總體隨高度的增加而增大,平均風速最大時氣象塔350 m處X向峰值加速度為19 cm/s2,Y向峰值加速度為20 cm/s2。下文統一選取2021年10月12日22:00—24:00持續大風時段數據作為分析樣本。

2.2 模態參數識別

圖5為各實測高度X向和Y向截止頻率為1.0 Hz的加速度功率譜曲線。由圖5可見,桅桿結構在X向和Y向分別呈現出近10個較為明顯的能量峰值。結合表2和3的有限元結果,圖5中桅桿的加速度功率譜中必然同時混雜了纖繩和塔身的振動信息。考慮到塔身振動的整體性和纖繩振動對塔身影響的局部性,利用穩定圖[14]進一步剔除纖繩局部振動引起的虛假模態,由5個高度的加速度響應得到的頻率穩定圖如圖6所示。

由圖6可見,X向、Y向分別有6個穩定的頻率極點。對于X向、Y向0.58 Hz以內的頻率,其能量峰僅在2,3個高度位置較為明顯,進一步通過不同高度信號的互功率譜相位信息及振型特征,并結合表2和3中的有限元結果,推斷出X向的0.182,0.238和0.576 Hz,Y向的0.210,0.239和0.573 Hz分別為不同錨固高度纖繩的頻率。同理,可推斷0.614,0.728和0.837 Hz分別為塔身X向前3階模態;0.603,0.713和0.822 Hz分別為塔身Y向前3階模態。圖7給出了實測和有限元的振型對比,由圖7可見,X向、Y向的前3階實測振型與有限元結果具有較好的一致性。

表4中匯總了氣象塔X向、Y向的前3階彎曲模態頻率實測結果,表中同時給出了基于隨機子空間法識別的塔身模態阻尼比。由表4可見,塔身的前3階阻尼比在1%~2%間變化,具有一定的離散性,這一結果與Harikrishna等[15]針對50 m高桅桿實測的1%~3%阻尼比相近,也接近中國規范[16]中鋼塔架2%的阻尼比建議值。

另外,表4中對比了氣象塔有限元模型與實測頻率的誤差及振型MAC(Modal Assurance Criterion)。由表4可見,有限元模型的前3階模態均與實測存在較大差異,振型MAC值均在0.75以下,其X向、Y向最大頻率誤差分別為18.08%和19.24%,均出現在1階模態。這些差異的出現可能是模型中各項設計參數與結構實際服役狀態不符的緣故。為了得到一個比較理想的桅桿塔身動力模型,需要對初始有限元模型進行修正。

3 有限元模型修正

3.1 參數敏感性研究

由于桅桿結構的非線性特征,其基本動力特性與結構組成形式、纖繩分布、結構剛度等相互關聯[17]。結合桅桿結構的實際構造和工程實際,初步選定的影響參數主要有:塔身彈模Et、塔身密度ρ、第1~5層纖繩線質量ρ1~ρ5、第1~5層纖繩預張力F1~F5和纖繩彈模Er。

首先分析塔身模態頻率和振型對以上參數的敏感性。參數的敏感性研究實質上是函數對自變量求導的數學問題,有限元差分法是常用的靈敏度數值分析方法。利用攝動法使修正參數發生微小擾動,再根據一階差分公式近似計算參數靈敏度, 一階差分公式為:

(1)

式中 模態特征量f可以是頻率、振型MAC等;xi為修正參數;Δxi為攝動量,取設計值的10%;各參數的設計值如第1節所述。

模態特征量對各參數的靈敏度分析如圖8和9所示。圖中1,3,5階模態分別對應Y向的1,2,3階彎曲模態;2,4,6階模態分別對應X向的1,2,3階彎曲模態。從圖8中頻率的敏感性可以看出:1)塔身模態頻率對塔身彈模和密度、纖繩彈模的敏感性最高,纖繩質量次之,纖繩張力最小;2)塔身密度和纖繩彈模對塔身X向、Y向的前3階彎曲模態頻率均有明顯影響,但塔身彈模僅對X向、Y向的第3階模態頻率存在較大影響;3)3~5層纖繩的線質量和預張力對塔身模態頻率影響較大,底部兩層纖繩的影響很小。

由圖9中對模態振型的敏感性可以看出:1)塔身振型對塔身密度、3~5層纖繩線質量及纖繩彈模均具有較高的敏感性,塔身彈模和1層、2層纖繩線質量對塔身振型的影響不大;2)3~5層纖繩預張力對模態振型也有一定影響,底部兩層纖繩的影響很小。

依據圖8和9中靈敏度分析結果,下文選取塔身彈模、塔身密度、纖繩彈模及第3~5層纖繩的線質量和預張力作為修正參數,以實測結果為目標,對塔身X向、Y向的前3階模態進行修正。

3.2 NSGA?Ⅱ算法與修正結果

NSGA?Ⅱ算法降低了非劣排序的復雜性,具有運算速度快,解集收斂性好的優點[13]。本文通過設計MATLAB?ANSYS聯合仿真優化程序,利用全局搜索和局部優化實現模型修正,每次迭代的參數由MATLAB導入ANSYS中,目標函數值通過MATLAB讀取ANSYS輸出文件進而分析優化。算法參數設置:初始種群500個,最大迭代次數20代,交叉概率0.7,變異概率0.02。根據工程經驗并經初步試算,表5給出了各參數的初始值與變化范圍。

通過定義頻率誤差平方和Q1、振型MAC之和Q2建立目標函數,構造的目標函數如下:

(2)

(3)

式中 ,分別為X向、Y向的第i階有限元計算頻率;,分別為X向、Y向的第i階實測頻率;,分別為X向第i階有限元振型向量和實測振型向量(由于表4中Y向振型MAC差別很小,目標函數Q2中僅考慮X向前3階振型)。

NSGA?Ⅱ算法對目標函數的優化過程如圖10所示。迭代至第9代滿足收斂條件,對應目標函數Q1收斂至0.009,目標函數Q2收斂至2.8(各階頻率誤差小于6%,各階振型MAC大于0.9)。

表6為修正后的參數值及變化率,表7為修正后塔身X向、Y向的前3階彎曲模態。由表6可見,纖繩線質量、張力及彈模變化最顯著,塔身彈模和密度較參數初始值變化較小。另外,修正后塔身彈模和密度較初始值略有增大;對于纖繩,除第5層的線質量比初始值大外,其他參數均有較大程度的減小。纖繩參數的變化率表明實際服役纖繩的一些性能較設計值可能存在一定衰減。由表7可見,修正后有限元模型X向、Y向的前3階頻率與實測的誤差均在5.31%以內,振型MAC不低于0.93,較表4中初始有限元模型有較大程度的改善,其相對誤差可以接受[18]。

3.3 有限元和實測結果的對比

為了進一步驗證模型修正的合理性,本節對梯度塔有限元風振結果與實測數據進行對比。由于2021年10月12日多個風速儀的數據存在異常,本節選用13日12:00—12:10的大風數據及梯度塔加速度實測數據。將梯度塔不同高度13臺超聲風速計(0.1 Hz采樣頻率)的平均風速進行指數律擬合,得到梯度塔對應時刻的風速沿高度的分布,如圖11所示。由于所選時段臺風已登錄海南,擬合到的剖面指數明顯偏大。通過對10 Hz采樣的實測風速數據進行頻譜分析發現,實測風譜與von Karman譜較為吻合,如圖12所示,圖中Lu為湍流積分尺度,U為對應高度處風速。通過以上實測風場的分析,采用諧波疊加法,基于實測風速擬合剖面每10 m高度生成一個風速時程。風速譜采用von Karman譜,僅考慮塔身高度方向脈動風的相關性,相關函數采用Simiu建議值(Cz=10),脈動風間隔取0.1 s。

采用完全法進行瞬態動力分析,阻尼比按表4中X向1階模態0.0209取值。通過分析,得到有限元模型修正前、后塔身響應的加速度均方根,其X向結果與實測值的對比如圖13所示。

由圖13可見,修正后模型加速度均方根要高于修正前,且更接近實測值。修正后上部四個測點值與實測值誤差均小于5%,但在50 m處兩者差異較為明顯,其誤差為-35.35%,這可能是因為實測中近地面風場更易受周邊干擾,高湍流導致結構脈動響應偏大。

4 結 論

(1)深圳氣象塔模態非常密集,且纖繩模態參與程度較為顯著。X向、Y向的基頻分別為0.614 Hz和0.603 Hz,結構阻尼比為1%~2%。

(2)塔身密度、纖繩彈模對動力特性有顯著影響,塔身彈模、高層纖繩的線質量也存在一定程度的影響,纖繩張力對其影響較低。

(3)結合NSGA?Ⅱ算法修正了有限元模型的材料參數,修正后的塔身彈模和密度略有增大,纖繩張力、彈模等參數明顯減小,說明實際服役纖繩的一些性能較設計值可能存在一定衰減。

(4)氣象塔有限元模型修正后的風致響應高于修正前,并更接近實測結果,驗證了修正模型的準確性。

參考文獻:

[1]"""" Asgarieh E, Moaveni B, Barbosa A R, et al. Nonlinear model calibration of a shear wall building using time and frequency data features[J]. Mechanical Systems and Signal Processing, 2017, 85: 236-251.

[2]"""" Yuan Z X, Liang P, Silva T, et al. Parameter selection for model updating with global sensitivity analysis[J]. Mechanical Systems and Signal Processing, 2019, 115: 483-496.

[3]"""" Berman A, Nagy E J. Improvement of a large analytical model using test data[J]. AIAA Journal, 1983, 21(8): 1168-1173.

[4]"""" Mazzotti M, Mao Q, Bartoli I, et al. A multiplicative regularized Gauss-Newton method with trust region sequential quadratic programming for structural model updating[J]. Mechanical Systems and Signal Processing, 2019, 131: 417-433.

[5]"""" 楊朋超, 薛松濤, 謝麗宇. 結構動力模型的改進直接修正方法及工程應用[J]. 建筑結構學報, 2021, 42(3): 34-40.

YANG Pengchao, XUE Songtao, XIE Liyu. An improved direct method for dynamic model updating and its practical engineering applications[J]. Journal of Building Structures, 2021, 42(3): 34-40.

[6]"""" 張建, 錢孝文, 楊娜, 等. 大跨度圓竹拱結構動力測試及有限元模型修正[J]. 建筑結構學報, 2022, 43(2): 94-104.

ZHANG Jian, QIAN Xiaowen, YANG Na, et al. Dynamic test and finite element model modification of long-span round bamboo arch structures[J]. Journal of Building Structures, 2022, 43(2): 94-104.

[7]"""" Saudi G. Structural assessment of a guyed mast through measurement of natural frequencies[J]. Engineering Structures, 2014, 59: 104-112.

[8]"""" Ni Y Q, Xia Y, Lin W, et al. SHM benchmark for high-rise structures: a reduced-order finite element model and field measurement data[J]. Smart Structures amp; Systems, 2012, 10(4): 411-426.

[9]"""" 葉錫鈞, 安關峰, 周朝陽, 等. 高聳柔性結構的模態參數識別及模型修正研究[J]. 建筑結構學報, 2014, 35(5): 33-39.

YE Xijun, AN Guanfeng, ZHOU Chaoyang, et al. Modal parameter identification and model updating of high-rise flexible structure[J]. Journal of Building Structures, 2014, 35(5): 33-39.

[10]""" Foti D, Diaferio M, Giannoccaro N I, et al. Ambient vibration testing, dynamic identification and model updating of a historic tower[J]. NDT amp; E International, 2012, 47: 88-95.

[11]""" Ren W X, Chen H B. Finite element model updating in structural dynamics by using the response surface method[J]. Engineering Structures, 2010, 32(8): 2455-2465.

[12]""" Brownjohn J M W, Xia P Q, Hao H, et al. Civil structure condition assessment by FE model updating: methodology and case studies[J]. Finite Elements in Analysis and Design, 2001, 37(10): 761-775.

[13]""" Deb K, Pratap A, Agarwal S, et al. A fast and elitist multiobjective genetic algorithm: NSGA-II[J]. IEEE Transactions on Evolutionary Computation, 2002, 6(2): 182-197.

[14]""" Reynders E, De Roeck G. Reference-based combined deterministic-stochastic subspace identification for experimental and operational modal analysis[J]. Mechanical Systems and Signal Processing, 2008, 22(3): 617-637.

[15]""" Harikrishna P, Annadurai A, Gomathinayagam S, et al. Full scale measurements of the structural response of a 50 m guyed mast under wind loading[J]. Engineering Structures, 2003, 25(7): 859-867.

[16]""" 中華人民共和國住房和城鄉建設部,國家市場監督管理總局. 高聳結構設計標準:GB 50135—2019[S]. 北京: 中國計劃出版社, 2019.

Ministry of Housing and Urban-Rural Development of the People’s Republic of China, State Administration for Market Regulation. Standard for design of high-rising structures: GB 50135—2019[S]. Beijing: China Planning Press, 2019.

[17]""" 王肇民. 桅桿結構[M]. 北京: 科學出版社, 2001.

WANG Zhaomin. The Guyed Masts[M]. Beijing: Science Press, 2001.

[18]""" 丁一凡, 劉宇飛, 樊健生, 等. 基于空間位形的在役索膜結構有限元模型修正與安全評估[J]. 工程力學, 2022, 39(5): 44-54.

DING Yifan, LIU Yufei, FAN Jiansheng, et al. Finite element model updating and safety assessment of in-service cable-membrane structure based on spatial configuration[J]. Engineering Mechanics, 2022, 39(5): 44-54.

Finite-element model updating and parameter sensitivity research of high-rise guyed mast

LIU Mu-guang1,2, QIAO Lei1, WANG Lei3, YU Xian-feng1, ZHANG Chun-sheng4, XIE Zhuang-ning1, ZHANG Li4

(1.School of Civil Engineering amp; Transportation, South China University of Technology, Guangzhou 510641, China; 2.State Key Laboratory of Subtropical Building and Urban Science, Guangzhou 510641, China; 3.Guangdong Communication Planning amp; Design Institute Group Co., Ltd., Guangzhou 510507, China; 4.Shenzhen National Climate Observatory, Shenzhen 518040, China)

Abstract: Based on the measured acceleration response under Typhoon Kompasu, the modal parameters of the 356 m Shenzhen Meteorological Gradient Tower (SMGT) are identified. The Non-dominated Sorting Genetic Algorithm Ⅱ(NSGA-Ⅱ), which is a fast and elitist genetic algorithm, is applied to update the finite element (FE) model of SMGT. The results show that the vibration modes of SMGT are very dense, and the involvement of cable vibration modes is obvious. The fundamental frequencies of SMGT in X and Y directions are 0.614 Hz and 0.603 Hz, respectively, and the damping ratio of the first 3 order bending modes are about 1%~2%. The tower density and cable elastic modulus have a significant effect on the modal frequency and mode shape of SMGT, the lineic mass of high-rise cable and tower elastic modulus also have a certain influence, while the cable tension has a relatively low influence on the modes of SMGT. The wind-induced response of the updating FE model is higher than that of initial model, and closer to the actual measurement, which verifies the accuracy of the updating FE model.

Key words: model updating;wind vibration response;guyed mast;modal identification;field measurements

作者簡介: 劉慕廣(1981―),男,博士,副教授。電話:(020)87110615;E-mail:liumg@scut.edu.cn。

通訊作者: 余先鋒(1985―),男,博士,講師。電話:(020)87110615;E-mail:ctxfyu@scut.edu.cn。

主站蜘蛛池模板: 99国产在线视频| 日韩中文字幕亚洲无线码| 亚洲欧美另类色图| 毛片手机在线看| 欧美一级99在线观看国产| 国产不卡网| 色婷婷成人网| 九九九久久国产精品| 亚洲综合在线最大成人| 福利片91| 国产麻豆精品在线观看| 91在线免费公开视频| 亚洲第一视频免费在线| 欧美中日韩在线| 欧美精品亚洲精品日韩专区| 亚洲高清日韩heyzo| 免费国产黄线在线观看| 99久久亚洲综合精品TS| 国产精品一线天| 操国产美女| 精品国产91爱| 欧美国产在线看| 国产96在线 | 国产精品 欧美激情 在线播放| 亚洲专区一区二区在线观看| 就去色综合| 欧美激情视频二区三区| 99视频全部免费| 日韩av资源在线| 国产毛片不卡| 99re热精品视频国产免费| 日韩毛片免费| 国产十八禁在线观看免费| 亚洲天堂网在线视频| 一级毛片网| 青青操国产视频| 成人福利在线免费观看| 98超碰在线观看| 国产清纯在线一区二区WWW| 99激情网| 婷婷伊人久久| 免费精品一区二区h| 亚洲三级电影在线播放| 亚洲精品国产日韩无码AV永久免费网 | 日韩专区第一页| 91视频青青草| 国内精品一区二区在线观看| 亚洲乱码在线播放| 亚洲免费人成影院| 欧美在线黄| 97狠狠操| 91最新精品视频发布页| 18黑白丝水手服自慰喷水网站| 99精品视频九九精品| 国产福利在线免费观看| 亚洲AV永久无码精品古装片| 日本久久久久久免费网络| 少妇人妻无码首页| 91丨九色丨首页在线播放| 成年人国产网站| 欧美丝袜高跟鞋一区二区| 国产99在线| 亚洲人成人伊人成综合网无码| 亚洲无码在线午夜电影| 国产一二三区在线| 日韩毛片在线播放| 美美女高清毛片视频免费观看| 日韩高清无码免费| 天天躁日日躁狠狠躁中文字幕| 欧美成人看片一区二区三区 | 久久综合五月婷婷| 久久99精品久久久久久不卡| 欧美不卡二区| 欧美成a人片在线观看| 久热99这里只有精品视频6| 欧美日本一区二区三区免费| 欧美另类精品一区二区三区| 亚洲天堂自拍| 国产主播在线一区| 国产成人精品优优av| 欧美综合成人| 婷婷亚洲最大|