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

一種盤式制動器鉗體輕量化設計研究

2015-10-29 06:27:03唐進元趙國偉
中國機械工程 2015年17期
關(guān)鍵詞:有限元優(yōu)化分析

唐進元 趙國偉

中南大學高性能復雜制造國家重點實驗室,長沙,410083

一種盤式制動器鉗體輕量化設計研究

唐進元趙國偉

中南大學高性能復雜制造國家重點實驗室,長沙,410083

以體積和疲勞壽命為輕量化設計目標,采用計算機輔助集成技術(shù)對盤式制動器的鉗體進行輕量化設計。利用有限元分析方法計算制動盤制動力矩并將其與實驗測得的制動力矩進行比較,驗證了有限元分析的步驟與方法正確可行;結(jié)合有限元分析得到的鉗體的最大主應力,根據(jù)名義應力法和Smith方法計算鉗體疲勞壽命;通過ISIGHT集成CATIA、ABAQUS和MATLAB對鉗體進行輕量化設計。

制動鉗體;輕量化;有限元分析;制動力矩;疲勞壽命

0 引言

汽車制動器主要分為盤式制動器和鼓式制動器,是汽車制動系統(tǒng)的主要零部件。盤式制動器由于其熱穩(wěn)定性好、水穩(wěn)定性好、制動力矩平穩(wěn)、制動塊上壓力分布均勻、噪聲低等優(yōu)點在汽車行業(yè)中得到了廣泛應用[1]。由于汽車制動性能的好壞直接關(guān)系到駕駛員和乘客的生命安全,因此制動器的設計是汽車設計中一項非常重要的工作。

目前,盤式制動器的研究主要涉及制動器的設計、加工以及制動性能的分析與提高等方面[2]。Hohmann等[3]建立了盤式制動器的有限元模型,對制動過程中接觸壓力的分布情況進行了研究。Liu等[4]利用ABAQUS通過復特征值法對盤式制動器制動過程中的尖叫現(xiàn)象進行了研究。Belhocine等[5]通過ANSYS對不同型號和不同材料的制動盤在制動過程中的熱現(xiàn)象進行了仿真分析。有關(guān)制動器的輕量化研究還處于起步階段,隨著汽車輕量化設計的發(fā)展,有必要對制動器進行輕量化設計。

現(xiàn)有的輕量化設計方法大都是改變零件的材料,使用質(zhì)量更小、機械性能更好的材料來代替原有材料[6];或先用CAE軟件進行有限元分析,然后將應力較小的區(qū)域去掉[7];或利用多目標遺傳算法,對機械零部件的近似模型進行優(yōu)化[8]。上述三種方法中,第一種方法改變材料后會提高機械零件的制造成本;第二種方法需要進行多次建模和分析,工作量大且難以確定最佳的優(yōu)化結(jié)果;第三種方法需要對模型進行大量的簡化,難以實現(xiàn)對復雜模型的輕量化。

本文不改變零件的材料,在保證疲勞壽命的情況下對盤式制動器的制動鉗體進行輕量化設計。首先通過有限元分析和實驗相結(jié)合的方法對盤式制動器制動力矩進行分析驗證;然后根據(jù)有限元計算得到的應力結(jié)果和材料的性能參數(shù)來計算鉗體的疲勞壽命;最后使用計算機輔助優(yōu)化軟件ISIGHT集成CATIA、ABAQUS和MATLAB對制動鉗體進行參數(shù)優(yōu)化。優(yōu)化過程中計算機自動更改優(yōu)化變量的值,并自動調(diào)用相關(guān)軟件進行分析計算,既節(jié)省時間又減少設計者的工作量,同時還能得到較好的輕量化結(jié)果。

1 盤式制動器有限元分析

1.1有限元模型

圖1 盤式制動器有限元模型的約束與加載情況

文本所用盤式制動器的有限元模型由制動盤、鉗體、制動塊和活塞組成[3]。被簡化的零部件對當前零部件的作用通過約束、載荷或邊界條件的形式施加到模型中。為了便于優(yōu)化時對優(yōu)化參數(shù)進行更改,盤式制動器相關(guān)零件的三維模型在CATIA中創(chuàng)建,通過STEP文件將其導入到ABAQUS中進行有限元分析。在ABAQUS中通過映射網(wǎng)格劃分方法將制動盤、制動塊和活塞劃分為六面體網(wǎng)格,為了便于優(yōu)化設計將鉗體劃分為四面體網(wǎng)格,劃分網(wǎng)格后的有限元模型如圖1所示。1.2約束與加載

當盤式制動器制動時,鉗體內(nèi)部液壓缸中的液體壓力推動活塞運動,使與活塞連接的制動塊與制動盤接觸,其反作用力推動鉗體運動并將壓力傳遞到另一側(cè)的制動塊使其與制動盤接觸,以此來實現(xiàn)對汽車的制動[3]。

在有限元模型中,我們將液體壓力施加給鉗體和活塞的相應位置(圖1b)。活塞與制動塊和鉗體與制動塊相接觸的面之間建立綁定約束;制動塊和制動盤之間建立接觸約束,設置摩擦因數(shù)為0.4。鉗體與導向銷連接的孔與其中心的參考點耦合,制動盤內(nèi)孔表面與其對應的參考點建立耦合約束。

盤式制動器的分析分為兩步進行(圖1c),第一步給鉗體和活塞施加大小為p0的壓力載荷,使制動塊與制動盤充分接觸并將其壓緊;第二步通過參考點RP-1使制動盤沿其軸線方向轉(zhuǎn)動。在整個過程中通過邊界條件使鉗體和制動塊只能沿制動盤的軸線方向移動,使制動盤只能沿其軸線方向轉(zhuǎn)動。

1.3結(jié)果分析

根據(jù)有限元分析結(jié)果,盤式制動器中的最大主應力出現(xiàn)在鉗體上,其大小和分布情況如圖2所示。

圖2 制動鉗體的最大主應力分布

在ABAQUS的后處理模塊中可以得到制動盤受到的制動力矩。有限元分析和實驗獲得的制動力矩的值如圖3所示。

圖3 有限元分析和實驗得到的制動力矩

對比圖3中的兩條曲線可以看出,有限元分析的過程與制動器實際制動過程基本一致,即在制動塊與制動盤產(chǎn)生相對運動之前制動力矩為0,開始制動后制動力矩在很短的時間內(nèi)達到最大值,并保持在該值附近波動直到制動結(jié)束。

由圖3可知,有限元分析得到的制動力矩的大小為2106 N·m,實驗測得的制動力矩大小為1883 N·m。實驗測得的制動力矩偏小,與有限元分析結(jié)果的誤差為10.6%。引起實驗測得的制動力矩偏小的原因主要是在有限元分析中沒有考慮液壓泵的機械效率,而是直接將液壓泵的壓力p0代入進行計算和分析。實際過程中,汽車盤式制動器液壓泵的機械效率一般取0.85~0.95[9],所以引起實驗測得的制動力矩比有限元分析結(jié)果小是正常的。結(jié)合上述有限元分析和實驗測得的制動力矩的對比結(jié)果,可以說明有限元分析得到的結(jié)果與真實情況接近,驗證了有限元分析結(jié)果的準確性。

圖3中有限元分析值相對于實驗值有一些滯后是因為該分析步的初始時間的取值偏大(這樣做可以提高分析速度)。減小初始時間可以減小兩者之間的差異,但分析速度過慢不利于之后的輕量化設計。

2 鉗體的疲勞壽命計算

2.1SN曲線的修正

根據(jù)SN曲線來判斷零件的疲勞壽命時,為了準確地計算出零件的疲勞壽命,首先需要對SN曲線進行修正。SN曲線通常都是在一定實驗條件下獲得的,當機械零件的實際工作條件和實驗條件不一致時需要對其SN曲線做相應的修正[10]。

SN曲線修正的關(guān)鍵是對疲勞極限進行修正。對疲勞壽命的修正主要根據(jù)不同的應力比、尺寸、表面精度和溫度等因素的影響進行修正。盤式制動器制動時的外載荷為液壓泵提供的液體壓力。本文所用載荷的施加方式與對盤式制動器進行疲勞測試時的加載方式相同。在疲勞測試時,車輪制動時液壓泵的壓力從0增大到p0,然后再卸載。因此,盤式制動器在實驗工況下的應力循環(huán)可以近似為脈動循環(huán),記脈動疲勞極限為σ0。

2.1.1疲勞極限的修正

鉗體的材料為QT450-10,材料參數(shù)如表1所示[11-12],其中σ-1為對稱循環(huán)應力下指定存活率為95%時的疲勞極限。根據(jù)相關(guān)手冊無法查閱該材料σ0的值,可以根據(jù)Smith圖來近似計算σ0[13]。

表1 QT450-10的材料參數(shù)

Smith圖是以平均應力σm為橫坐標,最大應力σmax和最小應力σmin為縱坐標的疲勞極限線圖,如圖4所示。圖4中,A點表示對稱疲勞極限σ-1,D點表示脈動疲勞極限σ0,C點表示強度極限σb,曲線ADC為最大應力線,曲線BEC為最小應力線。曲線ADC與曲線BEC所包圍的面積表示不產(chǎn)生破壞的應力水平。

圖4 Smith圖

由圖4可以看出,D點的橫坐標為BEC線與橫軸的交點(E點)。連接AC和BC(圖4中虛線所示),可得到BC線與橫軸交于點E′,通過E′作橫軸的垂線交AC于點D′,D′對應的縱坐標即為所求得的σ0。將表1中的數(shù)據(jù)代入圖4計算得σ0=300 MPa。上述方法計算得到的σ0的值比實際情況要小,這樣得到的SN曲線偏于保守。

機械零件的尺寸、表面精度和環(huán)境溫度等因素對疲勞極限的影響通過相應的影響系數(shù)來衡量[14]。本文考慮鉗體的尺寸和表面精度對疲勞極限的影響。記尺寸系數(shù)為Csize,表面加工系數(shù)為Csurf。根據(jù)相關(guān)手冊查得Csize和Csurf的值如表2所示。修正后的脈動疲勞極限為

(1)

相對于對稱疲勞極限σ-1的綜合修正系數(shù)為

(2)

表2 QT450-10材料的疲勞參數(shù)

2.1.2SN曲線的修正

圖5中,曲線Ⅰ為QT450標準試樣在對稱循環(huán)載荷下存活率為95%的SN曲線[11]。運用疲勞綜合修正系數(shù)對材料的SN曲線進行修正可以得到零件的SN曲線。由式(2)可得

σ′=σ/kσ

(3)

式中,σ′為作用于零件上的名義應力;σ為作用于標準試樣上的名義應力。

圖5 修正前后的SN曲線

在高周疲勞區(qū)域,零件的SN曲線可表示為

(σ′kσ)mN=σmN=C

(4)

式中,m為材料常數(shù)。

由式(4)可得修正后的SN曲線,即圖5中的曲線Ⅱ。

2.2疲勞壽命計算

當計算出零件的應力后,將其代入上述SN曲線中即可算得零件的疲勞壽命。通過SN曲線計算得到的是零件斷裂失效的疲勞壽命。

據(jù)統(tǒng)計,裂紋萌生壽命占零件總壽命的90%以上[15],因此也可以使用零件的裂紋萌生壽命來估算機械零件的疲勞壽命。Smith等提出了一種計算裂紋萌生疲勞壽命的方法[16](以下簡稱SWT公式):

(5)

材料疲勞參數(shù)的取值如表2所示。

上述SN曲線和材料的疲勞參數(shù)都是在單軸拉壓的條件下測得的,因此本文根據(jù)鉗體受到的最大主應力來計算疲勞壽命。將有限元分析得到的最大主應力代入SN曲線和SWT公式中計算得到的輕量化前的鉗體的疲勞壽命分別為2.32×106次和2.07×106次。在對制動器進行疲勞測試時,只要經(jīng)過2×105次循環(huán)后制動器沒有發(fā)生破壞即可視為合格。因此鉗體還有一定的輕量化空間。

3 輕量化設計

3.1輕量化流程

本文所用輕量化設計的思路是通過計算相關(guān)優(yōu)化變量在不同尺寸下的體積和疲勞壽命來確定最終的優(yōu)化尺寸。圖6為本文所用輕量化設計的流程圖。由圖6可以看出,在輕量化開始之前首先要確定輕量化對象并確定需要優(yōu)化的位置。

圖6 輕量化設計流程圖

本文的輕量化對象為鉗體,要優(yōu)化的位置為鉗體內(nèi)側(cè)凹槽(圖2)。確定好輕量化對象后在CATIA中建立用于有限元分析的模型,并通過宏錄制記錄修改優(yōu)化變量的過程,得到相關(guān)的宏命令;然后使用ABAQUS對模型進行有限元分析,輸出鉗體受到的最大主應力和體積,然后將有限元分析得到的最大主應力代入預先編寫的MATLAB程序中計算鉗體的疲勞壽命;最后在ISIGHT中調(diào)用CATIA、ABAQUS和MATLAB對鉗體進行輕量化分析;待優(yōu)化完成后再對優(yōu)化結(jié)果進行處理,若結(jié)果可行則可以根據(jù)輕量化結(jié)果試制產(chǎn)品,若結(jié)果不可行則尋找產(chǎn)生問題的原因,待問題解決后重新進行上述操作。

3.2ISIGHT集成CATIA、ABAQUS和MATLAB

圖7所示為在ISIGHT中集成CATIA、ABAQUS和MATLAB的優(yōu)化框架。使用Simcode組件通過批處理命令集成CATIA、ABAQUS、MATLAB和Delete模塊。其中CATIA模塊通過宏命令更改優(yōu)化變量的值,并為有限元分析提供模型。ABAQUS模塊通過python命令對更改后的模型進行有限元分析,并輸出鉗體的最大主應力和體積。MATLAB模塊用于計算并輸出鉗體的疲勞壽命。Delete模塊由自編程序組成,主要用于刪除當前循環(huán)中產(chǎn)生的中間文件,以保證下一個循環(huán)的順利執(zhí)行。Calculator組件用于將ABAQUS模塊中計算得到的最大主應力賦值給MATLAB中的相應變量。

圖7 ISIGHT集成CATIA、ABAQUS和MATLAB的優(yōu)化框架

3.3設置優(yōu)化參數(shù)

鉗體的優(yōu)化變量為內(nèi)側(cè)凹槽的深度,如圖2所示。選擇凹槽的半徑為優(yōu)化變量,記為R。在Loop組件中設置R由131 mm變化到135 mm,增量為0.02 mm。優(yōu)化目標為鉗體的體積和疲勞壽命,將體積小且疲勞壽命高的解作為最優(yōu)解。

3.4輕量化結(jié)果

優(yōu)化完成后可以在ISIGHT的History菜單中以列表的形式查看優(yōu)化結(jié)果,也可以在Graphs菜單下以圖表的形式查看。圖8~圖10所示分別為在ISIGHT中優(yōu)化得到的鉗體的體積、應力和疲勞壽命隨R變化的情況。

圖8 體積隨R的變化關(guān)系

圖9 最大主應力隨R的變化關(guān)系

圖10 疲勞壽命隨R的變化關(guān)系

由圖8可知,鉗體的體積隨著R的增加逐漸減小,由此說明通過改變R的值能對鉗體實現(xiàn)輕量化,且R的值越大輕量化效果越明顯。

圖9表示的是優(yōu)化過程中鉗體上的最大主應力隨R的變化情況。隨著R的變化,最大主應力的值在一定范圍內(nèi)波動,但仍有逐漸增大的趨勢,如圖中虛線所示。這表明隨著R的增加鉗體上的最大主應力會逐漸增大。

圖10所示為根據(jù)最大主應力計算得到的疲勞壽命隨R的變化關(guān)系。本文根據(jù)材料的SN曲線和SWT公式兩種方法計算得到鉗體的疲勞壽命。由圖10可見,兩種方法計算得到的疲勞壽命大小較為接近。總體而言,根據(jù)SWT公式計算得到的疲勞壽命比根據(jù)SN曲線計算得到的疲勞壽命更為保守。因此,本文根據(jù)SWT公式計算得到的疲勞壽命來近似作為零件的疲勞壽命。

結(jié)合圖8和圖10選取R的最終優(yōu)化尺寸為133.40 mm,優(yōu)化前后鉗體的體積和疲勞壽命如表3所示。由表3可以看出,通過優(yōu)化,鉗體的質(zhì)量減小了5%,此時鉗體仍具有較高的疲勞壽命。

表3 鉗體輕量化結(jié)果

4 結(jié)論

(1)通過有限元分析得到制動過程中制動盤受到的制動力矩,并將其與實驗測得的制動力矩進行比較,比較結(jié)果表明,有限元分析結(jié)果與真實情況一致性好。

(2)根據(jù)鉗體的尺寸和表面質(zhì)量對鉗體材料的疲勞極限和SN曲線進行修正,結(jié)合修正后的SN曲線和SWT公式對鉗體的疲勞壽命進行估算,根據(jù)優(yōu)化結(jié)果,選擇SWT公式計算結(jié)果作為鉗體的疲勞壽命。

(3)建立了ISIGHT集成CATIA、ABAQUS和MATLAB的優(yōu)化流程進行輕量化設計計算,設計計算結(jié)果表明:在保證鉗體的疲勞壽命的同時使其質(zhì)量減小了5%,在一定程度上實現(xiàn)了輕量化目標。

(4)本文僅對一個尺寸進行輕量化使鉗體的質(zhì)量減小了5%,可以對鉗體零件進行綜合分析確定出更多的優(yōu)化變量,進一步減小其質(zhì)量,本文方法同樣適用于其他機械零件,對進行輕量化設計有參考價值。

[1]蔡運迪, 唐文獻, 黃秋蕓, 等. 水冷盤式制動器熱疲勞失效有限元分析[J].中國機械工程, 2012, 23(22):2726-2731.

Cai Yundi, Tang Wenxian, Huang Qiuyun, et al. Finite Element Analysis(FEM) on Thermal Fatigue Failure of Water-cooling Disc Brake[J]. China Mechanical Engineering, 2012, 23(22):2726-2731

[2]K?semann C P, Huart M, Stobbe F, et al. Mechanical Braking System for the Pulsed Power Supply System of ASDEX Upgrade[J]. Fusion Engineering and Design, 2013, 88(9): 1491-1494.

[3]Hohmann C, Schiffner K, Oerter K. Contact Analysis for Drum Brakes and Disk Brakes Using ADINA[J]. Computers & Structures, 1999, 72(1): 185-198.

[4]Liu P, Zheng H, Cai C. Analysis of Disc Brake Squeal Using the Complex Eigenvalue Method[J]. Applied Acoustics, 2007, 68(6): 603-615.

[5]Belhocine A, Bouchetara M. Thermal Analysis of a Solid Brake Disc[J]. Applied Thermal Engineering, 2012, 32: 59-67.

[6]Duraiselvam M, Valarmathi A, Shariff S M. Laser Surface Nitrided Ti-6Al-4V for Light Weight Automobile Disk Brake Rotor Application[J]. Wear, 2014, 309(1): 269-274.

[7]Kim J G, Jang G W. Development of a Lightweight Frame for a 40-foot Flatbed Trailer by Using CAE-based Structural Optimization[J]. Proceedings of the Institution of Mechanical Engineers, Part D: Journal of Automobile Engineering, 2011, 225(5): 643-652.

[8]張勇, 李光耀, 王建華. 多目標遺傳算法在整車輕量化優(yōu)化設計中的應用研究[J]. 中國機械工程, 2009, 10(4): 500-503.

Zhang Yong, Li Guangyao, Wang Jianhua.Design Optimization on Lightweight of Full Vehicle Based on Muti -objective Genetic Algorithm[J]. China Mechanical Engineering, 2009, 10(4): 500-503.

[9]Jegadeeshwaran R, Sugumaran V. Comparative Study of Decision Tree Classifier and Best First Tree Classifier for Fault Diagnosis of Automobile Hydraulic Brake System Using Statistical Features[J]. Measurement, 2013, 46(9): 3247-3260.

[10]崔泗鵬, 姚衛(wèi)星, 夏天翔. 連接件振動疲勞壽命分析的名義應力法[J]. 中國機械工程, 2014, 25(18): 2519-2522.

Cui Sipeng, Yao Weixing, Xia Tianxiang. Nominal Stress Approach for Fatigue Life Prediction of Mutifastener Joints under Vibration Loading[J]. China Mechanical Engineering, 2014, 25(18): 2519-2522.

[11]朱森第. 機械工程材料性能數(shù)據(jù)手冊[M]. 北京: 機械工業(yè)出版社, 1994.

[12]秦大同, 謝里陽. 現(xiàn)代機械設計手冊(第6卷)[M]. 北京: 化學工業(yè)出版社, 2011.

[13]趙少汴, 王忠保. 抗疲勞設計——方法與數(shù)據(jù)[M]. 北京: 機械工業(yè)出版社, 1997.

[15]Fajdiga G, Sraml M. Fatigue Crack Initiation and Propagation under Cyclic Contact Loading[J]. Engineering Fracture Mechanics, 2009, 76(9): 1320-1335.[16]Xia Z, Kujawski D, Ellyin F. Effect of Mean Stress and Ratcheting Strain on Fatigue Life of Steel[J]. International Journal of Fatigue, 1996, 18(5): 335-341.

(編輯郭偉)

Research on Lightweight Design for Disc Brake Caliper

Tang JinyuanZhao Guowei

State Key Labratory of High Performance Complex Manufacturing,Changsha,410083

Taking volume and fatigue life as lightweight design goal, an integrated computer-aided technology was used for disc brake caliper’s lightweight design. First, through comparing the barking torque which got by FEA and experiments it is proved that FEA’s steps and methods are correct and feasible. Then, combining the maximum principal stress which got by FEA with nominal stress method and Smith method was used to calculate the fatigue of caliper. Finally, Through ISIGHT integrating CATIA、ABAQUS and MATLAB, the lightweight of caliper was performed.

brake caliper; lightweight; finite element analysis(FEA); brake torque; fatigue life

2014-12-02

國家自然科學基金資助項目(51275530)

TP302DOI:10.3969/j.issn.1004-132X.2015.17.002

唐進元,男,1962年生。中南大學機電工程學院教授。研究方向為機械設計及理論。趙國偉,男,1989年生。中南大學機電工程學院碩士研究生。

猜你喜歡
有限元優(yōu)化分析
超限高層建筑結(jié)構(gòu)設計與優(yōu)化思考
民用建筑防煙排煙設計優(yōu)化探討
關(guān)于優(yōu)化消防安全告知承諾的一些思考
一道優(yōu)化題的幾何解法
隱蔽失效適航要求符合性驗證分析
電力系統(tǒng)不平衡分析
電子制作(2018年18期)2018-11-14 01:48:24
電力系統(tǒng)及其自動化發(fā)展趨勢分析
磨削淬硬殘余應力的有限元分析
基于SolidWorks的吸嘴支撐臂有限元分析
箱形孔軋制的有限元模擬
上海金屬(2013年4期)2013-12-20 07:57:18
主站蜘蛛池模板: 欧美性爱精品一区二区三区| 狠狠色香婷婷久久亚洲精品| 97人人做人人爽香蕉精品| 啪啪永久免费av| 国产福利一区在线| 亚洲午夜福利精品无码不卡| 永久免费av网站可以直接看的 | 亚洲欧洲自拍拍偷午夜色无码| 亚洲浓毛av| 日本人又色又爽的视频| 国内精品视频| 国产精品嫩草影院视频| 日韩福利在线视频| 亚洲第一页在线观看| 久久国产高潮流白浆免费观看| 婷婷久久综合九色综合88| 国内嫩模私拍精品视频| 久久久久青草线综合超碰| 91精品免费久久久| 天堂av综合网| 美女毛片在线| 日韩在线2020专区| 国产美女视频黄a视频全免费网站| 国产中文一区二区苍井空| 久久青草热| 久久久波多野结衣av一区二区| 1024国产在线| 91成人免费观看| 欧美一级夜夜爽www| 超碰aⅴ人人做人人爽欧美 | 久爱午夜精品免费视频| 精品国产www| 亚洲欧美不卡中文字幕| 一本久道久综合久久鬼色| 欧美日本一区二区三区免费| 国产成人精品在线| 亚洲精品色AV无码看| 亚洲天天更新| 国产精品大尺度尺度视频| 国产亚洲欧美日韩在线一区二区三区| 精久久久久无码区中文字幕| 亚洲精品无码日韩国产不卡| 91在线无码精品秘九色APP| 男女精品视频| 99热这里只有免费国产精品 | 黄色网站不卡无码| 免费激情网站| 亚洲精品大秀视频| 国产尤物视频在线| 99热国产这里只有精品无卡顿"| 91成人在线免费视频| 欧美有码在线观看| 成人91在线| 为你提供最新久久精品久久综合| 亚洲激情区| 久久成人国产精品免费软件 | 国产无套粉嫩白浆| a天堂视频| 欧美翘臀一区二区三区| 亚洲男人的天堂久久香蕉| 亚洲最黄视频| 东京热一区二区三区无码视频| 午夜在线不卡| 国产精品嫩草影院视频| 精品国产欧美精品v| 噜噜噜综合亚洲| 91最新精品视频发布页| 在线欧美日韩国产| 亚洲乱亚洲乱妇24p| 色九九视频| 91欧美亚洲国产五月天| 伊人中文网| 国产高清不卡视频| 四虎免费视频网站| 真人免费一级毛片一区二区| 最新亚洲人成无码网站欣赏网 | 国产无吗一区二区三区在线欢| 在线视频亚洲欧美| 国产日韩欧美视频| 精品久久人人爽人人玩人人妻| 九色国产在线| 国产男女免费视频|