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

離心載荷對軸類部件過盈配合的影響研究

2023-09-01 07:38:02舒易亮劉志明楊廣雪高敬宇刁曉明
中南大學學報(自然科學版) 2023年7期
關鍵詞:有限元

舒易亮,劉志明,楊廣雪,高敬宇,刁曉明

(1. 北京交通大學 機械與電子控制工程學院,北京,100044;2. 中國鐵道科學研究院集團有限公司 機車車輛研究所,北京,100081)

過盈量是保證軸類過盈配合結構傳遞扭矩和載荷的重要參數,在工程結構件正常工作時,過盈量將直接影響結構的過盈配合性能。大多數工程機械緊固連接配合通過壓裝法或溫差法實現過盈裝配。然而,軸類過盈配合部件往往服役于高轉速環境,如鐵路機車輪軸[1]、高速旋轉主軸與轉子或軸承[2-4]、齒輪與軸[5]、渦輪發動機轉子[6]等。高速旋轉的軸類過盈配合部件在離心載荷的作用下,過盈配合的2 個部件在接觸面產生徑向位移差,從而導致旋轉部件過盈量減小[7-9]。而軸類過盈配合部件是通過過盈連接傳遞扭矩及軸向載荷的,當過盈量減小到機械設計允許的最小過盈量[10]時,將直接導致軸類過盈配合部件的性能失效。

國內外學者采用有限元方法對軸類過盈配合結構接觸面間的應力、位移,接觸壓應力及過盈量等參量的變化規律進行了大量的研究。ZHAO[11]通過有限元方法對實心和空心滾子軸承進行了多體接觸分析,研究了載荷在滾動軸承內部的分布規律。ZHANG 等[12]基于有限元方法對環形齒輪-輪轂過盈連接部件開展仿真研究,得到了比傳統Lame 方程更加完整、更加精確的過盈配合部位的應力。張松等[2]采用有限元方法對高速旋轉主軸的過盈連接特性進行了研究,發現主軸與轉子間的過盈連接特性與初始過盈量和旋轉角速度有關。當主軸轉速高于松脫速度時,主軸與轉子間處于分離狀態,兩者之間的接觸應力和扭矩下降為零。單寶峰等[13]對影響高速主軸與軸承之間過盈量的因素展開了分析,指出影響主軸軸承過盈量的因素主要有溫度、離心力及初始配合過盈量等。

由于過盈量是影響軸類過盈配合部件連接性能的一個關鍵參數,當軸類過盈配合部件承受周期性的旋轉彎曲載荷時,靠近接觸邊緣的局部接觸區域就會產生周期性的往復滑移,從而引起微動疲勞損傷[14]。楊廣雪等[15]指出,接觸壓力、微動滑移幅值、摩擦剪切應力等參數是影響微動疲勞損傷的重要參數,在離心載荷作用下,高速列車輪軸過盈量的變化將直接導致過盈配合面的微動損傷參量發生變化。ZHANG等[16-18]基于旋轉彎曲微動疲勞試驗和微動疲勞定量仿真模型,研究了過盈配合試樣的微動疲勞損傷,揭示了微動疲勞裂紋的表面形貌演化、磨損特征及損傷分布變化特征。ZENG等[19]提出了一種全尺寸過盈配合輪軸裂紋萌生的有限元預測模型,考慮了微動磨損引起的應力重分布對微動疲勞的影響,并提出了緩解輪軸過盈配合接觸邊緣部位應力集中的方法。陳剛等[20-24]開展了列車車軸過盈配合部位的微動疲勞試樣或者全尺寸實物輪軸的微動疲勞試驗及仿真研究,指出輪軸配合過盈量作為一個關鍵參量時刻影響車軸的微動疲勞的損傷。然而,實際輪軸過盈量在不同加載循環次數下的變化情況以及過盈量沿配合區域的軸向分布特征不能觀測到,而且現有的研究并沒有考慮過盈量在不同旋轉角速度下的減小效應對車軸微動疲勞的影響。因此,研究軸類過盈配合部件過盈量在不同離心力載荷作用下的減小效應,同時考慮初始過盈量和離心載荷的影響,建立過盈量隨角速度變化的定量分析公式模型是非常有必要的。

軸類部件過盈配合接觸仿真研究是屬于高度非線性的復雜問題,配合面間呈現復雜的接觸狀態和應力狀態。隨旋轉角速度增大,離心載荷增大,將導致軸類過盈配合部件的過盈量的減小效應逐漸增強,甚至可能引起軸類過盈配合部件的性能失效。孫林平等[8]對離心載荷作用下實心簡化輪軸的過盈量的變化進行了有限元仿真分析,并結合有限元仿真結果分析研究了接觸壓應力的變化規律,指出在離心載荷對輪軸過盈配合有重要影響,當轉速達到一定量級時,會出現輪軸分離現象。龐強宏等[9]采用Ansys Workbench 軟件分析了圓盤-軸和齒輪-軸在離心力和高溫同時作用下的接觸應力和等效應力分布規律,綜合考慮過盈連接件可能面臨的復雜工況,使設計的過盈配合結構性能更加可靠。盧萍等[7]仿真分析了列車輪軸在高速運行情況下的輪軸連接性能,指出離心載荷會降低輪轂孔面的接觸壓力。然而,研究者僅僅對軸類過盈配合部件在離心力載荷作用下過盈量的減小效應進行了研究,推導公式時未考慮初始過盈量的影響,也未將仿真分析結果和理論計算結果進行詳細對比,更未分析離心載荷作用下輪軸過盈量的定量變化規律和特征。

本文作者將離心載荷對軸類過盈配合部件連接性能的影響研究簡化成平面應力問題,推導含初始過盈量的軸類過盈配合部件在不同旋轉角速度下的過盈量變化公式。考慮初始過盈量的影響,定量地分析研究離心載荷對軸類過盈配合部件過盈量的影響。基于簡化的輪軸實物等效模型,將推導的理論公式進行數值模擬,結合有限元仿真分析軸類過盈配合部件在不同旋轉角速度下的過盈量及接觸壓應力的變化特征,以驗證所推導理論公式的準確性。

1 輪軸過盈配合計算方法

1.1 過盈配合接觸壓力

輪軸過盈配合仿真計算屬于接觸非線性問題,傳統設計方法以厚壁圓筒為模型,采用彈性力學的方法計算過盈連接所需要的過盈量,保證在合理的過盈量下能傳遞扭矩和軸向力。一般可將軸類過盈配合部件看作組合厚壁圓筒,進而分析厚壁圓筒過盈面間變形和應力的關系。如圖1 所示。假設組合圓筒中的內圓筒的內半徑為R1,共同配合半徑為R2,外圓筒的外半徑為R3,過盈配合的長度無限長,不考慮軸向的邊界條件,組合圓筒的半徑過盈量為δ,內、外圓筒的彈性模量分別為E1和E2,泊松比分別為υ1和υ2。根據過盈量δ和輪軸配合前、后的變形協調條件,輪軸過盈配合的接觸應力P可由Lame方程[25]求出:

圖1 組合厚壁圓筒簡圖Fig. 1 Sketches of combined thick-walled cylinder

若輪軸材料屬性相同時,則輪軸彈性模量E=E1=E2,泊松比υ=υ1=υ2,接觸應力P可進一步簡化為

當內圓筒為實心軸時(即R1=0時),接觸應力P可進一步簡化為

然而,在實際工程問題中,軸類過盈配合部件的輪軸的材料往往不一致,材料力學性能也各有差異。為準確獲得輪軸在初始過盈量下的接觸壓應力,有必要采用式(1)進行精確計算。若軸類過盈配合部件的輪軸的材料一致或者屬性相近,可采用式(2)進行計算,也可以得到接近實際的應力。

1.2 離心載荷下過盈量的減小效應

如圖1所示,若2個裝配的圓柱體處于含初始過盈量的接觸狀態,則在高速旋轉離心載荷作用下,內圓筒和外圓筒將出現徑向位移。對于軸類過盈配合部件,軸的長度遠比輪的厚度大,當輪軸結構在高速旋轉時,可以簡化成含初始接觸壓應力的2個等厚旋轉圓盤。對于任意的等厚旋轉圓盤,密度為ρ,若以均勻角速度ω旋轉,則圓盤內半徑為r的單位圓柱的離心力Fd大小相同,方向沿著各自的徑向方向,且與φ無關,可以看作軸對稱問題[26],如圖2 所示。圖2 中,u為微小單元的徑向位移,r為微小單元相對圓心位置的半徑,φ為力的夾角,σφ為旋轉角為φ時的應力,σr為徑向應力。

圖2 等厚旋轉圓盤受力示意圖Fig. 2 Schematic diagrams of force of rotating disc of equal thickness

2 輪軸等效模型建立及數值模擬

2.1 實物輪軸模型參數

高速列車輪軸往往依據普速既有列車輪軸的設計經驗選擇輪軸過盈配合量,標準TB/T 1718.3—2018[27]中規定動車組輪軸間的過盈量為0.10%dm~0.15%dm+0.06 mm(dm為車軸輪座直徑),故將輪軸間的半徑過盈量取為0.15 mm,接近過盈量的中間值。真實車軸輪座區域的直徑為200 mm,軸身直徑為173 mm,車軸輪座直徑與軸身直徑之比為1.156。空心車軸內孔直徑為30 mm,車軸輪座區域左右兩側的突懸量分別為1.635 mm 和3.5 mm,兩側端部過渡圓弧弧長分別為25 mm 和16 mm,卸荷槽半徑和深度分別為16 mm和1 mm,實物輪軸模型如圖3所示。考慮到模型的對稱性,選取高速列車實物輪軸模型的1/4模型,通過有限元仿真軟件計算不同初始過盈量下接觸壓應力在車軸輪座區域的軸向分布,為了獲得接觸邊緣區域的精確應力,在接觸邊緣區域進行局部細化網格,采用有限元數值求解進行網格劃分的最小單元長度和寬度均為400 μm[28],有限元網格模型如圖4 所示,仿真計算時在對稱面施加對稱邊界條件。

圖4 輪軸1/4對稱有限元網格模型Fig. 4 Symmetric finite element mesh model of 1/4 wheel axle

2.2 輪軸等效模型建立及仿真方法

采用有限元軟件計算輪軸實物1/4模型在不同過盈量下的接觸壓應力沿軸向相對位置的分布特征,如圖5 所示,軸向正方向定義為有A端指向B端,并以A端為橫坐標零點位置。由圖5可知:輪軸過盈配合區域的接觸壓應力與過盈量呈線性正相關關系,相同過盈量下分布呈現中部平緩、數值小,接觸邊緣應力急劇增大的特征,與文獻[15]中試樣軸過盈配合部件的接觸壓應力分布特征一致。實物輪軸過盈配合區域接觸邊緣左側的接觸壓應力比右側邊緣的接觸壓應力小,這是因為車軸輪座區域左側突懸量比右側突懸量小,且左側過渡圓弧比右側過渡圓弧大。車軸輪座區的微動疲勞損傷與接觸區域的接觸壓應力、相對滑移幅值及外界振動載荷密切相關,而車軸輪座區內側邊緣的接觸壓應力遠比輪座區外側邊緣的接觸應力大,故從接觸壓應力角度看,車軸輪座區域內側比車軸輪座區域外側更容易出現微動疲勞損傷[28]。

圖5 不同過盈量下接觸壓應力沿軸向分布Fig. 5 Distributions of contact compressive stress along axial direction with different interferences

將實物輪軸簡化成等效輪軸模型,通過式(1)結合仿真結果中車軸輪座中間區域的接觸壓應力數據,取平均得到的等效接觸壓應力Pave,進而求解出車輪的等效外徑R3。考慮到輪軸過盈配合中間區域的接觸壓應力與Lame方程的壓應力計算值較接近[29],故選取過盈量為150 μm、軸向相對位置為60~120 mm 的模型有限元節點的接觸壓應力進行求解,如下式所示:

由于輪軸材料彈性模量相差不大,為了計算方便,可認為輪軸彈性模量相等。取E=E1,υ=υ1,R1和R2分別參考圖3 中圖4 中輪軸1/4 對稱有限元網格模型實物輪軸參數,分別取值為15 mm 和100 mm。根據式(2)可計算車輪的等效外徑R3,如式(21)所示,同時考慮到建模方便,將計算得到的R3向上取整為230 mm。

為驗證車輪等效外徑的準確性,基于式(1)分別計算過盈量分別為100、115、130、150、180和210 μm 時的接觸壓應力,與仿真求解車軸輪座中間區域的節點接觸壓應力平均值進行對比,如表1所示,實物輪軸過盈配合中間區域的有限元求解節點平均值與采用等效車輪外徑基于式(1)計算的結果相對誤差均小于1%,證實了車輪等效外徑R3計算結果的合理性。

表1 不同過盈量下仿真平均值與Lame公式計算值比較Table 1 Comparisons of simulation average value and Lame formula calculation value with different interferences

基于高速列車實物輪軸模型過盈配合區域的尺寸參數及求解的車輪等效外徑R3,通過有限元前處理軟件Hypermesh進行幾何前處理建模,建立輪軸過盈配合等效模型。軸部件的內徑R1為15 mm,輪部件的外徑R3為230 mm,輪軸共同接觸區域的半徑R2為100 mm,模型參數如圖6所示。輪軸過盈配合部件的軸的彈性模量為E1為2.05×105MPa,泊松比為υ1為0.3,輪的彈性模量E2為2.1×105MPa,泊松比為υ2為0.3。采用1/4 軸對稱模型進行仿真分析,同時為了在接觸區域獲得精確的有限元仿真結果,采用精細化網格技術,在接觸區域的單元長度為ESmin為1 mm。

圖6 輪軸等效有限元模型Fig. 6 Equivalent finite element model of wheel and axle

在輪軸過盈配合部件的輪軸間建立接觸對,軸的接觸面為接觸從面,輪的接觸面為接觸主面,法向接觸行為采用硬接觸模擬,主從面節點保證一一對應,接觸面間的摩擦因數取0.6,切向行為采用罰函數接觸算法模擬,輪軸過盈配合部件的主從面間的過盈量取為0.15 mm。仿真求解在ANSYS軟件中求解,共分為2個載荷步:第1個載荷步施加初始過盈量CONF為0.15 mm,同時為避免輪軸主從接觸面節點存在初始的間隙,設定初始間隙調整量為0.001 mm,時間步長為2 s;第2 個載荷步施加全局的旋轉角速度0~1 207.73 rad/s,時間步長為20 s,表示在20 s時間內,旋轉角速度由0 rad/s線性增大到1 207.73 rad/s。

2.3 數值模擬分析

根據前文對離心載荷下過盈配合輪軸過盈量的變化理論公式推導結果,在MATLAB 進行數值模擬仿真。理論計算的過盈接觸壓應力、軸的徑向位移(uaxle)和輪的徑向位移(uwheel)及輪軸接觸徑向位移差隨旋轉角速度的變化如圖7所示。由圖7可知:當初始過盈量0.15 mm 時,接觸壓應力為124.62 MPa,隨旋轉角速度增大,接觸應力不斷減小,當旋轉角速度增至966 rad/s 時,接觸壓應力減小至0 MPa。

圖7 理論仿真接觸參數隨旋轉角速度的變化Fig. 7 Theoretical simulation of contact parameters as a function of rotational angular velocity

由含初始過盈的輪軸理論公式可知,軸部件外側接觸面和輪部件內側接觸面在徑向方向會發生相對位移,其相對徑向位移之差與過盈量相同。在離心載荷作用下,軸外表面和輪內表面的徑向位移及接觸過盈量的變化如圖7 所示。從圖7 可見:對于初始過盈量為0.15 mm的輪軸部件,當軸部件外側接觸面的徑向過盈位移為-0.05 mm,輪部件內表面的徑向過盈位移為0.10 mm。旋轉角速度小于966 rad/s 時,輪軸接觸面的徑向過盈位移和接觸過盈隨旋轉角速度增加不斷減小,呈非線性特征;當旋轉角速度等于966 rad/s 時,軸部件和輪部件的徑向過盈位移及接觸過盈位移均變為0 mm,接觸壓應力減小為0 MPa,表明此時輪軸過盈配合剛好由過盈配合狀態轉變為間隙配合狀態,隨旋轉角速度繼續增大,輪軸接觸間的間隙不斷增大,間隙配合狀態時的輪軸間隙的增長速率比過盈配合狀態的輪軸過盈量的減小速率大,這是由于在過盈配合狀態時,輪軸間存在相互作用約束,而在間隙配合狀態時,輪軸接觸面相互分開,輪軸間不再存在相互作用。

3 有限元仿真分析

3.1 含初始過盈量的接觸壓應力仿真

采用有限元軟件ANAYS對含初始過盈量下輪軸間接觸壓應力進行仿真分析,對1/4輪軸等效模型施加平面對稱約束,之后施加過盈量0.15 mm。由式(1)可計算出過盈配合輪軸的在初始過盈量下的接觸壓應力理論值約為124.62 MPa。含初始過盈量的接觸壓應力有限元仿真結果云圖如圖8 所示,仿真的過盈接觸壓應力沿軸向分布特征和公式理論計算結果比較如圖9 所示。由圖8 和圖9 可以看出:理論公式計算的接觸壓應力與仿真結果近似相等,且在接觸邊界處偏差較大,這是因為有限元仿真在接觸邊界存在一定的應力集中,故可以認為有限元仿真的結果與理論值比較吻合。

圖8 過盈接觸壓應力仿真結果Fig. 8 Simulation results of interference contact pressure stress

圖9 過盈接觸壓應力沿軸向分布Fig. 9 Distributions of compressive stress of interference contact along axial direction

3.2 離心載荷下的過盈量仿真

采用ANAYS有限元軟件對離心載荷作用下輪軸過盈量的減少量進行仿真分析,旋轉角速度ω的取值范圍為0~1 207.73 rad/s,并在仿真時間步20 s內呈線性增大。選取輪軸軸向位置中間處的節點作為研究過盈接觸參數的對象,ANSYS 中接觸參數隨旋轉角速度的變化如圖10 所示,其中uaxle、uwheel及Δuaxle_wheel分別為軸的徑向位移、輪的徑向位移及輪軸徑向位移差。初始接觸壓應力和初始過盈量分別為131.14 MPa和0.15 mm,且接觸壓應力及基礎過盈量隨角速度增大而減小,呈現非線性特征。當旋轉角速度為953 rad/s 時,接觸壓應力減小至0 MPa,此時,過盈量減小至0 mm。若角速度繼續增大,則輪軸配合由過盈配合狀態轉變為間隙配合狀態,與MATLAB 中的理論數值模擬結果基本一致。

圖10 仿真接觸參數隨旋轉角速度的變化Fig. 10 Changes of simulated contact parameters with rotational angular velocity

ANSYS 仿真計算的過盈量在旋轉角速度處于0~953 rad/s變化期間,輪軸接觸面的徑向位移差一直為0.15 mm,數值上與過盈量相等,這與MATLAB 中數值模擬過盈量隨角速度增大而減小的現象有所不同。這是因為在ANSYS有限元軟件中,含初始過盈量的過盈配合是基于有限元軟件算法通過輪軸接觸面的徑向位移施加實現的,所以,當施加初始過盈量時,輪的徑向位移為-0.05 mm,軸的徑向位移為0.10 mm。在旋轉角速度從0 rad/s增大到953 rad/s時,輪軸的徑向位移呈非線性同步增加。但是,徑向位移差一直保持為0.15 mm 不變,在此過程中輪軸處于過盈配合狀態。接觸面的位置是輪軸過盈配合面的共同接觸面處,軸部件的外表面和輪部件的內表面一直處于接觸狀態,并未發生分離。故此階段的相對位移差一直與初始過盈量的狀態一致。當角速度大于953 rad/s 時,輪軸過盈配合由過盈配合狀態轉變為間隙配合狀態,輪軸徑向位移差將繼續增大,增長斜率基本和輪的內表面徑向位移增長斜率基本一致。

3.3 有限元分析與理論模擬對比

含初始過盈量的輪軸過盈配合部件在離心力載荷作用下的接觸壓應力變化,如圖11 所示。由圖11 可知:隨旋轉角速度增大,理論數值模擬和有限元仿真計算得到的接觸壓應力隨旋轉角速度的變化曲線基本一致,都呈現非線性特征。理論數值模擬及有限元仿真計算的徑向位移差變化曲線如圖12 所示。為了將理論數值模擬和有限元仿真計算的過盈量變化結果進行對比,將理論數值模擬的過盈量的變化曲線沿縱軸正向增大0.15 mm得到過盈量偏置曲線,該曲線在過盈量大于0.15 mm后的變化趨勢和ANSYS仿真曲線變化趨勢一致。

圖11 數值模擬及仿真接觸壓應力變化對比Fig. 11 Comparisons of numerical simulation and simulation contact pressure stress change

圖12 數值模擬過盈量變化及仿真徑向位移差對比Fig. 12 Numerical simulations of interference change and comparison of simulated radial displacement difference

由圖12 可見:有限元軟件采用給定接觸主從面徑向位移差的方式施加初始過盈量。在過盈配合狀態轉變為間隙配合狀態之前,輪軸接觸面一直處于零間隙狀態,故認為接觸面間的徑向位移差不變。當輪軸接觸狀態轉變為間隙配合狀態時,輪軸過盈配合的間隙在0.15 mm 的基礎上繼續增加。即當輪軸配合處于間隙配合狀態時,有限元仿真得到的徑向位移差為初始過盈量和間隙值之和。綜上所述,對于含過盈量的輪軸過盈配合結構,在離心力載荷作用下仿真結果與理論數值模擬結果一致,驗證了含初始過盈量的軸類過盈配合部件在離心載荷作用下的過盈量減小量公式推導的準確性。

4 結論

1) 不同旋轉角速度下的過盈量及接觸壓應力的理論數值模擬結果與有限元仿真分析結果一致,驗證了含初始過盈量的軸類過盈配合部件的過盈量減少量公式的準確性。

2) 隨旋轉角速度從0 rad/s 增大到1 207.73 rad/s,輪軸間接觸壓應力及過盈量不斷減小,呈現非線性特征。當旋轉角速度增大到960 rad/s 左右時,接觸壓應力及過盈量減小為零,此時,輪軸配合狀態由過盈配合狀態轉變為間隙配合狀態,輪軸間接觸壓應力及過盈量的理論數值模擬結果與有限元仿真結果一致。

3) 由于有限元軟件采用對接觸主從面施加徑向位移差的方式施加初始過盈量,在過盈配合狀態轉變為間隙配合狀態之前,輪軸接觸面一直處于零間隙狀態。在過盈量減小為零之前,輪軸過盈配合部件的接觸面間的徑向位移差一直不變。

猜你喜歡
有限元
基于擴展有限元的疲勞裂紋擴展分析
非線性感應加熱問題的全離散有限元方法
TDDH型停車器制動過程有限元分析
新型有機玻璃在站臺門的應用及有限元分析
上海節能(2020年3期)2020-04-13 13:16:16
基于I-DEAS的履帶起重機主機有限元計算
基于有限元模型對踝模擬扭傷機制的探討
10MN快鍛液壓機有限元分析
磨削淬硬殘余應力的有限元分析
基于SolidWorks的吸嘴支撐臂有限元分析
箱形孔軋制的有限元模擬
上海金屬(2013年4期)2013-12-20 07:57:18
主站蜘蛛池模板: 国产精品网址在线观看你懂的| 九九视频免费看| 国产成人1024精品下载| 亚洲国产AV无码综合原创| 人妻无码中文字幕一区二区三区| 亚洲av色吊丝无码| 色AV色 综合网站| 香蕉国产精品视频| 国产精品久线在线观看| 有专无码视频| 日韩小视频在线观看| 国产福利一区视频| 国产精品亚洲片在线va| 欧美一级高清免费a| 亚洲国产系列| 亚洲午夜片| 成人福利在线免费观看| 亚洲成a人在线播放www| 日韩天堂在线观看| 老熟妇喷水一区二区三区| 欧美国产三级| 亚洲天堂成人在线观看| 日本在线免费网站| 久久99国产综合精品1| 色噜噜综合网| 成人免费网站在线观看| 无码一区18禁| 色天天综合久久久久综合片| 亚洲国产天堂久久九九九| 黄色网页在线观看| 亚洲一区二区视频在线观看| 日韩毛片免费视频| 国产精品成人免费视频99| 免费 国产 无码久久久| 日本午夜精品一本在线观看| 中文字幕乱妇无码AV在线| 蜜臀av性久久久久蜜臀aⅴ麻豆| 亚洲免费黄色网| 久久 午夜福利 张柏芝| 99久久精品美女高潮喷水| 欧美日韩精品在线播放| 久视频免费精品6| 亚洲精品第五页| 色亚洲激情综合精品无码视频| 亚洲日韩久久综合中文字幕| 亚洲91在线精品| 熟妇人妻无乱码中文字幕真矢织江| 婷婷六月在线| 99视频在线免费| 日本黄色a视频| 中文无码精品A∨在线观看不卡 | 一级毛片网| 成人综合网址| 日韩二区三区无| 国产啪在线91| m男亚洲一区中文字幕| 亚洲欧美在线综合一区二区三区 | 污污网站在线观看| 在线免费无码视频| 日韩高清中文字幕| 久久99久久无码毛片一区二区| 2021最新国产精品网站| 91日本在线观看亚洲精品| 亚洲欧洲日韩久久狠狠爱| 五月婷婷中文字幕| 97人妻精品专区久久久久| 白丝美女办公室高潮喷水视频| 无码aaa视频| 色婷婷天天综合在线| 波多野结衣久久精品| 日本人又色又爽的视频| 亚洲一区二区精品无码久久久| 青青热久免费精品视频6| 欧美日韩导航| 91蝌蚪视频在线观看| 亚洲激情99| 伊人中文网| 欧美成人精品高清在线下载| 日韩av无码DVD| 日韩欧美中文字幕一本| 亚洲欧洲日韩综合色天使| 无码精品一区二区久久久|