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

雙穩態壓電振動能量采集器的時-頻域動力學特性及實驗研究*

2019-09-21 08:00:18趙澤翔王光慶譚江平崔素娟武海強
傳感技術學報 2019年8期
關鍵詞:振動

趙澤翔,王光慶,譚江平,崔素娟,武海強

(浙江工商大學信息與電子工程學院,杭州 310018)

雙穩態壓電振動能量采集器是一種將環境中振動能采集并轉化成電能的非線性機電耦合器件,具有工作頻帶寬、轉換效率高等優點,在無線傳感器網絡節點和微電子設備的供電系統中具有良好的應用前景[1-6]。

國、內外對雙穩態壓電振動能量采集器進行了大量研究,并取得了較豐富的研究成果。Minghui Y A[7]等人提出了雙穩態L形梁結構,理論分析得到該系統的勢能函數具有兩個對稱的穩定勢阱。此外,當考慮重力勢能時,勢能函數則出現不對稱的雙勢阱。李海濤[8]等人利用Melnikov方法建立了雙穩態壓電模型并分析了系統的同宿分岔和混沌等非線性動力學行為,研究結果表明Melnikov方法可為能量采集器的參數設計提供有效的理論依據。Sovan Sundar Dasgupta[9]等人利用Fokker-Planck-Kolmogorov方程,研究了高斯白噪聲激發的非線性雙穩態能量采集器的系統參數對振蕩器的位移和速度的概率分布以及輸出電壓的均方的影響,結果表明振蕩器的位移和速度的聯合概率分布的最大峰值隨著噪聲強度值的增大減小。

盡管如此,雙穩態壓電振動能量采集器還有很多科學問題需要深入研究,如雙穩態壓電振動能量采集器的時-頻域動力學響應特性對系統參數的敏感性問題,它是優化雙穩態壓電振動能量采集器結構以及提升其綜合輸出性能的關鍵因素。目前,針對雙穩態壓電振動能量采集器動力學響應特性的研究要么側重于從時域角度去分析系統參數對其動力學特性的影響,要么側重于從頻域角度分析雙穩態壓電振動能量采集器在系統參數作用下的頻率響應特性。這些都不能全面揭示雙穩態壓電振動能量采集器的非線性振動激勵和動力學響應特性,以至于無法對雙穩態壓電振動能量采集器進行結構優化參數設計和綜合性能提高。本文分別從時域和頻域兩個角度綜合研究分析系統參數對雙穩態壓電振動能量采集器動力學響應特性的影響,為揭示雙穩態壓電振動能量采集器的動力學特性和非線性作用機制提供全方位的視角。本文考慮振型和分布參數對輸出性能的影響,利用拉格朗日方程和Raleigh-Ritz方法建立了雙穩態壓電懸臂梁振動能量采集器分布參數非線性機電耦合動力學模型,并利用龍格庫塔算法(ODE45)對雙穩態壓電懸臂梁振動能量采集器進行仿真,研究激振加速度、磁鐵對間距、激振頻率對系統動力學特性的影響;此外,利用諧波平衡法對雙穩態壓電懸臂梁振動能量采集器進行頻域分析,研究激振加速度、磁鐵對間距、激振頻率、負載阻抗對系統采集輸出電特性(電壓、電功率)的影響,比較了不同激勵頻率、磁鐵對間距對最優阻抗的影響,并利用實驗進行驗證。

1 雙穩態壓電振動能量采集器建模

圖1 雙穩態壓電振動能量采集器結構示意圖

考慮圖1所示的雙穩態壓電振動能量采集器結構為主要研究對象,它主要由壓電雙晶片懸臂梁和兩個磁鐵A、B構成。兩片長度均為L1且極化方向相反的壓電片對稱粘結在長度為L的懸臂梁根部上、下表面,壓電片串聯連接負載電阻R;壓電雙晶懸臂梁一端固定在基座左側,另一端連接磁鐵A(簡稱末端磁鐵);外部磁鐵B粘結在基座右側,并與末端磁鐵A水平相距為d。

1.1 非線性磁力模型

圖1所示壓電雙晶懸臂梁末端磁鐵A受到磁鐵B的非線性排斥力可以通過點磁荷模型計算得到,如圖2所示為磁鐵A與磁鐵B之間的幾何關系,圖中L為懸臂梁長度,a為磁鐵長度的一半,w(x,t)為懸臂梁x處的振動位移,α=?w(L,t)/?x為磁鐵A的偏轉角,這里假設磁鐵的偏轉角很小。同時,假設磁鐵A、B中磁場分布是均勻的,且磁極本身的幾何線度遠小于磁鐵間的距離,則磁鐵A、B可以看成是點磁荷,采用磁偶極子模型計算得到磁鐵B和磁鐵A之間的勢能為[9-10]

(1)

式中:μ0=4π×10-7H/m為真空磁導率,為向量梯度,mA為偶極子A的磁矩,mB為偶極子B的磁矩,且mB=MBVB,MB為磁鐵B的磁化強度,VB為磁鐵B的體積,rBA為磁鐵B到磁鐵A的方向向量。

圖2 非線性磁力模型

由圖2可以得到:

rBA=-d·i+w(L,t)·j+0.k
mA=MAVAcosα·i+MAVAsinα·j+0·k
mB=-MBVB·i+0·j+0·k

(2)

又α=arctan(w′(L,t)),則

(3)

將式(2)、式(3)代入式(1)得非線性勢能:

(4)

則磁鐵A受到磁鐵B作用的非線性磁力為:

(5)

1.2 壓電懸臂梁振型

由于壓電懸臂梁通常作低頻振動,其一階彎曲振動模態對系統動力學響應特性影響最大。因此,根據Raleigh-Ritz方法可將圖1所示壓電懸臂梁的振動位移w(x,t)寫成:

w(x,t)=φ1(x)q1(t)

(6)

式中:φ1(x)為一階彎曲模態振型,q1(t)為一階彎曲廣義模態坐標。

式中:

(7)

且振型函數可表示為:

(8)

為計算振型系數C1、C2、C3、C4、D1、D2、D3、D4,可列出以下邊界條件:

(9)

變截面梁兩部分的模態頻率關系為:

(10)

式中:ω1為一階模態振型角頻率。

1.3 壓電能量采集器動力學方程

壓電能量采集器的動力學方程采用拉格朗日方程建立,系統拉格朗日函數為[11-12]:

La(x,t)=Ts+Tp+TM+We-Us-Up-Um

(11)

式中:Ts、Tp分別為金屬基板和壓電片的動能,TM為懸臂梁末端磁鐵的動能,Us、Up分別為金屬基板和壓電片的勢能,We為壓電片的電能,Um為磁鐵對之間的勢能。它們的具體表達式如下:

(12)

(13)

(14)

(15)

(16)

(17)

式(11)中的非線性勢能Um可將式(4)在q1(t)=0處進行泰勒展開得到,即

Um=a0+0.5a1q1(t)2+0.25a2q1(t)4+o[q1(t)4]

(18)

式中:

將式(12)~式(18)代入式(11)所示的拉格朗日函數,并考慮系統廣義耗散能量得到系統動力學方程為

(19)

將式(19)轉換成狀態空間方程形式得:

(20)

1.4 諧波平衡法求解動力學方程

(21)

將式(21)代入式(19)中,令方程兩邊的sin(ωt)、cos(ωt)、常數項系數相同,并且忽略高次諧波項,則有:

C(k+sC2+1.5s(A2+B2))=0

(22)

B(k-ω2+0.75s(A2+B2))+2ξωA-θE+ΓA0=0

(23)

A(k-ω2+0.75s(A2+B2))-2ξωB-θD=0

(24)

E+θωRA+0.5ωCpRD=0

(25)

D-θωRB-0.5ωCpRE=0

(26)

r2=A2+B2

(27)

式中:r為一階模態振動位移的幅值。

由式(25)、式(26)得

(28)

(29)

將式(28)、式(29)代入式(23)、式(24)后平方相加得:

(30)

由式(22)解得,

(31)

將C代入式(30)后可以得到能量采集器的振動位移幅值r,并得到能量采集器輸出的電壓幅值和一個周期內的平均輸出功率分別為

(32)

(33)

2 采集器時-頻域動力學特性仿真

為了研究雙穩態壓電振動能量采集器動力學特性對系統參數敏感性,首先研究了系統靜態分岔特性以及勢能特性。圖3所示為能量采集系統靜態分岔分析結果,當末端磁鐵A和外部磁鐵B的中心距d大于分岔點間距d0=36.6 mm時,非線性磁力較小,系統只有一個穩定解,系統變成單穩態系統。當d小于分岔點間距d0時,能量采集器出現分岔現象,有兩個穩定的解和一個不穩定的解,且兩個穩定解的距離隨d的減小先增大后減小。圖3所示為能量采集器的勢能結果,當d<36.6 mm時,如圖4中d=35 mm、30 mm和25 mm時所示,采集器勢能曲線出現兩個對稱的勢阱和一個勢壘,表現出雙穩態特性,并且隨著d的減小,勢壘逐漸增大。當d大于36.6 mm,如圖4中d=40 mm所示,采集器的勢能曲線由雙勢阱退化成一個勢阱,表現出單穩態運動特性。

圖3 采集器靜態分叉特性

圖4 采集器勢能曲線

2.1 系統時域動力學特仿真分析

2.1.1 加速度幅值A0對動力學特性的影響

圖5 d=30 mm,f=5 Hz時采集器的時域動力學特性

根據靜態分叉特性和勢能特性的分析結果,取磁鐵間距d=30 mm,激振頻率f=5 Hz,研究激振加速度幅值A0對系統時域動力學特性的影響,結果如圖5所示。圖5(a)所示為不同激振加速度A0時能量采集器的分岔圖,由圖可知當A0小于8m/s2的時,系統作單穩態運動,并在某個勢能阱內作小幅周期振動,此時采集器的相圖為一個圓,且圓心位于某個勢能阱中心,如圖5(b)所示,采集器的振動位移幅值僅為1 mm;當激勵加速度幅值A0增大到8 m/s2~10.5 m/s2,采集器獲取的動能增加,系統開始越過勢壘,并由單周期小幅值振動進入三周期大幅值振動狀態,其相圖如圖5(c)所示,此時采集器的振動位移幅值和振動速度幅值都急劇增大;繼續增大激勵加速度幅值A0到10.5 m/s2~14 m/s2,采集器開始作大幅值的混沌運動,其相圖如圖5(d)所示;當激勵加速度幅值A0在14 m/s2~17 m/s2范圍內,采集器又進入三周期大幅值振動狀態,如圖5(e)所示。繼續增大激勵加速度A0,采集器最后進入大幅值雙穩態周期振動,如圖5(f)所示。總之,在A0<8 m/s2范圍內,采集器在某個勢能阱內作小幅值的單周期振動,在A0>8 m/s2時,采集器開始在兩個勢能阱之間來回作大幅值的振蕩,這有利于提高振動能量采集器的輸出性能,但對外部激勵條件要求較高,需要有足夠的激勵強度使采集器克服勢能阱之間的勢壘阻礙作用。此外,圖5所示運動現象說明在一定的激勵條件下混沌吸引子和周期吸引子可以互相轉化,從而是采集器表現出大幅值的振動特性,提高能量采集輸出能力。

2.1.2 磁鐵間距d對系統動力學特性的影響

取f=5 Hz,A0=8 m/s2,研究磁鐵間距對系統分叉特性和動力學特性的影響,如圖6所示。由圖6所示的分岔圖可知,當磁鐵間距d<23 mm時,由于磁力勢能太大,系統無法越過勢壘,采集器被限制某一個勢阱內作小幅運動,其相圖如圖6(b)所示;隨著d的不斷增大,磁力勢能不斷減小,系統開始越過勢壘在兩個勢能阱之間來回作大幅的多周期[d=30 mm,圖6(c)]和雙穩態[d=35 mm,圖6(d)]運動;繼續增大d,磁力勢能繼續減小,由于磁力作用減弱,采集器開始進入單穩態周期運動,懸臂梁末端振幅位移幅值也逐漸較小,如圖6(e)所示。

圖6 A0=8 m/s2,f=5 Hz時采集器的時域動力學特性

2.1.3 激振頻率f對系統動力學特性的影響

取磁鐵間距d=30 mm,激振加速度A0=8 m/s2,研究激振頻率f對系統分岔特性和動力學特性的影響,分別如圖7所示。由圖7(a)可知采集器系統頻率響應向特性非常復雜;在f=5 Hz以內,采集器作周期的大幅值振動,如圖7(b)所示為f=4 Hz時作單周期的大幅值振動,圖7(c)所示為f=5 Hz時作三周期的大幅值振動;在f=5 Hz~10.5 Hz范圍內,采集器作雙穩態大幅值振動,如圖7(d)所示為f=9 Hz時的雙穩態振動相圖;當f=10.5 Hz~21 Hz范圍內,系統進入混沌狀態,如圖7(e)和圖7(f)所示分別為f=11 Hz和f=21 Hz時的混沌振動相圖;當f大于21 Hz時,系統作小幅振動,如圖7(g)為f=22 Hz時的二周期小幅振動相圖;當f大于32 Hz時,采集器系統作小幅值單穩態周期運動,如圖7(h)所示為f=35 Hz時的小幅值單穩態相圖。從圖7所示結果可以看出,壓電能量采集器作大幅值振動的頻帶范圍為4 Hz~21 Hz。

圖7 A0=8 m/s2,d=30 mm時采集器的時域動力學特性

2.2 系統頻域動力學特仿真分析

本節利用諧波平衡法進行仿真,研究改變系統外部激勵(激振頻率、激振加速度)、磁鐵之間的距離、負載電阻時對系統輸出功率的影響。

2.2.1 采集器幅值-頻率特性

圖8所示時A0=10 m/s2時不同磁鐵間距d=[20 26 30 45]mm下的采集器末端振動幅值-頻率響應r-f曲線。圖中r小于0.5 mm為小幅阱內運動,r大于0.5 mm為大幅阱間運動。由圖8可知采集器的大幅運動和小幅運動可以共存,大幅運動振幅隨著f的增加而增加,且向右傾斜,表現出軟特性;小幅阱內運動主要集中在f=30 Hz~40 Hz,并且由于磁力引起的非線性剛度,頻率響應曲線向左傾斜,表現出硬特性。此外,隨著磁鐵間距d的增大,大幅阱間運動頻率范圍減小,振幅最大值增大,小幅阱間運動頻率范圍左移。

圖8 采集器幅值-頻率響應(實線為穩定解,虛線為不穩定解)

2.2.2 采集器輸出功率特性

圖9所示為f=5 Hz,R=7 MΩ,d=25 mm時激振加速度幅值與采集器輸出功率的關系。由圖9(a)可知當激振加速由小變大時功率緩慢增長,當激振加速到50 m/s2時,輸出功率突然變大,之后又緩慢增長,這是因為在剛開始振動幅值很小,懸臂梁振動位移和振動速度小,系統動能小于勢能,懸臂梁無法越過勢壘;當激振及速度不斷增大,其動能不斷增大,待動能大于勢能時才越過勢壘進入大幅振動狀態;反之,當激振加速度由大變小時,輸出功率先緩慢變小,隨著激振加速度的不斷減小,當A0大約為1 m/s2時,功率急劇下降,因為剛開始時激振加速度比較大,其動能也大,輕松越過勢壘作大幅阱間運動,因為其振動速度大,導致輸出功率也大,當A0<1 m/s2時,因動能太小無法越過勢壘,從而作小幅阱內運動。圖9(b)和圖9(c)分別給出了在大幅振蕩的情況下,不同磁鐵間距d=[20,25,30]mm時的輸出功率P隨激振加速A0的變化曲線。由圖9可知在阱間運動情況下輸出功率增長的速率隨d的增大而增大。

圖9 采集器輸出功率隨激振幅值的變化情況(實線為穩定解,虛線為不穩定解)

圖10是在大幅振蕩的情況下,激振頻率f=[4,8,12]Hz,A0=4 m/s2,d=30 mm時,輸出功率隨電阻的變化曲線,圖10可知,采集器輸出功率隨電阻的增大先增大后減小;當頻率增大,最優阻抗隨激振頻率增大而減小,最大輸出功率隨激振頻率增大而增大。

圖10 不同激勵頻率下輸出功率隨負載電阻的變化

圖11是在大幅振蕩的情況下,激振加速度A0=[4,10,15]m/s2,d=30 mm,f=8 Hz時輸出功率隨電阻的變化曲線,由圖11可知,最優阻抗并不會因為激振加速度改變而改變,但是最大輸出功率卻隨激振加速度增大而增大。

圖12是在大幅振蕩的情況下,磁鐵間距d=[10,20,30]mm,A0=4 m/s2,f=8 Hz時輸出功率隨電阻的變化曲線,由圖可知,在R=7 MΩ時采集器輸出功率達到最大,最優阻抗并不會因為磁鐵距離的改變而改變,但是最大輸出功率卻隨磁鐵之間距離的增大而減小。

圖11 不同激勵幅值下輸出功率隨負載電阻的變化

圖12 不同d下輸出功率隨負載電阻的變化

圖13 實驗測試系統

3 實驗驗證與結果分析

3.1 樣機與實驗系統

為了驗證雙穩態壓電振動能量采集系統參數對其時-頻域動力學特性的影響的正確性,搭建了實驗測試平臺,如圖13所示。實驗樣機的懸臂梁尺寸為70 mm×10 mm×0.15 mm,如圖14所示,懸臂梁末端用高強度膠粘結一永磁鐵,另一個外部磁鐵也通過高強度膠粘結在基座上;兩片壓電片(PZT-5A)通過高強度膠經高溫固化后粘結在金屬基體的上、下表面,PZT表明通過金屬銅電極串聯連接負載電阻R。實驗中,雙穩態壓電能量采集器樣機通過夾具安裝在激振器(JZK-5)的測試臺上,信號發生器(AFG3102C)產生的正弦信號經功率放大器(YE5871A)放大后,使激振器產生正弦振動;基礎加速度由加速度傳感器(YJ9A)測量得到,懸臂梁末端位移由激光位移傳感器(LK-G80)測量得到通過軟件在電腦上顯示,懸臂梁末端振動速度由激光測速儀(Polytech OFV-505)測得,采集輸出電壓由示波器測得和顯示。

圖14 懸臂梁樣機

3.2 實驗結果分析

根據仿真分析結果,選取磁鐵間距d=[26、30]mm,f=4.95 Hz,分別逐漸增大(正向掃描)和減小(反向掃描)激振加速度幅值,測量采集器的輸出功率有效值,實驗結果如圖15所示。由圖15可知,當d=26 mm時,正向增大激勵加速度幅值到51 m/s2,采集器由小幅值的單穩態運動跳轉到大幅值的雙穩態運動,這與圖9(a)所示仿真計算結果50 m/s2是吻合的,說明仿真計算結果是正確的。當反向減小激勵加速度幅值時,采集器輸出功率逐漸減小,當到加速度幅值減小到5 m/s2時,采集器由大幅值振動狀態跳轉到小幅值振動;同理,當d=30 mm時,采集器由小幅值振動跳轉到大幅值振動的激勵加速度幅值為19 m/s2,這與圖9(b)所示的仿真計算結果20 m/s2是吻合一致的。此外,由圖13所示的實驗結果還可以看出,增大磁鐵間距d,采集器由小幅值振動跳轉到大幅值振動所需要的激勵加速度幅值也減小,增大激振加速度采集器的突躍點左移,減小激振加速度,采集器的突躍點右移,這與圖9所示數值仿真的趨勢相似,仿真得以驗證。

圖15 采集器幅值-功率特性實驗結果

圖16是激振頻率f=4.95 Hz,激振加速度為A0=15 m/s2時,不同磁鐵之間的距離d=[26、28、30]mm下采集器輸出功率Prms隨負載電阻R的變化情況實驗結果。由圖可知,隨著電阻的增大,功率先增加到一個峰值,然后持續下降;功率隨著d的增大而減小,但是最優阻抗在d的改變下不變,為1.1 MΩ,趨勢與圖12所示數值仿真一致。

圖16 不同d下輸出功率隨負載電阻的變化實驗結果

圖17所示為磁體間距離d=28 mm,A0=10 m/s2,采集器在頻率f=[13 16]Hz時的輸出功率Prms隨負載電阻R的實驗變化結果。由圖可知,功率隨著激振頻率f的增大而增大,且最優阻抗隨著激振頻率f增大而減小,實驗結果趨勢與圖10仿真結果趨勢相似。

圖17 不同激勵頻率下輸出功率隨負載電阻變化的實驗結果

圖18是磁體間距d=28 mm,f=4.95 Hz,采集器不同激振加速度下的輸出功率Prms隨負載電阻R的變化實驗結果。由圖可知,功率隨著A0的增大而增大,且最優阻抗不變,實驗趨勢與圖11所示數值仿真結果相似。

圖18 不同激勵幅值下輸出功率隨負載電阻的變化實驗結果

為了進一步驗證仿真結果的正確性,取激振頻率為f=11 Hz,d=26 mm,逐漸增大激振加速度,測量雙穩態壓電發電系統不同激振加速度下的運動特性,結果如圖19所示。

圖19 d=26 mm,f=11 Hz時采集器的時域動力學特性實驗結果

圖19(a)為激振加速度幅值A0=10 m/s2時的相圖實驗結果,此時為小幅周期振動;圖19(b)為A0=20 m/s2時的相圖,此時采集器先作大幅值的雙穩態運動后進入混沌運動,繼續增大A0=30 m/s2,如圖19(c)所示,此時采集器則完全進入混沌運動狀態;繼續增加A0=45 m/s2,采集器則作大幅周期振動,如圖19(d)所示。上述實驗得到的運動特性變化趨勢與圖5所示仿真結果趨勢是相似的。

4 結論

建立了雙穩態壓電振動能量采集器的非線性磁力模型和機電耦合動力學方程,并通過諧波平衡法求解動力學方程的解,實驗結果驗證了仿真結果的正確性,得到以下結論:①改變磁鐵的間距會使系統由單穩態變為雙穩態;雙穩態系統在低頻率范圍內會出現大幅運動,且響應頻帶很寬。②系統存在大幅運動吸引子,一定情況下改變激振頻率以及激振加速度可以實現混沌吸引子和周期吸引子之間的轉化。③最優阻抗主要由激振頻率決定,激振加速度和磁鐵間距對其影響甚微。

猜你喜歡
振動
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
某調相機振動異常診斷分析與處理
大電機技術(2022年5期)2022-11-17 08:12:48
振動與頻率
天天愛科學(2020年6期)2020-09-10 07:22:44
This “Singing Highway”plays music
具非線性中立項的廣義Emden-Fowler微分方程的振動性
中立型Emden-Fowler微分方程的振動性
基于ANSYS的高速艇艉軸架軸系振動響應分析
船海工程(2015年4期)2016-01-05 15:53:26
主回路泵致聲振動分析
UF6振動激發態分子的振動-振動馳豫
計算物理(2014年2期)2014-03-11 17:01:44
帶有強迫項的高階差分方程解的振動性
主站蜘蛛池模板: 亚洲综合片| 国产尤物视频网址导航| 中文字幕人妻av一区二区| 精品一区二区三区水蜜桃| 在线观看无码a∨| 国产午夜人做人免费视频中文| 久久毛片网| 国产在线八区| 九九热这里只有国产精品| 美女扒开下面流白浆在线试听| 国产精品太粉嫩高中在线观看| 久夜色精品国产噜噜| 五月激激激综合网色播免费| 二级特黄绝大片免费视频大片| 久久婷婷国产综合尤物精品| 熟妇人妻无乱码中文字幕真矢织江| 广东一级毛片| 国产成人8x视频一区二区| 99尹人香蕉国产免费天天拍| 国产精品自在拍首页视频8| 亚洲天堂.com| 亚洲AV人人澡人人双人| 国产喷水视频| 亚洲国产精品日韩av专区| 国产免费看久久久| 国产精品无码作爱| 天天干伊人| 国产成人精品一区二区| 中文字幕人妻av一区二区| 无码国内精品人妻少妇蜜桃视频| 无码日韩人妻精品久久蜜桃| 毛片大全免费观看| 毛片免费网址| 人妻少妇乱子伦精品无码专区毛片| 国内熟女少妇一线天| 国产免费久久精品44| 免费a在线观看播放| 成人国产免费| 色AV色 综合网站| 欧美色视频日本| 国产亚洲精品va在线| 国产国产人成免费视频77777 | 美女无遮挡免费网站| 久久国产av麻豆| 中文字幕欧美日韩| 国产aⅴ无码专区亚洲av综合网| 亚洲成人精品| 欧美成人精品一区二区| 婷婷综合在线观看丁香| 秋霞一区二区三区| 91成人精品视频| 亚洲性影院| 色天天综合久久久久综合片| 国产综合色在线视频播放线视| 久久精品人人做人人爽电影蜜月| 国产精品无码AV片在线观看播放| 91国语视频| 99成人在线观看| 美女无遮挡被啪啪到高潮免费| 91福利在线看| 成人午夜免费视频| 天天躁夜夜躁狠狠躁图片| 色妞永久免费视频| 婷婷久久综合九色综合88| 东京热av无码电影一区二区| 久久性视频| 欧美综合在线观看| 青青草国产一区二区三区| 色吊丝av中文字幕| 亚洲国产精品久久久久秋霞影院| 曰韩免费无码AV一区二区| 国产精品乱偷免费视频| 中文字幕日韩丝袜一区| 超清无码熟妇人妻AV在线绿巨人| 亚洲成综合人影院在院播放| jizz国产在线| 亚洲伦理一区二区| 成人国产精品一级毛片天堂| 亚洲伦理一区二区| 国产在线精彩视频二区| 天堂网国产| hezyo加勒比一区二区三区|