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

散熱筋結(jié)構(gòu)盤式制動器溫度場有限元分析與優(yōu)化

2017-08-09 02:05:04潘公宇闞云峰
關(guān)鍵詞:有限元優(yōu)化分析

潘公宇,闞云峰

(江蘇大學(xué) 汽車與交通工程學(xué)院, 江蘇 鎮(zhèn)江 212013)

?

散熱筋結(jié)構(gòu)盤式制動器溫度場有限元分析與優(yōu)化

潘公宇,闞云峰

(江蘇大學(xué) 汽車與交通工程學(xué)院, 江蘇 鎮(zhèn)江 212013)

結(jié)合熱分析理論,采用數(shù)值法代入溫度邊界條件,對帶有散熱筋結(jié)構(gòu)的通風(fēng)盤式制動器進行了溫度場分析和優(yōu)化。建立制動盤的三維模型,在此基礎(chǔ)上對制動盤的溫度特性進行瞬態(tài)溫度場分布的仿真分析,得到制動盤在制動行程中的溫度場、內(nèi)外部溫度梯度分布以及溫度變化的時間曲線。使用Isight進一步對制動盤進行了結(jié)構(gòu)優(yōu)化,改善溫度場分布并返回到有限元模型進行仿真對比論證,為制動器的設(shè)計優(yōu)化提供了相應(yīng)的理論基礎(chǔ),也在很大程度上為由溫度場引起的制動器失效判斷提供了分析改進手段。

有限元法;溫度場;瞬態(tài)分析;結(jié)構(gòu)優(yōu)化

隨著交通運輸?shù)陌l(fā)展,越來越多的制動器廠商對制動器的設(shè)計制造提出了更加嚴(yán)格的質(zhì)量要求,制動器的規(guī)格尺寸逐步降低,而相應(yīng)的制動器部件所承載的傳遞熱量大大增加。

在盤式制動器制動過程中,摩擦片夾緊制動盤產(chǎn)生制動力矩,將車輛行駛動能轉(zhuǎn)化為摩擦產(chǎn)生的熱能消散。在制動過程中,制動盤因為受到摩擦而產(chǎn)生溫度變化,由此而產(chǎn)生熱衰退現(xiàn)象以及摩擦片偏磨現(xiàn)象[1],以及由溫度場分布中熱應(yīng)力過高引起的制動盤裂紋、翹曲等,嚴(yán)重影響了行車安全。因此,在對制動盤制動過程中溫度場分析的基礎(chǔ)上,針對制動過程中制動盤溫度場的改善以及優(yōu)化措施研究具有十分重要的意義。

1 能量傳遞分析

對制動過程中制動熱量的傳遞過程進行分析,以熱傳遞方式的形式計算能量的傳遞,并以能量折算的方式將其作為邊界條件加載在制動盤上以模擬制動過程中的溫度變化。

1.1 制動熱量及其傳輸

采用熱力學(xué)第一定律即能量守恒定律進行熱力學(xué)分析。封閉系統(tǒng)沒有質(zhì)量的流入和流出,熱能和機械能在轉(zhuǎn)移或轉(zhuǎn)換時,能量的總量必定守恒:

(1)

式中:Q為熱量;W為做功; ΔU為系統(tǒng)的內(nèi)能; ΔKE為系統(tǒng)的動能; ΔPE為系統(tǒng)的勢能。

對于瞬態(tài)熱分析,熱傳遞速率q=du/dt,即流入或流出的熱傳遞速率q等于系統(tǒng)內(nèi)能單位時間內(nèi)的變化[2]。

在汽車制動過程中,汽車將制動器運動部件的動能轉(zhuǎn)化成內(nèi)能[3],即

(2)

式中:M為整車質(zhì)量(kg);v0為車輛制動初速度(m/s);v1為制動行程末車輛速度(m/s);g為重力加速度,9.8 m/s2;i為道路坡度系數(shù),下坡路段大于0,上坡路段小于0;f為滾動阻力系數(shù);s為車輛制動距離(m)。

在車輛制動過程中,熱傳遞方式包括熱傳導(dǎo)、熱對流和熱輻射等3種,而其中熱輻射在實際分析過程中權(quán)重較小,可以忽略[4]。

1.2 能量折算法

采用能量折算法計算摩擦片與制動盤摩擦生熱產(chǎn)生的熱流密度,同時假設(shè)摩擦產(chǎn)生的熱量在摩擦表面上均勻分布,根據(jù)能量轉(zhuǎn)化以及能量守恒定律,從能量折算的角度分析車輛的制動行程,將制動過程中的動能全部轉(zhuǎn)化為制動盤摩擦產(chǎn)生的熱能。

熱流密度與對流換熱系數(shù)分別以一組與時間t相關(guān)的函數(shù) q表示,實際分析中將t轉(zhuǎn)化為更直觀的車輛行駛速度v代入計算。

2 有限元模型及其參數(shù)設(shè)置

從熱量傳遞上分析,制動盤吸收了在制動期間制動器產(chǎn)生的總熱量的 95%,而摩擦片或制動襯塊僅吸收剩余熱量的5%[3,5]。

在制動過程中,制動盤受到摩擦片的法向和切向作用力,同時受到熱載荷的作用。吳佳偉等[5]以提高制動盤散熱性能為目的,在制動盤盤體中間鑄有供氣流通過的散熱筋結(jié)構(gòu),即為雙層盤結(jié)構(gòu)的通風(fēng)盤式制動器,增加了散熱面積,有效地降低了溫升。

本文設(shè)計具有螺旋散熱風(fēng)道的通風(fēng)盤式制動器進行瞬態(tài)溫度場分析。螺旋式風(fēng)道制動盤縱向總厚度為 36mm,按順時針螺旋狀分布34條散熱筋,散熱筋厚度為 10mm,單個散熱筋寬度為9.8mm,兩側(cè)圓盤厚度均為 13mm。制動盤內(nèi)圓孔直徑為 175mm,外圓盤直徑為315mm。

使用三維軟件對實車制動盤建模。建模過程中做出以下簡化措施:

1) 制動盤所受壓力處均勻分布。

2) 制動過程中材料屬性不發(fā)生變化。

3) 忽略制動盤小特征影響。

4) 在制動過程中動能完全轉(zhuǎn)化為熱能。

5) 制動盤表面摩擦因數(shù)恒定。

制動盤模型如圖1所示。制動盤材料特性如表1所示。

圖1 制動盤三維模型表1 制動盤材料特性

分析對象材料熱傳導(dǎo)系數(shù)/(W·(m2·K)-1)比熱容/(J·(kg·K)-1)密度/(t·mm-3)制動盤鑄鐵405107.1×10-9

3 模型熱傳遞參數(shù)的構(gòu)建

整個模型的熱傳遞參數(shù)由熱流密度、對流換熱系數(shù)以及環(huán)境初始溫度確定,以時間為參考系加載至Ansys熱分析單元中,載荷步設(shè)置為0.2s。

3.1 熱流密度

制動盤溫度場分析中的熱載荷為熱流密度,以數(shù)值方式加載于制動盤面與制動閘片的摩擦表面上。假定整個制動過程為勻減速過程,可得熱流密度q與時間t的關(guān)系:

(3)

式中:q(t)為t時刻加載于制動盤表面的熱流密度(kW/m2);k為軸質(zhì)量(kg);a為制動加速度(m/s2);v0為制動初速度(m/s);n為每根軸上裝配的制動盤個數(shù);R和r分別為制動盤閘片與盤面摩擦的環(huán)形區(qū)域的外徑和內(nèi)徑(m)。丁群等[6]根據(jù)以往文獻的經(jīng)驗,認(rèn)為車體動能轉(zhuǎn)化為熱能的效率Z為0.9。

3.2 對流換熱系數(shù)

制動盤的對流換熱是指制動過程中制動器與周圍接觸的空氣之間由于溫差的存在而引起的熱量交換。制動盤的對流換熱系數(shù)近似為[7]

(4)

式中:v為車速(m/s);α是經(jīng)驗公式系數(shù),前輪制動取0.7,后輪制動取0.3。

3.3 環(huán)境初始溫度及熱輻射

環(huán)境初始溫度設(shè)置為22 ℃,而熱輻射在熱量交換中所占權(quán)重較低,所以忽略不計。

4 制動盤溫度場仿真分析

使用Hypermesh對模型進行網(wǎng)格劃分,制動盤有限元網(wǎng)格模型如圖2所示。

圖2 制動盤有限元網(wǎng)格劃分

以某型號汽車在高速公路上制動為初始工況。制動初速度為120km/h,制動時間為4.2s,之后汽車完全停止。在此過程中制動的加速度為恒值。

根據(jù)整車相關(guān)數(shù)據(jù)得到的熱流密度q(t)與對流換熱系數(shù)hr(t)如表2所示。

設(shè)置初始溫度為22 ℃,通過函數(shù)加載的方式將熱流密度以及對流換熱系數(shù)加載到制動盤相應(yīng)區(qū)域,通過仿真計算,得到的制動盤溫度場云圖如圖3所示。相對應(yīng)的溫度變化曲線如圖4所示。

表2 熱流密度與對流換熱系數(shù)

圖3 制動盤溫度場云圖

圖4 制動盤溫度變化曲線

從溫度云圖可以看出:在制動盤與摩擦片接觸區(qū)域溫度上升幅度最大,沿著制動盤徑向向中心輻射溫度梯度逐漸遞減。從最高溫度變化曲線可以看出:整個制動過程的溫度梯度變化呈現(xiàn)快速上升和緩慢下降2個階段,在制動過程中溫度達到最高值,并且在之后的制動過程中溫度逐漸降低。

沿制動盤徑向截取剖面得到的溫度分布以及溫度變化曲線如圖5、圖6所示。

圖5 制動盤徑向溫度分布

圖6 制動盤徑向溫度變化曲線

通過截取的制動盤剖面溫度云圖可以得到:制動開始后,溫度由制動盤摩擦表面逐漸向制動盤內(nèi)部傳導(dǎo);在厚度方向,溫度在從制動盤摩擦表面到散熱筋板的方向上逐漸減小。

另外可以看出:沿著制動盤徑向向外方向溫度上升速度呈現(xiàn)逐漸加快的趨勢,在制動盤邊緣溫度上升速度最高。

處于邊緣處的線速度由于半徑較大,大于靠近制動盤中心方向的其他位置,所以在單位時間內(nèi)由于摩擦產(chǎn)生的熱量也隨之變大,導(dǎo)致制動盤邊緣溫度上升要大于內(nèi)側(cè)。這一部分區(qū)域?qū)?yīng)的摩擦片外側(cè)區(qū)域接觸溫度也要高于摩擦片內(nèi)側(cè),在一定程度上加劇了摩擦片在內(nèi)外側(cè)方向上的偏磨現(xiàn)象[1]。

截取制動盤與散熱筋片接觸區(qū)域得到的溫度分布以及溫度變化曲線如圖7、圖8所示。

圖7 散熱筋溫度分布

圖8 散熱筋溫度變化曲線

由圖7、圖8可以看出:散熱筋片根據(jù)接觸制動盤內(nèi)表面區(qū)域的不同,溫度分布也呈現(xiàn)梯度分布,但與制動盤表面溫度變化趨勢不同的是,散熱筋片的最高溫度變化并沒有呈現(xiàn)快速上升和緩慢下降2個階段,而是隨著制動過程一直存在溫度上升的現(xiàn)象,也就是在整個制動過程中一直處于吸收制動盤外表面向內(nèi)部傳遞的熱量的狀態(tài)。散熱筋片的設(shè)計有效地改善了制動盤在制動過程中的熱傳遞問題,這也從一個側(cè)面反映出散熱風(fēng)道的設(shè)計對整個通風(fēng)盤式制動器散熱性能的改善有著不可忽視的影響。

5 通風(fēng)盤式制動器制動盤優(yōu)化設(shè)計

5.1 優(yōu)化方式及設(shè)計變量

通過改變制動盤厚度(x1)、散熱筋厚度(x2)、散熱筋縱向?qū)挾?x3)來改變制動盤散熱性能(這里主要指制動盤最高瞬態(tài)溫度)。

相較于有限元法的優(yōu)化過程,將有限元模型擬合成無限近似的數(shù)學(xué)模型能夠準(zhǔn)確、快捷地求解最優(yōu)解。求解優(yōu)化近似模型的過程包括樣本數(shù)據(jù)的采集、樣本模型擬合近似數(shù)學(xué)模型、近似模型的誤差分析以及對近似模型的精度優(yōu)化等。

5.2 優(yōu)化近似模型

為了最大程度地擬合近似數(shù)學(xué)模型,需要使用試驗設(shè)計的方法求解30組仿真解,而仿真解的采樣使用最優(yōu)拉丁超立方設(shè)計法,它比普通的正交采樣覆蓋率和均衡性都更高,具有非常優(yōu)秀的精確度,有效地覆蓋了采樣空間,能夠逼近至少2階的非線性關(guān)系[8](圖9)。

圖9 近似模型采樣方法

試驗設(shè)計所選擇的樣本設(shè)計變量以及變量的變化范圍圍繞優(yōu)化設(shè)計方式以及設(shè)計變量來確定,即在制動盤厚度、散熱筋厚度以及散熱筋縱向?qū)挾冗@3個設(shè)計變量上分別選取合適的范圍進行樣本采樣。具體采樣范圍如表3所示。

表3 設(shè)計變量以及變量的變化范圍

本文使用Isight中提供的最優(yōu)拉丁超立方試驗設(shè)計的采樣方法進行采樣,將選取的樣本點代入到有限元模型中進行仿真計算,得出最高溫度y。

樣本點的采樣以及計算結(jié)果如表4所示。

表4 樣本采樣數(shù)據(jù)

在30個采樣數(shù)據(jù)完成后,利用Isight的二次響應(yīng)面法來構(gòu)建制動盤的優(yōu)化數(shù)學(xué)模型,就是將這30組數(shù)據(jù)擬合成與變量x1、x2、x3有關(guān)的數(shù)學(xué)模型[9],最高溫度y= f (x1, x2, x3),即:

y=-24.685x1-10.651 4x2+17.408 15x3-0.101x1x2+0.443 061x1x3-1.528 31x2x3+0.794 672x12+1.316 931x22-0.474 3x32+323.327 8

(5)

因為建立的近似模型并不是響應(yīng)變量和設(shè)計變量之間的真正模型,而是近似的關(guān)系式,所以存在一定的誤差,本研究使用R2來進行誤差分析的評估。如果誤差較小則可以進一步使用模型進行優(yōu)化分析;如誤差較大,則繼續(xù)增加樣本點進行模型的重新構(gòu)建。

誤差分析采用5個樣本點(表5)以及計算值進行分析,對應(yīng)的R2結(jié)果和誤差分析如圖10所示。

表5 誤差分析樣本點

圖10 近似模型誤差分析

這樣,優(yōu)化的數(shù)學(xué)模型就可以用如下公式代替原來的復(fù)雜有限元模型計算:

最小厚度

最小厚度

最小間距

y=制動盤最高溫度擬合模型

即:

0.1

0.06

0.08

y=-24.685x1-10.651 4x2+ 17.408 15x3-0.101x1x2+ 0.443 061x1x3-1.528 31x2x3+ 0.794 672x12+1.316 931x22- 0.474 3x32+323.327 8

(6)

5.3 制動盤優(yōu)化過程

賴宇陽[10]利用Isight的Point算法來進行迭代運算,求解最優(yōu)解。Point 算法由4種算法組成,包括遺傳算法、最速下降法、線性單純形法以及序列二次規(guī)劃法。具體流程在Isight中的顯示如圖11所示。

在經(jīng)過一定次數(shù)的迭代運算之后,包括可行解、最優(yōu)解等在內(nèi)的可能結(jié)果會顯示在程序中。迭代過程中各項變量的可行解。

圖11 迭代算法流程框圖

圖12 近似模型參數(shù)迭代過程

最終的優(yōu)化結(jié)果如圖13所示。

Optimumdesignpoint:Run#=19736Objective=168.00976415641023Penalty=0.0ObjectiveAndPenalty=168.00976415641023x1=12.917281150817871x2=11.502419471740723x3=12.0y=168.00976415641023

圖13 迭代求解的優(yōu)化結(jié)果

5.4 有限元分析驗證

由最終優(yōu)化結(jié)果的各項變量值可以得出優(yōu)化后最高溫度168.1 ℃。在求出最優(yōu)解后將其轉(zhuǎn)化為有限元模型再代入有限元軟件中進行仿真分析,比較優(yōu)化前后溫度場分析結(jié)果的變化,評估優(yōu)化結(jié)果。優(yōu)化前后溫度場分析結(jié)果如圖14所示。

由圖14可知:溫度場在優(yōu)化后最高溫度下降了8 ℃左右,達到了169.97 ℃,說明優(yōu)化變量的選擇是有效的。另外溫度場仿真的結(jié)果與優(yōu)化近似模型的計算結(jié)果非常接近,進一步證明了優(yōu)化近似模型的準(zhǔn)確性。

圖14 優(yōu)化前后溫度場分析結(jié)果對比

6 結(jié)論

基于帶有散熱筋結(jié)構(gòu)的制動器三維模型,通過構(gòu)建有限元模型并導(dǎo)入進行溫度場仿真分析,研究在緊急制動的情況下通風(fēng)盤式制動器的溫度場分布,并根據(jù)優(yōu)化擬合模型研究優(yōu)化措施,得出以下結(jié)論:

1) 在緊急制動工況下,高溫區(qū)域主要集中在制動盤摩擦區(qū)域表面。制動摩擦產(chǎn)生的熱量在縱向沿著制動盤厚度方向向制動盤內(nèi)部傳遞,并且制動盤邊緣也就是摩擦片外側(cè)溫度上升較快,一定程度上加劇了摩擦片內(nèi)外側(cè)偏磨現(xiàn)象的產(chǎn)生。

2) 散熱筋構(gòu)造對制動盤散熱性能有著一定的影響,其設(shè)計優(yōu)劣勢必影響制動盤制動性能的穩(wěn)定。在某些情況下,良好的散熱結(jié)構(gòu)設(shè)計有助于改善制動盤摩擦表面溫度場分布的不均勻以及接觸溫度過高的情況,而摩擦表面的接觸溫度以及溫度梯度集中反映了載荷、速度、摩擦因數(shù)、材料的熱物理特性以及耐久性等因素的影響。

3) 對制動盤溫度場分布提出了更為準(zhǔn)確、快捷的結(jié)構(gòu)優(yōu)化方案,為相應(yīng)的熱應(yīng)力集中、材料的熱疲勞耐久問題研究提供了參考。

[1] 張光榮,謝敏松,黎軍,等.摩擦片偏磨引起的汽車制動低鳴噪聲[J].機械工程學(xué)報,2013,49(9):81-86.

[2] 臧國群,何忠韜,丁立利.基于ANSYS的高速客車車頂穩(wěn)態(tài)傳熱溫度場分析[J].機械工程與自動化,2013(1):79-81.

[3] ASHWORTH R J,NEWCOMB T P.Temperature distributions and thermal distortions of brake drums[J].Proc I Mech E,1977(191):169-177.

[4] 趙華,蔣卓良.空間薄壁圓管瞬態(tài)溫度場有限元分析[J].重慶理工大學(xué)學(xué)報(自然科學(xué)),2010,24(4):44-48.

[5] 吳佳偉,楊志剛.氣流方向?qū)νL(fēng)制動盤散熱性能的影響[J].汽車工程學(xué)報,2014,4(6):418-423.

[6] 丁群,謝基龍.基于三維模型的制動盤溫度場和應(yīng)力場計算[J].鐵道學(xué)報,2002,24(6):34-38.

[7] 馬迅,朱前進.蹄鼓式制動器瞬態(tài)溫度場的仿真分析[J].機械設(shè)計與制造,2008(6):71-73.

[8] 何為,薛衛(wèi)東,唐斌.優(yōu)化試驗設(shè)計方法及數(shù)據(jù)分析[M].北京:化學(xué)工業(yè)出版社,2012.

[9] 閔亞能.試驗設(shè)計(DOE)應(yīng)用指南[M].北京:機械工業(yè)出版社,2011.

[10]賴宇陽.Isight參數(shù)優(yōu)化理論與實例詳解[M].北京:北京航空航天大學(xué)出版社,2012.

(責(zé)任編輯 林 芳)

Temperature Finite Element Analysis and Optimization of the Radiating Rib Structure Disc Brake

PAN Gongyu, KAN Yunfeng

(School of Vehicle and Traffic Engineering, Jiangsu University, Zhenjiang 212013, China)

Combined with thermal analysis theory, numerical methods are used to substitute temperature boundary conditions, and the ventilated disc brakes with cooling ribs structure temperature field were analyzed and optimized. 3D model of the brake disc is built. Based on the temperature characteristics of the brake disc, we have the simulation analysis of transient temperature field distribution analysis obtained in the brake disc brake stroke in temperature, external temperature and temperature gradient distribution time variation curve. And using Isight, the brake discs is further structural optimized to improve the temperature field and return the results to finite element model simulation comparison argument. It provides a brake corresponding theoretical basis for design optimization, but also it provides the improved analysis tools to the brake failure judgment causedby the temperature fieldto a large extent.

finite element method; temperature field; transient analysis; structure optimization

2016 -08-09

國家自然科學(xué)基金資助項目(51375212),江蘇省道路載運工具應(yīng)用新技術(shù)重點實驗室開放基金資助項目(201509)

潘公宇(1965—),男,博士,教授,主要從事車輛系統(tǒng)動力學(xué)、車輛動態(tài)設(shè)計理論、車輛振動控制技術(shù)等方面的研究,E-mail:pangongyu@hotmail.com。

潘公宇,闞云峰.散熱筋結(jié)構(gòu)盤式制動器溫度場有限元分析與優(yōu)化[J].重慶理工大學(xué)學(xué)報(自然科學(xué)),2017(7):12-19.

format:PAN Gongyu,KAN Yunfeng.Temperature Finite Element Analysis and Optimization of the Radiating Rib Structure Disc Brake[J].Journal of Chongqing University of Technology(Natural Science),2017(7):12-19.

10.3969/j.issn.1674-8425(z).2017.07.002

U270.35

A

1674-8425(2017)07-0012-08

猜你喜歡
有限元優(yōu)化分析
超限高層建筑結(jié)構(gòu)設(shè)計與優(yōu)化思考
民用建筑防煙排煙設(shè)計優(yōu)化探討
關(guān)于優(yōu)化消防安全告知承諾的一些思考
一道優(yōu)化題的幾何解法
隱蔽失效適航要求符合性驗證分析
電力系統(tǒng)不平衡分析
電子制作(2018年18期)2018-11-14 01:48:24
電力系統(tǒng)及其自動化發(fā)展趨勢分析
磨削淬硬殘余應(yīng)力的有限元分析
基于SolidWorks的吸嘴支撐臂有限元分析
箱形孔軋制的有限元模擬
上海金屬(2013年4期)2013-12-20 07:57:18
主站蜘蛛池模板: 99在线视频精品| 亚洲第一区精品日韩在线播放| 国产肉感大码AV无码| 国产欧美日韩18| 午夜天堂视频| 99精品国产自在现线观看| 一级爆乳无码av| 9966国产精品视频| 91福利片| 国产精品一区二区不卡的视频| 欧美另类视频一区二区三区| 久久久久国产一区二区| 国产高清免费午夜在线视频| 国产午夜无码片在线观看网站| 欧美性猛交一区二区三区| 亚洲色图在线观看| 五月婷婷导航| 丝袜亚洲综合| 黄色网在线| 国产va在线观看免费| 亚瑟天堂久久一区二区影院| 亚洲av无码成人专区| 国产成人一区在线播放| 在线观看免费黄色网址| 香蕉色综合| 亚洲成a人在线观看| 色欲色欲久久综合网| 国产99精品视频| 日本欧美在线观看| 中文字幕无码av专区久久| 国产成人无码播放| 亚洲资源站av无码网址| 亚洲午夜天堂| m男亚洲一区中文字幕| 亚洲欧美h| 亚洲成aⅴ人在线观看| 97久久精品人人做人人爽| 国产成人AV综合久久| 国产特一级毛片| 国产成人av一区二区三区| 伊人色天堂| 丁香五月激情图片| 99青青青精品视频在线| 亚洲成年网站在线观看| 日韩美毛片| 国产成人一区二区| 91毛片网| 色婷婷国产精品视频| 国产第八页| 亚洲乱码视频| 欧美黄网在线| 亚洲第一福利视频导航| 国产精品吹潮在线观看中文| 久操中文在线| 欧美色综合网站| 亚洲综合在线网| 在线va视频| 国产日韩丝袜一二三区| 精品一区国产精品| 午夜日韩久久影院| 男人天堂伊人网| 高清免费毛片| 午夜国产理论| 亚洲欧美国产五月天综合| 国产剧情伊人| 午夜欧美在线| 亚洲天堂日本| 午夜成人在线视频| 米奇精品一区二区三区| 久久综合国产乱子免费| 好紧太爽了视频免费无码| 欧美日本在线播放| 欧美精品在线视频观看| 91丝袜乱伦| 午夜一区二区三区| 92精品国产自产在线观看| 国产性精品| 1级黄色毛片| 在线a网站| 欧美一区精品| 国产精欧美一区二区三区| AV在线麻免费观看网站|