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

摩擦力作用下電液伺服系統(tǒng)非線性動力學行為

2015-12-19 00:55:32朱勇姜萬錄鄭直
北京航空航天大學學報 2015年1期
關(guān)鍵詞:振動系統(tǒng)

朱勇,姜萬錄*,鄭直

(1.燕山大學 河北省重型機械流體動力傳輸與控制重點實驗室,秦皇島066004;2.燕山大學 先進鍛壓成形技術(shù)與科學教育部重點實驗室,秦皇島066004)

伺服系統(tǒng)中,實現(xiàn)運動轉(zhuǎn)換與傳遞的活塞和缸壁之間普遍存在著摩擦力,它在低速時表現(xiàn)出的強非線性作用,會導致爬行和振蕩現(xiàn)象的出現(xiàn),使得系統(tǒng)的控制效果劣化[1].因此,摩擦力對伺服系統(tǒng)動態(tài)特征的影響不容忽視.

目前,針對非線性摩擦力作用下電液伺服系統(tǒng)動態(tài)特征的研究尚不多見.對液壓系統(tǒng)動態(tài)特性的研究一般采用系統(tǒng)建模和數(shù)值仿真手段[2-3],所依據(jù)的理論多是經(jīng)典控制理論和線性動力學理論[4-5],較少運用非線性動力學理論[6-7],且多數(shù)在建模時對非線性因素進行線性化處理[8-9],研究結(jié)論與實際情況有較大差異,很難解釋實際動態(tài)測試中出現(xiàn)的時域波形復雜、頻域尖峰繁雜等異常現(xiàn)象.然而,電液伺服系統(tǒng)是非線性動力學系統(tǒng),執(zhí)行機構(gòu)運行過程是典型的非線性過程.故應引入非線性動力學理論和方法進行分析研究.

本文重點探究非線性摩擦力對電液伺服系統(tǒng)動態(tài)特征的影響規(guī)律.通過理論研究和數(shù)值試驗分析,深入揭示系統(tǒng)內(nèi)在的分岔現(xiàn)象及典型非線性動力學行為.用非線性動力學研究方法對實測的動態(tài)數(shù)據(jù)進行深入分析,以揭示非線性摩擦力引起的“極限環(huán)型振蕩”現(xiàn)象.旨在揭示伺服系統(tǒng)非線性振動的機理及誘因,使綜合分析系統(tǒng)的動態(tài)特征變得更接近實際.

1 執(zhí)行機構(gòu)的動力學模型

電液伺服系統(tǒng)的執(zhí)行機構(gòu)為伺服液壓缸,本文以雙作用單活塞桿液壓缸為例進行分析,其工作原理如圖l所示,其動力學方程為

式中,m為活塞及慣性負載的折合質(zhì)量;x為活塞位移;Fc為黏性力;Fs為彈性力;Ff為摩擦力;FL為負載力;p1,p2分別為無桿腔和有桿腔壓力;A1,A2分別為無桿腔和有桿腔活塞有效作用面積.

圖1 雙作用單活塞桿伺服液壓缸工作原理Fig.1 Working principle of double-acting single piston rod servo cylinder

2 非線性摩擦力

圖2所示為體現(xiàn)摩擦力與速度關(guān)系的Stribeck曲線,數(shù)學模型可由式(2)和式(3)描述[10].

圖2 摩擦力與速度的關(guān)系曲線Fig.2 Relation curve of friction and velocity

式中,F(xiàn)e為外力;Fj為最大靜摩擦力;Fk為庫侖摩擦力;vs為Stribeck速度;B為黏滯摩擦系數(shù);τ為經(jīng)驗參數(shù),一般在0.5~2.0之間取值.

摩擦力的作用效果隨工作點在Stribeck曲線上所處區(qū)段不同而異[1]:①當工作點位于區(qū)域Ⅰ時,Stribeck曲線斜率為很大的正數(shù),其作用效果使系統(tǒng)阻尼瞬間陡增,遏制系統(tǒng)的動態(tài)響應的快速性.②當工作點位于區(qū)域Ⅳ時,摩擦力特性為線性的,其作用效果相當于增加了系統(tǒng)的阻尼.③當工作點位于區(qū)域Ⅱ或區(qū)域Ⅲ時,摩擦力呈現(xiàn)非線性時變特性,當處于負阻尼狀態(tài)時會產(chǎn)生具有極限環(huán)的自激振動.

3 摩擦力非線性動態(tài)特征

3.1 非線性動力學模型

為了集中研究非線性摩擦力對系統(tǒng)動態(tài)特征的影響,暫不考慮液壓彈簧力、黏性力等非線性因素.當工作點位于區(qū)域Ⅱ或Ⅲ時,有[1]

則式(1)在工作點x附近的特性可表達為

式中,c0為結(jié)構(gòu)阻尼系數(shù);k為線性彈簧剛度.

移項整理,得

由于油源壓力脈動、閥口流量-壓力非線性等因素的影響,進入液壓缸的油壓有微觀波動,基本服從簡諧振動規(guī)律,故式(7)右邊的輸入項可近似表示為 Fsin(ωt+φ0),是系統(tǒng)的激振源[11].

據(jù)上述分析,式(7)可化為

式中,F(xiàn)為激振力;ω為激振角頻率;φ0為激振力的初相角.

進一步整理,得

由Van Der Pol方程的結(jié)構(gòu)形式可知,式(10)是受迫振動Van Der Pol方程.該方程為研究伺服系統(tǒng)非線性摩擦力的動態(tài)特征提供了結(jié)構(gòu)模型.

3.2 解析解及分析

把摩擦力的非線性動態(tài)特征歸結(jié)為Van Der Pol方程,就可以借助該方程的特性來揭示系統(tǒng)內(nèi)在的基本規(guī)律.用多尺度法可求得式(10)的近似解,其自由振動的振幅如式(11)所示[12]:a0為積分常數(shù).分析式(11)發(fā)現(xiàn):隨著t→∞,當η<0(即)時,a→0,自由振動趨于衰減,式(10)受激勵后的穩(wěn)態(tài)運動是頻率為ω的受迫振動;當η>0(即)時,,穩(wěn)態(tài)運動中除包含頻率為ω的受迫振動外,還含有頻率為ωn的自由振動,通常ω與ωn不可共約,此穩(wěn)態(tài)運動是非周期的.可見,大激勵引發(fā)自由振動衰減,小激勵引發(fā)非周期運動.

3.3 穩(wěn)定性分析

當F0=0時,設(shè)狀態(tài)變量y1=y,y2=,則式(10)的狀態(tài)方程為

在原點附近加以線性化,得

特征方程為

可解出特征根:

分析式(16),當 μ>0時,特征根的實部Re(λi)>0,由系統(tǒng)穩(wěn)定性與特征方程根的關(guān)系可知,工作點是不穩(wěn)定的平衡點,從工作點近旁出發(fā)的軌線是發(fā)散的,但由于穩(wěn)定極限環(huán)的存在,限制了軌線的無限發(fā)散;當μ<0時,特征根的實部Re(λi)<0,工作點是穩(wěn)定的平衡點,在極限環(huán)內(nèi)部的軌線收斂于平衡點[13].

4 數(shù)值試驗

以系統(tǒng)方程式(10)的具體算例式(17)為例,進行數(shù)值試驗研究,以探索系統(tǒng)阻尼系數(shù)μ和外加激振力F0對系統(tǒng)動態(tài)特征的影響.

4.1 分岔特性研究

F0取不同值時,以μ為分岔參數(shù)的分岔圖如圖3所示.圖中μ為單位質(zhì)量上的系統(tǒng)阻尼系數(shù),y為振動位移.

圖3 分岔參數(shù)為μ的非線性方程分岔圖Fig.3 Bifurcation diagrams with μ as parameter

由圖3可知,當參數(shù)μ和F0取不同值時系統(tǒng)發(fā)生了不同程度的分岔現(xiàn)象:①系統(tǒng)方程存在單解、多解和無窮多解,反映在分岔圖上表現(xiàn)為單值曲線、多值曲線和涂黑區(qū)等不同的區(qū)段,分別對應于單周期、多周期和混沌等不同的運動狀態(tài).②隨著參數(shù)的變化,系統(tǒng)會發(fā)生運動狀態(tài)突然變化的動態(tài)分岔現(xiàn)象.③隨著阻尼系數(shù)μ的增大,混沌區(qū)會突然消失,系統(tǒng)出現(xiàn)周期3、周期5或周期7的穩(wěn)定狀態(tài),接著再次進入混沌運動狀態(tài).

4.2 運動形態(tài)仿真

為了形象地體現(xiàn)系統(tǒng)在不同參數(shù)下的運動形態(tài),在MATLAB中建立仿真模型,對其典型的非線性動力學行為進行仿真.仿真中采用Runge-Kutta算法,采樣頻率100 Hz,遠大于外控力頻率fP=ω/2π =0.16 Hz,終了時間1000 s.

當 μ=2 N·s/(mm·kg),F(xiàn)0=2 N/kg時,仿真結(jié)果如圖4所示.由圖4可知,時間歷程呈周期重復;功率譜在基頻fP及其倍頻處出現(xiàn)尖峰,所有譜峰對應的頻率可共約;相軌跡在有限的區(qū)域內(nèi)重復,呈一封閉曲線,即有極限環(huán)存在;Poincaré圖在一定的區(qū)域上只有1個孤立點存在.說明此時系統(tǒng)處于極限環(huán)型振蕩狀態(tài).

圖4 極限環(huán)型振蕩形態(tài)Fig.4 Forma of limit cycle oscillation

當μ=1.7 N·s/(mm·kg),F(xiàn)0=1.2 N/kg時,仿真結(jié)果如圖5所示.由圖5可知,時間歷程呈周期重復;功率譜在分頻fP/3及其倍頻處存在尖峰,所有譜峰對應的頻率可共約;相軌跡在有限的區(qū)域內(nèi)重復,呈封閉曲線;Poincaré圖在一定的區(qū)域上有3個孤立點存在.說明此時系統(tǒng)處于3倍周期運動狀態(tài).

圖5 3倍周期運動形態(tài)Fig.5 Forma of triple periodic motion

當 μ=0.8 N·s/(mm·kg),F(xiàn)0=0.8 N/kg時,仿真結(jié)果如圖6所示.由圖6可知,時間歷程無規(guī)律;功率譜出現(xiàn)噪聲背景和寬峰;相軌跡在有限的區(qū)域內(nèi)不重復;Poincaré圖有無限個孤立點存在,且分布在有限的區(qū)域上.說明此時系統(tǒng)處于混沌運動狀態(tài).

由以上數(shù)值試驗分析可知,當系統(tǒng)阻尼系數(shù)μ和外加激振力F0取不同值時,系統(tǒng)在運行過程中蘊含豐富的非線性動力學行為.系統(tǒng)可能做極限環(huán)型振蕩、倍周期運動,進而通向混沌運動.

圖6 混沌運動形態(tài)Fig.6 Forma of chaotic motion

5 電液伺服系統(tǒng)動態(tài)試驗

本節(jié)將利用非線性動力學研究方法對實測的電液伺服系統(tǒng)的動態(tài)數(shù)據(jù)進行深入分析,以揭示非線性摩擦力引起的極限環(huán)型振蕩現(xiàn)象.

5.1 試驗系統(tǒng)組成

按圖7所示的系統(tǒng)原理圖搭建振動測試試驗臺[7].該試驗系統(tǒng)可實現(xiàn)對電液伺服系統(tǒng)在不同供油壓力和負載壓力下的狀態(tài)數(shù)據(jù)進行采集.系統(tǒng)中通過調(diào)節(jié)溢流閥閥口開度大小來改變系統(tǒng)供油壓力;通過調(diào)節(jié)節(jié)流閥閥口開度大小來改變負載壓力,以實現(xiàn)系統(tǒng)外加阻尼大小的調(diào)整;用精密壓力表對系統(tǒng)進、回油路壓力進行監(jiān)測;用振動傳感器對系統(tǒng)執(zhí)行機構(gòu)軸向振動信號進行監(jiān)測;用位移傳感器對執(zhí)行機構(gòu)實時位置進行監(jiān)視;用數(shù)據(jù)采集卡采集傳感器輸出信號,并送基于虛擬儀器的計算機系統(tǒng)進行分析處理.

圖7 狀態(tài)監(jiān)測及試驗系統(tǒng)原理圖Fig.7 Schematic of condition monitoring and experiment

5.2 振動信號的采集及處理

5.2.1 振動信號的采集

按表1設(shè)置的工況對系統(tǒng)工作狀態(tài)進行動態(tài)測試.其中,輸入信號為計算機控制系統(tǒng)給伺服放大器的電壓信號,以控制伺服閥的閥口開度大小.根據(jù)供油壓力和負載壓力大小將外加阻尼大小界定為4類:無、較小、適中、較大.輸入信號為0.2V時,調(diào)整主溢流閥及節(jié)流閥閥口開度,使系統(tǒng)分別在表1設(shè)定的12種工況下運行.同時用振動加速度傳感器對液壓缸整個運行過程中的軸向振動信號進行采集,采樣頻率為10 kHz.

表1 不同供油壓力情況下執(zhí)行機構(gòu)工作狀態(tài)Table1 Actuator working condition under different supply pressure

5.2.2 振動信號的處理

采用圖8中的數(shù)據(jù)處理方案對采集的振動加速度信號進行預處理,并采用非線性動力學研究方法中的時間歷程、頻閃采樣、功率譜分析等有效方法對預處理數(shù)據(jù)進行分析研究[15].

圖8 數(shù)據(jù)處理方案Fig.8 Program of data processing

5.3 試驗結(jié)果分析

圖9為供油壓力8 MPa時采用圖8中的數(shù)據(jù)處理方案對采集的振動加速度信號進行處理所得到的振動位移信號時域波形.比較4種工況發(fā)現(xiàn),在整個運行過程中,振動幅值存在顯著波動,隨著活塞位移的變化而變化.由此可知,執(zhí)行機構(gòu)在運行過程中,系統(tǒng)動態(tài)性能隨活塞位移的變化而變化,其變化規(guī)律隨工況不同而存在明顯差異.

圖10為供油壓力8 MPa時4種工況的全程頻閃采樣圖,它的覆蓋區(qū)域與Van Der Pol極限環(huán)極為相似.工況1和2的輪廓邊界由許多離散點構(gòu)成,發(fā)生了較強烈的振蕩;而工況3和4的輪廓邊界較清晰,振蕩程度相對較弱,每種工況都有一個極限環(huán)存在,說明“極限環(huán)型振蕩”現(xiàn)象存在的普遍性.工況1的面積最大,工況2和3的次之,工況4的最小,說明極限環(huán)面積的大小由系統(tǒng)外加阻尼的大小來決定.

圖9 振動位移信號時域波形Fig.9 Waveforms of vibration displacement signals

圖11~圖14為供油壓力8MPa時4種工況的分段功率譜圖.根據(jù)執(zhí)行機構(gòu)運行總時間長度,將其分成等分的4段,分別定義為:始段、中前段、中后段、終段.比較4種工況發(fā)現(xiàn),振動能量值隨負載壓力的增大而逐漸降低,說明隨系統(tǒng)阻尼的增大,振動幅值被抑制.同一工況下,不同分段中相同的頻率點處有相應的譜峰存在,進一步說明了摩擦極限環(huán)型振蕩現(xiàn)象的存在.

圖10 全程頻閃采樣圖Fig.10 Diagrams of entire stroboscopic sampling

圖11 負載壓力0 MPa時分段功率譜圖Fig.11 Power spectrums at load pressure of 0 MPa

圖12 負載壓力2 MPa時分段功率譜圖Fig.12 Power spectrums at load pressure of 2 MPa

圖13 負載壓力4 MPa時分段功率譜圖Fig.13 Power spectrums at load pressure of 4 MPa

圖14 負載壓力6 MPa時分段功率譜圖Fig.14 Power spectrums at load pressure of 6 MPa

為了驗證上述所得結(jié)論的普遍性,采用與供油壓力為8 MPa時相同的數(shù)據(jù)處理方法,分別對供油壓力為6 MPa和4 MPa時的試驗數(shù)據(jù)進行了進一步的分析研究.通過比較分析,不同供油壓力下同樣可以得到非線性摩擦力引起的極限環(huán)型振蕩現(xiàn)象存在的普遍性的結(jié)論.

6 結(jié)論

本文以電液伺服系統(tǒng)為研究對象,根據(jù)非線性動力學原理,重點探究了摩擦力的非線性作用對系統(tǒng)動態(tài)特征的影響規(guī)律.通過理論研究、仿真分析及試驗驗證,得到了如下結(jié)論:

1)電液伺服系統(tǒng)執(zhí)行機構(gòu)在往復運動過程中,摩擦力變化規(guī)律服從Stribeck曲線,在速度趨近于零時取值具有不確定性,在低速運動區(qū)段呈現(xiàn)非線性規(guī)律,在高速運動區(qū)段呈現(xiàn)線性規(guī)律.

2)摩擦力的非線性作用隨工作點在Stribeck曲線上所處區(qū)段不同而異.在靜摩擦區(qū),其作用效果使系統(tǒng)阻尼瞬間陡增,遏制了系統(tǒng)動態(tài)響應的快速性;在流體動壓潤滑區(qū),其作用效果增加了系統(tǒng)的阻尼;在邊界潤滑區(qū)和混合潤滑區(qū),其動態(tài)特征可以用Van Der Pol方程來描述,其作用效果使系統(tǒng)阻尼減小,有可能使系統(tǒng)變成負阻尼系統(tǒng),產(chǎn)生極限環(huán)型振蕩.

3)系統(tǒng)阻尼系數(shù)和外加激振力的大小影響系統(tǒng)的運動狀態(tài).當二者參數(shù)取不同值時,系統(tǒng)可能做極限環(huán)型振蕩、倍周期運動,進而通向混沌運動.

4)非線性摩擦力引起的極限環(huán)型振蕩會使系統(tǒng)響應穩(wěn)定區(qū)域變得復雜,進而造成執(zhí)行機構(gòu)在工作過程中發(fā)生非線性振動及動態(tài)特性變得復雜和多變.因此在系統(tǒng)建模與動態(tài)特性研究時應該將摩擦力的非線性作用考慮在內(nèi),其對系統(tǒng)動態(tài)特征的影響不容忽視.

References)

[1] 王林鴻,吳波,杜潤生,等.液壓缸運動的非線性動態(tài)特征[J].機械工程學報,2007,43(12):12-19.Wang L H,Wu B,Du R S,et al.Nonlinear dynamic characteristics of moving hydraulic cylinder[J].Chinese Journal of Mechanical Engineering,2007,43(12):12-19(in Chinese).

[2] Dasgupta K,Murrenhoff H.Modelling and dynamics of a servovalve controlled hydraulic motor by bondgraph[J].Mechanism and Machine Theory,2011,46(7):1016-1035.

[3] Lan Z K,Su J,Xu G,et al.Study on dynamical simulation of railway vehicle bogie parameters test-bench electro-hydraulic servo system[J].Physics Procedia,2012,33:1663-1669.

[4] Milic V,Situm Z,Essert M.Robust H∞position control synthesis of an electro-hydraulic servo system[J].ISA Transactions,2010,49(4):535-542.

[5] 郭棟,付永領(lǐng),盧寧,等.自抗擾控制技術(shù)在電液力伺服系統(tǒng)中的應用[J].北京航空航天大學學報,2013,39(1):115-119.Guo D,F(xiàn)u Y L,Lu N,et al.Application of ADRC technology in electrohydraulic force servo system[J].Journal of Beijing University of Aeronautics and Astronautics,2013,39(1):115-119(in Chinese).

[6] Chen C T.Hybrid approach for dynamic model identification of an electro-hydraulic parallel platform[J].Nonlinear Dynamics,2012,67(1):695-711.

[7] 朱勇,姜萬錄,王夢,等.非線性時變力作用下液壓缸爬行機理與抑制方法研究[J].農(nóng)業(yè)機械學報,2014,45(3):305-313.Zhu Y,Jiang W L,Wang M,et al.Creeping mechanism and suppression methods of hydraulic cylinder under nonlinear time-varying force[J].Transactions of the Chinese Society for Agricultural Machinery,2014,45(3):305-313(in Chinese).

[8] Tang R,Zhang Q.Dynamic sliding mode control scheme for electro-hydraulic position servo system[J].Procedia Engineering,2011,24:28-32.

[9] Seo J,Venugopal R,Kenné J P.Feedback linearization based control of a rotational hydraulic drive[J].Control Engineering Practice,2007,15(12):1495-1507.

[10] 張從鵬,劉強.直線電機定位平臺的摩擦建模與補償[J].北京航空航天大學學報,2008,34(1):47-50.Zhang C P,Liu Q.Friction modeling and compensation of positioning stage driven by linear motors[J].Journal of Beijing University of Aeronautics and Astronautics,2008,34(1):47-50(in Chinese).

[11] 楊安元,楊雪.液壓系統(tǒng)的減振方法研究[J].液壓與氣動,2004(2):51-53.Yang A Y,Yang X.Research on approaches to weakening vibration of hydraulic system[J].Chinese Hydraulics& Pneumatics,2004(2):51-53(in Chinese).

[12] 陳予恕,唐云,陸啟韶,等.非線性動力學中的現(xiàn)代分析方法[M].北京:科學出版社,2000:1-50.Chen Y S,Tang Y,Lu Q S,et al.Modern analysis methods from nonlinear dynamics[M].Beijing:Science Press,2000:1-50(in Chinese).

[13] 師漢民.機械振動系統(tǒng)[M].武漢:華中科技大學出版社,2004:210-212.Shi H M.Vibration system[M].Wuhan:Huazhong University of Science & Technology Press,2004:210-212(in Chinese).

[14] 劉延柱,陳立群.非線性振動[M].北京:高等教育出版社,2001:60-63.Liu Y Z,Chen L Q.Nonlinear vibration[M].Beijing:Higher Education Press,2001:60-63(in Chinese).

[15] 姜萬錄,劉思遠,張齊生.液壓故障的智能信息診斷與監(jiān)測[M].北京:機械工業(yè)出版社,2013:108-180.Jiang W L,Liu S Y,Zhang Q S.Intelligent information diagnosis and monitoring on hydraulic fault[M].Beijing:China Machine Press,2013:108-180(in Chinese).

猜你喜歡
振動系統(tǒng)
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
Smartflower POP 一體式光伏系統(tǒng)
噴水推進高速艇尾部振動響應分析
WJ-700無人機系統(tǒng)
ZC系列無人機遙感系統(tǒng)
北京測繪(2020年12期)2020-12-29 01:33:58
This “Singing Highway”plays music
基于PowerPC+FPGA顯示系統(tǒng)
半沸制皂系統(tǒng)(下)
振動攪拌 震動創(chuàng)新
中國公路(2017年18期)2018-01-23 03:00:38
中立型Emden-Fowler微分方程的振動性
主站蜘蛛池模板: 久久精品国产亚洲AV忘忧草18| 国产一区免费在线观看| 欧美一级在线| 日韩一二三区视频精品| 日本成人一区| 久996视频精品免费观看| 亚洲高清日韩heyzo| 色屁屁一区二区三区视频国产| 国产美女一级毛片| 国产精品色婷婷在线观看| 成人无码区免费视频网站蜜臀| 国产人成网线在线播放va| 国产成人综合在线观看| 国产农村精品一级毛片视频| 亚洲精品无码不卡在线播放| 综1合AV在线播放| 亚洲国内精品自在自线官| 综1合AV在线播放| 国产黑丝一区| 一级毛片网| 亚洲无码在线午夜电影| 97久久免费视频| 一级片一区| 日韩在线视频网站| 国产激爽大片在线播放| 欧美成人怡春院在线激情| 国产资源站| 欧美成人午夜影院| 日韩高清中文字幕| 国产高清免费午夜在线视频| 日韩视频精品在线| 日韩 欧美 小说 综合网 另类| 尤物亚洲最大AV无码网站| 国产欧美日韩资源在线观看| 国产成人精品在线| 男女猛烈无遮挡午夜视频| 99青青青精品视频在线| 亚洲国产精品无码AV| 国产精品丝袜视频| 69免费在线视频| 国产精品思思热在线| 国产精品尤物铁牛tv | 在线精品自拍| 精品色综合| 欧美色综合网站| 97亚洲色综久久精品| 欧美一级黄色影院| 亚洲成在人线av品善网好看| 国产精品一区二区不卡的视频| 亚洲高清在线播放| 国产精品欧美日本韩免费一区二区三区不卡 | 在线无码私拍| 中文字幕无线码一区| 亚洲无码高清一区二区| 91国内视频在线观看| 国产精品极品美女自在线| 亚洲精品久综合蜜| 四虎国产在线观看| 免费国产一级 片内射老| 亚洲天堂首页| 色综合中文| 亚洲欧美成人影院| 国产免费久久精品44| 久久婷婷色综合老司机| 色悠久久久久久久综合网伊人| 亚洲美女AV免费一区| 国产日韩欧美在线视频免费观看| 日本三级精品| 女人毛片a级大学毛片免费| 久久人妻xunleige无码| 国模私拍一区二区| 久久久久亚洲AV成人网站软件| 欧美亚洲第一页| 午夜不卡福利| 麻豆AV网站免费进入| 亚洲91精品视频| 99在线免费播放| 午夜a视频| 国产91丝袜| 日韩成人在线视频| 久久夜夜视频| 国产白浆一区二区三区视频在线|