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

150MD24Z7.5高速電主軸多場耦合模型與動態性能預測

2016-07-26 02:21:34張麗秀吳玉厚
振動與沖擊 2016年1期

張麗秀, 閻 銘, 吳玉厚, 陸 峰

(沈陽建筑大學 交通與機械工程學院,沈陽 110168)

?

150MD24Z7.5高速電主軸多場耦合模型與動態性能預測

張麗秀, 閻銘, 吳玉厚, 陸峰

(沈陽建筑大學 交通與機械工程學院,沈陽110168)

摘要:電主軸是機電一體化產品,充分考慮并預測其動態特性是機床主軸系統優化設計的前提。基于電主軸內部磁場、電場、溫度場、結構場間的耦合關系,建立了150MD24Z7.5高速電主軸多場耦合有限元模型,通過電主軸電機電磁損耗及軸承摩擦生熱計算,仿真電主軸溫度場及結構場變化,討論電主軸熱態特性與振動特性之間耦合關系,分析電主軸溫升熱膨脹后氣隙變化對振動特性的影響并通過實驗加以驗證。研究結果表明,電主軸溫升形變對振動幅值影響較大,其中由氣隙變化引起的電磁力幅值增加12.1%。利用該多場耦合模型可預測電主軸振動幅值,預測誤差為10.2%。

關鍵詞:電主軸;耦合;預測;動態性能

評價電主軸質量的指標中有很大一部分屬于其動態性能。電主軸的動態性能包含溫升、熱變形、熱應力、振動、噪聲及動態剛度等。若能在設計過程中提前預測電主軸動態性能并進行合理優化,對提高電主軸整體質量具有重要意義。電主軸是典型的機電一體化產品,其動態特性又與電磁損耗及軸承摩擦密不可分[1],因此,電主軸的優化必然涉及機、電、熱、磁等多種學科,采用有限元方法分析其動態特性是現代學者常用的方法之一。Holkup等[2]研究了電主軸中電機與軸承引起的瞬間溫度變化可能導致軸承等部件的損壞。由于電主軸的轉子和轉軸之間的過盈配合量直接影響機床的加工精度,芮執元等[3]研究了熱變形對電主軸轉子與轉軸過盈配合的影響。Huang等[4]設計了轉軸熱變形補償控制裝置,提高加工公差,減少機床加工誤差和時間。Creighton等[5]分析了高精度電主軸的熱誤差對工件進行微加工時的影響。鑒于上述原因,于翔[6]通過對影響轉軸熱變形的主要熱源及因素分析,并得出電主軸生熱散熱的重要公式為減少轉軸熱變形和提高加工精度提供理論依據。Abele等[7]提出關于電主軸的熱-動力學研究,同時對電主軸的發展等各方面做了總結。并且何俊等[8]進一步利用有限元分析軟件研究了電主軸系統熱-結構耦合模型,得到電主軸系統的溫度場分布和熱變形,通過實驗驗證模型的正確性。謝黎明等[9]采用仿真軟件研究電主軸在轉動中產生的熱量導致的轉軸變形量。文獻[10]通過分析劃片機氣靜壓電主軸的發熱和傳熱,用 ANSYS 軟件建立熱-應力有限元分析模型并計算主軸的溫度場分布及轉軸的熱變形值。Zverev 等[11]針對主軸系統熱-結構特性的研究,建立高速電主軸的彈塑性變形單元模型。文獻[12-14]運用有限元模型法來預測高速電機主軸熱特性。但關于熱量的輸入采用傳統的理論估計算法,并未考慮到電主軸實際溫度場分布不均勻的情況,考慮因素單一,因此對于轉軸形變的計算并不準確。本文綜合考慮電機損耗生熱以及軸承摩擦生熱兩大因素,結合電主軸的獨特散熱機構,建立電主軸熱態模型,計算結構模型,分析電主軸電磁場,溫度場,熱變形場以及熱應力場,并由此預測耦合因素下電主軸動態性能。

1電主軸耦合特性

圖1 電主軸內部的機-電-熱-磁耦合關系Fig.1 Mechanical-electric-thermal-magnetic coupling relationship of motorized spindle

高速電主軸運行過程中,會由于各部分的損耗而產生大量的熱,這些熱量的產生將使電主軸產生熱效應并嚴重影響電主軸的剛度、壽命及精度特性。Jedrzejewski等[15]研究的高速加工中心電主軸由于溫度場導致的主軸前端變形情況,研究結果表明,當電主軸轉速增大時,主軸前端的變形可達80 μm左右,達到穩定至少需要200 s時間。可以看出,電主軸溫升導致的熱效應不可忽視。由于高速旋轉的電主軸轉子系統周圍存在溫度場、磁場、電場及力場的多場耦合,電主軸內部存在復雜的機電熱磁耦合問題。電主軸內部的機-電-熱-磁耦合關系如圖1所示。從圖中可以看出,由于電主軸運行過程中發熱,導致內部零件產生熱變形,即電主軸的結構發生變化,主軸電機內電磁特性更新,最終影響電主軸的動力學行為。

2電主軸在機-電-熱-磁耦合狀態下動態性能預測方法

根據電主軸的機-電-熱-磁耦合關系,采用有限元方法,建立電主軸的機-電-熱-磁耦合模型,是分析預測其動態性能的有效方法。本文在電主軸驅動系統模型、電機模型、結構模型、溫度場模型及熱膨脹模型的基礎上,通過有限元計算方法,預測電主軸動態性能。首先,建立驅動系統與電機模型相耦合的時步有限元模型,計算電主軸電磁損耗;通過電主軸結構模型,分析其軸承在相應轉速下的摩擦損耗。其次,建立電主軸溫度場模型,結合電磁損耗與軸承摩擦損耗,對電主軸溫度場進行預測分析。第三,建立電主軸熱膨脹模型,計算電主軸零件各部分熱變形及熱應力。最后,在有限元軟件中導入電主軸結構變形,利用電主軸結構模型繼續進行振動特性分析。圖2為電主軸在機-電-熱-磁耦合狀態下動態性能預測策略。

圖2 機-電-熱-磁耦合狀態下動態性能預測策略Fig.2 Dynamic performance prediction strategy

3電主軸動態預測有限元模型

電主軸在運轉過程中由于損耗的作用會使其生熱,進而溫升形變會影響電主軸內部結構,其中最為主要的兩大生熱源為電機的電磁損耗熱以及軸承的摩擦熱。由于電主軸獨特結構,無法采用傳統風扇散熱。電主軸散熱主要通過是冷卻液對流換熱以及定轉子之間氣隙對流換熱。圖3為電主軸生熱和換熱結構圖。本文以150MD24Z7.5電主軸為研究對象建立電主軸的機-電-熱-磁有限元模型,電主軸的基本參數如表1所示,電主軸前軸承型號為7009C,后軸承型號為7007C。

表1 電主軸參數

圖3 電主軸生熱和換熱結構圖Fig.3 Heat and heat transfer structure

3.1電主軸軸承摩擦生熱有限元模型

摩擦生熱是高速電主軸軸承生熱的主要原因,除了球與軌道的差動摩擦生熱,還有球的陀螺轉動摩擦生熱以及球與保持架的滑動摩擦生熱[16]。在有限元軟件中在耦合場中用三維六面體二十節點Solid226單元建立有限元模型,仿真時模型中電主軸在轉速10 000 r/min空載運行,電主軸仿真的初始溫度20℃。將轉速轉化為切向位移載荷施加在轉軸上,軸承外圈固定約束,后軸承端面施加250 N預緊力。 電主軸網格劃分如圖4所示,將軸承轉速以載荷的形式施加到有限元模型上得到圖5基于軸承摩擦熱的轉軸溫度所示的電主軸轉軸溫升,將仿真計算所得轉軸的溫度作為電主軸溫度場仿真分析的一項熱源載荷。

圖4 軸承摩擦熱加載后網格圖Fig.4 Grid under load of bearing friction heat

圖5 基于軸承摩擦熱的轉軸溫度場Fig.5 Shaft temperature based on bearing friction heat

3.2電主軸電磁耦合與損耗計算

3.2.1時步有限元仿真模型

由文獻[17-18]可知SPWM供電對電機損耗有直接影響。為了精確計算電主軸正常工況下的電磁損耗,將在有限元軟件中建立的電主軸模型導入建立的變頻器外電路中,形成如圖6所示的時步有限元仿真模型。

圖6 電主軸時步有限元模型Fig.6 Time-step finite element model

3.2.2電主軸電磁損耗計算

利用時步有限元模型計算得到如圖7所示的電主軸磁場計算電磁損耗,圖8為電主軸10 000 r/min空載運行時各項電磁損耗。

圖7 電主軸電磁場Fig.7 Magnetic field

圖8 電主軸電磁損耗Fig.8 Electromagnetic loss

3.3電主軸傳熱系數

3.3.1定轉子氣隙對流換熱系數計算

定轉子氣隙對流換熱系數計算較復雜,當定轉子氣隙處于純層流狀態時,則為純導熱不具備散熱功能,換熱與電主軸的轉速無關。本文考慮定轉子氣隙非純層流狀態,電主軸定轉子氣隙對流換熱系數按照式(1)計算[19]:

hg=28(1+ωg0.5)

(1)

式中:ωg為氣隙平均風速(m/s),ωg為0.5轉子圓周速度(m/s)。

由式(1)可知電主軸定轉子氣隙的換熱系數與轉速有關,定轉子氣隙換熱系數hg在10 000 r/min轉速下計算得148 W/(m2·℃)。

3.3.2冷卻液對流換熱系數計算

對流換熱系數與流體物性、流動形態及換熱面的表面狀況等因素有關。根據流體相似性準則,有

Nμ=0.023Re0.8Ρr0.4(紊流狀態)

(2)

式中:ν為流速(1.4 m/s);D為特征尺寸,

μ為流體運動粘度;Re為流體雷諾數;λf為流體的導熱系數,W/(m·K);η為流體的動力粘度, Pa·s;Cp為恒壓熱容, J/(Kg·K);Pr為流體的普朗特數;Ρμ為流體的努塞爾特數。

按式(2)計算電主軸在溫度20℃時,定子鐵心與冷卻液的對流換熱系數α計算約為2 600 W/(m2·K)。

4電主軸動態性能仿真結果

4.1電主軸溫度場仿真分析

4.1.1溫度仿真條件

假設仿真計算所得損耗全部轉化為熱量,不考慮電主軸整體輻射散熱微弱影響;電主軸初始溫度20℃;電主軸轉子端面與外表面絕對光滑處理,所以不考慮空氣阻力摩擦損耗熱。最終由上面仿真出單純軸承摩擦熱影響下的轉軸溫度(66℃~72℃)作為電主軸一項熱源;將穩態下運轉下39.5~40.0 ms時間段的電磁損耗導入電主軸溫度場模型作為另一熱源。傳熱機制為理想傳熱狀態,定子絞線繞組與定子鐵心通過空氣接觸導熱;轉子導條與轉子鐵心是接觸導熱;轉子鐵心和軸承與轉軸為接觸導熱;定子鐵心與轉子鐵心除了通過氣隙接觸導熱外,還進行對流散熱,并且定子鐵心與外面冷卻液也進行對流散熱。

4.1.2溫度仿真結果

把電磁損耗及軸承摩擦生熱的仿真結果與換熱系數加載到有限元軟件中進行電主軸溫度場分析,獲得電主軸穩態溫度場如圖9所示,穩態熱流量場如圖10所示。從圖9電主軸溫度場分布看出最高溫度分布在電主軸轉軸部分,達到68℃,外表面由于冷卻液作用結果,溫度最低只有23.8℃,接近仿真溫度20℃。整個電主軸溫度場剖視度情況看溫度分布均勻,從內至外溫度呈現下降趨勢。從圖10電主軸熱流量場看出電主軸的定轉子部分與繞組絞線的傳熱最大,其中熱流量最大部分發生在定子齒部,最小部分是電主軸轉軸。這是由于定轉子氣隙換熱作用的結果。

圖9 電主軸溫度場Fig.9 Temperature field

圖10 電主軸熱流量場Fig.10 Heat flow field

4.2電主軸結構場仿真分析

4.2.1結構仿真條件

因為電主軸最外面是機殼,電主軸最外面采用固定約束,把電主軸理想為整體單元,根據電主軸力學材料特性如表2所示,在有限元軟件中繼續進行電主軸結構分析,將上面得到的熱載荷施加到電主軸模型中。

表2 電主軸力學材料特性

4.2.2結構仿真結果

得到電主軸溫升形變如圖11所示,電主軸熱應力如圖12所示。從圖11電主軸形變場可以看出,由于電主軸溫升引起的電主軸內部形狀變化導致電主軸的定轉子變化最嚴重,最大變形發生在定轉子齒頂部分,這是由于電主軸綜合生熱效果使電主軸定轉子氣隙兩側達到相對較高的溫度,而且由于定轉子氣隙兩側齒槽的特殊結構,使其容易發生形變,進而達到電主軸溫升形變后的最大變形量。

圖11 電主軸形變場Fig.11 Deformation field

圖12 電主軸應力場Fig.12 Stress field

從圖12電主軸熱應力場可以看出,熱應力集中在電主軸定轉子齒槽以及轉軸部分,這是因為定轉子齒槽特殊結構且變形量較大的原因,轉軸的溫度最高,并且因為電主軸在轉軸中磁場熱作用的綜合效果,最終產生不規則的熱應力。

4.3電主軸溫升形變前后電磁力分析

陳小安等[20]為了研究高速電主軸軸向振動對加工件品質的影響,采用有限元法通過建立高速電主軸轉子-軸承動力學模型分析系統的振動特性,但是忽視了熱變形對振動特性的影響,因為定轉子間氣隙作為電主軸電能與機械能轉換的唯一通道,其大小直接影響磁場的分布情況,進而引起不同階次的力波,導致電主軸產生電磁振動及噪聲。Chang等[21]研究了電主軸轉子熱膨脹引起氣隙動態變化對振動性能的影響。本文針對電主軸熱變形導致的氣隙形變量達到了0.023 2 mm,導致電主軸定轉子間氣隙磁密發生變化,最終影響電主軸的徑向電磁力。

根據Maxwell定律,由電主軸氣隙磁場產生的,作用在定子內表面單位面積上的徑向力的表達式為

(3)

通過有限元軟件仿真電主軸形變前后氣隙磁密空間分布如圖13所示。從圖中可知,由于熱膨脹導致定轉子間氣隙變小,使定轉子氣隙磁密增大。通過式(3)計算電主軸電磁力并進行FFT變換得到電主軸形變前后電磁力頻譜如圖14所示,從圖中看出電主軸溫升前后各階頻率的電磁力幅值均變大,在1 000 Hz時電磁力幅值增長了12.1%。

圖13 形變前后磁密對比圖Fig.13 Magnetic flux density contrast

圖14 形變前后電磁力頻譜對比圖Fig.14 Electromagnetic spectrum contrast

5電主軸動態性能實驗測試

5.1電主軸溫度場實驗分析

電主軸溫度測試系統如圖15所示,采用TC-2008多路溫度測試儀在電主軸10 000 r/min空載時運行5 300 s記錄電主軸外殼和電主軸定子繞組的溫升結果分別如圖16和圖17所示。從圖16外殼溫度監測數據曲線看出,在21.5℃的環境溫度條件下,電主軸外殼溫升4℃,這與圖9的3.8℃仿真溫升結果誤差為5%。從圖17繞組溫度監測數據曲線看出,在21.5℃的環境溫度條件下,電主軸繞組溫升28.5℃,這與圖9的28.357℃仿真溫升結果誤差為0.5%。溫升結果與仿真結果相符。

圖15 電主軸溫度測試系統Fig.15 Temperature testing system

圖16 電主軸溫度監測曲線Fig.16 The curve of temperature monitoring

圖17 電主軸繞組溫度監測曲線Fig.17 The curve of winding temperature

5.2電主軸振動實驗分析

在監測電主軸溫度變化情況的同時,電主軸振動測試系統如圖18所示,利用INV3018振動數據采集儀采集數據,其中壓電式振動傳感器測試點分布與溫度傳感器測試點相對應,利用加速度幅值衡量振動劇烈程度,并且使用DAXP軟件對電主軸運行中采集的數據進行頻譜分析,并得到電主軸形變前后頻譜分析對比圖如圖19所示。 從圖中可以清晰看出電主軸各階頻率的振動幅值在電主軸溫升形變后有所上升,一方面是因為溫升導致電主軸內部零件發生熱變形,即各零件的形狀發生變化,這可能影響其自身的固有頻率,另一方面形狀變化后的電磁力有所增加導致的,這與仿真出的電主軸形變后各階頻率電磁力幅值變大是吻合的。由于電主軸氣隙偏心引起的電磁振動頻率為2~5倍頻。實驗測得1 000 Hz附近的最大振動幅值增加了13.47%,與仿真在1 000 Hz時電磁力幅值增長了12.1%的結果相比,利用該方法預測電主軸振動幅值,預測誤差為10.2%。

圖18 電主軸振動測試系統Fig.18 Vibration test system

圖19 電主軸形變前后振動頻譜分析對比圖Fig.19 Vibration spectrum contrast

6結論

采用機-電-熱-磁耦合法對150MD24Z7.5高速電主軸進行分析,得到電主軸內部電磁場、溫度場,結構場以及熱應力場,并對比分析了電主軸溫升形變后的振動變化,結論如下:

(1) 高速旋轉電主軸溫度隨運轉時間增加逐漸上升。其中轉軸溫度最高,定轉子齒頂形變量最大,且最大熱應力集中在定轉子齒槽處。

(2) 電主軸溫升形變導致氣隙變化,電主軸氣隙磁密變大,電主軸各階頻率電磁力幅值變大,由氣隙變化引起的電磁力幅值增加,振動幅值加大。

(3) 采用機-電-熱-磁耦合模型預測150MD24Z7.5高速電主軸振動特性,預測振動幅值誤差為10.2%。

參 考 文 獻

[1] Demetriades G D, De La Parra H Z, Andersson E,et al. A real-time thermal model of a permanent-magnet synchronous motor[J]. IEEE Transactions on Power Electronics,2010,25(2): 463-474.

[2] Holkup T, Cao H, Kola ′ r P, et al. Thermo-mechanical model of spindles[J]. CIRP Annals-Manufacturing Technology,2010, 59: 365-368.

[3] 芮執元,王志強,周亞敏. 熱變形對電主軸電機轉子與主軸過盈配合的影響[J].機械與電子,2011(7): 3-5.

RUI Zhi-yuan, WANG Zhi-qiang, ZHOU Ya-min. Effect of thermal deformation to the interference fit between the motor spindle’s rotor and spindle[J]. Machinery & Electronics,2011(7): 3-5.

[4] Huang Kun-fang, Juan Yu-ling, Tang Chia-hui, et al. A control scheme of thermal growth compensator for motorized spindles[J]. ICIC Express Letters,2011, 5(9A): 3101-3108.

[5] Creighton E, Honegger A, Tulsian A, et al. Analysis of thermal errors in a high-speed micro-milling spindle[J]. International Journal of Machine Tools & Manufacture,2010, 50: 386-393.

[6] 于翔. 主軸系統熱變形分析[J]. 應用能源技術,2008(10):47-50.

YU Xiang. Research and analysis of the thermal deformation of the spindle system[J]. Applied Energy Technology,2008(10):47-50.

[7] Abele E, Altintas Y, Brecher C. Machine tool spindle units[J]. CIRP Annals-Manufacturing Technology,2010, 59: 781-802.

[8] 何俊,賴玉活,羅錫榮,等. 基于 ANSYS Workbench 的數控車床主軸系統熱-結構耦合分析[J]. 組合機床與自動化加工技術, 2011(7): 19-22.

HE Jun, LAI Yu-huo, LUO Xi-rong, et al. Thermal-structural coupling analysis of CNC lathe spindle system based on ANSYS Workbench. [J]. Modular Machine Tool & Automatic Manufacturing Technique, 2011(7): 19-22.

[9] 謝黎明,邵寬平,靳嵐,等.高速電主軸的熱變形研究[J].機械制造,2011, 49(559): 41-43.

XIE Li-ming, SHAO Kuan-ping, QIN Lan, et al. Research on thermal deformation of motorized spindle[J]. Machinery,2011, 49(559): 41-43.

[10] 王明權,易傳云. 劃片機氣靜壓電主軸熱變形的有限元分析[J]. 電子工業專用設備,2007(147): 39-44.

WANG Ming-quan, YI Chuan-yun. FEA of thermal deformation for air spindle of dicing saw[J]. Equipment for Electronic Products Manufacturing,2007(147): 39-44.

[11] Zverev I A, Eun I U, Hwang Y, et al. An elastic deformation model of high speed spindle units[J]. International Journal of Precision Engineering and Manufacturing,2006(7): 39-46.

[12] Bossmanns B, Jay F T. Thermal model for high speed motorized spindles[J]. International Journal of Machine Tools and Manufacture, 1999,39: 1345-1366.

[13] Kim J D, Zverv I, Lee K B. Thermal model of high-speed spindle units[J]. Intelligent Information Management,2010, 2: 306-315.

[14] Uhlmanna E, Hua J. Thermal modeling of a high speed motor spindle[J]. CIRP Conference on High Performance Cutting,2012, 1: 313-318.

[15] Jedrzejewski J, Kowal Z, Kwasny W, et al. Hybrid model of high speed machining centre headstock[J]. CIRP Annals Manufacturing Technology, 2004, 53(1): 285-288.

[16] 陳觀慈,王黎欽,古樂,等.高速球軸承的生熱分析[J].航空動力學報,2007(1): 163-168.

CHEN Guan-ci, WANG Li-qin, GU Le, et al. Heating analysis of the high speed ball bearing[J]. Journal of Aerospace Power,2007(1): 163-168.

[17] Bradley K, Cao W, Clare J, et al. Predicting inverter induced harmonic loss by improved harmonic injection[J]. IEEE Transactions on Power Electronics, 2008, 23(5): 2619-2624.

[18] Liu R T, Mi C C, Gao D W Z, et al. Modeling of eddy-current loss of electrical machines and transformers operated by pulse-width-modulated inverters[J].IEEE Transactions on Magnetics, 2008, 44(8): 2021-2028.

[19] 王寶民,胡赤兵,孫建仁,等.高速電主軸熱態特性的ANSYS仿真分析[J].制造技術與機床,2009, 35(1): 28-31.

WANG Bao-min, HU Chi-bing, SUN Jian-ren, et al. Simulation analysis of thermal characteristics of high-speed motorized spindle by using ANSYS[J]. Journal of Lanzhou University of Technology,2009, 35(1): 28-31.

[20] 陳小安, 張朋, 合燁, 等.高速電主軸軸向振動研究[J].振動與沖擊,2014,33(20):70-74.

CHEN Xiao-an, ZHANG Peng,HE Ye, et al. Axial vibration of high-speed motorized spindles[J]. Journal of Vibration and Shock,2014,33(20):70-74.

[21] Chang Ching-Feng, Chen Jin-Jia. Vibration monitoring of motorized spindles using spectral analysis techniques[J]. Mechatronics, 2009, 19(9): 726-734.

基金項目:國家自然科學基金(51375317);遼寧省科技創新重大專項(201301001);遼寧自然科學基金項目:(2014020069)

收稿日期:2014-06-24修改稿收到日期:2014-12-26

通信作者閻銘 男,碩士,研究員,1988年生

中圖分類號:TM3;V414.3+3

文獻標志碼:A

DOI:10.13465/j.cnki.jvs.2016.01.011

Multi-field coupled model and dynamic performance prediction for 150MD24Z7.5 motorized spindle

ZHANG Li-xiu, YAN Ming, WU Yu-hou, LU Feng

(Shenyang Jianzhu University, Shenyang 110168, China)

Abstract:Motorized spindles are products of electromechanical integration. Their dynamic performance analysis and prediction are the preconditions for optimization design of machine tool spindle systems. Based on the coupling relationship among magnetic field,electric field,temperature field and structural field in a motorized spindle, the finite element model with multi-field coupled for 150MD24Z7.5 motorized spindle was established here. By calculating electromagnetic loss and bearing friction heat, the temperature field and structural field of the motorized spindle were simulated. The coupling relationship between thermal characteristics and vibration characteristics of the motorized spindle was discussed. The effects of air-gap variation after thermal deformation on the spindle’s vibration characteristics were analyzed and verified with tests. The results showed that thermal deformation has larger influences on the vibration amplitude, electromagnetic force amplitude increases about 12.1% due to the change of air-gap; with this model, the vibration amplitude of the motorized spindle can be predicted with an error of 10.2%.

Key words:motorized spindle; coupling; prediction; dynamic performance

第一作者 張麗秀 女,博士,副教授,1970年生

主站蜘蛛池模板: 亚洲无码不卡网| 亚洲国产成人精品一二区| 国产成人做受免费视频| 亚洲综合一区国产精品| 国产成人乱码一区二区三区在线| 丁香婷婷在线视频| 国产成人亚洲精品无码电影| 波多野结衣无码中文字幕在线观看一区二区 | 九色视频一区| 国产不卡网| 99精品热视频这里只有精品7| 欧美日韩国产在线观看一区二区三区 | 亚洲精品欧美重口| 国产无码高清视频不卡| 性欧美在线| 成人伊人色一区二区三区| 手机在线国产精品| 国产精品久久久久鬼色| 国产在线观看第二页| 精品91视频| 超清无码一区二区三区| 中文天堂在线视频| 亚洲成人播放| 成人国产一区二区三区| 亚洲一区二区黄色| 99re这里只有国产中文精品国产精品| 久久久久免费看成人影片 | 久久久久久久蜜桃| 亚洲AV成人一区二区三区AV| 国产在线观看人成激情视频| 亚洲国产日韩视频观看| www.精品视频| 国产成人AV综合久久| 热这里只有精品国产热门精品| 99久久精品国产精品亚洲| 人妻无码中文字幕第一区| 亚洲精品中文字幕午夜| 国产精品浪潮Av| 国产美女人喷水在线观看| 国产网站黄| 亚洲久悠悠色悠在线播放| 亚洲高清无在码在线无弹窗| 亚洲人成影视在线观看| 色呦呦手机在线精品| 97狠狠操| 午夜无码一区二区三区| 亚洲一区二区日韩欧美gif| 亚洲第一成网站| 亚洲成肉网| 精品无码国产一区二区三区AV| 亚洲成人网在线播放| 欧美成人精品欧美一级乱黄| 亚洲综合中文字幕国产精品欧美| 欧美日韩精品在线播放| 国产亚洲精| 99精品国产自在现线观看| 久草热视频在线| 日本免费一级视频| 激情网址在线观看| 国产免费a级片| 国产精品2| 国产日韩久久久久无码精品| 2020国产精品视频| 狠狠色噜噜狠狠狠狠色综合久| 国产日本欧美在线观看| 69视频国产| 国产69精品久久久久孕妇大杂乱 | 波多野结衣一区二区三区四区视频| 无码日韩人妻精品久久蜜桃| 蜜桃视频一区二区三区| 国产精品女主播| 97精品国产高清久久久久蜜芽| 欧美国产日产一区二区| 凹凸精品免费精品视频| 日本伊人色综合网| 久久久精品久久久久三级| 久久人人妻人人爽人人卡片av| 亚洲丝袜第一页| 91精品情国产情侣高潮对白蜜| 99热这里只有精品在线播放| 亚洲色精品国产一区二区三区| 激情国产精品一区|