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ì)稱簡支梁的抗彎剛度
主站蜘蛛池模板: 日本成人在线不卡视频| 美美女高清毛片视频免费观看| 亚洲中字无码AV电影在线观看| 麻豆AV网站免费进入| 欧美色视频网站| 亚洲综合精品香蕉久久网| 日本欧美视频在线观看| 香蕉99国内自产自拍视频| 中国国产一级毛片| 另类专区亚洲| 国产在线精品人成导航| 亚洲精品动漫| 99久久精品国产综合婷婷| 亚洲天堂日韩在线| 人妻精品全国免费视频| 青青草原偷拍视频| 国产精彩视频在线观看| 五月婷婷丁香色| 国产乱人激情H在线观看| 99久久无色码中文字幕| 精品少妇人妻无码久久| 欧美性精品不卡在线观看| 免费看的一级毛片| 亚洲 成人国产| 国产区人妖精品人妖精品视频| 亚洲无线观看| 国产毛片高清一级国语 | 国产微拍一区| 国产一级做美女做受视频| 毛片网站免费在线观看| 91精品日韩人妻无码久久| 久久久久久久久18禁秘| 五月天在线网站| 日韩精品资源| 国产成人高清精品免费5388| 亚洲色图欧美| 欧洲亚洲一区| 国产欧美自拍视频| 日韩精品一区二区深田咏美| 成人精品在线观看| 真实国产精品vr专区| 老司机久久99久久精品播放| 欧美有码在线| 久久人妻xunleige无码| 无码精品国产VA在线观看DVD| 国产精品深爱在线| 成年人国产网站| 亚洲免费毛片| 18禁影院亚洲专区| 1769国产精品免费视频| 天天做天天爱夜夜爽毛片毛片| 国产小视频a在线观看| 亚洲激情区| 国产粉嫩粉嫩的18在线播放91| 久久亚洲中文字幕精品一区| 亚洲bt欧美bt精品| 精品一区二区三区水蜜桃| 国禁国产you女视频网站| 亚洲黄网视频| 国产欧美视频综合二区| 国产麻豆91网在线看| 亚洲一区无码在线| 波多野结衣久久精品| 国产精品一老牛影视频| 亚洲美女一区二区三区| 再看日本中文字幕在线观看| 亚洲精品不卡午夜精品| 国产亚洲视频中文字幕视频| 日本91视频| 久久久久国产精品熟女影院| 日本三区视频| 曰韩人妻一区二区三区| 2021国产精品自产拍在线观看| 不卡网亚洲无码| 婷五月综合| 久久久久国产一级毛片高清板| 在线另类稀缺国产呦| 欧美亚洲中文精品三区| www.国产福利| 国产 日韩 欧美 第二页| 熟女成人国产精品视频| 久久综合色天堂av|