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

波浪錐型風力俘能結構能量轉換效率

2023-08-30 13:29:20閆豫龍柳迪偉鄭云龍
上海交通大學學報 2023年8期
關鍵詞:振動效率結構

鄒 琳, 閆豫龍, 陶 凡, 柳迪偉, 鄭云龍

(武漢理工大學 機電工程學院,武漢 430079)

傳統的水平軸風力發電機存在結構復雜、安裝維護困難、承受載荷混雜等缺點[1],通常放置在沿海、近海和其他高風速地區,低風速條件下風力發電機轉速會發生驟降,導致功率輸出與部件可靠性降低[2].Adel[3]提出了一款名為Vortex的無葉片風力發電機,利用結構周圍空氣形成的漩渦晃動柱體,繼而將柱體振動機械能轉化為電能.無葉片風力機基于渦激振動的原理,能在低風速環境下有效捕獲漩渦能量,提高對風能的利用效率[4].雖然目前基于流致振動的俘能設備其風能轉換效率仍弱于水平軸風力機,例如Hemon等[5]設計的風力俘能結構,風能轉換效率不到1%;Badhurshah等[6]提出的雙穩態彈簧的渦激振動俘能結構,最大轉換效率為10%,但這類系統因其適應各種情況的潛力而仍然令人感興趣.相較于水平軸風力發電機的高風速要求,無葉片風力發電機則可進行小型化設計并可用于低風速地區分布式風力發電布局,適合作為無線傳感器網絡系統的能源利用和區域小型電力系統[7].增強渦激振動有利于無葉片風力機采集風能并轉換為電能,而質量比、阻尼比等控制參數對無葉片風力機俘能結構的渦激振動能量轉換有重要影響.對此,國內外學者開展了諸多研究.

Zhang等[8]研究表明高阻尼比顯著抑制圓柱俘能結構渦激振動能量轉換效率.Kumar等[9]和李小超等[10]發現圓柱俘能結構出現最佳能量轉換特性的折合流速與阻尼比有關.Zheng等[11]數值研究證明質量比會對渦激振動能量轉換造成影響.Zhao等[12]、Xu等[13]和Wang等[14]均設計出改型圓柱結構的渦激振動俘能裝置,其俘能功率和能量轉換效率較圓柱俘能結構顯著提升,并提出俘獲能量與阻尼比呈正相關.分析上述文獻可以看出,繞流結構的截面形狀對渦激振動俘能裝置的能量轉換也起著不可忽視的作用.Lam等[15]給出合適的波長比、波幅比參數使波浪型圓柱具有明顯的減阻效果.Chizfahm等[16]研究證明圓錐柱型無葉片風力俘能結構在高風速下具有良好的能量轉換性能.鄒琳等[17]在波浪圓柱中引入斜率得到一種波浪錐型圓柱,發現適合的波長比、波幅比、斜率組合可以增強柱體振動.

目前,渦激振動能量轉換的研究大多集中在水生能源圓柱俘能結構,本文提出的波浪錐型風力俘能結構的能量轉換特性尚無人研究.通過數值模擬分析質量比、阻尼比對波浪錐型俘能結構渦激振動特性、能量轉換效率的影響.尋求合適的質量阻尼參數,以增大振幅、提高振動頻率和拓寬鎖頻區間為目的,為波浪錐型風力俘能結構能量轉換效率的提升提供理論支持.

1 數學模型

1.1 波浪錐型俘能結構-發電機-負載耦合模型

俘能結構模型示意圖如圖1所示.圖1(a)為基于渦激振動的波浪錐型無葉片風力發電機物理模型,主要由4個部分組成:① 波浪錐型俘能結構;② 彈簧-阻尼系統;③ 傳動機構,包括連桿、直線軸承、萬向節聯軸器等;④ 直線發電機.波浪錐型俘能結構-發電機-負載耦合模型可以簡化為質量-彈簧-阻尼模型,如圖1(b)所示.圖中:Fflu為流體力;k為彈簧剛度;c,c0,cg分別為結構阻尼、發電機內能損耗阻尼、負載電阻發電阻尼;Fres為發電機對俘能結構的阻力;R0為發電機的內阻;RL為外界負荷電阻.對渦激振動而言,X方向的振幅往往遠小于Y方向振幅,且單自由度運動更便于能量收集,故俘能系統僅被允許在Y方向做單自由度運動.風力作用迫使俘能結構發生Y方向渦激振動,振動機械能經傳動機構傳遞至永磁發電機動子,進而動子切割磁力線產生電能.

圖1 俘能結構模型示意圖Fig.1 Schematic diagram of energy harvesting structure model

耦合模型發生渦激振動時,俘能結構在Y方向受流體力Fflu和發電機阻力Fres,波浪錐型俘能結構的運動方程可由下式表示:

(1)

(2)

(3)

假設俘能結構位移y與發電機動子位移yg存在倍數關系,傳遞效率為α,則俘能結構與發電機動子之間存在關系:

(4)

(5)

俘能結構受到的發電機阻力為

(6)

將式(1)、(2)、(3)、(6)整理化簡得耦合模型運動方程:

(7)

式中:mt為考慮流體附加質量的俘能結構總質量;ct為考慮發電機和傳動機構后的俘能結構總阻尼;F(t)為流體對俘能結構作用力,也被稱為升力.

1.2 能量轉換效率

由于渦激振動具有周期性,俘能結構一個周期T范圍內的平均有效功率如下式表示:

(8)

式中:T為渦激振動周期.于是,俘能結構從流體中獲取的平均有效功率為

(9)

式中:fosc為俘能結構振動頻率;A為俘能結構振幅;Cl為升力系數幅值;φ為相位角,即升力與波浪錐型俘能結構振動位移之間的相位差.

流體的動壓為1/2ρv2,流體掃過波浪錐型俘能結構的面積簡化為DmL,則施加在俘能結構上的力為1/2ρv2DmL.波浪錐型俘能結構的流體功率Pflu為流體作用力與來流速度的乘積:

(10)

故波浪錐柱俘能結構發生渦激振動時對流體能量的轉換效率為

(11)

2 計算模型及驗證

2.1 控制方程

將波浪錐型俘能結構的渦激振動視為非定常和不可壓縮過程.研究發現,SSTk-ω湍流模型能很好地預測反向壓力梯度流動,在總體預測能力方面似乎更好,結果更加真實[20].采用近壁面的k-ω模型和遠壁面的k-ε模型,用黏度限制器構建湍流黏度μt是SSTk-ω模型的最大優勢.因此采用雷諾平均Navier-Stokes (RANS)方程結合SSTk-ω湍流模型模擬俘能結構的渦激振動響應過程.

無量綱不可壓縮RANS方程為

(12)

(13)

(14)

(15)

μt為湍流黏度[21],

(16)

2.2 俘能結構

波浪錐型俘能結構如圖2所示,其中DH和DL分別為波浪錐柱對應斜率為k的直錐柱最大直徑和最小直徑,結構的幾何特征表達式如下:

圖2 波浪錐型結構示意圖Fig.2 Schematic diagram of wave conical structure

Dz=Dm+2acos(2πz/λ)+k(z-L/2)

(17)

式中:Dz表示波浪錐柱對應高度Z處的直徑;L=7Dm;a、λ分別表示波浪錐柱的波幅與波長.根據楊耀宗[22]的研究,選取振動性能最佳的一組參數,即k=0.05,a/Dm=0.1,λ/Dm=1.75.

2.3 計算模型

計算域大小為30Dm×20Dm×10Dm,笛卡爾坐標系原點位于圓柱中心,圓柱中心距入口10Dm, 距出口20Dm,x正方向為順流向,y方向為渦激振動的方向,z正方向沿俘能結構中心向上,如圖3所示.設置來流為空氣,入口邊界采用均勻速度入口(Velocity inlet),在雷諾數Re=3 900的條件下,入口速度為5.772 m/s;出口邊界為壓力出口(Pressure outlet),相對壓力為0 MPa;上下面及側面設置為對稱邊界(Symmetry);俘能結構表面設置為無滑移壁面(No slip wall).

2.4 可靠性驗證

2.4.1計算策略驗證 采用瞬態CFD方法,為保證計算策略的準確性,需要對網格劃分數量N和時間步長Δt進行無關性驗證. 如圖4所示,驗證對比折合流速范圍為Ur=4.8~13.2,質量比與阻尼比固定為m*=10,ζ=0.01下的波浪錐形俘能結構的振幅比.圖4(a)對3種網格密度(904 337/1 382 078/1 499 493)進行了網格敏感性研究,圖4(b)則對3個時間步長(0.001/0.000 5/0.000 25)進行了時間步長獨立性驗證.從圖中可以看出,網格數量達到 1 382 078 后振幅幾乎不再變化,而時間步長設定Δt=0.000 5 時柱體振幅已和Δt=0.000 25 時近似相同.進一步地,以升力系數、阻力系數和振幅作為對比因素,選取Ur=7.2時波浪錐型俘能結構進行渦激振動數值模擬的驗證結果,如表1所示.可以認定Case4網格密度N=1 380 278、無量綱時間步長Δt=0.000 5 的計算精度與速度正確合理,平均阻力系數Cd,m、脈動升力系數Cl,r均已近似不再變化,最大振幅Amax與文獻[17]中相差不到0.3%,可用于波浪錐型俘能結構仿真實驗.

表1 無關性驗證結果對比Tab.1 Comparison of verification results

圖4 計算策略無關性驗證Fig.4 Calculation strategy independence verification

2.4.2實驗驗證 渦激振動實驗布置如圖5(a)所示,試驗段4個側面采用可拆卸亞克力擋板,實驗模型上下兩端開孔方便炭纖維桿支撐固定并連接到尼龍線.波浪錐型俘能結構如圖5(b)所示,采用8200Pro材料3D打印制作.位移測量原理如圖 5(c)所示, 圖中U∞為風洞來流. 風洞來流繞過實驗模型在其后方產生渦脫,渦脫頻率接近俘能結構-尼龍線-彈簧系統固有頻率時發生共振,通過試驗臺下方激光位移器發射和回收激光,由采集設備輸出渦激振動位移時程曲線.整體測位移裝置包括示波器、電源、控制器和激光頭.

圖5 渦激振動實驗裝置Fig.5 Vortex-induced vibration experimental device

振幅比驗證對比如圖6所示.圖中:Num表示數值仿真;Exp表示實驗測試.由圖可見,質量比m*=10,阻尼比ζ=0.01時仿真結果與文獻[22]中結果非常吻合,鎖頻區間和振幅比基本一致.因為風洞實驗箱和俘能結構模型尺寸并不完全匹配,所以采用等比例設計放大俘能結構,其質量比達到m*=46.2,在與其相同的質量比、阻尼比參數下進行仿真對比,發現鎖頻區間與振幅比近似相同,曲線變化趨勢較接近,因此認定本文的計算策略和用戶自定義函數(UDF)程序正確可靠.

圖6 振幅比驗證對比Fig.6 Comparison of amplitude ratio verification

3 結果分析

3.1 方案設計

由上文可知,基于渦激振動的小型風能采集器可以在相對較低的風速下產生運動,可以從各種地區頻繁和廣泛的風速中俘獲能量,如Vortex無葉片風力發電機在風速區間為1~9 m/s的地區進行模型試驗.本文設計的波浪錐型俘能結構希望應用于城市環境下的小型區域分布式無葉片風力發電機發電系統,而城市環境風速一般在1~7 m/s范圍內[23],因此確定Re=3 900 也即風速v=5.772 m/s條件下進行仿真模擬,從而切合城市環境常見風速條件.因此,探究合適的波浪錐型俘能結構質量比、阻尼比組合,有效地將低速環境風的動能轉換為電能至關重要.表2為質量比、阻尼比組合的實驗方案設計,探討質量阻尼參數對波浪錐型俘能結構渦激振動響應及能量轉換的影響規律.

表2 實驗方案設計Tab.2 Experimental design

3.2 m*和ζ對振幅比、頻率比的影響

對振幅A、振動頻率fosc進行無量綱化處理:A*=A/Dm,f*=fosc/fn.其中fn為結構固有頻率,A*和f*為振幅比、頻率比.A*和f*直接影響波浪錐型俘能結構捕獲能量的大小.

圖7(a)為波浪錐型俘能結構振動響應振幅比A*隨折合流速Ur的變化曲線.可以看出質量比和阻尼比對波浪錐型俘能結構鎖頻區間有顯著影響.當阻尼比ζ=0.05不變的情況下,鎖頻折合流速區間會隨著質量比的增加而縮短.阻尼比對鎖頻區間也有著相似的影響,由m*=5的4組方案可以看出,隨著阻尼比的增加鎖頻區間減小.m*=5,ζ=0.01的鎖頻區間最廣,在6.6

圖7 無量綱振幅比、頻率比曲線Fig.7 Dimensionless amplitude ratio and frequency ratio curves

圖7(b)給出了頻率比f*隨折合流速的變化曲線.Ur<6.0時,頻率比趨于0,此時俘能結構未進入鎖頻區間,當m*=2,ζ=0.05,Ur=6.0時,頻率比變化到0.821,隨后頻率比隨著折合流速增加而逐漸增加,當Ur=10.8時,f*增加到1.029,波浪錐型俘結構的振動頻率大于固有頻率,但此時振動幅值已經較小即將退出鎖頻區間.其他質量比、阻尼比條件下頻率比f*的變化規律類似,進入鎖頻區間時,頻率比迅速從0上升至趨于1附近,隨后在整個鎖頻區間,頻率比f*都隨折合流速的增加而增加并在退出鎖頻區間前達到最大,隨后退出鎖頻區間f*降到0附近.

可以看出質量比m*對頻率影響比f*的影響大,以ζ=0.05為例,同一折合流速下,較大的質量 比對應的頻率比也越大.而阻尼比ζ對頻率比的影響較小,當m*=5不變而ζ改變時,進入鎖頻區間后各折合流速下的f*基本重合.

3.3 m*和ζ對升力系數幅值、相位角的影響

過往的研究中發現,除了振幅和頻率之外,發生渦激振動時的升力系數幅值Cl和相位角正弦值sinφ同樣會影響流體和俘能結構之間的能量傳遞.圖8(a)、8(b)分別為升力系數幅值Cl、相位角正弦值sinφ隨折合流速的變化曲線.

圖8 升力系數幅值、相位角變化曲線Fig.8 Lift coefficient amplitude and phase angle change curves

從圖8(a)可以看出,當Ur=6.0時,m*=2,ζ=0.05的俘能結構升力系數幅值Cl迅速從0增加到0.537,隨后隨著折合流速的增加而逐漸減小,當Ur=10.8時,Cl降低到0.08,直至退出鎖頻區間減小到0.其他質量比、阻尼比條件下也有類似規律,升力系數幅值在進入鎖頻區間之后迅速達到最大值,隨后隨著折合流速的增加而逐漸減小,直至退出鎖頻區間趨于0.同時可以發現質量比、阻尼比對Cl也有較大影響,質量比、阻尼比越大,同一折合流速下對應的升力系數幅值越小.同一質量比、阻尼比條件下升力系數幅值對應的折合流速與振動幅值對應折合流速并不重合,說明升力越大并不意味著升力對渦激振動的作用越大.

學者們認為相位角φ在0°~180°范圍內時流體可將能量傳遞給俘能結構,結構發生周期性運動;而當相位角在 -180°~0°范圍內時,能量從俘能結構傳輸給流體,此時結構不會發生振動.從圖8(b)中可以看出,在鎖頻區間內sinφ均為正值,相位角在0°~180°之內,此時波浪錐型俘能結構發生渦激振動,俘能結構可以從流體中俘獲能量.未進入鎖頻區間前,相位角均在0°附近,剛進入鎖頻區間sinφ較小,之后sinφ隨著折合流速的增加而逐漸增加,在退出鎖頻區間之前達到最大值,直至退出鎖頻區間后,sinφ又降到0附近.質量比和阻尼比對sinφ也存在著影響,隨著質量比或者阻尼比的增加,在同一折合流速下,對應的sinφ先增加后減小,這也說明質量比和阻尼比對能量轉換的影響并不是單純的越大越好或者越小越好,存在一組合適的質量比、阻尼比組合使振動響應更好.

Cl和sinφ可以表征升力的大小和升力對振動響應的作用,升力系數幅值和相位角正弦值的乘積Clsinφ越大,表示波浪錐型俘能結構與來流的相互作用越強,如圖9所示.從圖中可以看出在進入鎖頻區間后,Clsinφ迅速達到最大值,隨后隨著折合流速增加而逐漸減小.當質量比不變時,阻尼比越大,俘能結構與流體間的相互作用也越大,在阻尼比不變的情況下,改變質量比也有類似的規律.質量比、阻尼比對能量轉換的影響是非線性的,Clsinφ隨著質量比、阻尼比增加而變大,但同時振幅與鎖頻區間也會減小.

圖9 Clsin φ隨折合流速變化曲線Fig.9 Clsin φ versus reduced velocity

為了更好地展示相位角對振動的影響,以m*=2,ζ=0.05的升力系數幅值與振動位移時程圖為例,如圖10所示.折合流速Ur=6.0,7.2時,波浪錐型俘能結構的升力系數和位移呈現“同相位”特征,此時φ值較小;當Ur=9.6時,升力系數與位移除呈現出“同相位”特征之外,升力系數還呈現出“拍”的特征;當Ur=10.8時,升力系數與位移交替出現“同相位”與“反相位”特征.結合圖8,Ur=6.0,7.2時sinφ較小,位移與升力系數相位差較小,表現為時程圖“同相位”,而當Ur=9.6,10.8時相位差較大,時程圖出現“拍”的特征和正反相位交替出現的現象.

圖10 m*=2,ζ=0.05的升力系數幅值與振動位移時程圖Fig.10 Time chart of Cl and vibration displacement at m*=2 and ζ=0.05

3.4 質量阻尼參數對能量轉換效率的影響

根據式(11)推導可以計算出各折合流速下的能量轉換效率η,不同質量比m*、阻尼比ζ、質量阻尼比m*ζ下能量轉換效率變化曲線如圖11所示.圖11(a)中同一折合流速下的能量轉換效率隨阻尼比的增加先增大后減小,最大能量轉換效率對應的折合流速并不與最大振幅相同.以ζ=0.05,m*=2,3,4為例,m*=2條件下Clsinφ未達到最大值,但此時振幅較大,能量轉換效率最高.說明能量轉換是各項參數綜合作用的結果.同時可以發現質量阻尼比m*ζ相同的情況下,其能量轉換效率曲線也十分接近.隨著阻尼比的增加,鎖頻區間內無量綱振幅值和能量轉換效率的增加趨勢變緩.可以看出Ur=7.8時m*=2,ζ=0.05的最大能量轉換效率略低于Ur=8.4時m*=4,ζ=0.05的轉換效率,但m*=2,ζ=0.05條件下鎖頻區間更廣,在整個折合流速范圍內俘能效率最高,當Ur=7.8時風能轉換效率達到2.8%.因此認為m*=2,ζ=0.05時,波浪錐形俘能結構的能量轉換效率最高.

圖11 能量轉換效率曲線Fig.11 Energy conversion efficiency curves

本次實驗使用的質量阻尼比m*ζ=0.05,0.10,0.15,0.20,計算出每個m*ζ下的最大能量轉換效率ηmax及平均能量轉換效率ηmea,整理出如圖11(b)所示的曲線.可以看出質量阻尼比m*ζ影響著能量轉換效率,而且存在著最優的m*ζ組合使得能量轉換效率最高.m*ζ=0.10的情況下能量轉換效率最高,平均能量轉換效率較m*ζ=0.05提高了近31%.當應用渦激振動發電時,引入發電機會造成阻尼比較大的結果,這也提供一種思路:如果想要提高能量轉換效率,可以在m*ζ一定的情況下,將俘能結構輕量化的同時提高發電阻尼在系統總阻尼中的比例.

3.5 fn對能量轉換效率的影響規律

以質量比、阻尼比組合m*=2,ζ=0.05為對象,研究固有頻率fn對波浪錐型俘能結構能量轉換的影響,圖12給出了固有頻率fn=5,10,20 Hz下振幅比隨折合流速、風速的變化曲線.從圖12(a)中可以看出,改變固有頻率不影響俘能結構發生渦激振動的鎖頻區間,同一折合流速下俘能結構的振幅隨固有頻率的增加而變大.

圖12 固有頻率對振幅比的影響Fig.12 Influence of natural frequency on amplitude ratio

本文的波浪錐型俘能結構能量轉換效率研究限定在Re=3 900情況下,風速v固定為5.772 m/s.而在實際應用中來流速度會不斷變化,因此有必要研究實際風速下的能量轉換規律.根據Ur=v/(fnDm)可以得出不同固有頻率下折合流速Ur對應的實際風速,圖12(b)為振幅比隨風速的變化曲線.從圖中可以看出,固有頻率增加不但可以提高振幅,同時又可以拓寬渦激振動鎖頻區間,波浪錐型風力俘能結構可以在更大的風速范圍內俘獲能量.

固有頻率對能量轉換效率的影響如圖13所示.從圖13(a)中可以看出,能量轉換效率η在進入鎖頻區間之后迅速達到峰值,隨后η隨著Ur的增加逐漸降低,在Ur不變的條件下,固有頻率越大能量轉換效率越高.在質量比阻尼比組合為m*=2,ζ=0.05,固有頻率fn=20 Hz條件下,波浪錐型俘能結構的風能轉換效率達到5.1%,并且風能轉換效率會隨著固有頻率的提高而繼續增大.從圖13(b)中可以看出,fn的增加也會使得波浪錐型俘能結構俘能的實際風速范圍增加.綜上所述,在運用渦激振動進行風力發電時,適當的提高波浪錐型風力俘能結構固有頻率可以提高能量轉換效率,拓寬俘能風速范圍,有效提升其能量轉換特性.

圖13 固有頻率對能量轉換效率的影響Fig.13 Influence of natural frequency on energy conversion efficiency

4 結論

(1) 質量比m*和阻尼比ζ對振動響應有較大影響,m*和ζ越大,渦激振動振動幅值和鎖頻區間越小.質量比相較于阻尼比對頻率比影響較大.

(2) 質量阻尼比m*ζ相同的情況下,不同參數組合其能量轉換效率曲線變化趨勢基本一致.m*=2,ζ=0.05的波浪錐型俘能結構俘能的鎖頻區間更廣,能量轉換效率更高.

(3) 折合流速不變的情況下,固有頻率fn越大,能量轉換效率越高,適當的提高波浪錐型風力俘能結構的固有頻率可以有效提升其能量轉換特性.

猜你喜歡
振動效率結構
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
提升朗讀教學效率的幾點思考
甘肅教育(2020年14期)2020-09-11 07:57:42
振動與頻率
天天愛科學(2020年6期)2020-09-10 07:22:44
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
中立型Emden-Fowler微分方程的振動性
論《日出》的結構
跟蹤導練(一)2
創新治理結構促進中小企業持續成長
現代企業(2015年9期)2015-02-28 18:56:50
“錢”、“事”脫節效率低
中國衛生(2014年11期)2014-11-12 13:11:32
主站蜘蛛池模板: 在线无码av一区二区三区| 97久久人人超碰国产精品| 国产成人福利在线| 欧美成人区| 国产精品一线天| 欧美激情,国产精品| 亚洲Aⅴ无码专区在线观看q| 国产极品美女在线播放| 波多野结衣的av一区二区三区| 久草青青在线视频| 久久黄色影院| 欧美在线精品一区二区三区| 噜噜噜久久| 久久亚洲国产最新网站| 国产精品网曝门免费视频| 国产性猛交XXXX免费看| 精品视频在线一区| 国产青榴视频在线观看网站| 青青青国产视频| 国产99视频精品免费视频7| 国产毛片片精品天天看视频| 日韩一区二区在线电影| 9久久伊人精品综合| 久久狠狠色噜噜狠狠狠狠97视色| 怡春院欧美一区二区三区免费| 久久久久亚洲AV成人人电影软件 | 免费视频在线2021入口| 国产精品无码一区二区桃花视频| 国产在线高清一级毛片| 中文字幕在线欧美| 久久亚洲黄色视频| 综合天天色| 久久99热这里只有精品免费看| 尤物精品视频一区二区三区| 日本精品αv中文字幕| 999精品色在线观看| 成人小视频在线观看免费| 日韩经典精品无码一区二区| 亚洲无码日韩一区| 九九九国产| 国产成本人片免费a∨短片| 久久情精品国产品免费| 国产成人综合亚洲欧美在| 亚洲成人精品| 无码国产偷倩在线播放老年人| 免费a级毛片视频| 久久99国产综合精品1| 欧美成人精品一级在线观看| 丰满少妇αⅴ无码区| a级毛片在线免费| 无码高潮喷水专区久久| 日韩毛片视频| 亚洲无码熟妇人妻AV在线| 国产成人高清精品免费软件| 亚洲小视频网站| 色噜噜久久| 国产成人a在线观看视频| 中文字幕一区二区人妻电影| 国产精品任我爽爆在线播放6080| 666精品国产精品亚洲| 亚洲 成人国产| 在线精品亚洲国产| 2022国产91精品久久久久久| 一区二区日韩国产精久久| 久久综合五月| 免费福利视频网站| 久久国产av麻豆| 在线观看亚洲人成网站| 4虎影视国产在线观看精品| 五月激情综合网| 欧美啪啪网| 久久精品人人做人人爽| 亚洲人成电影在线播放| 久久国产成人精品国产成人亚洲 | 色偷偷一区二区三区| 911亚洲精品| 欧美在线中文字幕| 国产亚洲欧美日韩在线一区二区三区| 精品中文字幕一区在线| 欧美日韩激情在线| 国产精品尤物在线| 美女啪啪无遮挡|