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

非線性振蕩的周期態(tài)與混沌態(tài)

2016-04-15 09:52:16韋忠海黃蕙玲黃炳華
太原理工大學(xué)學(xué)報 2016年1期

韋忠海,黃蕙玲,黃炳華

(廣西大學(xué) 電氣工程學(xué)院, 南寧 530004)

?

非線性振蕩的周期態(tài)與混沌態(tài)

韋忠海,黃蕙玲,黃炳華

(廣西大學(xué) 電氣工程學(xué)院, 南寧 530004)

摘要:證明了具有多頻激勵源混頻構(gòu)成的一階微分電路也能誕生混沌。一階非自治動態(tài)電路,當(dāng)僅加入一個激勵源時,受迫振蕩呈現(xiàn)周期態(tài);當(dāng)加入多個不同頻激勵源時,誕生的受迫振蕩呈現(xiàn)混沌態(tài);以三個不同頻的激勵源為例,證明混沌振蕩的誕生,僅不過是振蕩周期的充分延長。可以用諧波平衡原理與功率平衡定理求微分方程的主諧波解。求解結(jié)果的正確性可用仿真相圖驗(yàn)證。

關(guān)鍵詞:混頻;混沌態(tài);周期態(tài);諧波解;空間曲線

通常可用非線性微分方程描寫系統(tǒng)的動態(tài)過程,多數(shù)的方程無法求出其解析,只能畫出其圖形解又稱相圖,方程的解包含非振蕩解與振蕩解兩種類型,前者指動態(tài)變量經(jīng)暫態(tài)后是單調(diào)變化的,或趨于無窮或趨于穩(wěn)定平衡點(diǎn),沒有持續(xù)振蕩特性;后者又可分為周期態(tài)與混沌態(tài),都是有界的。周期態(tài)表現(xiàn)為在仿真時間內(nèi),相軌線明顯的在一個單循環(huán)或多循環(huán)的閉合軌上不斷重復(fù),混沌態(tài)表現(xiàn)為在仿真時間內(nèi),相軌線在一個有限的區(qū)域內(nèi)不斷環(huán)繞,而無法看到軌線的閉合與重復(fù),延長仿真時間(要延長多少時間)會不會看到軌線的閉合,這是一個遙遙無期的答案。

混沌科學(xué)尚未誕生前,在兩變量構(gòu)成的坐標(biāo)相平面上,畫出相點(diǎn)連續(xù)運(yùn)動的相軌跡,顯示兩個動態(tài)變量相互間的非線性函數(shù)關(guān)系,做為方程的求解結(jié)果。這在20世紀(jì)50年代就已經(jīng)在非線性微分方程中廣泛應(yīng)用。例如范德堡振蕩的相圖是一個閉合的周期軌。

混沌相圖是求不出解析表達(dá)式微分方程的圖形解,是一條三維的空間曲線,這條曲線顯然要比范德堡振蕩形成的極限環(huán)復(fù)雜得多。它是三個動態(tài)變量相互關(guān)系的多元非線性函數(shù),這就是混沌函數(shù)。人類依靠這個圖形解推進(jìn)了非線性動態(tài)科學(xué)的發(fā)展,這就是混沌科學(xué)。用三維空間表示三個動態(tài)變量的相互函數(shù)關(guān)系,其有關(guān)的定義,概念,推理和結(jié)論可以推廣到多維的歐氏空間。

當(dāng)今用數(shù)值仿真畫出的相圖,都是用有理數(shù)近似的代替無理數(shù)進(jìn)行運(yùn)算與仿真,因而嚴(yán)格的說,畫出的相圖全部是周期態(tài)。當(dāng)非線性振蕩的周期較短,在仿真時間內(nèi)能完成一個周期,相圖呈現(xiàn)閉合軌,振蕩呈現(xiàn)周期態(tài);當(dāng)非線性振蕩的周期較長,在仿真時間內(nèi)不能完成一個周期,相圖呈現(xiàn)非周期的混沌態(tài)。由此可見混沌是振蕩周期的無限(或充分)延長,相圖的軌線在顯示的時段內(nèi),始終不重復(fù)。周期態(tài)與混沌態(tài)只不過是周期的長短。

1多諧波激勵引起的混沌態(tài)

例1電路如圖1所示。

圖1 混頻電路Fig.1 Mixing circuit

根據(jù)KCL與KVL可列式(1),電路參數(shù)如式(2),激勵源參數(shù)如式(3),式中w0表示公共基頻,是三個激勵頻率的最大公約數(shù)。電路中全部諧波成分均為w0的整數(shù)倍。

(1)

(2)

(3)如果電路的三個激勵源全部置零,并設(shè)電容的起始電壓u(0)=U0=1V,求零輸入解畫出相圖2-a。

如果設(shè)電容起始電壓u(0)=U0=0V,求零狀態(tài)解,當(dāng)E1m=10,E2m=0,E3m=0時,畫相圖2-b,當(dāng)E1m=10,E2m=8.3,E3m=0時畫相圖2-c;當(dāng)E1m=10,E2m=8.3,E3m=11時畫相圖2-d。

圖2 例1的相圖Fig.2 Phase diagram of example 1

(4)

由圖2-b可見,當(dāng)uF由單個正弦激勵源組成時誕生周期軌。由圖2-c與圖2-d可以看出,隨著加入激勵源個數(shù)的增加,相圖性狀逐步趨于復(fù)雜而形成混沌。

2公共基頻ω0高低和振蕩形態(tài)的關(guān)系

激勵值參數(shù)如式(3),設(shè)置較低的電壓,是為了便于用實(shí)驗(yàn)驗(yàn)證,制作一個實(shí)驗(yàn)裝置,即可以誕生周期態(tài)的信號,也可以誕生混沌態(tài)的信號,兩者的區(qū)別僅僅不過是周期的長短。因?yàn)橐话汶娮有盘柊l(fā)生器的電壓不宜過高。故三個激勵電壓按式(3)設(shè)置。只要三個激勵頻率公共基頻由式(5)變化,相圖形態(tài)隨著按圖3-a逐步變化為圖3-f。

(5)

圖3 公共基頻ω0逐步減低,振蕩由周期態(tài)演變成混沌態(tài)Fig.3 Public fundamehtal frequency steps down,oscillation changes from periodicstate to chaotic state is the abscissa,is the ordinate

由6個相圖的變化可見,隨公共基頻ω0的減低,振蕩周期延長,從開始是明顯閉軌的圖3-a,演變成混沌如圖3-f.只要ω0<10 000,即相圖均為充滿整個顯示范圍。如果三個頻率變?yōu)棣?=1 234 567;w2=2 211 111; ω3=3 322 111;ω0=1涂滿整個顯示范圍的顏色會更為加深,軌線環(huán)繞程度更加稠密。

3算法誤差嚴(yán)重影響相圖性狀

由以上論述可見,多頻激勵混頻誕生的振蕩,即可以是周期態(tài)也可以是混沌態(tài),取決于各激勵頻率公共基頻ω0的高低,當(dāng)基頻ω0很低,周期T很長時,在一定的仿真時間TsT的間隔內(nèi),相圖表現(xiàn)為閉合的周期軌。

如例1電路,隨激勵頻率的變化,可以有以下兩種振蕩性狀。

第一種性狀,如果取三個頻率ω1=106,ω2=3×106,ω3=7×106則ω0=106,穩(wěn)態(tài)振蕩周期T=6.28μs在Ts=62.8μs時段內(nèi)畫出的相圖如圖4,其相軌線很明顯形成閉合軌。說明當(dāng)Ts>T+Th時相圖是周期的,其中Th是進(jìn)入穩(wěn)態(tài)振蕩前暫態(tài)過程的時間。

第二種性狀,如果取三個激勵頻率ω1=1 234 567,ω2=3 322 111,ω3=7 654 321則公共基頻ω0=1是三個激勵頻率的最大公約數(shù),全體頻率成份(包括組合頻率ωeq)都是ω0的整數(shù)倍。因而這個周期函數(shù)的周期T=6.283 2s,在Ts=62.8μs的時段內(nèi)畫出的相圖如圖5顯然是非周期,說明當(dāng)Ts

當(dāng)三個頻率取ω1=(106+1),ω2=3×106,ω3=7×106則嚴(yán)格求其公共基頻ω0=1,按嚴(yán)格精確的算法,其周期T=6.283 2s,屬于第二種性狀,但畫出的相圖保持如圖4,這是由于算法誤差引起的。其中ω1=(106+1)和ω1=106在數(shù)值仿真中被視為相等。因而相圖畫出的結(jié)果,和第一種性狀完全相同。

圖4 當(dāng)T=6.28 μs時,Ts=62.8 μs的相圖       圖5 當(dāng)T=6.28 μs時,Ts=62.8 μs的相圖      Fig.4 Phase diagram of Ts=6.28 μs          Fig.5 Phase diagram of Ts=62.8 μs

4求三個主諧波解

4.1求三個主諧波的迭加解

電路圖1如果三個激勵源分別單獨(dú)求各自的諧波解,用MATH程序small-1a/2b/3c.nb可求得

當(dāng)僅有E1激勵源時,E1r=6, E1x=8, E2r=E2x=E3r=E3x=0,得

u1=3.427 67cos(1 010 000 t)+

2.819 44sin(1 010 000 t)

(6a)

當(dāng)僅有E2激勵源時,E1r=E1x=0, E2r=20/3, E2x=5, E3r=E3x=0,得

u2=2.327 38cos(2 550 000 t)+

3.965 89sin(2 550 000 t)

(6b)

當(dāng)僅有E3激勵源時,E1r=E1x=E2r=E2x=0, E3r=6.6,E3x=8.8,得

u3=-3.431 95cos(3 150 000 t)-

3.412 38sin(3 150 000 t)

(6c)

三個解迭加的結(jié)果可得

(6d)

4.2求主諧波的耦合解.

首先求非線性支路的耦合關(guān)系,采用相量法用下標(biāo)ν表示相量,下標(biāo)m表示幅值, r表示實(shí)部x表示虛部,文中的虛數(shù)單位采用電工學(xué)的符號j,角頻率記為ω=2πf。以變量u為例,u1表示u的ω1諧波分量的瞬時值,U1v表示u1的相量值,U1m表示幅值;又如非線性支路電流用in表示瞬時值,in1表示in的ω1諧波分量,In1v表示in1的相量,In1m表示幅值,In1r表示實(shí)功成份,In1x表示虛功成份,各值關(guān)系如式(7)。

u1=U1rsinω1t+U1xcosω1t,

(7)

電路圖1包含三個不同頻的諧波源,可以劃分成三個分部網(wǎng)絡(luò)。但三個諧波分部網(wǎng)絡(luò)相互之間有非線性耦合,表現(xiàn)在ω1諧波分部網(wǎng)絡(luò)的gn1eq包含有另兩個諧波U2m,U3m的貢獻(xiàn),余類推。也就是說三個分部網(wǎng)絡(luò)不能獨(dú)立求解后迭加;而必須首先求耦合關(guān)系而后共同聯(lián)合求耦合解。

(8)

u=u1+u2+u3=U1rsinω1t+U1xcosω1t+

(9)

(10)

壓控非線性電導(dǎo)gn的伏安關(guān)系如式(8),電流in包含有眾多諧波,其中imain代表主諧波,inon代表非主諧波。設(shè)gn的控制電壓包含u1+u2+u3三個主諧波,其瞬時值如式(9),相量值和幅值如式(10)。以式(9)代入式(8),在本例證的計(jì)算中暫時忽略非主諧波inon的影響,設(shè)in≈imain,只考慮主諧波的計(jì)算。MATH程序coupl.nb解出三個諧波電流(in1,in2,in3)與三個諧波電壓(u1,u2,u3)的關(guān)系如式(11)與式(12)。

(11)

(12)

三個主諧波的電流相量與同諧波電壓相量同相位,其中比值gn1eq,gn2eq,gn3eq,稱等效基波電導(dǎo)是實(shí)常數(shù)如式(11)和式(12),它體現(xiàn)三個分部網(wǎng)絡(luò)是相互耦合關(guān)聯(lián),不是分成相互無關(guān)的三個孤立部分。要注意式(11),式(12) 與k2無關(guān),因?yàn)榕即雾?xiàng)的展開對基波沒有貢獻(xiàn)。更要注意不同頻率的諧波成份的相量不能相加。

(13)

(14)

三個激勵源的相量值如式(13).用相量法列出三個分部網(wǎng)絡(luò)節(jié)點(diǎn)電壓U的三個電流平衡方程如式(14).其求解變量是相量U1v,U2v,U3v包含有6個未知實(shí)數(shù),三個方程的實(shí)部與虛部必須各自等于零,因而可以建立6個實(shí)數(shù)平衡等式,解出6個未知實(shí)數(shù),用MATH程序small-123.nb求其結(jié)果為:

(15)

4.3比較迭加解與耦合解的差別

比較式(6)的迭加解與式(15)的耦合解,兩者的差別很微小,其結(jié)果幾乎一致,有下列的原因。

因?yàn)橹磺笕齻€主諧波,求解中沒有包括主諧波的倍頻與組合頻率成份。三個主諧波相互之間的耦合量比較少。從這個觀點(diǎn)考慮,這是一個弱非線性系統(tǒng):況且偶次項(xiàng)展開對主諧波沒有貢獻(xiàn),在求主諧波時偶次項(xiàng)沒有發(fā)揮非線性的作用可忽略不計(jì)。三次項(xiàng)展開對主諧波的貢獻(xiàn),僅僅體現(xiàn)在非線性因子k3,而k3=1/106=k1/20;只有線性項(xiàng)因子k1的1/20。

三個激勵值較小,如果激勵量充分小,可視為在微變范圍內(nèi)電路是線性的,因而求迭加解會得到正確的結(jié)果。激勵值參數(shù)設(shè)置較低的電壓,是為了便于用實(shí)驗(yàn)驗(yàn)證非線性振蕩的兩種振蕩形態(tài),可以用同一實(shí)驗(yàn)裝置實(shí)現(xiàn)。因?yàn)橐话汶娮有盘柊l(fā)生器電壓不宜過高。故式(3)設(shè)置較低電壓,但在這樣低電壓下,當(dāng)三個激勵頻率公共基頻的變化如式(5)變化時,相圖3變化可印證有兩種振蕩形態(tài)。

再次還要看非線性元件gn在電路網(wǎng)絡(luò)中的影響力。如果視gn和電容c以左部分是一個有源網(wǎng)絡(luò),gn和c是這個有源網(wǎng)絡(luò)的端口負(fù)載,則當(dāng)負(fù)載的導(dǎo)納值增大時,負(fù)載電流增大端口電壓即我們的求解變量u必然降低,這種影響力還要看并聯(lián)電導(dǎo)gn在整個導(dǎo)納中的比重。當(dāng)控制電壓變量U1m,U2m,U3m增大時,gn的等效電導(dǎo)值gn1eq,gn2eq,gn3eq增大,僅受三次項(xiàng)系數(shù)k3的影響,按本例提供的電路參數(shù),由于k3較小,其他諧波成份通過耦合,促使待求諧波成份的基波電導(dǎo)值增大,從而引起求解變量u的降低量很微小,幾乎可以忽略。例如ω2,ω3諧波成份通過耦合,促使ω1諧波成份的等效基波電導(dǎo)值gn1eq增大,從而引起求解變量u1的降低量很微小。

4.4改變電路參數(shù)

為了說明耦合解與迭加解有明顯差別,可以改變電路參數(shù),如果電路參數(shù)改成k3=8/105=k1,三個激勵值各放大10倍,則由于考慮其他諧波對本諧波耦合量的增加,促使本諧波成份的基波電導(dǎo)值增大,從而引起求解變量u的降低量就非常明顯。

令E1r=60, E1x=80, E2r=200/3, E2x=50;E3r=66, E3x=88,程序large-1a/2b/3c.nb求得:

當(dāng)僅有E1激勵源時,E1r=60, E1x=80, E2r=0, E2x=0, E3r=0, E3x=0,得

u1=22.44cos[1 010 000 t]+

17.87sin[1 010 000 t],

(16a)

當(dāng)僅有E2激勵源時,E1r=0, E1x=0, E2r=200/3, E2x=50, E3r=0, E3x=0,得

u2=15.87cos[2 550 000 t]+

24.67sin[2 550 000 t]

(16b)

當(dāng)僅有E3激勵源時,E1r=0, E1x=0, E2r=0, E2x=0, E3r=66, E3x=88,得

u3=-22.54cos[3 150 000 t]-

20.18sin[3 150 000 t],

(16c)

(16d)

式(16d)是沒有考慮耦合影響,用迭加法求出的電壓總量(假設(shè)各諧波的電壓峰值有重迭的機(jī)會)。如果考慮耦合影響用程序large-123.nb求聯(lián)合解。

當(dāng)E1r=60, E1x=80, E2r=200/3, E2x=50, E3r=66, E3x=88,得

(17)

式(17)是考慮耦合影響,聯(lián)解求出的電壓總值(假設(shè)各諧波的電壓峰值有重迭的機(jī)會)。比較式(16)與(17)可見,考慮耦合影響后,輸出電壓三個主諧波的總 值有明顯減小。但電子實(shí)驗(yàn)室的信號發(fā)生器一般達(dá)不到這樣高的電壓,不便于實(shí)驗(yàn)驗(yàn)證。

5結(jié)論

(18)

(19)

含有激勵源一階微分方程的普遍形式表為式(18),式中uF是激勵源組成的激勵成份,可表為uF=uF(t)等效于一個非正弦的周期源;f(u)是僅含u的非線性函數(shù),u可表為u=u(t)。用式(18)對照式(1)可得

uF=e1g1/c+e2g2/c-e3g1/c.

(20)

圖3仿真解橫座標(biāo)代表求解變量,其最大值和主諧波解的式(6)或式(15)是大體近似的。

由圖3可見,uF由三個激勵源組成,當(dāng)公共基頻較高時誕生多循環(huán)的周期軌如圖3a或圖4。當(dāng)公共基頻較低時誕生混沌如圖3f或圖5。由此可見,只要混頻激勵組成的uF有足夠多的諧波源,并且公共基頻足夠低, 就能很普遍的在一階微分電路誕生混沌,那么二階以上的微分電路產(chǎn)生混沌的可能就更加廣泛了。 因而說,不久的將來,千千萬萬的混沌函數(shù)會相繼涌現(xiàn)出來。當(dāng)今出現(xiàn)的一些混沌系統(tǒng)是千萬個混沌函數(shù)先行的例子。

參考文獻(xiàn):

[1]黃炳華,劉慧杰,梁永清.一階微分電路構(gòu)成的混沌[J].ModernPhysics,2014,4(5):86-99.

[2]黃炳華,梁永清,韋忠海.多激勵源混頻構(gòu)成的混沌[J].ModernPhysics,2014,4(6):147-159.

[3]HUANGBinghua,LIGuangming,LIUHuijie.Powerbalancetheoremoffrequencydomainanditsapplication[J].JournalofModernPhysics,2014,5:1097-1108.

[4]HUANGBinghua,HEXiaoyang.Powerbalanceofmulti-harmoniccomponentsinnonlinearnetwork[J].JournalofModernPhysics,2014,5:1321-1331.

[5]HUANGBinghua,WEIYafen,HUANGYing.Describingchaosofcontinuoustimesystemusingboundedspacecurve[J].JournalofModernPhysics,2014,5:1489-1501.

[6]梁永清,黃炳華.非線性電路頻域的功率平衡[J].太原理工大學(xué)學(xué)報,2014,45(3):328-333.

[7]馮久超,李廣明.功率平衡理論在研究非線性電路與混沌中的進(jìn)展[J].華南理工大學(xué)學(xué)報,2012,40(11):13-18.

[8]HUANGBinghua,HUANGXinmin,LIHui.Maincomponentsofharmonicsolutionsofnonlinearoscillations[J].ProcediaEngineeringVolume,2011,16:325-332.

[9]HUANGBinghua,YANGGuangsong,WEIYafen,etal.Harmonicanalysismethodbasedonpowerbalance[J].AppliedMechanicsandMaterials,2013:325-326.

[10]Mickens R E.Comments on the method of harmonic balance[J].Journal of Sound Vibration,1984,94(3):456-460.

[11]Atadan A S,Huseyin K.An intrinsic method of harmonic analysis for non-linear oscillations[J].Journal of Sound Vib,1984,95(4):525-530.

[12]馮朝文,蔡理,康強(qiáng),等.一種新的三維自治混沌系統(tǒng)[J].物理學(xué)報,2011,60(3):030503.

[13]唐良端,李靜,樊冰,等.新三維混沌系統(tǒng)及其電路仿真[J].物理學(xué)報,2011,58(2):785-793.

[14]許喆,劉崇新,楊韜.一種新型混沌系統(tǒng)的分析及其電路實(shí)現(xiàn)[J].物理學(xué)報,2010,59(1):131-139.

[15]王杰智,陳增強(qiáng),袁著祉.一個新混沌系統(tǒng)及其性質(zhì)研究[J].物理學(xué)報,2006,55(8):3956-3963.

[16]劉凌,蘇燕辰,劉崇新.一個新混沌系統(tǒng)及其電路仿真實(shí)驗(yàn)[J].物理學(xué)報,2006,55(8): 3933-3937.

[17]YU Simin,QIU Shuisheng,LIN Qinghun.New results of study on generating multiple-scroll chaotic attractors[J].Science in China,2003,46(2):104-105.

[18]YU S M,TANG W K S,LU J,et al.Generating 2n-wing attractors from Lorenz-like systems[J].International Journal of Circuit heolications,2010,38(3):243-258.

(編輯:劉笑達(dá))

Periodic and Chaotic Status of Nonlinear Oscillations

WEI Zhonghai,HUANG Huiling,HUANG Binghua

(ElectricalEngineeringSchoolofGuangxiUniversity,NanningGuangxi530004,China)

Abstract:This paper proves that the first order differential circuit which is constituted by mixing of multi-excitation source with different frequency can also produce chaos. Regarding the first order non-autonomous circuits,the forced oscillation shows periodic status when only one excitation source is added.The forced oscillations display non-periodic status when multi-excitation source are added. Taking three excitation source as example, this paper proves that the births of chaotic oscillations are from sufficient extension of oscillation cycle. The main harmonic components in the differential equations can be solved by using the principle of the harmonics balance and power balance theorem.The correctness of solving results can be verified by phase portrait plotted by simulation.

Key words:mixing;chaotic status;periodic status;harmonic solution;space curve

中圖分類號:TN711.4;TN752

文獻(xiàn)標(biāo)識碼:A

DOI:10.16355/j.cnki.issn1007-9432tyut.2016.01.010

作者簡介:韋忠海(1964-),男,廣西平果人,實(shí)驗(yàn)師,主要從事電力電子及運(yùn)動控制研究,(E-mail)wshqfy@gxu.edu.cn通訊作者:黃炳華,教授,主要從事電路與系統(tǒng)研究,(E-mail)gxuhbh@163.com

基金項(xiàng)目:國家自然科學(xué)基金資助項(xiàng)目:基于非線性微分方程基礎(chǔ)上的功率平衡的研究(60662001)

收稿日期:2015-05-26

文章編號:1007-9432(2016)01-0046-07

主站蜘蛛池模板: 一级黄色网站在线免费看| 国内嫩模私拍精品视频| 欧美久久网| 欧美日韩中文国产| 九色综合视频网| 天天综合天天综合| 在线视频一区二区三区不卡| 亚洲综合18p| 亚洲精品少妇熟女| 波多野结衣亚洲一区| 久青草免费视频| 欧洲一区二区三区无码| 国产精品19p| 蜜臀AV在线播放| 国产微拍一区| 91精品啪在线观看国产91| 天天色综合4| 丰满少妇αⅴ无码区| 污视频日本| 亚洲欧洲自拍拍偷午夜色无码| 国产黄色片在线看| 久久黄色视频影| 国产亚洲欧美日韩在线观看一区二区| 国产精品丝袜视频| 国产精品永久在线| 欧美影院久久| 69国产精品视频免费| 欧美一区二区精品久久久| 91最新精品视频发布页| 日韩一级二级三级| 尤物特级无码毛片免费| 国产美女精品在线| 午夜无码一区二区三区在线app| 高清色本在线www| 免费又爽又刺激高潮网址| 亚洲视频免费播放| 高清欧美性猛交XXXX黑人猛交| 露脸国产精品自产在线播| www.91在线播放| 99视频精品在线观看| 看国产毛片| 国产一级二级在线观看| 色婷婷亚洲十月十月色天| 无码中文字幕精品推荐| 91青草视频| 国产91丝袜| 国产99视频免费精品是看6| 2021精品国产自在现线看| 国产超碰在线观看| 亚洲中文字幕国产av| 亚洲综合专区| 欧美亚洲另类在线观看| 欧美乱妇高清无乱码免费| 欧美一区二区人人喊爽| 日韩一级毛一欧美一国产| 97国产成人无码精品久久久| 欧美a网站| 538精品在线观看| 婷婷亚洲最大| 国产在线观看成人91| 国产亚洲精品97AA片在线播放| 97狠狠操| 中文字幕永久在线观看| 久久久久人妻一区精品色奶水 | 2022精品国偷自产免费观看| 国产一级小视频| 中国毛片网| 精品自拍视频在线观看| 国产一级一级毛片永久| 欧美成人一区午夜福利在线| 亚洲永久视频| 久久不卡国产精品无码| 国产精欧美一区二区三区| 最新亚洲人成无码网站欣赏网 | 国产成人精品视频一区视频二区| 欧美一级色视频| 久久国产精品夜色| 精品国产亚洲人成在线| 韩日无码在线不卡| 欧美色丁香| 天天婬欲婬香婬色婬视频播放| 少妇精品网站|