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

斜齒行星傳動多體動力學建模與分析

2014-09-05 09:57:58劉先增
振動與沖擊 2014年7期
關鍵詞:模型

張 俊, 劉先增

(安徽工業(yè)大學 機械工程學院,安徽 馬鞍山 243032)

行星齒輪傳動的振動和噪聲是影響系統(tǒng)可靠性、壽命及操作環(huán)境的關鍵因素。圍繞該類傳動的動力學問題,學術界開展了大量研究,內(nèi)容涉及動力學建模、固有特性分析、動態(tài)響應求解、振動噪聲抑制等多個方面[1]。這其中,動力學建模和分析是進行后續(xù)動力性能研究及減振降噪的理論基礎。相比直齒行星傳動,斜齒行星傳動的結構、受力更為復雜,其動力學建模和分析也更具挑戰(zhàn)性,故迄今為止針對斜齒行星傳動動力特性的研究較少[2-4]。

按建模方法和考慮因素的不同,大致可將行星傳動的動力學模型分為解析模型、有限元模型和多體模型三類。其中,解析模型因建模簡單、求解容易等優(yōu)點而被廣泛采用。Kahraman[2]建立了斜齒行星傳動的三維解析模型,并依托所建模型分析了行星輪嚙合相位等參數(shù)對輪系動力特性的影響。該解析模型中,各構件均被視為具有6自由度的剛體,各剛體間的嚙合和支承簡化為具有集中效應的彈簧阻尼單元。延續(xù)這一思路,Lin和Parker等通過建立計入構件平移、扭轉(zhuǎn)運動的解析模型獲得了直齒行星傳動固有特性的解析表達式,并將系統(tǒng)的自由振動歸納為扭轉(zhuǎn)(rotational)、平移(translational)和行星輪(planet)三種模式[5]。采用類似解析模型,學者們進一步研究了參數(shù)靈敏度[6]、模態(tài)躍遷[7]、參數(shù)穩(wěn)定性[8]、非均布行星輪系統(tǒng)[9]、復合行星輪系[10]及人字齒、錐齒行星輪系模態(tài)[11]等問題。

不同于解析模型將輪系視為具有集中效應的彈簧質(zhì)量系,有限元模型可有效計入系統(tǒng)各環(huán)節(jié)的影響。Kahraman等運用有限元法建立了行星輪系的準靜態(tài)受力模型,分析了內(nèi)齒圈柔性對齒輪應力和行星輪載荷分配的影響[12]。采用類似手段,該文作者又進一步研究了內(nèi)齒圈柔性對系統(tǒng)動態(tài)特性的影響[13]。文獻[12-13]的研究表明,構件柔性有助于補償因齒輪和系桿制造安裝誤差引起的非均載效應,且對系統(tǒng)動態(tài)特性具有重要影響。此外,文獻[14-16]也采用類似的有限元模型對行星輪系的動態(tài)特性進行了研究,也得出了一些重要結論。

盡管有限元模型相較于解析模型具備更高的分析精度,但由于該類模型的建立和求解均較為費時,故不適用于需要反復迭代的初始設計階段。相比之下,介于解析模型和有限元模型之間的多體動力學模型既能很好地反映行星輪系的質(zhì)量集中特點,又能有效規(guī)避有限元模型的冗余繁雜,為行星輪系的動力學分析提供了新的解決方案。特別是隨著近年來一些商用軟件如Simpack、Virtual.lab、ADAMS、RecurDyn等的成功應用,采用多體動力模型對行星輪系進行動力學分析已成為行之有效的手段[17-19]。

有鑒于此,本文擬運用多體動力學方法,建立斜齒行星傳動的多體動力學模型,進而依托所建模型對輪系進行自由振動分析,并將分析結果與前人提出的解析模型的結果作對比,以驗證多體動力學模型的可靠性。在此基礎上,進一步開展斜齒行星輪系的動態(tài)響應和參數(shù)影響分析,以期為后續(xù)的動態(tài)設計與性能優(yōu)化提供理論依據(jù)。

1 多體動力學建模

為方便建模并不失一般性,作如下假設:

(1) 忽略各齒輪體、系桿和箱體的柔性,將其按剛體處理;

(2) 計入各齒輪副的嚙合變形和各構件的支承變形,以具有等效剛度的線性時變/不變彈簧代替;

(3) 各齒輪副的嚙合力始終作用在嚙合平面內(nèi),并與理論接觸線垂直;

(4) 各嚙合輪齒間不存在脫齒和嚙入/嚙出沖擊;

(5) 不考慮原動機和負載的慣性;

(6) 各行星輪的質(zhì)量、慣量和平均嚙合剛度均相同。

基于上述假設,建立如圖1所示的多體動力學模型。

圖1 斜齒行星傳動多體動力學模型

圖中,除箱體按固定件處理,其余構件均擁有6自由度。太陽輪和行星輪、行星輪和內(nèi)齒圈之間以線性時變彈簧聯(lián)接,彈簧剛度取為相應齒輪副的嚙合剛度(詳見GB/T 3480—1997斜齒輪平均嚙合剛度的計算)。太陽輪、內(nèi)齒圈及系桿與箱體之間、行星輪與系桿之間的聯(lián)接均以6自由度襯套代替,襯套各方向數(shù)值的計算可參見文獻[20]。篇幅所限,上述各剛度數(shù)值的計算不再詳列。

圖2為斜齒行星傳動的機構示意圖,為清晰計未示出系桿。定義全局坐標系O-XYZ,其中原點O取為輪系的幾何中心,X軸取為水平向右,Y軸取為垂直向上,Z軸由右手定則確定。在該坐標系下,各構件的彈性位移以qi表示,且有qi=[xi,yi,zi,uxi,uyi,uzi]T(i=s,r,c,1,2,…,N),下標s、r、c、n(n=1,2,…,N)分別表示太陽輪、內(nèi)齒圈、系桿和第n個行星輪。各符號的含義可參見文獻[3]。

圖2 斜齒行星傳動嚙合關系示意圖

根據(jù)圖2所示的幾何關系,不難推導出內(nèi)、外嚙合副沿理論嚙合線的相對位移δsn、δrn。

(1)

(2)

(3)

(4)

(5)

(6)

式中,ψsn=ψn-αs,ψrn=ψn+αr,ψn、αs、αr分別為第n個行星輪的位置角、外嚙合副嚙合角和內(nèi)嚙合副嚙合角。

由多體動力學理論,系統(tǒng)的動能可寫成如下形式

(7)

式中,mi為各構件在全局坐標系下的質(zhì)量。該式可進一步寫成如下矩陣形式

(8)

同樣不難給出系統(tǒng)的勢能表達式

(9)

上式同樣可寫成矩陣形式如下

(10)

進一步運用拉格朗日方程,可推導出系統(tǒng)的運動微分方程如下

(11)

式中,φ為系統(tǒng)的總體約束矩陣,F(xiàn)為系統(tǒng)外載荷列陣,φq為各構件彈性變形的Jacobian矩陣,μ、Ω為介于1~10間的罰系數(shù)[21]。限于篇幅,略去上述各矩陣元素及標量的具體推導過程。

2 動態(tài)特性仿真

不失一般性,以表1所示的參數(shù)為例進行斜齒行星輪系的動態(tài)特性仿真。設系統(tǒng)的輸入構件為太陽輪,輸出構件為系桿,內(nèi)齒圈與機殼制為一體。

表1 斜齒行星傳動算例系統(tǒng)計算參數(shù)

2.1 自由振動分析

求解系統(tǒng)運動微分方程的特征值問題,即可獲知該類傳動的自由振動特性。取系統(tǒng)動力學參數(shù)如表1,行星輪數(shù)目分別取3、4、5時,可得系統(tǒng)各階固有頻率如表2所示。

根據(jù)系統(tǒng)特征值的重根數(shù)、中心構件的振型坐標以及各行星輪間振型坐標的比例等特點,可將系統(tǒng)振型劃分為3類,即:①軸向平移—扭轉(zhuǎn)振動模式;②徑向平移—扭擺振動模式;③行星輪振動模式。各類振型的特點簡述如下:

(1) 軸向平移—扭轉(zhuǎn)振動模式: ①各中心構件即太陽輪、系桿和內(nèi)齒圈只有軸向平移振動和繞z軸的扭轉(zhuǎn)振動,其他方向上的振動均為零;所有行星輪的振型相同。 ②有11個特征值與此模式對應,且相應的特征值均為單根; ③各階固有頻率受行星輪個數(shù)影響,其中除零頻外,1、2階固有頻率隨行星輪個數(shù)增加而單調(diào)遞減,其他階次固有頻率均隨行星輪個數(shù)增加而單調(diào)遞增。

(2) 徑向平移—扭擺振動模式: ①各中心構件即太陽輪、系桿和內(nèi)齒圈只有徑向的平移振動和繞x、y軸徑向扭擺振動,其他方向上的振動均為零; ②有12個特征值與此模式相對應,且相應的特征值均為2重根; ③ 各階固有頻率同樣受行星輪個數(shù)影響,其中除第1階固有頻率隨行星輪個數(shù)的增加呈先增后減趨勢外,第4階固有頻率隨行星輪個數(shù)增加而單調(diào)遞減,其他階次固有頻率均隨行星輪個數(shù)增加而單調(diào)遞增。

(3) 行星輪振動模式: ①各中心構件即太陽輪、系桿和內(nèi)齒圈的扭轉(zhuǎn)、平移運動均為零,且各行星輪的位移為第1個行星輪位移乘以一個系數(shù);②有6個特征值與此模式相對應,其中有5個只在行星輪個數(shù)N>3時出現(xiàn),且重復率為N-3,另一個特征值在行星輪個數(shù)N=3時也存在,重復率為N-2。以上各特征值都不受行星輪個數(shù)影響。

表2 斜齒行星輪系固有頻率

為直觀計,圖3進一步給出了上述3種振動模式的振型示意圖。

圖3 斜齒行星傳動振動模式

將表2的仿真結果與采用文獻[2-3]的集中參數(shù)模型對比,可以發(fā)現(xiàn)在同等條件下,兩類模型所得的各階固有頻率和振型完全一致,表明本文所建模型的正確性,進而可依托該模型預估斜齒行星輪系的動態(tài)特性。

2.2 穩(wěn)態(tài)動力學響應

不妨以表1所示的3行星輪系統(tǒng)為例,采用上述多體動力學模型分析系統(tǒng)的穩(wěn)態(tài)動力學響應。設系統(tǒng)的輸入轉(zhuǎn)速恒為1 500 r/min,負載轉(zhuǎn)矩恒為400 N·m。

首先考察無誤差情況下,系統(tǒng)的動態(tài)嚙合力。圖4和5分別示出了行星輪系各外、內(nèi)嚙合力及其對應的頻譜。

圖4 無誤差情況下輪系外嚙合動態(tài)載荷

圖5 無誤差情況下輪系內(nèi)嚙合動態(tài)載荷

由圖4、圖5可知,太陽輪、內(nèi)齒圈與3個行星輪的動態(tài)嚙合力變化規(guī)律大致相同,只是在數(shù)值和相位上略有差別。一個嚙合周期內(nèi),內(nèi)、外嚙合力的變化顯著,說明系統(tǒng)在嚙合過程中存在一定的沖擊振動。進一步分析可知,上述動態(tài)嚙合力均圍繞831這一均值上下震蕩,而這一均值恰為均載條件下輪系靜態(tài)理論嚙合力。

而從嚙合力對應的頻譜圖可知,在1 152 Hz處有一條幅值突出的譜線,它恰好對應于系統(tǒng)的嚙頻。除此之外,其他幾條幅值較大的譜線正好對應嚙頻的倍頻。此點說明嚙頻激勵是引起斜齒行星輪系的主要激勵源,這也與前人的研究結論相吻合[3-4]。

3 參數(shù)影響分析

斜齒行星輪系的結構復雜,設計參數(shù)多,使得該類傳動的動態(tài)設計和性能優(yōu)化較為困難。為此,采用靈敏度分析方法,考察系統(tǒng)設計參數(shù)對斜齒行星輪系動態(tài)特性的影響,希冀為后續(xù)系統(tǒng)的性能改善提供一定的理論參考。篇幅所限,下文僅給出構件支承剛度和行星輪周向誤差對輪系動態(tài)特性的影響規(guī)律。

3.1 構件支承剛度

為方便分析,不妨定義一個嚙合周期內(nèi)各嚙合力的最大值與靜態(tài)均值的比例為載荷波動系數(shù),并以符號χi表示。顯然,χi越大,表明輪系嚙合過程中的振動越強烈。

采用前述的多體動力學模型,進一步分析各構件支承剛度對載荷波動系數(shù)的影響。在進行靈敏度分析時,僅改變待分析參數(shù),而保持其他參數(shù)不變。

圖6 太陽輪支承剛度對載荷波動系數(shù)的影響

分析結果表明,各中心構件的支承剛度對系統(tǒng)的均載特性影響顯著。隨著太陽輪、行星架、內(nèi)齒圈支承剛度的增加,行星系統(tǒng)的載荷波動系數(shù)變大,表明系統(tǒng)的均載特性變差,振動加劇。這其中,太陽輪支承剛度對系統(tǒng)動態(tài)性能的影響最大,行星架次之,內(nèi)齒圈影響最小。由此可見,降低中心構件的支承剛度有利于實現(xiàn)系統(tǒng)均載進而抑制振動。正是基于這一判斷,在斜齒行星輪系的工程應用中,可使太陽輪、行星架、內(nèi)齒圈其中之一浮動,或同時浮動其中的二者,以降低傳動系統(tǒng)嚙合過程中的振動。限于篇幅,下文僅給出太陽輪支承剛度對載荷波動系數(shù)的影響規(guī)律。

3.2 行星輪周向安裝誤差

再來分析行星輪安裝誤差對系統(tǒng)動態(tài)特性的影響。行星輪安裝誤差可分為徑向誤差和分度誤差。研究表明,與分度誤差相比,行星輪徑向誤差對系統(tǒng)動態(tài)特性的影響很小[22]。故下文僅考慮行星輪分度誤差對系統(tǒng)動態(tài)特性的影響。為便于分析,不妨將行星輪分度誤差換算成與其等價的周向誤差來表征。以符號eci表示周向誤差,并將該誤差納入前述的多體動力學模型。通過求解含誤差的動力學模型,可獲知相關動力參數(shù)。篇幅所限,僅給出其中某一個行星輪(本例為行星輪2)含有35 μm的周向安裝誤差時(即ec2=35 μm)系統(tǒng)的動態(tài)嚙合力及其頻譜圖,其結果如圖7、圖8所示。

圖7 輪系外嚙合動態(tài)載荷(ec2=35 μm)

表3 周向安裝誤差對嚙合力的影響(嚙合力單位N)

圖8 輪系內(nèi)嚙合動態(tài)載荷(ec2=35 μm)

進一步分析可知含不同誤差工況下系統(tǒng)的動態(tài)特性。為方便比較,將各工況下的動態(tài)嚙合力列表處理,其結果如表3所示。

由上述分析結果可知,行星輪周向安裝誤差極大地影響斜齒行星輪系的動態(tài)特性。就本例而言,周向安裝誤差的引入,改變了輪系內(nèi)各嚙合副的動態(tài)嚙合力及各功率流的分配,使得含有誤差的行星輪所承擔的內(nèi)、外嚙合力均值變小,同時使得其他行星輪的內(nèi)、外嚙合力均值變大;由于周向安裝誤差使得各行星輪的嚙合力均值進一步偏離了其靜態(tài)理論值,從而使得系統(tǒng)內(nèi)各嚙合力副的載荷波動系數(shù)增大,降低了系統(tǒng)的均載性能。進一步的數(shù)值分析表明,隨著誤差量的增大,存在周向安裝誤差的行星輪所在的功率支路的內(nèi)、外嚙合力均值單調(diào)遞減,而其他功率支路的內(nèi)、外嚙合力均值單調(diào)遞增。不僅如此,行星輪周向安裝誤差還會影響內(nèi)、外嚙合副在嚙頻倍頻處的嚙合力。周向安裝誤差使各內(nèi)、外嚙合力在部分嚙頻倍頻處的幅值增加,而在部分嚙頻倍頻處的幅值減??;且隨著周向安裝誤差量的增大,各功率支路的內(nèi)、外嚙合力在嚙頻倍頻處的變化規(guī)律不盡相同。由此可見,斜齒行星輪系的動態(tài)特性受行星輪周向安裝誤差的影響顯著,在進行輪系設計、制造和安裝時,必須嚴格控制這一誤差環(huán)節(jié),以降低系統(tǒng)振動提示傳動的動態(tài)性能。

4 結 論

(1) 建立了計入多種影響因素的斜齒行星傳動的多體動力學模型,并據(jù)此分析了傳動系統(tǒng)的自由振動特性,其仿真結果與前人的集中參數(shù)模型所得結果吻合,表明所建多體動力學模型能正確揭示斜齒行星傳動的動態(tài)特性。

(2) 快速求解了系統(tǒng)的穩(wěn)態(tài)動力學響應,獲得了輪系各環(huán)節(jié)的動態(tài)載荷。仿真結果表明,無誤差情況下,輪系各嚙合副的動態(tài)載荷圍繞靜態(tài)理論嚙合力上下波動,且嚙頻激勵是引起系統(tǒng)振動的主要原因。

(3) 中心構件支承剛度對斜齒行星輪系的動態(tài)特性影響明顯。隨著中心構件支承剛度的增加,行星傳動的載荷波動系數(shù)單調(diào)遞增,系統(tǒng)動態(tài)特性變差。這其中,太陽輪支承剛度對系統(tǒng)動態(tài)性能影響最大,行星架次之,內(nèi)齒圈影響最小。

(4) 斜齒行星傳動的動態(tài)特性對行星輪周向安裝誤差較為敏感。行星輪周向安裝誤差不僅改變各功率支路嚙合力的大小和分配情況,還改變了各嚙合副在嚙頻倍頻處的幅值。

參 考 文 獻

[1]YANG Jian-ming, DAI Li-ming. Survey of dynamics of planetary gear trains[J]. Int. J. Materials and Structural Integrity, 2008, 1(4): 302-322.

[2]Kahraman A. Natural modes of planetary gear trains[J]. Journal of Sound and Vibration, 1994, 173(1):125-130.

[3]楊通強, 宋軼民, 張策, 等. 斜齒行星齒輪系統(tǒng)自由振動特性分析[J]. 機械工程學報, 2005, 41(7): 50-55.

YANG Tong-qiang, SONG Yi-min, ZHANG Ce, et al. Propety analysis of free vibration of helical planetary gear trains[J]. Chinese Journal of Mechanical Engineering, 2005, 41(7): 50-55.

[4]Kahraman A. Planetary gear train dynamics[J]. Journal of Mechanical Design, 2002, 38(3): 6-9.

[5]Lin J, Parker R G. Analytical characterization of the unique properties of planetary gear free vibration[J]. Journal of Vibration and Acoustics, 1999, 121(2): 316-321.

[6]Lin J, Parker R G. Sensitivity of planetary gear natural frequencies and vibration modes to model parameters[J]. Journal of Sound and Vibration, 1999, 228(1): 109-128.

[7]王世宇, 宋軼民, 沈兆光, 等. 行星傳動系統(tǒng)的固有特性及模態(tài)躍遷研究[J]. 振動工程學報, 2005, 18(4): 412-417.

WANG Shi-yu, SONG Yi-min, SHEN Zhao-guang, et al. Research on natural characteristics and loci veering of planetary gear transmissions[J]. Journal of Vibration Engineering, 2005, 18(4): 412-417.

[8]Lin J, Parker R G. Planetary gear parametric instability caused by mesh stiffness variation[J]. Journal of Sound and Vibration, 2002, 249(1): 129-145.

[9]Lin J, Parker R G. Structure vibration characteristics of planetary gears with unequally spaced planets[J]. Journal of Sound and Vibration, 2000, 233(5): 921-928.

[10]Kiracofe D, Parker R G. Structured vibration modes of general compound planetary gear systems[J]. ASME Journal of Vibration and Acoustics, 2007, 129(2): 1-16.

[11]BU Zhong-hong, LIU Geng, WU Li-yan. Modal analyses of herringbone planetary gear train with journal bearings[J]. Mechanism and Machine Theory, 2012, 54: 99-115.

[12]Kahraman A, Vijayakar S. Effect of internal gear flexibility on the quasi-static behavior of a planetary gear set[J]. ASME Journal of Mechanical Design, 2001, 123: 408-415.

[13]Kahraman A, Kharazi A A, Umrani M. A deformable body dynamic analysis of planetary gears with thin rims[J]. Journal of Sound and Vibration, 2003, 262: 752-768.

[14]Parker R G. Dynamic response of a planetary gear system using a finite element/contact mechanics model[J]. Journal of Mechanical Design, Transaction of the ASME., 2000, 122: 304-310.

[15]Parker R G. A physical explanation for the effectiveness of planet phasing to suppress planetary gear vibration[J]. Journal of Sound and Vibration, 2000, 236(4): 561-573.

[16]Abousleiman V, Velex P, Becquerelle S. Modeling of spur and helical gear planetary drives with flexible ring gears and planet carries[J]. ASME Journal of Mechanical Design, 2007, 129: 95-106.

[17]Qin D T, Wang J H, Lim T C. Flexible multibody dynamic modeling of a horizontal wind turbine drivetrain system[J]. ASME Journal of Mechanical Design,2009,131: 14501-14508.

[18]Jan H, Frederik V, Ben M, et al. Multibody modeling of varying complexity for modal behavior analysis of wind turbine gearboxes[J]. Renewable Energy, 2011,36:3098-3113.

[19]Lethe G, Guyper , Kang J, et al. Simulating dynamics, durability and noise emission of wind turbines in single case environment[J]. Journal of Mechanical Science and Technology,2009,23:1089-1093.

[20]萬長森. 滾動軸承的分析方法[M]. 北京: 機械工業(yè)出版社, 1987.

[21]Shabana A A. Dynamics of multibody systems (third edition)[M]. Cambridge University Press, 2005.

[22]Cheon G J, Parker R G. Influence of manufacturing errors on the dynamic characteristics of planetary gear systems[J]. KSME International Journal, 2004, 18(4): 606-621.

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數(shù)模型及應用
p150Glued在帕金森病模型中的表達及分布
函數(shù)模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 午夜精品福利影院| 99re精彩视频| 色丁丁毛片在线观看| 青青草原偷拍视频| 亚洲乱亚洲乱妇24p| 国产免费久久精品99re不卡 | 日本久久久久久免费网络| 亚洲免费播放| 72种姿势欧美久久久大黄蕉| av一区二区无码在线| 国产精品亚洲一区二区三区在线观看| 欧美国产成人在线| AV不卡国产在线观看| 亚洲成A人V欧美综合天堂| 女同久久精品国产99国| 国产在线拍偷自揄拍精品| 久久亚洲综合伊人| 日韩在线视频网| 色香蕉网站| 亚洲精品午夜天堂网页| 亚洲有无码中文网| 久久久亚洲色| 久久毛片基地| 国产成人精品三级| 欧美亚洲欧美| 欧美在线国产| 人妻中文字幕无码久久一区| 国产呦视频免费视频在线观看| 怡春院欧美一区二区三区免费| 亚洲香蕉久久| 国产精品久久久久久影院| 欧美在线网| 中字无码av在线电影| 最新无码专区超级碰碰碰| www.av男人.com| 高清无码手机在线观看| a级免费视频| 精品无码一区二区三区在线视频| 亚洲国产成人麻豆精品| 国产乱人免费视频| 视频在线观看一区二区| 日韩欧美国产精品| 国产自在自线午夜精品视频| 亚洲成A人V欧美综合天堂| 国产三级韩国三级理| 99精品国产自在现线观看| 不卡视频国产| 91蜜芽尤物福利在线观看| 在线精品视频成人网| 日本黄色a视频| 乱码国产乱码精品精在线播放| 日日拍夜夜操| 99资源在线| 亚洲国产日韩欧美在线| 欧美第二区| 欧美日韩综合网| 狠狠色丁香婷婷| 亚洲国产成人精品无码区性色| 亚洲欧美天堂网| 亚洲精品在线影院| 99在线观看国产| 萌白酱国产一区二区| 国产欧美日韩91| 国产一级在线观看www色 | 久久国产精品国产自线拍| 久久精品人妻中文视频| 国产制服丝袜91在线| 免费av一区二区三区在线| 美女被躁出白浆视频播放| 免费一级大毛片a一观看不卡| 成人国产免费| A级全黄试看30分钟小视频| 88国产经典欧美一区二区三区| 中文字幕一区二区人妻电影| 国产第八页| 欧美日一级片| 欧美色视频日本| 久久婷婷色综合老司机| 久草热视频在线| 欧美一区二区三区不卡免费| 制服丝袜一区| 精品自拍视频在线观看|