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

基于遺傳算法的多自由度波浪能裝置浮體形狀優化

2022-07-05 03:41:22譚銘楊宇軒岑雨昊司玉林錢鵬張大海
中國艦船研究 2022年3期
關鍵詞:優化

譚銘,楊宇軒,岑雨昊,司玉林,錢鵬,張大海

浙江大學 海洋學院,浙江 舟山 316021

0 引 言

波浪能的高效利用需要波浪能發電裝置在各級能量轉換系統上都具有較高的效率。根據剛體在波浪中的振動理論,只有當波浪能發電系統與海浪之間產生共振時,海浪能才能最大程度地轉化為電能。然而,由于海浪的不規則性和非線性特性,不能簡單地將波浪能發電裝置(也稱波浪能裝置)視作質量彈簧阻尼系統,因此高效的電能轉換成為了一個難題。近年來,波浪能裝置的優化設計與波浪能控制已成為研究熱點。波浪能裝置優化是指對其各級能量轉換系統的傳動結構和元器件參數進行優化設計,例如Penalba 等[1]提出的波浪能裝置能量轉換方案的設計方法。波浪能控制是指在能量轉換系統已確定的情況下,根據實際海況中波浪條件的變化,通過引入先進控制系統,實時調整波浪能裝置某一環節的關鍵參數或者采取主動的驅動或約束動作,來改變波浪能裝置的運動響應模式和固有頻率,實現共振條件和最優電能轉換,例如Wang 等[2]提出的波浪能控制方法。

浮體作為一級能量捕獲機構,其效率決定了是否能夠高效驅動后端的動力輸出(power take off, PTO)系統。對于波浪能裝置而言,波浪能轉換至浮體機械能的過程涉及了浮體的水動力運動響應,而浮體形狀對浮體在波浪激勵下的水動力性能和能量轉換則有著很大影響。

目前,對于波浪能裝置浮體形狀的選取和優化一般基于枚舉法[3-7],但此方法很大程度上受人為因素的影響,所得浮體形狀并不是嚴格意義的最優解。近年來,迭代算法在波浪能裝置系統參數設計優化方面得到了廣泛應用。例如,Babarit等[8]結合梯度法與遺傳算法,以單位排水量下的吸收功率最大為目標對SEAREV 波浪能裝置的浮體形狀開展了多目標優化。McCabe 等[9-10]基于B 樣條曲面和頻域模型,應用遺傳算法計算了在垂蕩、縱蕩和縱搖方向同時運動的波浪能裝置的最優浮體形狀。Kurniawan 等[11]研究繞水下固定軸振蕩的均勻截面浮體,計算了給定頻率范圍內由主要參數定義的浮體最優形狀。Mohamed 等[12]采用進化類算法,選取切向力系數和透平效率為目標函數,基于FLUENT 計算流體力學軟件,對振蕩水柱式波浪能裝置的透平翼型優化進行了研究。Colby 等[13]采用神經網絡函數代替完整水動力模型求解吸收功率,計算波浪能裝置的最優壓載幾何,與無壓載情況相比,可將輸出功率提升84%。Goggins 等[14]將確定海域的年平均波浪能資源作為輸入,對不同形狀和半徑的浮體的目標函數進行了優化。Esmaeilzadeh 等[15]研究了水下平面壓差型波浪能裝置,結果表明在不規則波條件下,與同等面積的圓形相比,最優形狀浮體在某些工況下的平均轉換功率可提升數倍。

上述研究絕大多數基于頻域模型進行簡化計算,雖然可節省計算成本,但忽視了時域模型的非線性效應對計算精度的影響,進而對形狀優化結果的有效性產生影響。此外,優化過程中浮體形狀定義方法和算法設置對優化結果的影響規律鮮有深入研究。本文擬研究的多自由度波浪能裝置具有更復雜的動力學行為,其浮體的水動力性能對整體裝置的能量轉換效率有顯著影響。為此,本文首先建立該波浪能裝置的時域數值模型,采用遺傳算法搜尋浮體最優形狀,分別將平均發電功率與排水體積和濕表面面積之比作為目標函數,基于B 樣條曲面選用2 種形狀定義方案進行浮體形狀的參數化描述,獲得并分析優化結果。在數值模擬基礎上,開展縮比實驗研究,將優化的形狀與半球形和圓柱形2 種基準形狀進行對比,研究浮體形狀優化設計方法,以提高浮體的能量捕獲效率。

1 多自由度波浪能裝置簡述

目前,波浪能開發領域已基本形成了多種技術共存的格局。按照波浪能裝置的能量捕獲方式分類,有浮子式、振蕩水柱式和越浪式等。浮子式波浪能發電裝置作為較早開發的一種技術,具有結構簡單、技術發展成熟、較易維護等諸多優勢。但是,現有的浮子式波浪能發電裝置受限于自身的結構形式,往往設計成僅將浮體在一個運動方向上的機械能轉換為電能,例如以垂蕩運動為主的裝置只通過浮體與參照物的相對垂蕩運動來吸收垂蕩方向的浮體機械能,而其他方向(縱蕩和橫蕩,縱搖和橫搖)的浮體機械運動卻受到結構設計約束,導致這一部分浮體機械能無法被吸收轉換為電能,大大減少了理論上利用振蕩浮子原理可捕獲的波浪能總量。

本文研究的對象是多自由度波浪能裝置[16],如圖1 所示。在入射波浪的作用下,與浮體連接的機械傳動系統沿軸承方向直線運動,以及繞水平面2 個相互正交的軸線轉動,類似于陀螺儀的運動原理。其中直線運動轉換的是浮體在垂蕩運動模態上的機械能,轉動轉換的是浮體在縱蕩、橫蕩等其他運動模態上的機械能。實際上,該多自由度波浪能裝置各方向的機械傳動機構均可安裝直線或旋轉發電機輸出電能,因此該裝置具有多個PTO 輸出軸。本文的研究重點是浮體形狀優化,在后續仿真和實驗研究中,各自由度的PTO 力均用線性阻尼值代替。

圖1 多自由度波浪能裝置[16]Fig. 1 The multi-DOFs wave energy converter[16]

多自由度波浪能裝置的運動響應包括3 個自由度,如圖1 所示,其運動響應位移分別定義為平動PTO 軸的直線位移z,和2 個轉動PTO 軸的旋轉位移 γ和θ。基于勢流理論和牛頓第二定律,對這3 個自由度方向的運動進行動力學分析,可以得到裝置的時域運動方程為:

式中:m和ma(∞)分別為裝置浮體的質量和無窮大周期下的附加質量矩陣;K為靜水恢復剛度矩陣;τ為計算輻射力的卷積積分自變量;R(t)為遲滯函數,表示為

其中,λ(ω)為與頻率相關的輻射阻尼系數。

2 浮體形狀優化方法

浮體形狀優化包括3 個主要步驟:首先,基于B 樣條曲面的形狀描述方法,采用一系列控制點坐標定義浮體外形曲面形狀;隨后,加入遺傳算法,將所有控制點坐標作為種群中的個體變量;最后,將每一代生成的所有浮體形狀依次應用到基于勢流理論建立的裝置時域數值模型中,計算獲得每一種浮體形狀所對應的目標函數值,通過多次迭代獲得最優控制點坐標集合,確定浮體形狀的最優結果。

2.1 遺傳算法

遺傳算法由Holland 于1975 年提出[17],并由Goldberg 等[18]進一步完善。該算法是一種基于生物遺傳學的自適應全局優化搜索算法,使用群體搜索技術,通過對當前群體施加選擇、交叉、變異等一系列遺傳操作,從而產生新一代的群體,并逐步使群體進化到包含或接近最優解的狀態。該算法的優勢在于:1)同時使用多個搜索點的信息,可實現并行計算;2)以決策變量的編碼作為運算對象,適用于無數值概念或數學建模較困難的優化問題;3)其直接以目標函數值作為尋優依據,無需目標函數的導數值,不要求函數在定義域內連續可微。

將選取的浮體輪廓面控制點坐標值作為變量,采用實數編碼。由于目標函數的計算量遠遠超過遺傳算法本身的計算量,需要減少對目標函數的調用次數,使用變異率較高而數量較小的種群來增加種群多樣性。通常推薦的種群大小是搜索空間維度的1~2 倍,因此在本研究中,將種群大小NP設為20。選擇和交叉操作采用君主方案,即在根據目標函數值高低對種群進行排序的基礎上,將最優個體與其他偶數位的所有個體進行交叉,交叉率設為0.6,每次交叉產生2 個新個體。在交叉過后,對新產生的種群進行多點變異產生子種群,再計算其目標函數值,然后與父種群合并,并根據目標函數值進行排序,取前NP個個體為新種群,進行下一次遺傳操作,最大遺傳代數設為50。通常,建議的變異率為種群大小的倒數,但為了增加本算法中種群的多樣性,將變異率設為0.3。

迭代算法流程如圖2 所示。具體步驟為:

圖2 算法流程圖Fig. 2 Algorithm flow chart

1) 種群初始化,隨機生成NP個個體;

2) 評估個體的目標函數;

3) 確定個體的選擇概率,采用隨機遍歷抽樣法確定用于交叉操作的父代個體;

4) 采用中間重組將父代個體進行隨機交配生成子代;

5) 變異操作;

6) 采用精英主義策略,保留每一代種群中目標函數值最優的個體,與生成的子代共同構成下一代種群;

7) 判斷是否滿足收斂條件,若不滿足,從步驟2)起循環迭代,直到達到收斂,輸出最優解。

在遺傳算法中,迭代產生的每一個浮體形狀的好壞通過計算比較目標函數來判斷,目標函數的選擇對于優化結果的影響十分顯著。

波浪能發電技術的瓶頸之一是度電成本很高,難以達到商業化發展的成熟度。因此,對波浪能發電裝置進行研發優化的目標之一是在保證輸出發電功率和發電可靠性的同時,盡量降低發電成本。浮體作為波浪能裝置的一級獲能機構,浮體質量是其主要成本,因此浮體優化的目標可設定為單位質量下的發電量。浮體成本的評估依據包括浮體尺寸、浮體排水體積/濕表面面積。分別采用排水體積、濕表面面積來表征發電裝置浮體的成本,優化目標是提高波浪能裝置單位排水體積以及單位濕表面面積的發電功率。在此,設定目標函數f1為波浪能裝置的平均發電功率與排水體積之比,目標函數f2為波浪能裝置的平均發電功率與濕表面面積之比,即

2.2 時域數值模型

使用AQWA/WEC-Sim 搭建時域數值仿真模型,WEC-Sim 是模擬波浪能裝置的開源代碼[19]。代碼是在Matlab/Simulink 軟件中使用多體動力學求解器Simscape Multibody 開發。WEC-Sim 需要的唯一外部輸入是來自WAMIT,AQWA,Capytaine,HAMS 等軟件的邊界元水動力數據。基于WEC-Sim 使用這些數據在時域中模擬實際裝置,然后再與控制系統、PTO、其他物體和作用力相互耦合。WEC-Sim 工具箱的代碼結構如圖3所示。

圖3 WEC-Sim 代碼結構圖[17]Fig. 3 WEC-Sim code structure[17]

本文研究的多自由度波浪能裝置的2 個旋轉PTO 軸是對稱的。在仿真研究中,為簡便處理,可以設定波浪的入射角方向與裝置的1 個旋轉PTO軸的軸線朝向相同,則多自由度波浪能裝置在平移PTO 軸和其中1 個旋轉PTO 軸上輸出電能,據此建立如圖4 所示的仿真模型。

圖4 多自由度波浪能裝置的WEC-Sim 仿真模型Fig. 4 WEC-Sim simulation model of the multi-DOFs wave energy converter

2.3 浮體形狀定義

選擇恰當的幾何參數化描述方法來定義算法搜索空間,對于波浪能發電裝置形狀優化的效率是極其重要的。本文優化算法是通過在所允許的范圍內操縱變化定義外形輪廓來獲得各種形狀的波浪能裝置浮體,這一過程需要選擇具有代表性的幾何參數來表征形狀。B 樣條曲面是一種曲面造型參數化設計方法,通過確定控制點集合就能夠生成多種樣條曲面,通過修改控制點位置便可控制整個曲面的形狀,同時曲面具有良好的光滑度。本文基于B 樣條曲面采用2 種方法對浮體形狀進行幾何參數化定義。

B 樣條曲線方程可表示為:

基于B 樣條曲面,根據控制點空間的不同采用2 種形狀定義方案(以下稱形狀方案),且均采用11 個坐標可變的控制點。其中,形狀方案-1 中還有一個坐標恒定不變的控制點P12,用于使控制點連接生成的曲線與坐標軸相交實現閉環。如圖5 所示,形狀方案-1 中浮體外形曲面的控制點分布用圓柱體空間表示,各控制點坐標在z軸方向等距排列,僅在x軸方向變化,其所有坐標值作為基因組成遺傳算法中的基因型,即將控制點P1~P11的x軸坐標值作為遺傳算法的變量組,共11 個變量。形狀方案-2 中浮體外形曲面的控制點分布用多面體空間表示,各控制點坐標分別在z軸與x軸方向變化,其所有坐標值作為基因組成基因型,即將控制點P1~P11的x軸和z軸坐標值作為遺傳算法的變量組,共22 個變量。

圖5 形狀定義方案:控制點空間和示例曲面Fig. 5 Shape definition scheme: control point space and examples

3 數值模擬與分析

本文選取皮爾遜?莫斯柯維奇譜(簡稱P-M 譜)描述的不規則波開展數值模擬,它是由莫斯柯維奇對北大西洋1955-1960 年的海浪觀測資料進行譜分析得到,適用于充分成長的海浪。波譜公式為:

式中:fp為譜峰頻率,此處fp=0.2 Hz,Hm0為有義波高;f為波浪頻率。第i個波浪頻率fi對應的波面高程幅值為:

將不規則波看作由大量不同頻率、不同波幅、不同入射方向的規則波線性疊加而成,則波面高程可表示為:

采用仿真波浪工況生成的波譜及波面高程時歷曲線如圖6 所示。

圖6 仿真波浪工況Fig. 6 Simulated wave conditions

3.1 形狀方案對優化結果的影響分析

在設置目標函數f1(單位排水體積的平均發電功率)的優化仿真算例中,對比分析應用2 種形狀方案的迭代收斂情況和優化結果。圖7 所示為2 種方案各次迭代的最優目標函數值的變化趨勢。由圖可見,2 種方案下的目標函數值進化曲線在數代內均保持不變,其間偶爾會出現較大的跳躍,達到一個更高的目標函數值。之所以會出現目標函數值保持不變的現象,是由于算法采用了君主方案,從父代中保留了最優解,同時在交叉和變異過程中并沒有得到更優解。而曲線中目標函數值出現跳躍的現象則是因隨機的基因變異或交叉導致更優解的出現所致。可見,2 種方案下的目標函數值在有限迭代次數下的收斂效果均較好。

圖7 目標函數值進化曲線Fig. 7 Evolution curves of objective function value

另外,在迭代50 次后,形狀方案-1 與形狀方案-2 的最優目標函數值分別為108.37,112.58,由此得到2 種方案最優浮體形狀的平均發電功率分別為14 725.12 和16 453.69 W。選取相同體積、吃水與高度的圓柱形浮體作為基準形狀,計算獲得其目標函數值,并與優化形狀的目標函數值進行對比,如表1 所示。

表1 基準形狀與優化形狀性能對比Table 1 Comparison of the performance between benchmark shape and optimal shape

由表1 可知,形狀方案-1 與形狀方案-2 的優化形狀相比基準形狀在能量捕獲特性上有了顯著提高,平均發電功率可提高30%以上。同時,在目標函數值、平均發電功率及增幅方面,后者相比前者都較優,說明后者在以考慮體積為浮體成本的前提下能量捕獲特性優于前者。

3.2 目標函數對優化結果的影響分析

應用2 種形狀方案,分別以目標函數f1(單位排水體積的平均發電功率)和目標函數f2(單位濕表面面積的平均發電功率)開展優化仿真,得到的優化形狀分別如圖8 和圖9 所示。

圖8 目標函數f1 下的優化形狀Fig. 8 Optimal shape under objective function 1

圖9 目標函數f2 下的優化形狀Fig. 9 Optimal shape under objective function 2

由圖8 可見,通過形狀方案-1 優化得到的浮體形狀趨向于漸變的規則幾何形狀,浮體寬度在靜水面附近達到控制點坐標優化約束范圍的最大寬度,從靜水面開始隨吃水的增大而減小,一段距離后形狀收斂為與圓柱體相似,再往下時浮體寬度基本保持不變。分析這一段與圓柱體相似的形狀曲線幾何控制點的坐標結果可知,控制點坐標均收斂到坐標約束范圍的下限,即最小值附近。這樣的形狀呈現的規律為:在盡可能獲得更大的水線面面積的同時,浮體的排水體積盡可能趨小。通過形狀方案-2 得到的結果則呈現了“倒駝峰”幾何形狀特征,幾何體兩側向后傾斜形成側翼,首部較為突出。

由圖9 可見,通過形狀方案-1 得到的優化形狀同樣較為規則,但與目標函數f1的結果不同的是,浮體寬度的變化隨吃水深度的增大逐步減小,且變化較平緩。另外,可以發現在目標函數f2下采用2 種形狀方案優化獲得的浮體外形曲面均較平滑,截面面積自上而下逐漸減小,呈現出“倒金字塔”的幾何形狀特征。

將相同條件下的圓柱形和半球形作為基準形狀,計算其目標函數值,對比優化形狀。為保證在同一尺寸規模下對比不同形狀的性能,在目標函數f1下,基準形狀與優化形狀具有相同的排水體積和吃水深度;而對于目標函數f2,除吃水深度相同外,基準形狀與優化形狀的濕表面面積也一樣。表2 所示為得到的各形狀目標函數值。由表2 可以發現,無論采用形狀方案-1 還是形狀方案-2,得到的優化形狀的目標函數值均明顯優于基準形狀,其中,采用形狀方案-2 得到的結果均遠遠大于圓柱形的目標函數值,充分證明了這種形狀定義方案的有效性。

表2 基準形狀與最優形狀對比Table 2 Comparison between benchmark shape and optimal shape

綜上所述,可以發現:無論是形狀方案-1 還是形狀方案-2,優化形狀都呈現出上寬下窄的“倒金字塔”的幾何形狀特征,水線面以下體積的占比較大;形狀方案-1 的結果較為規則,易于制造,成本較低,形狀方案-2 的結果則具有更好的能量捕獲性能。另一方面,不同目標函數下2 種方案的優化形狀存在差異,說明浮體的優化形狀會隨著形狀方案和目標函數變化導致算法中的基因型復雜程度或搜索空間維度的不同而變化。

4 實驗研究

為驗證本文提出的多自由度波浪能裝置形狀優化方法的有效性,在浙江大學海洋學院小型波浪水槽中開展了縮比例實驗研究。實驗布置方案如圖10 所示。波浪水槽長25 m,水槽寬0.7 m,最大深度0.7 m,最大實驗允許水深0.5 m。在滿足造波質量的情況下,水槽造波波浪周期范圍為0.9~4.0 s,波高范圍為0.02~0.10 m。

圖10 實驗布置圖Fig. 10 Experimental layouts

將目標函數f1獲得的優化形狀浮體作為實驗組,將圓柱形浮體和半球形浮體這2 種基準形狀作為對照組。根據水槽尺寸,小比例浮體模型的寬度均為300 mm,比例尺為0.15,采用三維樹脂打印技術加工制作。在裝置各個PTO 軸方向上布置有減速電機來模擬負載阻尼力,負載力可通過電機輸出側的電阻值來調節。實驗中使用0.25,0.50,0.75,1.00 和1.25 Ω 這5 種負載電阻來產生負載阻尼力。浮體在平移PTO 軸上的直線相對運動通過滾珠絲桿和聯軸器等機械傳動部件轉換為旋轉運動再連接電機,浮體在旋轉PTO 軸上的旋轉相對運動通過增速齒輪箱連接電機。圖11 所示為實驗結果。

由圖11 可以看出,在所有負載情況下,優化得到的浮體運動響應比基準形狀更優,驗證了本文形狀優化方法的有效性。在負載電阻配置為0.75 Ω 時平均發電功率達到最大值,在配置電阻值遠離0.75 Ω 時平均發電功率降低,這可能是因為配置0.75 Ω 電阻時產生的負載阻尼力落在了實驗波況條件下模型裝置的最佳阻尼附近,因此具有最好的俘能效率。同時,當電阻值偏小或者偏大時,也就是負載阻尼相對于最優負載阻尼偏小或者偏大時,優化形狀相比基準形狀在轉換效率上的性能提升不明顯。這表明優化結果同時受到PTO 負載力的影響。在實際海況中,波浪能裝置的PTO 負載阻尼力不是簡單線性的,甚至是時變的,需在優化計算中增加考慮PTO 負載阻尼力對優化結果的影響,這是后續將開展的工作。

圖11 實驗結果Fig. 11 Experimental results

5 結 論

本文研究了不規則波下多自由度波浪能裝置浮體形狀的優化設計問題,設計了一種基于遺傳算法和B 樣條曲面的優化設計方法。該方法將平均發電功率與排水體積和濕表面面積之比作為目標函數,基于B 樣條曲線選用了2 種形狀方案來進行浮體外形的參數化描述。通過數值和實驗研究將優化得到的浮體形狀與基準形狀的目標函數值進行對比,分析了不同形狀方案和目標函數對優化結果的影響。可得到以下結論:

1) 優化形狀的能量捕獲特性比基準形狀有顯著提高;與基準形狀的平均發電功率相比,2 種形狀方案所得優化形狀的平均發電功率增幅都超過了30%。

2) 2 種形狀方案得到的優化浮體形狀都呈現出了上寬下窄的幾何形狀特征,水線面以下體積的占比較大;形狀方案-1 的優化結果浮體曲面較規則,形狀方案-2 的優化結果浮體曲面較平滑。

3) 與形狀方案1 相比,形狀方案-2 在目標函數值、平均發電功率及增幅方面都較優,說明形狀方案-2 在保證成本的前提下其能量捕獲特性優于形狀方案-1。同時,形狀方案-1 獲得的優化形狀較為規則,具有更易于加工的優點。

猜你喜歡
優化
超限高層建筑結構設計與優化思考
房地產導刊(2022年5期)2022-06-01 06:20:14
PEMFC流道的多目標優化
能源工程(2022年1期)2022-03-29 01:06:28
民用建筑防煙排煙設計優化探討
關于優化消防安全告知承諾的一些思考
一道優化題的幾何解法
由“形”啟“數”優化運算——以2021年解析幾何高考題為例
圍繞“地、業、人”優化產業扶貧
今日農業(2020年16期)2020-12-14 15:04:59
事業單位中固定資產會計處理的優化
消費導刊(2018年8期)2018-05-25 13:20:08
4K HDR性能大幅度優化 JVC DLA-X8 18 BC
幾種常見的負載均衡算法的優化
電子制作(2017年20期)2017-04-26 06:57:45
主站蜘蛛池模板: 欧美精品亚洲二区| 中文字幕不卡免费高清视频| 尤物午夜福利视频| 欧美日韩第三页| 九九九精品成人免费视频7| 亚洲色图欧美| 亚洲一区二区三区香蕉| 尤物成AV人片在线观看| 日韩人妻无码制服丝袜视频| 久久综合九色综合97网| 无码啪啪精品天堂浪潮av| 久久人人爽人人爽人人片aV东京热| 精品亚洲国产成人AV| 亚洲色图欧美激情| 日韩国产无码一区| 国产欧美日韩18| 97国产在线视频| 欧美国产视频| 666精品国产精品亚洲| 青青操国产| 国产丝袜一区二区三区视频免下载| 国产精品自在拍首页视频8| 亚洲第一色网站| 亚洲a级在线观看| 免费女人18毛片a级毛片视频| 亚洲国产精品国自产拍A| 狠狠色噜噜狠狠狠狠色综合久 | 亚洲最猛黑人xxxx黑人猛交 | 特级毛片8级毛片免费观看| 欧美视频免费一区二区三区| 操国产美女| 性欧美在线| 波多野结衣的av一区二区三区| 国产成人亚洲毛片| 精品少妇人妻av无码久久| 国产精品思思热在线| 青青久久91| 91精品综合| 无码乱人伦一区二区亚洲一| 呦女亚洲一区精品| 人妻丰满熟妇av五码区| 91麻豆精品国产91久久久久| 日韩性网站| 亚洲AV无码不卡无码| 九九精品在线观看| 免费看黄片一区二区三区| 国产日韩久久久久无码精品| 久久成人国产精品免费软件| 99re在线观看视频| 在线99视频| 国产精品永久免费嫩草研究院| 成人久久精品一区二区三区| 国产91丝袜| 色男人的天堂久久综合| 热这里只有精品国产热门精品| 中文字幕在线播放不卡| 成AV人片一区二区三区久久| 狠狠综合久久| 色综合中文字幕| 国产精品亚洲一区二区三区z | 2019年国产精品自拍不卡| 国产区精品高清在线观看| 五月激情婷婷综合| 日韩色图区| 波多野结衣一二三| 午夜欧美理论2019理论| 美女一级免费毛片| 在线国产三级| P尤物久久99国产综合精品| 亚洲精品另类| 97国产在线视频| 日本三区视频| 亚洲国产成人超福利久久精品| 国产一区在线视频观看| 亚洲中文字幕av无码区| 青青草原国产精品啪啪视频| 99尹人香蕉国产免费天天拍| 亚洲婷婷丁香| 日韩毛片免费| 青青操国产| 日本成人在线不卡视频| 久久伊人操|