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

基于Gam bit前處理的氣膜冷卻火焰筒壁溫分析

2017-11-13 04:53:49扈鵬飛
航空發(fā)動(dòng)機(jī) 2017年2期
關(guān)鍵詞:有限元

程 明,扈鵬飛,萬 斌,常 峰

(中國航發(fā)沈陽發(fā)動(dòng)機(jī)研究所,沈陽110015)

基于Gam bit前處理的氣膜冷卻火焰筒壁溫分析

程 明,扈鵬飛,萬 斌,常 峰

(中國航發(fā)沈陽發(fā)動(dòng)機(jī)研究所,沈陽110015)

為更好實(shí)現(xiàn)航空發(fā)動(dòng)機(jī)燃燒室氣膜冷卻結(jié)構(gòu)火焰筒的壁溫分析和冷卻結(jié)構(gòu)優(yōu)化,針對(duì)原有的火焰筒2維壁溫計(jì)算程序開展了2次開發(fā)工作,形成了基于G A M BIT前處理的火焰筒壁溫分析程序。新的計(jì)算過程采用G ambit軟件對(duì)幾何模型進(jìn)行前處理,生成三角形網(wǎng)格、指定邊界條件分組,通過編制前處理模塊代碼,對(duì)導(dǎo)出的網(wǎng)格文件進(jìn)行解析,進(jìn)一步將各種信息導(dǎo)入已有的有限元壁溫計(jì)算程序,完成壁溫計(jì)算分析。給出了采用以上方法進(jìn)行火焰筒壁溫計(jì)算和優(yōu)化的實(shí)例,結(jié)果表明:該方法較為有效地克服了原方法的各項(xiàng)缺點(diǎn),其交互性強(qiáng),大幅度提高了工作效率。

火焰筒;壁溫分析;有限元G ambit;前置處理;航空發(fā)動(dòng)機(jī)

0 引言

火焰筒是航空發(fā)動(dòng)機(jī)最重要的受熱部件之一。隨著燃燒室溫升的不斷提高,火焰筒的工作環(huán)境變得極為惡劣,在極大的熱負(fù)荷以及嚴(yán)重的熱沖擊下,火焰筒承受著很大的熱應(yīng)力、蠕動(dòng)應(yīng)力和疲勞應(yīng)力。因此,必須采取有效措施,對(duì)火焰筒壁面進(jìn)行冷卻,以降低火焰筒壁溫水平及溫度梯度,保證火焰筒的使用壽命[1-2]。

氣膜冷卻是1種主動(dòng)防護(hù)性冷卻技術(shù)[1-2]。典型的燃?xì)廨啓C(jī)燃燒室氣膜冷卻結(jié)構(gòu)中,冷卻氣流沿壁面切線方向以一定的入射角向高溫燃?xì)庵袊娚洌谥髁鞯膲毫湍Σ亮ψ饔孟拢鋮s氣流黏附在壁面附近形成溫度較低的冷氣膜,將高溫燃?xì)馀c壁面隔離,避免了高溫燃?xì)庵苯訉?duì)壁面進(jìn)行對(duì)流換熱,并能將一部分高溫燃?xì)饣虬l(fā)光火焰對(duì)壁面的輻射熱量帶走,從而對(duì)壁面起到良好的防護(hù)作用。通常沿著火焰筒軸向大約30~50 mm的間距就要引入下一段冷卻進(jìn)氣。

葛紹巖、劉登瀛[3]給出了氣膜冷卻結(jié)構(gòu)火焰筒內(nèi)外表面對(duì)流換熱系數(shù)、氣膜冷卻效率等的經(jīng)驗(yàn)準(zhǔn)則,使火焰筒計(jì)算域的能量方程得到封閉,可用于壁溫計(jì)算。在火焰筒壁溫計(jì)算中,一般采用有限元分析方法[4-5]。將火焰筒內(nèi)/外壁面抽象為2維軸對(duì)稱計(jì)算域,對(duì)其進(jìn)行三網(wǎng)格劃分,加載對(duì)流換熱、輻射換熱等邊界條件,可獲得詳細(xì)的壁面溫度分布。這種方法已經(jīng)在工程實(shí)踐中得到了應(yīng)用[6-7]。

楊志民等[8]對(duì)氣膜冷卻結(jié)構(gòu)火焰筒開展了3維壁溫計(jì)算,考慮了復(fù)雜幾何結(jié)構(gòu)及火焰筒內(nèi)進(jìn)氣不均的影響。張勃、吉洪湖[9]采用流熱耦合的方法進(jìn)行實(shí)際燃燒室內(nèi)燃燒場(chǎng)和氣膜冷卻火焰筒壁面溫度場(chǎng)的3維耦合計(jì)算,可獲得更為詳盡的信息,但是由于工作量較大,很顯然這種方法暫時(shí)還不適用于方案設(shè)計(jì)階段。另外,隨著燃燒室設(shè)計(jì)技術(shù)的發(fā)展,更先進(jìn)的多斜孔、層板等冷卻結(jié)構(gòu)已經(jīng)在航空發(fā)動(dòng)機(jī)上得到了應(yīng)用[2,10],但目前仍有相當(dāng)數(shù)量的在役航空發(fā)動(dòng)機(jī)采用機(jī)加環(huán)氣膜冷卻結(jié)構(gòu)。這樣,針對(duì)機(jī)加環(huán)氣膜冷卻結(jié)構(gòu)火焰筒、基于經(jīng)驗(yàn)換熱公式的有限元壁溫計(jì)算方法仍有進(jìn)一步研究的需求。

在以前的有限元壁溫計(jì)算中,針對(duì)三角形有限元網(wǎng)格的劃分往往采用手工操作的方式,邊界單元換熱條件的加載也是手工進(jìn)行的。常用做法是:(1)將壁面2維模型放大5~10倍,并繪制在坐標(biāo)紙上。在紙上人工劃分為三角形網(wǎng)格,并讀出網(wǎng)格節(jié)點(diǎn)坐標(biāo),人工進(jìn)行節(jié)點(diǎn)/線元/單元編號(hào)和組合,并錄入數(shù)據(jù)文件;(2)人為指定計(jì)算域內(nèi)的邊界單元的邊界條件類型和計(jì)算換熱量的特征參數(shù)。該做法的缺點(diǎn)為:一方面工作效率很低,劃分的網(wǎng)格數(shù)量/尺度不可控,點(diǎn)坐標(biāo)讀取精度較差/易出低級(jí)錯(cuò)誤;另一方面整個(gè)壁溫分析的過程不能實(shí)現(xiàn)交互式處理,不適應(yīng)需要多方案優(yōu)化時(shí)的工程需求。

本文開發(fā)了1種基于Gambit前處理的氣膜冷卻結(jié)構(gòu)火焰筒壁溫分析方法,可以對(duì)火焰筒冷卻方案快速實(shí)現(xiàn)幾何建模、網(wǎng)格生成、邊界條件加載和壁溫的數(shù)值模擬,成熟可靠、交互性強(qiáng),有效提高了火焰筒冷卻結(jié)構(gòu)優(yōu)化設(shè)計(jì)工作的效率。

1 計(jì)算方法和計(jì)算模型

在不考慮火焰筒頭部/旋流器的情況下,將火焰筒內(nèi)/外壁面抽象為2維軸對(duì)稱計(jì)算域。火焰筒壁溫計(jì)算幾何模型如圖1所示。

在幾何建模和網(wǎng)格劃分時(shí)不考慮氣膜孔、主燃孔和摻混孔的存在,但是在邊界條件加載時(shí),要考慮開孔流量對(duì)火焰筒燃?xì)鈧?cè)、冷氣側(cè)熱力參數(shù)的影響。對(duì)應(yīng)基于2維軸對(duì)稱穩(wěn)態(tài)導(dǎo)熱微分方程[3]

圖1 火焰筒壁溫計(jì)算幾何模型

其邊界條件為

第Ⅰ類:T|Γr=T0

第Ⅱ類:-λ/(?T/?n)|Γq=q

第Ⅲ類:-λ/(?T/?n)|Γh=h(T-T0)。式(1)的泛函表達(dá)式為

式中:x、r為軸向、徑向坐標(biāo);λ為導(dǎo)熱系數(shù);T為溫度;T0為周圍介質(zhì)溫度;q為熱流密度;h為換熱系數(shù);ΓT、Γq、Γh分別為第Ⅰ、Ⅱ、Ⅲ類邊界。

將計(jì)算對(duì)象劃分為Ne個(gè)3節(jié)點(diǎn)三角形環(huán)元素。元素的3個(gè)節(jié)點(diǎn)分別為i、j、m,則單元內(nèi)溫度插值函數(shù)可表示為

其中 [N]=[Ni,Nj,Nk],{T}e={Ti,Yj,Tk}τ

溫度插值函數(shù)為

式中:t=s/L,s為距節(jié)點(diǎn)的距離;邊界線段i、j的長(zhǎng)度L=(xi-xj)2+(ri-rj)2。

將單元內(nèi)及邊界線段上的溫度分布公式帶入泛函表達(dá)式(2),可得到Je及(?Je)/?{T}e的表達(dá)式。對(duì)于計(jì)算域內(nèi)的任一點(diǎn)M,把所有與M點(diǎn)有關(guān)的元素進(jìn)行單元分析,將其系數(shù)和常數(shù)部分按 ?J/?TM=Σ(?Je/?TM)建立代數(shù)方程。

對(duì)所有節(jié)點(diǎn),得到代數(shù)方程組

采用Cholesky法求解方程組(6)可得出整個(gè)計(jì)算域上的溫度分布。

2 邊界條件

針對(duì)圖1中火焰筒壁面模型,火焰筒壁面換熱表面分類如圖2所示。圖中給出了典型氣膜段及其各換熱表面的特征[2]:

(1)第 1類:壁面接受高溫燃?xì)獾妮椛浜蜌饽だ鋮s的對(duì)流換熱;

(2)第 2類:壁面向燃燒室機(jī)匣輻射熱,同時(shí)與2股腔道的冷卻空氣之間進(jìn)行對(duì)流換熱;

(3)第3、4類:壁面承受冷卻空氣與壁的對(duì)流換熱及二者之間的輻射換熱的復(fù)合作用;

(4)第5類:近似認(rèn)為與外界的熱交換為0。

對(duì)以上各類表面,其具體的換熱方程見文獻(xiàn)[3-6]。在計(jì)算過程中所采用的開孔流量分配、火焰筒燃?xì)鈧?cè)/冷氣側(cè)沿程熱力參數(shù)使用1維流量分配程序的計(jì)算結(jié)果。

圖2 火焰筒壁面換熱表面分類

3 半自動(dòng)化計(jì)算過程的前置處理和壁溫計(jì)算

3.1 模型導(dǎo)入和網(wǎng)格劃分過程

(1)在幾何建模軟件中,對(duì)壁面型線進(jìn)行處理(圖2),并以acis文件的格式導(dǎo)出。

(2)打開Gambit軟件[11]:導(dǎo)入幾何模型(使用前一步驟生成的acis文件);對(duì)型線進(jìn)行必要的切割、合并處理;指定網(wǎng)格劃分的尺度,完成三角形網(wǎng)格劃分;選定求解器類型為polyflow;為后續(xù)計(jì)算加載邊界條件方便,對(duì)如圖2所示的5類表面,按特定的規(guī)則指定所有的邊界線的名稱;導(dǎo)出網(wǎng)格文件(其擴(kuò)展名為“neu”)。

火焰筒壁溫計(jì)算的網(wǎng)格示例如圖3所示。

3.2 “neu”網(wǎng)格文件解析

“neu”文件中包含了所有網(wǎng)格的節(jié)點(diǎn)、線元、三角形單元的信息以及邊界線元的邊界條件分組等信息。“neu”文件所含數(shù)據(jù)信息解析見表1。同時(shí),該文件中還包含了一些與本文計(jì)算無關(guān)的信息,表中未全部列出。

圖3 火焰筒壁溫計(jì)算的網(wǎng)格示例

3.3 有限元壁溫計(jì)算

表1 "neu"文件所含數(shù)據(jù)信息解析

在有限元壁溫計(jì)算程序中,讀入前述的neu文件中的數(shù)據(jù)及由1維性能計(jì)算程序獲得的相關(guān)開孔流量分配、沿程熱力參數(shù)數(shù)據(jù),即可通過迭代計(jì)算獲得火焰筒壁溫分布。

在火焰筒流路方案發(fā)生變化、導(dǎo)致幾何模型變化時(shí),只需重新運(yùn)行第3.1節(jié)中的“模型導(dǎo)入和網(wǎng)格劃分過程”,就可迅速得到新模型對(duì)應(yīng)的有限元網(wǎng)格及邊界分組信息,結(jié)合前期已實(shí)現(xiàn)交互式自動(dòng)化運(yùn)行的“1維性能計(jì)算程序”,可獲得新方案的火焰筒壁溫分布。

采用上述方法,1個(gè)全新方案的氣膜冷卻火焰筒的壁溫計(jì)算時(shí)間可以由以前的1~2天縮短為現(xiàn)在10~30 min,從而有效提高了工作效率。

4 改進(jìn)方案壁溫計(jì)算實(shí)例

圖4 火焰筒外壁初始方案和改進(jìn)方案

火焰筒外壁的初始方案和改進(jìn)方案如圖4所示。根據(jù)燃燒室總體設(shè)計(jì)要求,與初始方案相比,改進(jìn)型燃燒室火焰筒的長(zhǎng)度和出口外徑略有增加,而火焰筒進(jìn)口尺寸和冷卻空氣量基本不變。對(duì)此,以燃燒室集成設(shè)計(jì)系統(tǒng)為平臺(tái),開展了多方案燃燒室的優(yōu)化設(shè)計(jì),主要工作內(nèi)容包括流路和流量分配設(shè)計(jì)、火焰筒開孔方案及氣膜段設(shè)計(jì)、2維壁溫計(jì)算等。如圖4所示的改進(jìn)方案為最終的外壁優(yōu)化設(shè)計(jì)結(jié)果。其特征是,除滿足前述燃燒室總體設(shè)計(jì)要求以外,仍采用氣膜冷卻結(jié)構(gòu),且氣膜段數(shù)保持不變,但是各氣膜段的長(zhǎng)度與初始方案相比有所變化,各段氣膜的冷卻空氣量在初始方案基礎(chǔ)上進(jìn)一步優(yōu)化、調(diào)整,以達(dá)到火焰筒壁面溫度低、溫度梯度小的目的。

地面起飛狀態(tài)初始方案和改進(jìn)方案火焰筒外壁壁溫計(jì)算結(jié)果如圖5所示。

圖5 初始方案和改進(jìn)方案火焰筒外壁壁溫計(jì)算結(jié)果

從圖中可見:

(1)在每個(gè)氣膜段內(nèi),火焰筒壁溫呈現(xiàn)先低后高的過程。“先低”是由于初始段氣膜冷卻的效率最高,“后高”是由于隨著壁面氣膜與火焰筒內(nèi)燃?xì)獾膿交欤瑲饽だ鋮s效率逐漸降低。

(2)與初始方案相比,由于改進(jìn)方案火焰筒氣膜段長(zhǎng)度的調(diào)整以及冷卻空氣量分配的優(yōu)化,其壁溫最大值有明顯下降。同時(shí),各氣膜段的溫度最高值(在該氣膜段末尾處)的離散度大幅度減小。

5 結(jié)束語

(1)針對(duì)氣膜冷卻火焰筒壁溫計(jì)算所采用的傳統(tǒng)有限元方法,開展了二次開發(fā)工作,形成了基于商業(yè)軟件前處理的火焰筒壁溫分析軟件。通過Gambit軟件幾何模型前處理、網(wǎng)格生成、邊界條件分組、網(wǎng)格文件解析、前處理模塊代碼的編制等,進(jìn)一步完成了壁溫計(jì)算分析。通過文中給出的火焰筒壁溫計(jì)算和優(yōu)化的實(shí)例,表明該方法較為有效地克服了原方法的各項(xiàng)缺點(diǎn),成熟可靠、交互性強(qiáng),大幅度提高了工作效率。

(2)對(duì)文中涉及的火焰筒外壁初始設(shè)計(jì)方案,通過流路和流量分配設(shè)計(jì)、火焰筒開孔方案及氣膜段設(shè)計(jì)調(diào)整,在冷卻空氣總量不變的設(shè)計(jì)要求下可實(shí)現(xiàn)火焰筒壁面溫度低、溫度梯度減小的目標(biāo)。

[1]金如山.航空燃?xì)廨啓C(jī)燃燒室 [M].北京:中國宇航出版社,1988:340-342.JIN Rushan.Aero-turbine combustion [M].Beijing:China Astronautic Publishing House,1988:340-342.(in Chinese)

[2]Lefebvre A H,Ballal D R.Gas turbine combustion-alternative fuels and emissions[M].Third edition,Boca Ratou:CRC Press,Taylor&Francis Group,2010:20-33,315-355.

[3]葛紹巖,劉登瀛,徐靖中,等.等,氣膜冷卻[M].北京:科學(xué)出版社,1985:25-61.GE Shaoyan,LIU Dengying,XU Jingzhong,et a1.Film cooling[M].Beijing:Science Press,1985:25-61.(in Chinese)

[4]韓振興,朱谷君,冀守禮,等.氣膜冷卻燃燒室火焰筒2維壁溫分布計(jì)算[J].航空動(dòng)力學(xué)報(bào),1995,10,(1):83-86.HAN Zhenxing,ZHU Gujun,JI Shouli,et al.Numerical calculation of 2 dimension wall temperature distribution of flame tube of film cooling combustion[J].Journal of Aerospace Power.(in Chinese)

[5]朱谷君.熱傳導(dǎo)及流體流動(dòng)的數(shù)值解法[D].北京:北京航空航天大學(xué),1992.ZHU Gujun.Numerical method for heat conduction and fluid flow[D].Beijing:Beihang University,1992.(in Chinese)

[6]蘇克.某型發(fā)動(dòng)機(jī)燃燒室火焰筒壁溫分析[J].航空發(fā)動(dòng)機(jī),1993(3):22-31.SU Ke.Wall temperature analysis of an aeroengine flame tube[J].Aeroengine,1993(3):22-31.(in Chinese)

[7]朱長(zhǎng)青,董志銳.鉆孔式氣膜冷卻火焰筒壁溫計(jì)算 [J].航空學(xué)報(bào),1991,(03):206-209 ZHU Changqing,DONG Zirui.Numerical analysis of wall temperature of film cooling combutor flame tube[J].Aeroengine,`1991(03):206-209(in Chinese)

[8]楊志民,韓振興.燃燒室火焰筒壁溫三維數(shù)值分析[J].航空發(fā)動(dòng)機(jī),1998,24(3):23-27.YANG Zhimin,HAN Zhenxing.Three-dimensional numerical analysis of wall temperature of combustor flame tube[J].Aeroengine,1998,24(3):23-27.(in Chinese)

[9]張勃,吉洪湖等.多斜孔壁與機(jī)加環(huán)氣膜冷卻燃燒室的壁面換熱特性數(shù)值研究 [J].航空動(dòng)力學(xué)報(bào),2012(4):832-836.ZHANG Bo,JI Honghu.Wall heat exchange property analysis of effusion cooling and machining ring film cooling structure flame tube[J].Journal of Aerospace Power.2012(4):832-836.(in Chinese)

[10]左渝鈺,張寶誠.航空發(fā)動(dòng)機(jī)主燃燒室火焰筒壁冷卻的研究[J].航空發(fā)動(dòng)機(jī),2002,28(4):38-43.ZUO Yuyu,ZHANG Baocheng.Investigation on cooling of flam tube wall of aeroengine combustor[J].Aeroengine,2002,28(4):38-43.(in Chinese)

[11]Fluent Inc.GAMBIT 2.0 User's Guide[M].Fluent Inc.,2003.

Wall Temperature Analysis of Film Cooling Flame Tube Based on Gambit

CHENG M ing,HU Pengfei,WAN Bin,CHANG Feng
(AECC Shenyang Engine Research Institute,Shengyang 110015,China)

In order to effectively realize wall temperature analysis and cooling structure optimization of film cooling flame tube for aeroengine combustor,secondary development work was carried out on an existing two-dimensional temperature calculation program for flame tube wall.A new flame tube wall temperature analysis program generates based on Gambit.The new calculation procedure deal with geometric model,form a triangle mesh and specify the boundary conditions using Gambit.Then,use a developed preprocessing module code to explain the exported mesh file,and various needed information was imported into the existing wall temperature calculation program.An example was given to calculate and optimize the flame tube wall temperature using the method.The result shows that the method is more effective and overcomes the shortcomings of original method,with better efficiency and strong interaction.

flame tube;wall temperature analysis;finite element;Gambit;aeroengine

V 231.1+3

A

1 0.1 3477/j.cnki.aeroengine.201 7.02.002

2016-09-30 基金項(xiàng)目:國家重大基礎(chǔ)研究項(xiàng)目資助

程明(1971),男,博士,自然科學(xué)研究員,主要從事航空發(fā)動(dòng)機(jī)主燃燒室技術(shù)研究工作;E-mail:cm.ln@hotmail.com。

程明,扈鵬飛,萬斌,等.基于 Gambit前處理的火焰筒壁溫分析[J].航空發(fā)動(dòng)機(jī),2017,43(2):6-9.CHENG Ming,HU Pengfei,WAN bin,etal.Wall temperature analysis offilm cooling flame tube based on Gambit[J].Aeroengine,2017,43(2):6-9.

(編輯:張寶玲)

猜你喜歡
有限元
基于擴(kuò)展有限元的疲勞裂紋擴(kuò)展分析
非線性感應(yīng)加熱問題的全離散有限元方法
TDDH型停車器制動(dòng)過程有限元分析
新型有機(jī)玻璃在站臺(tái)門的應(yīng)用及有限元分析
基于I-DEAS的履帶起重機(jī)主機(jī)有限元計(jì)算
基于有限元模型對(duì)踝模擬扭傷機(jī)制的探討
10MN快鍛液壓機(jī)有限元分析
磨削淬硬殘余應(yīng)力的有限元分析
基于SolidWorks的吸嘴支撐臂有限元分析
箱形孔軋制的有限元模擬
上海金屬(2013年4期)2013-12-20 07:57:18
主站蜘蛛池模板: 国产精品密蕾丝视频| 免费看的一级毛片| 成年人国产网站| 亚洲综合极品香蕉久久网| 亚洲AV无码久久天堂| AV不卡国产在线观看| 91精品免费高清在线| 日本国产在线| 亚洲 欧美 偷自乱 图片| 国产JIZzJIzz视频全部免费| 国产超薄肉色丝袜网站| www.91中文字幕| 久青草国产高清在线视频| 中文字幕 日韩 欧美| 久久a级片| 亚洲欧美精品日韩欧美| 久久国产亚洲偷自| 国产成人亚洲精品色欲AV | 久久香蕉国产线看观看式| 美女被躁出白浆视频播放| 国产美女在线观看| 国产毛片久久国产| 国产美女精品一区二区| 怡红院美国分院一区二区| 国产亚洲视频中文字幕视频| 日韩av无码DVD| 国产麻豆精品久久一二三| 精品欧美日韩国产日漫一区不卡| 国产精品网址你懂的| 亚洲国产中文在线二区三区免| 婷婷六月天激情| 欧美一级视频免费| 91在线精品麻豆欧美在线| 亚洲无码高清免费视频亚洲| 日韩国产高清无码| 国产靠逼视频| 丁香综合在线| 亚洲综合18p| 久久青草精品一区二区三区| 国产资源免费观看| 99精品福利视频| 欧美精品一二三区| 亚洲欧美日韩综合二区三区| 九月婷婷亚洲综合在线| 欧美成人午夜影院| 91极品美女高潮叫床在线观看| 日本成人不卡视频| 国产在线自揄拍揄视频网站| 成人看片欧美一区二区| 国产成人av一区二区三区| 国产菊爆视频在线观看| 日韩国产欧美精品在线| 亚洲AⅤ综合在线欧美一区| 日韩欧美国产综合| 国产亚洲欧美日韩在线一区| 97se亚洲综合| 久久久精品无码一二三区| 久久午夜夜伦鲁鲁片无码免费| 69免费在线视频| 性69交片免费看| 日本午夜精品一本在线观看| 精品伊人久久大香线蕉网站| 国产电话自拍伊人| 免费看一级毛片波多结衣| 久久免费成人| 伊人中文网| 亚洲91在线精品| 女人毛片a级大学毛片免费| 国产高清在线丝袜精品一区| 亚洲色图欧美视频| 十八禁美女裸体网站| 国产9191精品免费观看| 亚洲第一中文字幕| 5388国产亚洲欧美在线观看| 久久精品无码中文字幕| 日韩天堂网| 在线精品欧美日韩| 日本欧美在线观看| 国产永久无码观看在线| 全午夜免费一级毛片| 亚洲欧美在线精品一区二区| 亚洲国产清纯|