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

基于ANSYS CFX的魚雷渦輪機(jī)通流流場及性能計(jì)算方法

2012-05-28 03:02:54伊進(jìn)寶師海潮
關(guān)鍵詞:模型

陳 剛, 伊進(jìn)寶, 師海潮

?

基于ANSYS CFX的魚雷渦輪機(jī)通流流場及性能計(jì)算方法

陳 剛1, 伊進(jìn)寶2, 師海潮2

(1. 海軍裝備部駐西安地區(qū)軍事代表局, 陜西 西安, 710054; 2. 中國船舶重工集團(tuán)公司第705研究所, 陜西 西安, 710075)

運(yùn)用計(jì)算流體力學(xué)(CFD)軟件ANSYS CFX對(duì)魚雷渦輪機(jī)通流部分流場及性能進(jìn)行了數(shù)值計(jì)算, 介紹了CFD技術(shù)應(yīng)用到魚雷燃?xì)鉁u輪機(jī)通流設(shè)計(jì)中的一般步驟, 旨在提供一種燃?xì)鉁u輪機(jī)通流部分流場及性能的預(yù)測方法, 探索CFD技術(shù)在魚雷燃?xì)鉁u輪機(jī)動(dòng)力系統(tǒng)設(shè)計(jì)中的應(yīng)用途徑。研究結(jié)果表明, 利用CFD技術(shù)計(jì)算魚雷燃?xì)鉁u輪機(jī)通流流場, 預(yù)測通流氣動(dòng)性能與試驗(yàn)數(shù)據(jù)吻合較好。本文的研究成果對(duì)魚雷渦輪機(jī)工程研制具有一定的參考價(jià)值。

魚雷渦輪機(jī); ANSYS CFX; 數(shù)值仿真

0 引言

隨著計(jì)算機(jī)技術(shù)的飛速發(fā)展和計(jì)算方法的不斷完善, 數(shù)值仿真作為一種獲得魚雷渦輪機(jī)內(nèi)部復(fù)雜流動(dòng)細(xì)節(jié)的手段體現(xiàn)出其有效性和優(yōu)越性, 為研究者們分析渦輪機(jī)內(nèi)部精細(xì)流場結(jié)構(gòu)提供了一個(gè)強(qiáng)有力的工具。

目前, 在大多數(shù)研究機(jī)構(gòu)中, 已把通過計(jì)算流體力學(xué)(computational fluid dynamics, CFD)技術(shù)來優(yōu)化渦輪機(jī)設(shè)計(jì)作為一個(gè)必不可少的手段, CFD已經(jīng)進(jìn)入燃?xì)鉁u輪機(jī)氣動(dòng)熱力學(xué)每個(gè)領(lǐng)域, 在流場分析、氣動(dòng)熱力學(xué)設(shè)計(jì)、性能分析、發(fā)現(xiàn)和驗(yàn)證新的流動(dòng)現(xiàn)象等方面均發(fā)揮著重要作用。使用經(jīng)過技術(shù)驗(yàn)證的全3D粘性流場數(shù)值計(jì)算軟件, 不僅推動(dòng)了燃?xì)鉁u輪機(jī)氣動(dòng)熱力學(xué)自身的迅速發(fā)展, 而且促進(jìn)了燃?xì)鉁u輪機(jī)研制的革命性變化, 使其逐漸擺脫了耗資多、周期長、風(fēng)險(xiǎn)大、主要依靠完備實(shí)驗(yàn)數(shù)據(jù)庫的傳統(tǒng)設(shè)計(jì)方法, 開始向“預(yù)測設(shè)計(jì)”過渡, 從而使先進(jìn)渦輪機(jī)的研制周期大大縮短[1]。

本文利用CFD軟件ANSYS CFX作為平臺(tái), 主要介紹了CFD技術(shù)應(yīng)用到魚雷渦輪機(jī)設(shè)計(jì)中的一般步驟。

1 渦輪機(jī)通流流場及性能數(shù)值方法

ANSYS CFX是解決流體動(dòng)力學(xué)分析的專業(yè)軟件, 由前處理、求解器和后處理3個(gè)子模塊組成, 具有以下特點(diǎn): 1) 提供了比較精確的數(shù)值方法; 2)提供了快速穩(wěn)定的求解技術(shù); 3) 提供了先進(jìn)的湍流模型等[2]。

1.1 數(shù)值方法計(jì)算流程

如圖1所示, CFD工作流程分為以下幾個(gè)步驟: 1) 利用CAD軟件(如CATIA、UG 等)建立幾何模型; 2) 導(dǎo)入網(wǎng)格劃分軟件ICEM CFD, 生成計(jì)算網(wǎng)格; 3) 將生成的計(jì)算網(wǎng)格導(dǎo)入Ansys CFX-pre模塊, 設(shè)置邊界條件, 如速度、壓力、溫度、流量等流體進(jìn)口和出口, 輸入工質(zhì)參數(shù), 選擇數(shù)值仿真參數(shù), 如穩(wěn)態(tài)或瞬態(tài), 可壓縮流體或不可壓縮流體; 4) 進(jìn)入Ansys CFX-slover模塊, 開始求解運(yùn)算, 并動(dòng)態(tài)監(jiān)視迭代殘差; 5) 收斂后導(dǎo)入后處理模塊, 進(jìn)行結(jié)果處理, 評(píng)估設(shè)計(jì)方案。

圖1 CFD工作流程

1.2 建立渦輪機(jī)通流部分CFD仿真模型

仿真之前必須要定義所要解決的問題, 即仿真的對(duì)象。仿真對(duì)象要涉及到具體幾何結(jié)構(gòu), 即所仿真的計(jì)算區(qū)域。

算例渦輪機(jī)為部分進(jìn)氣、沖動(dòng)式、軸流式渦輪機(jī)。研究對(duì)象是魚雷渦輪機(jī)通流部分, 其結(jié)構(gòu)簡圖見圖2, 主要包括噴管、葉輪等[3-4]。

1.3 網(wǎng)格劃分

網(wǎng)格劃分采用ICEM CFD處理器。ICEM CFD是一款世界頂級(jí)的CFD/CAE前處理器, 為各流行的CFD/CAE軟件提供高效可靠的分析模型。ICEM CFD 提供了廣泛的CAD 軟件接口, 當(dāng)幾何實(shí)體模型建好以后, 可以通過ICEM CFD 中的File/Import Geometry 直接導(dǎo)入CATIA V4, STEP, IGES, UG, Pro/E 等格式的CAD 幾何模型文件。ICEM CFD具有豐富的網(wǎng)格類型, 如四面體網(wǎng)格(Tetra Meshing)、三棱柱網(wǎng)格(Prism Meshing)、六面體網(wǎng)格(Hexagon Meshing)、棱柱形網(wǎng)格(Pyramid Meshing)、O型網(wǎng)格(O-grid Meshing)、自動(dòng)六面體網(wǎng)格(Auto Hexa)等。ICEM CFD 本身附帶了多種有限元求解軟件的輸入輸出轉(zhuǎn)換器。用戶可以通過Output/Select Solver 選擇求解器并通過Output/Write Input 來輸出文件[2]。本文選擇的求解器是ANSYS CFX, 輸出文件格式默認(rèn)為*.CFX5。處理器輸入和輸出流程如圖3所示。

圖2 渦輪機(jī)通流部分示意簡圖

圖3 ICEM CFD輸入與輸出

在對(duì)魚雷燃?xì)鉁u輪機(jī)通流部分計(jì)算域劃分網(wǎng)格時(shí), 需考慮葉輪徑向間隙和輪盤前后間隙。針對(duì)渦輪機(jī)通流部分結(jié)構(gòu)形式, 計(jì)算網(wǎng)格采用H-O-H多塊結(jié)構(gòu)化網(wǎng)格, 葉片壁面及輪轂壁面網(wǎng)格局部加密以捕捉邊界層, 葉尖間隙、葉輪與噴管環(huán)軸向間隙也適當(dāng)加密, 算例的網(wǎng)格點(diǎn)數(shù)約為106。由于渦輪機(jī)采用部分進(jìn)氣方式, 造成葉輪進(jìn)口流場不均勻, 無法采用周期性邊界條件, 需對(duì)渦輪機(jī)進(jìn)行全周數(shù)值仿真。渦輪機(jī)通流部分計(jì)算域網(wǎng)格示意圖見圖4。

圖4 渦輪機(jī)通流部分計(jì)算網(wǎng)格

1.4 網(wǎng)格導(dǎo)入ANSYS CFX

把劃分好的計(jì)算域網(wǎng)格導(dǎo)入到ANSYS CFX -pre; 可導(dǎo)入不同格式的網(wǎng)格, 如CFX Mesh(gtm), Ansys 及Nastran 等, 從ICEM CFD 中導(dǎo)入網(wǎng)格直接選擇ICEM CFD 網(wǎng)格格式。

1.5 設(shè)置工質(zhì)流體物性參數(shù)

CFX中包括了理想氣體、水和金屬等常用流體和固體的物理參數(shù), 也可設(shè)置新的流體的密度、粘度、導(dǎo)熱系數(shù)和比熱容等。

本文計(jì)算時(shí)假定通流部分燃?xì)饬鲃?dòng)是定常流動(dòng), 燃?xì)鉃槔硐霘怏w, 遵守理想氣體定律, 給定燃?xì)夥肿恿?、定壓比熱? 動(dòng)力粘度隨溫度的變化遵從Sutherland方程。

1.6 噴管盒與葉輪交界面處理

由于噴管盒固定不動(dòng), 而葉輪高速旋轉(zhuǎn), 噴管盒和葉輪部分必須分開建模, 而在計(jì)算時(shí)噴管盒計(jì)算域和葉輪計(jì)算域使用交接面相互連接, 如圖5所示。交接面采用多參考系模型和凍結(jié)轉(zhuǎn)子模型來建立噴管盒計(jì)算域和葉輪計(jì)算域之間的連接, 并實(shí)現(xiàn)數(shù)據(jù)相互傳遞。

圖5 噴管盒與葉輪交界面

1.7 選擇湍流模型

基于雷諾平均N-S方程對(duì)粘性流場的數(shù)組計(jì)算的一個(gè)核心問題是湍流模型的選擇和使用。對(duì)于工程流動(dòng)問題, 目前已經(jīng)發(fā)展了多種具有相當(dāng)經(jīng)驗(yàn)的湍流模型[5-7]。由于魚雷渦輪機(jī)通流在變工況下會(huì)出現(xiàn)激波/邊界層的相互干涉等跨音速流動(dòng)現(xiàn)象, 鑒于SST湍流模型具有較好的模擬強(qiáng)逆壓力梯度的能力, 在跨聲速/邊界層干涉流動(dòng)中取得了很大的成功[8], 本文研究選擇SST湍流模型。

1.8 邊界條件

所謂邊界條件, 是指求解域邊界上所求解的變量或其1階導(dǎo)數(shù)隨地點(diǎn)及時(shí)間變化的規(guī)律。在CFX中邊界條件包括: 流動(dòng)進(jìn)出口邊界、給定壓力邊界、壁面邊界、對(duì)稱邊界和周期性邊界等。不同的求解模型, 邊界條件的選擇也不同。

針對(duì)魚雷渦輪機(jī)通流部分流動(dòng)特點(diǎn), 基于對(duì)當(dāng)?shù)卣挥谶吔绲?D流動(dòng)分析, 對(duì)亞音速進(jìn)口流動(dòng), 3個(gè)進(jìn)口氣流參數(shù)需要給定。本文進(jìn)口邊界采用壓力進(jìn)口, 給定氣流的總壓、總溫和氣流方向, 其他的氣流參數(shù)由流場外插可得; 出口采用壓力出口, 給定靜壓, 其他的氣流參數(shù)由流場外插可得; 對(duì)于湍流模型進(jìn)口給定中等湍流度和旋渦粘性比; 假定通流部分壁面是絕熱、光滑壁面。

1.9 設(shè)置求解參數(shù)

通過Solver Control 設(shè)置對(duì)流項(xiàng)離散格式和時(shí)間步長以及收斂精度。通過Write Solver File進(jìn)行求解計(jì)算。

本文求解的是3D定常粘性的雷諾平均NAVIER-STOKES方程組。控制方程的離散采用基于有限元的有限體積法, 耦合隱式格式的時(shí)間推進(jìn)算法, 對(duì)控制方程對(duì)流項(xiàng)的離散采用2階迎風(fēng)格式。

1.10 收斂標(biāo)準(zhǔn)

計(jì)算實(shí)時(shí)監(jiān)測NAVIER~STOKES方程的迭代殘差、通流部分進(jìn)出口流量以及葉輪扭矩。收斂標(biāo)準(zhǔn): 1) 流量方程、動(dòng)量方程、湍流方程的迭代殘差要小于10-4或出現(xiàn)震蕩; 2) 計(jì)算進(jìn)出口流量近似平衡, 誤差不超過0.5%; 3) 計(jì)算葉輪扭矩基本不隨計(jì)算時(shí)間變化。

只要滿足上述收斂標(biāo)準(zhǔn), 如圖6所示, 再計(jì)算100步, 可視為計(jì)算收斂。

圖6 計(jì)算監(jiān)測窗口

1.11 后處理及結(jié)果評(píng)價(jià)

ANSYS CFX提供了強(qiáng)大的后處理功能, 把結(jié)果. res文件導(dǎo)入到CFX-POST, 可對(duì)結(jié)果進(jìn)行后處理。

計(jì)算完成以后, CFX會(huì)自動(dòng)生成一個(gè).out文件和.res文件。其中.out文件中包含求解信息, 如模型設(shè)定、求解過程中解得狀態(tài)以及求解工作量統(tǒng)計(jì)等, 在其中可把渦輪機(jī)進(jìn)出口流量、葉片負(fù)載等性能參數(shù)摘取出來, 也可把結(jié)果.res文件導(dǎo)入CFX-POST中進(jìn)行后處理, 利用Report功能得出相應(yīng)性能參數(shù)。

2 算例結(jié)果分析

以魚雷燃?xì)鉁u輪機(jī)通流部分為研究對(duì)象, 利用上面介紹的數(shù)值計(jì)算方法, 開展通流流場及性能的數(shù)值計(jì)算分析。

2.1 計(jì)算方法試驗(yàn)驗(yàn)證

應(yīng)用魚雷渦輪機(jī)熱吹風(fēng)試驗(yàn)測量數(shù)據(jù)對(duì)計(jì)算方法進(jìn)行驗(yàn)證。計(jì)算結(jié)果與試驗(yàn)數(shù)據(jù)比較情況見表1。

表1 相對(duì)有效效率和流量的計(jì)算結(jié)果與試驗(yàn)數(shù)據(jù)比較

Table 1 Comparison of CFD results and test data of re- lative effective efficiency and flow rate

由表1可知, 數(shù)值計(jì)算的渦輪機(jī)性能與試驗(yàn)數(shù)據(jù)吻合較好, 最大偏差不超過6%, 這表明應(yīng)用計(jì)算流體動(dòng)力學(xué)軟件ANSYS CFX所構(gòu)建的魚雷燃?xì)鉁u輪機(jī)通流部分流場及性能數(shù)值計(jì)算模型是可信的。

2.2 噴管內(nèi)部流場細(xì)節(jié)

圖7分別給出了噴管內(nèi)部燃?xì)怦R赫數(shù)、靜壓云圖分布及沿流向的變化趨勢??芍? 氣體在噴管中流動(dòng)時(shí), 由于發(fā)生膨脹, 其壓力很快降低, 而絕對(duì)速度急劇增加。馬赫數(shù)由進(jìn)口的0.06增大到2.23, 在喉部位置馬赫數(shù)為1.0。靜壓由入口2.966 MPa降低到0.21 MPa, 喉部為1.74 MPa。噴管處于超臨界狀態(tài), 沒有激波。只有氣流與壁面的摩擦損失[9]。

圖7 噴管內(nèi)部燃?xì)鈪?shù)沿流向變化曲線

2.3 葉輪內(nèi)部流場細(xì)節(jié)

圖8給出了噴管環(huán)所對(duì)應(yīng)的葉輪內(nèi)部氣流參數(shù)沿流向的變化以及葉中位置相對(duì)馬赫數(shù)云圖分布[10]??傮w上來講, 葉輪入口燃?xì)庀鄬?duì)速度要比出口大。葉輪入口燃?xì)庀鄬?duì)馬赫數(shù)超音, 燃?xì)馀龅饺~片前緣, 分開向吸力面和壓力面流動(dòng), 在吸力面氣流繞過前緣以后先有一段加速過程。相對(duì)于燃?xì)鈬娙肴~輪方向, 流道入口為一斜切入口, 是一個(gè)收縮的通道, 超音速燃?xì)庠谶@里受到阻滯, 產(chǎn)生正激波, 激波以后燃?xì)馑俣燃眲∠陆? 變?yōu)閬喴羲賉11]。正激波打在相鄰葉片的吸力面上, 形成激波邊界層干涉, 引起吸力面邊界層增厚, 嚴(yán)重時(shí)會(huì)產(chǎn)生邊界層分離[12-13]。通道中間處燃?xì)饩S持亞音速流動(dòng)。葉輪通道出口也為一斜切口, 相對(duì)于燃?xì)鈦碇v是擴(kuò)張型, 亞音速燃?xì)庥珠_始膨脹,產(chǎn)生膨脹波, 燃?xì)庀鄬?duì)速度增大, 在葉片尾緣處燃?xì)庀鄬?duì)馬赫數(shù)又有超音; 此時(shí)與葉片吸力面上的低速邊界層氣流相遇, 相當(dāng)于產(chǎn)生內(nèi)折, 在葉片尾緣處產(chǎn)生斜激波, 成燕尾波, 后燃?xì)庵匦伦優(yōu)閬喴羲賉12-13]。

圖8 葉輪內(nèi)部流場細(xì)節(jié)

3 結(jié)束語

本文應(yīng)用ANSYS CFX軟件對(duì)魚雷燃?xì)鉁u輪機(jī)通流部分流場及性能進(jìn)行建模和數(shù)值仿真, 介紹了計(jì)算流體動(dòng)力學(xué)技術(shù)應(yīng)用到魚雷燃?xì)鉁u輪機(jī)設(shè)計(jì)中的一般步驟。研究表明, 所構(gòu)建的魚雷燃?xì)鉁u輪機(jī)通流流場及性能計(jì)算模型是可信的, 預(yù)測結(jié)果與試驗(yàn)數(shù)據(jù)吻合較好。本文的研究成果對(duì)魚雷渦輪機(jī)工程研制具有較強(qiáng)的參考價(jià)值。

[1] 劉大響. 對(duì)加快發(fā)展我國航空動(dòng)力的思考[J]. 航空動(dòng)力學(xué)報(bào), 2001, 16(1): 1-7. Liu Da-xiang. Deliberation Upon Advance of Aeroengine Development in China[J]. Journal of Aerospace Power, 2001, 16(1): 1-7.

[2] ANSYS 公司. ANSYS CFX產(chǎn)品簡介[EB/OL]. http://www. ansyscom.cn.

[3] 趙寅生, 錢志博. 魚雷渦輪機(jī)原理[M]. 西安: 西北工業(yè)大學(xué)出版社, 2002.

[4] 查志武, 史小峰, 錢志博. 魚雷熱動(dòng)力技術(shù)[M]. 北京: 國防工業(yè)出版社, 2006.

[5] Wilcox D C. Simulation of transition with a Two- Equation turbulence model[J]. AIAA Journal, 1994, 32(2): 247-253.

[6] Menter F R. Two-Equation Eddy-Viscosity Turbulence Models for Engineering Applications[J]. AIAA Journal, 1994, 32(8): 1598-1605.

[7] Langtry R B, Menter F R. Transition Modeling for General CFD Applications in Aeronautics[R]. AIAA Paper 2005-522,Reno, Nevada.

[8] Thivet F. Lessons learned from RANS simulations of shock- wave/boundary-layer interactions[R]. AIAA paper 2002-0583.

[9] 伊進(jìn)寶, 趙衛(wèi)兵, 師海潮. 魚雷渦輪機(jī)斜切噴管內(nèi)流場數(shù)值模擬[J]. 魚雷技術(shù), 2010, 18(3): 223-227. Yi Jin-bao, Zhao Wei-bing, Shi Hai-chao. Numerical Simulation of Flow Field in Scarfed Nozzle of Torpedo Turbine[J]. Torpedo Technology, 2010, 18(3): 223-227.

[10] 伊進(jìn)寶, 趙衛(wèi)兵, 師海潮. 部分進(jìn)氣燃?xì)鉁u輪機(jī)葉輪流場數(shù)值模擬[J]. 魚雷技術(shù), 2010, 18(6): 456-460. Yi Jin-bao, Zhao Wei-bing, Shi Hai-chao. Numerical Simulation on Inner Flow Field of Gas Turbine with Partial Inlet Flow[J]. Torpedo Technology, 2010, 18(6): 456-460.

[11] 潘錦珊. 氣體動(dòng)力學(xué)基礎(chǔ)[M]. 西安: 西北工業(yè)大學(xué)出版社, 1995.

[12] 王保國, 劉淑艷, 黃偉光. 氣體動(dòng)力學(xué)[M]. 北京: 北京理工大學(xué)出版社, 2006.

[13] 伍貽兆, 楊岞生. 跨聲速空氣動(dòng)力學(xué)[M]. 北京: 國防工業(yè)出版社, 2004.

Numerical Simulation of Torpedo Turbine Flow Passage Field and Performance Using ANSYS CFX

CHEN Gang1, YI Jin-bao2, SHI Hai-chao3

(1. Xi′an Representative Bureau, Naval Armament Department, Xi′an 710054, China; 2. The 705 Research Institute, China Shipbuilding Industry Corporation, Xi′an 710075, China)

Torpedo turbine flow passage field and performance are simulated by using computational fluid dynamics(CFD) software ANSYS CFX, and the general steps of designing torpedo gas turbine flow passage with CFD are introduced in order to predict the gas turbine flow passage field and performance and to promote application of CFD technique to design of torpedo turbines. The results show that the predicted turbine efficiency agrees well with the test data. This study may offer a reference for development of new torpedo turbines.

torpedo turbine; ANSYS CFX; numerical simulation

TJ630.32

A

1673-1948(2012)04-0285-05

2011-03-06;

2011-11-07.

陳 剛(1974-), 男, 工程師, 主要從事魚雷總體及質(zhì)量方面的工作.

(責(zé)任編輯: 陳 曦)

猜你喜歡
模型
一半模型
一種去中心化的域名服務(wù)本地化模型
適用于BDS-3 PPP的隨機(jī)模型
提煉模型 突破難點(diǎn)
函數(shù)模型及應(yīng)用
p150Glued在帕金森病模型中的表達(dá)及分布
函數(shù)模型及應(yīng)用
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計(jì)的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 波多野结衣久久精品| 伊人久热这里只有精品视频99| 在线国产三级| 一区二区三区成人| 夜夜高潮夜夜爽国产伦精品| av在线手机播放| 久久综合丝袜长腿丝袜| 成人午夜福利视频| 欧美成人影院亚洲综合图| 亚洲无码免费黄色网址| 亚洲无码在线午夜电影| 亚洲欧洲日产国码无码av喷潮| 精品视频在线一区| 国内视频精品| 在线看国产精品| 99无码中文字幕视频| 一级毛片免费不卡在线| 国内精自线i品一区202| 五月激情综合网| 欧洲免费精品视频在线| 午夜精品福利影院| 亚洲第一视频免费在线| 午夜在线不卡| 亚洲无线视频| 中文字幕在线播放不卡| 日本一本正道综合久久dvd| 玖玖精品视频在线观看| 国产精品久久久久久久久kt| 无码视频国产精品一区二区| 综合天天色| 国产新AV天堂| 国产成人啪视频一区二区三区| 91久久国产成人免费观看| 97在线视频免费观看| 亚洲女同一区二区| 亚洲AⅤ无码日韩AV无码网站| 中文纯内无码H| 久久精品女人天堂aaa| 国禁国产you女视频网站| 国产正在播放| 久久久久青草大香线综合精品| 欧美激情网址| 国产高潮视频在线观看| 日韩精品免费一线在线观看| 欧美全免费aaaaaa特黄在线| 欧美a在线| 日本人妻一区二区三区不卡影院| 97人妻精品专区久久久久| 亚洲天堂网在线播放| 国产一区二区网站| 午夜日b视频| 亚洲无码91视频| 婷婷亚洲最大| 欧美国产精品不卡在线观看 | 在线观看热码亚洲av每日更新| 露脸一二三区国语对白| 国产成人高清精品免费软件 | 亚洲国产成人综合精品2020 | 四虎永久在线视频| 91亚瑟视频| 日韩精品欧美国产在线| 四虎国产精品永久在线网址| 国产一级片网址| 91麻豆精品国产高清在线| 日本欧美一二三区色视频| 色哟哟色院91精品网站| a色毛片免费视频| 欧美www在线观看| 国产高潮视频在线观看| 黄色a一级视频| 伊人久久综在合线亚洲2019| 18禁不卡免费网站| 国产尤物在线播放| 全部免费特黄特色大片视频| 黄片在线永久| 五月天天天色| 日本午夜影院| 欧美视频免费一区二区三区| 国产精品99久久久久久董美香| 97久久精品人人| 国产午夜人做人免费视频中文| 国产欧美日本在线观看|