亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        分段式結(jié)構(gòu)對風(fēng)力機(jī)塔架振動特性的影響分析

        2017-06-28 16:24:20程友良薛占璞渠江曼
        噪聲與振動控制 2017年3期
        關(guān)鍵詞:模態(tài)振動結(jié)構(gòu)

        程友良,薛占璞,渠江曼,蔣 衍

        (華北電力大學(xué) 能源動力與機(jī)械工程學(xué)院,河北 保定 071003)

        分段式結(jié)構(gòu)對風(fēng)力機(jī)塔架振動特性的影響分析

        程友良,薛占璞,渠江曼,蔣 衍

        (華北電力大學(xué) 能源動力與機(jī)械工程學(xué)院,河北 保定 071003)

        針對分段式結(jié)構(gòu)引起的風(fēng)力發(fā)電塔架失穩(wěn)問題,利用剛體有限元法建立分段式結(jié)構(gòu)三維動力學(xué)模型,提取模態(tài)分析的前五階非零模態(tài)的振型,研究在湍流工況下分段式塔架連接處內(nèi)外法蘭的振型特性,通過內(nèi)外法蘭的振型與頻譜響應(yīng)對比,分析分段結(jié)構(gòu)對塔架振動特性的影響。結(jié)果表明:分段式的外法蘭比內(nèi)法蘭各階振動位移量小,振動變化幅度小,不易發(fā)生共振,外法蘭的第1階振型起主導(dǎo)作用,非線性的扭轉(zhuǎn)變形小,該結(jié)果可為塔架在運(yùn)行過程中避免發(fā)生共振及失穩(wěn)現(xiàn)象提供一定的參考。

        振動與波;塔架;分段式結(jié)構(gòu);頻譜;扭轉(zhuǎn)變形

        風(fēng)力機(jī)塔架為分段式結(jié)構(gòu)。塔架是整個風(fēng)力發(fā)電機(jī)組的支撐部件,分段的塔架連接處分為內(nèi)法蘭和外法蘭結(jié)構(gòu),法蘭主要承受風(fēng)的橫向剪切及塔架的扭轉(zhuǎn)等載荷作用,其振動特性直接關(guān)系到塔架結(jié)構(gòu)可靠性和運(yùn)行穩(wěn)定性。張豐豪等利用剛體有限元方法對塔架進(jìn)行三維結(jié)構(gòu)動力學(xué)建模,通過浮動坐標(biāo)系和坐標(biāo)變換的基本原理,得出塔架1階彎曲振動隨著結(jié)構(gòu)阻尼的增大主導(dǎo)作用減弱[1]。許斌等對混凝土填充鋼箱的連接段進(jìn)行彈性和彈塑性分析,發(fā)現(xiàn)連接段受力性能優(yōu)于傳統(tǒng)的法蘭盤連接,傳遞載荷比較均勻[2]。王振宇等利用葉素理論,應(yīng)用該風(fēng)向突然偏轉(zhuǎn)的工況,對塔架的應(yīng)力分布進(jìn)行分析,發(fā)現(xiàn)風(fēng)向突然偏轉(zhuǎn)90°時,塔架16.8 m高處發(fā)生屈服破壞[3]??率捞玫壤弥C波疊加法和改進(jìn)的葉素動量理論并結(jié)合風(fēng)振精細(xì)化頻域計算方法有效模擬了葉片和塔架的脈動風(fēng)速時程,發(fā)現(xiàn)離心力效應(yīng)減小了塔架的風(fēng)振系數(shù)[4]。應(yīng)有等出了一種塔架主動阻尼控制方案,并通過技術(shù)仿真和現(xiàn)場試驗(yàn)對其控制效果進(jìn)行了驗(yàn)證分析[5]。趙榮珍等運(yùn)用次空間迭代法對塔架進(jìn)行數(shù)值模擬,發(fā)現(xiàn)擺振是塔架的主要振動形式,1階固有頻率高于葉片通過頻率[6]。程友良等對分段式塔架進(jìn)行了數(shù)值模擬,發(fā)現(xiàn)分段式比錐筒式的應(yīng)力變化幅度小[7]。黃珊秋等通過模態(tài)分析使塔架的尺寸和材料得到了優(yōu)化優(yōu)化設(shè)計[8]。張羽等從工程實(shí)用角度出發(fā),總結(jié)了塔架的結(jié)構(gòu)選型、屈曲特性及靜動態(tài)特性等,并對塔架的未來研究趨勢作出了展望[8]。呂鋼根據(jù)多自由度模態(tài)分析理論,提出了彈性鉸力學(xué)模型,發(fā)現(xiàn)振動的能量主要集中在第1、2階頻率處[9]。

        在以往的研究中,通常只關(guān)注風(fēng)力機(jī)塔架整體的風(fēng)振效應(yīng),而塔架分段連接處的法蘭振動情況研究甚少,法蘭盤是各段連接的部件,對其進(jìn)行振動特性研究直接關(guān)系到整個塔架的穩(wěn)定性[10–13],受實(shí)驗(yàn)條件所限,數(shù)值模擬成為研究分段式內(nèi)外法蘭的有效手段。研究一定工況下內(nèi)外法蘭的振動特性,結(jié)合模態(tài)分析,分析了分段式內(nèi)外法蘭對塔架的影響,該結(jié)果在一定程度上為避免塔架的共振及失穩(wěn)現(xiàn)象提供參考。

        1 數(shù)學(xué)模型

        1.1 塔架分段式內(nèi)外法蘭數(shù)學(xué)模型

        以某額定功率2 MW風(fēng)力發(fā)電機(jī)塔架為研究對象,塔架采用內(nèi)、外法蘭兩種連接方式,分別建立兩個模型,內(nèi)法蘭外徑4 837 mm,內(nèi)徑4 737 mm,厚度90 mm,外法蘭與內(nèi)法蘭的幾何尺寸相同,如圖1為塔架內(nèi)法蘭,圖2為塔架外法蘭。材料為16Mn,密度為7 900 kg﹒m-3,彈性模量206 GPa,泊松比0.26。塔架內(nèi)、外法蘭三維結(jié)構(gòu)動力學(xué)模型建立時,一般采用模態(tài)分析法或者有限元法,但這兩種方法各有其不足之處,模態(tài)分析法建模時需要事先計算各階固有頻率及振型,一般考慮前2階彎曲模態(tài)振型,而忽略了扭振及其他方向的振型情況。

        對于小變形假設(shè)的有限元法不考慮各階振型位移之間耦合作用。鑒于此,本文利用剛體有限元法對塔架分段式內(nèi)外法蘭進(jìn)行建模,該方法采用浮動坐標(biāo)和坐標(biāo)變換,可以進(jìn)行柔性部分大變形及各方向耦合的非線性計算。該方法的基本原理是是將柔性結(jié)構(gòu)分割為由彈性阻尼節(jié)點(diǎn)連接的剛體單元,分別計算剛體單元的動能、重力勢能以及彈性阻尼節(jié)點(diǎn)的彈性勢能和阻尼耗散能,利用拉格朗日方程進(jìn)行動力學(xué)建模[1]。其優(yōu)點(diǎn)在于可計算任意變形范圍內(nèi)的振動,統(tǒng)一采用拉格朗日方程建模,還可以節(jié)省計算時間,在相同條件下,只需要有限元法的三分之一。圖1為剛體有限元法示意圖,圖2為風(fēng)力機(jī)塔架動力學(xué)模型,圖3、圖4分別為塔架外、內(nèi)法蘭結(jié)構(gòu)動力學(xué)模型。

        圖1 剛體有限元法示意圖

        圖2 風(fēng)力機(jī)塔架動力學(xué)模型

        圖3 塔架外法蘭

        圖4 塔架內(nèi)法蘭

        1.2 基本控制方程

        塔架內(nèi)外法蘭滿足的基本控制方程,圖2、圖3中模型的3個平動位移(x1,x2,x3)和3個轉(zhuǎn)動位移(φ1,φ2,φ3),將每個剛體單元的自由度統(tǒng)一用qi表示

        剛體單元中任意一點(diǎn)在參考坐標(biāo)系中的位移可寫為

        式中Ti——從坐標(biāo)系{}′向坐標(biāo)系{}的變換矩陣,

        r′——固結(jié)坐標(biāo)系{}′中的位置用矢量

        r′表示,r——在坐標(biāo)系{}中的位置用r表示,

        剛體單元的動能為

        式中tr——矩陣的秩;B——坐標(biāo)變換矩陣;變量上方加“.”——變量對時間的導(dǎo)數(shù)。

        模態(tài)是機(jī)械結(jié)構(gòu)的固有振動特性,每一個模態(tài)具有特定的固有頻率、阻尼比和模態(tài)振型[14–17]。具有n個自由度的無阻尼振動系統(tǒng)的振動微分方程可表示為

        式中M——塔架結(jié)構(gòu)的整體質(zhì)量矩陣,

        K——為塔架結(jié)構(gòu)的整體剛度矩陣。

        非零矢量?i的比例解為

        將n個特征矢量?i按列排成n×n階陣,可得系統(tǒng)的特征矢量矩陣

        此時特征矢量為模態(tài)矢量或模態(tài)矩陣,即為模態(tài)振型。

        2 計算工況及方法

        2.1 邊界條件及方法

        為了對比塔架內(nèi)外法蘭振動特性,在塔高35 m法蘭連接處施加50年一遇良態(tài)風(fēng),該脈動風(fēng)風(fēng)速波動值基本在±8 m/s范圍內(nèi)如圖5。

        塔架分為兩段,法蘭位于35 m處,塔高70 m,計算時間200 s,時間步長為10-5s,將塔架內(nèi)外法蘭的模態(tài)分析結(jié)果進(jìn)行對比,分析其中的影響因素。

        2 計算工況及方法

        2.1 邊界條件及方法

        為了對比塔架內(nèi)外法蘭振動特性,在塔高35 m法蘭連接處施加50年一遇良態(tài)風(fēng),該脈動風(fēng)風(fēng)速波動值基本在±8 m/s范圍內(nèi)如圖5。

        圖5 塔高35 m處脈動風(fēng)速(50年一遇)

        塔架分為兩段,法蘭位于35 m處,塔高70 m,計算時間200 s,時間步長為10-5s,將塔架內(nèi)外法蘭的模態(tài)分析結(jié)果進(jìn)行對比,分析其中的影響因素。

        2.2 數(shù)值計算的準(zhǔn)確性驗(yàn)證

        為了保證數(shù)學(xué)模型的正確性與可靠性,進(jìn)行對比驗(yàn)證。阻尼比為0.025時,將本文模型分段式的內(nèi)、外法蘭計算結(jié)果與獲得認(rèn)證的風(fēng)力機(jī)結(jié)構(gòu)計算程序FAST結(jié)果進(jìn)行對比驗(yàn)證[18–20]。圖6分別給出了內(nèi)、外法蘭橫向位移對比結(jié)果,可得出本文模型與FAST計算程序的結(jié)果相吻合。塔架分段處的橫向位移量約為0.1 m,與現(xiàn)實(shí)相符。

        圖6 塔架內(nèi)外法蘭橫向位移對比

        對分段式內(nèi)法蘭振動位移曲線進(jìn)行快速傅里葉變換(FFT),利用獲得的振動頻率與有限元程序Bmode的計算結(jié)果進(jìn)行對比,提取塔架的橫向前五階振動頻率,如圖7可知,本文內(nèi)法蘭分段式的振動頻率小于計算機(jī)有限元程序Bmode的計算結(jié)果,誤差較小,各階振動頻率相對誤差在4%以內(nèi),充分證明了本文模型進(jìn)行數(shù)值研究的可行性。

        圖7 振動頻率對比

        3 計算結(jié)果與討論

        3.1 分段式內(nèi)法蘭模態(tài)振型圖分析

        在塔架結(jié)構(gòu)動力響應(yīng)中,前5階非零模態(tài)振型是所有振型的基礎(chǔ),與其對應(yīng)的模態(tài)頻率是振動特性分析的重點(diǎn),確定塔架結(jié)構(gòu)的模態(tài)頻率直接決定塔架在什么范圍內(nèi)振動幅度較大。由圖8可知,內(nèi)法蘭1階模態(tài)振型(f=0.76 Hz)至5階模態(tài)振型(f= 11.16 Hz),最大位移量為7.2 mm,前2階振型中,看不出明顯的扭轉(zhuǎn)變形,第3、第4階模態(tài)振型顯示徑向擠壓振動并伴隨扭轉(zhuǎn)變形,第5階顯示橫向扭轉(zhuǎn)變形。由于在振動過程中,塔架的能量體現(xiàn)在前兩階振型對應(yīng)的頻率中,所以擺振是該內(nèi)法蘭的主要振型。當(dāng)外加載荷頻率和固有頻率相同時,即發(fā)生共振現(xiàn)象。圖7表明各階振型位移量將提供共振參考,如,塔架內(nèi)法蘭振動測量的位移達(dá)到1.008 mm,對應(yīng)的頻率為0.76 Hz,此時將發(fā)生共振,應(yīng)避免在該頻率下運(yùn)行。

        3.2 分段式外法蘭模態(tài)振型圖分析

        圖9給出了分段式外法蘭前5階振型圖,1階模態(tài)振型(f=0.57 Hz)至5階模態(tài)振型(f=12.24 Hz),最大位移量為6.06 mm,第3、第4階模態(tài)振型顯示徑向擠壓振動并伴隨扭轉(zhuǎn)變形,第5階顯示橫向扭轉(zhuǎn)變形,擺振依然是該外法蘭的主要振型。對于塔架內(nèi)法蘭而言,外法蘭在前5階的頻率分布及振動位移量均有優(yōu)勢。內(nèi)法蘭前5階對應(yīng)頻率平均間隔為1.301 Hz,而外法蘭為2.31 Hz,且各階振動位移量比內(nèi)法蘭小。因此,外法蘭塔架相對于內(nèi)法蘭塔架更不易發(fā)生共振,且振動變化小。

        3.3 分段式結(jié)構(gòu)對塔架振動特性的影響

        根據(jù)圖5邊界條件,結(jié)合內(nèi)外法蘭結(jié)構(gòu)動力響應(yīng),圖10、圖11給出了在不同頻率下的頻譜局部放大圖,可以得出外法蘭彎曲振動對應(yīng)的高階頻率比內(nèi)法蘭少,并且變化幅度小,說明塔架在實(shí)際運(yùn)行中,前3階非零模態(tài)振型及頻率對塔架的正常運(yùn)行影響極大,甚至無需考慮第4、第5階彎曲振動造成的影響。

        圖8 內(nèi)法蘭前5階振型圖

        圖9 外法蘭前5階振型圖

        圖10 內(nèi)法蘭頻譜局部放大

        圖11 外法蘭頻譜局部放大

        從圖10發(fā)現(xiàn)塔架隨著頻率的逐漸增大,內(nèi)法蘭塔架第1階振動的主導(dǎo)作用逐漸減弱,振動特性主要源于外界湍流風(fēng)的變化。在外湍流風(fēng)變化極大的情況下,內(nèi)法蘭塔架的振動位移比對應(yīng)的外法蘭塔架位移量大,并且非線性作用占主導(dǎo)地位。內(nèi)法蘭極有可能在1階振動頻率發(fā)生共振,引起風(fēng)力發(fā)電系統(tǒng)失穩(wěn)。而外法蘭相對于內(nèi)法蘭在1階振動頻率發(fā)生共振的幾率小,因此應(yīng)該避免該頻率范圍內(nèi)引起的風(fēng)力發(fā)電系統(tǒng)失穩(wěn)現(xiàn)象發(fā)生。

        4 結(jié)語

        利用剛體有限元法建立塔架分段處內(nèi)外法蘭三維數(shù)學(xué)模型,進(jìn)行結(jié)構(gòu)動力學(xué)分析,研究了內(nèi)外法蘭在一定工況下的振動特性,并重點(diǎn)分析了內(nèi)外法蘭的振型變化情況,結(jié)論如下:

        (1)外法蘭的振動幅值相對于內(nèi)法蘭變化較小,振型中的位移量較小,擺振是內(nèi)外法蘭的主要振動形式,并且伴隨輕微的扭轉(zhuǎn)振動和彎曲振動,以前的研究注重塔架的彎曲振動,但隨著柔性系統(tǒng)的發(fā)展,應(yīng)該注重扭轉(zhuǎn)振動的研究。

        (2)內(nèi)外法蘭的1階振型在塔架振動過程中起重要作用,其中非線性位移量是導(dǎo)致失穩(wěn)的重要因素,外法蘭相對于內(nèi)法蘭發(fā)生共振的幾率小。

        (3)內(nèi)法蘭彎曲振動的高階頻率比外法蘭多,應(yīng)增加以柔性系統(tǒng)為結(jié)構(gòu)的減震裝置,避免發(fā)生失穩(wěn)現(xiàn)象。內(nèi)法蘭是比較傳統(tǒng)的塔架連接方式,但是外法蘭的振動特性優(yōu)于內(nèi)法蘭,外法蘭的結(jié)構(gòu)優(yōu)化及性能分析是未來研究的重點(diǎn)。

        [1]張豐豪,何榕.結(jié)構(gòu)阻尼對風(fēng)力機(jī)塔架振動特性的影響[J].太陽能學(xué)報,2015,36(10):2467-2473.

        [2]許斌,李正超,謝詠劍.組合風(fēng)電塔架混凝土填充鋼箱連接段數(shù)值模擬[J].機(jī)械設(shè)計與制造,2016,(3):30-33.

        [3]王振宇,張彪,趙艷,等.臺風(fēng)作用下風(fēng)力機(jī)塔架振動響應(yīng)研究[J].太陽能學(xué)報,2013,34(8):1434-1442.

        [4]柯世堂,王同光,曹九發(fā),等.考慮葉片旋轉(zhuǎn)和離心力效應(yīng)風(fēng)力機(jī)塔架風(fēng)振分析[J].太陽能學(xué)報,2015,36(1):33-40.

        [5]應(yīng)有,朱重喜,楊帆,等.大型風(fēng)電機(jī)組塔架主動阻尼控制技術(shù)研究[J].太陽能學(xué)報,2015,36(1):54-60.

        [6]趙珍榮,呂鋼.水平軸風(fēng)力發(fā)電機(jī)塔架的振動模態(tài)分析[J].蘭州理工大學(xué)學(xué)報,2009,35(2):33-36.

        [7]程友良,薛占璞,渠江曼.新型分段式風(fēng)力發(fā)電塔結(jié)構(gòu)改進(jìn)及性能研究[J].制造業(yè)自動化,2016,38(9):45-49.

        [8]張羽,蔡新,高強(qiáng),等.風(fēng)力機(jī)塔架結(jié)構(gòu)研究概述[J].工程設(shè)計學(xué)報,2016,23(2):108-123.

        [9]呂鋼.基于有限元法的水平軸風(fēng)力機(jī)塔架動態(tài)響應(yīng)與優(yōu)化問題研究[D].蘭州:蘭州理工大學(xué),2009.

        [10]HUANG LINHONG,SONG LILI,LI GANG,et al. Variation Characteristics of Regional Synchronous Wind in Hami[J].Journal of Meteorological Research,2014, 29:344-357.

        [11]WITTBRODT E,ADAMIEC-WóJCIK I.Dynamics of flexible multibody systems:Rigid finite element method [M].Berlin:Springer,2006:5-82.

        [12]黃珊秋,陸萍.ZONDZ-40風(fēng)力機(jī)塔架的模態(tài)分析[J].太陽能學(xué)報,2001,22(2):153-156.

        [13]龍凱,賈嬌.大型水平軸風(fēng)力機(jī)塔筒渦激振動焊縫疲勞分析[J].太陽能學(xué)報,2015,36(10):2455-2459.

        [14]程友良,薛占璞,楊國寧.一種風(fēng)力發(fā)電機(jī)塔架[P].中國專利:ZL201420777451.1,2015-06-17.

        [15]戴建鑫.風(fēng)力機(jī)塔架的有限元建模及靜動態(tài)特性分析[D].蘭州:蘭州理工大學(xué),2011.

        [16]趙世林,李德源,黃小華.風(fēng)力機(jī)塔架在偏心載荷作用下的屈曲分析[J].太陽能學(xué)報,2010,31(7):901-906.

        [17]王印軍,任勇生,孫丙磊,等.基于Ansys的垂直軸風(fēng)力機(jī)塔架的力學(xué)分析及結(jié)構(gòu)優(yōu)化[J].山東科技大學(xué)學(xué)報,2011,30(5):96-102.

        [18]陸萍,秦惠芳,欒芝云.基于有限元法的風(fēng)力機(jī)塔架結(jié)構(gòu)動態(tài)分析[J].機(jī)械工程學(xué)報,2002,38(9):127-130.

        [19]劉雄,李鋼強(qiáng),陳嚴(yán),等.水平軸風(fēng)力機(jī)筒型塔架動態(tài)響應(yīng)分析[J].太陽能學(xué)報,2010,31(4):412-417.

        [20]李源霜,李汶柏,張文明.大型橢圓振動篩空心裂紋轉(zhuǎn)軸非線性振動特性分析[J].噪聲與振動控制,2016,36(2):46-51.

        Analysis of the Influence of Sectional Structure on Vibration Characteristics of Wind Turbine Towers

        CHENG You-liang,XUE Zhan-pu,QU Jiang-man,JIANG Yan
        (School of Energy Power and Mechanical Engineering,North China Electric Power University, Baoding 071003,Hebei China)

        The instability problem of the wind power towers induced by the sectional structure is studied.The rigid finite element method is used to establish the three-dimensional dynamic model of the sectional structure.The first five nonzero modals are extracted.The modal characteristics of the inner and outer flanges of the joints of the tower are studied under the turbulent conditions.The results show that the vibration displacements of the outer flange are all smaller than those of the inner flange.The magnitudes of vibration variation are small,so the resonance is unlikely to occur.The first order vibration of the outer flange plays a dominant role in nonlinear torsional deformation.The results provide a reference for the towers to avoid the resonance in the operation process.

        vibration and wave;tower;sectional structure;vibration spectrum;torsional deformation

        TK83

        :A

        :10.3969/j.issn.1006-1355.2017.03.007

        1006-1355(2017)03-0037-05

        2016-11-11

        國家自然科學(xué)基金重點(diǎn)基金(11232012);中央高?;究蒲袠I(yè)務(wù)費(fèi)專項資金項目(2016XS107)

        程友良(1963-),男,博士,教授,博士生導(dǎo)師,主要研究方向?yàn)榱黧w動力學(xué)理論及其應(yīng)用、流體設(shè)備與節(jié)能、可再生能源理論及其應(yīng)用、清潔能源利用技術(shù)與設(shè)備。

        薛占璞(1983-),博士研究生,專業(yè)為流體機(jī)械及工程,研究方向?yàn)榍鍧嵞茉蠢眉夹g(shù)與設(shè)備。E-mail:shenghuo166@163.com

        猜你喜歡
        模態(tài)振動結(jié)構(gòu)
        振動的思考
        《形而上學(xué)》△卷的結(jié)構(gòu)和位置
        振動與頻率
        論結(jié)構(gòu)
        中華詩詞(2019年7期)2019-11-25 01:43:04
        中立型Emden-Fowler微分方程的振動性
        論《日出》的結(jié)構(gòu)
        國內(nèi)多模態(tài)教學(xué)研究回顧與展望
        創(chuàng)新治理結(jié)構(gòu)促進(jìn)中小企業(yè)持續(xù)成長
        基于HHT和Prony算法的電力系統(tǒng)低頻振蕩模態(tài)識別
        UF6振動激發(fā)態(tài)分子的振動-振動馳豫
        計算物理(2014年2期)2014-03-11 17:01:44
        蜜桃臀av一区二区三区| 久久av无码精品人妻糸列| 日韩欧美精品有码在线观看| 成人全部免费的a毛片在线看| 在线观看av网站永久| 少妇被猛男粗大的猛进出| 亚洲国产美女在线观看| 一区二区三区在线观看视频免费 | 国产97在线 | 亚洲| 黄色毛片视频免费| 亚洲国产免费一区二区| 手机看片自拍偷拍福利| 亚洲国产免费不卡视频| 中文字幕日韩三级片| 久久人人97超碰超国产| 国产香蕉尹人在线视频你懂的| 日本在线观看一二三区| 久久精品国产久精国产果冻传媒| 狠狠爱无码一区二区三区| 午夜天堂精品一区二区| 激情文学婷婷六月开心久久 | 国产乱子乱人伦电影在线观看| 东北无码熟妇人妻AV在线| 一本色道久久综合中文字幕| 亚洲av高清一区二区在线观看| 成人午夜视频精品一区| 国产手机在线αⅴ片无码| 国产精品女人一区二区三区 | 日韩爱爱网站| 蜜臀一区二区av天堂| 国产日产欧产精品精品蜜芽| 中文人妻av久久人妻18| 亚洲片在线视频| 国产精品久久婷婷六月| 亚洲国产a∨无码中文777| 好吊色欧美一区二区三区四区 | 人妻少妇偷人精品免费看| 亚洲色www成人永久网址| 国产综合激情在线亚洲第一页| 宅男久久精品国产亚洲av麻豆| 亚洲一区二区国产激情|