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

基于離散伴隨求解器的鐵路下導風工程外形優化研究

2019-09-11 12:49:38陳柏羽程建軍辛林桂丁泊淞
鐵道科學與工程學報 2019年8期
關鍵詞:風速鐵路優化

陳柏羽,程建軍,辛林桂,丁泊淞

基于離散伴隨求解器的鐵路下導風工程外形優化研究

陳柏羽,程建軍,辛林桂,丁泊淞

(石河子大學 水利建筑工程學院,新疆 石河子 832003)

為提高鐵路沿線前傾式下導風工程輸導效果,基于離散伴隨求解器對前傾式下導風工程導風板進行局部外形優化,并運用數值模擬分析優化前后下導風工程控制下的風沙流場。研究結果表明:導風板形狀敏感性存在由中部向兩端逐漸增大的分布規律,其數值大小與風速呈正相關,通過優化導風板外形可以有效降低導風板阻力,減少氣流經過導風板時損耗的動能,氣流通過保護區時具有更高的運動風能;導風板外形優化可以增加導風板有效輸導高度,擴大保護區范圍,提高保護區風速,導風板聚集壓縮加速氣流的效果更好;導風板外形優化可以增強導風板輸沙效果,減少鐵軌附近的積沙量,增大鐵路與鐵路背風側積沙范圍之間的距離,延長下導風工程使用年限。

下導風工程;離散伴隨求解器;外形優化;數值模擬

下導風工程又稱聚風板工程,由柵欄工程發展而來,廣泛應用于防治風沙、風吹雪對鐵路和公路的危害[1]。下導風工程的排列和結構形式多樣,根據導風板傾角不同分為直立式(傾角=90°)、前傾式(傾角<90°)和后傾式(傾角>90°)3種[2]。導風板寬度、傾角及下口高度作為下導風工程的主要結構特征,直接影響著下導風工程的輸導效果,已有學者對此進行了大量研究:劉賢萬等[1]通過風洞試驗研究導風板傾角對輸沙效果的影響,發現導風板采用偏離直立的傾角可以有效減少板前和板后沙粒堆積。程建軍等[3]通過三維數值模擬對比分析直立式和前傾式導風板的輸沙效果,結果表明前傾式導風板的輸沙效果明顯優于直立式導風板。辛林桂等[4]通過數值模擬對下導風工程的主要影響參數進行了優化分析,發現導風板與風向夾角為60°、下口高度為2 m時輸沙效果較好。LIU等[5]通過數值模擬研究具有不同下口高度的直立式下導風工程,發現導風板背風側保護區范圍與下口高度呈正相關。相較于導風板寬度、傾角及下口高度等,導風板外形通常被認為是影響下導風工程輸導效果的次要因素,因此關于導風板外形優化的研究較少。另外,導風板氣動外形優化涉及諸多變量,傳統試錯法在處理這類問題時存在一定局限性,往往需要反復實驗才能達到設計目標[6]。而伴隨方程法是一種基于梯度算法的優化方法,其計算量與目標變量數目無關,可大幅降低計算成本[7]。1974年Pironneau[8]首次將伴隨方程法應用于流體動力學,隨后Jameson[9]將其擴展至航空領域并用于飛行器翼型優化,此后伴隨方程法一直作為空氣動力學領域的主要優化方法[10]。目前已有學者基于伴隨方程法開發出可用于工程實際的軟件工具[11?14]。其中,離散伴隨求解器作為一種基于離散伴隨方程法開發的CFD優化工具,旨在進一步開發和改進流體系統的性能。它擴展了傳統流體求解器的分析范疇,僅通過一次計算就能提供流體系統詳細的性能敏感性數據,例如阻力、升力和壓力等,敏感性數據的可視化可指導設計人員對目標幾何形狀進行正確的改進,同時離散伴隨求解器內置的網格變形工具也可根據性能敏感性數據進行網格自由變形,從而能實現智能化設計。目前,隨著CFD技術的不斷發展和成熟,基于伴隨方程法的氣動外形優化設計技術日趨完善,已在空氣動力學領域發揮了重要作用[15?17],然而尚未有基于伴隨求解技術進行風沙防護工程外形優化方面的研究。本文基于離散伴隨求解器對前傾式下導風工程導風板進行局部外形優化,在保證導風板風向投影面積、傾角及下口高度一定的情況下,通過降低導風板阻力來減少氣流經過導風板時產生的動能損耗,進一步提高導風板輸導效果。為下導風工程的優化設計和應用提供理論依據。

1 數值模擬方法

1.1 幾何建模及網格劃分

運用AutoCAD對前傾式下導風工程進行幾何建模。由于下導風工程在結構上呈幾何對稱,因此可將其簡化為二維模型計算。計算域長度為150 m,高度為30 m,導風板放置在距入口60 m處,導風板寬度為3 m,與水平風向夾角60°,下口高度為1.5 m。計算域模型如圖1所示。

流體計算域模型網格劃分類型采用四邊形主導(Quadrilateral Dominant),并對導風板附近的網格進行加密,增長率為1.1,網格質量平均值為0.983 8。計算域模型網格數量約為5×104個。計算域網格劃分結果如圖2所示。

圖1 計算域模型及邊界條件示意圖

圖2 計算域網格劃分示意圖

1.2 計算參數及控制方程

根據空氣動力學原理,風沙流馬赫數均小于0.3時被認為是不可壓縮流動[17],模型入口邊界條件為速度入口(Velocity-inlet)。由于計算域中流態充分發展條件下才采用自由出流,而此模型中出口不能保證自由出流,故模型出口邊界條件為壓力出口(Pressure-outlet),壓力差為0。計算域模型上下邊界條件為Wall,下邊界粗糙度為0.001,上邊界粗糙度為0。空氣密度k=1.225 kg/m3,空氣動力黏度=1.789×10?5Pa·s,壓力為常壓。風沙流中沙粒粒徑通常為0.075~0.25 mm,本文取沙粒粒徑s=0.15 mm,沙粒密度s=2 650 kg/m3,黏度=0.004 7 Pa·s[18],初始沙粒體積分數為1%,類型為FLUID。入流風速選取10,15和20 m/s 3種。

計算域求解模型采用歐拉雙流體模型,附加-湍流模型,湍流強度=0.05,湍流半徑=0.5,并選取Syamlal-O’Brien 曳力模型[19]。所需控制方程主要有連續性方程、動量方程和-湍流方程[20]。方程組求解計算方法采用SIMPLEC算法,此算法適用于不可壓縮流動,并且可以加快迭代過程的收斂,各分量值收斂標準為10?6量級。

2 離散伴隨方法

2.1 優化流程

導風板外形優化系統主要包括幾何建模工具、網格劃分工具、CFD程序及離散伴隨求解器。其中離散伴隨求解器的工作流程與傳統CFD程序類似,即參數設置、初始化伴隨場并迭代至收斂,其計算量也與傳統CFD程序相當。導風板外形優化設計流程如下:1) 計算域幾何建模及網格劃分;2) 利用CFD程序求解計算域模型并得到穩定收斂的流場解;3) 離散伴隨求解器基于梯度算法求解伴隨方程;4) 輸出模型的敏感性數據;5) 離散伴隨求解器內置的網格變形工具根據敏感性數據進行網格自由變形;6) 重復上述過程直至目標變量收斂,設計目標達成后,將最終模型網格化并通過CFD程序重新評估該幾何形狀。基于離散伴隨求解器進行導風板外形優化的流程如圖3所示。

2.2 計算參數及網格變形控制

以導風板所受阻力為目標變量,阻力敏感性取向(Sensitivity orientation)設定為最小化(Minimize)。離散格式采用二階迎風格式(Second order upwind),各分量值收斂標準為10?6量級。

圖3 優化流程示意圖

網格變形對象為導風板,網格變形比例因子(Freeform scale factor)為0.1。單次優化的阻力降低幅度為2%。網格變形區域高度為1 m,寬度為5 m。同時,網格變形區域內和坐標方向控制點數量各為20個。網格變形區域如圖4所示。

圖4 網格變形區域示意圖

3 導風板外形優化結果與分析

下導風板主要借助風的動力作用,使風沙流通過導風板下口時,風道斷面減小導致風速增大,從而使風沙流加速通過線路并清除線路上的積沙[3]。然而,導風板在聚集加速氣流的過程中會將導風板所在流層的流體動能和位能轉化為加強導風板上、下端流體運動的動能和位能[1]。這一過程勢必會損耗氣流的部分動能,而損耗動能的大小可通過導風板所受阻力來體現[2]。因此,在保證導風板風向投影面積、傾角和下口高度一定的情況下,對導風板局部外形進行正確的減阻設計,可以有效降低氣流經過導風板時產生的動能損耗,增加導風板下端氣流的運動風能,進一步提高導風板輸導效果。

3.1 形狀敏感性分析

基于穩定收斂的流場解,離散伴隨求解器可從中推導出流體系統詳細的敏感性數據。值得注意的是,敏感性數據僅針對于指定的目標變量及當前的入流條件,伴隨優化結果隨目標變量和入流條件的變化而變化。因此,在使用伴隨方法進行氣動外形優化時通常選取某一風速下求解出的流場信息作為伴隨優化依據[6]。圖5為不同風速下導風板形狀敏感性示意圖。可以看出,不同風速下導風板形狀敏感性均存在由中部向兩端逐漸增大的分布規律,且對同一區域而言,形狀敏感性的分布和方向不受風速影響,僅在數值大小上與風速呈正相關。由于風速會影響敏感性數據,從而影響伴隨優化結果,所以,為便于研究,本文以風速15 m/s時計算出的流場信息作為伴隨優化依據,進行導風板外形局部優化。

從圖5可知,改變導風板上下端外形可以顯著降低導風板阻力,但導風板上端變形會向后彎曲,導致更多氣流從導風板上端分離,經導風板聚集加速后通過下口的氣流減少,這降低了導風板輸沙效果。因此本文通過約束網格變形區域,僅對導風板下部進行減阻外形優化。

(a) 風速:10 m/s;(b) 風速:15 m/s;(c) 風速:20 m/s

圖6 優化前后導風板外形對比圖

優化前后導風板外形對比如圖6所示。可以看出,優化后導風板下部與導風板主體呈一定角度,轉角處近似由一段弧形平滑過渡,導風板下部與氣流夾角變小。導風板下部與氣流夾角變小使得導風板下端對氣流剪切作用減弱,氣流損耗的運動風能減少,此外,轉角處由弧形過渡也使氣流沿導風板流動更加平順,進一步減少了氣流流動過程中的動能損耗。

圖7為優化過程中導風板阻力變化曲線。可以看出,在最初的幾輪優化循環中,阻力降低顯著,隨優化循環次數的增加,導風板單次優化降低的阻力值逐漸減少并最終趨于收斂。經30輪優化循環后,導風板阻力由最初的660.68 N降低至543.01 N,導風板阻力降低17.81%。這表明氣流經過導風板時損耗的動能顯著降低,氣流能以更高的動能通過導風板下口。同時阻力降低也表明下導風工程在極端風環境下的抗傾覆能力增強,下導風工程結構可靠度相對較高。

圖7 導風板阻力變化曲線

3.2 流場分析

導風板所在流層氣流沿導風板流動時會在板面某一點產生分離,即氣流分離點,該點以上氣流沿板面向上流動并加速導風板上端,而該點以下氣流沿板面向下流動并加速通過下口。如圖8所示。據此,本文定義氣流分離點至導風板下端點之間的豎直高度為導風板有效輸導高度。在導風板風向投影面積、傾角及下口高度一定的情況下,導風板有效輸導高度越高,表明經導風板聚集壓縮加速的氣流越多,導風板輸導效果越好。以風速15 m/s為例,優化前導風板有效輸導高度為1.74 m,而優化后導風板有效輸導高度為2.05 m。相較而言,氣流分離點向上提升0.31 m,導風板有效輸導高度增加17.82%。由此可見,優化導風板外形可以將其所在流層中更多的氣流聚集壓縮并加速通過下口,導風板輸導效果更加顯著。

導風板外形決定著導風板周圍氣流的流向及流速,是影響導風板輸導效果的結構特征之一,其周圍流場分區與普通擋墻類似,也有減速區,渦流區及加速區生成。特別的是,導風板的加速區可根據所處位置不同分為上下兩部分,下方加速區通常稱為保護區,其范圍和風速大小可以反映導風板輸導效果。如圖9所示。形成流場分區的原因在于導風板是一種帶尖緣的鈍體,當氣流受擠壓繞過導風板兩端時會在上下端點處產生分離形成強剪切層,強剪切層兩側的壓力差迫使流線向導風板中部彎曲,從而在導風板背風側形成渦流區,在導風板上下端形成加速區[21]。

(a) 優化前;(b) 優化后

從圖9可知,優化前后的流場分區范圍明顯不同:相較于優化前,優化后導風板保護區長度明顯增加,而迎風側減速區、上方加速區及背風側渦流區范圍則明顯減小。下方保護區變化趨勢表明優化導風板外形有利于增大導風板吹刮范圍,提高導風板輸導效果。而上方加速區及背風側渦流區變化趨勢表明優化導風板外形可以減弱導風板背風側氣流對列車及受電弓等設施的干擾,更利于列車高速通行。為進一步比較優化導風板外形對保護區風速的影響,繪制不同風速下導風板背風側0.5 m高度處風速變化曲線,如圖10所示。可以發現,優化后的保護區風速明顯高于優化前,并且隨風速增大,外形優化體現出的優勢越顯著。例如,在3種風速條件下,優化前背風側0.5 m高度處風速極大值分別為13.06,20.51和26.31 m/s,而優化后0.5 m高度處風速極大值分別為13.68,21.70和27.91 m/s,較優化前分別提高4.75%,5.80%和6.08%。可見,優化導風板外形可以增大導風板保護區范圍,同時也可以提高保護區內氣流流速,導風板輸導效果進一步增強。

(a) 優化前;(b) 優化后

圖10 導風板背風側0.5 m高度處風速變化曲線

3.3 鐵路積沙分析

圖11為不同風速時導風板外形優化前后鐵路積沙分布情況。可以發現,不論風速如何,2種導風板控制下的鐵路附近均有積沙分布。鐵路附近容易積沙的主要原因是由于枕木、鐵軌等結構凹凸不平,對鐵路附近流場產生擾動,加速區無法完全覆蓋鐵路,另外,軌道凹形槽結構也導致鐵軌附近存在小范圍渦流區,易使沙粒在鐵軌附近堆積[4]。

從圖11可知,風速為10 m/s時,在鐵路兩側及鐵軌之間均勻積沙分布,導風板輸導效果不太理想;風速為15 m/s時,鐵路迎風側基本無積沙,鐵軌之間積沙有所減少,鐵路背風側積沙范圍與鐵路之間的距離增大,導風板輸導效果有所提高;風速為20 m/s時,鐵路迎風側與鐵軌之間積沙進一步減少,且鐵路背風側積沙范圍與鐵路之間的距離進一步增大,導風板輸導效果較好。可見,鐵路周圍積沙分布有隨風速的增大而減少的變化規律,這也證實了風速是影響導風板輸導效果的主導因素。另外,鐵路積沙分布也受導風板外形影響:在相同風速下,優化前鐵軌附近積沙量較多,鐵路背風側積沙范圍與鐵路之間的距離相對較近;而優化后鐵軌附近積沙量明顯減少,鐵路背風側積沙范圍與鐵路之間的距離相對較遠。由此可見,優化導風板外形可增強下導風工程輸導效果,延長下導風工程使用年限。

圖11 下導風工程控制下鐵路積沙分布圖

注:圖中“A”代表優化前,“B”代表優化后。

4 結論

1) 前傾式導風板形狀敏感性存在由中部向兩端逐漸增大的分布規律,其數值大小與風速呈正相關。通過對導風板下部外形優化,導風板所受阻力由660.68 N降低至543.01 N,阻力值降低17.81%,氣流經過導風板時損耗的動能顯著降低,導風板輸導效果增強。

2) 優化與否,導風板周圍均有明顯的減速區、渦流區及加速區生成。但優化導風板外形可以減小導風板迎風側減速區、背風側回流區及上方加速區范圍,增大下方保護區范圍及風速,并且可以增加導風板有效輸導高度,導風板聚集壓縮加速氣流效果更加顯著。

3) 鐵路結構易使沙粒堆積在鐵軌附近及線路背風側,通過優化導風板外形可以進一步減少鐵軌之間的積沙量,增大鐵路背風側積沙范圍與鐵路之間的距離,延長下導風工程使用年限。

[1] 劉賢萬, 凌裕泉, 賀大良, 等. 下導風工程的風洞實驗研究——[1]平面上的實驗[J]. 中國沙漠, 1982, 2(4): 14?21. LIU Xianwan, LING Yuquan, HE Daliang, et al. Studies on the effects of fence in wind tunnel. part 1. Experiment with level snuface[J]. Journal of Desert Research, 1982, 2(4): 14?21.

[2] 王中隆. 中國天山艾肯達坂透風式下導風防雪工程[J].山地學報, 1994(4): 193?200. WANG Zhonglong. Application of blower fences with a permeability new technology on Aiken pass in Tianshan mountain of China[J]. Mountain Research, 1994(4): 193?200.

[3] 程建軍, 智凌巖, 薛春曉, 等. 鐵路沿線下導風板對風沙流場的控制規律[J]. 中國鐵道科學, 2017, 38(6): 16?23. CHENG Jianjun, ZHI Lingyan, XUE Chunxiao, et al. Control law of lower air deflector for sand flow field along railway[J]. China Railway Science, 2017, 38(6): 16?23.

[4] 辛林桂, 程建軍, 王連, 等. 基于Design Exploration方法對鐵路下導風工程關鍵設計參數的優化[J].中國沙漠, 2018(6): 1193?1199. XIN lingui, CHENG Jianjun, WANG Lian, et al. Optimization of key design parameters of the aviation baffle engineering of railway based on design exploration method[J]. Journal of Desert Research, 2018(6): 1193? 1199.

[5] LIU D, LI Y, WANG B, et al. Mechanism and effects of snow accumulations and controls by lightweight snow fences[J]. Journal of Modern Transportation, 2016, 24(4): 261?269.

[6] Munoz-Paniagua J, Garcia J, Crespo A, et al. Aerodynamic optimization of the nose shape of a train using the adjoint method[J]. Journal of Applied Fluid Mechanics, 2015, 8(3): 601?612

[7] 陳頌, 白俊強, 史亞云,等. 民用客機機翼/機身/平尾構型氣動外形優化設計[J]. 航空學報, 2015, 36(10): 3195?3207. CHEN Song, BAI Junqiang, SHI Yayun, et al. Aerodynamic shape optimization design of civil jet wing-body-tail configuration[J]. Acta Aeronautica et Astronautica Sinica, 2015, 36(10): 3195?3207.

[8] Pironneau O. On optimum design in fluid mechanics[J]. Journal of Fluid Mechanics, 1974, 64(1): 97?110.

[9] Jameson A. Aerodynamic design via control theory[J]. Journal of scientific computing, 1988, 3(3): 233?260.

[10] Tzanakis A. Duct optimization using CFD software ANSYS fluent adjoint solver’[D]. Goteborg, Sweden: Chalmers University of Technology, 2014.

[11] Nielsen E J, Diskin B, Yamaleev N K. Discrete adjoint-based design optimization of unsteady turbulent flows on dynamic unstructured grids[J]. AIAA Journal, 2009, 48(6): 1195?1206.

[12] Mader C, Martins J, Marta A. Towards aircraft design using an automatic discrete adjoint solver[C]// 18th AIAA Computational Fluid Dynamics Conference. 2007: 3953.

[13] Martín M, Andrés E, Wildham M, et al. CAD-based aerodynamic shape design optimization with the DLR tau code[C]// 27th International Congress of the Aeronautical Sciences, France. 2010.

[14] Carrier G, Destarac D, Dumont A, et al. Gradient-based aerodynamic optimization with the elsA software[C]// 52nd Aerospace Sciences Meeting, 2014: 0568.

[15] Petrone G, Hill C, Biancolini M. Track by track robust optimization of a F1 front wing using adjoint solutions and radial basis functions[C]// 32nd AIAA Applied Aerodynamics Conference. 2014: 3174.

[16] Dhert T, Ashuri T, Martins J R R A. Aerodynamic shape optimization of wind turbine blades using a Reynolds- averaged Navier–Stokes model and an adjoint method[J]. Wind Energy, 2017, 20(5): 909?926.

[17] 程建軍, 龐巧東. 戈壁強風區擋風構筑物限制下列車動力學特性分析[J]. 鐵道標準設計, 2013(1):1?5. CHENG Jianjun, PANG Qiaodong. Analysis on train aerodynamics characteristics under different types of windbreak structures at strong wind zone in Gobi[J]. Railway Standard Design, 2013(1): 1?5.

[18] HUANG N, ZHENG X J, ZHOU Y H. A multi-objective optimization method for probability density function of lift-off speed of wind-blown sand movement[J]. Advances in Engineering Software, 2006, 37(1): 32?40.

[19] 辛國偉, 程建軍, 楊印海. 鐵路沿線掛板式沙障開孔特征與風沙流場的影響研究[J]. 鐵道學報, 2016, 38(10): 99?107. XIN Guowei, CHENG Jianjun, YANG Yinhai. Study on effect of characteristics of hanging-type concrete sand barrier opening and wind-sand field along railway[J]. Journal of the China Railway Society, 2016, 38(10): 99?107.

[20] Lee S J, Park K C, Park C W. Wind tunnel observations about the shelter effect of porous fences on the sand particle movements[J]. Atmospheric Environment, 2002, 36(9): 1453?1463.

[21] 高永平. 蘭新鐵路的防風工程[J]. 路基工程, 2010(增1): 100?103. GAO Yongping. Windproof engineering of Lanzhou- Xinjiang railway[J]. Subgrade Engineering, 2010(Suppl 1): 100?103.

Shape optimization of blower fences of railway based on discrete adjoint solver

CHEN Boyu, CHENG Jianjun, XIN Lingui, DING Bosong

(College of Water Resources and Architectural Engineering, Shihezi University, Shihezi 832003, China)

In order to reveal the transport effects of the blower fence along the railway, a discrete adjoint solver was used to optimize the local shape of the blower fence. The wind and sand flow field under the control of the blower fence before and after the optimization was studied by numerical simulation. The results show that: The shape sensitivity of the air deflector is gradually increased from the middle to the both ends of the air deflector, and its value is positively correlated with the wind speed. By optimizing the shape of the air deflector, the air deflector resistance can be effectively reduced, and the kinetic energy lost when the airflow passes through the air deflector is reduced, and the airflow passes through the protection zone to have higher wind energy; The shape optimization of the air deflector can increase the effective transmission height of the air deflector, expand the range of the protection zone, increase the wind speed of the protection zone. The effect of the air deflector gathering and compressing the accelerated airflow is better; The shape optimization of air deflector can enhance the effect of sediment transport, reduce the amount of sediment accumulated near the railroad, increase the distance between the range of sediment accumulated on the leeward side of the railway and the railway, and prolong the service life of the blower fence.

blower fence; discrete adjoint solver; shape optimization; numerical simulation

X169

A

1672 ? 7029(2019)08? 1923 ? 08

10.19713/j.cnki.43?1423/u.2019.08.008

2018?11?05

國家自然科學基金資助項目(51641808,51568057,51268050)

程建軍(1979?),男,河北衡水人,教授,博士,從事道路風沙防治研究;E?mail:chengdesign@163.com

(編輯 涂鵬)

猜你喜歡
風速鐵路優化
超限高層建筑結構設計與優化思考
房地產導刊(2022年5期)2022-06-01 06:20:14
民用建筑防煙排煙設計優化探討
關于優化消防安全告知承諾的一些思考
一道優化題的幾何解法
沿著中老鐵路一路向南
云南畫報(2021年12期)2021-03-08 00:50:54
基于Kmeans-VMD-LSTM的短期風速預測
基于最優TS評分和頻率匹配的江蘇近海風速訂正
海洋通報(2020年5期)2021-01-14 09:26:54
鐵路通信線路維護體制改革探索與實踐
基于GARCH的短時風速預測方法
無人機在鐵路工程建設中的應用與思考
主站蜘蛛池模板: 亚洲精品图区| av在线无码浏览| a毛片在线| 久久99精品久久久久久不卡| 国产亚洲欧美在线视频| 凹凸国产分类在线观看| 亚洲成人在线免费观看| 五月婷婷中文字幕| 久久狠狠色噜噜狠狠狠狠97视色 | 欧美国产精品不卡在线观看| 久久久久人妻一区精品色奶水| 香蕉视频国产精品人| 综合色区亚洲熟妇在线| 亚洲αv毛片| 国产高潮流白浆视频| 啪啪免费视频一区二区| 亚洲人妖在线| 黄色福利在线| 国产日本欧美亚洲精品视| 国内黄色精品| 国产精品福利尤物youwu| 中文字幕永久视频| 性视频久久| 91在线视频福利| 国产亚洲精品无码专| 美女国内精品自产拍在线播放| 色一情一乱一伦一区二区三区小说| 亚洲无码高清视频在线观看| 欧美日韩成人| 亚洲黄色高清| 婷婷色狠狠干| 亚洲不卡影院| 亚洲精品午夜无码电影网| 丰满的少妇人妻无码区| 国产毛片片精品天天看视频| 精品无码人妻一区二区| 亚洲精品在线影院| 中文字幕久久波多野结衣| 看看一级毛片| 美女无遮挡免费视频网站| 国产日本一线在线观看免费| 国产粉嫩粉嫩的18在线播放91 | 欧美日韩第三页| 白丝美女办公室高潮喷水视频| 18禁黄无遮挡免费动漫网站| 91久久精品国产| 中文字幕66页| 国产免费久久精品99re不卡| 免费看av在线网站网址| 亚洲国产成人麻豆精品| 国产成人亚洲毛片| 热伊人99re久久精品最新地| 国产肉感大码AV无码| 就去色综合| 久久一本日韩精品中文字幕屁孩| 国产日韩欧美在线视频免费观看| 动漫精品啪啪一区二区三区| 国产福利在线免费观看| 国产网站在线看| 99精品一区二区免费视频| 波多野吉衣一区二区三区av| 国产在线拍偷自揄观看视频网站| 日本国产精品| 国产精品蜜臀| 日本不卡视频在线| 国产午夜精品鲁丝片| 午夜啪啪福利| 毛片免费视频| 一区二区日韩国产精久久| 一区二区三区高清视频国产女人| 九九免费观看全部免费视频| 亚洲欧美成aⅴ人在线观看| 欧美精品一二三区| 三级国产在线观看| 九月婷婷亚洲综合在线| 亚洲国产成人精品无码区性色| 色香蕉网站| 青草视频网站在线观看| 日本精品视频一区二区| 91久久偷偷做嫩草影院电| 中文字幕在线播放不卡| 国产精品亚洲一区二区三区z|