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

        ?

        等效線化方法分析亞音速壁板非線性極限環(huán)顫振

        2015-01-07 08:42:32唐懷平楊翊仁
        振動(dòng)工程學(xué)報(bào) 2015年5期
        關(guān)鍵詞:分析質(zhì)量系統(tǒng)

        唐懷平,楊翊仁

        (西南交通大學(xué)力學(xué)與工程學(xué)院,四川成都610031)

        等效線化方法分析亞音速壁板非線性極限環(huán)顫振

        唐懷平,楊翊仁

        (西南交通大學(xué)力學(xué)與工程學(xué)院,四川成都610031)

        研究了受集中質(zhì)量與非線性運(yùn)動(dòng)約束聯(lián)合作用下的二維亞音速壁板的極限環(huán)顫振問題。采用Galerkin方法將非線性壁板運(yùn)動(dòng)方程離散為常微分方程組。分析了集中質(zhì)量大小及其位置對(duì)壁板系統(tǒng)失穩(wěn)特性的影響;采用等效線化方法研究了系統(tǒng)的分叉特性及極限環(huán)顫振穩(wěn)定性。結(jié)果表明:系統(tǒng)會(huì)產(chǎn)生顫振失穩(wěn),質(zhì)量塊的大小及其位置對(duì)顫振臨界速度有著重要的影響;系統(tǒng)會(huì)經(jīng)歷超臨界的 Hopf分叉而處于穩(wěn)定的極限環(huán)運(yùn)動(dòng);等效線化方法可在一定范圍內(nèi)較為精確地對(duì)極限環(huán)穩(wěn)定性及其幅值進(jìn)行判定。

        壁板;亞音速流;等效線化;Hope分叉;極限環(huán)顫振

        引 言

        近些年,隨著高速鐵路運(yùn)行速度的不斷提升,列車的氣動(dòng)彈性問題也越來越突顯,并已經(jīng)成為高速列車中亟待解決 的關(guān)鍵基 礎(chǔ)問題 之 一[1-2]。由 于高速列車采用流線型的設(shè)計(jì),因此列車車體中存在著大量的蒙皮等壁板結(jié)構(gòu)。這些壁板結(jié)構(gòu)在列車高速運(yùn)行時(shí)會(huì)產(chǎn)生明顯的振動(dòng),如武廣客運(yùn)專線試驗(yàn)中,當(dāng)列車運(yùn)行速度達(dá)到350 km/h時(shí),列車車身蒙皮和車窗的振動(dòng)非常顯著,并會(huì)產(chǎn)生很大的輻射噪聲。因此有必要對(duì)這類特殊壁板的氣動(dòng)彈性問題進(jìn)行相關(guān)研究。

        就現(xiàn)有高速列車的運(yùn)行速度而言,其基本上屬于低亞音速范圍(馬赫數(shù)約為0.3)。針對(duì)亞音速氣流中壁板的氣動(dòng)穩(wěn)定性問題,文獻(xiàn)[2]曾指出兩端簡(jiǎn)支的壁板會(huì)在馬赫數(shù)為0.125時(shí)發(fā)生顫振失穩(wěn),并得到了風(fēng)洞吹風(fēng)試驗(yàn)的證實(shí)。而文[3-6]卻指出,兩端固定支撐的壁板在亞音速氣流中并不會(huì)發(fā)生顫振失穩(wěn),而僅會(huì)出現(xiàn)發(fā)散失穩(wěn);而一端固支一端自由的壁板卻會(huì)出現(xiàn)顫振失穩(wěn)。除了壁板的氣動(dòng)失穩(wěn)問題,亞音速壁板結(jié)構(gòu)在失穩(wěn)后的復(fù)雜非線性運(yùn)動(dòng)特性也是學(xué)者們關(guān)注的焦點(diǎn)?,F(xiàn)有的研究較為廣泛地考慮壁板大變形而產(chǎn)生的幾何非線性對(duì)系統(tǒng)失穩(wěn)特性的影響。文[7-9]均針對(duì)該非線性作用下的壁板的穩(wěn)定性及極限環(huán)響應(yīng)進(jìn)行了研究。事實(shí)上,由于生產(chǎn)、安裝過程中的誤差,壁板結(jié)構(gòu)常常還會(huì)受到其他結(jié)構(gòu)非線性因素的作用。這些非線性因素會(huì)約束壁板結(jié)構(gòu)的位移并導(dǎo)致壁板呈現(xiàn)出復(fù)雜動(dòng)力學(xué)特性。Li等[10-11]考慮 支撐松動(dòng)產(chǎn)生的接觸非線 性因素,研究了亞音速壁板的極限環(huán)顫振及混沌響應(yīng)。相比于幾何非線性而言,針對(duì)壁板在位移約束下的研究還比較欠缺。另外,上述研究主要是以理論模型分析為主,均未涉及到實(shí)際的風(fēng)洞模型。事實(shí)上,對(duì)于實(shí)際風(fēng)洞實(shí)驗(yàn)中的壁板而言,不可避免地需要對(duì)壁板結(jié)構(gòu)施加某些必要的集中質(zhì)量,以滿足模型設(shè)計(jì)要求及數(shù)據(jù)測(cè)試要求,例如在壁板關(guān)鍵位置安裝有一定質(zhì)量的傳感器及某些必要的實(shí)驗(yàn)掛件等。這些額外的重量對(duì)壁板的動(dòng)力學(xué)行為,尤其是對(duì)非線性極限環(huán)運(yùn)動(dòng)有何影響,也是需要關(guān)注的一個(gè)重要問題。

        因此,本文綜合考慮非線性約束及集中質(zhì)量?jī)蓚€(gè)因素的影響,對(duì)壁板的穩(wěn)定性及極限環(huán)響應(yīng)進(jìn)行分析。文中采用Galerkin方法對(duì)非線性亞音速黏彈性壁板的運(yùn)動(dòng)方程進(jìn)行離散;采用等效線化方法研究壁板的非線性運(yùn)動(dòng)特性。著重考察集中質(zhì)量對(duì)系統(tǒng)穩(wěn)定性及極限環(huán)運(yùn)動(dòng)的影響,并采用數(shù)值方法進(jìn)行積分驗(yàn)證。

        1 壁板顫振模型

        考慮一端固支一端受位移約束的懸臂二維黏彈性壁板,如圖1所示。壁板長(zhǎng)度為l,厚度為h,且h?l,壁板單位長(zhǎng)度的質(zhì)量為ρs。壁板上表面作用有沿x方向的亞音速不可壓縮氣流,來流速度為U∞,空氣密度為ρ∞。壁板在lm處作用有一質(zhì)量為m的集中質(zhì)量塊;壁板在端部受到的非線性運(yùn)動(dòng)約束fnon可以表示為

        式中w為壁板的橫向振動(dòng)位移,K1和K3為非線性運(yùn)動(dòng)約束的控制參數(shù)。

        圖1 亞音速懸臂壁板的幾何模型Fig.1 Schematic diagram of a cantilevered plate in subsonic flow

        利用Hamilton原理可得壁板的橫向振動(dòng)方程[10-11]

        式中D=為板的彎曲剛度,E為板的彈性模量,gs為黏性阻尼系數(shù),ν為泊松比,ΔP為作用在壁板上表面的氣動(dòng)壓力載荷。

        而壁板的邊界條件為

        由文[3-4,10]可知作用在壁板單側(cè)的氣動(dòng)力

        引入如下無量綱參數(shù)

        可得壁板的無量綱運(yùn)動(dòng)方程

        在下面的計(jì)算中,選取如下的基本參數(shù):E= 70 GPa,ρs=2 750 kg/m3,ρ∞=1.25 kg/m3,l= 1.0 m,h=2.0 mm,υ=0.3,gs=0.000 5。

        文獻(xiàn)[12]的研究表明采用系統(tǒng)前兩階模態(tài)可以得到較為滿意的定性和定量的分析結(jié)果。本文的目的在于定性分析和展示系統(tǒng)所蘊(yùn)含的典型的非線性特性,因此選取懸臂梁的前兩階模態(tài)對(duì)方程(6)進(jìn)行離 散[12],即

        采用Galerkin方法對(duì)式(6)進(jìn)行離散化,可得

        引入如下變換

        式(8)變?yōu)?/p>

        在選定的基本參數(shù)情況下,式(8)的各系數(shù)為:

        式(10)的各系數(shù)為

        2 顫振邊界分析

        首先,考察集中質(zhì)量及運(yùn)動(dòng)約束剛度對(duì)系統(tǒng)顫振穩(wěn)定性的影響。為計(jì)算系統(tǒng)的臨界顫振速度,將式(10)寫作狀態(tài)空間形式

        采用數(shù)值方法計(jì)算式(11)在零平衡點(diǎn)的雅克比矩陣的特征值,并依據(jù)特征值得特性對(duì)系統(tǒng)的顫振臨界速度進(jìn)行判定。

        圖2給出了當(dāng)附加質(zhì)量塊位于端部,即em=1時(shí),不同端部支撐剛度k1對(duì)應(yīng)的臨界顫振速度λf與集中質(zhì)量ˉm的變化關(guān)系。由圖2可知,隨著集中質(zhì)量的增加,顫振臨界速度呈現(xiàn)下降趨勢(shì),系統(tǒng)的顫振臨界速度隨著剛度的增加而增加。有趣的是,當(dāng)集中質(zhì)量較小時(shí),剛度對(duì)顫振臨界速度影響較為明顯;而當(dāng)集中質(zhì)量增加至0.5左右時(shí),不同剛度對(duì)應(yīng)的系統(tǒng)顫振臨界速度將趨于相同的值。

        圖3給出了當(dāng)無量綱集中質(zhì)量=0.1,不同剛度k1時(shí),系統(tǒng)顫振臨界臨界速度與集中質(zhì)量的位置ξm之間的變化關(guān)系。圖3中的點(diǎn)劃線對(duì)應(yīng)的是無集中質(zhì)量時(shí)系統(tǒng)的顫振臨界速度。由圖3可知:不同剛度時(shí)系統(tǒng)的顫振臨界速度隨位置的變化呈現(xiàn)相

        圖2 不同集中質(zhì)量時(shí)系統(tǒng)顫振邊界Fig.2 Flutter boundary for different

        圖3 不同集中質(zhì)量位置時(shí)系統(tǒng)顫振邊界Fig.3 Flutter boundary for different

        似性;系統(tǒng)的顫振臨界速度隨位置的變化呈現(xiàn)非線性關(guān)系。特別值得指出的是:不同剛度對(duì)應(yīng)的顫振邊界與點(diǎn)劃線的交點(diǎn)有著相同的橫坐標(biāo),即

        當(dāng)<<時(shí),集中質(zhì)量會(huì)提高系統(tǒng)的臨界顫振速度并增加系統(tǒng)的穩(wěn)定性;而當(dāng)<或ξm>ξmB時(shí),集中 質(zhì) 量卻 起 到了相 反 的作 用;而 在ξm=ξmA及ξm=ξmB兩點(diǎn)處,集中質(zhì)量對(duì)系統(tǒng)的顫振穩(wěn)定性沒有影響。

        圖4給出了=0.1,k1=5.0,λ=49.29時(shí),不同集中質(zhì)量位置對(duì)應(yīng)的系統(tǒng)線性的響應(yīng)相圖。從圖可知:即使在相同的動(dòng)壓下,由于不同集中質(zhì)量的放置位置,系統(tǒng)也會(huì)呈現(xiàn)出截然不同的穩(wěn)定特性。

        圖4 不同集中質(zhì)量位置時(shí)線性系統(tǒng)的運(yùn)動(dòng)響應(yīng)相圖Fig.4 Phase plots for differentξm

        3 極限環(huán)顫振分析及討論

        在上一節(jié)中,已經(jīng)分析了系統(tǒng)的顫振失穩(wěn)。下面就針對(duì)系統(tǒng)在顫振失穩(wěn)后可能呈現(xiàn)的非線性極限環(huán)運(yùn)動(dòng)進(jìn)行分析。采用等效線化方法來分析極限環(huán)顫振(LCO),已有許多重要的成果發(fā)表。文獻(xiàn)[11,13-14]均基于等效線化方法對(duì)機(jī)翼系統(tǒng)中的極限環(huán)顫振運(yùn)動(dòng)進(jìn)行了研究。這種方法為定性及定量分析壁板系統(tǒng)的極限環(huán)顫振提供了指導(dǎo)。因此本文也采用等效線化方法對(duì)極限環(huán)運(yùn)動(dòng)進(jìn)行分析。在下面的計(jì)算中進(jìn)一步給定如下參數(shù):k1=5,k3=2,ˉm= 0.1。

        針對(duì)本文中的立方非線性而言,依據(jù)等效線化方法,式(11)中的非線性可以表示為

        式中Keq=,A為系統(tǒng)(11)的極限環(huán)運(yùn)動(dòng)分支p1的幅值。

        將式(13)代入至方程(11)后,非線性系統(tǒng)(11)將變成線性系統(tǒng),其狀態(tài)方程可以寫作

        下面考察ξm=0.7和ξm=1.0兩個(gè)特殊位置的顫振邊界曲線,結(jié)果分別如圖5(a)和(b)所示。從圖5可知,隨著等效剛度的增加兩組顫振臨界速度都呈現(xiàn)單調(diào)增加趨勢(shì)。

        圖5 系統(tǒng)顫振速度隨等效剛度的變化關(guān)系Fig.5 Curves of critical flutter dynamics pressure vs.equivalent linearized stiffness

        下面以ξm=1.0為例說明如何依據(jù)等效剛度-幅值-顫振臨界動(dòng)壓三者的耦合圖(如圖6所示)對(duì)系統(tǒng)極限環(huán)的穩(wěn)定性進(jìn)行判定。首先,假設(shè)動(dòng)壓λt對(duì)應(yīng)的極限環(huán)幅值為As,依據(jù)耦合圖可知,當(dāng)其有一個(gè)正的增量ΔAs時(shí),等效剛度及顫振臨界速度都將會(huì)增大,系統(tǒng)將呈現(xiàn)穩(wěn)定的收斂到幅值為As的運(yùn)動(dòng)。而當(dāng)幅值有負(fù)的增量——ΔAs時(shí),由于等效剛度及顫振臨界速度都將降低,系統(tǒng)也將呈現(xiàn)收斂于幅值為As的運(yùn)動(dòng)。因此,無論對(duì)于正的或者負(fù)的幅值增量,系統(tǒng)都將最終穩(wěn)定于幅值為As的穩(wěn)定極限環(huán)運(yùn)動(dòng)。因此,系統(tǒng)在λ=λf處產(chǎn)生超臨界的Hopf分叉。同理,針對(duì)ξm=0.7這一工況,也會(huì)得到相同的結(jié)論。

        圖6 極限環(huán)穩(wěn)定性判定耦合圖Fig.6 A coupled scheme for the stability of LCOs

        圖7(a)及(b)分別給出了當(dāng)動(dòng)壓λ=40(小于顫振速度λf=41.5)時(shí)及λ=43時(shí),系統(tǒng)的運(yùn)動(dòng)相圖。由圖7(a)及(b)可知,系統(tǒng)分別處于收斂運(yùn)動(dòng)和穩(wěn)定極限環(huán)運(yùn)動(dòng),這也與理論預(yù)測(cè)的一致。圖8給出了ξm=0.7時(shí)系統(tǒng)穩(wěn)定極限環(huán)幅值隨動(dòng)壓的變化曲線圖,圖8中圓圈表示數(shù)值模擬的結(jié)果,而實(shí)線代表等效線化的分析結(jié)果。由圖8可知:當(dāng)動(dòng)壓較小時(shí)(Ⅰ區(qū)域),等效線化的分析結(jié)果與數(shù)值模擬結(jié)果吻合較好;但當(dāng)動(dòng)壓較大(Ⅱ區(qū)域)時(shí),兩者的相差較大。這主要是由于在Ⅱ區(qū)域內(nèi)系統(tǒng)已經(jīng)處于非周期運(yùn)動(dòng),如混沌等復(fù)雜運(yùn)動(dòng)(如圖9所示)。而等效線化分析方法卻只能對(duì)周期運(yùn)動(dòng)進(jìn)行分析。對(duì)于混沌等運(yùn)動(dòng)的研究可參見文獻(xiàn)[10],而本文暫不涉及。

        圖7 不同動(dòng)壓時(shí)系統(tǒng)的運(yùn)動(dòng)相圖Fig.7 Phase plots for differentλ

        圖9 λ=55時(shí)系統(tǒng)運(yùn)動(dòng)相圖Fig.9 Phase plots forλ=55

        圖10(a)及(b)分別給出了ξm=0.7及ξm=1.0時(shí),系統(tǒng)處于穩(wěn)定的極限環(huán)運(yùn)動(dòng)時(shí)壁板的振動(dòng)形態(tài)圖。從圖10可知看出,雖然在兩種不同位置時(shí)壁板均呈現(xiàn)極限環(huán)運(yùn)動(dòng),但其運(yùn)動(dòng)的形態(tài)已經(jīng)有了很大的不同。相比于圖10(b)而言,圖10(a)中在ξ=0.7附 近 顯示 出 了明 顯 的節(jié) 點(diǎn)[10-11]。

        4 結(jié)束語

        本文利用等效線化方法分析了受集中質(zhì)量和位移約束作用的亞音速黏彈性壁板的極限環(huán)顫振運(yùn)動(dòng)。結(jié)果表明:系統(tǒng)的顫振臨界速度隨著集中質(zhì)量位置及位移約束剛度的增大而分別呈現(xiàn)減小和增大的趨勢(shì)(如圖2,3所示);不同的集中質(zhì)量位置會(huì)導(dǎo)致不同性質(zhì)的運(yùn)動(dòng)響應(yīng)(如圖4所示);集中質(zhì)量的位置對(duì)系統(tǒng)的極限環(huán)振動(dòng)的節(jié)點(diǎn)及其位置有著顯著的影響(如圖10所示);系統(tǒng)經(jīng)歷超臨界的Hopf分叉而處于穩(wěn)定的極限環(huán)運(yùn)動(dòng)(如圖7(b)所示);等效線化方法可在一定范圍內(nèi)(如在圖8中的Ⅰ區(qū)域)對(duì)極限環(huán)運(yùn)動(dòng)穩(wěn)定性(如圖6所示)及極限環(huán)的運(yùn)動(dòng)幅值進(jìn)行比較精確地分析。但需要特別指出的是,當(dāng)壁板處于幅值相對(duì)較大的極限環(huán)運(yùn)動(dòng)時(shí),壁板將不可避免地受到其自身大變形幾何非線性的約束并影響到其非線性運(yùn)動(dòng)特性,這也是后續(xù)理論及實(shí)驗(yàn)工作研究的重點(diǎn)。

        [1] Schetz J A.Aerodynamics of high-speed trains[J]. Ann.Rev.Fluid.Mech.,2001,33:371—414.

        [2] 李鵬,楊翊仁,魯麗.氣動(dòng)力作用下高速車輛橫向穩(wěn)定性分析[J].振動(dòng)與沖擊,2010,29(11):135—138. Li Peng,Yang Yi-ren,Lu Li.Lateral stability of a high-speed train under aerodyanmic force[J].Journal of Vibration and Shock,2010,29(11):135—138.

        [3] Dugundji J,Dowell E H,Perkin B.Subsonic flutter of panels on continuous elastic foundations[J].AIAA Journal,1963,5:1 146—1 154.

        [4] Kornecki A,Dowell E H,O’Brien J.On the aeroelastic instability of two-dimensional panels in uniform incompressible flow[J].Journal of Sound and Vibration,1974,47:163—178.

        [5] 李鵬,楊翊仁,魯麗.外激勵(lì)作用下亞音速二維壁板的分叉及響應(yīng)研究[J].力學(xué)學(xué)報(bào),2011,43(4):746—754. Li Peng,Yang Yi-ren,Lu Li.Bifurcation and responses analysis of two-dimensioan panel with external excitation in subsonic flow[J].Chinese Jouranl of Theoretical and Applied Mechancis,2011,43(4):746—754.

        [6] 李鵬,楊翊仁,魯麗.微分求積分析二維亞音速壁板的失穩(wěn)問題[J].動(dòng)力學(xué)與控制學(xué)報(bào),2012,10(1):11—14. Li Peng,Yang Yi-ren,Lu Li.Instability analysis of two-dimensional thin panel in subsonic flow with differential quadrature method[J].Jouranl of Dyanmics and Control,2012,10(1):11—14.

        [7] YAO G,LI F M.Chaotic motion of a composite lami-nated plate with geometric nonlinearity in subsonic flow[J].International Journal of Non-Linear Mechancis,2013,50:81—90.

        [8] Ellen C H.The non-linear stability of panels in incompressible flow[J].Journal of Sound and Vibration,1977,54:117—121.

        [9] Matsuzaki Y.Reexamination of stability of a two-dimensional finite panel exposed to an incompressible flow[J].ASME Journal of Applied Mechanics,1981,48:472—478.

        [10]Li P,Yang Y R.On the stability and chaos of a plate with motion constraints subjected to subsonic flow [J].International Journal of Non-Linear Mechancis,2014,59:28—36.

        [11]Li P,Yang Y R.Hopf and two-multiple semi-stable limit cycle bifurcations of a restrained plate subjected to subsonic flow[J].Journal of Sound and Vibration,2015,335:286—303.

        [12]Jin J D,Stability and chaotic motions of a restrained pipe conveying fluid[J].Journal of Sound and Vibration,1997,208:427—439.

        [13]Yang Z C,Zhao L C.Analysis of limit cycle flutter of an airfoil in incompressible flow[J].Journal of Sound and Vibration,1988,123(1):1—13.

        [14]Yang Y R.KBM method of analyzing limit cycle flutter of a wing with an external store and comparison with a wind-tunnel test[J].Journal of Sound and Vibration,1995,187(2):163—178.

        圖8 系統(tǒng)極限環(huán)幅值隨動(dòng)壓的變化關(guān)系
        Fig.8 The curve of the amplitudes of LCOs vs.dynamic pressure

        Analysis of nonlinear limit cycle flutter of a subsonic plate based on equivalent linearized method

        TANG Huai-pin,YANG Yi-ren
        (School of Mechanics and Engineering,Southwest Jiaotong University,Chengdu 610031,China)

        The limit cycle flutter of a subsonic plate with a concentrated mass and subjected to nonlinear motion constraints is addressed in this paper.The Galerkin method is used to transfer the partial differential equation of motion of the plate to a set of ordinary differential equations.The theoretical analysis of the bifurcations and the stabilities of the limit cycles are conducted with the help of the equivalent linearized method.Results show that the system undergoes flutter instability;the mass and its location has significant effect on the flutter critical dynamic pressure;the system undergoes stable limit cycle oscillations (LCO)due to the supercritical Hopf bifurcation;the results obtain by equivalent linearized method are in good agreement with the numerical ones.

        plate;subsonic flow;equivalent linearized method;Hope bifurcation;limit cycle flutter

        U270.1+1

        A

        1004-4523(2015)05-0748-06

        10.16385/j.cnki.issn.1004-4523.2015.05.009

        唐懷平(1967—),男,博士研究生,副教授。電話:(028)87600797;E-mail:thp-vib@163.com

        2014-11-28;

        :2015-04-28

        國(guó)家自然科學(xué)基金資助項(xiàng)目(11302183;11072204;11102170;11102172);中央高?;究蒲袠I(yè)務(wù)費(fèi)專項(xiàng)資金資助項(xiàng)目(2013TD0004)

        猜你喜歡
        分析質(zhì)量系統(tǒng)
        Smartflower POP 一體式光伏系統(tǒng)
        “質(zhì)量”知識(shí)鞏固
        WJ-700無人機(jī)系統(tǒng)
        隱蔽失效適航要求符合性驗(yàn)證分析
        ZC系列無人機(jī)遙感系統(tǒng)
        質(zhì)量守恒定律考什么
        做夢(mèng)導(dǎo)致睡眠質(zhì)量差嗎
        電力系統(tǒng)不平衡分析
        電子制作(2018年18期)2018-11-14 01:48:24
        連通與提升系統(tǒng)的最后一塊拼圖 Audiolab 傲立 M-DAC mini
        電力系統(tǒng)及其自動(dòng)化發(fā)展趨勢(shì)分析
        手机在线免费观看av不卡网站| 综合无码一区二区三区四区五区| 国产亚洲精品综合在线网址| 超碰青青草手机在线免费观看 | 蜜臀av中文人妻系列| 精品乱色一区二区中文字幕| 国产乱人偷精品人妻a片| 亚洲av无码一区二区二三区 | 亚洲综合av大全色婷婷| 又黄又硬又湿又刺激视频免费| 欧美精品在线一区| 亚洲一区二区三区乱码在线| 免费亚洲老熟熟女熟女熟女| 久久综合狠狠色综合伊人| 日本高清不卡二区| 国产精品国产三级在线专区| 日韩中文字幕有码午夜美女| 极品av麻豆国产在线观看| 无码免费人妻超级碰碰碰碰| 一本久道在线视频播放| 激情人妻另类人妻伦| 好大好深好猛好爽视频免费| 国产午夜激情视频自拍| 快射视频网站在线观看| 草草地址线路①屁屁影院成人| 97免费人妻在线视频| 国产颜射视频在线播放| 在线观看国产视频你懂得| 小鲜肉自慰网站| 在线观看网址你懂的| 国产午夜激情视频在线看| 亚洲一区二区三区四区五区黄| 无码h黄动漫在线播放网站| 漂亮的小少妇诱惑内射系列| 国产成人大片在线播放| 男女后进式猛烈xx00动态图片| 2022Av天堂在线无码| 国成成人av一区二区三区| 一色桃子中文字幕人妻熟女作品| 日本a级特黄特黄刺激大片| 成人国产自拍在线播放|