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

電動(dòng)汽車(chē)牽引用水冷異步電機(jī)耦合場(chǎng)分析*

2018-05-15 00:49:29王小飛
電機(jī)與控制應(yīng)用 2018年4期

王小飛, 代 穎, 羅 建, 黃 闖

(上海大學(xué) 機(jī)電工程與自動(dòng)化學(xué)院,上海 200072)

0 引 言

異步電機(jī)是電動(dòng)汽車(chē)驅(qū)動(dòng)電機(jī)的主要類(lèi)型之一。其成本低、結(jié)構(gòu)堅(jiān)固、可靠性高、調(diào)速范圍寬,特斯拉電動(dòng)車(chē)Model S、寶馬Mini-E及部分新能源客車(chē)采用了異步電機(jī)作為牽引電機(jī)[1]。電動(dòng)汽車(chē)牽引用電機(jī)電磁負(fù)荷大、工作環(huán)境惡劣,對(duì)電機(jī)的散熱設(shè)計(jì)提出了挑戰(zhàn),過(guò)高的溫升不但會(huì)破壞電機(jī)的定子繞組絕緣、降低籠型繞組的強(qiáng)度,還會(huì)造成導(dǎo)條斷裂[2-9],還可能導(dǎo)致鐵心與機(jī)殼的過(guò)盈配合由于熱膨脹系數(shù)不同引起的熱應(yīng)力而出現(xiàn)滑移或漲開(kāi)現(xiàn)象,降低電機(jī)結(jié)構(gòu)剛度,引發(fā)較大的電磁振動(dòng)。

本文以1臺(tái)額定功率40 kW的車(chē)用異步電動(dòng)機(jī)為研究對(duì)象,基于Ansys Workbench多物理場(chǎng)仿真平臺(tái),對(duì)電機(jī)運(yùn)行在額定工況(40 kW)和峰值工況(80 kW)進(jìn)行磁-熱-結(jié)構(gòu)的多物理場(chǎng)耦合的有限元仿真,分析其熱-結(jié)構(gòu)特性。

1 電磁場(chǎng)計(jì)算

1. 1 樣機(jī)基本參數(shù)

對(duì)樣機(jī)進(jìn)行電磁場(chǎng)仿真,樣機(jī)參數(shù)如表1所示,樣機(jī)如圖1所示。

表1 樣機(jī)參數(shù)

圖1 樣機(jī)

1. 2 損耗計(jì)算

對(duì)電機(jī)進(jìn)行電磁場(chǎng)仿真,計(jì)算導(dǎo)致電機(jī)發(fā)熱的電機(jī)運(yùn)行過(guò)程中產(chǎn)生的各種損耗,如式(1)所示[10]:

P∑=PCu1+PCu2+PFe+Pfw+PΔ(1)

式中:PΣ——總損耗;

PCu1——定子繞組的銅損;

PCu2——轉(zhuǎn)子銅損;

PFe——鐵心損耗;

Pfw——機(jī)械損耗;

PΔ——雜散損耗。

其中機(jī)械損耗和雜散損耗由經(jīng)驗(yàn)公式計(jì)算得到。損耗的仿真計(jì)算結(jié)果如表2所示。

表2 各部分損耗 W

2 水冷異步電機(jī)的三維瞬態(tài)溫度場(chǎng)

2. 1 求解域的模型

本文中的異步電機(jī)具有水冷套管結(jié)構(gòu),采用全封閉式外部冷卻系統(tǒng),無(wú)需內(nèi)部通風(fēng)。求解電機(jī)溫度場(chǎng)時(shí),為簡(jiǎn)化分析,作如下假設(shè):(1)繞組的絕緣漆分布均勻,繞組完全浸涂;(2)電機(jī)外殼的冷卻效果大致由水套的表面熱交換系數(shù)表示;(3)忽略輻射散熱以及轉(zhuǎn)子與轉(zhuǎn)軸之間的接觸熱阻;(4)電機(jī)各部分損耗不隨溫度變化而變化。

基于上述假設(shè),電機(jī)的求解域模型如圖2所示。

圖2 電機(jī)溫度場(chǎng)模型

為了簡(jiǎn)化分析,建立定子繞組等效模型,引入1個(gè)替代繞組銅線(xiàn)的等效導(dǎo)體和槽內(nèi)所有絕緣材料的等效絕緣層,如圖3所示。按參考文獻(xiàn)[11-12]可計(jì)算得槽絕緣的等效導(dǎo)熱系數(shù)為0.26 W/m·K。

圖3 等效繞組模型

2. 2 定轉(zhuǎn)子空氣間隙的有效熱導(dǎo)率

定轉(zhuǎn)子和氣隙之間的熱傳遞比較復(fù)雜,引入有效熱導(dǎo)率表示氣隙中流動(dòng)空氣的熱交換容量。有效熱導(dǎo)率λef計(jì)算如下[13]:

(2)

其中,η=R2/R1。

式中:Reg——空氣間隙的雷諾數(shù);

νa——空氣的運(yùn)動(dòng)黏度;

ωr——轉(zhuǎn)子的圓周速度;

δ——?dú)庀堕L(zhǎng)度;

R2——轉(zhuǎn)子的外半徑;

n——電機(jī)轉(zhuǎn)速;

Rec——臨界雷諾數(shù);

R1——定子的內(nèi)半徑。

2. 3 瞬態(tài)溫度場(chǎng)求解結(jié)果

2. 3. 1 額定工況運(yùn)行時(shí)溫度場(chǎng)仿真

環(huán)境溫度設(shè)為45 ℃。選定40 kW、3 600 r/ min額定工況持續(xù)運(yùn)行工作點(diǎn),對(duì)電機(jī)求解域模型進(jìn)行電磁場(chǎng)與溫度場(chǎng)的有限元耦合場(chǎng)仿真,按參考文獻(xiàn)[14-16]計(jì)算得出水冷機(jī)殼水套表面的等效換熱系數(shù)為15 508.8 W/(m2·K),按參考文獻(xiàn)[17]計(jì)算得定子繞組端的換熱系數(shù)為117.05 W/(m2·K)、定子鐵心端面的換熱系數(shù)77.64 W/(m2·K)、轉(zhuǎn)子鐵心端面的換熱系數(shù)51.61 W/(m2·K)、轉(zhuǎn)子端環(huán)的換熱系數(shù)205.36 W/(m2·K)、轉(zhuǎn)軸靠近轉(zhuǎn)子鐵芯端面換熱系數(shù)為27.32 W/(m2·K)。仿真得出電機(jī)的溫升分布,如圖4、圖5所示。

圖4 額定工況下運(yùn)行60 min時(shí)電機(jī)各部分溫度云圖

圖5 額定工況下運(yùn)行60 min時(shí)電機(jī)各部分溫度曲線(xiàn)

從圖4、圖5可知,電機(jī)在額定工況下持續(xù)運(yùn)行,溫升達(dá)到穩(wěn)態(tài)時(shí),繞組的溫度最高,繞組端部溫度明顯高于定子槽內(nèi)繞組,最高溫度為87 ℃。結(jié)構(gòu)上定子鐵心與機(jī)殼直接接觸,殼體中的冷卻水使定子鐵心軛部散熱優(yōu)于齒部,故定子鐵心齒部的溫度要高于軛部,最高為65.5 ℃。轉(zhuǎn)子部分溫度最高出現(xiàn)在端環(huán)與轉(zhuǎn)子導(dǎo)條連接處為56.7 ℃。由于機(jī)殼水道的原因,機(jī)殼溫度兩端明顯高于中部,轉(zhuǎn)軸的熱傳遞使得端蓋與軸承連接處溫升略高于端蓋其他區(qū)域,機(jī)殼最高溫度為70.8 ℃。由仿真結(jié)果可知,電機(jī)的最大溫升在合理范圍內(nèi),電機(jī)在額定運(yùn)行點(diǎn)可以安全可靠運(yùn)行。

2. 3. 2 峰值工況運(yùn)行時(shí)電機(jī)溫度場(chǎng)仿真

考慮實(shí)際電動(dòng)汽車(chē)可能運(yùn)行的工況,按照額定功率/額定轉(zhuǎn)速工作點(diǎn),峰值功率/峰值轉(zhuǎn)速工作點(diǎn)兩種工況依次進(jìn)行電機(jī)的溫度場(chǎng)仿真。電機(jī)在額定工況下60 min 后溫升基本達(dá)到穩(wěn)定,峰值工作點(diǎn)的溫度場(chǎng)仿真從電機(jī)額定工況運(yùn)行60 min后開(kāi)始,按技術(shù)指標(biāo)要求的60 s,校核電機(jī)峰值工作點(diǎn)的溫升。

圖6 峰值工況運(yùn)行60 s電機(jī)溫度云圖

圖7 峰值工況運(yùn)行60 s電機(jī)各部分溫度曲線(xiàn)

由圖6可知,電機(jī)峰值運(yùn)行60 s(總時(shí)間為3 660 s)時(shí),電機(jī)各部分的溫度分布與電機(jī)額定工況下溫度場(chǎng)分布相似。由圖7可知電機(jī)峰值工況運(yùn)行時(shí)繞組溫度上升最快,第60 s時(shí),繞組的最高溫度為107 ℃,定子鐵心的最高溫度為75.4 ℃,轉(zhuǎn)子鐵心的最高溫度為64.5 ℃,機(jī)殼的最高溫度為75 ℃,均在允許溫升限度以?xún)?nèi),滿(mǎn)足技術(shù)指標(biāo)的要求。

圖8 峰值工況運(yùn)行200 s電機(jī)溫度曲線(xiàn)

為了預(yù)測(cè)在此工況下電機(jī)峰值運(yùn)行的最長(zhǎng)允許時(shí)間,讓電機(jī)繼續(xù)在峰值工況下運(yùn)行一段時(shí)間。由圖8可知當(dāng)峰值工況下連續(xù)運(yùn)行200 s(總運(yùn)行3 880 s)時(shí),定子齒部最高溫度升高到74 ℃,轉(zhuǎn)子最高溫度升高到70 ℃,繞組的最高溫度達(dá)到123 ℃,超過(guò)了E級(jí)絕緣的最高允許溫度(120 ℃),接近繞組的最高溫升標(biāo)準(zhǔn)(80 K),所以電機(jī)在此環(huán)境下額定工況運(yùn)行一段時(shí)間后,峰值工況的運(yùn)行時(shí)間不宜超過(guò)200 s,否則易導(dǎo)致繞組絕緣損壞,影響電機(jī)運(yùn)行的安全性和可靠性。

3 三維熱應(yīng)力場(chǎng)的仿真研究

溫度變化產(chǎn)生的熱應(yīng)力會(huì)導(dǎo)致電機(jī)結(jié)構(gòu)發(fā)生形變,在電機(jī)的薄弱結(jié)構(gòu)部位較大的熱應(yīng)力可能導(dǎo)致電機(jī)結(jié)構(gòu)出現(xiàn)問(wèn)題。車(chē)用水冷異步電機(jī),轉(zhuǎn)子導(dǎo)條端部與端環(huán)的焊接部位和電機(jī)定子鐵心與機(jī)殼的過(guò)盈配合面由于材料熱膨脹系數(shù)不同在接合面產(chǎn)生的熱應(yīng)力可能會(huì)導(dǎo)致導(dǎo)條開(kāi)焊或過(guò)盈配合面漲開(kāi),產(chǎn)生安全隱患。

在溫度場(chǎng)分析的基礎(chǔ)上,借助相關(guān)彈性力學(xué)知識(shí),進(jìn)行熱應(yīng)力分析。力學(xué)彈性塑性問(wèn)題可以根據(jù)系統(tǒng)最小勢(shì)能的位移來(lái)解決:

Π=U+Ω(3)

式中:Π——系統(tǒng)總勢(shì)能;

U——系統(tǒng)內(nèi)部應(yīng)力變形勢(shì)能;

Ω——系統(tǒng)因外力作用而獲得的勢(shì)能。

對(duì)于電機(jī)熱應(yīng)力分析,可以認(rèn)為其變形遠(yuǎn)小于本身的幾何尺寸。故靜力學(xué)方程[18]如下:

(4)

式中:n——元素的數(shù)量;

δ——節(jié)點(diǎn)位移矢量矩陣;

D——彈性剛度矩陣;

B——基于單元形狀函數(shù)的應(yīng)變-位移矩陣;

ε0——熱應(yīng)變矢量矩陣;

αx、αy、αz——x、y、z方向上的材料熱膨脹系數(shù);

ΔT——溫升值;

σ——應(yīng)力矢量。

根據(jù)建立的數(shù)學(xué)模型,利用Ansys Workbench把峰值工況運(yùn)行60 s(總運(yùn)行3 660 s)的溫度場(chǎng)直接載荷到熱應(yīng)力場(chǎng)進(jìn)行仿真分析。

3. 1 轉(zhuǎn)子導(dǎo)條熱應(yīng)力場(chǎng)

電機(jī)轉(zhuǎn)子導(dǎo)條的熱應(yīng)力仿真結(jié)果如圖9所示。由圖9可以看出熱應(yīng)力最大值出現(xiàn)在轉(zhuǎn)子導(dǎo)條與端環(huán)連接的位置,最大值為31 MPa。圖10所示為轉(zhuǎn)子導(dǎo)條各部分熱應(yīng)力的軸向分布,在整體分布上同一截面導(dǎo)條底部熱應(yīng)力最大、頂部次之、中部最小,也可發(fā)現(xiàn)轉(zhuǎn)子導(dǎo)條與端環(huán)連接處的熱應(yīng)力明顯要比其他位置大。圖11所示為導(dǎo)條由于熱應(yīng)力產(chǎn)生的徑向位移分布圖,可以看出導(dǎo)條頂部的徑向位移比導(dǎo)條底部的徑向位移大,導(dǎo)條徑向形變的最大值也出現(xiàn)在端環(huán)與導(dǎo)條連接處(即導(dǎo)條端部),為38 μm。綜上可得轉(zhuǎn)子導(dǎo)條與端環(huán)連接處是轉(zhuǎn)子導(dǎo)條最容易斷裂的位置。

圖9 轉(zhuǎn)子導(dǎo)條熱應(yīng)力云圖

圖10 轉(zhuǎn)子導(dǎo)條熱應(yīng)力軸向分布曲線(xiàn)

圖11 導(dǎo)條徑向位移曲線(xiàn)

導(dǎo)條與端環(huán)的連接是焊接,按局部無(wú)損檢測(cè)最低標(biāo)準(zhǔn),焊接接頭強(qiáng)度與母材強(qiáng)度之比為0.8,可得焊接接頭處能承受的屈服強(qiáng)度為105 MPa,故此工況運(yùn)行下導(dǎo)條與端環(huán)連接處的焊接接頭不會(huì)開(kāi)焊,電機(jī)可安全運(yùn)行。

3. 2 鐵心與機(jī)殼過(guò)盈配合面量的選取

機(jī)殼與鐵心過(guò)盈配合面的熱應(yīng)力如圖12、圖13所示。由圖12、圖13可以看出過(guò)盈配合面熱應(yīng)力軸向分布不均,但同一截面相差不大,熱應(yīng)力最大處在端部,為32 MPa。

圖12 過(guò)盈配合面熱應(yīng)力云圖

圖13 過(guò)盈配合面熱應(yīng)力軸向分布曲線(xiàn)

為避免機(jī)殼與鐵心之間出現(xiàn)滑移或漲開(kāi)現(xiàn)象,需要合理選擇過(guò)盈配合量的大小,以保證接觸面的過(guò)盈應(yīng)力大于熱應(yīng)力。配合面的過(guò)盈應(yīng)力可采用式(5)進(jìn)行近似計(jì)算:

(5)

式中:σr——徑向應(yīng)力;

b——被包容件外壁半徑;

P——內(nèi)壁(或外壁)壓力;

a——被包容件內(nèi)壁半徑;

r——過(guò)盈配合面半徑;

σθ——切向應(yīng)力;

E——材料的彈性模量;

δ——有效過(guò)盈量;

c——包容件的外壁半徑。

由式(5)計(jì)算徑向過(guò)盈應(yīng)力等于熱應(yīng)力時(shí)的有效過(guò)盈量為0.018 mm,計(jì)算0.018 mm過(guò)盈量對(duì)應(yīng)的切向過(guò)盈應(yīng)力為88.49 MPa,遠(yuǎn)大于過(guò)盈配合面的切向熱應(yīng)力。

機(jī)殼與鐵心的過(guò)盈配合不需要傳遞扭矩,為保證配合安全可靠,過(guò)盈配合的過(guò)盈量往往取值很小,機(jī)殼與鐵心在過(guò)盈配合接觸面處于彈性范圍。按第四強(qiáng)度理論,過(guò)盈配合不產(chǎn)生塑性變形所容許的最大過(guò)盈量計(jì)算公式:

(6)

式中:pfmax——過(guò)盈配合聯(lián)結(jié)件不產(chǎn)生塑性變形容許最大結(jié)合壓力的最小值;

df——過(guò)盈配合面直徑;

Ea——被包容件彈性模量;

Ei——包容件彈性模量;

Ca——被包容件系數(shù);

Ci——包容件系數(shù)。

由式(6)計(jì)算得出使機(jī)殼不發(fā)生塑性變形的最大過(guò)盈量為0.22 mm。

綜上所述,保證電機(jī)在額定工況下連續(xù)運(yùn)行60 min后,又在峰值工況下運(yùn)行60 s的情況下定子結(jié)構(gòu)熱態(tài)剛度不發(fā)生較大改變的最小過(guò)盈量為0.018 mm,使機(jī)殼不發(fā)生塑性變形的最大過(guò)盈量為0.22 mm。

4 結(jié) 語(yǔ)

本文基于Ansys多物理場(chǎng)仿真平臺(tái)對(duì)1臺(tái)額定功率40 kW的車(chē)用異步電機(jī)進(jìn)行多物理場(chǎng)耦合的有限元仿真,分析了電機(jī)連續(xù)運(yùn)行在額定功率和峰值功率兩種熱負(fù)荷工況時(shí)電機(jī)結(jié)構(gòu)部件的溫升和熱應(yīng)力,校核了電機(jī)繞組絕緣的安全性,預(yù)測(cè)了電機(jī)峰值工況的最長(zhǎng)允許時(shí)間。通過(guò)有限元仿真分析了電機(jī)轉(zhuǎn)子導(dǎo)條端部焊點(diǎn)和機(jī)殼與定子鐵心過(guò)盈配合面的熱應(yīng)力,校核了電機(jī)轉(zhuǎn)子端部導(dǎo)條焊點(diǎn)是否存在開(kāi)焊問(wèn)題,計(jì)算了保證電機(jī)在最高溫升工況下過(guò)盈配合面不發(fā)生滑移或漲開(kāi)的過(guò)盈配合量取值。

【參 考 文 獻(xiàn)】

[1] EHSANI M, YIMIN G, EGAY S. 現(xiàn)代電動(dòng)汽車(chē)混合動(dòng)力電動(dòng)汽車(chē)和燃料電池車(chē):基本原理,理論和設(shè)計(jì)[M].倪光正,倪培宏,熊素銘,譯.北京:機(jī)械工業(yè)出版社,2008.

[2] LIU Y J, LEE Y S, JUNG H K, et al. 3D thermal stress analysis of the rotor of an induction motor[J]. IEEE Transactions on Magnetics,2000,36(4): 1394-1397.

[3] MROZ J. Temperature field of a double squirrel-cage motor during start up[J]. IEEE proceedings: Electric Power Applications,2005,152(6): 1531-1538.

[4] 錢(qián)榮超.復(fù)合籠條感應(yīng)電動(dòng)機(jī)三維熱應(yīng)力場(chǎng)與起動(dòng)性能計(jì)算與分析[D].哈爾濱:哈爾濱理工大學(xué),2012.

[5] 王艷武,楊立,陳翾,等.異步電機(jī)轉(zhuǎn)子三維溫度場(chǎng)及熱應(yīng)力場(chǎng)研究[J].電機(jī)與控制學(xué)報(bào),2010,14(6): 27-32.

[6] BANGURA J F, DEMERDASH N A. Diagnosis and characterization of effects of broken bars and connectors in squirrel-cage induction motor by time-stepping coupled finite element state space modeling approach[J]. IEEE Trans Energy Convers,1999,14(4): 1167-1176.

[7] 謝穎,王澤,單雪婷,等.基于多場(chǎng)量的籠型感應(yīng)電機(jī)三維瞬態(tài)磁熱固耦合計(jì)算分析[J].中國(guó)電機(jī)工程學(xué)報(bào),2016,36(11): 3076-3084.

[8] 湯蘊(yùn)繆,梁艷萍.電機(jī)電磁場(chǎng)的分析與計(jì)算[M].北京:機(jī)械工業(yè)出版社,2012.

[9] 趙博,張洪亮.Ansoft 12在工程電磁場(chǎng)中的應(yīng)用[M].北京:中國(guó)水利水電出版社,2010.

[10] LIU R, ZHENG P, XIE D, et al. Research on the high power density electromagnetic propeller[J]. IEEE Transactions on Magnetics,2007,43(1): 355-358.

[11] XYPTRAS J, HATZIATHANASSIOU V. Thermal analysis of an electrical machine taking into account the iron losses and the deep-bar effect[J].IEEE Transactions on Energy Conversion,1999,14(4): 996-1003.

[12] 翁榮周.傳熱學(xué)的有限元方法[M].廣州:暨南大學(xué)出版社,2000.

[13] LI W L, LI S F, XIE Y, et al. Stator-rotor coupled thermal field numerical calculation of induction motors and correlated factors sensitivity analysis[J]. Proceedings of the CSEE,2007,27(24): 85-91.

[14] STATON D, BOGLIETTI A, CAVAGNINO A. Solving the more difficult aspectsof electric motor thermal analysis[C]∥ IEEE International Electric Machinesand Drives Conference, 2003 IEMDC'03, 2003: 747-755.

[15] STATON D A, CAVAGNINO A. Convection heat transfer and flow calculations suitable for electric machines thermal models[J]. IEEE Transactions on Industrial Electronics,2008,55(10): 3509-3516.

[16] 郁嵐,衛(wèi)運(yùn)鋼,杜雅琴,等.熱工基礎(chǔ)及流體力學(xué)[M].北京:中國(guó)電力出版社,2006.

[17] 黃國(guó)治,傅豐禮.中國(guó)旋轉(zhuǎn)電機(jī)設(shè)計(jì)手冊(cè)[M].北京:中國(guó)電力出版社,2007.

[18] COOK R D. Finite element modeling for stress analysis[M]. New York:John Wiley & Sons,1995.

主站蜘蛛池模板: 成人福利在线免费观看| 亚洲热线99精品视频| 99在线观看视频免费| 少妇精品网站| 久久精品91麻豆| 午夜啪啪网| 日韩东京热无码人妻| 日韩欧美国产综合| 偷拍久久网| 日韩精品一区二区三区中文无码| 五月婷婷导航| 国产二级毛片| 777午夜精品电影免费看| 久久久噜噜噜| 国产污视频在线观看| 成人在线视频一区| 草草影院国产第一页| 无码精品福利一区二区三区| 亚洲日本在线免费观看| 国产亚洲高清视频| 国产精品午夜福利麻豆| 欧美自慰一级看片免费| 日韩AV手机在线观看蜜芽| 欧美激情视频一区| 伊人久久婷婷| 九九免费观看全部免费视频| 国产鲁鲁视频在线观看| 国产成人精品在线1区| 中国精品久久| 国产一线在线| 国产精鲁鲁网在线视频| 亚洲精品天堂在线观看| 国产精品网拍在线| 麻豆精品国产自产在线| jizz在线免费播放| 综合社区亚洲熟妇p| 99久久国产精品无码| 啪啪永久免费av| 中文字幕免费在线视频| 中文字幕亚洲另类天堂| 免费a级毛片视频| 日韩美一区二区| 亚洲中文字幕久久精品无码一区| av无码一区二区三区在线| 日本高清免费不卡视频| 91福利在线看| 国产精品男人的天堂| 国内精品自在欧美一区| 无码网站免费观看| 在线观看视频99| 欧美性久久久久| 国产青青操| 欧美人与牲动交a欧美精品| 十八禁美女裸体网站| 亚洲成人高清无码| 欧美亚洲网| 亚洲欧洲日产无码AV| 伊人网址在线| 91黄视频在线观看| 久草视频福利在线观看| a级毛片免费看| 免费毛片在线| 国产综合另类小说色区色噜噜| 在线看AV天堂| 伊人蕉久影院| 中文字幕亚洲乱码熟女1区2区| 精品国产福利在线| 国产成人精品一区二区免费看京| 播五月综合| 97视频精品全国免费观看| 日韩在线成年视频人网站观看| 无码精品国产dvd在线观看9久 | 人妻夜夜爽天天爽| 拍国产真实乱人偷精品| 色综合久久无码网| 国产亚洲欧美在线视频| 青青草综合网| 欧洲精品视频在线观看| 国产麻豆91网在线看| 国产玖玖玖精品视频| 国产欧美日韩视频怡春院| 亚洲免费福利视频|