亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        一類對(duì)稱雙勢(shì)阱Duffing系統(tǒng)周期解的衍生與演化

        2017-12-19 07:03:58肖世富陳紅永牛紅攀
        關(guān)鍵詞:勢(shì)阱振子非對(duì)稱

        肖世富,陳紅永,牛紅攀

        (中國工程物理研究院總體工程研究所,四川 綿陽 621999)

        一類對(duì)稱雙勢(shì)阱Duffing系統(tǒng)周期解的衍生與演化

        肖世富,陳紅永,牛紅攀

        (中國工程物理研究院總體工程研究所,四川 綿陽 621999)

        針對(duì)自旋單自由度三次強(qiáng)化彈簧—質(zhì)量振子系統(tǒng)建立的對(duì)稱雙勢(shì)阱Duffing方程,通過把數(shù)值計(jì)算與諧波平衡半解析分析相結(jié)合,系統(tǒng)分析了該類Duffing系統(tǒng)在諧波強(qiáng)迫激勵(lì)下周期解隨激勵(lì)頻率變化的衍生與演化現(xiàn)象,獲得了不同頻段諧波強(qiáng)迫激勵(lì)下系統(tǒng)周期解的類型、周期解的衍生模式與演化規(guī)律。分析結(jié)果表明,該類Duffing方程存在平衡點(diǎn)臨域局部周期解以及鞍結(jié)點(diǎn)分岔衍生的對(duì)稱、反對(duì)稱與非對(duì)稱等多種全局周期解;局部或無對(duì)稱性的全局周期解直接通過倍周期分岔通向混沌運(yùn)動(dòng);全局對(duì)稱周期一解和反對(duì)稱周期三次諧波解首先各自發(fā)生對(duì)稱和反對(duì)稱破缺,再通過倍周期分岔演化為混沌。研究有助于深化對(duì)Duffing方程非線性現(xiàn)象及其演化規(guī)律的認(rèn)識(shí)。

        Duffing方程;諧波平衡法;周期解;鞍結(jié)點(diǎn)分岔;對(duì)稱性破缺;混沌

        0 引言

        非線性因素為動(dòng)力系統(tǒng)帶來了一系列線性系統(tǒng)所沒有的復(fù)雜動(dòng)力學(xué)現(xiàn)象,包括多解、跳越與滯后、極限環(huán)、固有頻率漂移、次諧波共振、超諧波共振、組合共振、內(nèi)共振、倍周期分岔、概周期分岔、混沌運(yùn)動(dòng)等等。Duffing方程是一類典型的非線性動(dòng)力系統(tǒng),是Duffing于1918年描述機(jī)械系統(tǒng)中彈簧硬化效應(yīng)建立的一類方程,具有如下形式

        Duffing方程蘊(yùn)含著豐富的復(fù)雜非線性現(xiàn)象,包括周期振動(dòng)、次諧波共振、局部與全局倍周期分岔、對(duì)稱與反對(duì)稱破缺、概周期運(yùn)動(dòng)、混沌以及多種周期或混沌運(yùn)動(dòng)共存等。因此,在非線性振動(dòng)理論中研究Duffing方程具有解剖一只麻雀的作用[1]。目前,Duffing方程與van der Pol方程是非線性振動(dòng)問題研究中最常用的例子[1-3],在力學(xué)、電子學(xué)等領(lǐng)域具有廣泛的研究和應(yīng)用價(jià)值。

        劉延柱等[2]和聞幫椿等[3]系統(tǒng)地總結(jié)了求解弱非線性Duffing方程近似解析解的各種方法。Marina和Herisanu[4]給出了求解強(qiáng)非線性Duffing系統(tǒng)近似解析解的積分迭代方法。陳艷鋒等[5]對(duì)無阻尼Duffing方程自由振動(dòng)頻率的求解方法進(jìn)行了探討,應(yīng)用Gauss*Chebyshev求積公式計(jì)算了Duffing方程的自由振動(dòng)頻率,得到了精確解析解表達(dá)式。王坤等[6]獲得了一類Duffing振子系統(tǒng)具有唯一周期解的必要條件,給出了一定條件下Duffing振子系統(tǒng)精確周期信號(hào)的獲取方法。馮少東與陳立群[7]利用同倫分析方法求解了恢復(fù)力為有理函數(shù)的Duffing簡(jiǎn)諧振子的響應(yīng)和頻率的近似周期解,與精確解吻合?;煦缡欠蔷€性動(dòng)力學(xué)系統(tǒng)研究的重點(diǎn)之一,倍周期分岔被公認(rèn)為是通向混沌的道路之一,而在對(duì)稱非線性系統(tǒng)中,對(duì)稱破缺又是倍周期分岔的前兆[8-9],驅(qū)動(dòng)振子逃離勢(shì)阱[10-11]。諧波強(qiáng)迫激勵(lì)對(duì)稱Duffing方程存在倍周期分岔演化混沌現(xiàn)象,其對(duì)稱性破缺研究也獲得了學(xué)術(shù)界的關(guān)注。Novak和Frehlich研究表明Duffing系統(tǒng)是由對(duì)稱破缺到倍周期分岔進(jìn)入混沌的一個(gè)典型例子[8];畢勤勝等研究了單勢(shì)阱Duffing方程對(duì)稱破缺分岔轉(zhuǎn)遷集的解析表達(dá)式[9];張瑩等研究了諧波激勵(lì)雙勢(shì)阱Duffing方程混沌與周期吸引子的對(duì)稱破缺激變[12],以及隨機(jī)參數(shù)Duffing系統(tǒng)的對(duì)稱破缺分叉[13]。

        目前,對(duì)Duffing方程蘊(yùn)含的復(fù)雜性已經(jīng)有比較深入的認(rèn)識(shí),然而,也存在一些認(rèn)識(shí)不夠深入的問題。例如Duffing方程的周期解是如何衍生的?其反對(duì)稱和非對(duì)稱周期解是如何演化的?本文針對(duì)這些問題,以自旋單自由度三次強(qiáng)化彈簧—質(zhì)量振子系統(tǒng)建立的對(duì)稱雙勢(shì)阱Duffing方程為例,采用數(shù)值計(jì)算與諧波平衡半解析分析相結(jié)合的方法,系統(tǒng)分析該類Duffing系統(tǒng)在諧波強(qiáng)迫激勵(lì)下周期解隨激勵(lì)頻率變化的衍生與演化現(xiàn)象,進(jìn)一步提升對(duì)Duffing方程非線性復(fù)雜現(xiàn)象及其演化規(guī)律的認(rèn)識(shí)。

        1 雙勢(shì)阱Duffing方程數(shù)學(xué)模型與周期解

        圖1 轉(zhuǎn)動(dòng)單自由度彈簧-質(zhì)量振子系統(tǒng)示意圖Fig.1 Sketch of the rotational single degree of freedom spring-mass oscillator system

        F=Kx+εKx3

        (1)

        其中,K為彈簧剛度,ε為三次非線性強(qiáng)化系數(shù)。

        (2)

        (3)

        方程(3)對(duì)應(yīng)的勢(shì)能為

        (4)

        其靜態(tài)平衡解滿足

        (1-Ω2)x+εx3=0

        (5)

        當(dāng)Ω<1時(shí),方程(5)只有穩(wěn)定的平凡解零解;當(dāng)Ω≥1時(shí),平凡解發(fā)生叉式分岔而失穩(wěn),系統(tǒng)存在1對(duì)對(duì)稱穩(wěn)定的分岔解:

        (6)

        此情況下,系統(tǒng)的勢(shì)能函數(shù)與準(zhǔn)靜態(tài)分岔圖如圖2所示。

        圖2 對(duì)稱雙勢(shì)阱Duffing方程的勢(shì)能函數(shù)和準(zhǔn)靜態(tài)分岔圖Fig.2 Diagram of the potential energy and quasi-static bifurcation solution of the symmetric double-well Duffing equation

        當(dāng)Ω<1時(shí),方程(3)周期解的性質(zhì)研究已非常成熟,本文不再討論。對(duì)Ω>1于,方程(3)為負(fù)剛度強(qiáng)化系統(tǒng),是一類典型的強(qiáng)迫諧波激勵(lì)對(duì)稱雙勢(shì)阱Duffing方程,下面通過數(shù)值模擬分析不同頻段諧波激勵(lì)下該系統(tǒng)的周期解。為了分析的方便,確定參數(shù)ε=0.3,Ω=1.05,c=0.04,f=0.2,只考查不同激勵(lì)頻率ω下系統(tǒng)的周期響應(yīng)。

        對(duì)于方程(3),在取定的參數(shù)下,采用數(shù)值計(jì)算軟件對(duì)系統(tǒng)的周期解進(jìn)行分析,在不同激勵(lì)頻率和初值情況下,系統(tǒng)的周期解如圖3所示。

        圖3 不同頻率諧波激勵(lì)下對(duì)稱雙勢(shì)阱Duffing系統(tǒng)存在的部分周期解Fig.3 The partial periodic solutions for the symmetric doubkle-well Duffing systems under the harmonic excitation of different frequencies

        圖3表明,在不同頻段諧波激勵(lì)下,對(duì)稱雙勢(shì)阱Duffing方程存在多種周期解,包括靜態(tài)穩(wěn)定平衡點(diǎn)臨域的局部周期解(圖3a黑粗線所示)、全局反對(duì)稱次諧波周期解(圖3a虛線所示)、全局對(duì)稱周期解(圖3a幅值最大的周期解)、低頻非對(duì)稱周期解(圖3b所示)、低頻反對(duì)稱次諧波周期解(圖3c所示)、低頻反對(duì)稱周期解(圖3d、圖3e所示,表明不僅僅是次諧波是反對(duì)稱的),且多種周期解往往共存(圖3a所示,非線性動(dòng)力系統(tǒng)中,還可出現(xiàn)混沌與周期解共存的現(xiàn)象,本文不作分析)。圖3a共有5個(gè)周期解,包括兩個(gè)靜態(tài)平衡點(diǎn)附近的局部二周期解、兩個(gè)全局三周期次諧波解、以及1個(gè)全局大范圍的對(duì)稱一周期解,這些周期解存在各自的吸引域,將使Duffing系統(tǒng)重復(fù)正弦掃頻實(shí)驗(yàn)從本質(zhì)上出現(xiàn)不同的實(shí)驗(yàn)結(jié)果。

        Duffing系統(tǒng)的這種多周期解共存特性,可解釋即使系統(tǒng)本身沒有不確定性和再精細(xì)的控制,非線性系統(tǒng)振動(dòng)實(shí)驗(yàn)的實(shí)驗(yàn)結(jié)果也可能不重復(fù)現(xiàn)象的內(nèi)在機(jī)理。

        2 周期解的衍生

        強(qiáng)迫諧波激勵(lì)對(duì)稱雙勢(shì)阱Duffing方程的周期解分析表明,系統(tǒng)存在多種周期解,包括全局反對(duì)稱次諧波解、全局對(duì)稱周期解、全局非對(duì)稱周期解等,本文采用諧波平衡法[3]分析這些周期解衍生的分岔模式。

        設(shè)方程(3)的全局反對(duì)稱次諧波解具有以下形式

        x(τ)=A1cos(ωτ/3)+B1sin(ωτ/3)+A2cosωτ+B2sinωτ

        (7)

        將表達(dá)式(7)代入方程(3),應(yīng)用諧波平衡法,得到

        (8)

        (9.1)

        將代數(shù)方程(8)的后兩式平方求和,并由前兩式做變換消去交叉項(xiàng)有

        (9.2)

        聯(lián)立代數(shù)方程(9.1)、(9.2)進(jìn)行數(shù)值求解,結(jié)果表明次諧波由鞍結(jié)點(diǎn)分岔衍生,其臨界激勵(lì)頻率為ωc=1.33,其衍生的分岔模式如圖4所示。

        圖4 次諧波周期解的衍生模式Fig.4 Derivation pattern of sub-harmonic periodic solution

        設(shè)方程(3)的全局對(duì)稱周期解具有以下形式

        x(τ)=ρ1sin(ωτ+φ)

        (10)

        圖5 全局對(duì)稱周期解的衍生模式(ρ1~ω)Fig.5 Derivation pattern of global symmetric periodic solution(ρ1~ω)

        將表達(dá)式(10)代入方程(3),將由諧波平衡法得到的兩式求平方和有

        (11)

        代數(shù)方程(11)可解析求解,但其表達(dá)式非常復(fù)雜。仍然采用數(shù)值方法,分析表明全局周期解也由鞍結(jié)點(diǎn)分岔衍生而成,臨界激勵(lì)頻率為ωc=1.07,其衍生模式如圖5所示。

        設(shè)方程(3)全局非對(duì)稱周期解具有以下形式

        x(τ)=A0+ρ1sin(ωτ+φ1)+ρ2sin(2ωτ+φ2)

        (12)

        將表達(dá)式(12)代入方程(3),應(yīng)用諧波平衡法,并消去相位變量,有

        (13)

        數(shù)值計(jì)算(13)表明,非對(duì)稱周期解仍然通過鞍結(jié)點(diǎn)分岔衍生,其臨界激勵(lì)頻率為ωc=0.32,其衍生模式如圖6所示。

        同理可采用諧波平衡法分析低頻激勵(lì)下系統(tǒng)出現(xiàn)的反對(duì)稱周期解、反對(duì)稱次諧波周期解等周期解衍生的模式,分析的關(guān)健是諧波需要假設(shè)合理的形式。分析結(jié)果表明,本文所討論的強(qiáng)迫諧波激勵(lì)對(duì)稱雙勢(shì)阱Duffing方程初期演化的周期解都通過鞍結(jié)點(diǎn)分岔衍生而成。

        3 周期解的演化

        3.1 局部周期解的倍周期分岔

        雙勢(shì)阱Duffing方程靜平衡位置臨域的局部周期解通過倍周期分岔進(jìn)入混沌,其分岔圖如圖7a所示,相平面演化如圖7b~圖7f所示。

        圖6 全局非對(duì)稱周期解的衍生模式Fig.6 Derivation pattern of global asymmetric periodic solution

        圖7 局部周期解的分岔圖與倍周期分岔演化相圖Fig.7 Bifurcation diagram and Phase diagram of the evolution of local periodic solution

        需要說明的一點(diǎn)是,圖7b~圖7f周期解演化相圖僅給出了系統(tǒng)其中一個(gè)平衡點(diǎn)臨域局部周期解的相平面演化圖,其另一個(gè)平衡點(diǎn)臨域局部周期解同樣通過倍周期分岔進(jìn)入混沌,且與圖7所給的相應(yīng)結(jié)果關(guān)于速度軸對(duì)稱。

        3.2 全局反對(duì)稱次諧波解的對(duì)稱性破缺與三倍周期分岔

        數(shù)值分析表明,雙勢(shì)阱Duffing方程由鞍結(jié)點(diǎn)分岔衍生的全局次諧波解隨著激勵(lì)頻率的下降,相軌跡首先發(fā)生反對(duì)稱破缺,隨后以3的倍數(shù)發(fā)生分岔,演化到概周期直至混沌狀態(tài),其分岔圖如圖8a所示,相軌跡演化如圖8b~圖8f所示。然而,數(shù)值模擬中,次諧波共振周期解三倍周期分岔的概周期與混沌狀態(tài)比較難于捕捉。

        圖8 全局次諧波相軌跡的反對(duì)稱破缺及其分岔現(xiàn)象Fig.8 Anti-symmetry breaking and bifurcation of global subharmonic phase trajectory

        圖9的分析結(jié)果表明,雙勢(shì)阱Duffing方程全局次諧波周期解反對(duì)稱破缺后,解的周期數(shù)并沒改變,僅是解的個(gè)數(shù)加倍,有的文獻(xiàn)直接將對(duì)稱性破缺后的解認(rèn)為周期加倍是值得商榷的。

        4 結(jié)束語

        本文建立了自旋單自由度三次強(qiáng)化彈簧—質(zhì)量振子系統(tǒng)的對(duì)稱雙勢(shì)阱Duffing方程,采用數(shù)值計(jì)算與諧波平衡半解析分析相結(jié)合的方法,系統(tǒng)分析了該類Duffing系統(tǒng)在諧波強(qiáng)迫激勵(lì)下周期解隨激勵(lì)頻率變化的衍生與演化現(xiàn)象,獲得以下結(jié)果:

        1)Duffing方程在不同頻段諧波強(qiáng)迫激勵(lì)下,存在多種周期解:平衡點(diǎn)臨域局部周期解、全局對(duì)稱周期解、全局反對(duì)稱次諧波周期解、低頻全局非對(duì)稱周期解(成對(duì)出現(xiàn))、低頻全局反對(duì)稱周期解、低頻全局反對(duì)稱次諧波周期解等。

        2)Duffing系統(tǒng)周期解由兩類模式衍生:靜態(tài)平衡點(diǎn)臨域衍生的局部周期解;鞍結(jié)點(diǎn)分岔衍生全局對(duì)稱周期解、全局反對(duì)稱次諧波周期解、全局非對(duì)稱周期解等。

        3)Duffing系統(tǒng)平衡點(diǎn)臨域衍生的局部周期解直接通過倍周期分岔演化為混沌;鞍結(jié)點(diǎn)分岔衍生的全局對(duì)稱周期一解和反對(duì)稱周期三次諧波解首先發(fā)生對(duì)稱性破缺(對(duì)稱/反對(duì)稱破缺),再通過倍周期分岔演化為混沌、或通過三倍周期分岔演化為概周期運(yùn)動(dòng)。

        這些成果有助于深化對(duì)Duffing方程非線性現(xiàn)象及其演化規(guī)律的認(rèn)識(shí);同時(shí),也能為實(shí)際工程結(jié)構(gòu)振動(dòng)試驗(yàn)時(shí)出現(xiàn)的一些現(xiàn)象和問題提供機(jī)理性解釋。例如,非線性工程結(jié)構(gòu)振動(dòng)實(shí)驗(yàn)時(shí),經(jīng)常出現(xiàn)重復(fù)實(shí)驗(yàn)的實(shí)驗(yàn)結(jié)果不重復(fù)或偶爾出現(xiàn)控制超差等現(xiàn)象或問題,其原因除了研究對(duì)象或?qū)嶒?yàn)系統(tǒng)的不確定性以外,非線性系統(tǒng)的多解、突變等特性也是重要原因。

        [1]丁同仁. 常微分方程定性方法的應(yīng)用[M]. 北京: 高等教育出版社, 2004.

        [2]劉延柱, 陳立群. 非線性振動(dòng)[M]. 北京: 高等教育出版社, 2003.

        [3]聞邦椿, 李以農(nóng), 韓清凱. 非線性振動(dòng)理論中的解析方法及工程應(yīng)用[M].沈陽: 東北大學(xué)出版社, 2001.

        [4]Marinca V, Herisanu N. Periodic solutions of Duffing equation with strong nonlinearity[J]. Chaos, Soliton and Fractals, 2008, 37(1): 144-149.

        [5]陳艷鋒, 鄭建華, 吳新躍,等. 無阻尼Duffing方程高精度近似解研究[J]. 機(jī)械科學(xué)與技術(shù), 2008, 27(12): 1591-1594.

        ChenYanfeng, Zheng Jianhua, WuXinyue, et, al. On high-accuracy approximate solution of undamped Duffing equation[J]. Mechanical Science and Technology for Aerospace Engineering,2008, 27(12): 1591-1594.

        [6]王坤, 關(guān)新平, 丁喜峰,等. Duffing振子系統(tǒng)周期解的唯一性與精確周期信號(hào)的獲取方法[J]. 物理學(xué)報(bào), 2010, 59(10): 6859-6863.

        Wang Kun, Guan Xinping Ding Xifeng,et al. Acquisition method of precise periodic signal and uniqueness of periodic solutions of Duffing oscillator system[J]. Acta Physica Sinica, 2010, 59(10): 6859-6863.

        [7]馮少東, 陳立群. Duffing簡(jiǎn)諧振子同倫分析法求解[J]. 應(yīng)用數(shù)學(xué)和力學(xué), 2009, 30(9): 1015-1019.

        Feng Shaodong, Chen Liqun. Homotopy analysis approach to the Duffing-Harmonic oscillator[J]. Applied Mathematics and Mechanics, 2009, 30(9): 1015-1019.

        [8]Novak S, Frehlich R G. Transition to chaos in the duffing oscillator[J]. Physical Review, 1982, 26A: 3660-3663.

        [9]畢勤勝, 陳予恕. Duffing系統(tǒng)解的轉(zhuǎn)遷集的解析表達(dá)式[J]. 力學(xué)學(xué)報(bào), 1997, 29(5): 573-581.

        Bi Qinsheng, Chen Yushu. Analytical expression of transition boundaries of the solution of Duffing systems[J]. Acta Mechanica Sinica, 1997, 29(5): 573-581.

        [10] Bishop S R, Sofroniou A, Shi P. Symmetry-breaking in the response of the parametrically excited pendulum model[J]. Chaos, Solitons and Fractals, 2005, 25(2): 257-264.

        [11] Mann B P,Koplow M A. Symmetry breaking bifurcations of a parametrically excited pendulum[J]. Nonlinear Dynamics, 2006, 46(4): 427-437.

        [12] 張瑩, 雷佑銘, 方同. 混沌吸引子的對(duì)稱破缺激變[J]. 物理學(xué)報(bào), 2009, 58(6): 3879-3805.

        Zhang Ying, Lei You-Ming, Fang Tong. Symmetry breaking crisis of chaotic attractors[J]. Acta Physica Sinica, 2009, 58(6): 3879-3805.

        [13] Zhang Y, Du L, Yue X L, et al. Analysis of symmetry breaking bifurcation in duffing system with random parameter[J]. Computer Modeling in Engineering & Sciences, 2015, 106(1): 37-51.

        PeriodicSolutionsInitiationandEvolvingofaSymmetricalDouble-PotentialWellDuffingSystemwithForcedHarmonicExcitation

        XIAO Shifu, CHEN Hongyong, NIU Hongpan

        (Institute of Systems Engineering, CAEP, Mianyang 621999, China)

        For a Duffing equation with symmetrical double potential well, which is derived from a self-rotating mass-spring pendulum with three-order nonlinear stiffness, the periodic solutions type, initiation mode and evolving behavior adapting to the frequency of excitation are investigated meticulously by Harmonic Balance Method and numerical simulation. Some rewarding results are obtained. This Duffing system has local periodic solutions in the region near its equilibrium point and many global periodic solutions with symmetry, antisymmetry and dissymmetry. The global periodic solutions derive from saddle-node bifurcation. The local periodic solutions and the dissymmetrical global periodic solutions evolve directly into chaos passing through double-periodic bifurcation. The symmetrical global 1-periodic solutions evolve into chaos through passing the double-periodic bifurcation only after losting their symmetry by symmetry-breaking, and the antisymmetric global 3-periodic solutions evolve into quasi-periodic solutions through passing the triple-periodic bifurcation also after lost their symmetry by antisymmetry-breaking.The studies help to raise the understanding of the nonlinear phenomenon and the evolution of solutions of the Duffing equation.

        duffing equation; harmonic balance method; periodic solution; saddle-node bifurcation; symmetry-breaking; chaos

        1672-3813(2017)03-0103-08;

        10.13306/j.1672-3813.2017.03.011

        O322

        A

        2017-06-13;

        2017-08-24

        國家自然科學(xué)基金(11472256,11402244),中物院院長(zhǎng)基金(YZ2015011);中物院發(fā)展基金(2015B0201024)

        肖世富(1970-),男,四川中江人,博士,研究員,主要研究方向?yàn)槿嵝远囿w系統(tǒng)動(dòng)力學(xué)、計(jì)算固體力學(xué)。

        (責(zé)任編輯耿金花)

        猜你喜歡
        勢(shì)阱振子非對(duì)稱
        含有陡峭勢(shì)阱和凹凸非線性項(xiàng)的Kirchhoff型問題的多重正解
        分?jǐn)?shù)階量子力學(xué)下的二維無限深方勢(shì)阱
        時(shí)空分?jǐn)?shù)階量子力學(xué)下的δ勢(shì)阱
        對(duì)稱三勢(shì)阱玻色—愛因斯坦凝聚體的非線性效應(yīng)
        彈簧振子問題的分析與求解
        非對(duì)稱Orlicz差體
        非線性Duffing擾動(dòng)振子共振機(jī)制的研究
        點(diǎn)數(shù)不超過20的旗傳遞非對(duì)稱2-設(shè)計(jì)
        基于近似熵和混沌振子的電力諧波檢測(cè)與估計(jì)
        非對(duì)稱負(fù)載下矩陣變換器改進(jìn)型PI重復(fù)控制
        99精品人妻少妇一区二区| 成人国产av精品麻豆网址| 91色老久久偷偷精品蜜臀懂色 | 国产成人精品视频网站| av在线播放免费观看| 国产99视频精品免视看7| 午夜dj在线观看免费视频| 日本不卡一区二区高清中文| 在线观看亚洲视频一区二区| 女人张开腿让男人桶爽| 日韩内射美女人妻一区二区三区 | 日本人妻伦理片在线观看| 中文字幕精品人妻在线| 久久久久成人片免费观看蜜芽 | 亚洲午夜成人片| 久久亚洲宅男天堂网址| 色又黄又爽18禁免费网站现观看| 老熟妻内射精品一区| 亚洲网站免费看| 成人国产一区二区三区av| 国产av无码专区亚洲avjulia | av少妇偷窃癖在线观看| 精品国产一区二区三区av新片| 中文字幕精品一区二区精品| 亚洲av无码一区二区三区在线| 国产aⅴ丝袜旗袍无码麻豆| 国产性虐视频在线观看| 东北女人毛多水多牲交视频 | 国产精品,在线点播影院| 日韩人妻中文字幕专区| 国产免费人成视频在线观看| 无码中文日韩Av| 亚洲天堂男人的av天堂| 亚洲av成人网| 91av手机在线观看| 国产一区二区免费在线观看视频| av网站免费线看精品| 97久久天天综合色天天综合色hd| 无码精品国产午夜| 久久精品国产亚洲av四叶草| 在线不卡av片免费观看|