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

冷卻塔風振效應特征及影響因素分析

2018-09-28 02:27:50張軍鋒楊軍輝葛耀君
振動與沖擊 2018年18期
關鍵詞:風速

張軍鋒, 楊軍輝, 葛耀君, 趙 林

(1. 鄭州大學 土木工程學院,鄭州 450001;2. 同濟大學 土木工程防災國家重點實驗室,上海 200092)

風荷載是冷卻塔的設計控制荷載,其動力作用作為冷卻塔結構的關鍵問題之一,長期以來一直受到設計和研究人員的關注[1-2]。作為高聳空間薄壁結構,冷卻塔風振效應的影響因素眾多,包括紊流度、風速取值V、結構尺度L、結構基頻fmin以及阻尼比ζ等,但大多研究都是在特定參數下分析冷卻塔的風振效應[3-10],并未系統分析上述參數的影響。近年來,隨著我國冷卻塔建設高度的持續增加且壁厚相對更薄,結構愈加輕柔,冷卻塔對風荷載的動力作用將更加敏感,亟需明確其風致動力響應特征及影響因素。

實際上,已有學者對上述因素的影響展開研究,但分析方法、參數取值及具體結果等都仍有較大差異。文獻[11]僅分析了阻尼比ζ對共振分量的影響,但其計算風速高達100 m/s,現實中并不會出現。文獻[12]通過對三座冷卻塔的計算分析了結構基頻fmin對風振效應的影響,但其所得結果實際上還同時受結構尺度L的影響。另外,文獻[12]均選擇位移作為參考指標,但冷卻塔結構設計的控制指標是內力而非位移,且內力和位移的風振效應是不同的。文獻[13-14]雖考慮了多個因素,但每種因素僅考慮2~3個參數值,無法得到系統的變化規律。文獻[15-17]分別經試驗和計算分析,給出了綜合考慮V,L和fmin等多因素的表達式,并分別被英國和德國規范采納,但具體結果和應用形式亦不相同[18]。而我國規范[19]僅區分場地類型也即紊流度給出了風振系數取值,對其他因素并未涉及。這都說明對冷卻塔風振效應的影響因素還未有系統和明確的認識。

由于新建冷卻塔幾乎都位于城鎮郊區,也即我國規范中的B類地貌,其紊流度基本是確定的;結構尺度L則通過影響結構動力特性和風荷載表面相關性而最終影響風振效應,這三種因素實際難以獨立分析。因此本文以待建某內陸核電廠大型冷卻塔為例,通過剛體模型風洞試驗和結構動力響應時程計算,在前期風致動力響應特征的基礎上,通過調整阻尼比ζ、風速取值V和風壓子午向相關性,分析這三個因素對風振效應的影響。

1 工程背景與計算方法

本研究中的冷卻塔(見圖1)處B類場地,10 m高度基本風速V0=26.7 m/s,對應塔頂風速VH=40 m/s。對塔筒表面位置的描述區分環向和子午向兩個維度,環向角度θ以迎風點為0°,逆時針轉動為正;子午向以0≤hS/HS≤1表示塔筒相對高度位置。

試驗在同濟大學TJ-3風洞中進行,模擬B類地貌,幾何縮尺比λL=1∶ 200,采用適當措施進行雷諾數效應補償。塔頂試驗風速10 m/s,即風速比λV=1∶ 4,并由此可得時間比λt=1∶ 50。信號采樣頻率312.5 Hz,采樣時長19.2 s,每個測點6 000個數據。限于設備采集能力,僅測面A的9排斷面同步測量并用于時程計算,圖2給出了側面A斷面的平均風壓、脈動風壓以及依Davenport理論所得峰值因子g,圖2中各斷面的參考風壓均為本斷面高度的來流速壓。試驗及結果的詳細介紹可見文獻[20]。

圖2 不同斷面的環向風壓參數 Fig.2 Wind pressure parameters of different sections

采用ANSYS程序進行結構計算,塔筒采用Shell188單元模擬并劃分72×35(環向×子午向)個單元,下支柱和檐口采用Beam188單元模擬。結構前100階頻率分布在0.763~2.903 Hz,且各階頻率基本隨階數線性增加。此頻段范圍內的振型可分為三種[21](見圖3),即雙向諧波耦合振型、側彎振型和豎向伸縮振型,并以第一種為主導:三種振型分別有95個、4個和1個。動力計算采用直接積分法進行,風荷載時程采用風洞試驗結果,根據風速比λV和時間比λT換算到原型結構,并通過POD(Proper Orthogonal Decomposition)分解和預測,將9×36=324個風壓測點擴展到36×72=2 592個有限元節點;時間積分步長Δt=1/37.5 s,涵蓋前53階模態(f53=1.842 Hz),且即使將Δt降低一半,脈動響應的增幅亦不足2%;采用瑞利阻尼,各階模態阻尼比ζ約為1%;對于結構響應,直接以時程最值作為極值。除采用上述直接積分法進行計算外,還采用模態疊加法進行了校核計算,除FX外,兩者所得脈動響應誤差一般不超過1%。

圖3 結構典型模態Fig.3 Characteristic modes of the HCT

眾所周知,風對結構的動力效應體現在兩個方面,即自然風本身的脈動性以及其引發的結構慣性效應(包含阻尼效應)。結構風振效應在最初較長的時間里以頻域計算為主,并且脈動響應σ也相應地被近似劃分為背景分量σB和共振分量σR[22]:式(1)~式(3)以單自由度結構的位移響應為例給出了其計算表達式。其中k,f0為系統的剛度和頻率,S(f)為荷載P(t)的頻譜。

(1)

(2)

(3)

類似地,背景和共振響應在時域內也可進行劃分。式(4)以單自由度為例給出了時域內的動力學方程,其中m,c分別為系統的質量和阻尼,且c=2m2πf0ζ。如果忽略系統的慣性效應,則式(4)退化為式(5)。顯然,式(5)中的響應yB(t)完全是由自然風本身的脈動性引起的,也是精確意義的背景響應,并可直接通過擬靜力計算求得。此時將y(t)與yB(t)相減(見式(6)),則可得共振分量yR(t)。顯然,σ,σB和σR分別與y(t),yB(t)和yR(t)對應,下文對背景和共振分量的劃分也以時域方法為準。

(4)

kyB(t)=p(t)

(5)

yR(t)=y(t)-yB(t)

(6)

2 計算結果及分析

2.1 計算結果概述

圖4首先給出了基準工況下(ζ=1%且λT=1∶ 50)部分塔筒內力及位移的均值和陣風響應因子G沿塔筒高度的分布,包括子午向軸拉力FY,T,環向軸拉力FX,T,環向內凹和外凸彎矩MX,I、MX,O,內凹和外凸位移UI,UO,且均以幅值給出。可以看出,不同響應的均值響應和動力效應差別極為顯著,對其詳細分析參見張軍鋒等的研究。總的來說,除塔筒頂部0.2HS高度范圍,雙向軸力和位移的風振效應均較小而MX的風振效應顯著,在塔筒頂部,各內力均值響應下降而脈動響應提高,故風振效應均甚為顯著。需要注意的是,眾多文獻對冷卻塔的動力計算僅關注其FY,T,忽略了其環向受力,但計算表明,環向彎矩MX的G值明顯大于FY,T,所閱文獻中,只有文獻[23]在未考慮共振效應的計算中關注了MX且與張軍鋒等的研究有相同的結論。冷卻塔風荷載作用下的結構設計控制指標則有塔筒中下部的FY,T、塔筒上端局部區域的FX,T以及環向彎矩MX,下文則通過調整計算參數分析對上述響應的影響。

為進一步分析脈動響應特征,表1給出了ζ=1%時4個特征位置(A點為hS/HS=0.19,θ=0°;B點為hS/HS=0.49,θ=0°;C點為hS/HS=0.70,θ=-95°;D點為hS/HS=0.92,θ=-95°)內力的平均和脈動響應。其中:A點和B點位于塔筒中下部的迎風子午線,主要關注其FY,T;C點和D點位于塔筒中上部的側風區,主要關注其環向軸拉力FX,T;同時還關注4個位置的環向彎矩MX。

從表1可知,ζ=1%時,雙向軸拉力的總脈動響應σ都遠小于其均值E,σ/E僅為1/5.8~1/7.6,且FX,T的σ/E,或者說FX,T的脈動效應,略大于FY,T的σ/E;MX的E和σ則在同一量級,甚至在C、D兩點各個阻尼比下都有E<σ,這也說明MX的脈動效應較雙向軸力更為顯著,這也是圖4中彎矩的G極為顯著而雙向軸拉力的G相對偏小的原因。在脈動響應中,雙向軸拉力的σR也都明顯小于其σB,背景分量在總脈動響應中的比例σB/σ為85%~95%,且亦有FX,T的σB/σ略小于FY,T的σB/σ,說明FX,T的共振分量權重較FY,T略大;而MX的σR則和σB接近,甚至有σR>σB,說明共振分量在MX的脈動效應中占有較大的權重。也正因為MX的共振分量權重較大,本文MX的G值為3~6,而Abu-Sitta等的研究在忽略共振分量時所得MX的G值為2~3。

圖4 各響應的均值和陣風響應因子GFig.4 Mean responses and gust response factors

響應位置均值EσBζ=1%σRσζ=3%σRσζ=5%σRσFY,T/(kN·m-1)A471.3558.8319.6063.4411.5861.439.1061.00B346.2343.8517.1947.539.7545.367.3844.91FX,T/(kN·m-1)C37.835.902.436.451.536.161.246.09D49.307.023.818.132.457.582.027.44MX/(kN·m·m-1)A-1.440.931.311.580.751.160.561.05B-5.071.541.822.371.021.830.771.70C0.671.311.181.750.741.490.611.43D0.751.000.961.430.651.230.561.18

上述3個內力脈動效應(σ/E)的大小以及共振分量在總脈動分量中權重(σR/σ)的大小關系,實際上也是冷卻塔結構受力特性的表現[24-25]:MX屬于局部內力,主要受其附近位置風壓的影響,而FY則屬整體內力,與整個塔筒表面的風壓有關,顯然,大面積風壓的疊加組合使其脈動性下降。

另外,從圖9可知,如忽略共振分量,對于FY,hS/HS<0.5,0.5

2.2 阻尼比ζ的影響

圖5 脈動響應隨ζ的相對變化Fig.5 Gust responses for different damping ratios

另外,鑒于前文對背景和共振分量的劃分,兩者之間將存在一定的耦合分量。如果依然根據σB和σR以式(1)近似計算總脈動響應σ,表1中各內力σ的誤差均小于5%,說明實際的耦合分量極小。這也是因為冷卻塔表面的風荷載頻譜能量主要集中在0.5 Hz以下(見圖6),而結構基頻為0.763 Hz,所以背景和共振分量之間幾乎沒有耦合性(見圖7):在低于0.5 Hz的頻段,總脈動響應的功率普曲線與背景響應重合;而在高于0.7 Hz的頻段,總脈動響應的功率譜曲線則與共振響應重合,而在0.5~0.7 Hz的背景和共振分量雖有一定的耦合效應但幅值極小。同時,各分量的頻譜分布也是對表1脈動響應數值的一個驗證。另外,對比各內力共振分量的功率譜分布可知,同一振型對不同內力共振響應的貢獻并不一致:圖7中3個內力共振分量的貢獻分別主要來自0.7~0.8 Hz,0.7~0.9 Hz且1.2~1.4 Hz以及0.8~1.0 Hz之間的振型。并且從圖7還可知,各響應共振響應都主要分布在0.7~1.6 Hz,與本計算的截取頻段(<1.84 Hz)一致,也說明積分步長Δt取值合理。

圖6 脈動風風速頻譜Fig.6 Davenport wind velocity spectrum

圖7 特征位置內力響應的功率譜密度Fig.7 PSD of three internal forces

需要說明,由于FY的共振分量相對極小,所以在圖7(b)以對數坐標給出以便于辨識。圖7(d)還給出了模態疊加法所得的FX頻譜,與完全積分法所得結果(見圖7(a))對比可知,前者總脈動響應的頻譜在低頻段偏小,但其共振分量在0~0.5 Hz仍有分布,而這顯然是不合理的。這一現象也使模態疊加法所得σ小于完全積分法的結果,而σR大于完全積分法的結果。模態疊加法出現這種結果的具體原因有待進一步分析,所以本文所給均為完全積分法的結果。

2.3 風速取值V的影響

盡管本冷卻塔設計風速下因風譜能量集中的區段低于結構基頻,共振效應較弱,但如果設計風速增加,脈動風的頻譜漸趨高頻(見圖6),結構共振效應σR在總脈動響應σ中的權重也將增加:因為σR不僅與荷載幅值有關,還與荷載頻譜分布和結構動力特性有關。為得到不同風速下的響應,可在保持λL=1∶ 200不變的情況下對λV=1∶ 2.22~1∶ 13.33和λT=1∶ 90~1∶ 15進行調整,以分析風速變化對共振和總脈動響應σR和σ的影響。但這樣實際上同時調整了荷載幅值和頻譜分布,不便于獨立分析荷載頻譜分布的影響。

眾所周知,對于線彈性結構的動力響應,如果荷載p(t)的時頻特征不變而僅是幅值變化,則各響應的時頻特征亦保持不變,響應的各個特征值如E,σB,σR和σ均正比于荷載幅值,各響應的峰值因子g和最終的陣風響應因子G亦保持不變。所以,對λV=1∶ 13.33 &λT=1∶ 15的計算,可將其時程結果乘以系數(15/50)2,以屏蔽荷載幅值的影響而僅計入頻譜變化的影響。這也相當于在保持λV=1∶ 4不變的情況下僅調整λT=1∶ 90~1∶ 15來分析荷載頻譜變化對共振和總脈動響應σR和σ的影響,即將原型結構風壓時程在保持幅值不變的情況下,將時間坐標軸從960 s(λT=1∶ 50)逐漸拉伸至1 728 s(λT=1∶ 90)或壓縮至288 s(λT=1∶ 15),以分別得到相當于V0=14.8 m/s和V0=89 m/s的頻譜:顯然,此過程中的平均和背景響應保持不變。另外,我國規范的V0取值絕大多數在33.4 m/s(即λV=1∶ 5 &λT=1∶ 40)以下,過高V0下的結構響應計算并無實際意義,因此可僅關注響應隨荷載頻譜的變化而忽略荷載幅值的影響。

圖8給出了4個位置內力共振和總脈動響應σR和σ隨風壓頻譜的變化,其橫坐標僅表示此風速下的頻譜而不計其幅值。由圖8可知,脈動響應僅在V>22.3 m/s的區段有所變化,說明在V<22.3 m/s時,風壓頻譜在0.763 Hz以上區段不再隨風速變化,故脈動響應趨于穩定。在V>22.3 m/s時,隨著V的增加,因荷載頻譜趨向高頻,各內力的σR都顯著增加,但增幅并不完全一致,這也說明同一振型對不同內力共振響應的貢獻并不一致。由于σR對σ的貢獻有限,所以λT增加時,σ的增加相對較小,并且依然以雙向彎矩的增加最為顯著,FX其次,FY最小。

另外,即使是現階段高200 m的冷卻塔,其基頻也一般在0.7 Hz以上,而設計基本風速V0多在30 m/s以下,由此及風速頻譜分布可知,冷卻塔的共振分量在總脈動響應中的比重依然很小。

需要指出,在V>22.3 m/s時,隨著風速V增加所帶來的風壓頻譜的變化,各內力的σR基本隨風速V呈線性增加。這也是因為,在0.763 Hz以上的區段,每個頻率的能量也基本隨風速線性增加(見圖6)。如果再考慮風速增加所帶來的風壓幅值的增加,則各內力的σR將與V3成正比:本次計算中,塔筒4個高度位置(hS/HS=0.19,0.35,0.49,0.70)迎風點的FY的σR與V2.85~3.25成正比,與預期相符。然而,文獻[15-17]分別經計算、試驗和計算指出,對于塔筒中下部的FY,其σR分別與V3.4~3.75,V4,和V3.35成正比。其中,Armitt等的研究因模型基頻偏小且計入了群塔干擾效應,故其所得σR最為顯著,風速V的冪指數也最大;Hashish等和Niemann的研究只是簡單給出了σR與V的關系,并沒有詳細的計算結果。相比之下,本文所得冪指數雖略有偏小,但可從風譜能量隨風速的變化規律(見圖6)給予解釋。

2.4 風壓相關性的影響

部分學者根據實測或試驗風壓進行風振響應計算時,由于測點僅布置在一個斷面,所以不得不假定脈動風壓在子午向完全相關[27-28],或者采用假定的子午向相關系數進行計算,且完全相關與部分相關所得風振

響應亦有差異。又因風壓環向相關性的研究較多且結論較為一致,本節在保持環向相關性的基礎上,僅在子午向考慮完全相關和部分相關兩種情況進行分析。

圖8 脈動響應隨時間比λT的變化Fig.8 Gust responses for different λT

前述計算采用了整個塔筒表面的風壓,已計入了風壓在子午向實際的部分相關性。為實現子午向的完全相關,對試驗數據進行如下調整:選擇某一斷面z0的脈動風壓時程為基準,其環向θ位置測點的脈動風壓時程為cP(z0,θ,t),則任意位置(z,θ)的脈動風壓時程為

(7)

式中:σP(z,θ)和σP(z0,θ)分別為所求位置及參考位置的脈動風壓根方差。如此,既可保留各斷面測點風壓的根方差幅值,又可使不同斷面之間同一θ位置的測點風壓時程完全相關。由于塔筒表面各測點的平均風壓并未改變,結構的均值響應E亦保持不變。分別取多個斷面的風壓作為基準進行調整及計算,圖9給出了幾種風壓所得內力的根方差。另外,為了說明子午向相關性對脈動響應σ的影響,還給出了原始全斷面風壓和以斷面9為基準所得全相關風壓的背景分量σB。

圖9 子午向相關性對脈動響應的影響Fig.9 Influences of meridian correlation on gust responses

由于保留了各斷面脈動風壓的根方差幅值,所以選取不同基準斷面所得內力根方差基本一致,尤其是對雙向軸力;只是不同基準斷面內的風壓頻譜及斷面內的相關性存在一定差異,所以各自所得內力脈動響應僅略有差異,下文分析僅針對以斷面9為基準所得全相關的結果。

與全斷面風壓結果相比,子午向全相關的風壓使FY的脈動響應σ在各高度位置都有較為明顯的增加,尤其在塔筒中部區段:塔筒中下區段(hS/HS<0.7)σ的增幅在10%~50%,且增幅沿高度的增加而增加;上部區段σ的增幅又沿高度急劇回落至零。這是因為FY的脈動響應以背景分量σB為主,并且從內力影響線(見圖10)可知,FY受整個塔筒高度范圍內荷載的影響,所以當采用子午向全相關風壓時,其背景響應σB有明顯增加,從而使總脈動響應σ亦明顯增加。

圖10 B點位置的內力影響線Fig.10 Influence lines of internal forces at position B

脈動響應σ的增加也使FY的G值有顯著增加,尤其是在中部區段:在hS/HS=0.3和0.7兩個位置,G值分別從部分相關的1.48和1.73提高到約1.57和2.12。考慮到自重在子午向產生的壓力,風荷載子午向拉力的些許增加就將使荷載效應組合中的凈拉力顯著增加,因此,假定子午向風壓完全相關時,FY的脈動響應和最終配筋都將會有顯著增加。

對于環向內力,采用子午向全相關風壓時,塔筒不同高度位置FX和MX的脈動響應σ各有增減,且增減幅值并不顯著,塔筒中下區段FX的σ甚至下降約25%。這是由于FX的脈動響應仍以σB為主,但其主要受相鄰高度區域內荷載的影響,而此區域內的影響線存在反號,所以FX的σB和σ隨子午向相關性的增加而減小。MX依然主要受相鄰高度區域內荷載的影響且影響線符號恒定,而區域外的影響線幅值雖小但符號相反,所以MX的σB隨子午向相關性的增加無明顯變化,僅在中部區域下降約10%。另外,對于MX的σ,其受子午向相關性的影響亦不明顯。

需要說明,限于篇幅,上述內力σB與子午向相關性的關系僅是借助θ=0°子午線上的荷載與B點位置內力影響線的定性分析,實際上應根據影響面進行分析。對于結構設計來說,盡管MX的σ在塔筒下部區段有明顯增加且使此區域MX的陣風響應因子G也增加約10%,但考慮到FX的σ有所下降,且風荷載的環向內力在多荷載效應組合中的權重有限,其對結構配筋設計的影響也可以忽略。

3 結 論

以風洞試驗所得風壓為基礎,通過風振計算分析了塔筒各響應的風振特征,并探討了阻尼比ζ、風速取值V和風壓子午向相關性這3個參數對風振效應的影響,主結論要如下:

(1)除塔筒頂部0.2HS高度范圍,雙向軸力和位移的風振效應均較小而MX的風振效應顯著,在塔筒頂部,各內力的均值響應下降而脈動響應提高,故風振效應均甚為顯著。在總脈動響應σ中,各內力背景分類σB和共振分量σR的耦合性均極小,可以忽略。其中,雙向軸力的σR均明顯小于其σB,σB/σ為85% ~95%,即脈動響應以背景分量為主,而MX的σR則和σB接近,甚至有σR>σB,即共振分量貢獻顯著。

(2)各內力共振分量σR對總脈動響應σ的貢獻隨高度位置變化。忽略σR時,在塔筒中下區段,雙向軸力的σ下降不超過10%,而MX的σ下降約20%;在塔筒上部,各內力的σ均下降約20%~30%。對于結構配筋設計來說,FY的σR可以忽略而塔筒中上區段環向內力的σR不可忽略。

(3)隨阻尼比ζ的增加,各響應的共振分量σR基本按統一規律顯著下降,且近似符合式(3)的量化關系。但由于各內力σR對σ的貢獻不同,所以ζ增加至5%時,雙向軸力的σ僅下降約5%~7%,而MX的σ下降約25%。

(4)共振分量σR隨風速V的增加有兩方面的原因,荷載頻譜趨向高頻而導致的σR的增加,以及荷載幅值增加使各響應分量的增加。僅計第一個原因時,各響應σR基本隨V(或時間比λT)呈線性增加,所以同時考慮兩方面的原因時σR與V3成正比。

(5)子午向相關性主要是通過σB而影響σ,并且主要對雙向軸力的σ有影響。這也是由于雙向軸力的σ均以σB為主,子午向相關性的提高使FY的σ有所增加但使FX的σ有所增減小,這可通過內力影響面得以驗證。對于環向彎矩MX,由于其σB相對偏小,且受子午向相關性影響不大,所以子午向相關性對MX的σ影響亦不明顯。

猜你喜歡
風速
邯鄲市近46年風向風速特征分析
基于Kmeans-VMD-LSTM的短期風速預測
基于最優TS評分和頻率匹配的江蘇近海風速訂正
海洋通報(2020年5期)2021-01-14 09:26:54
基于時間相關性的風速威布爾分布優化方法
陜西黃土高原地區日極大風速的統計推算方法
陜西氣象(2020年2期)2020-06-08 00:54:38
基于GARCH的短時風速預測方法
快速評估風電場50年一遇最大風速的算法
風能(2016年11期)2016-03-04 05:24:00
考慮風切和塔影效應的風力機風速模型
電測與儀表(2015年8期)2015-04-09 11:50:06
GE在中國發布2.3-116低風速智能風機
考慮風速分布與日非平穩性的風速數據預處理方法研究
主站蜘蛛池模板: 亚洲91在线精品| 大香伊人久久| 动漫精品啪啪一区二区三区| 91在线播放免费不卡无毒| 日本妇乱子伦视频| 久久semm亚洲国产| 久久久久夜色精品波多野结衣| 国产免费羞羞视频| 91精品免费高清在线| 精品人妻AV区| 亚洲男人的天堂在线观看| a天堂视频| 亚洲综合精品香蕉久久网| 日本精品影院| 国产丝袜第一页| 亚洲综合第一页| 亚洲天堂区| 99热在线只有精品| 99热这里只有精品久久免费| 国产黄色视频综合| 亚洲国产成人久久77| 婷婷亚洲最大| 人妻精品全国免费视频| 人妻无码AⅤ中文字| 成人一区在线| 国产精品播放| 国产嫩草在线观看| 91精品专区| 女人毛片a级大学毛片免费| 中文天堂在线视频| 国产免费网址| 国产精品福利尤物youwu | 国产九九精品视频| 国产精品部在线观看| 三级毛片在线播放| 亚洲视频a| 国产免费看久久久| 亚洲天堂网视频| julia中文字幕久久亚洲| 色婷婷成人网| 伦伦影院精品一区| 国产国模一区二区三区四区| 一级毛片免费播放视频| 亚洲男人在线天堂| 国产精品大尺度尺度视频| 亚洲精品第一页不卡| 久久免费成人| 国产人成网线在线播放va| 人妻无码中文字幕一区二区三区| 欧美性久久久久| 国产精品13页| 91精品网站| 欧美激情福利| 最新加勒比隔壁人妻| 国产欧美又粗又猛又爽老| aa级毛片毛片免费观看久| 国产无遮挡猛进猛出免费软件| 欧美在线导航| 国产在线一区视频| 国产丰满大乳无码免费播放| 午夜激情婷婷| 国产精品无码一二三视频| 亚洲av无码久久无遮挡| 中文字幕在线视频免费| 亚洲女人在线| 五月天在线网站| 国产高清在线精品一区二区三区| 亚洲娇小与黑人巨大交| 国产视频久久久久| 天天综合色网| 美女潮喷出白浆在线观看视频| 国产美女在线观看| 99视频在线精品免费观看6| 四虎AV麻豆| 精品偷拍一区二区| 在线观看无码av免费不卡网站| 日韩a在线观看免费观看| 国产精品无码翘臀在线看纯欲| 亚洲国产精品久久久久秋霞影院| 精品视频一区在线观看| 久久久久亚洲Av片无码观看| 国产一区二区福利|