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

風(fēng)力發(fā)電機(jī)組塔筒渦致橫振研究

2012-07-02 10:48:40董占琢
東方汽輪機(jī) 2012年2期
關(guān)鍵詞:模態(tài)模型

董占琢 廖 暉

(東方汽輪機(jī)有限公司,四川 德陽,618000)

0 引言

高聳結(jié)構(gòu)橫風(fēng)向振動(dòng)的機(jī)理較為復(fù)雜,影響因素很多,在工程結(jié)構(gòu)中較為常見且機(jī)理相對(duì)清楚的橫向風(fēng)振內(nèi)容包括:渦激振動(dòng)[1]、馳振[2]、顫振[3]。本文所研究的塔筒橫截面為規(guī)則圓形,不存在攻角問題,風(fēng)繞流塔筒時(shí)不會(huì)發(fā)生馳振和顫振,主要是由卡門渦街的漩渦發(fā)放引起的垂直于來流方向的渦激振動(dòng)。風(fēng)繞流塔筒產(chǎn)生的卡門渦街以及升阻力方向如圖1所示。

圖1 卡門渦街與升阻力方向

1 塔筒渦激振動(dòng)CFD計(jì)算

1.1 雷諾相似準(zhǔn)則

風(fēng)繞塔筒的流動(dòng)主要受粘性力、壓力和慣性力的作用。從力學(xué)相似的觀點(diǎn)看,若兩個(gè)流場在對(duì)應(yīng)點(diǎn)作用的同種力方向相同、大小成同一比例,則滿足動(dòng)力相似。在幾何相似的前提下,兩個(gè)流動(dòng)只要在對(duì)應(yīng)點(diǎn)滿足代表粘性力與慣性力比值的雷諾數(shù)相等,則表示壓力與慣性力之比的歐拉準(zhǔn)則必然相等,因此風(fēng)繞流塔筒滿足雷諾相似條件。只需給出不同雷諾數(shù)下的力系數(shù)的大小即可表示不同直徑、不同風(fēng)速下的受力。

1.2 二維圓形CFD計(jì)算

圓柱繞流問題是典型的鈍體大分離問題,本文采用SST湍流模型進(jìn)行計(jì)算。

1.2.1 計(jì)算網(wǎng)格

采用圖2所示的計(jì)算網(wǎng)格。計(jì)算域大小為10D×20D,計(jì)算網(wǎng)格總數(shù)34302。

圖2 計(jì)算網(wǎng)格與邊界條件

1.2.2 邊界條件

入口設(shè)置為速度入口條件,按照不同的雷諾數(shù)進(jìn)行計(jì)算;出口設(shè)置為背壓出口,出口為大氣壓;壁面設(shè)置為無滑移壁面條件;計(jì)算域上下邊界為對(duì)稱邊界,用以模擬無窮大空間。

1.2.3 計(jì)算結(jié)果

主要計(jì)算塔筒表面受力情況,圖3為隨時(shí)間變化圓柱表面的升力系數(shù)與阻力系數(shù)。

圖3 Re=3.34×106時(shí)的圓形升、阻力系數(shù)變化

由圖3可以看出,升力系數(shù)和阻力系數(shù)都隨時(shí)間周期按正弦或余弦變化,阻力系數(shù)變化幅值很小,可視為不變。升力系數(shù)呈關(guān)于0線的余弦變化,升力系數(shù)以其幅值Cl給出,頻率以斯特勞哈爾數(shù)的形式給出。

1.2.4 CFD二維計(jì)算模型與三維計(jì)算模型的比較

我們對(duì)圓柱的三維和二維建模都進(jìn)行了分析,結(jié)果對(duì)比如表1所示:

表1 二維、三維計(jì)算結(jié)果的比較

三維和二維的計(jì)算結(jié)果中,阻力系數(shù)和斯特勞哈爾數(shù)基本相同,而升力系數(shù)中三維的結(jié)果明顯比二維要小[3]。二維結(jié)果更加保守[4],以二維圓形計(jì)算結(jié)果來作為塔筒載荷。

1.2.5 CFD計(jì)算結(jié)果同DIN1055-4標(biāo)準(zhǔn)的比較

GL規(guī)范[5]規(guī)定橫向振動(dòng)載荷可按DIN1055-4[6]來計(jì)算。圖4給出了DIN1055-4標(biāo)準(zhǔn)查表數(shù)值和本文CFD計(jì)算數(shù)值的比較。

圖4 升、阻力系數(shù)DIN標(biāo)準(zhǔn)和CFD計(jì)算結(jié)果比較圖

本文的計(jì)算忽略了轉(zhuǎn)捩等復(fù)雜因素的影響,因而趨勢(shì)與實(shí)驗(yàn)趨勢(shì)有所偏差。將在下一節(jié)具體分析偏差產(chǎn)生的原因及修正方法。

1.2.6 二維CFD計(jì)算結(jié)果的修正

旋渦脫落引起的力是復(fù)雜的流體力學(xué)現(xiàn)象的結(jié)果,對(duì)描述流體和結(jié)構(gòu)物理特性的許多參數(shù)敏感。下面對(duì)轉(zhuǎn)捩、表面粗糙度、來流湍流度、實(shí)驗(yàn)影響等因素對(duì)結(jié)果的影響進(jìn)行分析。

采用關(guān)聯(lián)轉(zhuǎn)捩的模型之后,計(jì)算所得的圓柱的阻力系數(shù)與Schlichting[7]實(shí)驗(yàn)曲線的趨勢(shì)能夠吻合,如圖5所示。但由于采用的是經(jīng)驗(yàn)關(guān)聯(lián)式,與實(shí)驗(yàn)值的絕對(duì)值有一定的差別。但從圖6中可以看到,考慮轉(zhuǎn)捩之后的結(jié)果比全湍流的結(jié)果偏小。因此,盡管全湍流計(jì)算模型有一定局限,為了分析的安全性,我們?nèi)匀徊捎萌牧鞯挠?jì)算結(jié)果作為以后的分析基礎(chǔ)。

圖5 考慮轉(zhuǎn)捩后的阻力系數(shù)計(jì)算結(jié)果與實(shí)驗(yàn)值比較

圖6 轉(zhuǎn)捩對(duì)圓柱阻力系數(shù)的影響

粗糙表面對(duì)計(jì)算結(jié)果產(chǎn)生了較大的影響,從圖7可以看出,總體趨勢(shì)是使阻力系數(shù)和升力系數(shù)提高,而使斯特勞哈爾數(shù)降低。

圖7 表面粗糙度對(duì)圓柱升、阻力系數(shù)及斯特勞哈爾數(shù)的影響

湍流度對(duì)圓柱繞流有一定影響,如圖8所示,湍流度增大會(huì)使阻力系數(shù)降低。

圖8 湍流度對(duì)圓柱阻力系數(shù)的影響

目前,有關(guān)圓柱繞流所受阻力的實(shí)驗(yàn)值一般取Schlichting[7]的數(shù)值,然而,來流湍流度、圓柱表面粗糙度、壓力測(cè)量方式等都會(huì)對(duì)實(shí)驗(yàn)結(jié)果產(chǎn)生影響,實(shí)驗(yàn)結(jié)果也有很大誤差。

考慮表面粗糙度、來流湍流度等的變化和影響,升、阻力系數(shù)的修正以偏安全計(jì)算作為準(zhǔn)則,以防止實(shí)際中由于各種參數(shù)等的變化產(chǎn)生危險(xiǎn)工況。由于本文后續(xù)諧響應(yīng)計(jì)算方法,對(duì)于斯特勞哈爾數(shù)的精度要求不高,因此對(duì)其不作修正。對(duì)于升、阻力系數(shù),本文推薦修正系數(shù)為1.5。經(jīng)過修正后的升阻力系數(shù)如表2所示 (篇幅原因,此表僅為部分計(jì)算結(jié)果)。

表2 修正后的圓柱升力系數(shù)、阻力系數(shù)與斯特勞哈爾數(shù)

2 塔筒模態(tài)分析

為了分析FD82E風(fēng)力機(jī)塔筒在正常運(yùn)行與吊裝過程中的氣流激振安全性,需對(duì)其正常運(yùn)行與各吊裝環(huán)節(jié)進(jìn)行模態(tài)分析。分為四種工況,具體定義見表3。

表3 FD82E風(fēng)力機(jī)塔筒模態(tài)分析計(jì)算工況

下面以工況1為例,介紹模態(tài)分析流程及方法。

2.1 工況1塔筒有限元模型及邊界條件

塔筒按圖紙建模。塔筒頂部用一集中質(zhì)量點(diǎn)模擬機(jī)艙、輪轂和葉片的總質(zhì)量,坐標(biāo) (m)為(-0.5957, 0, 67.867), 總質(zhì)量為 106566 kg, Z方向轉(zhuǎn)動(dòng)慣量Iz=3540000 kg·m2。塔筒底部同樣用一集中質(zhì)量點(diǎn)模擬地基, 坐標(biāo)為 (0,0,-1.5896),總質(zhì)量為870040 kg,X、Y方向轉(zhuǎn)動(dòng)慣量Ix=Iy=9165000 kg·m2。集中質(zhì)量點(diǎn)分別與塔筒頂部和底部面建立剛性區(qū)域。地基質(zhì)量點(diǎn)周圍建立X、Y方向兩個(gè)扭轉(zhuǎn)彈簧模擬土壤約束,彈簧剛度系數(shù)取3×1010N·m/rad,彈簧末端節(jié)點(diǎn)坐標(biāo)為(5, 0, -1.5895), (0, 5, -1.5895)。整個(gè)塔筒網(wǎng)格均為六面體,壁厚方向劃分三層網(wǎng)格。模型總單元數(shù)為38866個(gè),總節(jié)點(diǎn)數(shù)為189358個(gè)。有限元模型中焊縫壁厚突跳處以及門框和筒體之間采用bonded接觸。地基質(zhì)量點(diǎn)加UX,UY,UZ,ROTZ四個(gè)方向約束,扭轉(zhuǎn)彈簧末端采用全約束,如圖9所示。

圖9 工況1條件下塔筒模態(tài)分析有限元模型

2.2 工況1模態(tài)計(jì)算結(jié)果及安全性校核

工況1模態(tài)計(jì)算結(jié)果見表4,具體振型如圖10所示。

表4 塔筒前6階模態(tài)計(jì)算結(jié)果

圖10 工況1塔筒X、Y向前兩階彎曲振動(dòng)振型

已知發(fā)電機(jī)轉(zhuǎn)速為1100~2000 r/min,風(fēng)輪轉(zhuǎn)速為 10.58~19.23 r/min,葉片 1P轉(zhuǎn)動(dòng)頻率為0.176~0.321 Hz, 3P 轉(zhuǎn)動(dòng)頻率為 0.528~0.962 Hz。滿足GL規(guī)范[5]要求的5%避開率。

2.3 所有工況模態(tài)分析結(jié)果

工況2、3、4的模態(tài)結(jié)果如表5所示,無風(fēng)輪激振,風(fēng)致振的校核會(huì)在后續(xù)章節(jié)分析。

表5 所有工況塔筒前6階模態(tài)計(jì)算結(jié)果

3 塔筒渦激振動(dòng)諧響應(yīng)分析

3.1 諧響應(yīng)分析載荷簡化

實(shí)際中由于塔筒外徑及風(fēng)速隨高度不斷變化,因此雷諾數(shù)沿高度方向不為定值,造成實(shí)際載荷狀況十分復(fù)雜,難以分析。本文對(duì)風(fēng)繞流塔筒載荷做如下簡化假設(shè):即假設(shè)塔筒沿高度方向作用的激振力頻率、幅值和相位均相同,激振力幅值取共振頻率附近CFD計(jì)算得到的最大值。此種簡化方法使激振力對(duì)塔筒的作用放大。

實(shí)際塔筒在風(fēng)繞流作用下,某一高度塔筒橫截面壓力分布不均,且隨時(shí)間變化。CFD中的載荷結(jié)果是將橫截面的壓力換算成作用在中心軸上的集中力。考慮到載荷簡化,沿高度方向塔筒載荷相當(dāng)于作用于中心軸的均布載荷,進(jìn)一步簡化為塔筒頂部的集中載荷和彎矩。

圖11 簡化載荷雷諾數(shù)計(jì)算取值說明

3.2 工況1條件下塔筒諧響應(yīng)分析

此工況為風(fēng)力機(jī)吊裝完畢的正常運(yùn)行條件,包括所有風(fēng)力機(jī)部件。在考慮風(fēng)繞流塔筒的橫向激振力之外還需考慮風(fēng)輪轉(zhuǎn)動(dòng)對(duì)塔筒的激振作用。

由于CFD計(jì)算本身存在一些誤差,安全起見取與工況1固有頻率0.453Hz最為接近的較大結(jié)果 (參見表2), 此時(shí)雷諾數(shù)約為 1.95×106, 風(fēng)速約為7m/s,均布載荷為58.73N/m (依照GL規(guī)范,此載荷已乘以安全系數(shù)1.35)。假設(shè)風(fēng)來流方向在X向,則風(fēng)繞流塔筒產(chǎn)生的橫向激振力在Y向,由理論力學(xué)簡化得到塔筒頂部集中載荷幅值為Fy=3884.99N,Mx=128496.04N·m。從計(jì)算結(jié)果可以看出,風(fēng)力機(jī)正常運(yùn)行時(shí)塔筒共振頻率附近的風(fēng)繞流激振力很小 (僅有不到4kN),對(duì)塔筒影響也很小,可以忽略。

3.3 工況2條件下塔筒諧響應(yīng)分析

3.3.1 工況2塔筒諧響應(yīng)分析有限元模型

工況2塔筒渦激振動(dòng)諧響應(yīng)分析有限元模型除塔頂質(zhì)量點(diǎn)的設(shè)置與工況1模態(tài)分析的有限元模型有差別外,其他部分完全相同。工況2諧響應(yīng)分析塔頂質(zhì)量點(diǎn)用于模擬機(jī)艙質(zhì)量,坐標(biāo) (m)為 (-1.19,-0.02996,67.69), 質(zhì)量為69000kg,Z方向轉(zhuǎn)動(dòng)慣量Iz=359500 kg·m2,X方向轉(zhuǎn)動(dòng)慣量Ix=91270 kg·m2, Y 方向轉(zhuǎn)動(dòng)慣量 Iy=393600 kg·m2。塔筒頂部同時(shí)建立一無質(zhì)量點(diǎn),用于加載集中載荷,坐標(biāo)為(0,0,66.15)。網(wǎng)格均為六面體,塔筒壁厚方向劃分三層網(wǎng)格。模型總單元數(shù)為38867個(gè),總節(jié)點(diǎn)數(shù)為189359個(gè),如圖12所示。

圖12 工況2條件下塔筒諧響應(yīng)分析有限元模型

3.3.2 工況2加載方式與邊界條件

取與工況2固有頻率0.542Hz最為接近的較大結(jié)果 (參見表2),此時(shí)雷諾數(shù)約為2.5×106,風(fēng)速約為9m/s,均布載荷為95.18N/m(已考慮安全系數(shù)1.35)。經(jīng)簡化得到塔筒頂部集中載荷幅值為Fy=6295.83N,Mx=208234.45N·m,阻尼系數(shù)值取0.005,相位角均取零,激振頻率范圍為0.541~0.5425Hz,計(jì)算步數(shù)為15步。有限元模型中焊縫壁厚突跳處以及門框和筒體之間采用bonded接觸。地基質(zhì)量點(diǎn)加UX,UY,UZ,ROTZ四個(gè)方向約束,扭轉(zhuǎn)彈簧末端采用全約束。

3.3.3 工況2諧響應(yīng)計(jì)算結(jié)果

由于結(jié)構(gòu)中定義了阻尼,所以結(jié)構(gòu)響應(yīng)與激振力之間不同步,存在相位差[5]。ANSYS計(jì)算的結(jié)果會(huì)以復(fù)數(shù)形式存儲(chǔ),而實(shí)際的計(jì)算結(jié)果是由實(shí)部和虛部的SRSS值給出。

由振動(dòng)理論可知,諧響應(yīng)激勵(lì)作用下塔筒一階共振時(shí)塔筒頂部位移最大,因此取塔筒頂部節(jié)點(diǎn)為觀察節(jié)點(diǎn)。可觀察到的塔筒頂部節(jié)點(diǎn)Y向振幅隨激振頻率變化,如圖13所示。

圖13 工況2條件下塔筒頂部節(jié)點(diǎn)Y向振幅隨激振頻率變化

由圖13可知,在激振力頻率為0.5419Hz時(shí),塔筒頂部節(jié)點(diǎn)位移最大,最大位移為0.103304m。由相位角-90.7555°可判斷此時(shí)已基本與共振峰值點(diǎn)重合。進(jìn)一步觀察0.5419Hz下塔筒共振時(shí)整體位移和應(yīng)力分布情況,如圖14和圖15所示。

圖14 工況2條件下塔筒共振時(shí)應(yīng)力分布/Pa

圖15 工況2條件下塔筒共振時(shí)位移分布/m

由圖14和圖15可知,共振時(shí)實(shí)部計(jì)算結(jié)果為:最大應(yīng)力3.43MPa,最大位移0.001636m。虛部計(jì)算結(jié)果為:最大應(yīng)力25.2MPa,最大位移0.107364m。

綜合實(shí)部、虛部計(jì)算結(jié)果可知,塔筒產(chǎn)生的最大應(yīng)力為25.43MPa,位于壁厚發(fā)生突跳的焊縫處。塔筒頂部最大振幅為0.107376m。由此可知,工況2條件下風(fēng)繞流塔筒產(chǎn)生的橫向激振力對(duì)塔筒的影響很小,幾乎不對(duì)塔筒造成危害。

3.4 所有工況諧響應(yīng)計(jì)算結(jié)果

由于工況3和工況4計(jì)算過程與工況2類似,不再詳述過程。工況3條件下塔筒產(chǎn)生的最大應(yīng)力為131.535MPa,塔筒頂部最大振幅為0.5027m。工況4塔筒產(chǎn)生的最大應(yīng)力為265.55MPa,塔筒頂部最大振幅為0.5075m。塔筒仍符合安全性要求。

4 塔筒安裝建議

塔筒安裝主要關(guān)心高于多少風(fēng)速時(shí)不能繼續(xù)實(shí)施吊裝,主要考慮風(fēng)激振頻率與已安裝塔筒固有頻率是否會(huì)發(fā)生共振。此時(shí)研究的工況主要為工況2、3和4,因工況2條件下塔筒與風(fēng)激振力共振時(shí)影響很小,因此忽略此工況。圖16為激振力頻率隨風(fēng)速的分布。

圖16 激振力頻率隨風(fēng)速分布

風(fēng)力機(jī)吊裝能否實(shí)施取決于風(fēng)激振與這兩種工況塔筒固有頻率是否會(huì)發(fā)生共振,計(jì)算結(jié)果見表6。

表6 吊裝工況下塔筒振動(dòng)避開率

從表6可看出,塔筒在風(fēng)速為15m/s時(shí)不滿足避開率大于5%的要求,按照GL規(guī)范要求取安全系數(shù)1.1,此時(shí)吊裝風(fēng)速應(yīng)為14/1.1=12.73m/s。相關(guān)文獻(xiàn)[8,9]風(fēng)速大于12m/s時(shí),風(fēng)力機(jī)應(yīng)停止吊裝,以保證工程安全。結(jié)合以上計(jì)算結(jié)果,建議當(dāng)風(fēng)速大于12m/s時(shí)FD82E風(fēng)力機(jī)停止吊裝。

5 結(jié)論

通過對(duì)FD82E塔筒渦激振動(dòng)、固有頻率、諧響應(yīng)等方面分析計(jì)算,得出了以下結(jié)論:

(1)塔筒附近流動(dòng)符合雷諾相似準(zhǔn)則,升力系數(shù)、阻力系數(shù)和斯特勞哈爾數(shù)等可近似認(rèn)為僅與雷諾數(shù)相關(guān)。

(2)塔筒繞流的CFD計(jì)算中,采用二維方法比三維方法計(jì)算出的升力系數(shù)偏高,而阻力系數(shù)和斯特勞哈爾數(shù)基本相同,因而二維計(jì)算方法所得載荷更趨保守。

(3)采用轉(zhuǎn)捩計(jì)算模型可以得到與阻力系數(shù)實(shí)驗(yàn)值趨勢(shì)符合良好的結(jié)果,但轉(zhuǎn)捩模型計(jì)算值比全湍流模型計(jì)算值偏低。采用全湍流模型計(jì)算可得到更加保守的載荷。

(4)表面粗糙度、來流湍流度等對(duì)塔筒的升力系數(shù)、阻力系數(shù)有較大影響,綜合考慮這些因素,在實(shí)際計(jì)算中建議乘以修正系數(shù),本文建議取值為1.5。

(5)在塔筒正常運(yùn)行的工況1中,塔筒的固有頻率和風(fēng)輪振動(dòng)頻率滿足GL規(guī)范要求的5%避開率。由工況1條件下的諧響應(yīng)分析結(jié)果可知,風(fēng)力機(jī)正常運(yùn)行時(shí)塔筒共振頻率附近的風(fēng)繞流激振力很小 (僅有不到4kN),對(duì)塔筒影響也很小,可以忽略。

(6)風(fēng)繞流塔筒諧響應(yīng)計(jì)算中,當(dāng)激振力與塔筒產(chǎn)生共振時(shí),塔筒位移與激振力的相位角相差約π/2。由工況2、3、4條件下諧響應(yīng)分析結(jié)果可知,風(fēng)繞流塔筒產(chǎn)生激振力對(duì)塔筒影響較小(最大應(yīng)力265MPa,小于許用應(yīng)力336MPa),不會(huì)產(chǎn)生應(yīng)力破壞;

(7)通過分析計(jì)算不同風(fēng)速下風(fēng)繞流塔筒激振力與塔筒的共振情況,得出當(dāng)風(fēng)速大于12m/s時(shí),風(fēng)力機(jī)應(yīng)該停止吊裝。

[1]詹昊,李萬平,方秦漢,等.不同雷諾數(shù)下圓柱繞流仿真計(jì)算 [J].武漢理工大學(xué)學(xué)報(bào),2008,30(12):129-132

[2]謝宇新.馳振穩(wěn)定性分析及其工程應(yīng)用 [D].天津大學(xué):碩士學(xué)位論文,2003

[3]陳桂彬,楊超,鄒叢青.氣動(dòng)彈性設(shè)計(jì)基礎(chǔ) [M].北京:北京航空航天大學(xué)出版社,2010

[4]王亞玲,劉應(yīng)中,繆國平.圓柱繞流的三維數(shù)值模擬 [J].上海交通大學(xué)學(xué)報(bào),2001,35(10):1464-1469

[5]GL規(guī)范.GL Wind 2003

[6]DIN1055-4德國風(fēng)載荷標(biāo)準(zhǔn),NABau 2005

[7]何鴻濤.圓柱繞流及其控制的數(shù)值模擬研究 [D].北京交通大學(xué):碩士學(xué)位論文,2009

[8]廖明夫,R.Gasch,J.Twele.風(fēng)力發(fā)電技術(shù) [M].西安:西北工業(yè)大學(xué)出版社,2009

[9]何顯富,盧霞,楊躍進(jìn),劉萬琨.風(fēng)力機(jī)設(shè)計(jì)、制造與運(yùn)行 [M].北京:化學(xué)工業(yè)出版社,2009

猜你喜歡
模態(tài)模型
一半模型
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計(jì)的漸近分布
3D打印中的模型分割與打包
車輛CAE分析中自由模態(tài)和約束模態(tài)的應(yīng)用與對(duì)比
國內(nèi)多模態(tài)教學(xué)研究回顧與展望
FLUKA幾何模型到CAD幾何模型轉(zhuǎn)換方法初步研究
高速顫振模型設(shè)計(jì)中顫振主要模態(tài)的判斷
基于HHT和Prony算法的電力系統(tǒng)低頻振蕩模態(tài)識(shí)別
由單個(gè)模態(tài)構(gòu)造對(duì)稱簡支梁的抗彎剛度
主站蜘蛛池模板: 97亚洲色综久久精品| 91亚洲精选| 大陆精大陆国产国语精品1024| 91美女视频在线观看| 天堂成人av| 免费观看亚洲人成网站| 在线观看国产精品第一区免费| 欧美精品亚洲精品日韩专| 国产一区亚洲一区| 国产精品19p| 亚洲国产精品国自产拍A| 欧美成人区| 亚洲综合天堂网| 亚洲成人免费在线| 天天色综合4| 欧美一级片在线| 国产成人亚洲综合a∨婷婷| 一级全免费视频播放| 91麻豆国产视频| 毛片大全免费观看| 99久久精品国产综合婷婷| 97久久精品人人| 在线中文字幕日韩| 精品91视频| h网址在线观看| 99视频只有精品| 国产成人h在线观看网站站| 日韩欧美国产精品| 无码丝袜人妻| 91色爱欧美精品www| 99久久性生片| 免费国产高清精品一区在线| 国产精品亚洲一区二区三区在线观看| 毛片基地视频| 日本一区二区三区精品国产| 天堂va亚洲va欧美va国产 | 国产一区二区三区免费| 欧美三级自拍| 亚洲制服中文字幕一区二区| 国产亚洲成AⅤ人片在线观看| 亚洲无码视频一区二区三区| 亚洲不卡影院| 欧美国产日本高清不卡| 天天爽免费视频| 夜夜操天天摸| 免费看a级毛片| 国产免费黄| 欧美三級片黃色三級片黃色1| 免费一级大毛片a一观看不卡| 国产AV毛片| 免费一级无码在线网站| 欧美亚洲国产一区| 亚洲视频三级| 国产爽歪歪免费视频在线观看 | 欧美精品H在线播放| 亚洲二三区| 国产女人在线观看| 国产精品无码制服丝袜| 91在线精品免费免费播放| 一区二区影院| 999国内精品久久免费视频| 国产欧美又粗又猛又爽老| 国产成人精品一区二区三在线观看| 成人一级免费视频| 波多野结衣一区二区三区四区| 中国精品久久| 久久精品人妻中文系列| 日韩在线播放中文字幕| 免费国产福利| 国产一级在线观看www色| 国产欧美精品一区二区| 国内精品久久人妻无码大片高| 亚洲性视频网站| 2021精品国产自在现线看| 国产亚洲一区二区三区在线| 免费一级毛片不卡在线播放| 久久香蕉国产线看观看精品蕉| 91精品最新国内在线播放| 亚洲国产精品美女| 国产一区三区二区中文在线| 欧美一区二区精品久久久| 国产欧美日本在线观看|