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

航空發動機燃燒室參數化設計

2022-02-06 08:08:42賈文杰
航空發動機 2022年6期
關鍵詞:模型設計

賈文杰

(中國飛行試驗研究院,西安 710089)

0 引言

目前,航空發動機燃燒室設計主要包括方案設計、初步設計、細節設計和試驗驗證[1]。其中,初步設計階段主要采用經驗和半經驗方法[2-5],手工計算數目繁多的參數,工作效率低且容易出錯;詳細設計階段主要采用計算流體力學(Computational Fluid Dynamics,CFD)方法對初步計算得到的幾何結構進行流場分析和結構優化,在數值求解過程中需要無數次重復“修改幾何模型-生成網格-設置邊界條件-求解-后處理”過程[5]。隨著計算機能力的提高,前處理環節(包括幾何建模、網格生成和邊界條件設置)占整個數值計算周期的時間比達到40%~80%[6],成為制約數值計算能力提高的關鍵。

針對傳統燃燒室設計方法中存在的不足,國外航空發動機公司提出了參數化設計方法。該方法最先由Tangirala等[7]提出,用變量代替常量建立幾何模型和分析計算結果,通過簡單修改模型和網格中的某些參數值,建立和分析新的模型,同時保持其它參數不變。采用參數化數值仿真方法可以有效解決傳統方法迭代計算次數多、人工介入多等問題,顯著提高計算效率,節約時間成本和人工成本。

Dawes等[8]提出,在數值仿真技術求解精度逐漸提高的背景下,將CFD方法用于燃燒室設計的技術瓶頸已經不是求解器自身的模型選擇和計算結果的精度與可靠性,而是在“幾何建模-網格生成-求解”過程中所耗費的巨大的時間成本和人工成本,而采用參數化方法解決上述問題的最大難點在于幾何建模到網格生成過程中的數據傳遞(Geometry Modeling And Grid Generation,GMGG);Tangirala等[7]分析了GE公司某環形燃燒室的摻混孔和冷卻孔參數對燃燒室出口溫度分布的影響,表明采用參數化方法在單次數值仿真周期中所需時間是傳統方法的28%;Honeywell等[9-10]自主開發了燃燒室參數化設計工具(Advanced Combustion Tools,ACT),涵蓋了燃燒室數值仿真的全部過程;Pegemanyfar等[11-12]也開發出一種知識工程(Knowledge Based Engineering,KBE)系統,可以實現從輸入設計參數和性能指標到自動計算出幾何模型和生成網格、求解和后處理的全部過程,真正實現了燃燒室的一體化設計,顯著縮短了燃燒室設計周期,驗證了參數化設計方法的有益效果。目前國外針對燃燒室的參數化數值仿真和一體化設計開展了大量研究[13-15],而中國相關研究較少,仍處于起步階段。

本文針對中國在燃燒室參數化設計方面存在的不足,開發燃燒室一體化設計平臺,并對其有效性進行驗證。

1 參數化數值仿真的難點

航空發動機燃燒室的參數化主要包括參數化建模、參數化網格生成和參數化求解及后處理3部分。目前,實現參數化數值仿真的2大難點在于如何實現各部分的參數化以及各模塊之間的數據傳遞。

1.1 變量的參數化

隨著計算機輔助設計(Computer Aided Design,CAD)技術的發展,已經有很多成熟的軟件用于幾何建模和求解。在前處理環節,由于燃燒室存在旋流器、火焰筒開孔等復雜結構,需要軟件具備強大的幾何處理能力,包括燃燒室幾何模型的建立、修復、計算域的提取以及高質量結構/非結構網格生成能力;在求解環節,需要求解器提供豐富的物理化學模型,具有靈活的適應性和較高的求解精度。除上述要求外,還需要各部分軟件具備腳本運行的能力以實現整個參數化流程的自動進行。

參數化方法的另一關鍵在于控制參數的提取。燃燒室的幾何模型一般由數百個草圖點和曲線結合回轉、拉伸、陣列和布爾等命令生成,結構復雜且參數眾多,因此無法將所有的幾何尺寸都標記為控制變量。解決上述問題的辦法為:將幾何模型腳本文件中的草圖點坐標、曲線長度等參數由常量改為變量;用表達式將幾何尺寸中不相互獨立的量相互關聯,以減少控制參數的數量;從表達式中篩選并標記出可能需要優化和反復修改的變量,作為最終可用的控制參數。在完成幾何模型的參數化后,采用類似方法修改并標記出控制網格尺寸、邊界條件和求解模型的變量。

1.2 數據傳遞

燃燒室數值仿真過程需要CAD軟件、網格生成軟件和求解及后處理軟件協同進行,但不同軟件的語言類型和語法規則不同,使得實現各軟件之間的參數傳遞成為燃燒室參數化數值仿真的難點之一。為此,本文采用Python語言編寫了燃燒室初步設計與參數化數值仿真軟件,作為協同初步設計和詳細設計的一體化設計框架。針對各模塊生成了對應的GUI界面,通過回調函數將GUI界面上的對話框輸入參數和腳本文件中提取出的控制變量相關聯,在對話框中修改數據即可同步更新腳本語言中對應的控制參數。通過Python語言編寫文件格式轉換程序實現了Salome軟件的Python腳本語言、ICEM CFD軟件的Tcl/Tk語言以及Fluent的Scheme語言之間的相互調用,同時實現了軟件的初步設計模塊和參數化數值仿真模塊之間的數據傳遞。

2 燃燒室參數化設計

2.1 燃燒室初步設計

目前,航空發動機燃燒室初步設計主要通過經驗公式和氣動熱力計算得到燃燒室的幾何結構,在給定燃燒室的性能參數和進口條件后,確定燃燒室內各部分的壓力損失、空氣流量分配和火焰筒上各孔的流量系數,代入經驗和半經驗公式中計算得到燃燒室的幾何參數,如擴壓器尺寸、參考截面(機匣橫截面最大處)尺寸、旋流器尺寸、主燃孔、摻混孔和冷卻孔參數等。燃燒室初步設計流程如圖1所示。

圖1 燃燒室初步設計流程

以目前燃燒室中常用的突擴型擴壓器初步設計為例,其結構如圖2所示。突擴型擴壓器幾何參數主要包括面積比AR、入口段直徑W、前置擴壓器長度L、擴張角θ和突擴段長度D。其中,突擴段長度D無法在設計階段確定,需要在數值模擬階段反復修正。

圖2 突擴型擴壓器結構

其中,面積比AR為

根據連續性方程可得面積比AR與前置擴壓器進口速度u3及出口速度u3.1之間的關系

式中:A3.1和A3分別為前置擴壓器出口段和入口段面積。

在擴壓器設計過程中,當燃燒室進口空氣流速u3確定時,便可以根據式(2)得到面積比AR。設前置擴壓器的總壓損失為△Pdiff,根據伯努利方程

式中:p3和q3為擴壓器進口氣流的靜壓和動壓頭;p3,1和q3,1為擴壓器出口氣流的靜壓和動壓頭。

擴壓器的靜壓恢復系數Cp為

擴壓器靜壓恢復系數如圖3所示。曲線CP*表示在給定的無量綱參數N/R1下,靜壓恢復系數的最大值和AR之間的關系;曲線表示在給定的面積比下,靜壓恢復系數的最大值和無量綱參數N/R1之間的關系。其中,橫坐標中在設計階段,當AR確定后,便可以從圖中得到最優的無量綱參數N/R1,代入式(1)中即可得到擴壓器的幾何參數。

圖3 擴壓器靜壓恢復系數

2.2 參數化幾何建模

目前,中國在燃燒室參數化幾何建模方面主要應用的商業軟件為UG,但其生成的模型無法與網格生成軟件相互關聯,只能用于燃燒室幾何拓撲結構不發生任何變化時的參數化數值仿真,局限性很大,因此本文選擇開源軟件Salome作為參數化前處理工具。

Salome軟件具備通過Python腳本語言運行的能力,在建立幾何模型時,首先采用交互式圖形界面操作的方法生成對應幾何模型的腳本文件,然后修改腳本文件并標記出特征參數,該方法相比根據Salome軟件的語法規則手動編寫幾何模型可以顯著縮短建模時間,這也是該軟件優于其它CAD軟件的特征之一。參數化建模流程如圖4所示。

圖4 參數化建模流程

第1步,在軟件中采用交互式方法創建燃燒室的幾何模型。為了避免燃燒室結構復雜造成約束過度的問題,采用模塊化建模的方法,將燃燒室分為火焰筒、旋流器和機匣3個模塊單獨建模,分別錄制各模塊的腳本文件。燃燒室模塊化建模流程如圖5所示。采用這種方法還可以提高腳本文件的通用性,為各種常用結構的火焰筒和旋流器生成模板文件并保存在數據庫中,排列組合后可以快速生成多個燃燒室模型。

圖5 燃燒室模塊化建模流程

第2步,將腳本文件中的幾何尺寸參數化。腳本文件中包含數百個與燃燒室結構相關的參數,需要使用表達式和變量對燃燒室中的參數進行簡化。例如,在生成,2級旋流器葉片草圖時,用于確定某2級葉片尺寸和位置的4個草圖點在原始腳本中表示為

pt_1=geompy.MakeVertex(2,4.5,0)

pt_2=geompy.MakeVertex(10,-4.5,0)

pt_3=geompy.MakeVertex(10,-3.5,0)

pt_4=geompy.MakeVertex(2,3.5,0)

可見在原始腳本文件中,葉片是由4個相互獨立的草圖點約束,需要通過手動計算才能得到葉片的長度、寬度和安裝角等參數,不便于控制參數的修改,使用參數化方法將程序修改為

d_2=0.5

1_2=4

degree_2=-45

Ang_2=math.pi*(degree_2/180)

pt_blade_2_mid_point_2=geompy.MakeVertex(6,-d_2,0)

pt_blade_2_stg_2=geompy.MakeVertex(6+l_2,l_2*math.tan(Ang_2)-d_2,0)

pt_blade_2_stg_3=geompy.MakeVertex(6+l_2,l_2*math.tan(Ang_2)+d_2,0)

pt_blade_2_stg_4=geompy.MakeVertex(6-l_2,-1*l_2*math.tan(Ang_2)+d_2,0)

在上述語句中,2級旋流器葉片厚度D_2=2*d_2;2級旋流器葉片長度L_2=2*l_2;Ang_2和ang_2分別表示弧度制和角度制下的2級旋流器葉片安裝角,通過給d_2和l_2賦不同的值,即可完成對旋流器葉片參數的修改,然后重新運行腳本文件即可生成新的模型,部分可以修改為變量的幾何參數及其關聯的表達式見表1。

表1 部分可以修改為變量的幾何參數及其關聯的表達式

第3步,為生成的幾何體劃分邊界并導出幾何模型。本文主要通過Salome軟件的“explode”命令和“group”命令為幾何模型劃分邊界,“explode”命令的函數表述為geompy.SubShape AIIIDs(Shape,Type),用于將實體拆分為片體并編號分類;“group”命令的函數表達為geompy.CreateGroup(Shape Name,Shape Type[“FACE”]),用于將相同或不同性質的幾何體組合為1個組,具體實現思路為:使用“explode”命令將燃燒室實體幾何模型拆分為片體,然后使用“group”命令根據邊界條件將片體整合為多個組,不同的組對應不同的邊界,如:“air_inlet”、“symmetry”以及“pressure_outlet”等,最后將不同的組分別導出為.stl格式的幾何文件,從Salome中導出的組以及在ICEM CFD中自動識別出的邊界條件如圖6、7所示。從圖中可見,Salome中劃分的組可以與ICEM CFD中的邊界條件準確對應。

圖6 從Salome中導出的組

圖7 邊界條件

2.3 參數化網格生成

本文使用商業軟件ICEM CFD作為參數化網格生成工具。ICEM腳本文件即.rpl語言,用于記錄網格生成過程中的完整操作信息,其生成網格的部分語言表示為

ic_geo_set_family_params TOP_AND_BOTTOM.1 prism 1 emax 6.0 ehgt 0.2 hrat 1.2 nlay 4 erat 0 ewid 0 emin 0.0 edev 0.0 prism_height_limit 0 law-1 split_wall 0 internal_wall 0

腳本中包含了邊界選擇(TOP_AND_BOTTOM)、邊界層初始高度(ehgt 0.2)和邊界層生長率(hrat 1.2)等網格生成過程中的全部內容。在網格生成過程中,提取可能需要多次修改的變量設置為自定義變量,并在設計平臺中設置對應的輸入框,ICEM CFD中可變的控制參數見表2。當需要修改網格參數時,可以通過主程序與腳本之間的數據接口在主程序用戶界面上的對話框中修改,保存后運行腳本即可重新生成網格。

表2 ICEM CFD中可變的控制參數

2.4 參數化求解及后處理

本文中選擇的參數化求解及后處理工具為商業軟 件Fluent。Fluent可 以使用Scheme語言編寫的日志文件來自動執行運算,非常適用于執行一系列具有相似參數和邊界條件的模型計算。只需要手動設置某一工況的參數,在執行其余相似工況的燃燒室數值計算時,重新運行日志文件即可重復所有與求解器設置相關的操作。

除了執行運算,Fluent還具備強大的后處理功能,可以方便快捷地顯示流場區域,創建等值線云圖、矢量圖、動畫和視頻等。在完成計算后,錄制顯示計算結果和流場分析的日志文件,然后重新運行即可顯示并輸出計算結果。

3 一體化設計平臺

航空發動機燃燒室參數化設計平臺主要用于燃燒室幾何結構的初步設計,并結合參數化數值仿真對零維計算得到的幾何結構進行詳細優化設計,最終輸出可以用于試驗研究的燃燒室幾何結構。

3.1 軟件功能

根據燃燒室設計流程與要求,本軟件平臺可以實現的主要功能為:

(1)輸入燃燒室性能參數與進口條件,得到燃燒室部分氣動參數;

(2)燃燒室幾何結構設計。包括擴壓器、參考截面、旋流器和火焰筒開孔等設計;

(3)參數化數值仿真過程中幾何參數和網格參數的修改、軟件和程序的后臺運行和調用,計算結果的顯示;

(4)不同格式文件的數據傳遞、調用和保存。

3.2 軟件界面設計

燃燒室參數化設計平臺采用Python語言的Tkinter模塊定制生成。燃燒室一體化設計平臺界面(如圖8所示)主要由位于頂部的菜單欄、左側的導航窗口與中間的計算窗口3部分組成。

圖8 燃燒室一體化設計平臺界面

導航窗口主要分為燃燒室初步設計和參數化數值仿真2個模塊,各部分均按照設計流程的順序依次排布。為了便于使用、修改和維護,按模塊化編程的思路進行軟件GUI設計,將計算窗口分為性能參數計算、擴壓器設計、燃燒室參考截面設計、旋流器設計、火焰筒開孔設計、參數化幾何以及參數化網格與求解幾部分,在使用時分別按順序點擊左側導航欄,就會彈出對應的設計窗口,輸入對應的參數后,點擊“計算”按鈕得到對應的輸出參數,然后點擊“保存”按鈕將數據保存至EXCEl文件中。

3.3 不同模塊間的數據傳遞

燃燒室初步設計與參數化數值仿真方法的主要難點之一在于不同模塊之間的數據傳遞,數據傳遞流程如圖9所示。在實際的燃燒室設計過程中,需要各模塊和函數之間相互調用。例如,在燃燒室初步設計環節中,火焰筒開孔設計和旋流器設計均與參考截面參數相關,因此需要在程序中實現各設計模塊的參數傳遞。

圖9 數據傳遞流程

在程序設計過程中,使用global語句將函數中的局部變量放入全局作用域,然后使用import語句將包含全局變量的程序模塊導入程序中,導入模塊和全局變量的部分語句如下:

from parametric_1 import p1,p3,p4,q3#從parametric_1模塊導入全局變量

def calculate_3():#計算燃燒室參考截面的函數

global s1#將流阻系數定義為全局變量

s1=float(obstacle.get())

s2=float(refdiam.get())

s3=float(ang.get())

s4=float(arate.get())

global s5#Aref面積,將其定義為全局變量

s5=math.sqrt(287.04/2*(p1*math.sqrt(p4)/p3)**2*s1*(q3**-1))*1000000 global s6 # AL面積

4 軟件有效性驗證

在完成軟件設計后,需要對軟件的有效性進行驗證。燃燒室扇段主要性能參數見表3。

表3 燃燒室扇段主要性能參數

本文設計的燃燒室為帶3級直葉片軸向旋流器的環形燃燒室,無主燃孔和摻混孔,除冷卻氣體外的其余氣體全部從頭部進入。取燃燒室的1/12進行研究,燃燒室和旋流器的幾何結構分別如圖10、11所示。

圖10 燃燒室幾何結構

圖11 旋流器幾何結構(剖視)

為了驗證軟件的有效性,在計算得到燃燒室的幾何模型后,分別探究了第1級旋流器葉片安裝角、第3級旋流器葉片安裝角和冷卻孔直徑對燃燒室出口溫度分布系數(Overall Temperature Distribution Factor,OTDF)的影響,數值計算方案劃分如圖12所示。計算采用基于壓力的求解器,湍流模型為RealizableK-ε 2方程模型,燃料為航空煤油(C12H23)。燃燒模型為非預混燃燒中的混合分數/PDF模型,使用DPM模型模擬燃油噴射,粒徑分布服從R-R分布,使用質量流量入口和壓力出口邊界條件。經過網格無關性驗證,所有計算方案的網格數量均大于495萬。

圖12 數值計算方案劃分

4.1 第1級葉片安裝角對燃燒室出口參數的影響

不同第1級葉片安裝角方案的回流區輪廓如圖13所示。雖然根據理論計算,葉片安裝角為30°時第1級旋流器的旋流數為0.495,不滿足常規發動機燃燒室設計中要求各級旋流器旋流數≥0.6[5]的條件,在5種方案下都形成了穩定的回流區。文獻[16]分別研究了在冷態和熱態環境中,內級旋流器旋流數為0.5、外級旋流器旋流數分別為±0.56時,旋流器出口流場的流動特性,發現在2級旋流器反向時,在冷態和熱態中均出現了回流區,與本文中的計算結果相符,說明單級旋流器和多級旋流器的流動特征不同,旋流數不是判斷燃燒室中能否形成回流區的充要條件。因此在設計時,不僅需要考慮旋流數和進口氣流Re的影響,還需要考慮各級旋流器之間的相互作用,文獻[16]也印證了這一結論。

圖13 不同第1級葉片安裝角方案的回流區輪廓

不同第1級葉片安裝角方案下燃燒室出口參數如圖14所示。從圖中可見,保持第3級旋流器葉片安裝角為45°,冷卻孔直徑為2 mm不變,在第1級旋流器葉片安裝角為30°時,回流區軸向尺寸較大,使得燃燒不充分,在燃燒室出口產生局部高溫區(超過2400 K),造成溫度分布不均勻;當第1級旋流器葉片安裝角為60°時,回流區更靠近火焰筒入口,高溫區遠離燃燒室出口,導致出口平均溫度過低,OTDF較大。結合出口溫度分布要求,本文中第1級旋流器葉片安裝角最佳取值范圍為45°~50°。

圖14 不同第1級葉片安裝角方案下燃燒室出口參數

4.2 第3級葉片安裝角對燃燒室出口參數的影響

上述計算中得到的出口OTDF最小值為0.264,超過設計要求,需要對燃燒室幾何結構繼續優化,不同第3級葉片安裝角方案中燃燒室出口參數如圖15所示。

圖15 不同第3級葉片安裝角方案中燃燒室出口參數

從圖中可見,當第3級旋流器葉片安裝角為45°~60°時,燃燒室出口OTDF基本保持不變。不同方案中火焰筒內燃油質量Q分數分布如圖16所示。從圖中可見,隨著第3級葉片安裝角的增大,燃油質量分數沿徑向擴張角增大,沿軸向貫穿距離減小,霧化性能增強,燃油分布更加合理;當葉片安裝角為45°~50°時,燃油得到充分霧化擴散,且未打到火焰筒壁面上,燃油分布最為合理。因此本方案中第3級旋流器葉片安裝角的最佳范圍為45°~50°。

圖16 不同方案中火焰筒內燃油質量分數分布

4.3 冷卻孔徑對燃燒室出口參數的影響

在優化旋流器結構過程中,為了限制網格數量,選擇冷卻孔直徑為2 mm,大于常規燃燒室的氣膜冷卻孔直徑,導致最終計算得到的出口OTDF的最小值為0.264,不滿足燃燒室設計需求。為了進一步優化燃燒室結構,同時測試軟件對于冷卻孔的參數化能力,分析了冷卻孔直徑對燃燒室出口溫度的影響。在冷卻孔直徑分別為2.0、1.5和1.2 mm時,燃燒室出口參數如圖17所示。從圖中可見,隨著冷卻孔直徑減小,燃燒室OTDF接近設計要求,這是由于隨著冷卻孔直徑減小,數量增加,冷卻效果增強且與主流摻混較為充分,使得燃油燃燒更加充分,出口溫度分布得到顯著優化。

圖17 不同冷卻孔直徑時燃燒室出口參數

在驗證軟件有效性過程中,對比統計了傳統方法和參數化方法在數值求解過程中消耗的時間,見表4(表中為多次計算方案中消耗時間的均值)。從表中可見,采用參數化方法可以顯著縮短數值計算周期,驗證了本文提出的參數化方法的有益效果。

表4 傳統方法和參數化方法的時間對比

5 結論

(1)提出一種參數化燃燒室設計平臺,通過將燃燒室初步設計過程需要多次修改的參數改為變量并提取出來,使用Python編程實現了從幾何建模、網格生成到求解和后處理環節的完全自動化,建立了燃燒室的一體化設計平臺;

(2)通過該平臺設計了一種帶直葉片3級旋流器的燃燒室幾何模型,并使用參數化數值仿真方法對其幾何結構進行優化,經過優化后的燃燒室出口OTDF降為0.235,接近設計要求值0.2;

(3)通過計算分析,使用參數化數值仿真方法能夠將建模時間縮短為原來的1.4%,網格生成時間縮短為原來的25%,求解時間縮短為原來的13.3%,驗證了本文中提出的參數化數值仿真方法的有益效果和一體化設計平臺的可靠性。

猜你喜歡
模型設計
一半模型
重要模型『一線三等角』
何為設計的守護之道?
現代裝飾(2020年7期)2020-07-27 01:27:42
重尾非線性自回歸模型自加權M-估計的漸近分布
《豐收的喜悅展示設計》
流行色(2020年1期)2020-04-28 11:16:38
瞞天過海——仿生設計萌到家
藝術啟蒙(2018年7期)2018-08-23 09:14:18
設計秀
海峽姐妹(2017年7期)2017-07-31 19:08:17
有種設計叫而專
Coco薇(2017年5期)2017-06-05 08:53:16
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
主站蜘蛛池模板: 天天干天天色综合网| 国产91小视频| 四虎精品黑人视频| 97综合久久| 国产亚洲精品无码专| 人妻丰满熟妇av五码区| 亚洲 成人国产| 久久精品无码中文字幕| 国产精欧美一区二区三区| 亚洲午夜福利精品无码| 欧美一级高清免费a| 青青操国产| 久久永久视频| 国产精品乱偷免费视频| 久久99国产精品成人欧美| 极品国产一区二区三区| 无码'专区第一页| 54pao国产成人免费视频| 中文字幕无线码一区| 无码国内精品人妻少妇蜜桃视频| 国产成人亚洲精品蜜芽影院| 国产永久在线视频| av午夜福利一片免费看| 精品欧美一区二区三区久久久| 国产嫖妓91东北老熟女久久一| 国产乱人视频免费观看| 欧美性天天| 色婷婷色丁香| 国产97视频在线| 亚洲精品第五页| 最新日本中文字幕| 中文字幕啪啪| 国产精品美女在线| 国产丝袜一区二区三区视频免下载| 亚洲成a人片| 9966国产精品视频| 亚洲乱强伦| 99热最新网址| 亚洲欧洲自拍拍偷午夜色| 99久久精彩视频| 欧美日韩一区二区三区在线视频| 尤物亚洲最大AV无码网站| 91久久国产综合精品| 午夜精品区| 国产精品一区二区在线播放| 国产一区免费在线观看| 亚洲天堂视频网站| 午夜日b视频| 中文字幕在线免费看| 国产Av无码精品色午夜| 日韩视频免费| 91在线中文| 久久黄色一级片| 国内精品小视频福利网址| 国产av一码二码三码无码| 午夜性刺激在线观看免费| 欧美黄网站免费观看| 中文字幕av一区二区三区欲色| 亚洲国产精品无码AV| 亚洲制服丝袜第一页| 日韩不卡高清视频| 国产福利在线观看精品| 精品少妇人妻一区二区| 91在线一9|永久视频在线| 福利在线不卡| 亚洲高清资源| 亚洲美女AV免费一区| 成年人国产网站| 在线观看国产精品日本不卡网| 国产在线观看一区精品| 亚洲综合网在线观看| 欧美一级夜夜爽www| 色欲不卡无码一区二区| 青青青国产视频手机| 1024国产在线| 亚洲AⅤ无码国产精品| 亚洲欧美日韩久久精品| 无码国产伊人| 一区二区三区毛片无码| 国产av无码日韩av无码网站| 任我操在线视频| 全部无卡免费的毛片在线看|