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

多參量耦合的電主軸熱特性建模及分析

2016-12-23 00:47:29康躍然史曉軍高建民李法敬
西安交通大學學報 2016年8期
關鍵詞:變形模型

康躍然,史曉軍,高建民,李法敬

(西安交通大學機械制造系統(tǒng)工程國家重點實驗室,710054,西安)

?

多參量耦合的電主軸熱特性建模及分析

康躍然,史曉軍,高建民,李法敬

(西安交通大學機械制造系統(tǒng)工程國家重點實驗室,710054,西安)

為了深入研究電主軸熱特性,通過分析軸承非線性力學特性與熱效應間的耦合作用機制,建立了綜合考慮熱誘導預緊力和黏溫效應的軸承動態(tài)生熱計算模型;在引入接觸熱阻的基礎上,基于系統(tǒng)內部多參量耦合關系,建立了電主軸熱-結構耦合計算方法。以某型號電主軸為試驗平臺,進行了溫升及熱變形試驗,利用試驗數(shù)據(jù)對所建模型和分析方法進行了驗證。最后,采用所建模型,對電主軸進行了熱特性分析,結果表明:不同工況下電主軸仿真結果與試驗數(shù)據(jù)吻合度良好,所建模型和分析方法進一步提高了仿真計算精度;電主軸溫度場分布表現(xiàn)為外冷內熱;軸芯熱積聚導致的電主軸軸向熱伸長是影響加工精度的主要因素;受到配合方式和裝配位置的影響,軸承在溫升過程中預緊力不斷減小。

熱誘導預緊力;黏溫效應;熱-結構耦合計算;熱伸長

電主軸是現(xiàn)代高速、高精度數(shù)控機床的核心功能部件,其性能直接決定加工精度與效率。隨著電主軸零部件設計、制造和裝配等技術的不斷進步,諸如幾何裝配、刀具磨損及伺服控制等引起的誤差所占比例逐漸減小,而在高速、高精度極端加工條件下,電主軸熱變形導致的加工誤差則日益凸顯。大量研究與生產實踐表明,在高速、高精密加工中,熱變形引起的誤差占總誤差高達60%~80%[1]。因此,建立精確的電主軸熱特性分析模型,深入研究不同運行工況條件下電主軸的溫升和熱變形機理,對進一步發(fā)展電主軸技術無疑有著重要的理論意義和工程應用價值。

文獻[2]借助有限元分析軟件對電主軸進行了溫度場的仿真模擬,通過試驗驗證了有限元方法分析電主軸熱特性的可行性。文獻[3]通過對轉速為40 000 r/min的高速電主軸進行熱分析及試驗研究,獲得了在不同轉速下的溫度場及熱變形。文獻[4]建立了一種考慮系統(tǒng)熱響應和預緊方式影響的角接觸球軸承熱-機耦合動力學模型,分析了不同運行狀態(tài)下主軸軸承的摩擦損耗及動態(tài)支承剛度。文獻[5]通過構建電主軸瞬態(tài)熱網絡模型,獲得了軸系溫度場及熱參數(shù)的瞬態(tài)變化特性,并結合試驗驗證了瞬態(tài)模型的準確性。

國內外學者對電主軸熱特性分析模型進行了一定研究,但考慮電主軸溫度場與熱變形耦合機理的分析模型仍有待深入研究,以進一步提高模型準確性。為此,本文通過分析軸承非線性力學特性與熱效應間的耦合作用機制,建立了軸承動態(tài)生熱計算模型。在引入接觸熱阻模型的基礎上,基于系統(tǒng)內部多參量耦合關系,建立了電主軸熱-結構耦合計算方法,并通過試驗數(shù)據(jù)與有限元計算結果分析對比,對所建模型和分析方法進行了驗證。最后,采用所建模型對電主軸熱特性進行了分析。

1 電主軸熱特性建模

1.1 軸承動態(tài)生熱計算模型

電主軸軸承摩擦生熱受到軸承結構、外加負載和潤滑油黏度等多個復雜因素共同影響,在系統(tǒng)達到熱平衡之前,各因素間是動態(tài)耦合相互作用的。基于Palmgren整體生熱模型,將軸承摩擦力矩分為內、外圈分量Mij、Moj[6]

(1)

式中:Z為軸承滾動體個數(shù);Db為滾動體直徑;di、do分別為內、外圈滾道與滾動體接觸點處的直徑;dm為軸承節(jié)圓直徑;n為軸承轉速;f0為與軸承類型和潤滑方式有關的系數(shù);ν是潤滑油的運動黏度;P1i、P1o為軸承內、外圈摩擦力矩的計算負荷;f1為與軸承類型和所受負載有關的系數(shù)。

由于潤滑油本身的黏溫效應,其對軸承生熱的影響如圖1所示,隨著潤滑油溫度升高,其運動黏度逐漸降低,進而導致軸承生熱量不斷減少。以軸承常用的32#機油為例,其不同溫度下的運動黏度可表示為

(2)

式中:T為潤滑油溫度。

圖1 潤滑油黏溫效應對軸承生熱的影響

此外,P1i、P1o和f1均與軸承負載受力有關,以f1為例

(3)

式中:C0為軸承基本額定靜載荷;P0為軸承當量靜載荷,可由下式獲得

P0=0.5Fa+0.46Fr

(4)

其中Fa、Fr分別為軸承的軸向及徑向載荷,它們需要根據(jù)軸承不同的配置形式做具體的分析計算,如表1所示。

軸承組承受總載荷表達式為

(5)

式中:α為軸承實際接觸角;m為軸承個數(shù)。

對于采用定位預緊的軸承,其溫升過程中的熱變形會造成軸承內、外圈相對幾何關系的改變,由此產生的熱誘導預緊力,反過來會進一步影響電主軸熱特性。忽略軸承徑向和滾動體變形以及由于離心力帶來的軸承內、外接觸角不相等,假設軸承只發(fā)生剛體位移和局部彈性變形,可建立如下模型。

表1 軸承組受力分析

注:Fao為軸承運行過程中的實際預緊力;Fae、Fre分別為軸承組承受的軸向及徑向載荷;Fap為軸承組承受的總載荷。

(6)

圖2 軸承接觸角變化

參考軸承內部載荷計算理論,通過軸向位移與接觸角的關系式,采用Newton-Raphson法進行迭代,求解出軸承軸向位移δa與實際變形后軸承接觸角α的關系式[7]

(7)

式中:B為總曲率。將求得的接觸角α代入式(8),可求得軸承熱變形后的實際預緊力

(8)

軸承在運轉時,滾動體受到外滾道控制,滾動體與內圈之間存在自旋摩擦力矩[4]

(9)

式中:μSi為滾動體與內圈滾道的滑動摩擦因數(shù);aij為滾動體與內圈滾道的赫茲接觸橢圓長半軸;L(K)ij為滾動體與內圈滾道赫茲接觸橢圓的第二類完全積分。

由此可獲得單個滾動體與軸承內、外圈接觸區(qū)的摩擦生熱為

(10)

式中:ωsij為滾動體的自旋角速度。

進一步根據(jù)Burton的理論[8],忽略保持架和導引套圈對熱的影響,將上述軸承摩擦生熱量按1∶1的比例平均分配到滾動體和與其接觸的內、外圈接觸滾道上,從而有

(11)

式中:Hi、Ho分別為軸承內、外圈滾道表面生熱;Hb為滾動體表面生熱。

1.2 電主軸內部接觸熱阻

電主軸各零部件在微觀上并未完全接觸,存在著一定的空隙,由此產生的接觸熱阻會使通過的熱流產生一部分損失,進而影響電主軸溫度場分布及熱變形結果。通過在關鍵部位建立接觸熱阻計算模型可進一步提高所建模型精度。

(1)軸承滾動體與內外圈接觸熱阻R由接觸面的形狀和大小決定[9]

(12)

式中:ψ為與接觸面大小有關的幾何參數(shù);λ1、λ2分別為軸承內、外圈和滾動體的熱導率;a為軸承赫茲橢圓接觸區(qū)域的長半軸。

(2)軸承內圈與軸芯是過盈配合,其熱傳導率hs與零件及介質熱導率、結合面加工參數(shù)有關

(13)

式中:Lg為未接觸空間厚度;A為名義接觸面積;Ac為實際接觸面積;Av為未接觸部分的面積;λa及λf分別為軸芯熱導率和空隙間介質的熱導率。A、Ac、Av均需根據(jù)接觸理論,綜合考慮粗糙表面的微觀形貌與微凸體變形來求解獲得[10]。

(3)軸承外圈由于與軸承座采用間隙配合,其熱傳導率與軸承幾何結構、配合方式及材料物性參數(shù)等因素有關

(14)

式中:hring為軸承外圈的厚度;hgap為軸承外圈與軸承座的平均間隙;λair為空氣的熱導率。

1.3 電主軸熱-結構耦合計算方法

電主軸在達到熱平衡前,系統(tǒng)內部多參量間相互耦合,且各參量是隨溫度、時間變化的非線性函數(shù)。基于上述計算模型,通過有限元仿真軟件Ansys與Matlab的交互協(xié)同,建立電主軸熱-結構耦合計算方法,其算法流程如圖3所示。

圖3 電主軸熱-結構耦合計算方法流程圖

電主軸熱-結構耦合算法的思路為:基于上文所建軸承生熱及接觸熱阻計算模型,以及電主軸內置異步電機熱損耗[11]和散熱邊界條件進行計算[12],獲得初始熱源及熱邊界的解析值,然后加載至有限元模型;按一定規(guī)律設置電主軸瞬態(tài)溫升過程中的時間節(jié)點,先通過有限元分析求得第一個時間節(jié)點的溫度場分布和內、外圈的相對熱變形,進而在Matlab中對換熱系數(shù)、潤滑液黏度、預緊力及接觸熱阻等熱參量進行修正;然后,將修正后的熱參量作為下一個時間節(jié)點的熱源及邊界條件,加載至有限元分析模型上進行下一個時間節(jié)點的仿真計算;如此迭代計算,直到系統(tǒng)某一時刻的特征溫度與上一時刻的差值小于0.2 ℃,即可認定系統(tǒng)達到熱平衡。

2 電主軸熱特性模型驗證

在建立的電主軸試驗平臺上,測量某型號磨削電主軸在不同轉速工況條件下的溫升和軸端熱變形[13],利用試驗數(shù)據(jù)對所建模型和分析方法進行驗證。該電主軸額定功率為12.4 kW,最高轉速為24 000 r/min,其內置電機為三相異步交流電機,軸承串聯(lián)配置,前軸承B7008定位預緊,后軸承B7007定壓預緊,油霧潤滑。在電主軸前軸承1外圈布置pt100型溫度傳感器(見圖4),該傳感器精度為A級、誤差為(0.15±0.000 2|t|) ℃,通過預緊彈簧保證其與目標測試面的緊密接觸,保證測溫迅速準確。試驗時采集儀實時記錄軸承溫升數(shù)據(jù),直到系統(tǒng)達到熱平衡為止。在軸芯前端加裝電渦流位移傳感器來測量軸端面的熱變形(見圖5)。試驗前依據(jù)被測面材料特性對傳感器進行標定。由于被測端面存在一定的跳動量,所以需要對試驗數(shù)據(jù)進行高頻濾波處理,保留低頻下隨溫升變化的熱變形量,以保證試驗數(shù)據(jù)測量精度。

圖4 軸承溫度傳感器布置

圖5 軸端熱變形測量裝置

試驗工況如下:定子冷卻水套流量為16 L/min,環(huán)境溫度為23 ℃,系統(tǒng)軸端空載,軸承中度預緊力為300 N,油霧潤滑流量為10 L/s。共設定6 000、8 000、10 000、12 000、15 000 r/min 5組不同轉速工況。每組試驗結束后,確保將電主軸系統(tǒng)冷卻至室溫,再進行下一組試驗。

針對目標對象特征,建立三維簡化幾何模型,如圖6所示。在有限元仿真中,將采用定壓預緊方式的軸承外圈與軸承座內圈結合面設為摩擦接觸類型,軸承內、外圈與滾動體接觸設置為不分離狀態(tài),對于采用螺釘連接的結合面,算法設置為MPC。基于上文所建模型,采用熱-結構耦合計算方法對電主軸進行有限元分析。

圖6 電主軸簡化模型

圖7為電主軸不同轉速下軸承外圈瞬態(tài)溫升的計算結果和試驗數(shù)據(jù)對比,可見兩者溫升趨勢吻合較好。表2為所建模型與文獻[14]模型的計算誤差對比,本文所建模型將穩(wěn)態(tài)時溫升最大誤差由15.5%減小為6.4%。

圖7 不同轉速下軸系瞬態(tài)溫升仿真與試驗對比

轉速/r·min-1文獻誤差/%本文誤差/%600015.55.280008.16.4100008.61.01200011.61.3150008.46.3

圖8為不同轉速工況下電主軸軸端熱變形的計算和試驗結果對比。由圖可見,隨著轉速增加,電主軸軸端熱變形不斷增大,不同工況下試驗與仿真數(shù)據(jù)最大誤差為7.8%,與文獻[14]中的模型相比,本文所建模型的計算精度有了進一步的提升。

圖8 兩種模型電主軸熱變形仿真與試驗結果對比

3 電主軸熱特性分析

3.1 電主軸穩(wěn)態(tài)溫度場

采用所建模型和方法,分析了該電主軸在特定工況條件(轉速為10 000 r/min,空載)下的溫度場分布,如圖9所示。系統(tǒng)電機轉子處的溫度最高為63.42 ℃,其次是前、后軸承,而殼體部分由于靠近冷卻水套溫升最低,在7 ℃左右。由于冷卻水套只能作用于相鄰部分的電機定子和殼體,無法實現(xiàn)系統(tǒng)內部,特別是軸芯處積聚熱的快速引出和高效冷卻,從而使整個電主軸形成了外冷內熱的溫度分布格局。由于存在接觸熱阻,軸承外圈與軸承座間出現(xiàn)了明顯的溫度分布不連續(xù)。

圖9 電主軸內部溫度場

3.2 電主軸熱變形

將上文求解獲得的系統(tǒng)溫度場結果作為體載荷導入靜力學分析中。在零件螺栓連接面和軸承定位處施加位移約束,對整個系統(tǒng)施加重力場,求解獲得電主軸熱變形,結果如圖10、圖11所示。

圖10 電主軸軸芯軸向熱變形

圖11 電主軸徑向熱變形

在熱穩(wěn)態(tài)條件下,電主軸軸芯兩端軸向熱變形較大,后端熱伸長最大達19.2 μm,雖然熱特性分析顯示該區(qū)域的溫升并不是最大的,但由于軸芯熱積聚,使其軸向方向溫度梯度變化較大(最大溫差達23 ℃左右),再加之各零部件間的配合關系進一步限制了軸向方向的自由度,熱變形逐漸累積疊加至前、后端。電主軸徑向熱變形最大為15.6 μm,由于系統(tǒng)本身的熱對稱設計,變形沿中軸線對稱分布,對加工精度影響相對較小。

3.3 電主軸瞬態(tài)熱特性分析

圖12為電主軸各典型區(qū)域最高溫度隨時間的變化情況。在系統(tǒng)運行初始階段,各區(qū)域溫度急速上升,而由于彼此熱環(huán)境的差異,各部分溫升速率不同。待系統(tǒng)穩(wěn)定運行一段時間后,在油霧、冷卻水套等裝置連續(xù)、有效的運行作用下,系統(tǒng)生熱與換熱之間趨于平衡,最終達到熱穩(wěn)態(tài)。

圖12 電主軸各區(qū)域溫度變化

前軸承預緊力隨時間的變化如圖13所示。在電主軸運行初期,對于采用串聯(lián)配置的前軸承,由于僅考慮內、外圈的相對軸向熱位移,且其變化方向與軸承初始預緊方向相反,產生的熱誘導預緊力使軸承預緊力急劇減小,前軸承2由于裝配位置臨近電機,溫升及預緊力變化較大,而隨著系統(tǒng)逐漸達到熱穩(wěn)態(tài),軸承預緊力的變化也趨于平緩。

圖13 軸承預緊力隨時間變化

4 結 論

(1)建立了綜合考慮熱誘導預緊力和黏溫效應的軸承動態(tài)生熱計算模型,在引入接觸熱阻計算模型的基礎上,基于系統(tǒng)內部多參量耦合關系,建立了電主軸熱-結構耦合計算方法。

(2)有限元計算結果與試驗數(shù)據(jù)對比顯示:不同轉速工況下,兩者瞬態(tài)溫升趨勢一致;熱穩(wěn)態(tài)時軸承外圈溫升偏差最大為6.4%,軸向熱變形誤差最大為7.8%,所建模型的計算精度進一步提高。

(3)采用所建模型和計算方法,對電主軸進行了熱特性分析。結果表明:電主軸表現(xiàn)外冷內熱的溫度分布;軸芯熱積聚導致的電主軸軸向熱伸長是影響加工精度的主要因素;受到配合方式和裝配位置的影響,軸承在溫升過程中預緊力不斷減小。

[1] OKAFOR A C, YALCIN M. Derivation of machine tool error models and error compensation procedure for three axes vertical machining center using rigid body kinematics [J]. International Journal of Machine Tools & Manufacture, 2000, 40(8): 1199-1213.

[2] JIN K C, DAI G L. Thermal characteristics of the spindle bearing system with a gear located on the bearing span [J]. International Journal of Machine Tools & Manufacture, 1998,38(9): 1017-1030.

[3] JIAN L, DONG H K. A study on the thermal characteristics and experiments of high-speed spindle for machine tools [J]. International Journal of Precision Engineering and Manufacturing, 2015, 16(2): 293-299.

[4] 陳小安, 劉俊峰, 合燁, 等. 高速電主軸熱態(tài)性能及其影響 [J]. 機械工程學報, 2013, 49(11): 135-142. CHEN Xiao’an, LIU Junfeng, HE Ye, et al. Thermal properties of high speed motorized spindle and their effects [J]. Journal of Mechanical Engineering, 2013, 49(11): 135-142.

[5] 米維, 閆柯, 吳文武, 等. 考慮熱-變形耦合的主軸-軸承系統(tǒng)瞬態(tài)熱特性分析 [J]. 西安交通大學學報, 2015, 49(8): 52-57. MI Wei, YAN Ke, Wu Wenwu, et al. Transient thermal property analysis for spindle-bearing system considering thermo-deformation coupling [J]. Journal of Xi’an Jiaotong University, 2015, 49(8): 52-57.

[6] 康輝民, 陳小安, 陳文曲, 等. 高速電主軸軸承熱分析與實驗研究 [J]. 機械強度, 2011, 33(6): 797-802.

KANG Huimin, CHEN Xiao’an, CHEN Wenqu, et al. High speed motorized spindle bearing thermal analysis and experimental research [J]. Journal of Mechanical Strength, 2011, 33(6): 797-802.

[7] HARRIS T A, KOTZALAS M N. 軸承技術的基本概念 [M]. 羅繼偉, 馬偉, 等譯. 北京: 機械工業(yè)出版社, 2010: 138-140.

[8] BURTON R A, STAPH H E. Thermally activated seizure of angular contact bearing [J]. ASLE Trans, 1967, 10(4): 408-417.

[9] BERND B, JAY F. A thermal model for high speed motorized spindles [J]. International Journal of Machine Tools & Manufacture, 1999, 39(9): 1345-1366.

[10]劉志峰, 馬澄宇, 趙永勝, 等. 基于接觸熱阻的高速精密電主軸熱特性分析 [J]. 北京工業(yè)大學學報, 2016, 42(1): 17-23. LIU Zhifeng, MA Chengyu, ZHAO Yongsheng, et al. Analysis of high speed motorized spindle thermal characteristics based on thermal contact resistance [J]. Journal of Beijing University of Technology, 2016, 42(1): 17-23.

[11]BOSSMANNS B, TU J F. A power flow model for high speed motorized spindles: heat generation characterization [J]. Journal of Manufacturing Science and Engineering, 2001, 123(3): 494-505.

[12]MA Chi, YANG Jun, ZHAO Liang, et al. Simulation and experimental study on the thermally induced deformations of high-speed spindle system [J]. Applied Thermal Engineering, 2015, 86: 251-268.

[13]王夢茜. 機床電主軸熱特性分析技術研究 [D]. 西安: 西安交通大學, 2014.

[14]WANG Mengxi, HONG Jun, WU Wenwu, et al. An improved model of motorized spindle for transient thermo-structural analysis [C]∥IEEE International Symposium on Assembly and Manufacturing. Piscataway, NJ, USA: IEEE, 2013: 91-96.

(編輯 杜秀杰)

Modeling and Analyzing Multi-Variable Coupling Therma Characteristics for Motorized Spindle

KANG Yueran,SHI Xiaojun,GAO Jianmin,LI Fajing

(State Key Laboratory for Manufacturing Systems Engineering, Xi’an Jiaotong University, Xi’an 710054, China)

Considering the coupling behaviors between the non-linear mechanics character and the thermal effect, a dynamic heat calculation model for rolling bearings taking heat-induced preload and temperature-viscosity effect into account is established. A thermo-structural coupling calculation method introducing the thermal contact resistance is developed according to the coupling parameters in the motorized spindle system. The temperature rising and thermal deformation tests are carried out to verify the accuracy of this model and algorithm, and the thermal characteristics of the motorized spindle are analyzed. It is found that the simulation results coincide well with the experimental data under different working conditions to indicate the effect of the proposed model on improving simulation precision. The temperature field distribution of the spindle trends high inside and low outside. The thermal expansion of spindle caused by the heat accumulation of shaft mainly affects the manufacturing accuracy. During the temperature rising process, the preload of bearings decreases due to the fitting manner and the assembly position.

heat-induced preload; temperature-viscosity effect; thermo-structural coupling calculation; thermal expansion

10.7652/xjtuxb201608006

2016-03-06。 作者簡介:康躍然(1991—),男,碩士生;史曉軍(通信作者),男,副教授。 基金項目:國家機床重大專項資助項目(2014ZX04001191)。

時間:2016-05-17

http:∥www.cnki.net/kcms/detail/61.1069.T.20160517.1902.012.html

TH113

A

0253-987X(2016)08-0032-06

猜你喜歡
變形模型
一半模型
重要模型『一線三等角』
談詩的變形
中華詩詞(2020年1期)2020-09-21 09:24:52
重尾非線性自回歸模型自加權M-估計的漸近分布
“我”的變形計
變形巧算
例談拼圖與整式變形
會變形的餅
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
主站蜘蛛池模板: 国产欧美专区在线观看| 精品国产Av电影无码久久久| 99激情网| 国产午夜精品鲁丝片| 国产黑丝视频在线观看| 日本欧美一二三区色视频| 一级毛片免费不卡在线| 高清不卡毛片| 婷婷六月综合网| 中国国语毛片免费观看视频| 久久久精品无码一区二区三区| 91麻豆国产视频| www.国产福利| 精品国产网| 中文字幕1区2区| 香蕉精品在线| 黄色在线网| 国产69精品久久久久妇女| www.狠狠| 欧美a在线看| 免费aa毛片| 婷婷在线网站| 日本福利视频网站| 亚洲无码37.| 国产欧美亚洲精品第3页在线| 日本91视频| 欧美激情视频一区| 国产性精品| 99re热精品视频国产免费| 欧美人与牲动交a欧美精品 | 国产第一色| 欧美另类图片视频无弹跳第一页| 69av免费视频| 国产第一页屁屁影院| 国产中文一区二区苍井空| 四虎免费视频网站| 一级毛片无毒不卡直接观看| 国产精品免费入口视频| 亚洲性日韩精品一区二区| 在线观看国产黄色| 国产精品熟女亚洲AV麻豆| 国产在线观看人成激情视频| 日韩黄色精品| 午夜三级在线| 国产人成在线观看| 国产乱子伦一区二区=| 亚洲成aⅴ人在线观看| 国产成人无码AV在线播放动漫| 久久99国产综合精品1| 日韩无码视频网站| 欧美亚洲另类在线观看| 亚洲一级毛片免费观看| 国产成人麻豆精品| 久久精品视频一| 欧美日韩国产在线播放| 日韩欧美中文| 又爽又黄又无遮挡网站| 欧美在线一二区| 久久天天躁狠狠躁夜夜躁| 中文字幕亚洲综久久2021| 成人免费视频一区二区三区 | 波多野结衣一区二区三视频| 色天堂无毒不卡| 国产香蕉在线视频| 国产精品福利尤物youwu| 91在线播放免费不卡无毒| 激情国产精品一区| 国产在线观看高清不卡| 99re热精品视频国产免费| 久久精品无码国产一区二区三区| 亚洲天堂2014| 99资源在线| 中文字幕久久亚洲一区| 国产福利免费在线观看| 在线观看免费人成视频色快速| 久久男人资源站| 久久伊人操| 日韩av电影一区二区三区四区 | 99er这里只有精品| 国产精品对白刺激| 四虎影视无码永久免费观看| 亚洲高清中文字幕在线看不卡|