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

海上風力機多場平臺構(gòu)建及地震動力學研究

2022-05-18 08:25:44閆陽天繆維跑岳敏楠牛凱倫
動力工程學報 2022年5期
關(guān)鍵詞:結(jié)構(gòu)模型

閆陽天, 李 春,2, 繆維跑, 岳敏楠, 牛凱倫

(1. 上海理工大學 能源與動力工程學院, 上海 200093;2. 上海市動力工程多相流動與傳熱重點實驗室, 上海 200093)

針對海上風力機遭受地震及多種環(huán)境載荷作用的結(jié)構(gòu)安全問題,美國和歐洲海上風電發(fā)展較為迅速的國家和地區(qū)已制定了相關(guān)規(guī)范。挪威船級社(DNV)、丹麥Risφ國家實驗室、德國勞氏船級社(GL)[1]、美國風能協(xié)會(AWEA)[2]及美國土木工程協(xié)會(ASCE)[3]分別針對各國特有的氣候、地質(zhì)地理條件制定了相關(guān)設(shè)計規(guī)范。規(guī)范中主要涉及3個方向:計算方法、結(jié)構(gòu)模型及風-震-波浪場耦合模型。

海上風力機會遭受不同大小的湍流風、波浪及海流載荷。針對環(huán)境載荷多變的海上風力機,國內(nèi)外諸多學者采用不同的耦合模型對風力機正常運行時遭受突發(fā)性地震沖擊下的結(jié)構(gòu)響應(yīng)展開研究。

早期學者們基于建筑規(guī)范中常用的地震響應(yīng)譜分析法計算風力機地震響應(yīng)[4-5]。Ma等[6]對比了瞬態(tài)分析法與響應(yīng)譜分析法計算所得地震作用下的塔架結(jié)構(gòu)響應(yīng),發(fā)現(xiàn)響應(yīng)譜分析結(jié)果較瞬態(tài)分析結(jié)果更劇烈,按照建筑規(guī)范程序設(shè)計風力機基礎(chǔ)會導致對地震載荷評估過于保守。Witcher[7]認為針對傳統(tǒng)建筑地震評估采用的頻域法不能解決地震載荷與氣動載荷同時作用的氣彈問題。Santangelo等[8]認為響應(yīng)譜分析法在解決氣彈和土壤柔性這2個關(guān)鍵問題上存在缺陷。楊陽等[9-10]基于FAST軟件建立湍流風-地震時域耦合模型,模擬風力機在不同土質(zhì)和地震強度下的塔架動力學響應(yīng)。上述研究基于FAST構(gòu)建了考慮氣動載荷、波浪載荷及土-構(gòu)耦合(SSI)效應(yīng)的多體動力學耦合模型,且運算準確高效,但計算結(jié)果比較宏觀,針對結(jié)構(gòu)局部承載規(guī)律研究能力薄弱。有限元方法可以有效彌補這一不足。若采用有限元殼模型進行研究,基于計算流體動力學-計算結(jié)構(gòu)動力學(CFD-CSD)的全耦合時域方法計算昂貴且低效。Santangelo等[8]采用非耦合時域法降低計算消耗量,研究了5 MW風力機在地震及風載荷作用下的結(jié)構(gòu)響應(yīng),發(fā)現(xiàn)與全耦合仿真結(jié)果相比,時域非耦合分析亦可得到較為準確的結(jié)果。

風、浪載荷通常是海上風力機結(jié)構(gòu)設(shè)計中最重要的環(huán)境載荷,因此在研究風力機結(jié)構(gòu)地震響應(yīng)時有必要考慮[11]。目前的研究方法主要依賴海上風力機的準靜態(tài)載荷計算[8]。Sun等[12]研究風力機在風、波和地震載荷下結(jié)構(gòu)響應(yīng)和振動控制,其中對于風載荷,采用葉素動量理論計算作用在旋轉(zhuǎn)葉片上的氣動載荷,固定樁基承受的波浪載荷采用莫里森方程計算。李穎等[13]基于FAST探討了研究地震響應(yīng)時考慮風、浪載荷的必要性。Wang等[14]研究了風、浪、流和地震聯(lián)合作用下的結(jié)構(gòu)響應(yīng),結(jié)果表明土-構(gòu)耦合效應(yīng)及多種環(huán)境載荷缺一不可。但在研究海上風力機地震動力學時,針對是否考慮環(huán)境載荷的必要性探討尚顯不足。

針對上述問題,筆者以10 MW單樁海上風力機為對象,采用有限元方法構(gòu)建海上風力機非線性Winkler土-結(jié)構(gòu)耦合模型,并在此基礎(chǔ)上建立了復雜的土-樁-塔架非線性殼模型。通過對模態(tài)計算結(jié)果的校核,確保構(gòu)建有限元模型的有效性,并基于上述模型進行從局部到整體的瞬態(tài)分析,研究其在風-浪-海流載荷作用下的地震動力學響應(yīng),為大型海上風力機結(jié)構(gòu)設(shè)計和優(yōu)化提供參考。

1 研究對象

2012年,丹麥科技大學(DTU)與丹麥Vestas合作設(shè)計了一款10 MW及以上容量風力機[15-16]。筆者以DTU 10 MW風力機為研究對象,假設(shè)風力機安裝海域平均水深20 m;埋土樁基38 m,根據(jù)不同土壤性質(zhì)將其分為非等長的6段;塔架壁厚隨塔高線性減小;塔基外徑為7 665 mm,壁厚為34 mm;塔頂外徑為5 500 mm,壁厚為26 mm。塔架選用密度為7 850 kg/m3、彈性模量為210 GPa、泊松比為0.3、材料屈服極限為380 MPa的A709圓截面高強鋼。考慮到結(jié)構(gòu)油漆、法蘭、螺栓和焊接質(zhì)量,計算密度修正為8 500 kg/m3。

1.1 有限元方法及模型

在求解瞬態(tài)問題時,采用有限差分法,通過對動力作用和反應(yīng)過程的時間離散來求解瞬態(tài)問題。對于具有線性阻尼和剛度的單自由度系統(tǒng),控制彈簧-阻尼-質(zhì)量系統(tǒng)運動的動力學方程為:

(1)

解決結(jié)構(gòu)系統(tǒng)動力學響應(yīng)的最常用方法是動力平衡方程的直接數(shù)值積分方法。求解微分方程的方法較多,從這些方法獲得的解被假定為具有連續(xù)高階導數(shù)的光滑函數(shù),但是非線性結(jié)構(gòu)的精確解要求加速度和位移的二階導數(shù)不是光滑函數(shù)。因此,針對加速度不連續(xù)性的求解,Newmark開發(fā)了一系列單步積分方法,用于解決爆炸和地震載荷下的結(jié)構(gòu)動力學問題,稱為Newmark時間積分方法。

對風力機的支撐結(jié)構(gòu)建立3D有限元模型。采用殼單元模型,并針對連接位置局部網(wǎng)格進行細化,具體支撐結(jié)構(gòu)的有限元模型如圖1所示。

圖1 支撐結(jié)構(gòu)及有限元模型

1.2 土-構(gòu)耦合模型

建立土-構(gòu)耦合模型時,結(jié)合p-y曲線法與Winkler土-結(jié)構(gòu)耦合模型[17],采用非線性彈簧描述埋土樁基與土壤間的相互作用,基于p-y(土壤反力-橫向位移)曲線及Q-z(土壤反力-縱向位移)曲線確定橫向、縱向的非線性彈簧剛度。土壤參數(shù)如表1所示,其中E為彈性模量,γ為有效土重度,fc為單位表面壓力,ft為表面拉力,q為樁底單位阻力。

景觀娛樂用水區(qū),指以滿足景觀、療養(yǎng)、度假和娛樂需要為目的的江河湖泊等水域。景觀娛樂活動不得危及景觀娛樂用水區(qū)的水質(zhì)控制目標。

表1 土壤參數(shù)

p-y曲線是非線性的,并且沿單樁深度(H)變化。API規(guī)范中根據(jù)最終的橫向承載力pu來得到p-y關(guān)系,這取決于該深度的土壤屬于淺層土還是深層土,深淺的劃分由樁基外徑及土壤力學性質(zhì)共同決定。砂性土的p-y曲線表達式為:

式中:A為經(jīng)驗調(diào)整系數(shù),A=3.0-0.8H/D;y為橫向位移;m;pu為樁深至H時的極限承載力,kN/m;pus為淺層土時的極限承載力,kN/m;pud為深層土時的極限承載力,kN/m;k為地基反力系數(shù),kN/m3;D為樁基外徑,m;C1、C2與C3為系數(shù),由內(nèi)摩擦角φ確定[17],分別取3.2、3.6和60。

經(jīng)計算,不同土壤深度處的p-y曲線如圖2所示。

圖2 不同土壤深度處的p-y曲線

2 多場耦合模型構(gòu)建

筆者構(gòu)建了考慮湍流風、波浪、海流及地震的多場非全耦合模型,見圖3,其中K1~K7為非線性彈簧剛度。將250 m×236 m的湍流風風場通過以輪轂為中心的13×13個節(jié)點進行計算,每根葉片由37個質(zhì)點的多體模型描述,氣動載荷采用FAST中的AeroDyn求解。波浪和海流載荷通過在AQWA中定義P-M譜進行求解。地震載荷將實測地表位移響應(yīng)作用于非線性土-構(gòu)耦合模型,地震攜帶能量通過土-構(gòu)耦合模型傳至樁基。

圖3 多場耦合模型

2.1 湍流風

建立湍流風風場模型的關(guān)鍵是選取合適的風譜模型,筆者采用IEC Kaimal風譜模型,湍流風載荷的計算結(jié)果如圖4所示,可以看出因變槳效應(yīng)導致切出風速時結(jié)構(gòu)承受的風載荷更低。

圖4 湍流風載荷時域圖

2.2 波浪與海流

使用Morison方程可以有效計算出振蕩流中作用在圓柱體上的力。該方程將單位長度的總力描述為阻力和慣性力的總和:

(6)

2.3 地震運動

圖5 各地震載荷x、y和z方向位移譜

(a) 監(jiān)測站CHY036

3 結(jié)論與分析

為更準確地研究地震響應(yīng)時湍流風、波浪及海流載荷同時作用對結(jié)構(gòu)的影響,針對不同地震發(fā)生時湍流風、波浪及海流載荷的作用構(gòu)建13個工況,計算載荷如表2所示,考慮波浪載荷的工況時,海流載荷固定為4 m/s。總計算時長為100 s,載荷步為0.02 s,根據(jù)瞬態(tài)交變載荷增量的不同,將各載荷步劃分子步,每個算例計算步長不超過15 000步,每個算例計算時長約28 h。

表2 計算載荷

分析各種環(huán)境載荷變化對地震動力學響應(yīng)的影響。在時域內(nèi)給出了整個風力發(fā)電系統(tǒng)的塔頂位移和應(yīng)力應(yīng)變能。為了研究不同高度時沿塔身周向的應(yīng)力分布,按照圖7所示定義周向探針分布,設(shè)置湍流風來流方向為0°探針,每間隔30°布置一個探針。

圖7 塔架周向應(yīng)力探針分布

3.1 地震載荷

針對變地震載荷條件下結(jié)構(gòu)塔頂位移、結(jié)構(gòu)應(yīng)力及應(yīng)變能進行時域響應(yīng)分析,結(jié)果如圖8所示。

(a)

從圖8可以看出,工況1~工況4為4種變地震載荷工況,地震均發(fā)生自風力機啟動運行后的第40 s,各對比算例的時域響應(yīng)在前40 s保持一致,由無地震發(fā)生的工況1可知,風力機啟動的前20 s各響應(yīng)波動幅度較大,隨后隨湍流風風速的變化而變化,但整體趨于平穩(wěn)。隨著地震的發(fā)生,結(jié)構(gòu)各時域響應(yīng)產(chǎn)生不同程度的增幅,塔頂位移變化更為顯著。由圖8(a)~圖8(c)可知,塔頂在x方向發(fā)生偏心往復運動,風、浪載荷造成塔頂位移在y和z方向位移分量遠小于地震載荷。風、浪載荷雖然方向與大小不穩(wěn)定,但主要以來流方向為主。而地震載荷各方向分量大小由橫波、縱波的大小決定,無主流方向。x和y方向的位移增幅明顯大于z方向,說明結(jié)構(gòu)橫向遭受地震載荷沖擊對地震載荷更敏感。工況2~工況4選取的地震PGD依次降低,但工況4所受應(yīng)變能與應(yīng)力峰值最大,工況2其次。結(jié)合圖6發(fā)現(xiàn),產(chǎn)生這一現(xiàn)象的原因是與工況2相比,工況4的PGD數(shù)值較小,但持續(xù)時間長,且近地位移響應(yīng)在x和y方向分量較大。工況4的位移峰值出現(xiàn)在80~90 s,原因是在該時刻結(jié)構(gòu)因風、浪載荷導致局部響應(yīng)峰值,因此在風力機結(jié)構(gòu)地震研究時,需重視橫波導致的結(jié)構(gòu)橫向響應(yīng),且需要同時考慮風、浪載荷。

對工況 2~工況 4地震發(fā)生后結(jié)構(gòu)應(yīng)力達到峰值時的表面應(yīng)力依據(jù)圖7的思路進行提取,結(jié)果見圖9,工況后面給出的是應(yīng)力峰值時間。應(yīng)力響應(yīng)達到峰值時支撐結(jié)構(gòu)壁面應(yīng)力在塔高方向的分布均呈現(xiàn)自上而下的非線性遞增現(xiàn)象。湍流風來流方向為0°,即探針Path0的位置。因此,此刻應(yīng)力峰值出現(xiàn)在Path0與Path6,即0°與180°附近。由圖9(a)可以看出,該現(xiàn)象是由強勁的湍流風載荷導致的,隨著地震載荷的加入,以工況4為例,由圖9(d)可知應(yīng)力最大區(qū)域在周向上分布的最大位置在0°~50°及200°~250°,峰值區(qū)域向逆時針方向轉(zhuǎn)動。

(a) 工況1-83.62 s

圖10為變地震載荷條件下峰值時刻塔架位移云圖,其中最左邊的圖是未發(fā)生地震載荷時的位移云圖。圖11為變地震載荷條件下峰值時刻應(yīng)變能云圖。

圖11 變地震載荷塔架應(yīng)變能云圖

結(jié)合峰值瞬間支撐結(jié)構(gòu)的位移及應(yīng)變能進行分析。風力機在遭受湍流風及地震作用時,結(jié)構(gòu)產(chǎn)生的響應(yīng)以一階模態(tài)為主,主要為塔頂?shù)那昂蠹皺M向位移。應(yīng)變能集中在26 m左右塔高位置。從圖10可以看出,風力機啟動造成的塔頂位移幅值高于地震及湍流風所造成的位移響應(yīng)幅值。湍流風作用對結(jié)構(gòu)塔頂運動起主導作用,地震造成塔架底部應(yīng)變能的聚集,對比工況 1與工況 4發(fā)現(xiàn),地震造成的塔架位移及應(yīng)變能增幅相較于風、浪造成的結(jié)構(gòu)位移及應(yīng)變能增幅較小,所以在地震發(fā)生時控制風、浪載荷可實現(xiàn)減少地震威脅結(jié)構(gòu)安全的可能。

圖10 變地震載荷塔架位移云圖

3.2 湍流風載荷

依據(jù)位移、應(yīng)力及應(yīng)變能的時域和瞬間響應(yīng),分析變湍流風載荷工況的影響。

圖12為相同地震載荷及波浪載荷作用下,3種變湍流風載荷工況的結(jié)構(gòu)時域響應(yīng),所選取的工況 2、工況 9和工況 10分別為額定風速(11.4 m/s)、無風載荷及切出風速(25.0 m/s)作用下支撐結(jié)構(gòu)時域響應(yīng)。從圖12可以看出,0~40 s時,因采用變槳系統(tǒng),額定風速下結(jié)構(gòu)時域響應(yīng)幅值最高,說明結(jié)構(gòu)承受的風載荷也更劇烈。當?shù)卣鸢l(fā)生后,使得結(jié)構(gòu)動力學響應(yīng)瞬間急劇增加,但三者風載荷的不同導致各響應(yīng)增幅差異較大。工況9比另外2個工況的響應(yīng)波動頻率更高。根據(jù)上述現(xiàn)象發(fā)現(xiàn),地震載荷造成的結(jié)構(gòu)響應(yīng)頻率較風載荷更高,額定風速下結(jié)構(gòu)承受風載荷更強,反而與地震載荷同時作用產(chǎn)生的結(jié)構(gòu)動力學響應(yīng)更小,說明湍流風載荷產(chǎn)生了氣彈阻尼效應(yīng),導致結(jié)構(gòu)受地震沖擊產(chǎn)生的響應(yīng)幅值增加,加速了地震傳遞至結(jié)構(gòu)的能量耗散。

(a)

圖13給出了變湍流風載荷時的周向應(yīng)力分布。可以看出,工況 2、工況 9和工況 10下的應(yīng)力響應(yīng)分布規(guī)律基本一致,結(jié)構(gòu)應(yīng)力產(chǎn)生峰值的時間較為接近,因此支撐結(jié)構(gòu)應(yīng)力達到峰值的時間主要取決于地震載荷。工況10應(yīng)力幅值明顯大于工況2和工況9,而工況9與工況2的應(yīng)力幅值較為接近。高風速下,結(jié)構(gòu)承受湍流風載荷較小,低風速下反之。但地震發(fā)生時,高的風載荷起到的氣彈阻尼效應(yīng)可有效降低結(jié)構(gòu)應(yīng)力的激增,基于此可開展有效的結(jié)構(gòu)抗震設(shè)計。

(a) 工況2-61.60 s

圖14給出了變湍流風載荷條件下峰值時刻塔架位移云圖,其中左邊的3個圖是未發(fā)生地震載荷時的位移云圖。圖15給出了變湍流風載荷條件下峰值時刻應(yīng)變能云圖。

圖14 變湍流風載荷塔架位移云圖

對塔架在地震發(fā)生前后的峰值位移云圖進行對比分析,并針對應(yīng)變能峰值進行研究。地震發(fā)生前,工況2的峰值位移明顯高于工況9和工況10,說明低風速湍流風造成的結(jié)構(gòu)位移響應(yīng)更劇烈;地震發(fā)生后,工況2的峰值位移響應(yīng)減小,而工況 9和工況10的峰值位移響應(yīng)明顯增加。工況10的峰值位移響應(yīng)明顯高于工況2和工況9,說明高氣動載荷造成的氣彈阻尼效應(yīng)可有效減小地震載荷導致的塔架位移響應(yīng)幅值。由圖15可知,氣彈阻尼效應(yīng)同樣可以有效緩解應(yīng)變能堆積,加速結(jié)構(gòu)遭受地震載荷后蓄積于結(jié)構(gòu)形變中的能量耗散。

圖15 變湍流風載荷塔架應(yīng)變能云圖

3.3 波浪載荷

對變波浪載荷結(jié)構(gòu)時域響應(yīng)結(jié)果進行對比分析,結(jié)果如圖16所示。從圖16可以看出,沒有地震發(fā)生時,波浪載荷的變化對結(jié)構(gòu)產(chǎn)生影響較小,湍流風載荷對結(jié)構(gòu)各響應(yīng)起主導作用。地震發(fā)生后,波浪載荷變化導致地震作用下結(jié)構(gòu)動力學響應(yīng)增加。工況2與工況13各響應(yīng)相差較小,說明波浪載荷較小時,可以忽略不計,但隨著波浪載荷的加劇,地震作用下的結(jié)構(gòu)響應(yīng)加劇。

(a)

分析圖17可知,各工況應(yīng)力峰值時刻并無明顯差異,但與其他2個算例相比,工況 6的應(yīng)力值有明顯變化,其節(jié)點應(yīng)力峰值為236.4 MPa;而工況 2和工況 13分別為191.2 MPa和196.6 MPa,相差約23.6%。考慮到不同支撐結(jié)構(gòu)對波浪載荷的敏感性不同,以及海水深度的增加所造成的影響,對固定式海上風力機各結(jié)構(gòu)制定詳細的標準時,波浪載荷不可忽略。

(a) 工況2-61.60 s

圖18給出了變波浪載荷條件下峰值時刻塔架位移云圖,其中左邊3個圖是未發(fā)生地震載荷時的位移云圖。圖19給出了變波浪載荷條件下峰值時刻應(yīng)變能云圖。由圖18和圖19可知,塔架的位移響應(yīng)依然保持以一階模態(tài)為主的塔頂俯仰運動,而波浪載荷的變化僅對工況 6在61 s峰值時刻的位移產(chǎn)生一定影響,但影響較小,說明波浪載荷的變化化對結(jié)構(gòu)位移影響較小。應(yīng)變能蓄積位置依然保持在塔底向上26 m的位置附近,但幅值分布差異十分明顯,工況6的應(yīng)變能幅值高于1 980 J的區(qū)域明顯大于工況2和工況13,說明波浪載荷對塔基應(yīng)變能影響較大。因此,波浪載荷對結(jié)構(gòu)的影響集中于結(jié)構(gòu)底部受波浪及海流沖擊較近的位置。地震能量由樁基自下而上傳遞,所以波浪載荷變化在一定程度上影響地震能量自樁基向結(jié)構(gòu)上部傳遞的過程,從而導致結(jié)構(gòu)動力學響應(yīng)的變化。

圖18 變波浪載荷塔架位移云圖

圖19 變波浪載荷塔架應(yīng)變能云圖

4 結(jié) 論

(1) 由于風、浪、流載荷的瞬態(tài)特性,風力機結(jié)構(gòu)的動態(tài)響應(yīng)隨環(huán)境載荷變化波動較大。支撐結(jié)構(gòu)動態(tài)響應(yīng)規(guī)律主要呈現(xiàn)塔頂前后及橫向的往復運動。

(2) 風、浪載荷導致結(jié)構(gòu)橫向和縱向各響應(yīng)增幅小于地震載荷。結(jié)構(gòu)橫向響應(yīng)受風、浪、流載荷和氣彈效應(yīng)影響更顯著。地震沖擊下,各方向響應(yīng)迅速增加,同時會放大風、浪、流載荷對結(jié)構(gòu)的作用。

(3) 在承受風、浪、流及地震載荷時,結(jié)構(gòu)塔基處應(yīng)力激增和應(yīng)變能的蓄積說明此處結(jié)構(gòu)失效風險較高。

(4) 不應(yīng)僅從塔頂或塔基的宏觀時域響應(yīng)判斷是否有考慮該環(huán)境載荷的必要性。動力學響應(yīng)分析應(yīng)考慮氣彈效應(yīng)對結(jié)構(gòu)位移響應(yīng)的影響和浪、流載荷引起的支撐結(jié)構(gòu)局部應(yīng)力增幅。湍流風載荷和浪、流載荷在研究風力機地震動力學響應(yīng)時不可忽視。

猜你喜歡
結(jié)構(gòu)模型
一半模型
《形而上學》△卷的結(jié)構(gòu)和位置
哲學評論(2021年2期)2021-08-22 01:53:34
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計的漸近分布
論結(jié)構(gòu)
中華詩詞(2019年7期)2019-11-25 01:43:04
新型平衡塊結(jié)構(gòu)的應(yīng)用
模具制造(2019年3期)2019-06-06 02:10:54
論《日出》的結(jié)構(gòu)
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉(zhuǎn)換方法初步研究
創(chuàng)新治理結(jié)構(gòu)促進中小企業(yè)持續(xù)成長
主站蜘蛛池模板: 亚洲国产欧美国产综合久久 | 手机在线国产精品| 国产一区二区精品福利| 一级做a爰片久久免费| 国产精品黑色丝袜的老师| 日本久久网站| 国内a级毛片| 国产性爱网站| 欧美一区国产| 无码AV日韩一二三区| 亚洲91精品视频| 制服无码网站| 99国产精品一区二区| 91精品免费高清在线| 亚洲欧美成aⅴ人在线观看| 国产手机在线ΑⅤ片无码观看| 2021国产v亚洲v天堂无码| 亚洲Va中文字幕久久一区| 亚洲天堂免费在线视频| 亚洲AⅤ永久无码精品毛片| 欧美亚洲日韩不卡在线在线观看| 国产精品免费福利久久播放| 精品欧美一区二区三区在线| 国产尹人香蕉综合在线电影| 欧美中文一区| 国产xx在线观看| 一本久道热中字伊人| 72种姿势欧美久久久久大黄蕉| 黄片一区二区三区| 亚洲欧美在线看片AI| 亚洲熟女中文字幕男人总站| 国产91精品久久| 欧美激情视频一区| 91美女视频在线观看| 国产一区二区三区视频| 午夜毛片免费看| 四虎影视库国产精品一区| 国产美女视频黄a视频全免费网站| 亚洲综合日韩精品| 久久久国产精品免费视频| 中文字幕有乳无码| 超碰免费91| 欧美成人日韩| 1769国产精品免费视频| 69视频国产| 精品国产黑色丝袜高跟鞋| 亚洲嫩模喷白浆| 蜜桃视频一区| 日本欧美中文字幕精品亚洲| 亚洲精品第一在线观看视频| 欧美午夜视频| 欧美一级片在线| www.亚洲国产| 亚洲一区波多野结衣二区三区| 国产一二三区视频| 亚洲码一区二区三区| 自拍偷拍一区| 国产精品免费福利久久播放| 97视频在线观看免费视频| 9cao视频精品| 黄色网站在线观看无码| 欧美国产日韩在线观看| 久久久久久久久久国产精品| 日韩中文精品亚洲第三区| 国产第八页| 亚洲精品制服丝袜二区| 精品国产成人av免费| 日韩小视频在线播放| 亚洲天堂在线免费| 日本道中文字幕久久一区| 91在线视频福利| 久久久久久国产精品mv| 国产精品美女网站| 国产精品亚洲专区一区| 亚洲无码A视频在线| 国产丝袜无码一区二区视频| 伦精品一区二区三区视频| a色毛片免费视频| 狠狠色丁香婷婷| 久久久亚洲国产美女国产盗摄| m男亚洲一区中文字幕| 亚洲一区二区三区在线视频|