• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    基于MATLAB的強(qiáng)迫振動(dòng)達(dá)芬方程的非線性幅頻響應(yīng)分析*③

    2016-09-05 03:42:15陳趙江
    物理通報(bào) 2016年6期
    關(guān)鍵詞:幅頻掃頻幅值

    陳趙江 陳 敏

    (浙江師范大學(xué)數(shù)理與信息工程學(xué)院 浙江 金華 321004)

    ?

    基于MATLAB的強(qiáng)迫振動(dòng)達(dá)芬方程的非線性幅頻響應(yīng)分析*③

    陳趙江陳 敏

    (浙江師范大學(xué)數(shù)理與信息工程學(xué)院浙江 金華321004)

    利用MATLAB研究了強(qiáng)迫振動(dòng)達(dá)芬(Duffing)方程的非線性幅頻響應(yīng)特性,分析了達(dá)芬方程非線性幅頻響應(yīng)近似解析求解和數(shù)值求解的方法和步驟,給出了相應(yīng)的MATLAB求解程序,并將解析解與數(shù)值解結(jié)果進(jìn)行了比較.仿真程序和結(jié)果能夠加深學(xué)生對(duì)非線性振動(dòng)相關(guān)知識(shí)的理解,提高大學(xué)物理及相關(guān)力學(xué)課程的課堂教學(xué)效果.

    達(dá)芬方程非線性振動(dòng)非線性幅頻響應(yīng)跳躍現(xiàn)象MATLAB

    強(qiáng)迫振動(dòng)達(dá)芬方程(Duffingequation)在非線性振動(dòng)領(lǐng)域是一個(gè)非常經(jīng)典的實(shí)例,被廣泛用來(lái)說(shuō)明跳躍現(xiàn)象(Jumpingphenomena)、頻響曲線彎曲和其他非線性行為.理解這個(gè)簡(jiǎn)單的低階非線性方程的振動(dòng)特性是學(xué)習(xí)更為復(fù)雜的非線性系統(tǒng)的基礎(chǔ).達(dá)芬方程解的跳躍現(xiàn)象與其非線性頻響特性有直接的關(guān)系,因此對(duì)頻響特性的研究非常重要.在目前大學(xué)物理以及相關(guān)力學(xué)課程的教學(xué)中,學(xué)生對(duì)非線性振動(dòng)知識(shí)了解甚少,一知半解.如果能利用軟件程序幫助學(xué)生了解掌握晦澀難懂的非線性振動(dòng)知識(shí),這將大大有益于相關(guān)教學(xué)工作的開(kāi)展.目前已往文獻(xiàn)對(duì)達(dá)芬方程非線性幅頻響應(yīng)的分析大多利用MAPLE和MATHEMATICA數(shù)學(xué)軟件且缺少數(shù)值解和解析解的比較分析.本文給出了達(dá)芬方程非線性幅頻響應(yīng)近似解析求解和數(shù)值求解的方法和步驟,利用MATLAB的強(qiáng)大計(jì)算功能編寫(xiě)了相應(yīng)的求解程序,并將解析解與數(shù)值解結(jié)果進(jìn)行了比較分析,以滿足學(xué)生對(duì)非線性振動(dòng)知識(shí)的進(jìn)一步理解.

    1 達(dá)芬方程非線性幅頻響應(yīng)的理論分析

    含阻尼和外部驅(qū)動(dòng)力的無(wú)量綱達(dá)芬方程一般具有如下形式[1]

    (1)

    其中u是位移,t是時(shí)間,ζ是阻尼比,γ是立方剛度系數(shù),F(xiàn)和Ω分別為激勵(lì)幅值和激勵(lì)頻率.需注意的是,式(1)中所有物理量均為無(wú)量綱量,因此Ω其實(shí)是激勵(lì)頻率和系統(tǒng)固有頻率之比,即

    Ω=1+εσ

    (2)

    則式(1)可改寫(xiě)為

    (3)

    我們采用多尺度法[1]對(duì)上式進(jìn)行近似解析求解,引入尺度T0=t和T1=εt

    (4)

    (5)

    (6)

    的式(3)的解.把式(4)~(6)代入式(3),并令ε的同次冪的系數(shù)相等,我們得到

    (7)

    (8)

    式(7)的通解可以表示為

    (9)

    (10)

    把A表示成如下的級(jí)數(shù)形式

    (11)

    則式(10)改寫(xiě)成

    (12)

    (13)

    將式(13)分成實(shí)部和虛部,得出

    (14)

    (15)

    式中a和β都是實(shí)數(shù),引進(jìn)變換

    (16)

    把式(14)和(15)化為自治系統(tǒng),由式(16)

    φ′=σ-β′

    (17)

    把式(16)代入式(14)和(15),并利用式(17)可得

    (18)

    (19)

    由于我們考慮穩(wěn)態(tài)運(yùn)動(dòng),即a和φ為常數(shù)值,設(shè)a′=0和φ′=0,可得到

    (20)

    (21)

    對(duì)式(20)和(21)取平方,將結(jié)果相加,得出

    (22)

    通常稱之為頻率響應(yīng)方程.利用式(2)我們將上式改寫(xiě)為Ω關(guān)于a的形式, 得到

    (23)

    從式(23)可知,對(duì)于某一幅值響應(yīng)a,其對(duì)應(yīng)兩個(gè)激勵(lì)頻率Ω1,2,但這兩個(gè)頻率所對(duì)應(yīng)的響應(yīng)a可能是不穩(wěn)定的.為了表征解的穩(wěn)定性,考慮如下的動(dòng)力學(xué)系統(tǒng)方程

    (24)

    其中,x1和x2為狀態(tài)量.則上述系統(tǒng)的雅可比矩陣為

    (25)

    根據(jù)式(18)和式(19),在本系統(tǒng)中狀態(tài)向量分別為a和φ,f1和f2的表達(dá)式如下

    (26)

    因而可得系統(tǒng)的雅可比矩陣為

    (27)

    穩(wěn)態(tài)運(yùn)動(dòng)的穩(wěn)定性依賴于雅可比矩陣的特征值,式(27)所對(duì)應(yīng)的特征方程如下

    (28)

    展開(kāi)此行列式,可得

    (29)

    本系統(tǒng)的本征值如下

    (30)

    (31)

    其中

    2 達(dá)芬振動(dòng)方程非線性幅頻響應(yīng)分析的MATLAB實(shí)現(xiàn)

    2.1MATLAB解析求解程序

    從式(22)可以看出為了得到頻率響應(yīng)曲線,可解出a2關(guān)于σ的函數(shù),或者可以解出σ關(guān)于a的函數(shù).顯然,后一種方法比較簡(jiǎn)單.因此,在編程時(shí)根據(jù)式(23)進(jìn)行計(jì)算.以下為MATLAB求解程序:

    %強(qiáng)迫振動(dòng)達(dá)芬方程幅頻響應(yīng)曲線的解析求解

    %參數(shù)設(shè)定

    epsilon=0.2; %小參數(shù)

    gamma1=4; %非線性參數(shù)

    F1=0.5; %激勵(lì)幅值

    zeta1=0.25;%阻尼項(xiàng)

    a=linspace(0.1,2.0,100); % 計(jì)算響應(yīng)幅值范圍

    forii=1:1:length(a)

    %根據(jù)(23)式,對(duì)應(yīng)一個(gè)響應(yīng)幅值,可能存在兩個(gè)激勵(lì)頻率值;

    %Omega1為較小的激勵(lì)頻率值

    Omega1(ii)=(1+3*epsilon*gamma1/8*a(ii)^2-sqrt((epsilon*F1/(2*a(ii)))^2-(epsilon*

    zeta1)^2));

    %Omega1對(duì)應(yīng)的本征函數(shù)值

    lamOmega11=sqrt(-((Omega1(ii)-1)/epsilon-3*gamma1/8*a(ii)^2)*((Omega1(ii)-1)/epsilon-9*gamma1/8*a(ii)^2))-zeta1;

    lamOmega12=-sqrt(-((Omega1(ii)-1)/epsilon-3*gamma1/8*a(ii)^2)*((Omega1(ii)-1)/epsilon-9*gamma1/8*a(ii)^2))-zeta1;

    %Omega2為較大的激勵(lì)頻率值

    Omega2(ii)=(1+3*epsilon*gamma1/

    8*a(ii)^2+sqrt((epsilon*F1/(2*a(ii)))^2-

    (epsilon*zeta1)^2));

    %Omega2對(duì)應(yīng)的本征函數(shù)值

    lamOmega21=sqrt(-((Omega2(ii)-1)/epsilon-3*gamma1/8*a(ii)^2)*((Omega2(ii)-1)/epsilon-9*gamma1/8*a(ii)^2))-zeta1;

    lamOmega22=-sqrt(-((Omega2(ii)-1)/epsilon-3*gamma1/8*a(ii)^2)*((Omega2(ii)-1)/epsilon-9*gamma1/8*a(ii)^2))-zeta1;

    iflamOmega11==conj(lamOmega21)

    %如果兩個(gè)本征值相等退出計(jì)算(即此時(shí)計(jì)算到頻響曲線最高點(diǎn))

    plot(Omega1(ii),a(ii),'k.','MarkerSize',15')

    holdon

    break

    else

    FG=((Omega1(ii)-1)/epsilon-3*gamma1/8*a(ii)^2)*((Omega1(ii)-1)/epsilon-9*

    gamma1/8*a(ii)^2)+zeta1^2;% 對(duì)應(yīng)上文中的(31)式

    ifFG<0

    plot(Omega1(ii),a(ii),'color',[.5 .5 .5],'MarkerSize',15) %不穩(wěn)定點(diǎn)

    else

    plot(Omega1(ii),a(ii),'k.','MarkerSize',15)%穩(wěn)定點(diǎn)

    holdon

    end

    GF=((Omega2(ii)-1)/epsilon-3*gamma1/

    8*a(ii)^2)*((Omega2(ii)-1)/epsilon-9*

    gamma1/8*a(ii)^2)+zeta1^2;% 對(duì)應(yīng)上文中的(31)式

    ifGF<0

    plot(Omega2(ii),a(ii),'.','color',[.5 .5 .5],'MarkerSize',15) %不穩(wěn)定點(diǎn)

    else

    plot(Omega2(ii),a(ii),'k.','MarkerSize',15)%穩(wěn)定點(diǎn)

    holdon

    end

    xlabel('Omega') %x軸標(biāo)題

    ylabel('a')%y軸標(biāo)題

    xlim([0.5,1.5]); %x軸范圍

    ylim([0,2]);%y軸范圍

    end%滿足條件退出循環(huán)(計(jì)算到幅頻響應(yīng)曲線頂點(diǎn))

    end%結(jié)束對(duì)不同a值的計(jì)算

    2.2MATLAB數(shù)值求解程序

    (32)

    在MATLAB中定義如下函數(shù)并保存為duffing.m

    functiondy=duffing(t,y)

    globalepsilongamma1F1zeta1Omega

    dy=zeros(2,1);

    dy(1)=-zeta1*y(1)+F1/2*sin(y(2));

    dy(2)=(Omega-1)/epsilon-3*gamma1/

    8*y(1)^2+F1/(2*y(1))*cos(y(2));

    end

    (33)

    在MATLAB中定義如下函數(shù)并保存為duffing1.m

    functiondydt=duffing1(t,y)

    globalepsilongamma1F1zeta1Omega

    dydt=zeros(2,1);

    dydt(1)=y(2);

    dydt(2)=-2*epsilon*zeta1*y(2)-y(1)-epsilon*gamma1*y(1)^3+epsilon*F1*cos

    (Omega*t);

    end

    達(dá)芬方程幅頻響應(yīng)數(shù)值求解主程序main.m如下:

    clc

    %closeall

    clearall

    globalepsilongamma1F1zeta1Omega樣% 定義全局變量

    %參數(shù)設(shè)定

    epsilon=0.2; %小參數(shù)

    gamma1 = 4; % 非線性參數(shù)

    F1=0.5;%激勵(lì)幅值

    zeta1=0.25; %阻尼項(xiàng)

    np=400;% 計(jì)算的頻率數(shù)目

    Omega1=linspace(.5,1.5,np); % 正向掃頻時(shí)的頻率

    Omega2=linspace(1.5,.5,np); % 反向掃頻時(shí)的頻率

    yy=[];

    Y0=[0.1 0.1];% 初始值

    fori=1:1:length(Omega1)

    Omega=Omega1(i);

    [T,Y] =ode45(@duffing,[0 400],Y0);

    %[T,Y] =ode45(@duffing1,[0 400],Y0);

    nn=length(Y(:,1));

    ymax=max(Y(nn-round(nn/2):nn,1));% 取穩(wěn)態(tài)幅值

    ymax1=max(Y(nn-round(nn/2):nn,2));

    Y0=[ymaxymax1];% 下一次計(jì)算的初始值

    yy=[yy;ymax];% 保存計(jì)算結(jié)果到y(tǒng)y向量

    end

    plot(Omega1,yy,'k','LineWidth',1.5)% 正向掃頻曲線用黑色表示

    holdon

    % 反向掃頻

    yy1=[];

    ymax1=max(Y(nn-round(nn/2):nn,1));% 反向掃頻的初始值

    ymax2=max(Y(nn-round(nn/2):nn,2));% 反向掃頻的初始值

    forj=1:1:length(Omega2)

    Omega=Omega2(j)

    [T,Y1] =ode45(@duffing,[0 400],[ymax1ymax2]);

    %[T,Y] =ode45(@duffing1,[0 400],Y0);

    nn1=length(Y1(:,1));

    ymax1=max(Y1(nn1-round(nn1/2):nn1,1));

    ymax2=max(Y1(nn1-round(nn1/2):nn1,2));

    yy1=[yy1;ymax1];

    end

    plot(Omega2,yy1,'color',[.5 .5 .5],

    'LineWidth',1.5)% 反向掃頻曲線用灰色表示

    xlabel('Omega') %x軸標(biāo)題

    ylabel('a')%y軸標(biāo)題

    xlim([0.5,1.5]); %x軸范圍

    ylim([0,2]); %y軸范圍

    3 求解結(jié)果分析

    (a)隨非線性參數(shù)的變化

    (b)隨激勵(lì)幅值的變化

    其次,我們利用2.2節(jié)的幅頻響應(yīng)數(shù)值求解程序?qū)?qiáng)迫振動(dòng)達(dá)芬方程的一階近似幅頻響應(yīng)進(jìn)行求解分析,其中使用的MATLAB子函數(shù)為duffing.m.圖2中黑色和灰色曲線分別為正向和反向掃頻計(jì)算結(jié)果.從圖2可以發(fā)現(xiàn)當(dāng)正向掃頻(即頻率比Ω增大)時(shí),響應(yīng)幅值一直增加.當(dāng)進(jìn)一步增加時(shí),就發(fā)生從A點(diǎn)到B點(diǎn)的跳躍并伴隨著響應(yīng)幅值的突然減小.隨后,響應(yīng)幅值隨Ω的增大而減小.對(duì)應(yīng)著A點(diǎn)的最大振幅只有從較低的頻率開(kāi)始增加才可能達(dá)到.在另一方面,反向掃頻(即頻率比Ω減小)時(shí),響應(yīng)幅值也增加.當(dāng)Ω進(jìn)一步減小時(shí)就發(fā)生從C點(diǎn)到D點(diǎn)的跳躍并伴隨著響應(yīng)幅值的突然增加.隨后,響應(yīng)幅值隨Ω的減小而減小.從圖中也可發(fā)現(xiàn)解析解中的不穩(wěn)定區(qū)域(從A點(diǎn)到C點(diǎn))在數(shù)值求解時(shí)是無(wú)法得到的,幅頻響應(yīng)一階近似解析解和一階近似數(shù)值解結(jié)果符合很好.

    圖2 強(qiáng)迫振動(dòng)達(dá)芬方程的幅頻響應(yīng)數(shù)值與解析解的比較

    最后,我們對(duì)強(qiáng)迫振動(dòng)達(dá)芬方程的幅頻響應(yīng)做進(jìn)一步的數(shù)值分析,其中使用的MATLAB子函數(shù)為duffing1.m.計(jì)算結(jié)果如圖3所示,從圖中可以發(fā)現(xiàn)數(shù)值解與近似解析解存在一定差異,且偏離Ω=1 越大兩者的差別越大.這是因?yàn)檫@里數(shù)值求解的是最初的不包含近似的微分方程式(3),而解析解是滿足Ω≈1的一階近似解.因此,當(dāng)Ω越接近1時(shí),兩者解的結(jié)果越接近.從上述分析也可知近似解析解式(22)是在弱阻尼、弱非線性系數(shù)并接近共振時(shí)得出的結(jié)果,并不適用于所有的情況,這一點(diǎn)在學(xué)習(xí)的時(shí)候要特別注意.

    圖3 強(qiáng)迫振動(dòng)達(dá)芬方程的幅頻響應(yīng)數(shù)值與解析解的比較

    4 總結(jié)

    本文對(duì)強(qiáng)迫振動(dòng)達(dá)芬方程的非線性幅頻響應(yīng)特性進(jìn)行了研究,分析了達(dá)芬方程非線性幅頻響應(yīng)解析求解和數(shù)值求解的方法和步驟,給出了相應(yīng)的MATLAB求解程序,并將解析解與數(shù)值解結(jié)果進(jìn)行了比較分析.仿真程序和結(jié)果能夠加深學(xué)生對(duì)非線性振動(dòng)的理解,提高大學(xué)物理以及相關(guān)力學(xué)課程的教學(xué)效果.此外,雖然本文只討論了達(dá)芬方程的非線性幅頻響應(yīng)特性,但仿照本文的求解步驟和程序也容易對(duì)達(dá)芬方程的非線性相頻特性進(jìn)行分析.

    1A·H·奈弗,D·T·穆克著. 非線性振動(dòng). 宋家骕, 羅惟德, 陳守吉譯. 北京:高等教育出版社, 1980

    2張志涌. 精通MATLABR2011a. 北京:北京航空航天大學(xué)出版社, 2011

    陳趙江(1980-),男,博士,主要從事大學(xué)物理教學(xué)和研究.

    2016-03-08)

    ③浙江師范大學(xué)“十二五”省級(jí)實(shí)驗(yàn)教學(xué)示范中心重點(diǎn)建設(shè)項(xiàng)目和2015年度浙江師范大學(xué)教改項(xiàng)目資助.

    猜你喜歡
    幅頻掃頻幅值
    杠桿型串聯(lián)非線性能量阱整星隔振系統(tǒng)的振動(dòng)控制
    幅頻電透視在探查煤層底板水及注漿檢驗(yàn)中的應(yīng)用
    煤炭與化工(2022年1期)2022-03-19 03:12:52
    正弦掃頻速率對(duì)結(jié)構(gòu)響應(yīng)的影響分析
    寬帶高速掃頻信號(hào)源的高精度功率控制設(shè)計(jì)
    帶電等效阻抗掃頻測(cè)試的互感器繞組及外絕緣隱患快速識(shí)別新技術(shù)的應(yīng)用研究
    電子制作(2017年8期)2017-06-05 09:36:15
    基于S變換的交流電網(wǎng)幅值檢測(cè)系統(tǒng)計(jì)算機(jī)仿真研究
    電子制作(2017年7期)2017-06-05 09:36:13
    一種線性掃頻干擾信號(hào)的參數(shù)估計(jì)方法
    正序電壓幅值檢測(cè)及諧波抑制的改進(jìn)
    探測(cè)器非線性對(duì)可見(jiàn)光通信系統(tǒng)幅頻響應(yīng)的影響
    低壓電力線信道脈沖噪聲的幅值與寬度特征
    亚洲午夜精品一区,二区,三区| 久久久久久免费高清国产稀缺| 亚洲第一欧美日韩一区二区三区 | 精品久久久精品久久久| 免费观看人在逋| 欧美av亚洲av综合av国产av| 人妻人人澡人人爽人人| 国产av国产精品国产| 亚洲情色 制服丝袜| 一本色道久久久久久精品综合| 后天国语完整版免费观看| 日本a在线网址| 国产免费现黄频在线看| 国产精品九九99| 19禁男女啪啪无遮挡网站| 国产野战对白在线观看| 法律面前人人平等表现在哪些方面 | 国产免费现黄频在线看| 日韩,欧美,国产一区二区三区| 国产男女超爽视频在线观看| 亚洲av欧美aⅴ国产| 日日爽夜夜爽网站| 亚洲精品自拍成人| 亚洲av日韩在线播放| 麻豆乱淫一区二区| 18禁观看日本| 国产一区有黄有色的免费视频| 欧美午夜高清在线| 岛国毛片在线播放| 永久免费av网站大全| 如日韩欧美国产精品一区二区三区| 18禁黄网站禁片午夜丰满| 大片电影免费在线观看免费| 黄色怎么调成土黄色| 在线观看舔阴道视频| 狂野欧美激情性bbbbbb| 狠狠精品人妻久久久久久综合| 男女无遮挡免费网站观看| 精品国产乱码久久久久久男人| 日本vs欧美在线观看视频| 亚洲中文日韩欧美视频| 好男人电影高清在线观看| 19禁男女啪啪无遮挡网站| 大陆偷拍与自拍| 午夜福利视频在线观看免费| 亚洲欧美色中文字幕在线| 在线永久观看黄色视频| 亚洲精品国产精品久久久不卡| 后天国语完整版免费观看| 我的亚洲天堂| 18禁观看日本| 亚洲人成77777在线视频| 午夜精品久久久久久毛片777| 天堂8中文在线网| 啦啦啦啦在线视频资源| 日韩熟女老妇一区二区性免费视频| 久久久国产成人免费| 精品少妇久久久久久888优播| 国产亚洲av高清不卡| 国产欧美日韩一区二区三 | 一级毛片电影观看| 蜜桃国产av成人99| 国产免费现黄频在线看| 中文精品一卡2卡3卡4更新| 两性午夜刺激爽爽歪歪视频在线观看 | 中亚洲国语对白在线视频| 日本猛色少妇xxxxx猛交久久| 老司机在亚洲福利影院| 日韩欧美一区视频在线观看| 国产精品.久久久| 亚洲全国av大片| 伦理电影免费视频| 2018国产大陆天天弄谢| 久久热在线av| 午夜91福利影院| 亚洲五月色婷婷综合| 午夜福利视频在线观看免费| 美女高潮喷水抽搐中文字幕| 黄频高清免费视频| 中文欧美无线码| 国产成人欧美在线观看 | 岛国毛片在线播放| 免费高清在线观看视频在线观看| 亚洲精品乱久久久久久| 国产1区2区3区精品| 久久久久久人人人人人| 99国产综合亚洲精品| 国产精品亚洲av一区麻豆| 高清欧美精品videossex| 免费少妇av软件| av网站免费在线观看视频| 一级毛片电影观看| 亚洲精品在线美女| 大香蕉久久网| e午夜精品久久久久久久| 我要看黄色一级片免费的| 国产在线观看jvid| 国产精品偷伦视频观看了| 久久99热这里只频精品6学生| 久久久久久久国产电影| 久久精品国产亚洲av高清一级| 亚洲国产中文字幕在线视频| 蜜桃国产av成人99| 国产av精品麻豆| 成年女人毛片免费观看观看9 | 国产高清videossex| 国产精品一二三区在线看| 日韩,欧美,国产一区二区三区| 午夜老司机福利片| 老司机影院成人| av视频免费观看在线观看| 免费人妻精品一区二区三区视频| 国产成人免费无遮挡视频| 国产深夜福利视频在线观看| 亚洲国产精品一区二区三区在线| 天堂8中文在线网| 欧美+亚洲+日韩+国产| 丝袜美腿诱惑在线| 自拍欧美九色日韩亚洲蝌蚪91| 手机成人av网站| 欧美变态另类bdsm刘玥| bbb黄色大片| 亚洲精品国产av成人精品| 男女无遮挡免费网站观看| 亚洲人成电影观看| 亚洲精品国产av成人精品| 午夜福利乱码中文字幕| 欧美日韩福利视频一区二区| 91九色精品人成在线观看| 国产在线观看jvid| tocl精华| 成年人午夜在线观看视频| 91麻豆精品激情在线观看国产 | 中文字幕精品免费在线观看视频| 老汉色∧v一级毛片| 高清视频免费观看一区二区| 亚洲国产av影院在线观看| av不卡在线播放| 久久精品aⅴ一区二区三区四区| 久久九九热精品免费| 日日爽夜夜爽网站| 午夜免费观看性视频| 欧美精品一区二区免费开放| 搡老熟女国产l中国老女人| 欧美在线黄色| 欧美亚洲 丝袜 人妻 在线| 窝窝影院91人妻| 涩涩av久久男人的天堂| 亚洲五月婷婷丁香| 中文字幕精品免费在线观看视频| 亚洲久久久国产精品| 久久这里只有精品19| 久久久久网色| 亚洲精品粉嫩美女一区| 夫妻午夜视频| 黑人巨大精品欧美一区二区蜜桃| 国产成人精品在线电影| 亚洲av男天堂| 大码成人一级视频| 久久香蕉激情| 国产亚洲欧美在线一区二区| 狠狠婷婷综合久久久久久88av| 日本av免费视频播放| 久久久精品94久久精品| 热99国产精品久久久久久7| 国产深夜福利视频在线观看| 成年人免费黄色播放视频| 国产成人影院久久av| 亚洲国产中文字幕在线视频| 免费在线观看视频国产中文字幕亚洲 | 天天躁狠狠躁夜夜躁狠狠躁| 最近中文字幕2019免费版| 亚洲人成电影免费在线| 日日摸夜夜添夜夜添小说| 不卡一级毛片| 亚洲国产精品999| 欧美成狂野欧美在线观看| 欧美日韩国产mv在线观看视频| 日韩电影二区| 亚洲伊人久久精品综合| 另类精品久久| cao死你这个sao货| 国产成人精品久久二区二区免费| 人妻 亚洲 视频| 一区在线观看完整版| 夜夜骑夜夜射夜夜干| 他把我摸到了高潮在线观看 | 美女福利国产在线| 18禁国产床啪视频网站| 狂野欧美激情性xxxx| 日韩中文字幕欧美一区二区| 亚洲av日韩在线播放| 国产精品一区二区免费欧美 | 正在播放国产对白刺激| 免费黄频网站在线观看国产| 亚洲人成电影免费在线| 亚洲一区中文字幕在线| 天天影视国产精品| 免费在线观看视频国产中文字幕亚洲 | 99九九在线精品视频| 波多野结衣av一区二区av| 十八禁人妻一区二区| 美女福利国产在线| 他把我摸到了高潮在线观看 | 亚洲欧洲精品一区二区精品久久久| 老熟妇乱子伦视频在线观看 | 国产亚洲午夜精品一区二区久久| 蜜桃国产av成人99| 男女边摸边吃奶| 国产成人a∨麻豆精品| 国产精品香港三级国产av潘金莲| 免费观看a级毛片全部| 色婷婷久久久亚洲欧美| 日韩视频一区二区在线观看| 脱女人内裤的视频| 亚洲视频免费观看视频| 一级片'在线观看视频| 日韩中文字幕欧美一区二区| av网站免费在线观看视频| 国产1区2区3区精品| 亚洲精品久久成人aⅴ小说| 99热国产这里只有精品6| a 毛片基地| 国产在线视频一区二区| 久久午夜综合久久蜜桃| 成年动漫av网址| 国产欧美日韩综合在线一区二区| 亚洲熟女毛片儿| 国产亚洲精品久久久久5区| 黄色a级毛片大全视频| 在线观看一区二区三区激情| 欧美人与性动交α欧美精品济南到| 热99国产精品久久久久久7| 黄色视频不卡| 精品第一国产精品| 国产精品久久久av美女十八| 亚洲欧美成人综合另类久久久| 一级a爱视频在线免费观看| 日韩视频一区二区在线观看| 日日摸夜夜添夜夜添小说| 午夜福利一区二区在线看| 亚洲伊人久久精品综合| 欧美成人午夜精品| 国产欧美日韩一区二区三 | 50天的宝宝边吃奶边哭怎么回事| 国产精品久久久久久精品古装| av有码第一页| 人妻久久中文字幕网| 黄网站色视频无遮挡免费观看| 久久久国产精品麻豆| 中文精品一卡2卡3卡4更新| 亚洲三区欧美一区| 日韩熟女老妇一区二区性免费视频| 成人国产av品久久久| 熟女少妇亚洲综合色aaa.| 丰满饥渴人妻一区二区三| 狠狠狠狠99中文字幕| 真人做人爱边吃奶动态| 国产伦理片在线播放av一区| 欧美少妇被猛烈插入视频| 亚洲精品一区蜜桃| 制服诱惑二区| 建设人人有责人人尽责人人享有的| 黄频高清免费视频| 搡老熟女国产l中国老女人| 十八禁网站免费在线| 国产激情久久老熟女| 人人妻,人人澡人人爽秒播| 两性夫妻黄色片| 精品福利观看| 香蕉丝袜av| 精品国产乱子伦一区二区三区 | 高潮久久久久久久久久久不卡| 国产真人三级小视频在线观看| 国产成人一区二区三区免费视频网站| 黄色 视频免费看| av超薄肉色丝袜交足视频| 夜夜骑夜夜射夜夜干| 亚洲中文av在线| av一本久久久久| 男女午夜视频在线观看| 精品一区二区三区四区五区乱码| 精品国产超薄肉色丝袜足j| 大型av网站在线播放| 亚洲国产欧美网| 中文字幕人妻丝袜一区二区| 国产一区二区 视频在线| 精品一区二区三卡| 各种免费的搞黄视频| 男人舔女人的私密视频| 亚洲人成电影免费在线| 久久久精品区二区三区| 纵有疾风起免费观看全集完整版| 在线十欧美十亚洲十日本专区| 久久久久精品人妻al黑| 成人国产av品久久久| 久久久精品国产亚洲av高清涩受| 精品国产乱码久久久久久小说| 亚洲精品久久午夜乱码| 天天躁夜夜躁狠狠躁躁| 免费观看人在逋| av天堂久久9| 少妇粗大呻吟视频| 国产精品av久久久久免费| 国产亚洲欧美精品永久| 麻豆乱淫一区二区| 午夜福利免费观看在线| 亚洲精品国产一区二区精华液| 热re99久久精品国产66热6| videosex国产| 国产激情久久老熟女| 欧美亚洲日本最大视频资源| 精品国产乱码久久久久久男人| 91大片在线观看| 亚洲 国产 在线| cao死你这个sao货| 脱女人内裤的视频| 欧美在线黄色| 香蕉国产在线看| 秋霞在线观看毛片| 汤姆久久久久久久影院中文字幕| 一边摸一边做爽爽视频免费| 最近最新中文字幕大全免费视频| 一本综合久久免费| 免费女性裸体啪啪无遮挡网站| 中文字幕高清在线视频| 一个人免费在线观看的高清视频 | 国产日韩欧美亚洲二区| 日本撒尿小便嘘嘘汇集6| 女人被躁到高潮嗷嗷叫费观| 日韩有码中文字幕| 欧美亚洲日本最大视频资源| 亚洲一卡2卡3卡4卡5卡精品中文| 美女视频免费永久观看网站| 久久久国产成人免费| 妹子高潮喷水视频| 最新的欧美精品一区二区| 777久久人妻少妇嫩草av网站| 亚洲av日韩在线播放| 中文字幕精品免费在线观看视频| 51午夜福利影视在线观看| 久久精品国产综合久久久| 黄色视频不卡| 精品国产乱子伦一区二区三区 | 99精品久久久久人妻精品| 国内毛片毛片毛片毛片毛片| 欧美亚洲 丝袜 人妻 在线| 久久ye,这里只有精品| 欧美成人午夜精品| 成人亚洲精品一区在线观看| 日本猛色少妇xxxxx猛交久久| 老司机靠b影院| 人成视频在线观看免费观看| 精品国内亚洲2022精品成人 | 亚洲av国产av综合av卡| 9191精品国产免费久久| 在线永久观看黄色视频| 中文字幕人妻熟女乱码| 美女国产高潮福利片在线看| 国产欧美日韩精品亚洲av| av网站在线播放免费| 亚洲 国产 在线| 这个男人来自地球电影免费观看| 欧美乱码精品一区二区三区| 国产又爽黄色视频| 国产男女内射视频| 国产精品久久久人人做人人爽| 国产精品免费大片| 国产精品国产av在线观看| √禁漫天堂资源中文www| 日韩 欧美 亚洲 中文字幕| 老司机午夜福利在线观看视频 | 亚洲国产精品999| 亚洲av电影在线观看一区二区三区| 国产日韩欧美在线精品| 99香蕉大伊视频| 自拍欧美九色日韩亚洲蝌蚪91| 日本av手机在线免费观看| 亚洲国产成人一精品久久久| 51午夜福利影视在线观看| 国产高清videossex| 男女午夜视频在线观看| av又黄又爽大尺度在线免费看| 亚洲第一av免费看| 午夜成年电影在线免费观看| 国产精品国产av在线观看| 欧美大码av| 一本一本久久a久久精品综合妖精| 淫妇啪啪啪对白视频 | 国产成人av教育| 嫩草影视91久久| 日韩一卡2卡3卡4卡2021年| 涩涩av久久男人的天堂| 大片电影免费在线观看免费| 久久毛片免费看一区二区三区| av一本久久久久| 久久亚洲国产成人精品v| 精品一区二区三区av网在线观看 | 免费观看人在逋| 免费日韩欧美在线观看| 电影成人av| 亚洲五月色婷婷综合| 汤姆久久久久久久影院中文字幕| 亚洲国产欧美日韩在线播放| 精品一区二区三区四区五区乱码| 国产成人系列免费观看| √禁漫天堂资源中文www| 1024香蕉在线观看| 亚洲 欧美一区二区三区| 国产免费现黄频在线看| 99国产精品一区二区三区| cao死你这个sao货| 天堂俺去俺来也www色官网| 国产av又大| 在线亚洲精品国产二区图片欧美| 亚洲av男天堂| 精品国产一区二区三区四区第35| 黑人猛操日本美女一级片| 热99久久久久精品小说推荐| 免费在线观看日本一区| 久久国产精品大桥未久av| 欧美 亚洲 国产 日韩一| 最近最新中文字幕大全免费视频| 黑人巨大精品欧美一区二区mp4| 青草久久国产| 99久久人妻综合| 可以免费在线观看a视频的电影网站| 亚洲国产看品久久| 成人18禁高潮啪啪吃奶动态图| 精品人妻熟女毛片av久久网站| 午夜福利视频精品| 免费女性裸体啪啪无遮挡网站| 曰老女人黄片| 大香蕉久久成人网| 好男人电影高清在线观看| 在线天堂中文资源库| 少妇猛男粗大的猛烈进出视频| 国产有黄有色有爽视频| 国产精品一区二区精品视频观看| 两性夫妻黄色片| 自线自在国产av| 中文字幕色久视频| 亚洲精品一二三| 午夜免费鲁丝| 免费高清在线观看日韩| 久久久久网色| 欧美成狂野欧美在线观看| 成年av动漫网址| 久久久久视频综合| 宅男免费午夜| 热re99久久国产66热| 日本a在线网址| 97精品久久久久久久久久精品| 欧美日韩视频精品一区| 国产在线一区二区三区精| 婷婷丁香在线五月| 亚洲av美国av| 9191精品国产免费久久| 美女大奶头黄色视频| 岛国在线观看网站| 国产亚洲精品久久久久5区| 妹子高潮喷水视频| 国产精品av久久久久免费| 免费看十八禁软件| 99国产精品一区二区三区| 国产成人精品在线电影| 免费在线观看日本一区| 18禁国产床啪视频网站| 多毛熟女@视频| 美女午夜性视频免费| 免费不卡黄色视频| 自线自在国产av| 99香蕉大伊视频| 久热爱精品视频在线9| 国产精品久久久久久人妻精品电影 | 亚洲熟女精品中文字幕| 性色av一级| 最近最新免费中文字幕在线| 国产精品久久久久成人av| 亚洲人成77777在线视频| 黑人欧美特级aaaaaa片| 国产91精品成人一区二区三区 | 窝窝影院91人妻| 一二三四社区在线视频社区8| 久久午夜综合久久蜜桃| 欧美午夜高清在线| 亚洲第一欧美日韩一区二区三区 | 麻豆乱淫一区二区| 别揉我奶头~嗯~啊~动态视频 | 久久影院123| 大香蕉久久网| 亚洲午夜精品一区,二区,三区| 国产一区二区三区av在线| 涩涩av久久男人的天堂| 搡老乐熟女国产| 成人亚洲精品一区在线观看| 美女国产高潮福利片在线看| 久久香蕉激情| 精品国产一区二区三区四区第35| 黄色怎么调成土黄色| 精品国产乱码久久久久久男人| 精品少妇一区二区三区视频日本电影| 国产极品粉嫩免费观看在线| 国产淫语在线视频| 中文字幕av电影在线播放| 国产免费福利视频在线观看| 999久久久精品免费观看国产| 一进一出抽搐动态| 亚洲国产看品久久| 亚洲色图综合在线观看| 日韩欧美免费精品| 亚洲欧美激情在线| 老司机午夜十八禁免费视频| 热99久久久久精品小说推荐| 99国产精品99久久久久| 久久免费观看电影| 91字幕亚洲| 欧美日韩亚洲国产一区二区在线观看 | 丰满人妻熟妇乱又伦精品不卡| 国产野战对白在线观看| 极品人妻少妇av视频| 国产在线一区二区三区精| 国产麻豆69| 久久精品熟女亚洲av麻豆精品| 中国国产av一级| 一个人免费在线观看的高清视频 | 狠狠婷婷综合久久久久久88av| 国产视频一区二区在线看| 亚洲国产av影院在线观看| 国产精品久久久人人做人人爽| 99国产精品99久久久久| 久久天躁狠狠躁夜夜2o2o| 亚洲 国产 在线| 国产精品 欧美亚洲| 国产三级黄色录像| 亚洲欧美精品综合一区二区三区| 久久这里只有精品19| 他把我摸到了高潮在线观看 | av在线老鸭窝| 午夜免费鲁丝| 纵有疾风起免费观看全集完整版| 999久久久国产精品视频| 午夜成年电影在线免费观看| 国产成人精品在线电影| 国产视频一区二区在线看| 亚洲av男天堂| 欧美黑人欧美精品刺激| 日日爽夜夜爽网站| 国产精品一区二区在线观看99| 我要看黄色一级片免费的| 亚洲精品自拍成人| 涩涩av久久男人的天堂| 精品一区二区三区av网在线观看 | tocl精华| 美女中出高潮动态图| 国产1区2区3区精品| 亚洲激情五月婷婷啪啪| 97在线人人人人妻| 久久综合国产亚洲精品| 久久99一区二区三区| 三上悠亚av全集在线观看| 美女扒开内裤让男人捅视频| 久久精品人人爽人人爽视色| 又紧又爽又黄一区二区| 成人国语在线视频| www.999成人在线观看| 亚洲国产精品成人久久小说| 热99re8久久精品国产| 1024香蕉在线观看| 免费在线观看视频国产中文字幕亚洲 | 手机成人av网站| 精品国产一区二区久久| 欧美激情 高清一区二区三区| 两性午夜刺激爽爽歪歪视频在线观看 | 一边摸一边做爽爽视频免费| 国产精品亚洲av一区麻豆| 精品少妇黑人巨大在线播放| a 毛片基地| 免费不卡黄色视频| 91精品三级在线观看| 欧美日韩精品网址| 操美女的视频在线观看| 十八禁网站免费在线| 97人妻天天添夜夜摸| 日本撒尿小便嘘嘘汇集6| 高清视频免费观看一区二区| 美女国产高潮福利片在线看| 国产一区二区在线观看av| 亚洲男人天堂网一区| 色老头精品视频在线观看| 亚洲国产av新网站| 亚洲成av片中文字幕在线观看| 国产精品久久久久久精品电影小说| 精品久久久久久电影网| 国产精品二区激情视频| 老熟妇乱子伦视频在线观看 | 亚洲精品国产av蜜桃| 久久ye,这里只有精品| 日本a在线网址| 欧美+亚洲+日韩+国产| 两个人看的免费小视频| 欧美在线一区亚洲| 新久久久久国产一级毛片| 国产99久久九九免费精品| 亚洲精品自拍成人| 午夜影院在线不卡| 国产成人a∨麻豆精品| 国产精品偷伦视频观看了| 国产淫语在线视频| 天天躁日日躁夜夜躁夜夜| 久久久国产成人免费| 久久久欧美国产精品| 超碰成人久久| 久久久精品国产亚洲av高清涩受| 深夜精品福利|