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

某柴油發動機缸內燃燒的數值模擬及優化分析

2017-07-06 11:02:55賴晨光陳永燕段孟華周毓婷
關鍵詞:優化模型

賴晨光,陳永燕,王 媛,段孟華,周毓婷

(1.重慶理工大學 a.車輛工程學院; b.化學化工學院, 重慶 400054;2.日本東北大學 流體科學研究所, 仙臺 980-8577)

?

某柴油發動機缸內燃燒的數值模擬及優化分析

賴晨光1a,2,陳永燕1a,王 媛1a,段孟華2,周毓婷1b

(1.重慶理工大學 a.車輛工程學院; b.化學化工學院, 重慶 400054;2.日本東北大學 流體科學研究所, 仙臺 980-8577)

采用基于克里精(Kriging)代理模型的混合多目標梯度優化算法(HMGE),通過FLUENT軟件平臺,對柴油發動機燃燒過程進行了多目標優化,優化變量為縮口率、余隙高度、凸臺高度、噴油錐角,優化目標為平均溫度、平均壓力、碳煙(Soot)排放、氮氧化合物(NOx)排放。完成多目標優化計算后,應用數據挖掘的總變差分析方法(ANOVA)和自組織映射分析方法(SOM)對優化變量和優化目標函數進行定性和定量分析,揭示出優化變量與優化目標之間的相互關系。分析結果表明:較好的動力性能需要較小的縮口率和較小的余隙高度;較少的碳煙排放需要較大的余隙高度,較低的氮氧化合物排放需要較小的凸臺高度和較小的縮口率。

柴油發動機;燃燒;代理模型;優化算法

柴油發動機缸內氣流流動和燃燒過程對整機的性能有很大的影響,并且也是減少有害排放、提高經濟性和動力性的決定性因素[1]。燃燒室形狀以及某些噴油參數直接影響柴油機進氣流動、混合氣體的形成和燃燒[2],因此近年來對于柴油機燃燒室形狀的優化研究也越來越受重視。傳統研究大多以缸內壓力、溫度以及排放的Soot、NOx含量作為優化對象,采用傳統優化方法對優化目標進行逐一優化。但是,傳統的優化方法優化目標單一、計算量大,并且變量參數的選擇是人為提取,提取的參數點空間分布不均勻,不具有代表性,從而影響計算結果的準確性[3]。

本文首先對柴油發動機缸內的工作過程進行數值模擬分析,包括燃燒室缸內的流場特性、碳煙排放、NOx排放等。再對燃燒室形狀進行多目標多變量優化,并對優化結果進行了數值模擬驗證。

1 數值模擬

1.1 幾何模型的建立

采用三維建模軟件CATIA建立某單缸四沖程柴油發動機燃燒室幾何模型,其主要技術性能參數見表1。

1.2 網格劃分

網格的劃分是進行數值模擬的基礎,生成網格的好壞直接影響計算結果的準確性。本研究網格劃分通過ANSYS軟件的IC Engine模塊中的Meshing完成。IC Engine模塊是安世亞太公司在ANSYS軟件中新增的一個專門用于進行發動機數值模擬研究的模塊。該模塊集成有Design Model、Meshing、IC Engine Solver、Fluent等軟件,將建模、網格劃分、邊界條件賦值、求解計算等過程直接關聯起來,操作方便,并且由于該模塊主要就是針對發動機數值模擬,所以劃分出的網格質量很高。由于本研究的柴油發動機采用4孔噴油器且孔分布均勻,為節約數值模擬計算時間,采用1/4 的模型進行網格劃分。當發動機工作時,網格數量會隨著曲軸轉角的變化不斷增加或者減少,上止點網格數量為92萬,下止點網格數量為226萬。

表1 柴油機基本技術參數

1.3 初始條件與邊界條件

在通用流體計算軟件FLUENT中,對發動機缸內燃燒的數值模擬是以曲軸轉角為單位進行計算的。本次數值模擬的發動機工況如下:轉速為2 200 r/min,負荷為100%,計算范圍從進氣門關閉(563°CA)到排氣門開啟(856°CA),噴油時刻為710°~720°CA。計算開始時假設燃燒室內的流場均勻單一,缸內初始溫度(T0)為420 K,初始壓力(P0)為0.17 MPa,渦流比為1.6,每個噴孔每循環噴油量為13.356 mg,燃油噴射的溫度為353 K。溫度采用恒溫邊界條件,氣缸蓋底面溫度為520 K,活塞頂面溫度為560 K,氣缸壁面溫度為460 K。進氣門關閉時燃燒室流場的湍動能(TKE)和湍動能耗散率(TLS)根據以下公式計算。

TKE=(3/2)×u2

(1)

u=1.4×h×(n/60)

(2)

TLS=hv/2

(3)

式中:h為沖程長度(m);n為轉速(r/min);hv為氣門最大升程(mm);u為湍流脈動速度(m/s)。

邊界條件中湍流模型選擇標準k-ε雙方程模型;噴霧模型選擇WAVE模型;燃燒模型選擇渦團破碎模型;排放物模型選擇Zeldovich NO模型和Moss-Brookes模型;離散方法選擇有限體積法;流場計算方法選擇SIMPLE算法;計算時應用了動網格技術。

2 計算結果分析

2.1 缸內壓力及溫度

圖1是燃燒室內的平均壓力與平均溫度曲線。燃燒室內最大壓力發生在上止點后3°CA,其值為14.6 MPa;燃燒室內最高溫度發生上止點后7° CA,其最高溫度為2 897 K。從平均溫度曲線的發展趨勢來看:在上止點前10° CA 的地方有一段平緩區域,在上止點前8° CA 的地方迅速上升,與噴油提前角為10° CA吻合,從爆發壓力時到達到最高溫度時(即從上止點后3° CA 到7° CA)為主燃期,發生在噴油結束后。

2.2 排放物濃度

圖2為柴油發動機燃燒室內碳煙和氮氧化合物含量的變化趨勢。碳煙生成的基本條件為缺氧,其溫度要求比氮氧化合物的低,因此使得碳煙的生成時刻比氮氧化合物提前。在上止點前2° CA左右碳煙開始生成,并在上止點后25° CA(B點)時達到了最大值,隨后燃燒室內的高溫以及殘余氧氣作用使生成的碳煙隨后被氧化[4],其含量漸漸降低。氮氧化合物生成的基本條件為高溫、富氧和較長的反應時間。在滯燃期階段氮氧化合物生成量基本上為0,這是由于燃燒室的溫度較低,不符合氮氧化合物生成條件。當缸內燃油進入主燃階段時,溫度不斷升高,氮氧化合物的生成量相應地增加,在上止點后10° CA(A點)處達到最大,之后由于缸內溫度和含氧量的下降NOx含量保持不變。

圖1 缸內壓力、溫度曲線

圖2 Soot、NOx質量分數

3 燃燒的優化

3.1 優化模型

本次優化的變量是縮口直徑、余隙高度、凸臺高度、油束錐角。表2為柴油發動機燃燒室形狀優化空間中設計變量的取值范圍。

表2 設計變量的取值范圍

優化目標是得到缸內平均壓力盡可能高、平均溫度盡可能大、氮氧化合物和碳煙含量盡可能少的關鍵參數組合,如表3所示。

表3 優化目標

3.2 優化方法

本文通過CFD軟件平臺、采用混合多目標梯度優化算法對4個優化變量進行多目標優化。但是,柴油機燃燒過程的數值模擬是一個很耗時的計算,在本次優化中,完成1組數據的數值模擬用1臺48核的工作站就需要計算18 h,所以如果僅采用數值模擬和優化算法,會使得計算任務太過繁重,占用資源較多。為了解決這一難題,便在數值計算和優化算法中間增加了建立代理模型這一過程。本文選擇的是Kriging模型。

Kriging模型是由一個參數模型和非參數隨機聯合構成的,計算時不需要建立某個特定的數學模型,只需要通過部分已知的信息就可以去擬合某一點的未知信息,比單個參數化模型具有更強的預測能力和靈活性[5]。在優化時采用Kriging代理模型可以通過已知樣本點對未知樣本點的值進行預測,并且可對預測值進行初步誤差計算。這樣就可以不用對每一組樣本數據都進行數值模擬,從而大幅降低了計算量。

混合多目標梯度優化算法(hybrid multi-gradient explorer,HMGE)是一種基于遺傳算法和梯度算法的優化算法,它結合了遺傳算法全局性好的優點,也保留了梯度算法有效性[6]。該算法優化途徑:首先,采用遺傳算法獲得非支配解集,然后采用梯度搜索法尋找非劣解,即搜索時隨機選取目標函數的梯度,分別按其正方向和負方向搜索,得到2個子個體,合并成一個大種群,再用Pareto前沿的判斷方法尋找該種群中的非劣解[7];將采用遺傳算法得到的非支配解集與通過梯度算法得到的非劣解集進行對比,尋找精確收斂到局部的Pareto解集,通常10~20次迭代即可獲得一個Pareto 解。該算法與傳統的優化算法相比,擁有計算精度高、效率高、全局求解能力強等優點。

3.3 優化流程

圖3為本次優化計算的流程。本文首先用拉丁超立方的取樣方法,從優化變量空間里提取40個樣本點;然后通過數值模擬計算出這40個初始點對應的缸內平均溫度、平均壓力、NOx平均質量分數、Soot平均質量分數的值;接著通過40個初始點及其求解值建立Kriging代理模型,并基于建立的Krigine代理模型采用混合多目標梯度優化算法搜索全局最優解;完成算法尋優后,將尋找到的最優解集通過k-mean聚類的方法聚為4類,選取每類的中心點作為優化結果驗證的樣本點;通過數值模擬的手段對選取的最優解集中的樣本點進行誤差驗證,若誤差值在10%以內,則建立的Kriging模型是有效的,若誤差值超過了10%,則需要通過尋找EI最大值、添加初始樣本的個數來提高代理模型的精度,再進行尋優計算。

表4為此次運用k-mean聚類方法選取的4個樣本點的誤差驗證結果,最大誤差均在Kriging代理模型的精度允許誤差10%以內,所以本次建立的Kriging模型是有效的。

圖3 優化流程

目標函數平均溫度平均壓力NOx含量Soot含量最大誤差6.5%8.2%0.49%4.4%

4 數據挖掘

數據挖掘(data mining)[8-14]是一種數據分析的方法,它可以從擁有海量的、隨機的、模糊的、殘缺的數據中提取潛在的、有用的信息,并且就數據分析的本質來說,它不僅能進行定性分析還能定量分析數據之間的相互關系。本文就是應用數據挖掘方法中的總變差分析方法(analysis of variance,ANOVA)和自組織映射分析方法(self-organization mapping,SOM)來探索4個優化變量對4個目標函數的影響規律,為以后的優化研究提供一定的參考。

4.1 總變差分析方法

總變差分析方法是一種統計學的定量分析方法,可以用來揭示設計變量對設計目標的影響規律。

如圖3所示,圖中相應區域的百分數的大小表示相應設計變量對目標的影響大小。對缸內平均壓力影響較大的是縮口率與余隙高度,其比重值之和達到了73.8%。平均溫度受縮口率、余隙高度、凸臺高度的影響均較大,其中:平均溫度受縮口率影響最大,其所占比重為32.3%。對缸內碳煙含量影響較大的設計變量為縮口率和余隙高度,其所占比例分別為 43.9%和29.9%。氮氧化合物主要受到余隙高度和凸臺高度的影響,所占比例之和高達79%。噴油錐角的影響最小。

圖4 總變差分析結果

4.2 自組織映射分析

自組織映射是一種降維且能保留原始數據結構特征的研究方法,它將任意維的輸入信號模式轉變為二維的離散映射,然后獲取設計變量與響應變量的二維神經元網絡分布,從神經元網格上可以定性分析出設計變量與響應變量之間的內在關系,以及設計變量之間交互影響關系及與響應變量的關系[15]。形成神經元網絡的原理:輸出層上某一結點能對某一模式作出特別反應來代表該模式類,當某類數據模式輸入時,會對輸出層某一結點產生最大刺激,同時也給周圍結點帶來刺激,產生最大刺激的點成為獲勝結點,每次的訓練都會使獲勝結點及其鄰域結點的連接權值得到調整,如此反復,直至連接權值調整微小為止。

圖5為設計變量和響應變量的神經元網格,本次使用了用于訓練Kriging代理模型的40個初始樣本來訓練神經元。8張神經元網格圖都是來源于同一張自組織映射網格,圖中的顏色代表了該設計變量值或響應變量值的大小,越偏向藍色代表值越小,越偏向紅色代表值越大。

圖5 設計變量和響應變量神經元網格圖

4.2.1 設計變量與設計目標之間關系分析

從總差變分析分析已知:對缸內平均溫度影響最大的是燃燒室的縮口率,其次是余隙高度。對比缸內平均溫度神經元網格的顏色模式和燃燒室縮口率神經元網格及余隙高度神經元網格的顏色模式,發現縮口率與平均溫度沒有完全相同或相反的趨勢,但是只有縮口率取藍色區域的某些值才能使平均溫度達到最大,縮口率取紅色區域的某些值才能使平均溫度達到最小;而余隙高度與平均溫度神經元網格的顏色模式正好相反,平均溫度與余隙高度是負相關的關系。

同上分析,平均壓力和碳煙的趨勢一樣,它們均與縮口率成非線性關系,但是縮口率大的時候平均壓力和碳煙含量均較小,平均壓力和碳煙含量的最大值都出現在縮口率小的某一區域;另外,兩者皆與余隙高度成負相關的關系。氮氧化合物(NOx)排放量與燃燒室凸臺高度神經元網格呈現出與之相似的趨勢,這說明氮氧化合物的排放量與凸臺高度為正相關關系;NOx排放與燃燒室縮口率、凸臺高度、油束夾角大小為非線性關系。

4.2.2 設計目標之間相互關系分析

對比圖4中的4張設計目標自組織映射神經元網格可以發現,平均溫度、平均壓力、碳煙排放的顏色變化趨勢大致相同,說明此三者是正相關的關系。但是對于優化要求,需要碳煙排放最少,平均溫度和平均壓力最大,所以碳煙排放與平均溫度和平均壓力不能同時達到最優,而最高溫度與平均壓力可以達到最優。缸內溫度與氮氧化合物排放關系沒有與碳煙排放的關系那么明顯,但是從圖5還是可以看出:高溫時會出現大的氮氧化合物排放。

對柴油發動機的整體性能而言,動力性和排放性在自組織映射神經元網格中的分布呈對立形勢,即這兩種性能在設計的時候不能同時達到最優,設計人員需選擇一種折中的方案進行設計。

5 優化結果驗證

通過K-mean聚類分析,選出4個讓至少一個目標達到最優且兼顧其他3個設計目標的最優點,然后進行數值模擬計算,將數值模擬的值與其預測值進行誤差分析。選出的4個點處的燃燒室模型如圖6所示。最優解A偏好于動力性好;最優解B偏好于碳煙排放最低;最優解C偏好于碳煙和氮氧化物排放都低;最優解D同時兼顧動力性能和排放性能。4個驗證點的優化預測結果與仿真結果的誤差如表5所示。由表5可知:其誤差均小于5%,在可以接受誤差范圍以內,說明該優化方法具有可行性。

表5 優化結果的預測值與仿真值的誤差 %

圖6 優化后驗證的燃燒室模型

通過數值模擬,對比優化前后得到的燃燒室模型與原模型對設計目標的影響,結果如表6所示。與原模型相比可見:最優解A在排放性能不變差的情況下,缸內平均壓力和平均溫度分別提升了6.28%和7.85%;最優解B在動力性和NOx排放不變差的條件下,碳煙含量減少了8.03%;最優解C在動力性不變壞的情況下,碳煙和氮氧化合物排放含量分別減少了5.62%和4.70%;最優解D在同時兼顧動力性能和排放性能的情況下,缸內平均壓力和平均溫度分別提升了4.52%和3.13%,缸內NOx含量和碳煙含量排放分別減少了3.85%和2.54%。

表6 優化后與原樣模型的數值模擬結果對比

6 結束語

本文以縮口率、余隙高度、凸臺高度和燃油噴射時的油束錐角作為設計變量,以平均溫度、平均壓力、碳煙排放含量、氮氧化物排放含量作為目標變量,采用了基于kriging代理模型的混合多目標梯度優化算法(HMGE算法)進行優化,并通過數據挖掘中的總變差和自組織映射(神經網絡法)對結果進行了分析,最后從最優解集中選擇4個點進行了數值模擬,并與原模型分析結果進行對比。結果表明:優化后的模型各方面性能均得到提高。

在采用數據挖掘的方法對設計變量與設計目標、設計目標與設計目標之間的關系進行分析時,得出以下規律:缸內平均溫度與余隙高度呈負相關的關系;平均壓力與凸臺高度、碳煙與余隙高度均呈負相關關系;氮氧化合物的排放量與凸臺高度呈正相關關系。從優化目標來說,碳煙排放和氮氧化合物排放可以同時得到優化。動力性能和排放性能在設計中是具有矛盾關系的,只能選擇折中方案進行設計。

多目標優化方法的應用極大縮短了計算時間,并且也使得優化取得了極好的效果。此次多目標優化研究得出的結果以及形成的一套可行的優化方法對設計人員的開發研究具有一定的指導意義。

[1] 王欣.4D24柴油機燃燒過程的多維數值模擬[D].南昌:南昌大學,2013.

[2] 胡林峰.柴油機對燃油噴射系統的要求和噴油系統的發展趨勢[J].現代車用動力,2002(4):1-4.

[3] DUAN Menghua,LAI Chenguang.Design space exploration on Combustion Chamber international[J].Conference on Flow Dynamics,2015(11):701-702.

[4] 袁方恩,林學東,田維,等.縮口燃燒室中氣流特性與燃油噴霧匹配對柴油機燃燒及排放的影響[J].吉林大學學報(工學版),2011(5):629-634.

[5] 鄒林君.基于Kriging模型的全局優化方法研究[D].武漢:華中科技大學,2011: 16-17.

[6] KAVEH R, NAGARAJAN G.Optimization of Diesel Engine Operating Parameters Using a Response Surface Method[J].SAE Technical Paper,2010(1):1262.

[7] 賴晨光,陳小雄.基于遺傳算法某汽車外形空氣動力學優化[J].重慶理工大學學報(自然科學),2016,30(4):1-5.[8] GUO Z,ZHOU Z,SONG L,et al.Aerodynamic Analysis and Multi-Objective Optimization Design of a High Pressure Ratio Centrifugal Impeller[C]//2014:V02DT42A013,ASME Turbo Expo 2014.2014.

[9] 賴晨光,陸茂桂,張海林.基于數據挖掘的皮卡貨車導流板減阻研究[J].重慶理工大學學報(自然科學),2017,31(2):1-6.

[10]黃解軍,潘和平,萬幼川.數據挖掘技術的應用研究[J].計算機工程與應用,2003,39(2):45-48.

[11]張春華,王陽.數據挖掘技術、應用及發展趨勢[J].現代情報,2003,23(4):47-48.

[12]郭慧東,王瑋,夏明超.基于數據挖掘的風電機組變槳系統劣化狀態在線辨識方法[J].中國電機工程學報,2016,36(9):2389-2397.

[13]李蔚,俞蕓蘿,盛德仁,等.基于動態數據挖掘的熱力參數傳感器故障診斷[J].振動、測試與診斷,2016,36(4):694-699.

[14]王文正,鄭鹍鵬,陳功,等.數據挖掘技術在飛行試驗數據分析和氣動參數辨識中的應用研究[J].空氣動力學學報,2016,34(6).

[15]SUMEET P,NADER F.Self Organizing Maps (SOM) for Design Selection in Multi-Objective Optimization Using modeFRONTIER[J].SAE,2008(1):1-10.

(責任編輯 劉 舸)

Numerical Simulation and Optimization Analysis of Combustion in a Diesel Engine

LAI Chen-guang1a,2, CHEN Yong-yan1a, WANG Yuan1a, DUAN Meng-hua2, ZHOU Yu-ting1b

(1.a.College of Vehicle Engineering; b.College of Chemistry and Chemical Engineering, Chongqing University of Technology, Chongqing 400054, China;2.Institute of Fluid Science, Tohoku University, Sendai 980-8577, Japan)

This research uses hybrid multi-objective Gradient exploration algorithm (HMGE) which is based on Kriging surrogate model to optimize and analysis combustion process and performance of diesel engine, and all of the numerical simulation are completed by FLUENT software. Caliber reducing rate, clearance height, convex platform height, injection cone angle are choose as optimization variables, average temperature, average pressure, soot and NOxemissions are the optimization objectives. After achieving multi-objective optimization, two data mining methods, analysis of variance(ANOVA) and self-organizing feature map(SOM) are used to qualitatively and quantitatively analyze the design variables and objective functions, and reveals the influence of design variables to design objectives. The results show that a good dynamic performance demands smaller convex platform height and caliber reducing rate; less soot emissions need larger clearance, and lower emissions of nitrogen oxides needs smaller convex platform height and smaller caliber reducing rate.

diesel engine; combustion; surrogate model; optimization analysis

2017-03-18

國家自然科學基金資助項目(51305477);重慶理工大學研究生創新基金資助項目(YCX2015204)

賴晨光(1978—),男,江西贛州人,博士,教授,主要從事汽車與高速列車空氣動力學研究,E-mail:Chenguanglai@cqut.edu.cn;通訊作者 陳永燕(1991-),女,重慶巫溪人,碩士研究生,主要從事汽車氣動特性和汽車發動機缸內流場的研究工作,E-mail:1453293871@qq.com。

賴晨光,陳永燕,王媛,等.某柴油發動機缸內燃燒的數值模擬及優化分析[J].重慶理工大學學報(自然科學),2017(6):23-30.

format:LAI Chen-guang,CHEN Yong-yan, WANG Yuan, et al.Numerical Simulation and Optimization Analysis of Combustion in a Diesel Engine [J].Journal of Chongqing University of Technology(Natural Science),2017(6):23-30.

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

U464.12+3

A

1674-8425(2017)06-0023-08

猜你喜歡
優化模型
一半模型
超限高層建筑結構設計與優化思考
房地產導刊(2022年5期)2022-06-01 06:20:14
民用建筑防煙排煙設計優化探討
關于優化消防安全告知承諾的一些思考
一道優化題的幾何解法
由“形”啟“數”優化運算——以2021年解析幾何高考題為例
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
主站蜘蛛池模板: 亚洲天堂视频在线观看免费| 久久国产av麻豆| 精品视频福利| 国内精品九九久久久精品| 中文成人在线| 久久情精品国产品免费| yjizz国产在线视频网| 成年女人a毛片免费视频| 综合成人国产| 亚洲精品中文字幕无乱码| 中文字幕亚洲另类天堂| 思思热在线视频精品| 国产精品专区第1页| 色成人亚洲| 色屁屁一区二区三区视频国产| 国产精品手机在线播放| 亚洲第一页在线观看| 最新精品国偷自产在线| a天堂视频在线| 免费毛片全部不收费的| 国产精品爽爽va在线无码观看| 国产乱子伦精品视频| 亚洲va在线∨a天堂va欧美va| 不卡国产视频第一页| 99视频在线免费看| 午夜国产精品视频| 国内精品小视频在线| 最近最新中文字幕在线第一页| 日本午夜在线视频| 欧美不卡视频一区发布| 中文字幕久久亚洲一区| 欧美一级在线看| 国产在线精品99一区不卡| 99这里只有精品在线| 伊人久久精品亚洲午夜| 国产亚洲视频免费播放| 熟妇无码人妻| 东京热一区二区三区无码视频| 国产内射一区亚洲| 欧类av怡春院| 97成人在线观看| 丁香婷婷久久| 免费看av在线网站网址| 日韩在线第三页| 999精品在线视频| 欧美黄网站免费观看| 热伊人99re久久精品最新地| 特级欧美视频aaaaaa| 亚洲中文精品久久久久久不卡| 香蕉国产精品视频| 久久一本日韩精品中文字幕屁孩| 男女精品视频| 国产H片无码不卡在线视频| 1级黄色毛片| 亚洲AⅤ波多系列中文字幕| 亚洲精品无码抽插日韩| 久久精品国产精品一区二区| 亚洲AV人人澡人人双人| 精品免费在线视频| 欧美人与动牲交a欧美精品| 视频一区亚洲| 国产精品对白刺激| 亚洲免费三区| 性网站在线观看| 国产正在播放| 午夜啪啪网| 久久久四虎成人永久免费网站| 色首页AV在线| 精品一区二区三区四区五区| 性69交片免费看| 四虎国产在线观看| 国产精品三区四区| 色视频国产| 玖玖精品在线| 国产精品护士| 亚洲一道AV无码午夜福利| 国产日韩欧美一区二区三区在线 | 免费国产小视频在线观看| 亚洲精品男人天堂| 国产自无码视频在线观看| 乱码国产乱码精品精在线播放 | 国产亚洲欧美在线中文bt天堂|