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

漂浮式風力機結構動力學響應TMD控制及其參數優化研究

2018-12-21 12:11:12丁勤衛郝文星葉柯華王淵博
振動與沖擊 2018年23期
關鍵詞:優化結構

丁勤衛, 郝文星, 李 春, 葉柯華, 王淵博

(上海理工大學 能源與動力工程學院,上海 200093)

隨著陸上風電場可開發資源的減少,海上風能因能量密度高和湍流度低等優勢逐漸為世界各國重視,“由陸地向海洋”業已成為未來風電場發展的必然趨勢[1-2]。目前,海上風電場主要集中在淺水區域,對于風資源更豐富、風況更優的深水區域,從經濟性考慮必須采用漂浮式。海上風力機運行過程中需始終承受波浪載荷作用,這也是與陸上風力機的最大不同,因此,其載荷特性更加復雜。較之于固定式風力機,漂浮式風力機因底部基礎不固定而導致的特有搖蕩特性使得其始終處于受力不平衡、運動非定常狀態,似此非線性載荷不僅影響機艙內傳動系統的正常工作,亦會降低風機的發電效率,甚至可能導致塔架的屈曲、傾覆等事故[3-5]。此外,與傳統海工平臺相比,漂浮式風力機平臺重心位置更高、水線面更小,波浪激勵更加劇了平臺的運動甚至傾覆。因此,研究漂浮式風力機在復雜多變的在役海洋環境中的安全穩定運行具有理論和應用價值。

目前,國內外學者針對漂浮式風力機穩定性問題展開了一些研究。文獻[6]提出將螺旋側板應用于漂浮式風力機Spar平臺的結構設計中,并采用勢-黏結合方法研究螺旋側板對Spar平臺動態響應影響,發現螺旋側板可明顯增大其穩定性。文獻[7]針對漂浮式風力機對風浪激勵響應過大問題提出了共用系泊的大型漂浮式風電場概念設計,并基于AQWA研究其運動響應特性,發現共用系泊可明顯增大平臺穩定性。上述研究不足之處在于將葉片、塔架及平臺等結構簡化為剛體,風載荷簡化為軸向推力,顯然這種簡化無法精確描述非定常氣動載荷,更無法揭示風浪激勵下塔架、葉片等結構的變形、屈曲等非線性動力響應過程。

部分學者通過控制葉片變槳和電機變扭矩以降低葉片氣動載荷進而達到提高漂浮式風力機穩定的目的。例如,Namik等[8]提出針對漂浮式風力機的葉片獨立變槳技術,研究獨立變槳對平臺運動的控制效果,結果表明獨立變槳降低平臺運動的幅度可達39%。再如,Fischer等[9]提出一種基于加速度反饋的非線性控制方法,在理論上分析漂浮式風力機的穩定性。上述方法雖可在一定程度上增大漂浮式風力機穩定性,但會加劇葉片根部疲勞載荷。

此外,一些學者提出將主動結構控制技術應用于提高漂浮式風力機穩定性。例如,文獻[10]建立了風波聯合作用下Spar平臺漂浮式風力機AMD控制系統運動方程,并利用最優二次控制算法實現對風力機縱蕩運動的主動控制。主動控制技術雖可取得比被動控制更好的效果,但主動控制不易實現且成本高昂。此外,現代風力機作為人類目前建造的最大的旋轉機械其結構巨大,主動控制時所需控制力非常巨大,已有學者指出,采用被動控制技術最適合于漂浮式風力機結構控制[11]。

調頻減振裝置主要包括:調頻質量阻尼器(Tuned Mass Damper,TMD)和調頻液體阻尼器(Tuned Liquid Damper,TLD),其作為結構控制技術被廣泛應用于高聳柔性結構設計中。實驗及數值計算均表明其對高樓、電視塔等高聳柔性結構減振效果明顯[12]。風力機塔架屬典型高聳柔性結構,因此,部分學者提出將調諧減振裝置應用于風力機被動結構控制領域。國內外學者針對TMD/TLD對提高風力機穩定性效果亦開展了相關研究。文獻[13]將TMD應用于海上風力機塔架中,并采用有限元方法研究其對風力機塔架振動的減振效果,但其將風輪和機艙簡化為塔尖堆聚非旋轉集中質量,此簡化無法獲知葉片等結構的動力學響應,更無法描述葉片-塔架等結構的耦合振動效應。Gordon等[14-16]研究TMD對漂浮式風力機穩定性的影響,但其不僅未考慮外載荷作用,且其工作僅為自由振動分析。文獻[17]借助開源軟件FAST考慮風浪載荷影響研究了TMD對Spar平臺漂浮式風力機穩定性影響,其工作雖然考慮水動-氣動-結構耦合作用,但其并未對TMD參數進行優化分析。

綜上所述,旨在提高漂浮式風力機穩定性的研究雖在漂浮式平臺結構設計、系泊系統優化及TMD結構控制方面開展了一些工作,但是:①多數研究均做了適當簡化,漂浮式風力機系統是復雜的強耦合非線性系統,如葉片揮舞-塔架前后振動耦合,如葉片擺振-塔架左右振動耦合,對結構、自由度等的簡化無法描述氣動-水動-結構耦合作用;②顯然,TMD控制效果與其固有質量、彈簧剛度及阻尼器阻尼均有關,但上述研究并未涉及其結構參數優化,因而導致其對振動位移及載荷等指標抑制率并不明顯;③目前國內外學者針對提高漂浮式風力機穩定的研究大都基于張力腿平臺、Spar平臺等,基于Barge平臺漂浮式風力機的甚少。

因此,本文借鑒相關學者研究經驗,提出將TMD結構控制技術應用到提高Barge平臺漂浮式風力機穩定性的研究中,進一步采用多島遺傳算法(Multi-Island Genetic Algorithm,MIGA)對所配置TMD的參數進行全局優化,以期為海上漂浮式風力機的穩定性控制提供一定的理論參考。

1 研究對象

本文研究對象為基于ITI Barge平臺NREL 5 MW風力機,風力機參數[18]如表1所示。平臺參數[19]如表2所示。漂浮式風力機模型,如圖1所示。

2 研究方法及TMD控制

基于NREL開發的氣動彈性結構動力學仿真軟件FAST研究TMD對漂浮式風力機的結構動力學響應影響,相關學者采用FAST對風力機結構動力學響應進行仿真分析并將計算結果與GH-Bladed進行對比發現二者計算結果具有極高的吻合度,可驗證FAST計算結果準確可信[20]。

表1 NREL 5 MW風力機參數

表2 ITI Barge平臺參數

FAST考慮氣動-水動-結構-控制耦合基于時間推進方法實現對非線性運動方程的求解[21-22]。氣動模塊基于Pitt-Peters動態入流理論并考慮軸向及切向風誘導的影響求解風輪平臺誘導速度,進一步基于葉素動量理論考慮葉尖損失、輪轂損失及Beddoes-Leishman動態失速模型修正求解風輪氣動力;水動模塊計算漂浮式平臺的水動力載荷,包括水線面面積矩和浮力;針對大尺度結構采用輻射/繞射理論求解波浪力,針對小尺度結構采用Morison方程考慮黏性效應和附加質量效應影響,但無法考慮VIV(Vortex Induced Vibration)激振力;系泊系統采用懸鏈線模型。彈性模塊通過Kane方法結合模態法并考慮子結構剛柔特性建立風力機多體動力學模型,其中葉片、塔架及低速傳動軸為柔性體模型,具有分布的質量、剛度及模態,而將輪轂、高速傳動軸及機艙等視為剛性體模型。以氣動模塊求解氣動力及水動模塊求解的水動力作為輸入激勵,得到該時間步的結構動力學響應并反饋至控制模型、水動模塊及氣動模塊,控制模塊則根據結構模塊反饋信息作出相應的控制策略,如變槳控制、偏航控制及TMD控制等。具體仿真流程,如圖2所示。

圖1 Barge平臺漂浮式風力機示意圖

圖2 仿真流程圖

TMD由質量系統、彈簧系統、阻尼系統及支撐系統組成,其通過改變自身阻尼或剛度來達到調整自振頻率的目的,使其接近被減振結構的固有頻率或外部激勵頻率,當結構在外部激勵下產生振動時,將帶動TMD振動,TMD產生的調諧慣性力將反作用到結構上,并通過阻尼系統將能量耗散,以達到結構穩定性控制的目的。本文將TMD布置在機艙內部示意圖,如圖3所示。

圖3 機艙布置TMD示意圖

3 湍流風場及不規則波建模

3.1 風場模型建立

為更真實模擬時域高風速湍流風作用下漂浮式風力機的結構動力學響應特性,建立隨時間及空間均變化的湍流風是仿真首先需要解決的問題。目前,常見的風場建模方法有:①基于測風塔實測數據,通過模糊邏輯預測、神經網絡預測等方法得到一定空間范圍內風速分布,該方法針對小空間風場建模較為準確,若要針對大空間風場建模,需要極度豐富實測數據,成本較高;②基于LES方法,考慮大氣邊界層效應等實現風場建模,但該方法所需計算資源極大;③基于風電場實測數據,通過氣象分析方法獲得風場風速分開,但該方法同樣適合小尺度空間風場建模,針對大尺度空間風場建模誤差較大。鑒于此,本文基于Kaimal風譜并考慮空間相干性建立三維時變風場模型。

根據所研究風力機幾何參數,設定風場覆蓋區域為195 m(水平方向)×195 m(垂直方向),之所以風場區域遠遠超過風輪掃略區域,是因為本文研究對象為漂浮式風力機,漂浮式平臺六自由度運動導致葉片存在垂蕩、橫蕩運動,為保障所建三維風場始終完全覆蓋葉片,故將風場布置較大。進一步設定網格節點為15×15,如圖4所示。考慮空間相干性,通過Kaimal風譜獲得每一節點的風速分布,進一步通過空間相干模型獲得整個風場的風速分布。空間相干模型如下

(1)

式中:Si,j(f)為節點i,j的互功率譜;C(Δr,f)為空間相干大小,節點之間距離為Δr;Si,i(f)Sj,j(f)分別為節點i,j的功率譜。

以輪轂中心為參考點,以時歷平均風速11.4 m/s為參考風速,仿真時間為600 s,建立三維時變湍流風風場,輪轂高度處三維風速分布,如圖5所示。輪轂點時域風速分布,如圖6所示。由圖6可知,所建湍流風在u方向波動范圍為10.52~12.55 m/s,v方向波動范圍為-1.052~0.993 2 m/s,w方向波動范圍為-0.637 8~0.664 8 m/s。其中u方向為來流入射方向,v方向水平方向,w方向為垂直方向。

圖4 風場計算區域及網格劃分

圖5 輪轂高度處風速分布

圖6 輪轂點時域風速分布

3.2 波浪譜及不規則波浪建模

波浪譜以風因素和波浪因素為參量,通過定義有義波高、波浪周期及有限風區等參數即可得到波浪的大致形式。波浪譜是隨機波浪的重要統計信息,可直接給出波能相對頻率和方向的分布。國內外學者所做海洋結構的研究針對波浪載荷的考慮幾乎全部采用波浪譜方法[23-27]。目前,常見的波浪譜有P-M譜、Jonswap譜和布氏譜等[28]。其中P-M譜根據大西洋波浪統計數據分析得出,其適合描述充分發展的波浪[29],其形式為

(2)

式中:U為某參考高度處平均風速;g為重力加速度;ω為波浪圓頻率。

選取P-M波浪譜生成有義波高為5 m,譜峰周期為12.4 s的波浪,基于設定參數建立不規則波浪,如圖7所示。

4 TMD對漂浮式風力機穩定性影響

4.1 TMD對平臺搖蕩特性影響

TMD結構控制技術廣泛應用于高樓、電視塔等高聳結構,結構參數大都按照Den Hartog原則[30-31]確定,鑒于目前尚未有漂浮式風力機TMD結構控制參數設定的原則,故本文借鑒高樓、電視塔等設計經驗,采用Den Hartog原則設定。此處TMD質量、阻尼及剛度質量為20 000 kg,剛度為10 000 N/m,阻尼為50 000 N/(m/s)。漂浮式風力機在TMD控制下平臺搖蕩特性的時程圖,如圖8所示。其中橫坐標為時間,縱坐標為平臺各自由度運動響應。

圖7 不規則波浪模型

由圖8可知,TMD控制下,TMD控制前后,平臺縱蕩、垂蕩及縱搖響應改變甚微,幾乎未發生變化,但平臺橫蕩、橫搖及首搖響應變化較明顯。橫蕩波動幅值由-1.03~1.12 m降低到-0.65~0.85 m,波動幅度降低約30%,橫蕩最大值由1.12 m減小到0.85 m。橫搖波動幅值由-1.11~1.38°減小到-0.54~0.88°,波動幅度降低約40%,最大值由1.37°減小到0.88°。配置TMD前后,首搖運動雖發生較大變化,但運動幅度變化較小。由計算可知,TMD控制前后,平臺橫蕩標準差分別為0.477和0.352,平臺橫蕩穩定性提高約25%;平臺橫搖標準差分別為0.492和0.289,比較可知,平臺橫搖穩定性提高約40%。

由上述時域分析可知,平臺縱蕩、垂蕩、縱搖和首搖的變化不明顯。因此,此處僅給出橫蕩與橫搖的幅值譜,如圖9所示。其中橫坐標為頻率,縱坐標為幅值響應。

由圖9可知,平臺橫搖的頻譜峰值約為0.09 Hz,橫蕩峰值頻率約為0.01 Hz。此外,平臺的橫蕩和橫搖頻譜峰值減小非常明顯,進一步從頻域角度驗證了TMD對平臺橫蕩和橫搖控制效果顯著。

圖9 平臺橫蕩及橫搖運動幅頻特性曲線

4.2 TMD對塔架頂部位移影響

漂浮式風力機在TMD控制下的塔尖位移時程圖,如圖10所示。其中橫坐標為時間,縱坐標為塔尖前后及側向位移。

圖10 塔尖前后及側向位移時歷曲線

由圖10可知,TMD控制前后,塔尖前后位移均在-0.8~1 m之間無規律波動,TMD對塔尖前后位移影響甚微。配置TMD前后,前300 s差別不明顯,300 s之后,TMD控制后塔尖側向位移明顯降低,且減小幅度非常明顯。計算可知,TMD控制前后,塔尖側向位移的標準差分別為0.066和0.041,塔尖側向位移穩定性提高約38%。

對塔尖前后及側向位移時域數據進行快速傅里葉變化,得到其幅值譜,如圖11所示。其中橫坐標為頻率,縱坐標為幅值響應。

圖11 塔架前后及側向位移幅頻特性曲線

由圖11可知,塔尖前后及側向位移的頻譜峰值均約為0.09 Hz。比較有無TMD控制兩種情況可知,無TMD控制時,塔尖前后及側向位移幅值均大于受控時的幅值,但是塔尖前后位移幅值變化較小,而側向位移幅值減小非常明顯。

5 TMD參數優化

由前文時域及頻域結果分析可知,機艙布置TMD對平臺橫蕩、橫搖運動及塔尖側向運動起到明顯的抑制作用,但影響TMD控制效果的參數包括質量系統質量、彈簧系統剛度及阻尼系統阻尼。TMD結構參數組合對優化漂浮式風力機穩定性效果必然不同,因此,探究較優的TMD結構參數組合對于提高漂浮式風力機的穩定性具有重要意義。本文采用多島遺傳優化算法對TMD結構參數進行優化,以確定較優的TMD結構參數組合。

多島遺傳算法[32]是在傳統遺傳算法基礎上建立的一種基于群體分組的并行性遺傳算法。多島遺傳算法將整個種群分解為多個子群(“島”),并將各個子群互相隔絕于不同的“島嶼”上,對每個子群中的個體進行傳統遺傳算法操作(選擇、雜交、變異操作),各個子群獨立地進化,而不是全部種群采用相同的進化機制,并且各個“島嶼”間以一定的時間間隔進行“遷移”操作,使各個“島嶼”間進行信息交換。并且每隔幾代挑選子群中的個體進行交換(遷移操作),保證了進化過程中優化解的多樣性,從而有效抑制了早熟現象的出現,有利于找到全局最優解。

MIGA反復地使用算子和選擇原則,從親代到子代再到孫代直至重孫代不停地繁衍,從而種群對環境的適應性得到不斷地升高。流程如下:

第1步 初始化群體;

第2步 計算個體的適應度函數值;

第3步 按個體適應度值決定的某種規則選擇進入下一代的個體;

第4步 按概率Pc進行交叉操作;

第5步 按概率Pm進行突變操作;

第6步 若未滿足停止條件,則轉“第2步”,否則進入“第7步”;

第7步 輸出種群中適應度值最優的染色體作為問題的滿意解或最優解。

MIGA流程圖,如圖12所示。

目標函數:以漂浮式風力機塔架的塔尖側向位移標準差σ1及平臺的橫搖標準差σ2之和最小為目標函數。

約束條件:

質量m:10 000 kg≤m≤40 000 kg

剛度k:5 000 N/m≤k≤25 000 N/m

阻尼d:6 000 N/(m/s)≤d≤21 000 N/(m/s)

優化設計所使用的MIGA算法控制參數中種群數、交叉概率、變異概率、遷移概率及最大代數分別為50、0.85、0.02、0.3和20。

6 TMD參數優化結果與分析

依據所設計的優化目標及多島遺傳算法參數設置方法,采用MIGA算法對計算結果進行優化,對質量、剛度及阻尼進行搜索尋優。

圖12 多島遺傳算法優化流程圖

圖13為塔尖側向位移標準差σ1值隨質量m、剛度k及阻尼d樣本點分布的四維云圖。圖14為σ1隨質量m變化的剖面圖。由圖14可知,隨著TMD質量增加,σ1變化較明顯,呈現出先減小后增大的趨勢,在15 000~25 000 kg區間內,塔尖側向位移偏差σ1較小。進一步,將圖14質量為15 000~25 000 kg間的剖面圖局部細化,得到圖15所示部分剖面圖。由圖15可知,當質量為18 695~22 935 kg區間內,塔尖的側向位移偏差存在極小值;當阻尼取值不變時,σ1隨剛度變化趨勢較小,說明剛度變化對σ1影響不明顯;當剛度取值不變,阻尼取值在12 000~17 000 N/(m/s)區間時,σ1存在明顯變化。

圖13 σ1四維圖

圖16所示平臺橫搖標準差σ2值隨質量m、剛度k及阻尼d樣本點分布的四維云圖。圖17所示σ2隨質量m變化的剖面圖。由圖17可知,隨著TMD質量增加,σ2變化較明顯,呈現出先減小后增大的趨勢,在20 000~30 000 kg區間內,平臺的橫搖偏差σ2較小。將該區域細化,得到局部剖面圖,如圖18所示。由圖18可知,當質量取值在18 315~25 905 kg區間內,平臺橫搖的標準差σ2存在極小值;當阻尼取值不變時,剛度變化對σ2影響較小;當剛度取值不變,阻尼取值在12 000~17 000 N/(m/s)區間時,σ2存在明顯變化。

圖14 σ1剖面圖

針對平臺橫搖運動,由圖19可知,TMD控制后,平臺橫搖運動得到明顯抑制,優化后TMD控制效果更好。優化TMD控制下,平臺的橫搖范圍減小到-0.49~0.85°。由橫搖運動時間序列數據計算可知,優化TMD控制后,平臺橫搖運動標準差約為0.23,均小于前文無TMD控制時的0.492和普通TMD控制時的0.289,平臺穩定性提高約53%。結果表明,參數優化后TMD對平臺橫搖運動具有顯著的控制效果。針對塔尖側向位移,由圖19可知,TMD控制后,塔尖側向位移得到明顯抑制,與對平臺橫搖運動的控制效果類似,優化后的TMD對塔尖側向位移的抑制效果更好。優化TMD控制下,塔尖側向位移波動范圍為-1.53~0.071 m。由塔尖側向位移時間序列計算可知,優化TMD控制后,塔尖側向位移標準差約為0.034,亦均小于前文無TMD控制時的0.066和普通TMD控制時的0.041,塔尖側向位移穩定性提高約50%。結果表明,參數優化后TMD對塔尖側向位移同樣具有顯著的控制效果。綜上所述,優化后TMD對平臺橫搖運動、塔尖側向振動具有明顯的抑制效果,計算結果與分析表明本文所提出的優化方法、優化算法及計算結果的準確性。

表3為MIGA算法優化后得到的TMD設計參數。為驗證優化方法的可靠性及優化結果的準確性,依據優化結果重新配置TMD質量、剛度和阻尼并進行計算,結果如圖19所示。為更清晰展示優化效果,此處選取300~600 s區間數據。

m=16 575 kgm=18 695 kgm=20 815 kgm=22 935 kgm=25 055 kgm=27 175 kg

圖15σ1隨m變化局部剖面圖

Fig.15 A partial cross-section ofσ1

圖16 σ2四維圖

圖17 σ2剖面圖

7 結 論

本文以NERL的ITI Barge型漂浮式風力機為研究對象,采用在機艙配置TMD控制方法,考慮漂浮式風力機實際部署海域的風況與海況,研究環境載荷作用下漂浮式風力機無控和受控時的穩定性,后采用MIGA對TMD各參數進行全局優化,結論如下:

(1) TMD對漂浮式風力機塔尖側向位移控制效果較明顯,穩定性提高約38%;TMD控制下,漂浮式風力機平臺橫蕩和橫搖降幅明顯。其中,平臺橫蕩穩定性提高了25%,橫搖穩定性提高了40%。

(2) 通過MIGA算法對TMD參數進行全局優化,結果表明漂浮式風力機的塔尖側向位移及平臺的橫搖隨TMD參數變化趨勢相似。隨TMD質量增加,均呈現出先減小后增大的趨勢;當阻尼取值不變時,剛度變化影響不明顯;當剛度取值不變,阻尼存在一定變化。

(3) 通過對比分析漂浮式風力機無控、普通TMD及優化后TMD的控制效果,驗證了MIGA優化結果的有效性,可為海上漂浮式風力機配置TMD提供一定的理論參考。

m=15 605 kgm=18 315 kgm=21 025 kgm=23 735 kgm=25 905 kgm=28 615 kg

圖18 σ2隨m變化局部剖面圖

圖19 優化后TMD對平臺橫搖及塔尖側向位移影響

猜你喜歡
優化結構
超限高層建筑結構設計與優化思考
房地產導刊(2022年5期)2022-06-01 06:20:14
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
民用建筑防煙排煙設計優化探討
關于優化消防安全告知承諾的一些思考
一道優化題的幾何解法
由“形”啟“數”優化運算——以2021年解析幾何高考題為例
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
新型平衡塊結構的應用
模具制造(2019年3期)2019-06-06 02:10:54
論《日出》的結構
創新治理結構促進中小企業持續成長
現代企業(2015年9期)2015-02-28 18:56:50
主站蜘蛛池模板: 99福利视频导航| 中国一级特黄大片在线观看| 天堂成人av| 国产成人精品亚洲77美色| 中文成人在线| 婷婷色狠狠干| 麻豆国产精品一二三在线观看| 一级毛片不卡片免费观看| 色偷偷av男人的天堂不卡| 天天躁日日躁狠狠躁中文字幕| 欧美精品在线看| 视频一区亚洲| 永久天堂网Av| 国产丝袜无码精品| 99视频在线观看免费| 97久久人人超碰国产精品| 精品国产毛片| 亚洲色大成网站www国产| 亚洲综合久久成人AV| 国产jizz| 欧美在线伊人| 久久天天躁狠狠躁夜夜躁| 18禁黄无遮挡免费动漫网站| 97精品伊人久久大香线蕉| 国产永久免费视频m3u8| 91毛片网| 欧美一级特黄aaaaaa在线看片| 免费国产不卡午夜福在线观看| 久久人人妻人人爽人人卡片av| 一区二区三区四区日韩| 无码内射在线| 国产精品香蕉在线| 国产女人在线观看| 999国内精品视频免费| 激情综合婷婷丁香五月尤物| 亚洲欧洲日韩国产综合在线二区| 97综合久久| 精品一区二区三区波多野结衣 | 一边摸一边做爽的视频17国产| 91尤物国产尤物福利在线| 老司机精品一区在线视频 | 欧美日韩中文字幕在线| 激情無極限的亚洲一区免费| 精品国产美女福到在线不卡f| 精品丝袜美腿国产一区| 国产激爽爽爽大片在线观看| 久久久久亚洲AV成人人电影软件| 亚洲一级毛片免费看| 伦伦影院精品一区| 欧美精品亚洲二区| 人人91人人澡人人妻人人爽| 无码日韩人妻精品久久蜜桃| 欧美亚洲香蕉| 亚洲免费人成影院| 在线观看无码av五月花| 91精品国产无线乱码在线| 日韩欧美亚洲国产成人综合| 国产成人资源| 国产精品999在线| 国产美女无遮挡免费视频网站| 日韩色图区| 亚洲成人黄色在线| 亚洲开心婷婷中文字幕| av在线手机播放| 久久久久夜色精品波多野结衣| 久久人妻系列无码一区| 欧美久久网| 成年av福利永久免费观看| 亚洲男女在线| 国产成人免费| 亚洲精品成人片在线观看| 久久精品亚洲中文字幕乱码| 青青草91视频| 浮力影院国产第一页| 91青青草视频在线观看的| 日韩精品一区二区三区swag| 国产综合无码一区二区色蜜蜜| 国产成+人+综合+亚洲欧美 | 午夜三级在线| 高清国产在线| 亚洲精品视频在线观看视频| 亚洲无卡视频|