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

電機-變距螺旋槳動力系統功率優化控制

2021-03-27 04:47:22段登燕裴家濤祖瑞李建波
航空學報 2021年3期
關鍵詞:優化模型

段登燕,裴家濤,祖瑞,李建波

南京航空航天大學 直升機旋翼動力學國家級重點實驗室,南京 210016

電動螺旋槳無人機應用于航拍、農藥噴灑、電力巡線等越來越普及,但由于電池能量密度低,此類無人機續航時間一般較短、載重能力低,限制了其應用領域的擴展[1-4]。在期待電池技術進一步突破的同時,不斷挖掘電機和螺旋槳的性能,減少功率消耗,提高飛行性能,是目前螺旋槳電動無人機普遍面臨的技術問題。電機-變距螺旋槳動力系統(以下簡稱變距電動力系統)可同時改變轉速、槳距兩個量,存在槳距和轉速的最佳組合,使系統功率最小[5-6]。相比電機-定距螺旋槳動力系統,其在耗能方面具有特殊優勢,但如何達到最小功率點,目前研究較少。

針對變距電動力系統功率優化控制問題,國內研究較少,國外處于理論研究及試驗階段。Cazenave等設計了一種基于卡爾曼濾波估計的變距電動力系統最小功率控制策略,該方法在模型具有不確定性及外界存在小擾動時仍適用[7-8],而采用的卡爾曼濾波算法不具備自適應性,不能對模型參數和噪聲特性進行在線估計及修正,可能會導致誤差較大。Cohen等提出了等步長、變步長兩種變距電動力系統最小功率控制方案并進行了對比,結果表明變步長方案更優[9],但所提變步長策略過于簡單,對不同系統不易自適應地改變步長,可能導致前期收斂速度慢的問題。Fresk[10]、Sheng[11]等針對四軸無人機開展了變距電動力系統最小功率控制工作,但螺旋槳均處于零入流狀態下,不適合目前研究較多的固定翼電動無人機、小型電動傾轉旋翼機、分布式電動無人機等螺旋槳存在較大入流狀態的飛行器。Henderson和Papanikolopoulos[12]針對某小型固定翼電動無人機,螺旋槳處于某一固定入流速度時,采用擾動與觀測(P&O)算法實現槳距優化,以PID控制調整螺旋槳轉速,設計了功率優化控制策略,但沒有對螺旋槳入流發生變化、拉力需求改變或外部存在擾動時,檢驗所提算法的有效性。

此外,在變距電動力系統建模中,電機建模以往包括一階電動機模型和電機損耗模型兩種方式[13]。等效回路模型較為直觀,但沒有計入電機的風損、摩擦損失、鐵損等,不能較好地反映電機實際效率。電機損耗模型可以較好地反映電機實際效率,但其建模依據的多個參數需由經驗或實驗獲得[14-16]。螺旋槳一般采用葉素理論、渦流理論開展建模[7-12],但均包含迭代過程,計算速度慢,不利于控制策略的研究及優化設計工作的開展。針對變距電動力系統功率優化控制問題,本文首先將電機一階電動機模型、損耗模型結合,建立了無刷直流電機模型;依據Goldstein渦流理論建立了螺旋槳模型;在此基礎上,為提高計算效率,基于改進天牛須算法的BP(Back Propagation)神經網絡辨識得到了變距電動力系統神經網絡代理模型。接著提出了一種變距電動力系統功率優化控制策略:在一定入流速度、拉力需求下,基于自適應擴展卡爾曼濾波-牛頓法實時優化槳距,并在一定槳距下利用模糊PID控制系統轉速以達目標拉力,實現目標拉力需求下的最小功率控制。仿真驗證結果表明,提出的功率優化控制策略魯棒性更強、優化速度更快、收斂效果更好。

1 變距電動力系統建模

1.1 無刷電機建模

一階電動機模型較為直觀,但沒有計入電機損失;電機損耗模型可以較好地反映電機實際工作效率,但其建模依據的多個參數需由經驗或實驗獲得。為了得到電機損耗模型的各個參數,本文將等效回路模型電機最大效率及效率為0.8時的計算數據,與電機損耗模型結合,開展電機的建模。

圖1 等效回路模型Fig.1 Equivalent circuit model

圖1為電機的等效回路模型,其中Vs為電池電壓,V;V為輸入電機的電壓,V;Vm為電機在相應轉速下的反電動勢,V。電調通過脈寬調制技術(PWM)調節V的大小,當β=100%時電機處于最大轉速。電機輸入電壓與電池電壓的關系為

V=βVs

(1)

電機有3個特征參數,Kv為電機轉速系數,r·min-1·V-1;R為內阻,Ω;i0為空載電流,A。由式(2)和式(3),可分別得到電機轉速ω、電機扭矩Q,即

(2)

(3)

由以上參數可得電機出軸功率Pshaft、入軸功率Pele、電機效率η和電機損失PL。

Pshaft=Qω=(V-iR)(i-i0)

(4)

Pele=Vi

(5)

(6)

PL=C0+C1ω+C2ω3+C3Q2

(7)

式中:Ci(i=0,1,2,3)為電機損失相關系數。

如式(7)所示,PL為電機損失,主要由摩擦損失、風阻損失、銅損、鐵損組成。其中摩擦損失和鐵損與ω成正比,見式(7)第2項;風阻損失與ω3成正比,見式(7)第3項;銅損與Q3成正比,見式(7)第4項。

由式(8)可求得電機效率η:

(8)

由式(8)分別對Q、ω求偏導可得

(9)

(10)

為了求最大效率,令偏導等于0,可得

(11)

(12)

(13)

(14)

將C1、C2、C3用C0表示,取電機效率為0.8時的轉速ω0.8、扭矩Q0.8及功率PL,0.8。由式(14) 可得C0,進而得到C1、C2、C3。

XM3電機特征參數如表1所示,從圖2電機效率隨電流的變化圖可看出,本文計算值相比一階電動機模型,與實驗值更接近,誤差更小。

以電機轉速Ω為x軸,電機扭矩為y軸,電機效率為z軸做XM3電機的效率圖如圖3。從圖中可以看出,電機在一定轉速及扭矩范圍內具有較大效率,超過了這個范圍,電機效率會有所下降。

表1 XM3電機特征參數Table 1 Characteristic parameters of motor XM3

圖2 電機效率隨電流的變化Fig.2 Motor efficiency with current change

圖3 XM3電機效率圖Fig.3 Efficiency map of motor XM3

1.2 螺旋槳建模

依據Goldstein渦流理論[17-18]對螺旋槳開展建模,圖4為作用在葉素上的速度和力。R為半徑,m;r為槳轂中心到槳葉剖面任一點處的距離,x為r的無量綱值。σ為槳葉實度;B為槳葉片數;b為葉素弦長,m;ω為螺旋槳轉速,rad/s;v為來流速度,m/s;λ為來流速度的無量綱值;φ為葉素來流角;φT為槳尖處的葉素來流角;αi為干涉角;vE為合速度;ωa、ωt分別為誘導速度的軸向和環向分量;α為葉素迎角;vR為來流速度v和旋轉速度ωr的合成量。

(15)

定義Γ為葉素環量,κ為Goldstein系數,以下計算中用Prandtl葉尖損失系數F來代替。

(16)

(17)

(18)

圖4 作用在葉素上的速度和力Fig.4 Velocity and force acting on blade element

由圖4的幾何關系,可得式(19),αi為干涉角。

(19)

(20)

Δα、Δθ需要引入以下修正,h為葉素厚度:

(21)

(22)

(23)

由式(19)和式(21)可得r處的迎角

α=θ+β-φ-αi-Δα

(24)

式中:θ為槳距;β為r處的安裝角。

拉力系數CT、功率系數CP可由式(25)和式(26) 得到。此外,翼型的雷諾數對升阻力系數的影響很大。本文采用二維差值方法計算不同雷諾數下的升阻力系數。

(CLcos(φ+αi)-CDsin(φ+αi))dx

(25)

(CLsin(φ+αi)+CDcos(φ+αi))dx

(26)

作電機螺旋槳匹配圖,見圖5,圖中3條粗實線代表不同螺旋槳不同工作狀態下,螺旋槳扭矩隨轉速的變化。可以看出,螺旋槳扭矩較小處,電機效率可能比較低;電機扭矩較高處,螺旋槳功率可能較高。而螺旋槳功率與扭矩、轉速均有關,這增加了電動力系統最小功率控制的難度。

圖5 電機螺旋槳匹配Fig.5 Motor propeller matching

1.3 變距電動力系統代理模型

由1.1節和1.2節電機、螺旋槳建模過程可以看出,變槳距電動力系統原始模型中存在迭代過程,計算速度慢,不利于后續開展控制及在線優化工作。因此,本節提出用神經網絡代理模型逼近原始模型的策略,并在傳統BP神經網絡的基礎上引入改進天牛須算法,以期達到全局收斂性更好、收斂速度更快的效果。

1.3.1 改進天牛須算法

天牛須(BAS)算法適合于單目標和多目標優化,與遺傳算法(GA)、粒子群優化(PSO)算法相似,不需要知道優化對象模型,可以實現自動尋優;且BAS只有一個個體,相比GA、PSO,尋優速度更快。天牛須算法原理歸納為:天牛覓食時,雖然不知道食物具體方位,但可以通過兩只長觸角辨別食物的氣味強度,如果左邊氣味強度大于右邊,天牛向左邊飛,反之向右飛。

傳統天牛須算法建模步驟如下:

步驟1創建天牛須朝向的隨機向量并作歸一化處理:

(27)

式中:rands(·)表示隨機函數;m為空間維度。

步驟2對天牛左、右須分別創建空間坐標:

(28)

式中:xr,t表示t次迭代時天牛右須的位置;xl,t表示t次迭代時天牛左須的位置;xt表示t次迭代時天牛質心的位置;d表示天牛左、右須的距離。

步驟3根據事先確定的適應度函數來模擬食物氣味的強度。

步驟4更新天牛質心位置:

xt+1=xt-δtbsgn(f(xr,t)-f(xl,t))

(29)

式中:δt表示t次迭代時的步長;sgn(·)為符號函數。

由步驟4可以看出,傳統天牛須算法中,天牛質心位置更新取決于當前迭代f(xr,t)、f(xl,t)的大小,容易偏離最優值。因此本文提出將歷史最優值xbest引入天牛質心位置更新的策略:

xt+1=xt-δtbsgn(f(xr,t)-f(xl,t))+

rands(·)xbest

(30)

1.3.2 BP神經網絡

BP神經網絡是典型多層前向網絡,如圖6所示,主要由信號正向傳播和誤差反向傳播組成。圖中:m為輸入量維度;l為隱含層神經元的個數;n為輸出層神經元個數;IW1,1、IW2,1分別為隱含層、輸出層神經元的權重;b1、b2分別為隱含層、輸出層神經元的閾值;p1為輸入量;a1為隱含層神經元的輸入量;a2為輸出層神經元的輸出量;f為激勵函數。正向傳播指信號依次經過輸入、隱層神經元,再由輸出神經元輸出的過程。反向傳播是指,當輸出值與目標值之間存在誤差時,將誤差反向傳播,并由梯度下降算法更新網絡權值、閾值的過程。當誤差足夠小時,網絡反向傳播停止,BP神經網絡訓練完成。

圖6 BP神經網絡Fig.6 BP neural network

1.3.3 基于改進天牛須算法的BP神經網絡系

統辨識

BP神經網絡在系統辨識、數據擬合、預測等多方面具有廣泛的應用。應用過程中,一般是先給定一組初始權值和閾值,然后由反向傳播不斷更新網絡的權值和閾值,直至達到網絡訓練目的。但這往往導致神經網絡的收斂過程依賴于初始權值和閾值,且容易陷入局部最優。以往有將GA、PSO等優化策略,用于優化初始權值、閾值的研究,并取得了較好的效果。研究發現用BAS代替GA、PSO后,收斂速度更快、全局收斂性一致、內存消耗更小[19]。本文將改進BAS算法用于BP神經網絡初始權值、閾值的優化,以期達到與BAS算法相比,全局收斂性一致、收斂速度更快的效果。基于改進天牛須算法的BP神經網絡系統辨識過程如圖7所示。

圖7 改進BAS-BP算法流程圖Fig.7 Flow chart of updated BAS-BP algorithm

1.3.4 變距電動力系統代理模型

建立代理模型前,需要獲得足夠描述變距電動力系統的數據點。本文以XM3電機-PROP1螺旋槳為例,開展代理模型的建立及驗證工作。PROP1螺旋槳的參數如表2所示。PROP1螺旋槳扭角、弦長、厚度沿槳葉徑向分布如圖8所示;0~0.5R、0.5R~0.75R、0.75R~R(R為半徑)翼型依次為S8037、NACA4412、ClarkY。以ClarkY翼型為例,給出了升阻力系數隨迎角、雷諾數(Re)的變化,如圖9所示。在螺旋槳建模過程中,采用二維插值的方法得到槳葉某處的升阻力系數。

螺旋槳槳距取-4°~10°,入流速度取0~20 m/s,轉速取2 000~8 000 r/min,計算各狀態下螺旋槳拉力、功率、扭矩,進而由螺旋槳轉速、扭矩得到電機效率。經過上述過程,共得到1 120個數據點,選取1 000個數據點用于訓練神經網絡,120個數據點用于驗證工作。

基于改進天牛須算法的BP神經網絡系統辨識,開展變距電動力系統代理模型的建立。螺旋槳代理模型建立時,輸入量為螺旋槳轉速、槳距、入流速度,輸出量為螺旋槳拉力、功率、扭矩。因此,選取輸入層神經元個數為3,輸出層神經元個數為3。依據模型復雜度,選取隱含層神經元個數為10。

表2 PROP1螺旋槳參數Table 2 Parameters of propeller PROP1

圖8 扭角、弦長、厚度分布Fig.8 Distribution of twist angle, chord length and thickness

圖9 升力、阻力系數Fig.9 Lift and drag coefficients

圖10為螺旋槳代理模型訓練時的迭代收斂過程,縱坐標自適應值表示上述1 000個數據點,經過訓練后的誤差平方和。從圖中可以看出,改進天牛須與天牛須相比,收斂速度更快,全局性基本一致。注意到自適應值沒有收斂到0,這與隱層神經元個數、迭代步數等有關。但對變距電動力系統而言,自適應值已足夠小。

圖11為以拉力為例,螺旋槳代理模型訓練完成后,用剩余120個數據點驗證的效果。從圖中可以看出,各數據點擬合較好。

電機代理模型建立時,輸入量為螺旋槳轉速、扭矩,輸出量為電機效率。因此,選取輸入層神經元個數為2,輸出層神經元個數為1。依據模型復雜度,選取隱含層神經元個數為10。圖12為電機效率神經網絡擬合效果圖。

圖10 迭代收斂過程Fig.10 Iterative convergence process

圖11 神經網絡擬合效果:拉力Fig.11 Fitting effect of neural network on force

圖12 神經網絡擬合效果:電機效率Fig.12 Fitting effect of neural network on motor efficiency

2 變距電動力系統最小功率控制策略

以XM3電機-PROP1螺旋槳,入流速度5 m/s、拉力需求20 N為例,在任一槳距下,調整螺旋槳轉速以滿足拉力需求,并記錄當前功率、電機效率,如圖13所示。從螺旋槳功率隨槳距的變化線可以看出,存在槳距和轉速的最優組合,使得螺旋槳功率具有最小值。從電機效率變化曲線可以看出,電機效率隨著槳距的增加而增加,但幅度變化不大,電機效率一直處于較高范圍內。電機功率是螺旋槳功率除電機效率所得的值,是電動力系統的實際消耗功率。從圖中可以看出,與螺旋槳功率曲線類似,電機功率存在最小值。因此,在變距電動力系統處于一定入流速度、拉力需求時,可以以系統消耗功率最小為目標,實時優化槳距并控制轉速,實現目標拉力需求下的功率優化控制。此外,從同等拉力螺旋槳轉速與槳距組合關系可以看出,隨著槳距的增加,為達到目標拉力,轉速減小;槳距超過10°后,轉速增加,是螺旋槳槳葉翼型達到失速迎角所致。

圖13 功率、電機效率、螺旋槳轉速圖Fig.13 Diagram of power, motor efficiency and rotational speed of propeller

為此,本文設計了如圖14所示的變距電動力系統最小功率控制策略,采用牛頓法在線優化槳距,同時利用自適應擴展卡爾曼濾波在線估計牛頓法中需要的梯度、海森矩陣;得到槳距后,輸送給變距電動力系統的神經網絡代理模型,得到螺旋槳拉力T和總消耗功率P;比較當前螺旋槳拉力T與目標拉力T0,采用模糊PID控制在線調整螺旋槳轉速,直至實際與目標拉力相等。

圖14 變距電動力系統最小功率控制策略Fig.14 Minimum power control strategy of variable-pitch propulsion system

2.1 基于自適應擴展卡爾曼濾波-牛頓法的槳距在線優化

針對變距電動力系統,由圖13可以看出,一定拉力需求、一定入流速度下電動力系統功率隨槳距的變化曲線,存在最小值。考慮到牛頓法收斂速度快,本文提出采用牛頓法在線優化槳距的策略。牛頓法公式為

(31)

將功率P寫為泰勒級數展開的形式:

(32)

式中:P′(θk)、P″(θk)、P?(θk)依次為P在θk處的一階、二階、三階導數,o為高階項。

從上述分析可以看出,牛頓法只需要一階、二階導數,此處引入三階導數是為了提高卡爾曼濾波估計的準確度。取θ=θk-1得到:

(33)

式中:δPk=Pk-Pk-1,δθk=θk-θk-1。

ΔPk=Hkxk+vk

(34)

(35)

觀測矩陣Hk展開為

Hk=

(36)

式中,為了保證泰勒展開式的有效性,δθk應取較小值;為了保證較好的觀測性,δθs,1、δθs,2應取較大值。

vk為測量噪聲,其協方差矩陣為R。將P′(θk)、P″(θk)、P?(θk)按泰勒級數展開得:

(37)

將式(37)寫成狀態方程的形式:

(38)

(39)

式中:ωk為測量噪聲,其協方差矩陣為Q。

(40)

(41)

式中:b為遺忘因子,取值在0.95~0.99之間。

2.2 基于模糊PID的轉速控制

變距電動力系統轉速控制時,螺旋槳入流速度、拉力需求、槳距會有變化,但傳統PID控制中參數不能實時調整,控制效果有時會較差。因此,本文采用模糊PID實現螺旋槳轉速控制,以達到某一槳距下,一定入流速度下的拉力需求,如圖15 所示。

圖15 模糊PID控制Fig.15 Fuzzy PID controller

模糊PID控制包括模糊化,確定模糊規則,解模糊等組成部分。以輸入量E為例,本文設計的隸屬度函數如圖16所示。其中,NB、NM、NS、ZO、PS、PM、PB分別表示“負大”“負中”“負小”“零”“正小”“正中”“正大”。

圖16 E的隸屬度函數Fig.16 Membership function of E

模糊規則在以往許多研究中均有提及,此處不再贅述。采用重心法解模糊:

(42)

式中:z0為模糊控制器輸出量解模糊后的精確值;zi為模糊控制量論域內的值;μ(zi)為zi的隸屬度值。

輸入螺旋槳模型中的控制量Ω由式(43)得到:

(43)

式中:Ω(k)為輸出量;e(k)為k時刻的偏差;Ts為采樣周期;KP、KI、KD依次為比例、積分、微分系數。

3 仿真結果

以XM3電機-PROP1螺旋槳,入流速度5 m/s、拉力需求50 s前為20 N,50 s后為30 N為例,開展變距電動力系統的仿真驗證。0~50 s需用拉力固定為20 N,采用上述算法實現功率最小,驗證算法的有效性;50 s后,拉力需求由20 N變為30 N,關注當拉力需求變化時,算法能否及時調整,及調整完成后能否再次實現功率最小,驗證本文算法的魯棒性。

3.1 槳距在線優化

基于自適應擴展卡爾曼濾波-牛頓法結合的槳距在線優化方法,優化螺旋槳槳距。

定義搜索步長如下:

(44)

為了說明基于自適應擴展卡爾曼濾波-牛頓法結合的槳距在線優化方法的有效性,針對以下3種方案進行了對比分析。

方案1step=0.5°,由槳距θ=-4°開始,逐步搜索最小功率點。

方案2初始步長step=0.5°,若找到最小點所在區間,則根據式(45)逐漸減小步長,直至找到最小功率點。

step=step·0.5

(45)

方案3擴展卡爾曼濾波(EKF)-牛頓法。

下面依次給出了在線優化過程中電動力系統總消耗功率、槳距、步長的變化過程。

由圖17可以看出,基于自適應擴展卡爾曼濾波-牛頓法(本文方案)、方案3收斂速度明顯大于方案1和方案2。在10 s、50 s左右,方案3有一定的波動,基于自適應擴展卡爾曼濾波-牛頓法收斂平緩。在最終收斂階段,自適應擴展卡爾曼濾波和方案2收斂性一致,方案1、方案3未收斂到最優點,且存在波動。這主要與步長變化有關,見圖18,基于自適應擴展卡爾曼濾波-牛頓法、方案3初始階段步長較大,所以收斂速度快;10 s、50 s左右,方案3步長有較大負值,造成方案3功率曲線的波動;方案1與方案2初始階段步長恒定,均為0.5,所以初始收斂階段相同,后期方案2步長逐漸減小,最終收斂效果大于方案1。

圖17 變距電動力系統總功率變化Fig.17 Total power change of variable-pitch propulsion system

圖19為優化過程中槳距的變化。與總消耗功率變化類似,基于自適應擴展卡爾曼濾波-牛頓法、方案3收斂速度明顯大于方案1和方案2;方案1最終收斂階段存在一定的波動。方案3最終收斂效果比基于自適應擴展卡爾曼濾波-牛頓法差。

圖18 步長變化Fig.18 Change of step length

圖19 槳距變化Fig.19 Change of pitch

3.2 轉速控制

為了更好地說明轉速控制的魯棒性,取3.1節仿真時間49~51 s為例,開展轉速控制仿真驗證。轉速控制仿真結果如圖20所示,從圖中可以看出,模糊PID控制與傳統PID控制相比,響應時間更短、穩態誤差更小。50 s時,目標拉力由20 N變為30 N,除初始階段有波動外,模糊PID仍具有很好的跟隨性。這是由于模糊PID中,其比例、積分、微分系數可以根據誤差、誤差變化率的大小,進行在線自適應調整。圖21以KP為例,給出了其隨時間的變化曲線。

圖20 轉速控制仿真結果Fig.20 Simulation results of rotational speed control

圖21 Kp隨時間的變化Fig.21 Change of Kp during simulation

4 結 論

1) 將歷史最優值引入天牛須算法的更新過程,提出了改進天牛須算法,并與BP神經網絡結合對變距電動力系統開展系統辨識工作,得到了系統神經網絡代理模型,代替了傳統計算方法中的迭代求解過程,計算效率更高,便于控制策略的仿真驗證。

2) 提出了一種變距電動力系統功率優化控制策略:基于自適應擴展卡爾曼濾波-牛頓法實時優化槳距,并在一定槳距下利用模糊PID控制系統轉速以達拉力需求,實現目標拉力需求下的最小功率控制。仿真結果表明,該策略魯棒性強,收斂速度快。

3) 本文提出的變距電動力系統功率優化控制策略有一定的實用性。應用于電動螺旋槳無人機上時,可以與機載飛控系統組合。變距電動力系統總消耗功率可由電機輸入電壓、電流相乘得到,但由于不易布置拉力傳感器,可以考慮將機載飛控系統得到的無人機速度與拉力聯系起來。根據變距電動力系統在無人機上的分布,由實際速度與目標速度的差異,實現螺旋槳轉速的調節。

猜你喜歡
優化模型
一半模型
超限高層建筑結構設計與優化思考
房地產導刊(2022年5期)2022-06-01 06:20:14
民用建筑防煙排煙設計優化探討
關于優化消防安全告知承諾的一些思考
一道優化題的幾何解法
由“形”啟“數”優化運算——以2021年解析幾何高考題為例
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
主站蜘蛛池模板: 精品一区二区三区自慰喷水| 午夜少妇精品视频小电影| 青草午夜精品视频在线观看| 国产欧美日韩在线一区| 免费毛片在线| 在线观看国产黄色| 久久久久久国产精品mv| 日韩精品成人在线| 午夜色综合| 亚洲天堂网2014| 免费99精品国产自在现线| 91成人试看福利体验区| 欧美五月婷婷| 国产喷水视频| 欧美一级大片在线观看| 日韩色图在线观看| yjizz国产在线视频网| 香蕉在线视频网站| 免费一级毛片在线播放傲雪网| www欧美在线观看| 中文字幕精品一区二区三区视频| 国产精品无码AⅤ在线观看播放| 免费播放毛片| 国产亚洲一区二区三区在线| 亚洲精品色AV无码看| 色综合五月婷婷| 成人永久免费A∨一级在线播放| 亚洲人成影视在线观看| 欧美性久久久久| 香蕉eeww99国产在线观看| 久夜色精品国产噜噜| 五月六月伊人狠狠丁香网| 国产一级小视频| 亚洲欧美精品在线| 最新亚洲人成无码网站欣赏网| 国产成人久视频免费| 成人国产三级在线播放| 日本久久免费| 亚洲第一色视频| 国产激爽大片高清在线观看| 亚洲男人的天堂在线观看| 久久久久无码精品| 亚洲国产精品人久久电影| 欧美日韩精品一区二区在线线 | 日韩二区三区无| 高清免费毛片| 自慰高潮喷白浆在线观看| 国产精品无码作爱| 亚洲欧洲美色一区二区三区| 国产亚洲欧美在线中文bt天堂| 亚洲欧洲日韩久久狠狠爱| 国产精品女人呻吟在线观看| 伊人国产无码高清视频| 国产一级无码不卡视频| 91毛片网| 国产精品毛片在线直播完整版| 免费看av在线网站网址| 好紧好深好大乳无码中文字幕| 欧美国产在线一区| 女人av社区男人的天堂| 试看120秒男女啪啪免费| 欧美日韩资源| 狼友视频一区二区三区| 国产91精品久久| 久久人妻xunleige无码| 国产色婷婷| 在线日韩日本国产亚洲| 特级aaaaaaaaa毛片免费视频| 国产精品天干天干在线观看| 蜜桃视频一区二区| 日韩欧美国产精品| 国产欧美日韩在线在线不卡视频| 精品国产自| 亚洲综合18p| 亚洲美女高潮久久久久久久| 国产女人综合久久精品视| 四虎影视无码永久免费观看| 国产精品视频导航| 97视频免费在线观看| 欧美日韩在线成人| 色男人的天堂久久综合| 国产微拍一区|