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

基于VOF-DPM的橫向射流霧化過程數值模擬

2024-01-01 00:00:00周濤濤唐志全王辰張郁
內燃機與動力裝置 2024年4期

摘 要:為高效精確模擬航空動力系統的液態燃料橫向射流多尺度霧化過程,分別采用離散相模型(discrete phase model, DPM)、流體體積(volume of fluid, VOF)法耦合DPM(VOF-DPM)對橫向射流霧化過程進行數值模擬,對比2種模型對橫向射流霧化過程的仿真結果,并研究模型轉換直徑與破碎模型對橫向射流霧化過程仿真結果的影響。仿真結果表明:相比DPM,VOF-DPM仿真得到的射流穿透深度更接近試驗結果,射流霧化過程更真實,并且能夠捕捉到更詳細的流場信息;當模型轉換直徑較小時,不能轉換為離散相顆粒的液滴相對較多,這些液滴仍由VOF求解,并阻擋氣流導致在其周圍產生小渦團;添加破碎模型對射流穿透深度和流場結構幾乎沒有影響,但導致離散相顆粒繼續破碎成更多更小的顆粒。

關鍵詞:橫向射流;霧化;VOF-DPM;模型轉換直徑

中圖分類號:TK401文獻標志碼:A文章編號:1673-6397(2024)04-0001-10

引用格式:周濤濤,唐志全,王辰,等.基于VOF-DPM的橫向射流霧化過程數值模擬[J].內燃機與動力裝置,2024,41(4):1-10.

ZHOU Taotao,TANG Zhiquan,WANG Chen,et al." Numerical simulation of the atomization process of liquid jet in a crossflow based on VOF-DPM[J].Internal Combustion Engine amp; Powerplant, 2024,41(4):1-10.

0 引言

橫向射流是一種典型的霧化方式,廣泛應用于航空、導彈動力系統中[1]。橫向射流中液體燃料的霧化效果是影響燃燒質量的重要因素,霧化效果越好,燃料與空氣混合越充分,越有助于提高燃燒效率和減少污染物排放[2]。對橫向射流霧化的研究主要分為試驗研究和數值模擬研究,由于霧化過程極快且物理過程非常復雜,雖然通過試驗可以得到一些宏觀參數,但難以測量流場細節;數值模擬研究能夠得到詳細的流場信息,成為研究霧化過程的重要手段[3]。

拉格朗日粒子追蹤法和界面捕捉法是用于霧化過程數值求解的2種典型方法。拉格朗日粒子追蹤法主要通過離散相模型(discrete phase model, DPM)實現,在DPM中將液相視為離散的顆粒,追蹤計算顆粒的速度、位置、動量等信息,并通過與液滴破碎模型、聚合模型等子模型的結合實現對霧化過程的模擬[4- 5]。該方法相對成熟且計算量較小,廣泛應用于射流霧化過程的模擬[6-7]。雖然DPM求解射流霧化過程計算量較小,但氣-液兩相的相互作用僅體現在歐拉方程的源項,不能真實還原霧化的全過程,尤其是忽略了一次霧化過程的具體信息。此外,早期研究指出采用DPM時液滴直徑需小于網格尺寸和Kolmogorov尺寸[8],近期也有研究提出液滴直徑小于網格尺寸的60%時才是合理的[9],這在一定程度上限制了拉格朗日粒子追蹤法的使用。界面捕捉法通過解析相界面獲得氣-液兩相界面的位置信息實現對射流霧化過程的直接求解,常用的界面捕捉方法有流體體積(volume of fluid, VOF)法[10-11]以及水平集(level set,LS)方法[12-13]。陳連華[14]采用界面捕捉法模擬了橫向射流的霧化過程,得到了詳細的霧化過程信息,研究了橫向射流一次霧化的機理,并進一步分析了不同無量綱參數對射流霧化過程的影響規律。界面捕捉法雖然能夠獲得詳細的霧化過程信息,但對網格要求較高,應足夠精細以捕捉霧化過程產生的液滴,這導致巨大的計算量。

鑒于2種方法的特點,有研究者研究結合VOF與DPM兩者優勢的新方法——流體體積法耦合離散相模型(VOF-DPM)方法[15-17]。該方法對連續的液柱、大液滴采用VOF求解,對尺寸較小的液滴通過DPM求解,既能在霧化初始階段捕捉到相界面的運動和變形,又根據一定的轉換機制將達到一定尺寸的小液滴轉換為拉格朗日粒子求解,從而減少計算量。Zhou等[18]采用VOF-DPM對不同入射角度的射流霧化進行數值模擬,并且對比了初始網格數量一致,分別采用VOF-DPM和VOF計算,自適應加密后的網格數量,僅計算48 μs后,采用VOF-DPM的網格數量比VOF的減少39.1%,表明采用VOF-DPM相比VOF可以顯著減少計算量。

合理的轉換機制是更好發揮VOF-DPM優勢的關鍵。Grosshans等[19]在研究VOF-DPM耦合模型的轉換機制時,提出了一種統計耦合方法,在靠近噴嘴處的濃密噴霧區域采用VOF求解,在下游稀疏噴霧區域采用DPM求解,在VOF中統計得到DPM中的液滴分布參數。統計耦合方法的噴霧濃密區與稀疏區的區分具有較大的隨意性,因此轉換機制更多選擇直接耦合方法,即在霧化過程模擬中,隨時識別流場中的小液滴,將滿足轉換機制的液滴轉換為拉格朗日粒子由DPM求解。直接耦合方法中,模型轉換機制通常為滿足模型轉換直徑條件和球形度條件。采用VOF-DPM對霧化過程進行數值模擬時,丁群[17]將轉換直徑設為網格尺寸的4倍,Zhou等[18]將轉換直徑設為0.1 mm,Shi等[20]將轉換直徑設為0.3 mm??梢?,不同研究人員設定的轉換直徑有一些差異,因此需要進一步研究轉換直徑對數值模擬結果的影響,以助于更合理地設定轉換直徑。

Mei等[21]采用VOF-DPM對橫向射流霧化進行了數值模擬,研究了添加液滴破碎模型的影響,研究表明添加泰勒類比破碎模型(Taylor analogy breakup model,TAB)對計算精度有負面影響;楊帆[22]采用VOF-DPM對燃氣渦輪發動機中低液氣密度比的橫向射流進行模擬時,研究了有無液滴破碎模型以及不同液滴破碎模型對模擬結果的影響,結果表明僅在氣-液界面生成液滴、沒有明顯二次霧化現象時,不宜添加液滴破碎模型。

本文中對噴嘴附近區域的液滴霧化過程進行研究,在靠近噴嘴附近以一次霧化為主,是否需要添加液滴破碎模型還需要進一步探究,因此本文中以Stenzler等[23]的試驗結果為參照,采用DPM與VOF-DPM分別進行數值模擬,對比DPM與VOF-DPM的區別,并基于VOF-DPM研究轉換直徑與破碎模型對模擬結果的影響。

2 計算域與算例設置

本文中數值模擬采用的幾何模型基于Stenzler[23]的試驗裝置進行外形尺寸的優化,計算域幾何模型如圖1所示。VOF-DPM中,幾何模型有實際射流入口,即底部的圓柱,液體燃料從底部圓柱口射入;DPM中,幾何模型中不需要實際的圓柱射流入口,在圓心處設置模型的射流入口即可。圓柱直徑d=0.254 mm,長L=1.270 mm,圓柱圓心在x方向上距離入口3 mm,距離出口10 mm,在y方向上居中。坐標系原點位于圓柱圓心處,氣體從左側流入,右側流出,其余面均為壁面;射流入口與空氣入口均為速度入口邊界,出口為壓力出口邊界,壁面邊界條件為不可滑移邊界條件。

采用ICEM軟件進行結構化網格劃分,用于VOF-DPM帶射流圓柱的幾何模型網格數為126萬,用于DPM的幾何模型網格數為128萬,除射流圓柱壁面外,其余的網格初始尺寸均為0.1 mm,網格橫截面圖如圖2所示。

進行VOF-DPM計算時開啟網格自適應加密功能,根據液相體積分數曲率進行網格細化與粗化。對液相體積分數曲率進行歸一化處理,最大值為1。細化閾值設為0.05,粗化閾值設為0.01,即在相鄰網格液體體積分數曲率大于0.05時,該處網格自動加密;在相鄰網格液體體積分數曲率小于0.01時,加密后的網格自動粗化。射流流體為水,水和空氣的物理性能參數如表1所示。

為對比VOF-DPM與DPM的區別,基于VOF-DPM研究模型轉換直徑以及破碎模型對模擬結果的影響,設置4個方案,各方案參數設置如表2所示。4個方案中,射流入口速度均為16.99 m/s,氣流入口速度均為116.00 m/s。

為了保證獲得的結果不受網格質量的影響,進行網格無關性檢驗。選用方案 2中的多相流模型及參數,分別在無網格自適應加密和2、3、4級網格自適應加密條件下進行計算,以中心對稱面x/d=-1.5位置處的氣流速度計算結果進行對比,結果如圖3所示,圖中z為笛卡爾坐標系z方向上的距離。由圖3可知:無網格自適應加密和添加網格自適應加密的結果相差較大,添加各級自適應加密的結果相差不大;z/d=4~9時,添加2級網格自適應加密的結果比3級小, 3、4級網格自適應加密的計算結果幾乎一致,說明3級自適應加密達到了計算要求,滿足了網格的無關性。因此,本文中VOF-DPM中的數值計算均采用3級網格自適應加密。

3 結果分析與討論

3.1 瞬時霧化結果

第0.45毫秒時,不同方案α=0.5相界面與顆粒分布圖像如圖4所示,紅色表示相界面,小圓點表示離散相顆粒。

由圖4a)可知:采用DPM得到的射流噴霧全部由離散相顆粒組成,在迎風面一側顆粒分布濃密且粒徑較大,在背風面一側顆粒稀疏且粒徑較?。槐緫沁B續的射流柱全部由顆粒表示,無法從模擬結果中觀察到射流柱的變形破裂及表面波現象,只能得到霧化后液滴的速度、大小和射流深度等特征信息。

由圖4b)、4c)、4d)可知:采用VOF-DPM可以更真實地再現射流霧化過程;當液體射流進入橫向氣流區時,液柱由圓柱狀變為扁平狀并逐漸彎曲,在氣動力的作用下在近壁面發生表面破碎,最終在較遠處發生柱狀破碎;同時,可觀察到射流液柱上的表面波,這有助于分析霧化的機理等。

由圖4b)、4c)可知:方案2中有更多的小液滴,這是因為方案 2、 3的模型轉換直徑不同;方案3中,滿足球形度條件的液滴的體積等效直徑小于0.25 mm時,便能夠轉換成離散相顆粒,方案2中滿足球形度條件的液滴的體積等效直徑小于0.15 mm時才轉換成離散相顆粒,即滿足球形度的體積等效直徑為0.15~0.25 mm的液滴在方案3中轉換成離散相顆粒,在方案2中依然還是由相界面表示的液滴,導致方案 2中的小液滴更多。由此可知,當模型轉換直徑設置得較小時,計算過程中更多的小液滴由VOF求解,但同時也帶來更大的計算量。

由圖4b)、4d)可知:在VOF-DPM模擬中添加破碎模型得到的結果中離散相顆粒數量更多。這是因為添加KH-RT破碎模型后,當滿足破碎模型中的破碎條件,通過VOF-DPM轉換生成的離散相顆粒繼續分裂成更多的離散相顆粒,在沒有添加破碎模型的方案2中,小液滴在滿足轉換機制轉換成離散相顆粒后,不繼續分裂,所以產生的離散相顆粒數量相應較少。

3.2 射流穿透深度

Stenzler[23]通過試驗總結出了不同條件橫向射流穿透深度關聯式,本文中的氣體為常溫,故選用適用于未加熱氣流中的射流穿透深度關聯式與模擬結果對比,關聯式如下:

zdmax=2.898q0.430xd0.384We-0.110,(14)

式中:q為動量通量比,q=ρlu2l/(ρgu2g);We為韋伯數。

數值計算的穿透深度結果以離散相顆粒位置分布圖表示, 4個方案模擬得到的射流穿透深度與試驗對比圖5所示。由圖5可知:采用DPM的方案 1得到的射流穿透深度比試驗結果偏大,方案 2、3、4得到的射流穿透深度與試驗結果十分接近;采用VOF-DPM的3個方案的射流穿透深度都是x方向上前半段略小于試驗結果,在后半段則幾乎與試驗結果重合并有大于試驗結果的趨勢,說明相比DPM,VOF-DPM可以得到與試驗結果更接近的射流穿透深度。3個方案得到的射流穿透深度基本沒有區別,說明模型轉換直徑和破碎模型對射流穿透深度幾乎沒有影響。這是因為不論液滴是否轉換為離散相顆?;蛘哳w粒破碎成更小的顆粒,液滴或者顆粒的動量都被保留繼承,因此基本不影響計算結果中的射流穿透深度。

3.3 液滴平均直徑

霧化后液滴大小是衡量霧化效果的關鍵因素。本文中采用索特平均直徑(SMD)表示霧化后的平均液滴大小。索特平均直徑

d32=∑(Nid3i)/∑(Nid2i),(15)

式中:Ni為直徑為di顆粒的數量。

由式(15)可得:方案1~4計算得到的索特平均直徑分別為60.87、71.99、77.31、65.53 μm。4個方案的結果相差不大, DPM計算得到的SMD最小。轉換直徑較大時,SMD較大。這可能是因為相比方案3,方案2中滿足球面度標準且體積等效直徑為0.15~0.25 mm的液滴仍由VOF求解,未能轉換為離散相顆粒,統計方法中只統計離散相顆粒,導致方案2中直徑較大的液滴未納入統計范圍,使得方案2中統計的SMD相對偏??;對比方案4和方案2、3的結果可知,添加破碎模型的方案4的SMD小于其他2個方案,這是因為添加破碎模型后,離散相顆粒滿足破碎條件,繼續分裂成更多的直徑更小的顆粒。

3.4 流場結構

在橫向射流中,橫向氣流遇到液體射流柱的阻擋改變流動方向,進而產生各種渦旋流場,這些漩渦流場又影響霧化液滴的分布混合。本文采用Q準則來識別流場中的渦,Q的定義為:

Q=Ω2-S2/2,

式中Ω和S分別是旋轉張量和應變張量。Q>0,表示區域中存在渦旋結構[23],Q準則的物理意義是流場旋轉部分的渦度大于變形而占主導地位。

4個方案在Q=7.5×107的渦量等值面與α=0.5相界面(以紅色表示)分布圖如圖6~9所示。由圖6~9可知:橫向氣流在射流液柱后產生劇烈的擾動,形成強烈的渦旋結構;DPM與VOF-DPM得到的渦量分布差異較大,DPM得到的渦量分布窄小且對稱,并且在后半部分的渦尺度較小,VOF-DPM由于計算得到射流分布矮小寬闊,得到的渦量分布也有同樣的特點;VOF-DPM得到的渦量在射流柱后底部有在y方向上相互交替的小尺度雙列渦,在中間有在y方向上相互交替的大尺度雙列渦,最上面是由于液團/液滴的阻擋形成的小尺度渦團。產生相互交替的雙列渦的原因可能是橫向氣流在遇到射流柱的阻礙后,產生了卡門渦街現象。4個方案的結果中都可以觀察到射流柱前產生了渦團,這是氣流受到阻礙后形成的馬蹄渦,VOF-DPM得到的馬蹄渦結構更詳細。因為方案 2、4中的轉換直徑更小,有更多的液滴未能轉換到離散相顆粒計算,因此產生了更多的小渦團。

方案1、3在x方向上x/d為2、5、10處橫截面下半平面的流線分布如圖10所示。由圖10可知:2個方案的流線分布差異較大;在x/d為2時,方案 1流線圖在靠近壁面處有一個反向渦旋對(counter-rotating vortex pair,CVP);x方向上,由于射流深度增加,CVP逐漸增大;不同x/d截面上,方案 3的渦流產生在背風面(低速區)的不同地方,并且沒有像方案1中那樣產生CVP,結合圖6~9分析,這是因為在方案3中射流柱后形成的是相互交替的雙列渦; VOF-DPM的背風面速度比DPM小的多,這也說明VOF-DPM中氣液剪切更加劇烈。

4 結論

分別采用DPM與VOF-DPM對橫向射流的霧化過程進行數值模擬,并基于VOF-DPM對比了模型轉換直徑以及破碎模型對模擬結果的影響。

1)相比DPM,VOF-DPM可以得到更真實的射流霧化過程、更準確的射流穿透深度,并能捕捉到更詳細的渦旋結構以及霧化后液滴周邊的小渦團,有利于下一步更準確地模擬霧化后液滴的蒸發過程。

2)當模型轉換直徑較小時,有相對較多的液滴不能滿足轉換機制,無法轉換為離散相顆粒,仍由VOF解析,相應的流場中產生更多的由液滴引起的小渦團;在可選擇的范圍內,轉換直徑設置越小,VOF-DPM求解中有更多的液相由VOF求解,對細節的捕捉更加詳細,但計算量的需求更大,因此應根據需要選擇合適的模型轉換直徑。

3)在靠近噴嘴的噴霧場采用VOF-DPM對霧化過程進行模擬,是否添加二次破碎模型對射流穿透深度以及流場結構幾乎沒有影響,但是添加破碎模型產生更多的離散相顆粒,得到的液滴SMD更小。

參考文獻:

[1] 常建龍,陳連華,趙永娟,等.橫向射流液滴霧化研究現狀分析[J].戰術導彈技術, 2022(2):29-36.

[2] LI X Y, SOTEROIU M, ARIENTI M, et al. High-fidelity simulation of atomization and evaporation in a liquid jet in cross-flow[C]//Proceedings of AIAA.49th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition. Reston,USA:AIAA,2011.

[3] 李佩波.超聲速氣流中橫向噴霧的混合及燃燒過程數值模擬[D].長沙:國防科技大學, 2019.

[4] DUKOWICZ J K. A particle-fluid numerical model for liquid sprays[J]. Journal of Computational Physics,1980,35(2):229-253.

[5] TSUJI Y, KAWAGUCHI T, TANAKA T. Discrete particle simulation of two-dimensional fluidized bed[J].Powder Technology, 1993,77(1):79-87.

[6] YOO Y L, HAN D H, HONG J S, et al. A large eddy simulation of the breakup and atomization of a liquid jet into a cross turbulent flow at various spray conditions[J].International Journal of Heat and Mass Transfer, 2017,112:97-112.

[7] IRANNEJAD A, JABERI F. Large eddy simulation of turbulent spray breakup and evaporation[J].International Journal of Multiphase Flow, 2014,61:108-128.

[8] FOX R O. Computational models for turbulent reacting flows[M].Cambridge ,UK:Cambridge University Press, 2003.

[9] ABOUELMAGD A, THVENIN D.Direct numerical simulation of spray evaporation and autoignition in a temporally-evolving jet[J].Proceedings of the Combustion Institute, 2017,36(2):2493-2502.

[10] SOH G Y, YEOH G H, TIMCHENKO V.Improved volume-of-fluid (VOF) model for predictions of velocity fields and droplet lengths in microchannels[J].Flow Measurement and Instrumentation, 2016,51:105-115.

[11] GROSSHANS H, MOVAGHAR A, CAO L, et al.Sensitivity of VOF simulations of the liquid jet breakup to physical and numerical parameters[J].Computers amp; Fluids, 2016,136:312-323.

[12] SHAO C X, YUAN S A, LUO K.A generalized coupled level set/volume-of-fluid/ghost fluid method for detailed simulation of gas-liquid flows[J].Journal of Computational Physics, 2023,487:112158.

[13] MNARD T, TANGUY S, BERLEMONT A.Coupling level set/VOF/ghost fluid methods: validation and application to 3D simulation of the primary break-up of a liquid jet[J].International Journal of Multiphase Flow, 2007,33(5):510-524.

[14] 陳連華.亞聲速氣流中液體橫向射流霧化數值模擬研究[D].太原:中北大學, 2022.

[15] STRM H, SASIC S, HOLM-CHRISTENSEN O, et al.Atomizing industrial gas-liquid flows:development of an efficient hybrid VOF-LPT numerical framework[J].International Journal of Heat and Fluid Flow, 2016,62:104-113.

[16] BHATIA B, JOHNY T, DE A.Understanding the liquid jet break-up in various regimes at elevated pressure using a compressible VOF-LPT coupled framework[J]. International Journal of Multiphase Flow, 2023,159:104303.

[17] 丁群. 基于OpenFOAM的VOF-DDM模型對霧化現象的研究[D].哈爾濱:哈爾濱工程大學, 2018.

[18] ZHOU W Y, CHEN B, ZHU Q B,et al.Numerical simulation of angled-injected liquid jet breakup in supersonic crossflow by a hybrid VOF-LPT method[J].International Journal of Multiphase Flow, 2023,166:104503.

[19] GROSSHANS H, SZSZ R Z, FUCHS L.Development of an efficient statistical volumes of fluid-Lagrangian particle tracking coupling method[J].International Journal for Numerical Methods in Fluids, 2014,74(12):898-918.

[20] SHI P, ZHU G Q, CHENG J M, et al.Simulation on atomization process of gas-liquid pintle injector in LRE under periodic conditions based on the VOF to DPM method[J].Aerospace Science and Technology, 2023,136:108222.

[21] MEI Y,ZHANG P,YAN Y W.Numerical simulation on spray characteristics of fuel jet in a crossflow[J].Transactions of Nanjing University of Aeronautics and Astronautics, 2020,37(Suppl.):18-27.

[22] 楊帆.低液氣密度比橫向射流一次霧化數值模擬研究[D].合肥:中國科學技術大學, 2021.

[23] STENZLER J, LEE J, SANTAVICCA D.Penetration of liquid jets in a crossflow[C]//Proceedings of 41st Aerospace Sciences Meeting and Exhibit.Reno,USA: AIAA,2003.

[24] WU P K, TSENG L K, FAETH G.Primary breakup in gas/liquid mixing layers for turbulent liquids[M]//Proceedings of 30th Aerospace Sciences Meeting and Exhibit. Reno,USA: AIAA,1995.

[25] BEALE J C, REITZ R.Modeling spray atomization with the Kelvin-Helmholtz/rayleigh hybrid model[J].Atomization and Sprays,1999,9(6):623-650.

Numerical simulation of the atomization process of liquid jet in

a crossflow based on VOF-DPM

ZHOU Taotao, TANG Zhiquan, WANG Chen, ZHANG Yu*

School of Automotive and Transportation Engineering,Hefei University of Technology, Hefei 230009,China

Abstract: In order to efficiently and accurately simulate the multi-scale atomization process of liquid jet in crossflow(LJIC) in aerospace power system, the LJIC atomization process is simulated using the discrete phase model (DPM) and coupled the volume of fluid (VOF) and DPM (VOF-DPM). The simulation results of the two models on the atomization process of LJIC are compared, and the effects of model conversion diameter and breakup model on the simulation results are investigated. The simulation results indicate that,compared with DPM, the jet penetration depth obtained by VOF-DPM is closer to the experimental results,the atomization process is more realistic and more detailed flow field information can be captured. When the model conversion diameter is smaller, there are more droplets that cannot be converted into discrete phase particles, which are still solved by VOF and block the flow causing small vortex clusters to form around them. Breakup model has little effect on the jet penetration depth and the flow field structure, but it causes the discrete phase particles to be further breakup into more and smaller particles.

Keywords: liquid jet in a crossflow; atomization; VOF-DPM; model transition diameter

(責任編輯:劉麗君)

收稿日期:2023-01-19

基金項目:國家自然科學基金青年基金資助項目(51906055,52106138);安徽省自然科學基金資助項目(2008085ME166);合肥工業大學學術新人提升B計劃資助項目(JZ2021HGTB0085)

第一作者簡介:周濤濤(1992—),男,山東菏澤人,工學博士,副教授,主要研究方向為多相流數值模擬,E-mail: zhoutt@hfut.edu.cn。

*通信作者簡介:張郁(1989—),男,河南洛陽人,工學博士,副研究員,主要研究方向為燃料液滴蒸發和燃燒,E-mail:yuzhang2020@sina.com。

主站蜘蛛池模板: 99草精品视频| 亚洲无线一二三四区男男| 91丝袜在线观看| 91在线播放国产| 伊人大杳蕉中文无码| 国产无码网站在线观看| 亚洲第一网站男人都懂| 大陆精大陆国产国语精品1024| 久久久精品国产SM调教网站| 国产在线拍偷自揄观看视频网站| 在线观看国产黄色| 亚洲视频色图| 日日拍夜夜操| 在线欧美一区| 精品视频福利| 在线日本国产成人免费的| 免费精品一区二区h| 亚洲成aⅴ人在线观看| 中文字幕久久波多野结衣 | a在线亚洲男人的天堂试看| 欧美成人精品在线| 日韩美毛片| 久久久久亚洲精品成人网| 国产精品无码一区二区桃花视频| 久久久黄色片| 色婷婷在线播放| 国产电话自拍伊人| 91网址在线播放| 97青草最新免费精品视频| 国产福利在线免费观看| 亚洲国产成人超福利久久精品| 国产免费一级精品视频| 国产99欧美精品久久精品久久| 国模视频一区二区| 亚洲国产精品日韩欧美一区| 在线观看精品自拍视频| 女同国产精品一区二区| 亚洲h视频在线| 亚洲欧美激情另类| 欧美亚洲欧美| 亚洲男人天堂久久| 久热99这里只有精品视频6| 区国产精品搜索视频| 亚洲第一视频区| 日韩精品一区二区三区视频免费看| 国产精品成人啪精品视频| 伊人久久久久久久| 国产精品污污在线观看网站| 在线色国产| 色视频国产| 又爽又黄又无遮挡网站| 网友自拍视频精品区| 亚洲成人播放| 日韩精品中文字幕一区三区| 国模视频一区二区| 亚洲精品你懂的| 综1合AV在线播放| 国产第一页屁屁影院| 国内精品免费| 国产91高跟丝袜| 欧美日韩另类在线| 亚洲色图另类| 亚洲人成在线精品| 自慰高潮喷白浆在线观看| 福利视频一区| 久久国产高清视频| 亚洲无线国产观看| 四虎国产精品永久一区| 国产在线无码av完整版在线观看| 四虎国产精品永久一区| 黄色网页在线观看| 99视频在线观看免费| 综合久久久久久久综合网| 日本一区二区三区精品AⅤ| 激情六月丁香婷婷| 亚洲欧美综合在线观看| 亚洲天堂日本| 国产裸舞福利在线视频合集| 亚洲一区无码在线| 日本福利视频网站| 久久午夜夜伦鲁鲁片不卡| 日韩无码精品人妻|