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

余維-1非光滑分岔下的簇發振蕩及其機理?

2017-08-01 00:35:02張正娣劉楊張蘇珍畢勤勝
物理學報 2017年2期
關鍵詞:界面區域系統

張正娣 劉楊 張蘇珍 畢勤勝

(江蘇大學理學院,鎮江 212013)

余維-1非光滑分岔下的簇發振蕩及其機理?

張正娣 劉楊 張蘇珍 畢勤勝?

(江蘇大學理學院,鎮江 212013)

(2016年8月15日收到;2016年11月2日收到修改稿)

不同尺度耦合會導致一些特殊的振蕩行為,通常表現為大幅振蕩與微幅振蕩的組合,也即所謂的簇發振蕩.迄今為止,相關工作大都是圍繞光滑系統開展的,而非光滑系統中由于存在著各種形式的非常規分岔,從而可能會導致更為復雜的簇發振蕩模式.本文旨在揭示存在非光滑分岔時動力系統的不同尺度耦合效應.以典型的含兩個非光滑分界面的廣義蔡氏電路為例,通過引入周期變化的電流源以及一個用于控制的電容,選取適當的參數使得周期頻率與系統頻率之間存在量級差距,建立了含不同尺度的四維分段線性動力系統模型.基于快子系統在不同區域中的平衡點及其穩定性分析,以及系統軌跡穿越非光滑分界面時的分岔分析,得到了不同余維非光滑分岔的存在條件及其分岔行為.重點探討了余維-1非光滑分岔下的簇發振蕩的吸引子結構及其產生機理,揭示了非光滑分岔下系統復雜振蕩行為的本質.

非光滑電路系統,不同尺度,簇發振蕩,非常規分岔

1 引 言

在實際工程系統中往往存在大量非光滑因素,如碰撞[1]、干摩擦[2]、開關[3]、脈沖控制[4]等.根據所建立的數學模型中的非光滑度的不同,非光滑系統大致可以分為以下三類:1)非光滑連續系統,系統的向量場連續,而其Jacobi矩陣不連續,如蔡氏電路系統[5];2)Filippov系統,系統的向量場和其Jacobi矩陣均不連續,但其狀態空間連續,如干摩擦系統[6];3)非光滑脈沖系統,系統的向量場、Jacobi矩陣及其狀態空間均不連續,如碰撞系統[7].

與光滑系統相比,非光滑動力系統存在許多特殊的動力學現象,如豐富的運動形式、分岔行為及通向混沌路徑的多樣性[8].同時,由于系統的非光滑性,傳統的分析分岔及復雜性的方法對非光滑系統不再適用,需要探索一些新的方法和手段[9],因此,非光滑系統動力學的理論分析、數值計算和應用研究具有一定的挑戰性.

針對非光滑動力系統,國內外許多學者開展了一系列研究.例如,Shaw和Holmes[10]對低維情形下的沖擊振子和分段線性振子系統展開了一系列研究,發現該類系統存在著次諧響應、倍周期分岔等非線性特征;Nordmark等[11]通過建立Poincaré-Nordmark映射,分析了沖擊振子發生碰撞時的運動形式,發現其中存在著一種被稱作擦邊碰撞的現象,且在周期運動狀態下的振子,通過擦邊分岔可以直接進入混沌運動;Hu[12]針對分段光滑的系統,分析了向量場的非光滑性對Poincaré映射可微性的影響;Xu[13]探究了該類系統復雜動力學行為的產生機理;陸啟韶和金俐[14]研究分析了在剛性約束下的n維非線性系統的動力學行為,建立了這一類型的系統在剛性約束附近具有的局部映射關系,給出了該系統局部映射Jacobi矩陣的計算方法;姜海波等[15]對脈沖作用下Chen系統進行了非光滑分岔分析,運用了Floquet理論揭示了系統周期解的分岔機理.

由于缺乏有效的分析手段,相關研究大都是圍繞著平衡點或周期軌道展開的[16],并通過不同的數值方法建立其相應的分岔圖,例如強制法,僅能給出穩定的周期軌道,又如軌道跟蹤法,雖然能夠給出相關的不穩定軌道,但僅能展現常規分岔模式.為解釋非光滑分岔機理,以Leine為代表的學者引入微分包含理論來分析系統在非光滑分界面上的分岔特性[17],通過分析廣義Jacobi矩陣的特征值隨輔助參數的退化情況,揭示各種非光滑分岔的動力特性及其相應的產生機理.

迄今為止,相關工作大都是針對同一尺度下的非光滑系統的分岔展開的,基本上沒有涉及不同尺度耦合下的非光滑動力系統的研究,而事實上,不同尺度耦合系統,涉及到科學和工程技術的各個領域.例如,飛行器中存在著快速的旋轉運動與相對較慢的平動之間的耦合[18],繩系衛星中系繩的縱橫向振動與衛星姿態動力學之間的耦合[19],催化反應中存在著不同量級的反應速率[20],而幾乎所有的神經元模型幾乎都包含不同尺度之間的耦合[21].

不同尺度耦合系統的研究最早可以追溯到Poincaré研究行星軌道時提出的奇異攝動方程,但是直到諾貝爾獎獲得者Hodgkin和Huxley[22]建立了快慢兩尺度神經元放電模型(H-H模型),成功地再現了其中的簇發放電行為,不同尺度耦合系統的復雜性才引起了學術界的高度重視.由于傳統的非線性理論無法解決不同尺度之間的相互作用,早期相關工作大都局限在近似求解、數值仿真和實驗分析.直到2000年,Izhikevich[23]引入了Rinzel的快慢分析法才將相關研究上升到機理分析的層次.雖然近年來在不同尺度耦合領域取得了一定的進展,但絕大部分結果均是針對光滑系統取得的,而非光滑系統由于存在著非光滑分岔,不僅會導致一些特殊的簇發振蕩現象,也會使得相應的簇發振蕩結果更加復雜.同時,針對時域上不同尺度耦合的工作開展較多,而頻域上不同尺度耦合由于不存在明顯的快慢子系統,相關工作開展較少.

自Chua和Lin[24]構建了存在混沌現象的蔡氏電路以來,其動力特性引起了國內外學者的廣泛關注并開展了大量的研究[25,26].以原始的具有分段線性特性的蔡氏電路為基礎,通過引起其他元器件或改變其中非線性阻尼的伏安特性,得到了各種形式的蔡氏電路,給出了許多典型的非線性現象,如倍周期分岔的混沌道路、多渦卷混沌吸引子等等[27,28].但是,相關工作均是在同一尺度上開展的,同時,很少考慮其中的非光滑分岔特性.本文旨在揭示存在非光滑分岔時系統頻域上的不同尺度效應.以典型的含兩個非光滑分界面的廣義蔡氏電路為例,通過引入周期變化的電流源以及用于控制的電容,適當選擇參數,使得激勵頻率遠小于系統的固有頻率,建立了含頻域兩尺度耦合的四維分段線性非光滑動力系統.討論了非光滑分界面上的非常規分岔條件,分析了余維-1非光滑分岔下的簇發振蕩及其產生機理.

2 數學模型

具有分段線性特性的蔡氏電路不僅在實驗室容易搭建,而且存在著諸如概周期振蕩、混沌振蕩等復雜動力學行為,因此常被作為基礎模型來研究各種非線性現象.本文在典型的非光滑廣義蔡氏電路的基礎上,引入控制電容C1和周期變化的電流源iG,如圖1所示,其中包含兩個電感L1和L2,兩個電容C1和C2以及一個分段線性的非線性電阻RG,同時并聯一個周期變化的電流源iG,其相應的動力學模型可以表示為

圖1 電路原理圖Fig.1.Circuit schematic diagram.

3 快子系統分岔分析

取定參數Ω=0.01,α=0.001,而其他參數為常規量,也即考慮在弱控制下系統激勵頻率與系統固有頻率ω之間存在量級差距時其動力學行為. 顯然,在固有頻率對應的任一周期T∈[τ0,τ0+2π/ω]內, 外激勵項w將在WA=Asin(Ωτ0)和WB=Asin(Ωτ0+2πΩ/ω)之間變化,由于Ω?ω,WA≈WB,也即雖然外激勵項會在±A之間變化,但在任一固有頻率所對應的周期內,w變化非常緩慢,因此,整個外激勵項w可視為慢變參數.因此快子系統可以表示為

而慢子系統則為

在區域D0中,快子系統具有如下惟一的平衡點E0:

而在區域D±中,快子系統分別存在惟一的平衡點E±:

其中k1=γ+δb-qδb+qδa,k2=γ-δγ+δγbqδγb+qδγa,k3=-qδγb+qδγa+δγb. 系統中可能會產生兩種不同形式的非光滑分岔,一是在分界面處發生fold分岔,其產生條件為

另一是在分界面處發生Hopf分岔,其產生條件為

圖2 (q,δ)平面上的非光滑分岔集Fig.2.The non-smooth bifurcation sets on the(q,δ)plane.

取定參數β=1.2,γ=0.6,a=-3.0,b=0.6,A=3.0,圖2給出了(q,δ)平面上的非光滑分岔圖,顯然當δ<0.9446時,隨著q在[0,1]之間變化,廣義Jacobi矩陣存在零特征值,因此會產生余維-1非光滑fold分岔,而當δ>0.9446時,不僅會存在零特征值,也會存在一對純虛根,因此存在二次穿越現象,導致余維-2非光滑fold-Hopf分岔.

4 非光滑余維-1周期簇發振蕩

4.1 周期簇發振蕩運動過程

由圖2可知,當δ=0.9時,在分界面上僅會產生余維-1非光滑fold分岔,此時快子系統在三個不同區域中分別存在著三個不同的平衡點,計算可知,E0為不穩定鞍點,其相應的特征值為λ1=2.8476,λ2,3=-0.3738±0.6551I;E±為穩定焦點,其相應的特征值為λ1=-1.0651,λ2,3=-0.0375±0.5503I.圖3給出了δ=0.9時系統在(y,v)和(u,v)平面上的相圖.

圖3δ=0.9時的相圖 (a)(y,v)平面的相圖;(b)(u,v)平面的相圖Fig.3.Phase portraits forδ=0.9:(a)On the(y,v)plane;(b)on the(u,v)plane.

我們以(y,v)平面上的相圖為例來說明軌跡的運動情況.假設軌跡從非光滑分界面v=+1上的A1點出發,由于在D+域內存在惟一的穩定焦點E+,此時軌跡將大致按照E+的特性,逐漸振蕩趨于E+.在此過程中,由于外激勵是一個慢變過程,因此系統軌跡主要由快子系統決定.當軌跡逼近于E+時,由于E+是穩定的焦點,快子系統只會使得系統軌跡駐留在E+點,因此慢變的外激勵的影響將會體現出來.在外激勵的作用下,其相應的平衡點E+的位置會發生變化,大致沿著E+A3向分界面移動,從而導致系統軌跡緩慢地由A2點運動到A3點.當軌跡到達分界面上的A3點,由于在D0域內存在惟一的不穩定鞍點,在慢變外激勵驅動下,系統軌跡快速地穿越D0域,到達分界面v=-1上的A4點,并從A4點出發,逐漸振蕩趨于D-區域中惟一穩定的焦點E-.經過與D+區域內相對稱的過程,抵達分界面v=-1上的A6點,并快速穿越區域D0,到達起始點A1,完成一個周期的簇發振蕩.

4.2 周期簇發振蕩的結構

為進一步說明周期簇發振蕩的結構,圖4給出了其相應的時間歷程.

從圖3(a)中可以看出,在簇發振蕩對稱的半個周期內,大致包含三個部分:1)從分界面上A1點逐漸趨于E+的逼近過程;2)從E+附近向分界面上A3點的運動過程;3)從A3穿越區域D0到A4的過程.其中第二過程占據了主要部分.

第一過程對應于從A1到E+的逼近過程(見圖4(c)),振蕩頻率主要受E+所對應的共軛復根的虛部的影響,通過計算,其振蕩周期的理論值為TS=2π/0.5503=11.4177,這與圖4(c)中的數值模擬值TP=11.4272相符合,其振蕩幅值的變化主要受該對共軛復根的實部影響,其在v方向上的變化可以用曲線AP來描述,近似為VP=VP0exp[-0.0375(t-tP0)],其中VP0和tP0分別對應于在A1處圍繞E+振蕩在方向上的幅值和時間.從圖4(c)可以看出曲線AP大致與從A1點到E+的振蕩幅值符合,因此,按照振幅衰減到起始振幅的1%為標準,從A1點到E+時間的理論值可以計算為T1=122.80,這與圖4(c)中的數值模擬值TP1=141.45也較為符合.

第二過程主要對應于隨慢變量x的變化平衡點向分界面移動的過程,由于受外激勵的影響,該過程從E+振蕩趨于A3點,其振蕩周期為TE=2π/Ω,這與圖4(b)中的數值模擬結果一致,而振蕩幅值主要受(3)式中的第一式所對應的特征值λT=0.0001(-β-1)=-0.00022的控制,其在v方向上隨時間的變化大致可以用曲線AT來描述,其相應的表達式為VE=V0exp[-0.00022(t-t0)],其中v0和t0對應于軌跡在E+處的v值及所對應的時間,因此,可以從理論上計算從E+到A3點的時間為T2=ln(10.8761)/0.00022=10848.0.

第三過程對應于穿越區域D0的過程,由于時間非常短,可以忽略不計,通過計算可以得到簇發振蕩周期的理論解為TJ=2(10848.0+122.80)=21941.6,與數值模擬結果T0=23140符合良好(見圖4(a)).

4.3 簇發振蕩機理分析

為揭示該簇發振蕩的機理,在此采用包絡分析法. 其主要思想是,分別考察外激勵項w=Asin(Ωt)在w=±A極值處的平衡態及其分岔模式,并結合系統的相圖,得到不同簇發振蕩的產生機理.圖5給出了在(x,v)平面上系統簇發振蕩的相圖與快子系統的包絡平衡線的疊加圖.

圖5 (x,v)平面上簇發振蕩相圖與快子系統(x,v)平衡曲線的疊加圖Fig.5.The overlap diagram between the phase portrait and the equilibrium branches of the fast subsystem on the(x,v)plane.

圖5中B±分別對應于w=±A=±0.5時三個不同區域的快子系統的平衡線.假設系統軌跡從分界面v=+1上的A1點(見圖3(a))出發,軌跡受D+區域內的系統控制,逐漸趨于D+區域內快子系統惟一穩定焦點E+,由于A1點與E+之間存在較大的距離,使得逼近過程存在大幅振蕩特性,也即處于激發態并大致按照E+的特性,振幅逐漸減小,振蕩趨于E+,進入沉寂態雖然此時快子系統的軌跡將穩定于E+,但耦合系統的軌跡由于受第二個尺度也即慢變周期激勵的影響,將按照外激勵周期做小幅周期振蕩,也即開始激發態振蕩.在該過程中,由于整個外激勵值在±A之間周期變化,導致耦合系統在該階段周期振蕩的軌跡在分別取w=±A所得到的快子系統的兩條平衡線之間來回變化,并受兩條平衡線的制約,使得軌跡逐漸振蕩趨于分界面v=+1.

當軌跡到達分界面v=+1上的A3點時,受平衡線的制約,軌跡將穿越分界面,而由于在D0域內存在惟一的不穩定鞍點,所以在慢變外激勵驅動下,系統軌跡快速地穿越D0域,到達分界面v=-1上的A4點,形成沉寂態并從A4點出發,逐漸振蕩趨于D-區域中惟一穩定的焦點E-.同樣,由于A4點和E-之間存在著較大的距離,使得該逼近過程存在著較大幅值的振蕩,導致激發態由于慢變激勵項w按照外激勵頻率Ω在w=±A之間周期變化,導致軌跡也按照頻率Ω在w=±A所對應的兩條平衡線之間周期來回變化,直至軌跡到達分界面v=-1上的A6點,并快速穿越區域D0,到達起始點A1,完成一個周期的簇發振蕩.

需要指出的是,由于D0區域中不穩定鞍點所對應的實特征值較大,使得軌跡會快速穿越D0區域,分別趨于D±中的穩定焦點,同時,當軌跡抵達D+或D-區域中w=±A時快子系統的兩條平衡線之間時,振蕩軌跡會在這兩條平衡線之間來回周期振蕩,也即軌跡為w=±A時快子系統的兩條平衡線包絡.

5 結 論

對于周期激勵下含兩非光滑分界面的廣義蔡氏電路,選取適當的參數使得周期頻率與系統頻率之間存在量級差,建立了一個含頻域兩時間尺度的四維分段線性系統.通過對快子系統進行平衡點以及非光滑分岔分析可知,在具體參數下,當系統軌跡穿越非光滑分界面時,會產生余維-1非常規fold分岔,從而導致了簇發振蕩現象的產生.這與光滑系統中的簇發振蕩現象不同,光滑系統中軌跡是在系統同時存在的不同穩定平衡點之間跳躍,而在該非光滑系統中,軌跡是在不同區域中惟一的平衡點之間的跳躍,從而揭示了非光滑因素下非常規fold分岔對系統簇發振蕩行為的影響機理.本文僅探討了非光滑系統中的余維-1非常規fold分岔導致的簇發振蕩,對于系統的更復雜的簇發振蕩還有待進一步研究.

[1]Siefert A,Henkel F O 2014Nucl.Eng.Des.269 130

[2]Duan C,Singh R 2005J.Sound Vib.285 1223

[3]Zhusubaliyev Z H,Mosekilde E 2008Phys.Lett.A372 2237

[4]Jiang H B,Li T,Zeng X L,Zhang L P 2013Acta Phys.Sin.62 120508(in Chinese)[姜海波,李濤,曾小亮,張麗萍2013物理學報62 120508]

[5]Galvenetto U 2001J.Sound Vib.248 653

[6]Carmona V,Fernández-García S,Freire E 2012Physica D241 623

[7]Dercole F,Gragnani A,Rinaldi S 2007Theor.Popul.Biol.72 197

[8]Zhusubaliyev Z T,Mosekilde E 2008Physica D237 930

[9]Rene O,Baptista M S,Caldas I L 2003Physica D186 133

[10]Shaw S W,Holmes P A 1983J.Sound Vib.90 129

[11]Nordmark A,Dankowicz H,Champneys A 2009Int.J.Non-Linear Mech.44 1011

[12]Hu H Y 1995Chaos,Solitons Fractals5 2201

[13]Xu H D 2008Ph.D.Dissertation(Sichuan:Southwest Jiaotong University)(in Chinese)[徐慧東2008博士學位論文(四川:西南交通大學)]

[14]Lu Q S,Jin L 2005Acta Mech.Sol.Sin.26 132(in Chinese)[陸啟韶,金俐 2005固體力學學報 26 132]

[15]Jiang H B,Zhang L P,Chen Z Y,Bi Q S 2012Acta Phys.Sin.61 080505(in Chinese)[姜海波,張麗萍,陳章耀,畢勤勝2012物理學報61 080505]

[16]Stavrinides S G,Deliolanis N C 2008Chaos,Solitons Fractals36 1055

[17]Leine R I 2006Physica D223 121

[18]Jia Z,Leimkuhler B 2003Future Gener.Comp.Syst.19 415

[19]Leimkuhler B 2002Appl.Numer.Math.43 175

[20]Gyorgy L,Field R J 1992Nature355 808

[21]Duan L X,Lu Q S,Wang Q Y 2008Neurocomputing72 341

[22]Hodgkin A L,Huxley A F 1952J.Physiol.117 500

[23]Izhikevich E M 2000Int.J.Bifurcation Chaos6 1171

[24]Chua L O,Lin G N 1990IEEE Trans.Circuits Syst.37 885

[25]Mkaouar H,Boubaker O 2012Commun.Nonlinear Sci.Numer.Simul.17 1292

[26]Kahan S,Sicardi-Schifino A C 1999Physica A262 144

[27]Baptist M S,Caldas I L 1999Physica D132 325

[28]Stavrinides S G,Deliolanis N C,Miliou A N,Laopoulos T,Anagnostopoulos A N 2008Chaos,Solitons Fractals36 1055

PACS:05.45.-a,05.45.Pq DOI:10.7498/aps.66.020501

Bursting oscillations as well as the mechanism with codimension-1 non-smooth bifurcation?

Zhang Zheng-DiLiu Yang Zhang Su-Zhen Bi Qin-Sheng?

(Faculty of Science,Jiangsu University,Zhenjiang 212013,China)

15 August 2016;revised manuscript

2 November 2016)

The coupling of different scales in nonlinear systems may lead to some special dynamical phenomena,which always behaves in the combination between large-amplitude oscillations and small-amplitude oscillations,namely bursting oscillations.Up to now,most of therelevant reports have focused on the smooth dynamical systems.However,the coupling of different scales in non-smooth systems may lead to more complicated forms of bursting oscillations because of the existences of different types of non-conventional bifurcations in non-smooth systems.The main purpose of the paper is to explore the coupling effects of multiple scales in non-smooth dynamical systems with non-conventional bifurcations which may occur at the non-smooth boundaries.According to the typical generalized Chua’s electrical circuit which contains two non-smooth boundaries,we establish a four-dimensional piecewise-linear dynamical model with different scales in frequency domain.In the model,we introduce a periodically changed current source as well as a capacity for controlling.We select suitable parameter values such that an order gap exists between the exciting frequency and the natural frequency.The state space is divided into several regions in which different types of equilibrium points of the fast sub-system can be observed.By employing the generalized Clarke derivative,different forms of non-smooth bifurcations as well as the conditions are derived when the trajectory passes across the non-smooth boundaries.The case of codimension-1 non-conventional bifurcation is taken for example to investigate the effects of multiple scales on the dynamics of the system.Periodic bursting oscillations can be observed in which codimension-1 bifurcation causes the transitions between the quiescent states and the spiking states.The structure analysis of the attractor points out that the trajectory can be divided into three segments located in different regions.The theoretical period of the movement as well as the amplitudes of the spiking oscillations is derived accordingly,which agrees well with the numerical result.Based on the envelope analysis,the mechanism of the bursting oscillations is presented,which reveals the characteristics of the quiescent states and the repetitive spiking oscillations.Furthermore,unlike the fold bifurcations which may lead to jumping phenomena between two different equilibrium points of the system,the non-smooth fold bifurcation may cause the jumping phenomenon between two equilibrium points located in two regions divided by the non-smooth boundaries.When the trajectory of the system passes across the non-smooth boundaries,non-smooth fold bifurcations may cause the system to tend to different equilibrium points,corresponding to the transitions between quiescent states and spiking states,which may lead to the bursting oscillations.

non-smooth circuit system,multiple scales,bursting oscillations,non-conventional bifurcation

:05.45.-a,05.45.Pq

10.7498/aps.66.020501

?國家自然科學基金(批準號:11472115,11472116)和江蘇省青藍工程資助的課題.

?通信作者.E-mail:qbi@ujs.edu.cn.

*Project supported by the National Natural Science Foundation of China(Grant Nos.11472115,11472116)and the Qinglan Project of Jiangsu Province,China.

?Corresponding author.E-mail:qbi@ujs.edu.cn.

猜你喜歡
界面區域系統
Smartflower POP 一體式光伏系統
工業設計(2022年8期)2022-09-09 07:43:20
WJ-700無人機系統
ZC系列無人機遙感系統
北京測繪(2020年12期)2020-12-29 01:33:58
國企黨委前置研究的“四個界面”
當代陜西(2020年13期)2020-08-24 08:22:02
基于FANUC PICTURE的虛擬軸坐標顯示界面開發方法研究
連通與提升系統的最后一塊拼圖 Audiolab 傲立 M-DAC mini
人機交互界面發展趨勢研究
關于四色猜想
分區域
手機界面中圖形符號的發展趨向
新聞傳播(2015年11期)2015-07-18 11:15:04
主站蜘蛛池模板: 色悠久久久久久久综合网伊人| 亚洲中文字幕av无码区| 免费在线不卡视频| 久久国产V一级毛多内射| 久久久久国色AV免费观看性色| 亚洲综合日韩精品| 亚洲天堂网在线观看视频| 国产SUV精品一区二区| 国产精品亚洲专区一区| 国内精自线i品一区202| 欧美亚洲中文精品三区| 99这里只有精品在线| 日韩精品亚洲一区中文字幕| 成人国产精品一级毛片天堂| 99精品欧美一区| 麻豆精品在线视频| 永久免费av网站可以直接看的| 久久婷婷五月综合97色| 国产爽妇精品| 中文字幕在线看| 亚洲午夜福利精品无码不卡| 国产精品成人免费视频99| 亚洲伊人久久精品影院| 沈阳少妇高潮在线| 波多野结衣一二三| 国产成人精品在线| 亚洲AⅤ综合在线欧美一区| 久久国产亚洲偷自| 色视频国产| 亚洲激情99| 亚洲国产亚洲综合在线尤物| 国产乱人激情H在线观看| 亚洲成人福利网站| 国产乱人伦AV在线A| 免费午夜无码18禁无码影院| 少妇精品在线| 国产香蕉一区二区在线网站| 欧美亚洲网| 99爱在线| 丝袜美女被出水视频一区| 中文无码伦av中文字幕| 久久精品无码中文字幕| 欧美日韩在线观看一区二区三区| 久久久国产精品无码专区| 精品国产中文一级毛片在线看| 青青草国产免费国产| 免费xxxxx在线观看网站| 中文无码精品a∨在线观看| 亚洲国产系列| 91视频日本| 日本一区二区三区精品国产| 免费在线看黄网址| 青草视频久久| 99精品在线看| 国产一在线观看| 国产精品私拍在线爆乳| AV不卡在线永久免费观看| 一本大道香蕉高清久久| 婷婷丁香在线观看| 中文字幕欧美成人免费| 日本在线国产| 青青青视频91在线 | 毛片基地美国正在播放亚洲| 日本三级欧美三级| 国产精品视频猛进猛出| 国产精品hd在线播放| 国产日韩av在线播放| 国产精品美女免费视频大全| jijzzizz老师出水喷水喷出| 亚洲国产综合自在线另类| 伊人久久福利中文字幕| 97青草最新免费精品视频| 国产精彩视频在线观看| 日韩欧美中文亚洲高清在线| 日本草草视频在线观看| 亚洲IV视频免费在线光看| 国产欧美成人不卡视频| 2020亚洲精品无码| 国产情侣一区二区三区| 97se亚洲综合在线天天 | 亚洲成年人网| 亚洲日韩欧美在线观看|