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

輸電線?絕緣子體系三維抖振響應的時/頻域理論方法研究

2022-11-14 01:08:24汪大海王濤汪偉徐康梁樞果
振動工程學報 2022年5期

汪大海 王濤 汪偉 徐康 梁樞果

摘要:強風作用下往往產生大變形的非線性動力抖振響應,是導致電力網絡風偏閃絡和支撐桿塔風災破壞的主要原因。為了揭示強風作用下輸電線路風效應形成機理,以典型兩跨輸電線?絕緣子為研究對象,以絕緣子的風偏位移和桿塔對絕緣子端部的空間支座反力響應為考察重點,基于索結構力學,給出了靜力平均風偏狀態的非線性解,并考察了體系在風偏狀態下模態、氣動阻尼等動力特性的變化;推導了響應的影響線函數及模態參與系數表達式,提出了線性時域動力響應的計算方法;進而依據風工程理論,給出了脈動風振響應背景分量和共振分量的頻域表達式;采用典型算例,通過與非線性有限元模型結果進行分析比較。結論表明,提出的三維振動時/頻域理論模型及參數計算取值方法具備足夠的工程效率和精度。研究為揭示輸電線路風災的破壞機理,完善輸電線路結構的抗風設計提供了重要的理論基礎和計算方法。

關鍵詞:非線性;抖振響應;輸電線;背景分量;共振分量

中圖分類號: O322;TU312+.1??? 文獻標志碼: A??? 文章編號:1004-4523(2022)05-1109-09

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

引言

高壓輸電線是典型的風荷載敏感結構,在各種風振響應中,強風作用下往往發生大變形的非線性動力抖振響應。這是導致電力網絡風偏閃絡和支撐桿塔風災破壞的主要原因。國內外學者從理論分析計算、風洞試驗和現場實測[1?7]等方面對大氣邊界層湍流作用下輸電線的抖振響應及其傳遞給桿塔的風振荷載進行研究,研究結果也被當前的國際設計規范和標準所采用[8]。

輸電線的風致響應可以分為靜風荷載作用下的靜力響應和靜平衡狀態處脈動風引起的動力響應。在較高的風速下,輸電線平面的靜力變形和風偏非常顯著,這也使得輸電線張力由初始狀態(只受重力)發生了相當大的變化,從而導致了系統動力特性的變化。這就需要一個用于靜態響應分析的非線性分析過程,而當張力變化可以被忽略時,線性理論更方便[9]。 Davenport[10]提出了基于傳統線性隨機振動理論的陣風響應因子方法。 Matheson 和 Holm ?es[11]通過使用有限差分法求解微分方程進行時域響應分析,與線性隨機振動理論進行比較,還基于準定常理論計算了考慮氣動阻尼效應的風致阻力。結果表明,輸電線的平均風偏對輸電線的豎向位移有明顯的影響,也證實了橫向支反力主要由背景分量貢獻。文章通過不考慮風偏角的線性隨機振動理論給出了較為準確的計算方法。Loredo ? Souza 和 Daven ?port[12]比較了使用基于擬靜力模型的線性隨機振動理論和氣彈風洞試驗得出的響應,結果證實了由于氣動阻尼較大,背景分量在橫向響應中占主導地位。輸電線的幾何非線性風振響應可用非線性有限元模型計算,也可用關于模態位移的耦合非線性運動方程進行分析,其動力響應可以由平面內和平面外模態位移來表示。由于動張力的作用,模態位移方程與各非線性項耦合,其中包含基于變形協調條件的二次項[13]。事實上,當靜力平衡狀態附近的動力變形較小時,由動力響應引起的動張力可以通過系統線性方程來確定,線性系統可以使用隨機振動的譜分析方法給出頻域動力響應解[6]。

綜上所述,目前的輸電線風振響應的理論方法主要集中在單跨鉸支輸電線的順風向振動方面。實際的輸電線路由輸電線與絕緣子連接而成。提出輸電線路三維抖振響應的理論模型是揭示輸電線路風災破壞機理的關鍵。本文基于索結構力學和風工程理論,以典型兩跨輸電線?絕緣子體系的絕緣子風偏位移以及端部對桿塔的空間支座反力響應為研究內容,推導并建立了體系強風抖振空間響應的時/頻域的理論計算方法。首先,推導給出了靜力平均風偏剛度和變形的非線性解,并考察了其對體系模態、氣動阻尼等動力特性的影響。進而在平均風偏狀態下,分別推導了兩種響應的影響線函數及模態參與系數,得到了時域線性動力響應理論計算方法,在此基礎上,給出了脈動風振響應背景分量和共振分量的頻域計算方法。最終,通過算例與非線性有限元模型結果分析進行比較,驗證了該理論方法的精確性和效率。研究為完善和修訂輸電線路的抗風設計提供了理論基礎。

1 靜風響應分析

考慮兩端固定在同一水平高度上的雙跨輸電線,中間通過長度為 l 的懸垂絕緣子與支座連接,垂跨比為 d0/L (d0為初始弧度)。初始狀態下,單位長度重力荷載為 mg=q0,水平張力為 H0,初始位置為 y0(x),垂度為 d0=q0L2/(8H0),如圖1( a )所示。不考慮輸電線的抗彎剛度、輸電線只承受拉力、輸電線橫截面積處處相等。本研究考慮風致響應最為不利的風向,即與輸電線垂直的方向。絕緣子風偏狀態示意圖如圖1(b)所示。單位長度上的平均風荷載由下式給出:

式中ρ為空氣的密度;為水平支撐以下 d0高度處的平均風速,其值看作沿輸電線跨度不變;CD 為阻力系數;D 為輸電線直徑。

帶絕緣子輸電線的微分方程與鉸支座單跨輸電線的靜力平衡方程相似。易求得靜風荷載作用下輸電線的豎向和橫向的靜風位移0( x ),0( x )分別為(其中,-L≤x≤L ):

式中? 0(0)= l (1- q0/q ),0(0)= D /q,分別為絕緣子端部豎向和橫向靜風位移,q =[ q 0(2)+fˉD2]0.5為風荷載與自重的合力,l 為絕緣子長度,L 為輸電線跨度。水平張力 H 由下面的非線性方程確定:

式中? E 為楊氏模量,A 為輸電線的橫截面積。新平衡狀態的 y( x )仍為拋物線,其弧垂為 d=qL2/(8H),輸電線平面的風偏角為θˉ= arctan ( fˉD /q0),絕緣子上端三個方向的反力分量為:

2 背景分量計算

背景分量可以直接由影響線函數確定。對于 t 時刻一個給定的脈動響應( t ),其背景分量 B 由下式給出:

式中μ( x )為影響線函數,表示當單位荷載作用在 x 處時引起 B 響應的增量。對于縱向脈動反力,其影響線函數可表示為:

其中,Irvine 常數λ2=(qL/H)2L/[HLe /(EA)];輸電線弧長 Le =L [1+8(d/L )2]≈L 。對應于橫向脈動反力的影響線函數為:

由于水平荷載作用在輸電線上不會引起支座處的豎向反力,因此,對于豎向反力的影響線函數縱向位移μu 或者橫向位移μw 的影響線函數可根據該方向的位移與反力的相似關系得到,可統一表示為:

依據公式(8)不難發現,由于橫向分力 T?z0( t )的影響線函數由平衡方程即可確定,因此,T?z0( t )的背景分量不受風偏狀態的影響。但是,絕緣子端部位移和縱向反力的影響線函數是水平張力 H 的函數,與風偏狀態下的剛度有關。

依據線性隨機振動理論,脈動風荷載作用下,背景分量的均方值可以計算為:

式中? IV =σV / V 為湍流強度,σV 和 S V ( f)分別為湍流的均方根和功率譜密度,f為頻率,coh ( x 1,x2,f )為脈動風相干函數??梢宰C明,直接使用影響線函數計算的背景分量與考慮所有背景模態的響應是等價的,因此影響線函數比模態分析方法效率更高,精度更好[6]。

3 共振分量計算

3.1? 模態位移響應

共振分量通過模態分析確定。當靜力平衡狀態附近的動力變形較小時,動力響應引起的動反力與靜反力相比也很小,動力方程可以基于準定常理論線性化,平面內和平面外模態位移的氣動阻尼耦合動力方程表示為[6]:

式中qiv (t)和qiw (t)分別表示面內和面外的模態位移;ξiv 和ξiw為結構模態阻尼比;固有頻率ωiv =2πfiv,ωiw =2πfiw,且fiv,fiw為輸電線自振頻率;Nv 和Nw分別為面內和面外模態截斷數。氣動阻尼比ξais1 和ξaijs1 s2( s1= v,w;s2= v,w )分別為:

其中,v 表示面內,w 表示面外,則有:

面內及面外模態的廣義力可分別表示為:

式中fD ( x,t)=ρDCD? V ( x,t )為脈動風荷載,V(x,t)是橫向脈動風速。φiv ( x )和φiw ( x )分別為平面內和平面外的振型。依據線性振動模態疊加理論,時域內,系統的任意響應均可采用模態位移和對應的模態參與系數ri ( x )表示為:

特別地,當所求 ( x,t )為位移響應時,ri ( x )即為模態φi ( x )。

3.2? 模態參與系數

模態參與系數可定義為模態慣性力下的靜態響應。引入響應的影響線函數有:

式中μv ( x ),μw ( x )為對應于該響應的面內及面外作用力影響線函數。考慮連續索結構的動力特征方程:

式中ΔHi 為水平附加張力。將公式(19)和(20)代入公式(17)和(18),可以得到:

對于本文重點考察的絕緣子端部支座處的反力響應,面內外反力T?y和T?z由模態位移和模態參與系數計算得到,而支座處縱向反力T?x可根據位移和反力的比例關系求出,如下所示:

表示x =0處 t 時刻的面外位移,動反力在 x0,y0,z0方向的分力為:

面內外反力T?y,T?z亦可以根據縱向張力定

3.3? 模態共振分量的頻域計算

基于隨機振動理論,對于非耦合模態振動,共振模態位移的均方根為:

其中,廣義力的功率譜密度為:

則有? vv = sin2θˉ,? ww = cos2θˉ。一旦確定了模態位移的均方根,可以由 SRSS 方法(總響應為背景分量與共振分量之和)計算( t )的共振分量的均方根。響應的極值為靜力響應加上峰值因子乘以響應的均方根[14]。依據現有研究,模態響應的背景分量與共振分量之間的相關性可以忽略不計[15].不難發現,由于風偏狀態下輸電線的頻率和氣動阻尼比發生了變化,共振分量也受到了靜力風偏的影響。

4 算例分析

4.1? 風場特性

考慮一根500 kV 的雙跨輸電線,它的參數如下:跨度 L=400 m,垂跨比 d0/L=1/30,線密度 m =2.39 kg/m,直徑 D=0.036 m,EA=48.8 MPa 。初始張力 H0=35.2 kN,Irvine 常數λ2=98.6。輸電線受到垂直于輸電線平面的風荷載作用,阻力系數CD =1.0。風速曲線為 =0( z/10)(1/6.5),0為10 m高度處平均風速,z 是地面以上高度,橫向脈動風的功率譜密度由Kaimal譜給出:

式中? u*為摩擦速度,σ V = u*。在參考高度即支撐水平以下的 d0處,z=70? d0=61 m 。

脈動風的相干函數由 Davenport 的指數函數模型給 c(出):

衰減因子為 C=16。除了上述模型之外,考慮到輸電線的幾何非線性,響應分析使用非線性有限元方法。

4.2? 靜風響應算例分析

圖2顯示了有限元和理論分析方法在不同風速下的垂度比 d/d0,水平張力比 H/H0,支座處反力比 y0/H0,面內外跨中位移比 w/d0和 v/d0,風偏角θ,結果表明,兩種分析方法中的靜力特性非常接近,驗證了靜風響應分析理論的正確性。從位移比和風速之間的斜率可以看出,由于幾何非線性的影響,輸電線的剛度隨著風速的增大而增大。風速較大時,靜張力隨著靜風荷載的增大而明顯增大。

只有當兩個振動模態頻率彼此非常接近時,平面內和平面外模態響應的耦合效應才很重要。事實上,反對稱的平面內和平面外的振動具有相同的模態頻率和形狀。一些對稱的平面內模態的頻率也可能非常接近具有對稱的平面外模態的頻率,需要考慮這些模態的耦合以及耦合氣動阻尼項的影響。雙跨輸電線?絕緣子體系的前幾階模態的示意圖如圖圖3所示,考慮到懸垂絕緣子軸向拉伸剛度較大,可忽略其軸向受拉變形,僅需考慮絕緣子對于面外對稱模態的影響。未耦合的模態響應可以看成單自由度系統進行計算。

4.3? 動力響應

采用譜表現方法,可以方便地模擬多個位置的平穩脈動風[16]。截止頻率fmax =5 Hz,頻率增量Δf=5/4096 Hz=0.0012 Hz 。脈動風持續時間為500 s,時間間隔為Δ t=0.125 s 。參考高度處的平均風速為 =30 m/s,湍流強度取為 IV =11%。

對于給定的脈動風時程樣本,采用有限元和模態分析方法計算風振響應時程。在模態分析方法中,包括在面內和面外方向上的前12階模態,即總共24階模態。基于耦合運動方程的 Newmark 方法確定模態分析響應,包括背景分量和共振分量。與氣動阻尼相比,結構模態阻尼很低,故可以忽略結構模態阻尼。首先消除了前100 s響應時程,即消除了

瞬態結構動力的影響。通過模態分析與非線性有限元分析結果進行了比較。圖 4和5分別給出了支座處位移時程和功率譜。圖 6和7給出了支座動反力響應時程及其功率譜。

除了進行以上模態響應時程分析外,下文還采用非線性有限元方法進行響應分析,并與響應譜理論分析模型的結果進行了比較。為了研究靜力平衡對動力響應的影響,基于理論公式的響應譜分析,其中背景分量使用影響線函數計算,共振分量通過模態分析確定,將是否考慮輸電線靜力風偏角的兩種計算結果與有限元方法比較。圖 8顯示了位移和反力分量的均方根隨風速的增加而增大,湍流強度 IV =11%。

可以發現,考慮輸電線平面靜力風偏的響應與非線性有限元分析結果吻合較好。由于氣動阻尼較大,動力響應主要來源于背景分量。如前所述,橫向動反力的背景分量不受輸電線平面靜力風偏的影響。然而如圖8( a )所示,當不考慮輸電線平面靜力風偏時,在支座處橫向位移的背景分量和總響應顯然被高估了。另外如圖8(b)和(d)所示,不考慮靜力風偏無法準確計算豎向位移和縱向反力分量。

對于上述響應,峰值因子范圍為2.68到4.03。如圖9所示,當 IV =11%時,隨著風速的增加,橫向動反力的陣風響應因子 GRF(Gust Response Fac? tor)從1.42略微增加到1.45。而橫向位移的 GRF 從1.42減小到1.15。由于脈動響應與湍流強度成正比,因此湍流強度越高,GRF 越大。

5 結論

本文以典型兩跨輸電線?絕緣子體系的絕緣子風偏位移以及端部對桿塔的空間支座反力響應為研究內容,將風效應分為平均風作用下的幾何非線性靜力效應和脈動風作用下的線性動力效應,基于懸索力學理論,通過推導體系的影響線函數和模態參與系數,分別計算了背景分量和共振分量。逐步建立了體系強風抖振空間響應的時/頻域的理論計算方法。具體結論如下:

(1)平均風荷載導致了幾何非線性的靜風響應,輸電線平面發生了風偏,水平風荷載還引起了平面內和平面外模態振動響應。與此同時,輸電線弦向張力也有了顯著提高,改變了系統的自振頻率和氣動阻尼等動力特性。

(2)利用影響線函數可方便高效地計算脈動風荷載的背景分量,且可以得出模態參與系數,用以計算共振分量。除了橫向支反力響應外,其他響應的影響線函數均會受到靜力風偏狀態的影響。由于較高的氣動阻尼,脈動響應以背景分量為主要貢獻,共振分量以一階振動為主。

(3)不考慮輸電線靜力風偏將顯著高估絕緣子橫向位移的背景分量及脈動響應,同時也無法計算豎向位移、縱向反力和豎向動反力分量。因此在動力響應分析中考慮靜力非線性變形對脈動風動力響應的影響是十分必要的。

(4)通過與非線性有限元分析比較,本文提出的三維振動理論模型具備足夠的工程效率和精度。

參考文獻:

[1] 鄧洪洲,朱松曄,陳曉明,等.大跨越輸電塔線體系氣彈模型風洞試驗[ J ].同濟大學學報(自然科學版),2003,31(2):132-137.

Deng? Hongzhou,Zhu? Songye,Chen? Xiaoming,et? al . Wind tunnel investigation on model of long span trans? mission? line? system [ J ]. Journal? of? Tongji? University (Natural Science),2003,31(2):132-137.

[2] 梁樞果,鄒良浩,韓銀全,等.輸電塔-線體系完全氣彈模型風洞試驗研究[ J ].土木工程學報,2010,43(5):70-78.

Liang? Shuguo, Zou? Lianghao, Han? Yinquan, et? al . Study of wind tunnel tests of a full aero-elastic model of electrical? transmission? tower-line? systems [ J ]. China Civil Engineering Journal,2010,43(5):70-78.

[3] 謝強,嚴承涌.1000 kV 特高壓交流同塔雙回輸電塔線耦聯體系風洞試驗[ J ].高電壓技術,2010,36(4):900-906.

XieQiang,Yan Chengyong . Wind tunnel test on 1000 kV UHV AC double circuit transmission tower conduc ? tor? coupling? system [ J ]. High? Voltage? Engineering,2010,36(4):900-906.

[4] 樓文娟,孫珍茂,許福友,等.輸電導線擾流防舞器氣動力特性風洞試驗研究[ J ].浙江大學學報(工學版),2011,45(1):93-98.

Lou Wenjuan,Sun Zhenmao,Xu Fuyou,et al . Experi? mental study? on? aerodynamic characteristics of air flow spoiler[ J ]. Journal of Zhejiang University (Engineering Science),2011,45(1):93-98.

[5] 李正良,任坤,肖正直,等.特高壓輸電塔線體系氣彈模型設計與風洞試驗[ J ].空氣動力學學報,2011,29(1):102-113.

Li Zhengliang,Ren? Kun,Xiao Zhengzhi,et? al . Aero ? elastic? model? design? and? wind? tunnel? tests? of? UHVtransmission? line? system[ J ]. Acta? Aerodynamica? Sini? ca,2011,29(1):102-113.

[6]? Wang? Dahai, Chen? Xinzhong, Li? Jie . Prediction? ofwind-induced buffeting response of overhead conductor:comparison? of linear? and? nonlinear? analysis? approaches [ J ]. Journal of Wind Engineering and Industrial Aerody? namics,2017,167:23-40.

[7]? Momomura Y,Marukawa H . Full-scale measurementsof wind-induced vibration of a transmission line system in a mountainous area[ J ]. Journal of Wind Engineering and Industrial Aerodynamics,1997,72:241-252.

[8]? Task Committee on Structural Loadings of the Commit?tee? on? Electrical Transmission? Structures? of the? Struc ? tural? Engineering? Institute? of the? American? Society? of Civil Engineers . Guidelines? for Electrical Transmission Line Structural Loading[M]. USA:American Society of Civil Engineers,2010.

[9]? Pasca? M ,Vestroni? F ,Gattulli? V . Active? longitudinalcontrol of wind-induced oscillations of a suspended cable [ J ]. Meccanica,1998,33(3):255-266.

[10] Davenport A? G . Gust response factors for transmissionline loading[ J ]. Wind Engineering,1980,2:899-909.

[11] Matheson,M J,Holmes J D . Simulation of the dynamic response of conductor in strong winds[ J ]. Engineer? ing Structures,1981,3(2):105-110.

[12] Loredo-Souza A M,Davenport A G . Wind tunnel aero ?elastic? studies? on? the? behaviour? of? two? parallel? cables [ J ]. Journal of Wind Engineering & Industrial Aerody? namics,2002,90(4):407-414.

[13] Di Paola M,Muscolino G,Sofi A . Monte Carlo simula ?tion for the response analysis of long-span suspended ca? bles under wind loads[ J ]. Wind & Structures,2004,7(2):107-130.

[14] Davenport A? G . Note on the distribution of the largestvalue? of? a? random? function? with? applications? to? gust loading[ J ]. Proceedings of the Institution of Civil Engi? neers,1964,28(2):187-196.

[15] Chen? X , Kareem? A . Evaluation? of? equivalent? staticwind loads on buildings[C]. Proceedings of 10th Ameri? cas? Conference? on? Wind? Engineering ,Baton? Rouge, LA,USA,2005.

[16] Peng L,Huang G,Kareem A,et al . An efficient space-time? based? simulation? approach? of? wind? velocity? field with? embedded? conditional? interpolation? for? unevenly spaced locations[ J ]. Probabilistic Engineering Mechan? ics,2016,43:156-168.

Frequency and time domain analytical methods for wind -induced buffeting response of overhead conductor

WANG Da-hai1,WANG Tao1,WANG Wei2,XU Kang3,LIANG Shu-guo4

(1.Department of Architectural Engineering,School of Civil Engineering and Architecture,Wuhan University of Technology,Wuhan 430070,China;2.Guangdong Electric Power Design Institute,Guangzhou 510799,China;3.State Grid Electric Power Research Institute Wuhan Nari Group Corporation,Wuhan 430072,China;4.School of Civil Engineering,Wuhan University,Wuhan 430070,China)

Abstract: The nonlinear dynamic buffeting response with large deformation often occurs under the action of strong wind,which will always lead to wind flashover of power network and the wind damage of supporting towers . In order to reveal the formation mechanism of wind effect on transmission lines under the action of strong wind,this paper takes the typical two-span transmission lines-insulators as the research object,and studies on the wind deflection displacement of insulators and the reaction response of towers to the spatial supports at the end of insulators . Based on the cable structure mechanics,the nonlinear solution of the static average wind deviation state is analyzed,and the changes of dynamic characteristics such as mode and aerodynamic damping of the system under wind deviation are investigated . The expressions of the influence line function and modal participation coefficient of the response are derived,and the calculation method of linear dynamic response in time domain is proposed . According to the theo? ry of wind engineering,the frequency domain expressions of background component and resonance component of fluctuating wind vibration response are given . A typical example is used to analyze and compare with the results of nonlinear finite element model . The results show that the three-dimensional vibration time/frequency domain theoretical model and the parameter calculation meth ? od proposed in this paper have sufficient engineering efficiency and accuracy . The study provides a theoretical basis and calculation method for revealing the damage mechanism of transmission line wind disaster and improving the wind-resistant design of transmis? sion line structure .

Key words : nonlinear;buffeting response;transmission line;background component;resonant component

作者簡介:汪大海(1975—),男,博士,教授。電話:13657236629;E-mail:wangdahai@whut .edu .cn。

主站蜘蛛池模板: 久久美女精品| 亚洲狠狠婷婷综合久久久久| 草草线在成年免费视频2| 国产中文一区a级毛片视频 | 国产乱人伦AV在线A| 国产激情国语对白普通话| 69精品在线观看| 久久国产精品电影| 国产福利影院在线观看| 99久久性生片| 美女视频黄频a免费高清不卡| 久久人妻xunleige无码| 激情综合婷婷丁香五月尤物| 无码人妻免费| 麻豆AV网站免费进入| 有专无码视频| 国产亚洲第一页| 成人午夜网址| 国产日韩久久久久无码精品| 亚洲国产精品一区二区第一页免| 老司机精品一区在线视频| 久久永久视频| 色偷偷av男人的天堂不卡| 青青青草国产| 色偷偷男人的天堂亚洲av| 先锋资源久久| 欧美精品二区| 亚洲aⅴ天堂| 免费精品一区二区h| 久久久久免费精品国产| 亚洲天堂网视频| 97se亚洲综合在线| 日韩毛片免费视频| 日韩少妇激情一区二区| 欧美色亚洲| 婷婷在线网站| 她的性爱视频| 性69交片免费看| 91精品国产综合久久香蕉922| 992tv国产人成在线观看| 国产精品区网红主播在线观看| 亚洲成肉网| 秘书高跟黑色丝袜国产91在线| 国产主播一区二区三区| 黄色国产在线| 国外欧美一区另类中文字幕| 国产高清在线观看| 亚洲香蕉伊综合在人在线| 潮喷在线无码白浆| 毛片视频网址| 国产网站在线看| 国产无遮挡猛进猛出免费软件| 日韩精品毛片| 国产精品久久久精品三级| 国产精品自拍露脸视频| 日韩第一页在线| 亚洲人成日本在线观看| 久久亚洲国产最新网站| 国产AV无码专区亚洲精品网站| 国产不卡在线看| 五月婷婷欧美| 亚洲中文字幕无码爆乳| 成人亚洲国产| 欧美一级在线看| 亚洲国产日韩欧美在线| 蜜桃视频一区| 干中文字幕| 国产精品视频猛进猛出| 久久夜色精品| 成人伊人色一区二区三区| 久久中文字幕2021精品| 久久福利片| 国产麻豆91网在线看| 国产成人麻豆精品| 在线毛片网站| 在线观看91香蕉国产免费| 精品人妻无码中字系列| 国产精品女在线观看| 日韩专区欧美| 美女啪啪无遮挡| 亚欧美国产综合| 91亚洲精品第一|