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

基于半剛性模型風洞試驗的鍋爐塔架風振分析

2015-04-14 08:42:27吳承卉黃銘楓陳勇軍
空氣動力學學報 2015年3期
關鍵詞:模態模型

吳承卉,黃銘楓,*,姜 雄,陳勇軍,邱 建

(1.浙江大學建筑工程學院結構工程研究所,浙江杭州 310058; 2.東方電氣東方鍋爐股份有限公司,四川成都 611731)

基于半剛性模型風洞試驗的鍋爐塔架風振分析

吳承卉1,黃銘楓1,*,姜 雄1,陳勇軍2,邱 建2

(1.浙江大學建筑工程學院結構工程研究所,浙江杭州 310058; 2.東方電氣東方鍋爐股份有限公司,四川成都 611731)

對于開展高頻天平測力試驗的半剛性模型,在進行傳統的頻域內結構風振分析之前,需要對高頻天平測力得到的試驗數據進行修正。從達朗貝爾原理出發,通過假定原型塔架振型與物理模型振型之間的轉換關系,運用模態疊加法,推導得出了半剛性模型天平測力數據的多階廣義風荷載譜修正公式。利用新得到的修正公式,發展了基于半剛性模型天平測力風洞試驗的塔架結構風振分析方法。開展了某實際超大型鍋爐塔架的天平測力風洞試驗,對天平測力數據進行了修正,結合鍋爐塔架的有限元模型,并完成了鍋爐塔架的風振分析工作,驗證了鍋爐塔架設計的抗風安全性。

鍋爐塔架;風洞試驗;半剛性模型;數據修正;廣義風荷載譜;風激振動

0 引 言

隨著我國經濟建設的高速發展,高層建筑及高聳結構在各地大量新建。由于該類建筑結構具有柔性高、阻尼低等動力特性,使得動力風荷載效應成為影響高層建筑結構安全性和使用性的主要因素之一。高頻天平風洞試驗技術[1-4]是目前分析和確定高層建筑動力風荷載效應的有效手段。高頻天平風洞試驗通過直接測量模型底部風致剪力、彎矩和扭矩,在線性振型假定的基礎上可以很好地估計原型結構的模態風力[5],使得結構風致振動分析工作能夠方便快捷地完成[6]。由于試驗中所用的測試模型只需模擬原型結構的幾何建筑外形,高頻天平風洞試驗技術具有模型制作簡單、試驗方便、省時省錢等優點,所以近年來被廣泛應用于高層及高聳結構的風效率研究中[7-10]。

為了保證高頻測力天平試驗的有效性,對風洞試驗的物理模型要求做到輕質高頻[11](即為剛性模型),使得模型的自振頻率遠離風洞試驗時風荷載譜密度的主要分布頻段范圍,避免產生風激共振,從而避免底部天平測力數據受到模型振動產生的慣性力干擾。但是一些特殊結構風洞試驗的模型由于受到原型結構幾何特征和模型制作材料、工藝等的限制,如格構式高聳結構,其風洞試驗剛性模型“輕質高頻”的要求往往難以達到,模型的自振頻率不夠高,在風洞試驗過程中仍然可能產生風激振動,這類模型可稱為半剛性模型[12]。該類模型天平試驗得到的風力數據由于含有慣性力成分,需要進行修正后才能用于通常的頻域內結構風振分析。

國內外關于半剛性模型風洞天平測力試驗技術研究的文獻目前可見的還較少。鄒良浩[12]等基于結構動力學及隨機振動理論對半剛性模型風洞試驗荷載譜的修正推導了一種實用計算公式,得到了消除共振影響的格構式塔架的一階廣義荷載譜;但由于給出的修正公式只能用于計算第一階振型為平動或者扭轉時的一階廣義荷載,并不適用于有多個振型參與風激振動的半剛性模型風洞試驗,應用范圍較狹窄,公式的推廣及應用受到限制。如果半剛性模型在風洞試驗中有多個振型參與振動,如何對高頻天平數據進行修正和開展原型結構風振分析,正是本文討論的主要問題。

本文從達朗貝爾原理出發,建立了半剛性模型的運動方程,通過假定原型塔架振型與物理模型振型之間的轉換關系,推導得出了半剛性模型天平測力數據的多階廣義風荷載譜修正公式。并利用新得到的修正公式,發展了基于半剛性風洞模型天平測力試驗的塔架結構風振分析方法。對某高139 m的萬噸級大型鍋爐塔架開展了半剛性風洞模型底部天平測力試驗。利用新得到的修正公式,對天平測力數據進行了修正,并結合鍋爐塔架的有限元模型,完成了鍋爐塔架的風振分析工作。

1 塔式鍋爐風洞試驗簡介

某萬噸級超超臨界機組鍋爐塔架為一格構式塔架,高139 m,平面尺寸為35.7 m×35.7 m,高寬比接近4∶1。利用ANSYS軟件[13-14]建立鍋爐塔架原型的有限元模型,可得該鍋爐塔架一階頻率為0.256 Hz,位于風荷載頻譜段內,屬于對動力風荷載作用敏感的高聳結構。由于我國規范對鍋爐塔架這種特殊結構風荷載沒有明確規定,因此為保證結構抗風設計的安全、經濟和合理,對該高聳結構開展了高頻天平風洞試驗,為該結構的抗風設計提供合理依據。

風洞試驗在浙江大學ZD-1邊界層風洞(試驗段尺寸為:4 m×3 m×18 m)進行,采用λL=1∶150模型縮尺比,物理模型高約0.95 m,平面尺寸約為0.24 m×0.24 m。由不銹鋼管及輕質薄型木板制作,其中不銹鋼的彈性模量為195 kN/mm2,高頻天平試驗主要相似比指標參見表1。圖1為風洞試驗模型在風洞中的照片。風洞流場為B類地貌的大氣邊界層氣流,地面粗糙度α=0.16,取離地面30 m處的大氣湍流度Iu近似與地面粗糙度α相等——即Iuh=30m≈α=0.16——風洞模擬流場的平均風速和湍流度剖面如圖2所示,Z表示高度,V為物理模型高度Z處風速,Vg為物理模型頂部風速,Iu為湍流強度。試驗在0°~360°范圍內每 隔 10°取 一 個風向角,共 36個 風 向角。風洞內1.5m高度處試驗風速約為8.0m/s,每個風向角下采樣長度為24 000個數據點,采樣時間為60 s,采樣頻率為400 Hz。

表1 高頻天平試驗主要相似比Table 1 Main scale ratios for high frequency force balance test

圖1 鍋爐塔架風洞試驗模型Fig.1 Wind tunnel test model of the tower installed with a boiler

圖2 風洞模擬的B類大氣邊界層風速與湍流度剖面Fig.2 Profiles of mean wind speed and turbulence intensities for category B

為保證高頻測力天平試驗的有效性,要求風洞試驗物理模型質量盡量小,頻率足夠高,使得物理模型自振頻率遠遠高于風荷載譜密度的主要分布頻段范圍,從而基本可以避免物理模型在風洞試驗過程中的風激振動。由于該鍋爐塔架結構及其使用要求的特殊性,中間鍋爐需要安裝到位,與塔架一起共同承受風荷載,相應風洞模型的制作比較復雜,在材料的選用上也受到一定的限制。對完工后的風洞試驗物理模型進行動力標定,發現前兩階自振頻率都約為12 Hz,物理模型剛度偏小,其自振功率譜見圖3(在物理模型頂部施加X向激勵時,底部FX的功率譜密度曲線)。圖4給出了0°風向角下物理模型基底彎矩和扭矩的功率譜密度曲線。從圖4可見由于風激共振產生的明顯譜峰。該風洞試驗模型底部天平的測力數據除包含平均風壓產生的平均風力與脈動風壓產生的隨機脈動風力外,還包含模型振動產生的慣性力成分,所以需要對半剛性模型的高頻天平數據進行頻響修正。由ANSYS有限元軟件分析得到原型塔架第一階模態振型為整體Y向平動,頻率為0.256 Hz;第二階為整體X向平動,頻率為0.276 Hz;第三階為整體扭轉,頻率為0.296 Hz。具體前十階模態頻率及振型見表2。從圖4的底部風荷載譜曲線,可以看出前兩階X向、Y向自振頻率十分接近,在12 Hz處有兩個緊密相鄰的共振峰;第三階扭轉自振頻率為40 Hz;前三階共振引起的譜峰值較大,而高階振型引起的譜峰值較小,故高頻天平數據修正只考慮前三階振動的影響,而忽略高階振型的影響。

圖3 模型自振功率譜圖Fig.3 Power spectrum of free vibration of the model

圖4 0°風向角下實測底部彎矩、扭矩譜Fig.4 Bending moment and torque spectra of the model under 0°wind attack

表2 塔式鍋爐前10階模態分析結果Table 2 The first ten vibration modes

2 半剛性物理模型頻響修正

在進行原型結構風致動力分析計算前,需要將上述物理模型振動慣性力對底部風力數據的干擾影響消除。基于隨機振動理論,可在頻域內對半剛性模型底部天平實測所得的風力譜進行頻響修正,從而消除慣性力對風力譜的影響。

2.1 公式推導

風洞試驗物理模型坐標系統如圖5所示。

記高頻測力天平試驗所測得的底部彎矩和扭矩如下:MX、MY、T。基于質量凝聚原理,鍋爐塔架風洞試驗物理模型的運動方程可依據達朗貝爾原理[15]給出如下:

圖5 模型坐標系統和風向角定義Fig.5 Coordinate system and wind angle definition of the model

式中:Z為N×1維高度向量ZT= [h1h2… hN], N為所取物理模型質點數;[FXFYTZ]T為作用在物理模型各質點上的風力分量;M、C分別為3N×3N質量矩陣、阻尼矩陣;[ X Y θ]T為物理模型各質點X、Y向平動及扭轉位移向量。

假定原型塔架基本振型為線性,并引入模態修正因子[16],則有:

其中,φiP為定義的原型塔架第i階模態振型,為3N×1維向 量,記 原 型 塔 架 振 型 為 ΦP= [φ1Pφ2P… φnP];Cix、Ciy、Ciθ為原型塔架第i階振型在頂層沿X、Y及扭轉向的分量,可用ANSYS[13-14]等有限元軟件建立原型結構的有限元模型分析計算得到;ηix、ηiy、ηiθ為模態修正因子,可取ηix=ηiy=1,ηiθ= 0.7;H、R風別為建筑高度及建筑平面回轉半徑。

ΦM= [φ1Mφ2M… φnM]為物理模型振型,φiM為模型第i階模態;Q= [q1Mq2M… qnM]T為各階模態位移響應,將式(4)代入式(3)有:

對于剛性模型,風洞試驗過程中基本可以認為其靜止不動,故相應的模態加速度響應、模態速度響應可取為0,第i階模態風力目標值Fi可以用風洞試驗所測底部彎矩、扭矩表示的模態力估計值近似得到:

對于半剛性模型而言,模型會在風洞測力試驗過程中出現明顯的振動,從而需考慮模型本身模態響應、的影響。此時將φiP表示為:

其 中 ΦP= ΦM[A1A2… An],A = [A1A2… An]為 振 型 轉 換 矩 陣,Ai= [a1ia2i… ani]T。將式(7)帶入式(5),并做功率譜變換有:

結構加速度是速度對時間的一階導數。從相平面上來看,加速度與速度矢量時間的相位角為90°,二者是完全不相關的,故速度和加速度時間的相關函數為0,所以對應的互功率譜耦合項可以完全忽略。此外,由于風荷載與速度、加速度的耦合項在整個功率譜能量中的貢獻較小,故也同樣進行省略的近似處理。故可以得到:

由于ΦP=ΦM[A1A2… An],取B=A-1,記為振型轉換矩陣逆矩陣,則物理模型振型可由塔架原型振型表示為:

則φiM=ΦPBi。那么對于物理模型其對應廣義模態力可表示為:

將式(11)兩邊做功率譜變換,并省略互譜項得物理模型的廣義模態力譜:

式中ω為圓頻率。基于隨機振動理論得:

式中ζj、ωj為物理模型第j階阻尼比、頻率。

這樣物理模型的廣義模態力譜Sqj(ω)可表示為實際各階模態力功率譜之和。類似的物理模型的模態速度譜、模態加速度譜也可由SFi(ω)表達。將其各自表達式代入式(9)中,得為:

若考慮結構n階振型修正,則有:

若只考慮第一階振型由式(16)得:

式(17)即與鄒良浩等人推導的公式[12]一致。

2.2 模態風力修正及風振分析結果

對風洞試驗數據采用式(16)進行修正,只考慮前三階振型,即取n=3。由于該塔式鍋爐為方形斷面,前兩階X向、Y向振型存在耦合作用,對物理模型需假設一種耦合模式,計算振型轉換矩陣A及其逆矩陣B。

假設物理模型振型為:

一階 二階 三階

根據試驗過程中觀測到的物理模型X、Y向的耦合振動模式,初略取a=、b=1。由有限元計算軟件可知該塔式鍋爐模型第一階振型為Y向平動,與X向平動及扭轉無耦合作用,根據式(7)可得:

可求得a11、a21、a31,同理求出A2、A3、B1、B2、B3,得:

帶入式(16)并只考慮前三階振型得:

由于規范對鍋爐類型的特殊結構有相關阻尼比的具體規定,考慮外部為鋼結構,內部為鍋爐體系的混合體系,取阻尼比為ζ=0.02。采用上述公式對塔式鍋爐風洞試驗數據進行處理,得到前三階模態風力譜如圖6所示。從圖中6可以看出,通過修正后,半剛性模型風激振動對前三階模態風力譜的共振峰影響得到了有效的消除。

圖6 修正后前三階模態風力譜圖Fig.6 Corrected model force spectra of three fundamental modes

2.3 原型結構風振計算結果

采用模態分析方法,將原型結構運動方程轉化為一組互不耦合的模態方程,每個模態方程可以獨立求解,而總的響應可以通過各個模態響應疊加得到。由隨機振動理論[17]得到原型塔架的模態位移響應譜如下:

式中CF為原型塔架風荷載位移譜相似常數,可由式(24)計算;λL=1∶150為模型幾何縮尺比例,λν為物理模型與原型塔架速度縮尺比,物理模型計算風速采用試驗1.5 m高度處的風速,平均值為8.0 m/s,并采用杭州地區百年一遇的基本風壓0.50 kPa折算,原型塔架設計風速為43.10 m/s,這樣風速縮尺比λν約為1∶5.4。

采用模態疊加方法,最后計算得到了原型結構在各風向角下的位移響應。圖7為各風向角下頂層位移變化曲線。從圖7可見,該鍋爐塔架在百年一遇風速條件下最大風振位移達40 cm,結構總的位移比為1/347,滿足高聳結構設計規范[18]對風荷載作用下結構水平變形的設計要求。

圖7 各風向角下頂層位移變化曲線Fig.7 Top displacements of the tower under each wind angle

3 結 論

本文首先從達朗貝爾原理出發建立運動方程,通過假定原型塔架振型與物理模型振型之間的轉換關系,運用模態疊加法,推導得出了半剛性模型天平測力數據的多階廣義風荷載譜修正公式。利用新得到的修正公式,發展了基于半剛性模型天平測力風洞試驗的鍋爐塔架風振分析方法。開展了某大型鍋爐塔架的半剛性模型高頻天平測力風洞試驗,采用考慮前三階模態影響的修正公式,對該大型塔架的高頻天平試驗底部三彎矩數據進行修正。修正達到所需效果:從修正后前三階模態風力譜圖上可知,由X和Y向耦合振動引起的譜峰得到了有效的消除。采用修正后風載荷數據進行鍋爐塔架的風振分析,計算得到結構總位移比為1/347,滿足高聳結構設計規范的要求。驗證了該超大型鍋爐塔架結構抗風設計的安全性和合理性。

本文提出和發展的復雜高聳結構半剛性模型風洞天平試驗和風振分析的方法,可為類似由于結構特殊及模型制作等原因形成的“半剛性”模型進行高頻天平風洞試驗及風振分析工作提供參照。

[1] Cook N J.A sensitive 6-component high-frequency-range balance for building aerodynamics[J].Phys.E.Sci.Instrum,1983,(116): 390-393.

[2] Tshanz T,Davenport A G.The base balance technique for the determination of dynamic wind loads[J].J.Wind Engineering&Industrial Aerodynamics,1983,13:429-439.

[3] Zhang Liangliang,Wu Yunfang,Yu Yang,et al.Experimental technology on dynamic wind loads on a tall building[J].Journal of Chongqing University,2006,29(2):99-102.

張亮亮,吳云芳,余洋,等.高層建筑風荷載高頻測力天平試驗技術[J].重慶大學學報,2006,29(2):99-102.

[4] Liang Shuguo,Zou Lianghao,Zhao Lin,et al.The investigation of 3-D dynamic wind loads on lattics towers by wind tunnel test[J].Acta Aerodynamica Sinica,2007,25(3):311-329.

梁樞果,鄒良浩,趙林,.格構式塔架三維動力風荷載的風洞試驗研究[J].空氣動力學學報,2007,25(3):311-329.

[5] Tang Yi,Gu Ming.Simulation method for spatial distribution of wind loads on tall buildings using HFFB data[C]//Shanghai International Industry Fair and the 3th Shanghai Engineering and Vibration Technology BBS,2007.

唐意,顧明.利用HFFB數據估計高層建筑多風荷載空間分布的方法[C]//上海市國際工業博覽會暨第三屆上海市“工程與振動”科技壇,2007.

[6] Yip D Y N,Flay R G J.A new force balance data analysis method for wind response predictions of tall buildings[J].Journal of Wind Engineering and Industrial Aerodynamics,1995,(54/55):457-71.

[7] Gu Ming,Zhou Yin,Zhang Feng,et al.Study on wind loads and wind-induced vibration of the Jin Mao building using high-frequencyforce-balance method[J].Journal of Building Structures.2000,21 (4):55-61.

顧明,周印,張鋒,項海帆,江歡成.用高頻動態天平方法研究金茂大廈的動力風荷載和風振響應[J].建筑結構學報,2000,21(4):55-61.

[8] Xie Zhuangning,Shi Biqing,Ni Zhenghua,et al.Experimental study on reduction of wind loads on the Shenzhen Kingkey Financial Tower by aerodynamic strategy[J].Journal of Building Structures.2010,31(10):1-7.

謝壯寧,石碧青,倪振華,等.深圳京基金融中心氣動抗風措施試驗研究[J].建筑結構學報,2010,31(10):1-7.

[9] Shi Biqing,Xie Zhuangning,Ni Zhenghua.Study on wind effects of Guangzhou West Tower using high-frequency-force-balance method[J].China Civil Engineering Journal,2008,41(2):42-48.

石碧青,謝壯寧,倪振華.用高頻底座力天平研究廣州西塔的風效應[J].土木工程學報,2008,41(2):42-48.

[10]Sun Yi,Li Zhengliang,Huang Hanjie,et al.Research on a super high-rise building based on the rigid model force measurement testing[J].Journal of Experiments in Fluid Mechanics,2010,24(4): 33-38.

孫毅,李正良,黃漢杰,等.某超高層建筑剛性模型測力試驗研究[J].實驗流體力學,2010,24(4):33-38.

[11]Cheng Feng,Hu Guofeng.Design and application of high frequency dynamic balance[J].Aerodynamic Experiment and Measurement&Control,1995,9(1):58-61.

程豐,胡國風.高頻底座天平的研制與應用[J].氣動實驗與測量控制,1995,9(1):58-61.

[12]Zou Lianghao,Liang Shuguo.A method to evaluate wind force spectra of semi-rigid model in wind tunnel tests[J].Journal of Experiments in Fluid Mechanics,2007,21(3):76-81.

鄒良浩,梁樞果.半剛性模型風洞試驗荷載譜的處理方法[J].實驗流體力學,2007,21(3):76-81.

[13]WangXinming.Numerical analysis of ANSYS engineering structure[M].Beijing:People’s trsffic press,2007.

王新敏.ANSYS工程結構數值分析[M].北京:人民交通出版社,2007.

[14]Yang Yongyi,Liao Haili,Li Yongle.Simplified buffeting analysis in time domain on the basis of ANSYS for cable-stayed bridge[J].Acta Aerodynamica Sinica,2004,22(4):457-460.

楊詠漪,廖海黎,李永樂.基于ANSYS的斜拉橋抖振時域實用分析方法[J].空氣動力學學報,2004,22(4):457-460.

[15]Li Xiqin.D’Alembert’s principle[J].Yunan Education,1980,12: 41-42.

李西秦.達朗貝爾原理[J].云南教育,1980,12:41-42.

[16]Huang Mingfeng,Tse K T,Chan C M,et al.Mode shape linearization and correction in coupled dynamic analysis of wind-excited tall buildings[J].The Structural Design of Tall and Special Buildings,2010,20:327-348.

[17] Cheng Daozhang.Random vibration theory and applications[J].Journal of UEST of China,1993,22(4):428-433.

陳道章.隨機振動理論及其應用[J].電子科技大學學報,1993,22(4):428-433.

[18]GB 50135-2006.Code for design of high-rising structures[S].The Ministry of Construction of the People’s Republic of China,Beijing:China Building Industry Press,2006.

GB 50135-2006.高聳結構設計規范[S].中華人民共和國建設部.北京:中國建筑工業出版社,2006.

Wind-induced vibration analysis of lattice-truss tower installed with a boiler based on semi-rigid model test

Wu Chenghui1,Huang Mingfeng1,*,Jiang Xiong1,Chen Yongjun2,Qiu Jian2
(1.College of Civil Engineering and Architechture,Zhejiang University,Hangzhou Zhejiang 310058,China; 2.Dongfang Boiler Group Co.,Ltd,Chengdu Sichuan 611731,China)

The model used in high-frequency-force-balance test must be light and rigid.However,for lattice-truss structures,it is difficult to ensure rigidity of their wind tunnel models due to structural features and model fabrication limits.Force data of a semi-rigid model measured by a force balance needs to be modified before applying them in vibration analysis of the prototype structure.Based on D'Alembert' s principle and mode superposition method,assuming the transformation relation between the prototype mode and the modal mode,a new formula for modifying modal force spectra of a semi-rigid model was obtained.With the aid of the newly derived formula,the traditional frequency domain analysis method can be used to predict wind-induced response of a lattice-truss tower structure.The wind tunnel test was conducted for a large-scale tower installed with a boiler.After applying the new formula,modal force spectra measured in the wind tunnel was corrected.Wind-induced response analysis was performed successfully for this tower structure.The results validated the adequacy of structural design of this tower in resisting wind storm.

lattice-truss tower;wind tunnel tests;semi-rigid model;data correction;modal wind force;wind-induced vibration

V211.7;TU312+.1

A

10.7638/kqdlxxb-2013.0041

0258-1825(2015)03-0353-07

2013-04-01;

2013-05-18

國家自然科學基金(51008275);浙江省公益性技術應用研究計劃(2012C21059);中央高校基本科研業務費專項資金(2012QNA4013)

吳承卉(1991-),女,浙江人,碩士生,研究方向:結構風工程.E-mail:sangyu1123@163.com.

黃銘楓*(1976-),男,浙江人,副教授,主要從事結構風工程研究.E-mail:mfhuang@zju.edu.cn

吳承卉,黃銘楓,姜雄,等.基于半剛性模型風洞試驗的鍋爐塔架風振分析[J].空氣動力學學報,2015,33(3):353-359.

10.7638/kqdlxxb-2013.0041 Wu C H,Huang M F,Jiang X,et al.Wind-induced vibration analysis of lattice-truss tower installed with a boiler based on semi-rigid model test[J].Acta Aerodynamica Sinica,2015,33(3):353-359.

猜你喜歡
模態模型
一半模型
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
車輛CAE分析中自由模態和約束模態的應用與對比
國內多模態教學研究回顧與展望
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
高速顫振模型設計中顫振主要模態的判斷
航空學報(2015年4期)2015-05-07 06:43:35
基于HHT和Prony算法的電力系統低頻振蕩模態識別
由單個模態構造對稱簡支梁的抗彎剛度
計算物理(2014年2期)2014-03-11 17:01:39
主站蜘蛛池模板: 少妇露出福利视频| 国产人在线成免费视频| 亚洲一区二区日韩欧美gif| 色丁丁毛片在线观看| 天天色综网| 72种姿势欧美久久久久大黄蕉| 国产区免费| 综合色天天| 91青青草视频| 日韩不卡高清视频| 亚洲人成色在线观看| av一区二区三区高清久久| 一本大道香蕉久中文在线播放| 国产一区二区三区日韩精品| 久久久久久国产精品mv| 欧美啪啪一区| 精品剧情v国产在线观看| 在线观看国产黄色| 91色在线观看| 国产91九色在线播放| 亚洲视频色图| 久久一级电影| 国产毛片网站| 午夜爽爽视频| 亚洲欧洲自拍拍偷午夜色| a毛片免费看| A级毛片无码久久精品免费| 伊在人亚洲香蕉精品播放| 夜夜操天天摸| 亚洲综合色区在线播放2019 | 久久综合国产乱子免费| 九九热视频精品在线| 草草线在成年免费视频2| 一级毛片免费的| swag国产精品| 国产精品亚洲精品爽爽| 免费国产在线精品一区| 99久久这里只精品麻豆| 91美女在线| 国产综合亚洲欧洲区精品无码| 亚洲大尺码专区影院| 中文字幕久久亚洲一区| 欧美综合区自拍亚洲综合天堂 | 在线中文字幕网| 欧美色丁香| 永久免费无码成人网站| 亚洲天堂.com| 色香蕉网站| 精品人妻AV区| 成人精品免费视频| 九色最新网址| 久久伊人久久亚洲综合| 国产三级国产精品国产普男人| 国产在线视频欧美亚综合| 欧美日韩亚洲综合在线观看| 97精品国产高清久久久久蜜芽| 亚洲人成网站在线播放2019| 在线播放国产99re| 亚洲成av人无码综合在线观看| 日韩精品一区二区三区中文无码 | 久视频免费精品6| 国产在线自乱拍播放| 国产在线观看高清不卡| 中文字幕乱码二三区免费| 亚洲精品大秀视频| 全裸无码专区| 男女性色大片免费网站| 久久香蕉欧美精品| 亚洲高清在线播放| 成年A级毛片| 亚洲另类色| 欧美色99| 免费看美女毛片| 国产超薄肉色丝袜网站| 一级毛片在线播放免费观看| 久久一级电影| 国产日本一线在线观看免费| 亚洲视频四区| 亚洲国语自产一区第二页| 五月天天天色| 国产成人无码AV在线播放动漫| 97国产精品视频人人做人人爱|