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

基于時變嚙合剛度的修形齒輪承載接觸分析*

2024-04-29 05:49:28韓振華周醞晨石萬凱單文桃曹忠亮
組合機床與自動化加工技術 2024年4期
關鍵詞:分配有限元

王 浩,韓振華,周醞晨,石萬凱,單文桃,曹忠亮

(1.江蘇理工學院機械工程學院,常州 213001;2.南京華秦光聲科技有限責任公司,南京 210034;3.重慶大學機械傳動國家重點實驗室,重慶 400044)

0 引言

輪齒承載能力是評價齒輪傳動品質的一項重要標準,齒間載荷分配系數和面接觸應力的模型精度是齒輪承載接觸力學計算結果精確性的前提與關鍵[1]。齒輪修形可以有效改善齒輪動載荷變化梯度,減少沖擊、振動與噪音,對提高齒輪傳動品質有十分重要的作用。然而,已有的齒輪承載接觸分析力學解析模型大都只能考慮鼓形的齒向修形,較難計入齒廓修形的影響[2-4]。

對于齒輪承載過程中產生的齒面接觸應力,學者常采用有限元法和解析法。其中,有限元法建模精度高,且計算結果精確,因而得到較為廣泛的應用。朱才朝[5]利用基于勢能法的三段函數式齒間載荷分配系數模型分析齒輪的彈流潤滑特性;AZNAR等[6]采用MPC算法與局部網格細化法優化齒輪的承載接觸有限元模型,提高了計算精度與效率;MAPER等[7]建立了包含齒廓修形的漸開線直齒圓柱齒輪承載接觸有限元模型;唐進元等[8]利用有限元分析了正交面齒輪的重合度、載荷分布、接觸應力等承載接觸性能參數;白恩軍等[9]建立了斜齒圓柱齒輪承載接觸有限元模型。然而,有限元算法計算精度依賴于模型的網格質量和局部網格細化程度,存在計算效率低、收斂難的問題。

相較于有限元法,基于赫茲彈性接觸理論的力學解析算法,計算快捷高效受到了越來越多的關注與研究。王會良等[10]基于齒輪TCA法和齒面柔度系數,分析齒輪承載接觸特性,建立了齒輪LTCA理論;曹雪梅等[11]提出一種簡化輪齒接觸分析中非線性方程數量的分解算法,提高了計算效率;王羽達等[12]基于輪齒TCA模型和赫茲接觸理論,計算了漸開線齒輪的動態接觸應力。

綜上所述,傳統齒間載荷分配系數解析建模未考慮齒廓修形及其引起的單雙齒接觸區域變化,無法準確計算承載嚙合齒對的齒間載荷分配系數與齒面接觸應力,同時,LTCA方法大都僅在TCA分析階段考慮齒廓修形,而在齒面柔度系數或承載接觸分析時采用有限元方法。為此,本文以齒廓修形的漸開線直齒輪為研究對象,提出修形齒輪單雙齒嚙合狀態的s判定方法與相應的嚙合剛度計算方法,建立更為精確的修形齒輪齒間載荷分配系數與齒面接觸應力解析算法,對不同修形高度、修形曲線與配對齒數條件下的齒輪算例進行承載接觸分析,并與同參數下的有限元結果進行對比,驗證解析算法的精度與穩定性。

1 齒廓修形的漸開線直齒輪模型

為建立齒廓修形漸開線直齒輪的齒廓曲線方程,需確定齒廓最大修形量Δmaxi、修形高度h、修形曲線,根據文獻[13]最大修形量Δmaxi為:

Δmaxi=ci+0.04FtB-1

(1)

式中:Δmaxi、ci分別為最大修形量與修形修正量(i為1、2,用以區分參與嚙合的主/從動輪),相應的主/從動輪修形修正量c1、c2分別取9與4[13],Ft為接觸點的圓周力,B為齒輪寬度。

建立坐標系XdOYd,如圖1所示,以Δx作為增量等距離散嚙合線,并在圓周方向上計算出輪齒1上的齒廓對應點,若該點位于修形齒廓區域內,嚙合線離散點的齒廓對應點處的修形量Δj為:

圖1 基于主動輪的嚙合線等距離散示意圖

Δj=Δmaxi(x/L)k

(2)

式中:修形曲線變化系數k由選取的修形曲線確定[14],Linear、Yoshio、Walker與拋物線修形曲線變化系數k分別為1.0、1.2、1.5與2.0,x為嚙合點M在嚙合線方向上距離單雙齒變化點A的距離,L為輪齒2的修形區域在高度方向上對應至嚙合線A點至C點的線段長度(長修形),C點為輪齒1的嚙合起始點;另外,根據文獻[15],正弦修形曲線的嚙合點修形量為:

Δj=Δmaxisin(πx/L)

(3)

確定齒廓對應點處的修形量Δj后,根據漸開線齒廓曲線展成原理,可建立修形齒廓曲線參數方程為:

(4)

式中:u為嚙合線離散點對應漸開線齒廓處的參角,u=θ+αt,θ為嚙合點處展角,αt為壓力角。

2 修形齒輪單雙齒嚙合狀態判定與嚙合剛度

標準漸開線齒輪副在承載嚙合傳動過程中齒輪嚙合剛度主要包括3部分:輪齒剛度(彎曲剛度kb、剪切剛度ks與徑向壓縮剛度ka)、齒輪接觸剛度kc與齒輪基體剛度kf。單齒嚙合剛度Ki為上述各剛度的串聯形式,雙齒嚙合剛度Kall為兩對輪齒嚙合剛度的并聯形式[16-17]。

(5)

(6)

式中:i為1、2,用于區分參與嚙合的不同齒對,下標g、e分別表示主動輪與從動輪。

當采用齒廓修形時,原單齒嚙合區間的齒輪嚙合剛度根據式(5)進行計算,原雙齒嚙合區間由于修形齒廓參與嚙合,嚙合狀態會發生變化,需要重新判定嚙合狀態:根據齒對嚙合力與嚙合剛度計算接觸點處沿嚙合線方向的變形量,然后與接觸點處修形量進行對比,基于二者的大小關系判定修形齒輪副的單雙齒嚙合狀態,并計算單雙齒的嚙合區間與嚙合剛度,其中,修形齒輪嚙合剛度定義為KTPM[18]。

如圖2a所示,當主動輪順時針轉動時,齒對1首先進入嚙合承擔載荷,由于齒對2主動輪處存在修形齒廓,當齒對1沿嚙合線方向的變形量δ1大于齒對2處主動輪的齒廓修形量Δ1時,齒對2將進入嚙合并參與傳動,此時兩對輪齒共同承擔負載扭矩Tc,并在嚙合線方向上的兩處作用點產生嚙合力F1、F2,根據力矩平衡和廣義胡克定律可得:

(7)

(a) 齒對2受力分析 (b) 齒對1受力分析

式中:K1、K2與δ1、δ2分別為齒對1、齒對2的單齒嚙合剛度與接觸點沿嚙合線方向的變形量,rb2為從動輪的基圓半徑。

與此同時,齒對1沿嚙合線方向的變形量δ1即為整體嚙合變形量δall,則δall=δ1>δ2,相應的修形齒輪嚙合剛度KTPM為:

KTPM=Fn/δ1,δall=δ1>δ2

(8)

式中:Fn為齒輪副法向力,根據式(7)將式(8)推導為:

KTPM=(K1δ1+K2δ2)/δ1,δall=δ1>δ2

(9)

同時,齒對2在嚙合點處嚙合變形量δ2與修形量Δj之和等于齒對1處嚙合變形量為δ1。

δ1=Δj+δ2

(10)

將式(10)帶入式(9),可得到:

KTPM=K1+K2-K2Δj(δ1)-1,δall=δ1>δ2

(11)

將式(8)轉化為關于變形量δ1的表達式,并代入式(11),可推導得出圖2a所示齒對承載接觸狀態下的修形齒輪嚙合剛度KTPM,即式(12)中條件為δall=δ1>δ2的表達式;同樣地,如圖2b所示的齒對承載接觸狀態,此時齒對2沿嚙合線方向的變形量δ2即為整體嚙合變形量δall,即δall=δ2>δ1,同理可推導出相應的修形齒輪嚙合剛度KTPM,即式(12)中條件為δall=δ2>δ1的表達式。因此,圖2所示齒對承載接觸狀態下的修形齒輪嚙合剛度KTPM可寫為:

(12)

式(12)即為齒廓修形漸開線直齒輪的嚙合剛度分析方法,下文將建立單雙齒嚙合狀態的判定模型,確定實際雙齒嚙合區間,然后依據式(12)計算雙齒區間內的齒輪嚙合剛度。

如圖3所示,基于主動輪齒廓,對齒輪副嚙合區間進行劃分,分析齒對1、2的接觸嚙合狀態及相應的嚙合剛度,點a1為嚙入點,點a3、a4為單雙齒嚙合變化點,點a2、a5為原雙齒接觸區域平分點,點a1~點a5將接觸區域劃分為A1、A2、A3、A4、A5。

(a) 接觸區域A1參與嚙合 (b) 接觸區域A2參與嚙合

如圖3a所示,當齒對1以a1點即將進入接觸區域A1參與嚙合時,齒對1中從動輪待接觸區域存在修形齒廓,若齒對2沿嚙合線方向的變形量δ2大于齒對1中從動輪修形量Δ2,則齒對1接觸參與嚙合,即為雙齒嚙合狀態,此時的修形齒輪嚙合剛度KTPM可依據式(12)計算;若齒對2沿嚙合線方向的變形量δ2小于齒對1中主動輪修形量Δ2,則齒對1不參與嚙合,此時為單齒嚙合狀態(齒對2單獨承載),齒對2的嚙合剛度K2依據式(5)進行計算。

如圖3b所示,當齒對1以a2點即將進入接觸區域A2參與嚙合時,齒對2中主動輪待接觸區域存在修形齒廓,若齒對1沿嚙合線方向的變形量δ1大于齒對2中主動輪修形量Δ1,則齒對2接觸參與嚙合,即為雙齒嚙合狀態,此時的修形齒輪嚙合剛度KTPM可依據式(17)計算;若齒對1沿嚙合線方向的變形量δ1小于齒對2中主動輪修形量Δ1,則齒對2不參與嚙合,此時為單齒嚙合狀態(齒對1單獨承載),齒對1的嚙合剛度K1依據式(5)進行計算。

如圖3c所示,當齒對1以a3點即將進入接觸區域A3參與嚙合時,此時接觸區域A3為單齒嚙合狀態,齒對1的嚙合剛度依據式(5)進行計算。

綜合圖3所述,即可得到齒廓修形齒輪在各嚙合位置的單雙齒嚙合狀態判定與相應嚙合剛度計算方法。

3 齒廓修形漸開線直齒輪承載接觸分析

3.1 齒間載荷分配系數

3.1.1 未修形齒輪的齒間載荷分配系數

根據式(6)確定的雙齒嚙合剛度Kall與齒輪整體變形量δall,法向力Fn可表示為:

Fn=Kallδall=(K1+K2)δall

(13)

根據式(7)中確定的載荷、剛度以及變形關系,帶入式(13)等式左側,可推導出:

K1(δall-δ1)=K2(δall-δ2)

(14)

由于嚙合齒對在嚙合線方向上變形量δ1、δ2中的較大值與雙齒嚙合整體變形量δall相等,此時無論δall=δ1≥δ2或δall=δ2≥δ1,根據式(14)皆可推導出各嚙合齒對在嚙合線方向上的變形量δ1、δ2與雙齒整體嚙合變形量δall相等,即:

δall=δ1=δ2

(15)

進一步可根據式(7)中嚙合力、剛度與變形的數學關系推導出:

Fi/Ki=Fn/Kall=δall

(16)

式中:Fi為接觸齒對的嚙合力(i為1、2)。

對式(16)進行比例互換,可得到未修形下的齒輪齒間載荷分配系數S為:

S=Fi/Fn=Ki/Kall

(17)

3.1.2 齒廓修形齒輪的齒間載荷分配系數

當齒廓修形后,原雙齒嚙合區域的嚙合狀態發生變化,根據標準漸開線齒輪時變嚙合剛度式(5)與齒廓修形齒輪嚙合剛度式(12),重新確定齒輪的單雙齒嚙合狀態與相應嚙合剛度,繼而根據式(17)計算修形齒輪不同齒廓區域參與嚙合的齒間載荷分配系數SLDF。根據圖3,以齒對1為例,分析齒對1承載嚙合下的齒間載荷分配情況,具體為:

(18)

式中:Ku為齒對1的嚙合剛度,當齒對1處于單齒嚙合狀態時,其嚙合剛度依據式(5)計算,即Ku=K1;當齒對1、2同時參與嚙合時,此時計算Ku需考慮齒對2嚙合剛度K2的疊加影響,即:

Ku=KTPM-K2

(19)

當齒對1在齒廓區域A1處嚙合時,由于齒對2接觸變形量δ2的影響,齒對1存在非接觸與接觸兩種狀態,若δ2>Δ2,則齒對1接觸,此時為雙齒嚙合狀態,將式(18)代入式(19)可得到齒對1的齒間載荷分配系數SLDF。

(20)

若δ2<Δ2,則齒對1為非接觸狀態,此時齒對1不分擔載荷,則齒間載荷分配系數為0,即:

SLDF=0

(21)

當齒對1在齒廓區域A2處嚙合時,由于齒對1接觸變形量δ1的影響,齒對2存在非接觸與接觸兩種嚙合狀態,若δ1>Δ1,則齒對1接觸,此時為雙齒嚙合狀態,齒對1的嚙合齒廓為非修形區域,其齒間載荷分配系數根據式(22)進行計算。

(22)

若δ1<Δ1,則齒對2為非接觸,此時齒對1承擔全部載荷,齒間載荷分配系數為1,即:

SLDF=1

(23)

當齒對1在齒廓A3處嚙合時,此時為單齒嚙合區域,齒間載荷分配系數為1,即:

SLDF=1

(24)

當齒廓進入A4、A5段參與嚙合時,齒對1皆為非修形齒廓參與嚙合,此時僅考慮齒對1的嚙合剛度,即Ku=K1,相應的齒間載荷分配系數為:

(25)

綜上所述,在一個嚙合周期內,根據齒對1在齒廓區域A1、A2、A3、A4、A5是否接觸參與嚙合及其嚙合狀態判定,計算參與嚙合齒對的嚙合剛度,根據式(18)~式(25),可得到修形漸開線直齒輪齒間載荷分配系數SLDF計算模型。

(26)

3.2 齒面接觸應力

根據上述建立的修形齒輪齒間載荷分配系數模型,可計算嚙合力Fi,繼而根據赫茲接觸理論計算齒面接觸應力。

首先,需確定不同齒廓赫茲接觸類型下接觸點處的綜合曲率半徑Re,針對未修形的齒廓區域,根據歐拉薩瓦利公式[19],通過計算沿嚙合線方向上基圓至嚙合點的距離求得其曲率半徑ρ1、ρ2,對于齒廓修形區域,以主動輪為例,根據微分幾何與修形齒廓方程式(4),可推導出修形齒廓曲率CTPM與曲率半徑ρTPM。

(27)

ρTPM=1/CTPM

(28)

β=u+tanu,rb1為主動輪基圓半徑。

針對齒廓修形齒輪承載嚙合時存在的不同赫茲接觸類型,可推導出相應的綜合曲率半徑Re。

(29)

通過式(26)確定的齒間載荷分配系數SLDF,可計算出參與嚙合齒對的接觸力Fi,根據式(29)所得綜合曲率半徑Re,基于赫茲彈性接觸理論,可推導出修形齒輪承載接觸時的齒面接觸應力σH與接觸半寬b。

(30)

式中:Ee為綜合彈性模量,具體表達式為:

(31)

式中:E1、E2與ν1、ν2分別為主、從動輪的彈性模量與泊松比。

如圖4所示為該解析算法的流程圖。

圖4 齒間載荷分配系數與齒面接觸應力計算流程圖

4 計算結果分析與討論

基于上述建立的算法,求解不同修形高度、修形曲線與配對齒數3種條件下不同算例的齒間載荷分配系數與齒面接觸應力,并與同參數下的有限元計算結果進行對比評價,驗證本文解析算法計算結果的準確性和穩定性。基于表1中拋物線修形齒輪的主要幾何參數,每種條件下選取5個算例求解以上不同齒廓修形齒輪算例的齒間載荷分配系數(圖5~圖7)和齒面接觸應力(圖12~圖14)。使用有限元法對比評價解析算法結果的準確性和穩定性,以驗證算法的正確性。

表1 齒輪參數

圖5 不同修形高度算例的齒間載荷分配系數

圖6 不同修形曲線算例的齒間載荷分配系數

圖7 不同配對齒數算例的齒間載荷分配系數

分析如圖5~圖7所示的齒間載荷分配系數結果可知:

(1)本文解析算法與有限元兩種計算方法的齒間載荷分配系數具有相同變化趨勢:由最小值逐漸增大至最大值1,此過程為雙齒嚙合區間;隨后進入單齒嚙合區間,即SLDF=1的水平線;再進入雙齒嚙合區,下降至最小值。同時嚙合線各點所對應的誤差數值較小。

(2)根據圖5~圖7數據結果,可得到不同齒廓修形算例在一個嚙合區間內的齒間載荷分配系數誤差均值,如圖8所示,每種條件下5個算例誤差均值變化范圍分別為7.5%~9.2%、7.5%~8.0%、5.6%~7.8%,表明誤差均值較小且不同算例的誤差均值相近。算例結果分析表明,對于不同算例情況下的齒輪副,本文解析算法的齒間載荷分配系數計算精度較高,且算法具有較好的誤差穩定性。

圖8 一個嚙合區間內的齒間載荷分配系數誤差均值

圖10 單齒嚙合區結束點位置誤差

(3)根據圖5~圖7數據結果,分析齒廓修形后每個算例的單齒嚙合區位置誤差,以此可進一步評價解析算法的計算精度。其中,齒廓修形后的單齒嚙合區起始點位置誤差、結束點位置誤差與嚙合區間位置誤差分別如圖9~圖11所示。不同修形高度、修形曲線與配對齒數算例的單齒嚙合區起始誤差變化范圍分別為5.8%~12.4%、5.1%~14.2%、5.8%~10.1%,單齒嚙合區結束點位置誤差變化范圍分別9.2%~13.1%、6.1%~13.3%、9.7%~12.6%,單齒嚙合區間位置誤差變化范圍分別為3.2%~12.5%、3.2%~13.0%、2.2%~8.8%,表明解析算法所得出的單雙齒變換點及其區間大小與有限元結果之間誤差較小,可較為準確地計算出單雙齒嚙合變換點以及嚙合區間位置。因此,本文解析算法中依據修形量與變形建立的單雙齒嚙合區間判定方法具有較高的計算精度,可較為準確的計算單雙齒變換點以及嚙合區間位置。

(4)另外,由于配對齒數的變化,齒輪副有效嚙合區間也會隨之變化,分析圖7中得到的橫坐標實際嚙合線長度可得到齒廓修形齒輪的有效嚙合區間,相對于有限元結果,5個配對齒數算例的有效嚙合區間誤差分別為1.3%、1.0%、3.0%、3.1%、1.1%,相對誤差較小,表明本文解析算法對不同配對齒數的齒輪副嚙合區間計算結果的準確度較高。

分析如圖12~圖14所示的4種條件下不同齒廓修形算例在一個嚙合區間內的齒面接觸應力結果可知:

圖12 不同修形高度算例的齒面接觸應力

圖13 不同修形曲線算例的齒面接觸應力

圖14 不同配對齒數算例的齒面接觸應力

(1)解析算法的齒面接觸應力結果與有限元結果在變化趨勢具有較好的一致性;根據圖12~圖14的齒面接觸應力結果,可得到3種條件下齒廓修形算例在一個嚙合區間內的齒面接觸應力誤差均值(圖15a)變化范圍分別為7.0%~9.0%、8.4%~10.4%、9.0%~10.7%,誤差均值較小,同時不同算例間誤差均值相近,表明解析算法具有較高的齒面接觸應力計算精度與誤差穩定性。

(a) 單嚙合區間 (b) 去除突變區域

(2)與此同時,在修形齒輪剛開始進入嚙合、即將退出嚙合以及修形曲線與漸開線過渡區域,有限元結果產生了局部應力突變,以圖12~圖14中的A、B區域為例,主要由以下兩個原因導致:①剛開始嚙入與即將退出嚙合時,接觸齒對嚙合力較小,需要極小網格尺寸(低于微米1~2個量級)才能有效捕捉到該處的接觸應力信息,建模困難,且嚙入時齒頂部位存在幾何尖點;②修形齒廓過渡區域對幾何曲率產生一定程度的嚙合畸變。因此,當求解步長位于上述位置附近時會直接影響有限元計算結果,導致有限元局部應力突變、計算誤差較大,這是不可避免的。除去由有限元結果突變區域所帶來的計算誤差后,解析算法的齒面接觸應力在一個嚙合區間內的誤差均值(圖15b)變化范圍分別為4.8%~6.4%、4.5%~6.8%、6.1%~8.6%,再次表明解析算法具有較高的齒面接觸應力計算精度與誤差穩定性。

5 結論

(1)本文針對齒廓修形的漸開線直齒輪,通過對比承載齒輪接觸點沿嚙合線方向的變形量與修形量,得到了修形齒輪單雙齒嚙合狀態的判定方法與相應嚙合剛度的計算方法,通過推導修形齒輪嚙合力、剛度與變形的力學關系,建立了修形齒輪各嚙合齒對精確的齒間載荷分配系數計算模型。

(2)基于齒間載荷分配系數計算了嚙合齒對實際承擔載荷,應用歐拉薩瓦利公式和赫茲接觸理論,建立了齒廓修形齒輪嚙合過程中的齒面接觸應力模型。

(3)不同修形高度、修形曲線與配對齒數條件下多種齒廓修形齒輪的算例結果表明:本文所提出解析算法的齒間載荷分配系數、齒面接觸應力與相應單雙齒嚙合區間結果與有限元方法結果變化規律相一致,相對誤差較小,具有較高的計算精度與穩定性,為齒廓修形漸開線直齒輪的承載接觸分析提供了一種高效便捷的新方法。

猜你喜歡
分配有限元
基于可行方向法的水下機器人推力分配
新型有機玻璃在站臺門的應用及有限元分析
上海節能(2020年3期)2020-04-13 13:16:16
應答器THR和TFFR分配及SIL等級探討
基于有限元的深孔鏜削仿真及分析
基于有限元模型對踝模擬扭傷機制的探討
遺產的分配
一種分配十分不均的財富
績效考核分配的實踐與思考
磨削淬硬殘余應力的有限元分析
基于SolidWorks的吸嘴支撐臂有限元分析
主站蜘蛛池模板: 99在线观看精品视频| 亚洲综合片| 国产微拍精品| 在线视频亚洲色图| 亚洲av日韩综合一区尤物| 欧美日韩一区二区三区四区在线观看 | 国产成人福利在线视老湿机| 国产在线观看第二页| 中字无码av在线电影| 国产va在线观看免费| 无码中字出轨中文人妻中文中| 2021亚洲精品不卡a| 色亚洲激情综合精品无码视频| 国产精品99在线观看| 国产麻豆va精品视频| 国产91特黄特色A级毛片| 伊人天堂网| 日韩小视频在线观看| 99国产精品免费观看视频| 久久这里只有精品国产99| 欧洲在线免费视频| 狠狠ⅴ日韩v欧美v天堂| 直接黄91麻豆网站| 国产成人综合亚洲欧美在| 国产精品无码作爱| 欧美一级夜夜爽| 国产丝袜无码精品| 免费在线视频a| 国产午夜一级毛片| 久久国产精品国产自线拍| 国产午夜一级毛片| 热久久这里是精品6免费观看| 欧美性久久久久| 亚洲一级色| 欧美黄网站免费观看| 精品少妇人妻无码久久| 91美女视频在线| 亚洲成肉网| 亚洲美女久久| 四虎成人在线视频| www.亚洲天堂| 国产99在线| 全裸无码专区| 国产成人精品亚洲77美色| 国产va欧美va在线观看| 亚洲日韩精品伊甸| 国产日韩欧美视频| 免费在线播放毛片| 亚洲黄色高清| 国产免费人成视频网| 国产精品3p视频| 日韩无码黄色网站| av一区二区三区在线观看 | 日韩少妇激情一区二区| 日本一区高清| 亚洲日本一本dvd高清| 人妻精品全国免费视频| www.99精品视频在线播放| 2021国产精品自拍| 波多野吉衣一区二区三区av| 国产高清在线丝袜精品一区| 国产乱子伦视频在线播放| 久久亚洲日本不卡一区二区| 欧美日韩福利| 国产亚洲精品资源在线26u| 欧美激情二区三区| 国产成人8x视频一区二区| 国产一级毛片高清完整视频版| 成人福利一区二区视频在线| 欧美激情首页| 青青青视频免费一区二区| 国产哺乳奶水91在线播放| 日本三级黄在线观看| 国产一级毛片yw| 色哟哟色院91精品网站| 亚洲人成网址| 亚洲热线99精品视频| 无码精油按摩潮喷在线播放 | 国产精品久久自在自2021| 在线精品亚洲一区二区古装| 制服丝袜亚洲| 国产无码高清视频不卡|