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

功率分流齒輪傳動(dòng)參數(shù)空間解域界及吸引域全局特性

2021-06-10 00:28:44林何洪靈劉霞胥光申
振動(dòng)工程學(xué)報(bào) 2021年2期

林何 洪靈 劉霞 胥光申

摘要: 為揭示齒輪系統(tǒng)動(dòng)力學(xué)全局解域特性,建立了功率分流直齒輪系統(tǒng)非線性振動(dòng)模型,結(jié)合胞映射與區(qū)域分解思想構(gòu)建了主要激振參數(shù)下的系統(tǒng)全局解域求解算法。求解了阻尼比分別與綜合傳動(dòng)誤差、嚙合頻率和齒側(cè)間隙等參數(shù)配置下的平面解域界結(jié)構(gòu),揭示了兩參量解域區(qū)內(nèi)的動(dòng)態(tài)行為潛在趨勢(shì),如倍周期分岔序列、混沌窗口帶等的遷移演化。借助Lyapunov指數(shù)對(duì)不同誤差幅值下分岔通道結(jié)構(gòu)進(jìn)行追蹤,驗(yàn)證了分岔節(jié)點(diǎn)與解域界子域臨界點(diǎn)相吻合。數(shù)值求解多組阻尼比下的吸引域全局性態(tài),發(fā)現(xiàn)混沌吸引域與周期1吸引子吸引域間相互擴(kuò)張與進(jìn)退演變激烈,二者胞域邊界處系統(tǒng)局部平衡態(tài)勢(shì)異常脆弱;當(dāng)阻尼比為0.04時(shí)吸引域行為相對(duì)敏感,多吸引子共存現(xiàn)象突出。其結(jié)果可為齒輪系統(tǒng)振動(dòng)優(yōu)化及參數(shù)全局設(shè)計(jì)提供參考。

關(guān)鍵詞: 機(jī)械振動(dòng); 功率分流傳動(dòng); 胞映射; 解域界; 吸引域

中圖分類號(hào): TH113.1; TH132.41; O322 ? ? 文獻(xiàn)標(biāo)志碼: A ? ?文章編號(hào): 1004-4523(2021)02-0235-08

DOI:10.16385/j.cnki.issn.1004-4523.2021.02.003

引 ?言

在裝備制造領(lǐng)域齒輪傳動(dòng)扮演了不可或缺的重要角色,特別在航空航天、新能源汽車、風(fēng)電機(jī)組及高速軌道交通裝備等國(guó)民經(jīng)濟(jì)重點(diǎn)行業(yè),其需求更為強(qiáng)烈,其中,功率分流齒輪傳動(dòng)形式因其獨(dú)特的結(jié)構(gòu)布局在功重比及緊湊性等方面展現(xiàn)出突出優(yōu)勢(shì),正越來(lái)越受到重視[1?2]。振動(dòng)和噪聲是衡量齒輪傳動(dòng)動(dòng)態(tài)品質(zhì)的顯著要素,通常當(dāng)系統(tǒng)運(yùn)行至“不良”參數(shù)激振區(qū),如間隙或時(shí)變剛度等的部分非線性激勵(lì)響應(yīng)帶,系統(tǒng)的瞬態(tài)碰振或混沌傾向?qū)@著增大[3]。實(shí)際中齒輪副的工作難以避免地遭受內(nèi)、外部激勵(lì),反映在振動(dòng)模型中即多參量間的作用最終誘發(fā)系統(tǒng)復(fù)雜動(dòng)態(tài)反饋,如分岔、跳躍和混沌等。因此,全面地開展齒輪副傳動(dòng)工況因素的監(jiān)控對(duì)合理預(yù)測(cè)或規(guī)避無(wú)規(guī)則振動(dòng)尤為必要。劉鎮(zhèn)星等[4]考慮了艦船搖擺環(huán)境下滑動(dòng)軸承?齒輪副系統(tǒng)的動(dòng)力學(xué)特性,觀測(cè)到較小齒隙下?lián)u擺引起的擠齒現(xiàn)象。魏靜等[5]建立了針對(duì)齒廓修形與齒向修形時(shí)斜齒輪傳動(dòng)剛度與誤差非線性耦合模型,分析了修形參數(shù)對(duì)系統(tǒng)動(dòng)態(tài)特性的影響。Yang等[6]從輪齒裂紋擴(kuò)展的角度,求解了裂紋尺寸對(duì)3自由度齒輪系統(tǒng)非線性振動(dòng)的影響。現(xiàn)有齒輪系統(tǒng)模型動(dòng)力學(xué)研究大都建立在單一參數(shù)激勵(lì)通道下,因分析手段的差異,對(duì)實(shí)際多源激勵(lì)的協(xié)同影響未能深刻揭示,同時(shí)響應(yīng)結(jié)果的完整性呈現(xiàn)也受到局限。

胞映射思想最初由Hsu提出,用于離散空間域胞化求解,目前已擴(kuò)展出插值胞映射、點(diǎn)映射?胞映射及廣義胞映射圖論法等諸多改進(jìn)版本,兼顧可行性的同時(shí)效率與精度也不斷提升。在非線性動(dòng)力學(xué)方面,胞化離散不僅可對(duì)系統(tǒng)狀態(tài)空間進(jìn)行求解,還可從不同維度方向自由配置受關(guān)注的參數(shù)區(qū),得到解域界細(xì)節(jié)[7?9]。解域界分析以區(qū)域離散技術(shù)為基礎(chǔ),結(jié)合胞映射揭示參數(shù)空間和狀態(tài)空間特定目標(biāo)域下的潛藏振動(dòng)信息,對(duì)指導(dǎo)齒輪系統(tǒng)減振控制與動(dòng)載優(yōu)化具有積極意義。Kahraman和Singh借助Poincaré映射和相空間離散追蹤了多組吸引子吸引域的共存行為[10]。Farshidianfar和Saghafi基于Melnikov方法數(shù)值計(jì)算了單對(duì)齒輪副在全局域中的同宿分岔和混沌遷移路徑[11?12],并對(duì)混沌參數(shù)帶進(jìn)行了綜合性預(yù)測(cè)。林何等[13]采用解域離散與Lyapunov指數(shù)等工具獲取了并車螺旋錐齒輪系統(tǒng)兩參量平面上的解域界和倍周期分岔分布形態(tài)。在齒輪動(dòng)力學(xué)領(lǐng)域[14?15],各類嚙合激勵(lì)參數(shù)是與系統(tǒng)動(dòng)態(tài)響應(yīng)和振動(dòng)行為息息相關(guān)的關(guān)聯(lián)要素,借助解域界和胞映射可將系統(tǒng)多源激振參數(shù)間共同作用下的動(dòng)力學(xué)行為整體性態(tài)從更廣層面高效、詳盡地完成解析。本文對(duì)功率分流齒輪系統(tǒng)開展全局特性研究以不同維度參數(shù)域和特定初值狀態(tài)域?yàn)檠芯繉?duì)象,解域界直觀呈現(xiàn)系統(tǒng)受激響應(yīng)振動(dòng)變化趨勢(shì),在二維視角下揭示全局分岔與混沌路徑,相空間吸引域信息反饋出嚙合振動(dòng)對(duì)初始狀況的敏感依賴性,揭示吸引子共存及吸引域演化等全局性態(tài)。

2 全局解域特性

假定從動(dòng)輪g_1,g_2幾何結(jié)構(gòu)與物理屬性完全相同;主、從動(dòng)輪齒數(shù)分別為32,105,轉(zhuǎn)動(dòng)慣量分別為0.003 kg·m2,0.197 kg·m2,模數(shù)為3,平均嚙合剛度為k=7.8×〖10〗^8 N/m;輸入扭矩T_0=188 N·m;阻尼比取ξ∈[0.01, 0.18];量綱一齒側(cè)間隙取b∈[0.1, 0.4];量綱一嚙合頻率取ω∈[0.1, 0.5]。

2.1 參數(shù)空間解域界

將ξ和E構(gòu)成的參數(shù)域劃分為一系列規(guī)則矩形胞,各維度方向分割300個(gè)胞,以參數(shù)胞幾何中心代表的數(shù)值實(shí)施計(jì)算。取ξ×E有限域[0.05,0.18]×[0.05,0.21],得到參數(shù)解域界如圖4所示,圖中Pi代表周期態(tài)為i。解域界內(nèi)部是由聚集的1,2,3,4和8周期的參數(shù)胞構(gòu)成,這些參數(shù)胞預(yù)示著系統(tǒng)的穩(wěn)態(tài)行為和態(tài)勢(shì)??梢钥闯鯬1和P2狀態(tài)參數(shù)胞占據(jù)了較大片域,具有良好的穩(wěn)定性,后續(xù)的4周期胞歷經(jīng)短暫參數(shù)帶直接轉(zhuǎn)為8周期響應(yīng)胞,再由8周期通道快速進(jìn)入混沌域,這一過(guò)程與倍周期分岔規(guī)律相一致,同時(shí)也揭示了分岔通道的分布結(jié)構(gòu)。當(dāng)阻尼比增加到0.08時(shí),解域界中出現(xiàn)了P3響應(yīng)區(qū),促使混沌區(qū)提前出現(xiàn)并呈擴(kuò)大趨勢(shì)。從ξ×E域中豎直方向可以看出阻尼比在從0.05增至0.18的過(guò)程中,周期狀態(tài)間的轉(zhuǎn)變更為顯著,表明系統(tǒng)振動(dòng)行為對(duì)阻尼比變化敏感。

從圖4中選取E=0.1,ξ∈[0.05, 0.18]工況參數(shù),由圖5可知嚙合運(yùn)動(dòng)在阻尼比減小時(shí)經(jīng)歷了倍周期分岔,依次出現(xiàn)1周期、2周期、4周期,…,混沌等轉(zhuǎn)變;此后,混沌域內(nèi)局部區(qū)再次生成倍周期分岔窗口,阻尼比的減小使嚙合運(yùn)動(dòng)復(fù)雜化?;煦鐓^(qū)響應(yīng)的無(wú)規(guī)則性還導(dǎo)致了嚙合副正反向振動(dòng),嚙合位移存在負(fù)值,X>-1表明未產(chǎn)生齒背沖擊。圖6中兩條最大Lyapunov指數(shù)曲線定量表征了周期區(qū)(λ_1<0)與混沌區(qū)(λ_1>0)的參數(shù)域范圍以及振動(dòng)的非線性強(qiáng)度,λ_1值越大表明非線性振動(dòng)越強(qiáng),λ_1曲線與分岔節(jié)點(diǎn)相吻合,二者共同驗(yàn)證了圖4中解域界結(jié)果的準(zhǔn)確性。

圖7中解域界由ξ×ω∈[0.05, 0.18]×[0.2, ? 0.6]構(gòu)成,域內(nèi)1周期和2周期為穩(wěn)定后的主要狀態(tài),其中ξ較大或ω較小時(shí)系統(tǒng)易趨向簡(jiǎn)單周期運(yùn)動(dòng)。1周期胞域內(nèi)出現(xiàn)不穩(wěn)定混沌,在阻尼比方向呈狹長(zhǎng)帶;隨著嚙合頻率的增大,系統(tǒng)運(yùn)動(dòng)狀態(tài)切換更頻繁并最終轉(zhuǎn)入混沌域,解域界直觀細(xì)致地展示了系統(tǒng)狀態(tài)的變化。圖8為ξ和b構(gòu)成的解域界,隨阻尼比增大嚙合振動(dòng)變化層次性較明顯,當(dāng)b>0.2時(shí),ξ在[0.07, 0.1]區(qū)間以混沌振動(dòng)為主,該參數(shù)域界上邊界胞較為光滑,阻尼比方向呈周期倍化演變且狀態(tài)轉(zhuǎn)換未出現(xiàn)周期性跨越,同時(shí)齒側(cè)間隙相比阻尼比對(duì)系統(tǒng)狀態(tài)改變的驅(qū)動(dòng)性要弱。從ξ×b中分別取胞值(0.15, 0.3),(0.12, 0.3),(0.115, 0.3)和(0.08, 0.3),在圖9中分別得到P1,P2,P4和混沌吸引子,與ξ×b解域界結(jié)果相一致。

2.2 狀態(tài)空間吸引域

狀態(tài)空間全局特性揭示系統(tǒng)在不同初值條件下長(zhǎng)期運(yùn)動(dòng)趨勢(shì)。在X×X ˙∈[-3, 3]×[-3, 3]區(qū)域內(nèi)取400×400個(gè)正規(guī)胞,在圖4混沌態(tài)子域中取ξ為0.1,E為0.205時(shí)測(cè)試吸引域性態(tài)(如圖10所示)。全局域主要由1周期吸引域(空白區(qū))和混沌吸引域(灰色區(qū))構(gòu)成,二者基本覆蓋整個(gè)求解區(qū),域邊界光滑連續(xù),域內(nèi)外層次分明,兩側(cè)存在跳變的可能,因外層為1周期域,故域外胞穩(wěn)定性更好。符號(hào)B表示域邊界,即1周期吸引域與混沌吸引域的邊界胞構(gòu)成,總數(shù)量為2479,此時(shí)中心區(qū)域吸引域演變更為激烈,兩吸引域間交錯(cuò)嵌入現(xiàn)象明顯。計(jì)算知混沌胞的占比為40.36%,對(duì)比解域界圖4知子域驟變會(huì)導(dǎo)致系統(tǒng)狀態(tài)急劇轉(zhuǎn)變,致使局部平衡態(tài)勢(shì)異常脆弱,吸引子在受碰激發(fā)時(shí)極易趨向更穩(wěn)定軌道上。

改變?chǔ)沃?.02時(shí)全局吸引域如圖11所示,相比ξ為0.1時(shí)吸引域間更加嵌入,但主要敏感空間有所縮小,覆蓋區(qū)由之前的[-3,3]×[-3,3]收縮至[-1,1]×[-1,1],全局域仍以1周期域和混沌域構(gòu)成,兩者狀態(tài)胞的數(shù)量分別為89760和71000,此時(shí),兩吸引域交錯(cuò)更深且接觸范圍更廣,表明低阻尼時(shí)系統(tǒng)響應(yīng)狀態(tài)是較敏感的,起始狀態(tài)受干擾時(shí)易從一個(gè)吸引域轉(zhuǎn)入另一個(gè)吸引域。

當(dāng)ξ增大為0.03時(shí)(如圖12所示),相空間中僅存在混沌域和P1吸引域,各吸引域均有向外發(fā)散擴(kuò)充態(tài)勢(shì),混沌吸引域螺旋結(jié)構(gòu)相較于ξ為0.02時(shí)更寬,結(jié)構(gòu)趨向簡(jiǎn)單,各吸引域胞集變得更為集中,此時(shí)P1吸引子與混沌吸引子間相互轉(zhuǎn)移敏感性下降,表明系統(tǒng)對(duì)初始條件具有更強(qiáng)的抗干擾能力。

進(jìn)一步增大ξ為0.04,由圖13和14可見混沌吸引域分布于4周期吸引域內(nèi)側(cè),此時(shí)域邊界分形特征清晰,這些鑲嵌于4周期吸引域內(nèi)部的初值將在穩(wěn)態(tài)時(shí)促使系統(tǒng)收斂入混沌態(tài)。狀態(tài)變化主要集中在4周期吸引域的外圍,1周期吸引域在整個(gè)域中所占比達(dá)40.77%,混沌胞全局占比約10.97%,從動(dòng)載特性角度知1周期具有更好的穩(wěn)定性,因?yàn)槠浒獦?gòu)成區(qū)域相對(duì)集中且未出現(xiàn)其他散亂的高周期干擾胞。4周期吸引域范圍內(nèi)還出現(xiàn)了少量不穩(wěn)定的2周期吸引胞,分布較為散亂,大部分散落在4周期域內(nèi)部。由于存在4周期、混沌和2周期吸引子共存聚集,此狀態(tài)域內(nèi)系統(tǒng)更易出現(xiàn)波動(dòng),少量2周期胞在受激時(shí)易被吸引到臨近的域內(nèi),促使原吸引子失穩(wěn)再遷移。繼續(xù)計(jì)算ξ為0.05時(shí),發(fā)現(xiàn)阻尼比的增大迫使4周期吸引域加速退化,1周期吸引域和混沌域不同程度地?cái)U(kuò)張。可見阻尼比越大系統(tǒng)振動(dòng)越易收斂到穩(wěn)定的1周期上,即增大阻尼比可在一定程度上抑制齒輪系統(tǒng)發(fā)生多周期或混沌振動(dòng)。

3 結(jié) ?論

(1)基于區(qū)域分解技術(shù)與胞映射思想構(gòu)建了齒輪系統(tǒng)參數(shù)空間解域界求解數(shù)值算法,并給出了其求解流程。

(2)解域界對(duì)混沌窗口帶、分岔節(jié)點(diǎn)及子域邊界等結(jié)構(gòu)信息給出了準(zhǔn)確預(yù)測(cè),綜合傳動(dòng)誤差與阻尼比構(gòu)成的解域界中模擬出了倍周期分岔序列演化發(fā)展歷程。齒側(cè)間隙與阻尼比構(gòu)成的解域界中齒側(cè)間隙相較于阻尼比對(duì)系統(tǒng)狀態(tài)改變的驅(qū)動(dòng)性要弱。

(3)解域邊界處系統(tǒng)狀態(tài)局部平衡較為脆弱,阻尼比增大,混沌吸引子吸引域不斷退化,1周期吸引子吸引域進(jìn)而擴(kuò)張,全局吸引域邊界分形特征趨向簡(jiǎn)單,增大阻尼比可抑制齒輪系統(tǒng)對(duì)初值條件的敏感依賴性。

參考文獻(xiàn):

[1] Roozegar M, Angeles J. A two-phase control algorithm for gear-shifting in a novel multi-speed transmission for electric vehicles[J]. Mechanical Systems & Signal Processing, 2018, 104:145-154.

[2] 向 ?玲,高 ?楠,唐 ?亮,等.支承剛度變化下風(fēng)電齒輪傳動(dòng)系統(tǒng)的非線性動(dòng)力學(xué)特性[J].振動(dòng)與沖擊, 2019, 38(01): 111-117.

Xiang Ling, Gao Nan, Tang Liang, et al. Nonlinear dynamic characteristics of wind turbine gear transmission system with varying support stiffness [J]. Journal of Vibration and Shock, 2019, 38(01): 111-117.

[3] 桂永方,朱如鵬,靳廣虎,等.間隙非線性圓柱齒輪分流傳動(dòng)系統(tǒng)動(dòng)力學(xué)與均載特性分析[J].振動(dòng)與沖擊, 2014, 33(18): 177-184.

Gui Yongfang, Zhu Rupeng, Jin Guanghu, et al. Dynamic and load sharing characteristic analysis of a nonlinear cylindrical gear split-torque transmission system with backlash[J]. Journal of Vibration and Shock, 2014, 33(18): 177-184.

[4] 劉鎮(zhèn)星,劉占生,于香宇,等.艦船搖擺作用下滑動(dòng)軸承-齒輪副系統(tǒng)動(dòng)力學(xué)分析[J].機(jī)械工程學(xué)報(bào), 2018, 54(17): 226-234.

Liu Zhenxing, Liu Zhansheng, Yu Xiangyu, et al. Dynamic analysis of journal bearing-gear system under swing movement of the ship[J]. Journal of Mechanical Engineering, 2018, 54(17): 226-234.

[5] 魏 ?靜,王剛強(qiáng),秦大同,等.考慮修形的斜齒輪系統(tǒng)非線性激勵(lì)與動(dòng)力學(xué)特性研究[J].振動(dòng)工程學(xué)報(bào), 2018, 31(4):561-572.

Wei Jing, Wang Gangqiang, Qin Datong, et al. Nonlinear excitation and dynamic characteristics of helical gear system with considering modification[J]. Journal of Vibration Engineering, 2018, 31(4):561-572.

[6] Yang Yi, Xia Wenkai, Han Jianming, et al. Vibration analysis for tooth crack detection in a spur gear system with clearance nonlinearity[J]. International Journal of Mechanical Sciences, 2019, 157: 648-661.

[7] Sun Jianqiao, Xiong Furui. Cell mapping methods beyond global analysis of nonlinear dynamic systems[J]. Advances in Mechanics, 2017, 47: 150-177.

[8] Li Zhigang, Jiang Jun, Li Jing, et al. A subdomain synthesis method for global analysis of nonlinear dynamical systems based on cell mapping[J]. Nonlinear Dynamics, 2018,95(1):715-726.

[9] 徐 ?偉,岳曉樂(lè),韓 ?群.胞映射方法及其在非線性隨機(jī)動(dòng)力學(xué)中的應(yīng)用[J]. 動(dòng)力學(xué)與控制學(xué)報(bào), 2017, 15(3): 12-20.

Xu Wei, Yue Xiaole, Han Qun. Cell mapping method and its applications in nonlinear stochastic dynamical systems[J]. Journal of Dynamics and Control, 2017, 15(3):12-20.

[10] Kahraman A, Singh R. Non-linear dynamics of a spur gear pair[J]. Journal of Sound and Vibration, 1990, 142(1): 49-75.

[11] Farshidianfar A, Saghafi A. Global bifurcation and chaos analysis in nonlinear vibration of spur gear systems[J]. Nonlinear Dynamics, 2014, 75(4): 783-806.

[12] Saghafi A, Farshidianfar A. An analytical study of controlling chaotic dynamics in a spur gear system[J]. Mechanism & Machine Theory, 2016, 96:179-191.

[13] 林 ?何,王三民,董金城.并車螺旋錐齒輪傳動(dòng)動(dòng)力學(xué)參數(shù)二維域界結(jié)構(gòu)分析[J].航空動(dòng)力學(xué)報(bào), 2017, 32(8): 2017-2024.

Lin He, Wang Sanmin, Dong Jincheng. Two-dimensional domain structure of dynamical parameters of combining spiral gear transmission[J]. Journal of Aerospace Power, 2017, 32(8):2017-2024.

[14] 唐進(jìn)元,熊興波,陳思雨.基于圖胞映射方法的單自由度非線性齒輪系統(tǒng)全局特性分析[J].機(jī)械工程學(xué)報(bào), 2011, 47(5): 59-65.

Tang Jinyuan, Xiong Xingbo, Chen Siyu. Analysis of global character of single degree of freedom nonlinear gear system based on digraph cell mapping method[J]. Journal of Mechanical Engineering, 2011, 47(5): 59-65.

[15] Gou Xiangfeng, Zhu Lingyun, Chen Dalin. Bifurcation and chaos analysis of spur gear pair in two-parameter plane[J]. Nonlinear Dynamics, 2015, 79(3):2225-2235.

Global characteristics of basin of attraction and parameterized solution domain of power-split gear transmission

LIN He1,2,3, HONG Ling2, LIU Xia1,3, XU Guang-shen1,3

(1. School of Mechanical and Electrical Engineering, Xi'an Polytechnic University, Xi'an 710048, China;

2. State Key Laboratory for Strength and Vibration of Mechanical Structures, Xi'an Jiaotong University, Xi'an 710049, China;

3. Xi'an Key Laboratory of Modern Intelligent Textile Equipment, Xi'an Polytechnic University, Xi'an 710600, China)

Abstract: In order to investigate the global characteristics of dynamic solutions of gear systems, a nonlinear dynamical model of the power-split spur gear transmission is established, and the calculation algorithm regarding global solutions under main excitation parameters is deduced based on cell mapping method (CMM) as well as domain decomposition method (DDM). Parametric planar solution domains constructed by damping ratio respectively with synthetically transmission error, mesh frequency and backlash are computed, and the potential global evolution behaviors within solution domains are exhibited, such as period-doubling bifurcation cascades, chaotic window zones. The bifurcation routes inside solution domain with respect to varied error magnitudes are tracked by applying largest Lyapunov exponent, which demonstrate that bifurcation nodes are in consistent with the subdomain boundary points presented in parameterized solution domain. By numerically calculating the global behaviors of the basin of attraction under multiple damping ratios, it is shown that mutual expansion and retrogression between the chaotic basin of attraction and the period 1 basin of attraction are remarkably, and the local equilibrium of the cells nearby the boundary is extremely unstable. The basin of attraction is sensitively while the damping ratio reaches 0.04, and multiple attractors coexisting phenomenon exhibits significantly. The result could provide references for vibration optimization or even global design of dynamic parameters for gear system.

Key words: mechanical vibration; power-split gear transmission; cell mapping method; parameterized solution domain; basin of attraction

作者簡(jiǎn)介: 林 ?何(1985-),男,副教授。電話:(029)62779107;E-mail: linhe@xpu.edu.cn

通訊作者: 胥光申(1964-),男,教授。電話:(029)62779109;E-mail: xugs988@126.com

主站蜘蛛池模板: 特级精品毛片免费观看| 58av国产精品| 亚洲狼网站狼狼鲁亚洲下载| 无码福利日韩神码福利片| 伊人91视频| 亚洲人妖在线| 日本不卡免费高清视频| 国产精品爽爽va在线无码观看 | 免费亚洲成人| 国产主播一区二区三区| 四虎永久免费地址在线网站| 欧美午夜在线播放| 欧美、日韩、国产综合一区| 一级毛片高清| 亚洲国产在一区二区三区| 高h视频在线| 在线另类稀缺国产呦| 国产成熟女人性满足视频| 国产黑丝视频在线观看| 国产成人欧美| 亚洲日韩精品伊甸| 久久夜色精品国产嚕嚕亚洲av| 欧美精品综合视频一区二区| 婷婷五月在线| 国产激情无码一区二区APP| av免费在线观看美女叉开腿| 白丝美女办公室高潮喷水视频| 久久久久免费精品国产| 91在线激情在线观看| 91成人在线免费观看| 亚洲国产成人久久77| 91伊人国产| 欧美特黄一级大黄录像| 亚洲an第二区国产精品| 九九九精品成人免费视频7| 国产成人精品优优av| 亚洲精品人成网线在线| 国产美女在线观看| 国产玖玖视频| 中文字幕亚洲专区第19页| 午夜视频www| 亚洲,国产,日韩,综合一区 | 夜夜操狠狠操| 国产原创演绎剧情有字幕的| 亚洲国产成人在线| 亚洲手机在线| 91www在线观看| 日韩欧美网址| 69av免费视频| 成人另类稀缺在线观看| 久久人搡人人玩人妻精品 | 国产精品浪潮Av| 亚洲视频四区| 亚洲AV电影不卡在线观看| 57pao国产成视频免费播放| 九九九精品成人免费视频7| 国产在线视频自拍| 成人伊人色一区二区三区| 9丨情侣偷在线精品国产| 成人在线天堂| 黄色网址免费在线| 亚洲国产欧洲精品路线久久| 国产美女一级毛片| 一本一本大道香蕉久在线播放| 一区二区在线视频免费观看| 欧美精品aⅴ在线视频| 国产美女在线观看| 午夜福利在线观看入口| 欧美日韩福利| 伊人狠狠丁香婷婷综合色| 色婷婷电影网| 国产日韩精品欧美一区灰| 久草视频福利在线观看| 伊人成色综合网| 成人国产精品2021| 伊人久久婷婷| 国产农村妇女精品一二区| 国产一区亚洲一区| 40岁成熟女人牲交片免费| 亚洲国产精品日韩av专区| 欧美午夜视频| 在线观看精品自拍视频|