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

正十二烷噴霧火焰大渦模擬研究

2023-05-26 08:06:12趙萬輝衛海橋
內燃機學報 2023年3期

趙萬輝,孫 婷,衛海橋,周 磊

(1. 中國民航大學 航空工程學院,天津 300300;2. 天津大學 內燃機燃燒學國家重點實驗室,天津 300350)

近年來,國內外各研究機構基于定容燃燒彈開展了大量的正十二烷噴霧火焰試驗.Gimeno等[1]對于正十二烷冷態噴霧及噴霧火焰的研究相對較早,他們開展了不同運行參數(環境溫度和密度等)對噴霧發展過程的影響研究.Sandia國家實驗室也進行了大量的試驗,并將標準化的試驗結果(氣相、液相噴霧貫穿距和火焰浮起長度等)公布到Engine Combustion Network(ECN)[2]網站上,方便全世界的研究學者訪問.國內方面,倪兆靜等[3]探究了液滴蒸發過程.賀鵬飛等[4]在試驗中探究了噴霧火焰燃燒和碳煙顆粒的形成,發現低氧體積分數條件下碳煙顆粒氧化過程明顯受到抑制,顆粒平均直徑和團聚物回轉半徑明顯增大;而碳煙的生成區域還會受到燃油溫度的影響.之前的研究有助于理解噴霧火焰中涉及的復雜的物理化學現象,并為高效數值仿真計算平臺的發展提供了試驗基礎.然而受限于拍攝精度、拍攝頻率等因素,試驗手段獲得的結果仍非常有限.計算流體力學方法有助于進一步理解噴霧火焰的時間和空間發展過程.

雷諾平均(RANS)和大渦模擬(LES)方法廣泛應用于噴霧火焰的模擬中.與RANS方法相比,LES方法能夠以可接受的計算量捕捉到更細的湍流渦結構、褶皺火焰結構、噴霧火焰中局部多點和多階段著火過程以及火焰穩定等現象,在模擬噴霧火焰燃燒過程方面極具潛力[5-7].噴霧火焰中的多階段著火現象可以利用關鍵組分(如CH2O、H2O2和OH等)的產生和消耗來進行區分[8].在噴霧中出現高溫著火之前,燃料分子發生脫氫裂解反應,隨后經歷加氧、分子異構化等過程后產生酮類過氧化物RO2,并伴隨著CH2O、H2O2等組分的快速累積,表明噴霧中發生一階段反應.受混合氣活性和初始工況的影響,在低溫條件下,一階段著火可以從稀薄混合氣區域過渡到濃混合氣區域,而在低氧體積分數、高溫或高密度條件下,油、氣混合受到限制,一階段著火只能出現在濃混合氣區域.由于大渦模擬方法需要采用更加精細的網格,對計算資源的需求顯著增大.自適應網格技術只需對局部網格進行加密,在一定程度上減少了計算量,然而總的計算量仍然非常巨大.另外,模擬湍流燃燒過程時需要采用合適的湍流燃燒模型對能量方程和組分方程封閉,這些模型進一步增加了計算量,并提高了湍流燃燒模擬過程的復雜性,導致大渦模擬方法的應用仍存在一定局限性,利用大渦模擬方法研究復雜幾何形狀燃燒室燃燒過程仍存在巨大挑戰.

正十二烷具有與柴油更加接近的碳鏈長度和沸點,以正十二烷為替代燃料能夠更加準確地貼近柴油的蒸發和混合過程.近年來對正十二烷噴霧火焰的模擬已然成為研究熱點[5,9].Zhang等[10]探究了湍流與化學反應相互作用對噴霧自燃和噴霧火焰結構的影響.Zhao等[11]探究了兩次噴射中噴霧與火焰之間的相互作用,發現兩次噴射的燃油與火焰存在較強的相互作用,其中預先噴射的燃油中發生的化學反應可以引起溫度升高,并產生不同中間組分,如OH和CH2O等.以上因素均會加速二次噴霧著火過程,并且初次噴霧引起的溫升發揮最主要的作用.Kaario等[12]以正十二烷噴霧為基礎,探究了甲烷氛圍下雙燃料噴霧火焰的發展過程,發現盡管初始環境中的甲烷會在不同程度上推遲噴霧著火,正十二烷噴霧火焰和雙燃料火焰中混合氣的著火過程均可以分為湍流混合、低溫燃燒、高溫燃燒和高溫擴散火焰等區域.目前國內對正十二烷噴霧火焰大渦模擬研究較少,王利民等[13]探究了不同溫度條件下正十二烷噴霧的燃燒過程,但仍缺少對低溫燃燒策略典型工況(低氧體積分數、高壓縮比(高環境密度))條件下噴霧著火和火焰發展過程的深入探究,同時對正十二烷噴霧火焰穩定過程的預測精度也需要做出進一步評價.因此,筆者旨在進一步闡明當前大渦模擬方法對寬工況范圍內正十二烷噴霧火焰的預測精度,并探究環境變量,如溫度、氧體積分數、環境密度及噴油壓力對噴霧著火和火焰穩定過程的影響規律.另外,文獻[7]中LESLEM模型被應用到正庚烷噴霧燃燒過程模擬研究中,并能夠較好地預測寬工況正庚烷噴霧混合氣形成、著火和燃燒過程.正庚烷噴霧工況與正十二烷有很大不同,筆者在該研究的基礎上展開,同時在不需要對噴霧模型和湍流燃燒模型相關參數進行修改的前提下進行模擬研究,并驗證了當前具有高階精度的仿真計算平臺——LES-LEM模型在廣泛工況下的適用性.

1 模型參數設置與驗證

1.1 參數設置及模擬工況介紹

計算域是一個高為100mm、直徑為30mm的圓柱,所采用的網格總數約為8×105.燃油從圓柱形計算域頂端中心位置向外噴射.試驗中,噴油參數的設置如表1所示;計算域中初始環境溫度、氧體積分數和環境密度等參數設置如表2所示,且與ECN網站試驗工況保持一致.表3為噴霧模型、差分格式和燃燒模型.其中壁面邊界條件為無滑移邊界條件.化學反應機理采用Yao等[14]的包含54種組分和269步化學反應的機理,該機理可以比較準確地預測噴霧著火過程,特別是在低溫、低氧體積分數條件下的模擬結果明顯優于其他機理,因而被廣泛應用于正十二烷噴霧火焰的模擬當中[15].

表1 噴油參數設置Tab.1 Parameters for spray modeling

表2 參數設置Tab.2 Numerical setting

表3 模型選擇Tab.3 Settings for numerical models

1.2 控制方程簡介

大渦模擬方法是在KIVALES模型[16]的基礎上發展而來.基本控制方程包括連續性方程、動量守恒方程和能量守恒方程等.

式中:ρ、u、t、e和分別為密度、速度、時間、內能和壓力;sgs為亞網格尺度變量;τij、和分別為黏性應力、熱通量、亞網格分子輸運和亞格子顯焓;和分別為噴霧產生的液滴阻力源項、化學反應源項及噴霧源項.為避免重復,只列出了部分控制方程,其他方程及求解過程參見文獻[7].

1.3 模型驗證

之前的研究[7]表明,當前大渦模擬方法適合模擬不同環境工況下的冷態噴霧以及噴霧火焰發展過程.對于未發生化學反應的正十二烷冷態噴霧,在不同噴油壓力條件下,當前方法得到的氣相噴霧貫穿距模擬結果與試驗結果非常接近,同時不同位置上燃油的濃度分布結果也與試驗結果非常接近[6].圖1給出了不同工況下著火延遲期和火焰浮起長度的模擬與試驗結果.著火延遲期(ID)定義為計算域中最高溫度達到著火溫度Tign對應的時刻,其中Tign定義為

圖1 不同工況下著火延遲期和火焰浮起長度結果Fig.1 Ignition delay times and flame lift-off lengths under different conditions

式中:Tamb與Tmax分別代表環境氣體的初始溫度以及噴霧火焰達到準穩態時計算域內的最高溫度[17].

火焰浮起長度(LOL)是指噴霧火焰上游位置溫度達到著火溫度Tign的點到達噴孔處的最小距離.LOL表示準穩態條件下,噴霧中高溫反應區最上游位置基本保持不變,此時燃油蒸發與燃燒過程相平衡,火焰達到準穩定狀態.從圖1中還可以看出增加氧體積分數或者提高環境密度均會導致LOL縮短.

當前大渦模擬方法整體上能夠捕捉正十二烷噴霧著火過程隨初始條件的變化趨勢,即溫度升高、氧體積分數升高或者環境密度升高都會導致著火提前、火焰浮起長度縮短.在高溫、高密度條件下著火延遲期和火焰浮起長度模擬結果與試驗值非常接近,而在不同氧體積分數條件(圖1a)下,預測結果比試驗值略低,這與該化學反應機理活性較強有關,混合氣更容易著火.而在低密度條件(圖1b)下,著火延遲期預測值比試驗值略高,該工況下環境密度比較低,混合氣活性下降.

圖2 對比了算例1基礎工況模擬得到的碳煙結果和噴霧內密度分布.相同工況下試驗結果中包含了CH2O和多環芳烴(PAH)的濃度分布,而在模擬采用的化學反應動力學機理中不包含PAH組分,因而利用C2H2來定性地比較碳煙圖像.從圖2中可以看到,C2H2主要分布在下游位置,與試驗結果中的組分分布趨勢保持一致.同時,密度分布的紋影結果來自ECN網站,試驗工況與模擬算例相同.從紋影圖中可以看出,火焰前端位置與試驗結果基本保持一致.因此,模擬結果較好地捕捉了組分分布、火焰前端位置.

圖2 在環境溫度為900K、0.7ms時刻算例1碳煙試驗和質量分數模擬結果及密度分布的紋影圖片和模擬結果對比Fig.2 Distributions of soot experiments and soot mass fraction numerical results,as well as density from schlieren images and simulations at 0.7ms,900K of Case 1

2 結果分析與討論

2.1 初始氧體積分數的影響

圖3展示了噴霧火焰發展過程.可以發現,在1.0ms時刻噴霧火焰已經達到準穩定燃燒狀態,火焰浮起長度與1.5ms時刻下的數值相差不大.之前的研究[18]表明,初始環境中氧體積分數決定著噴霧火焰中OH質量分數的最高值,即氧體積分數的升高有助于促進OH生成,并使得最高溫度明顯升高,混合氣活性增強.混合物分數體現了燃油分數分布,計算方法參見文獻[19].混合物分數為0代表氧化劑(空氣),混合物分數為1表示純燃料.

圖3 900K溫度條件下算例1火焰結構對比Fig.3 Comparison of flame structures of Case 1 at 900K

圖4 對比了不同氧體積分數條件下正十二烷噴霧火焰的著火過程.需要注意的是,ECN網站上僅給出了13%、15%和21%氧體積分數條件下的著火延遲期和火焰浮起長度數據,并未給出10%氧體積分數條件下的試驗結果.但由于發動機新型燃燒模式中多采用高EGR率,氧體積分數較低,為了保證燃料在低氧體積分數條件下能夠穩定著火,探究低氧體積分數條件下的著火過程十分必要.為此,圖4中增加了10%氧體積分數條件下的著火過程.該工況下溫度、壓力、環境密度和噴油壓力均與15%氧體積分數工況保持一致.但由于缺少著火延遲期試驗數據,因而未在圖1a中畫出.圖4中彩色散點代表H2O、H2、CO2和CO質量分數之和,可以用反應進度標量C來表示.隨著氧體積分數的升高,溫度開始升高以及高溫反應出現的時刻提前,火焰浮起長度縮短,導致濃混合氣參與高溫反應的量增多.

圖4 不同氧體積分數條件下的著火過程Fig.4 Ignition process under different oxygen volume fraction conditions

圖5顯示氧體積分數降低后噴霧火焰中反應進度緩慢,溫度升高速率明顯降低.低氧體積分數條件下火焰浮起長度變長,油、氣混合更佳,但是噴霧火焰中仍然存在著濃混合氣參與高溫反應的現象.當初始氧體積分數降低至8%時,混合氣活性較差,放熱過程趨緩,最高燃燒溫度明顯降低.但是噴霧火焰仍然呈現出準穩定燃燒狀態,火焰浮起長度(圖6中白色虛線)幾乎不隨時間變化.不同工況條件下,火焰浮起長度的變化反映了環境溫度、氧體積分數以及環境密度對火焰穩定過程的影響規律.在低溫、低氧體積分數或者低環境密度條件下,高溫反應難以在上游位置出現,火焰浮起長度更長,在這些工況下著火推遲,說明不利于噴霧火焰更早地達到穩定狀態.

圖5 不同氧體積分數條件下算例3~5火焰結構對比Fig.5 Comparison of flame structures of Cases 3—5 under different oxygen volume fraction conditions

圖6 8%氧體積分數條件下算例2火焰結構Fig.6 Flame structures of Case 2 under 8%O2 volume fraction condition

噴霧著火延遲期較短(小于噴油持續期)時,燃燒室內的燃油存在著邊蒸發邊燃燒的現象.燃油蒸發量可以反映出燃燒室內化學反應的程度和起始時刻.對于未發生化學反應的工況,燃燒室內各個單元燃油累計蒸發質量隨時間單調升高;當化學反應出現以后,燃油質量停止升高的點即代表了燃油消耗速率超過了蒸發速率的時刻.圖7所示氧體積分數降低導致著火推遲,因而噴油結束前燃油蒸發質量停止升高的時刻推遲;同時,在13%至21%氧體積分數條件下,燃油蒸發質量存在著超過1ms的穩定時刻,說明在此期間(0.2ms到1.5ms),燃油的蒸發過程與消耗過程存在著平衡關系,溫度較低的燃油蒸發之后迅速參與化學反應.

圖7 溫度為900K時不同氧體積分數下燃油蒸發質量隨時間的變化Fig.7 Temporal evolution of fuel evaporation mass at 900K under different oxygen volume fraction conditions

升高氧體積分數可以明顯提高噴霧火焰中的溫度,并有助于加快化學反應速率,使得碳煙氧化過程加快.Wang等[20]研究發現,21%氧體積分數條件下噴霧火焰中的碳煙出現時刻更早、峰值濃度更高,且分布在更細小的空間內.最終的碳煙生成量是其生成與氧化過程相互競爭的結果.高氧體積分數條件下噴霧火焰著火延遲期和火焰浮起長度均明顯縮短,有限的油、氣混合時間是造成碳煙排放較高的重要原因.圖8a顯示正十二烷噴霧燃燒過程中整個計算域內甲醛生成總量隨時間變化,甲醛質量的升高表明燃燒室內甲醛的生成速率大于消耗速率.10%氧體積分數條件下的著火延遲期為0.74ms,而甲醛質量開始下降的時刻大約發生在1.1ms時刻以后,此時噴霧內部出現大量的高溫反應,導致火焰內部溫度快速升高,使得化學反應速率加快,在1.1ms時刻以后,甲醛的消耗速率明顯高于生成速率,因而甲醛質量開始下降.圖8b顯示了1.5ms(噴油結束)時刻后碳煙質量與初始氧體積分數呈現出非單調變化趨勢,這主要是與噴油結束前燃油消耗的量、噴霧火焰與空氣的接觸面積有關.Chishty等[15]基于輸運概率密度函數探究正十二烷噴霧火焰碳煙生成過程時也發現,碳煙質量與氧體積分數的變化呈非單調變化趨勢.由于氧體積分數對碳煙生成的影響更加復雜,氧體積分數的變化不僅改變了混合氣的濃度,火焰的溫度和碳煙的氧化速率均隨著氧體積分數的升高而有所提升.混合氣活性的提高導致著火提前,隨后大量的濃混合氣參與高溫反應從而產生大量碳煙.與此同時,火焰溫度升高加快碳煙氧化速率.整體上碳煙的生成過程與氧化過程的競爭關系決定了最終的碳煙排放水平.另外,甲醛質量隨著初始氧體積分數的變化呈現出單調變化趨勢,并與碳煙的排放水平形成鮮明對比.對于13%、15%和21%氧體積分數條件下,著火延遲期相差不大,均小于0.5ms.盡管降低氧體積分數后導致混合氣活性下降,由于此時混合氣活性仍然較強,一階段著火出現時刻仍然比較早,造成CH2O大量累積.當氧體積分數降至10%時,著火過程明顯推遲,一階段著火過程受到抑制.

圖8 溫度為900K時不同氧體積分數下甲醛和碳煙質量隨時間的變化Fig.8 Temporal evolution of CH2O and soot mass at 900K under different oxygen volume fraction conditions

2.2 初始環境密度和噴油壓力的影響

高環境密度(或壓力)條件下,分子間的碰撞頻率提高,混合氣活性增強.冷態的燃油液滴噴入到高溫、高壓和高密度環境中,迅速蒸發形成可燃混合氣,并以較高速度向下游運動.圖9顯示了整個燃燒室內早期蒸發的燃油質量受環境密度的影響較小.高溫反應是造成燃油快速消耗的重要因素,并且在噴油結束前,燃油的消耗與蒸發基本保持平衡,蒸發的燃油迅速參與高溫反應,分子鏈斷裂產生大量的OH基團.在0.2ms時刻左右,高密度條件下燃油累積質量升高趨勢有所放緩;而在低密度條件下曲線的升高速率放緩的時刻略有推遲.圖10顯示了低環境密度(15.2kg/m3)條件下的火焰發展情況,高環境密度(22.8kg/m3)下的火焰發展參考圖3.對比圖3和圖10可以發現,低密度條件下噴霧火焰高溫反應區明顯位于下游位置,火焰浮起長度變長.

圖9 不同環境密度下算例1和6燃油蒸發質量隨時間的變化Fig.9 Temporal evolution of fuel evaporation mass under different ambient density conditions of Cases 1 and 6

圖10 低環境密度為15.2kg/m3工況下算例6火焰結構Fig.10 Flame structures of Case 6 under a low ambient density of 15.2kg/m3 condition

不同環境密度條件下正十二烷噴霧火焰著火過程如圖11所示.高溫著火位置首先出現在噴霧火焰的前端溫度和混合氣分數均比較合適的區域,高溫著火出現后,高溫反應區快速向上游方向發展,并逐漸形成穩定結構.此時由于火焰浮起長度縮短,濃混合氣開始參與高溫反應,并引起溫度升高.

圖11 不同環境密度條件下溫度-混合物分數散點Fig.11 Scatters of temperature-mixture fraction under different ambient density conditions

采用高壓噴射有助于改善發動機經濟性、熱效率和排放特性,因而高壓噴射受到了廣泛的關注.噴油壓力過高導致油束向下游運動速度過快,極易引起噴霧撞壁,造成發動機效率下降和排放變差.Liu等[21]研究表明,高壓噴射可以引起油、氣過度混合,從而導致著火失敗.因而有必要深入探究不同噴油壓力條件下噴霧著火過程,試驗工況設置如表2中算例8~10所示.圖12顯示不同噴油壓力下算例7~9的燃油噴射速率.筆者深入開展了初始環境溫度為800K、氧體積分數為21%條件下正十二烷噴霧火焰著火過程大渦模擬.

圖12 不同噴油壓力下算例7~9燃油噴射速率Fig.12 Rate of injection of Cases 7—9 with different injection pressures

圖13 對比了不同噴油壓力下的火焰著火延遲期(ID)和火焰浮起長度.可以發現,當前大渦模擬方法對著火過程的預測值偏低,而火焰浮起長度模擬結果與試驗值非常接近,特別是在噴油壓力為50MPa時,模擬結果與試驗值偏差的最小值約為5%.大渦模擬方法可以捕捉噴霧火焰瞬態燃燒特征,如早期噴霧尾端的自燃點與主火焰合并的瞬態過程,更容易造成對著火延遲期的預測偏低,綜合來看偏差在20%以內均是合理的[22].總之,當前模型可以捕捉噴油壓力對著火延遲期和火焰浮起長度的影響規律,即升高噴油壓力對著火過程的影響較小,卻能延長火焰浮起長度.

圖13 不同噴油壓力下算例8~10著火延遲期和火焰浮起長度Fig.13 Ignition delay and flame lift-off lengths of Cases 8—10 with different injection pressures

Chishty等[15]對比了不同初始環境溫度條件下的著火延遲期和火焰浮起長度,并探究了不同數值模擬方法(輸運概率密度函數和均質混合氣模擬性)以及不同化學反應動力學機理的影響表明,由于Yao等[14]的機理在較低溫度條件下反應活性較強,有助于改善對于較低溫度,如850K條件下著火過程的預測.然而對于溫度更低的800K工況,以上方法均難以給出令人滿意的結果.Blomberg等[23]對比了雷諾平均和大渦模擬方法的模擬結果,發現采用LES可以更早捕捉到噴霧中多點著火現象,在模擬噴霧著火過程方面更具優勢.

圖14 對比了1.5ms時刻不同噴油壓力條件下的火焰結構.在相同噴射時間內,高噴油壓力下噴入燃燒室的燃油更多,霧化后的油束動能更高,速度更快,同時可燃混合氣更多.因此,噴油壓力為150MPa條件下的噴霧火焰面積更大,高溫火焰區域到噴孔的距離更遠.圖15為不同噴油壓力條件下混合氣溫度與混合物分數分布.從噴孔噴出的大量液滴以非常高的速度向下游運動,在火焰浮起長度位置處燃油蒸發與消耗達到平衡狀態.在噴油壓力為150MPa條件下,噴出的燃油更多,蒸發過程吸收大量的熱,導致溫度降低,不利于早期氧化過程產物累積,并且高速運動的油束推動反應區向下游運動,使得濃油區(混合物分數Z大于0.15)參與高溫反應的量明顯減少,油、氣混合過程得到改善,有助于減少碳煙排放.

圖14 1.5ms時刻不同噴油壓力下火焰結構Fig.14 Flame structures with different injection pressures at 1.5ms

圖15 標記了反應進度標量C的T-Z散點Fig.15 Scatters in the T-Z space colored by progress variable C

3 結論

(1) 基于大渦模擬方法探究了不同氧體積分數、環境密度和噴油壓力條件下正十二烷噴霧著火和火焰發展過程,并闡明了火焰結構以及主要組分分布受不同工況的影響規律.

(2) 證實具有高階精度的LES-LEM模型在寬工況范圍內具有廣泛適用性;在高氧體積分數條件下,燃燒區域溫度升高,化學反應速率加快;同時高氧體積分數下混合氣活性更高,導致火焰的浮起長度縮短,限制了燃油與空氣的混合時間,存在大量的濃混合氣區域;最終碳煙的排放是其氧化過程和生成過程競爭的結果,高氧體積分數條件下碳煙氧化過程發揮的作用更大,整體碳煙排放水平下降.

(3) 降低環境密度會導致著火推遲,使得油、氣混合更加充分,但是當環境密度非常低時,高溫著火被推遲,造成燃燒不完全;另外,提高噴油壓力可以改善油、氣混合過程,高速運動的油束推動高溫反應區向下游移動,導致火焰浮起長度變長,但噴油壓力對正十二烷噴霧著火過程影響較小.

主站蜘蛛池模板: 72种姿势欧美久久久大黄蕉| 91亚洲精选| 浮力影院国产第一页| 真实国产乱子伦视频| 伊人成色综合网| 亚洲综合极品香蕉久久网| 青青青视频免费一区二区| 亚洲欧美在线精品一区二区| 久久伊人操| 国产精品片在线观看手机版 | 国产乱子伦手机在线| 青青青国产精品国产精品美女| 欧美区一区二区三| 欧美亚洲国产视频| a免费毛片在线播放| 国产精品浪潮Av| 成人免费午夜视频| 国产网站黄| 又黄又爽视频好爽视频| 老司国产精品视频91| 99久久精品无码专区免费| 青草国产在线视频| 亚洲第一视频网站| 欧美国产在线看| 国产麻豆另类AV| 国产人成网线在线播放va| 国产午夜在线观看视频| 91精品啪在线观看国产| 国产成人精品2021欧美日韩| 福利一区在线| 一区二区三区在线不卡免费| 亚洲婷婷丁香| 99视频在线免费| 91色在线观看| 日韩精品成人网页视频在线| 极品国产一区二区三区| AV不卡无码免费一区二区三区| 中文字幕av无码不卡免费| 一区二区理伦视频| 在线欧美一区| 免费一极毛片| 亚洲人成人无码www| 黄色三级网站免费| 精品国产成人三级在线观看| 亚洲免费福利视频| 亚洲中文字幕在线精品一区| 国产精品视频观看裸模| 国产chinese男男gay视频网| 亚洲第一香蕉视频| 国产午夜一级毛片| 免费中文字幕在在线不卡| 久久国产精品波多野结衣| 亚洲福利网址| 色综合激情网| 中文字幕乱码二三区免费| 欧美另类第一页| 欧美日韩中文国产va另类| 亚洲三级色| 亚洲精品制服丝袜二区| 中国成人在线视频| 国产午夜福利亚洲第一| 国产91成人| 日韩欧美中文在线| 久久久国产精品免费视频| 国产精品深爱在线| 波多野结衣一区二区三区四区| 欧美一区二区三区国产精品| 国产成人精品亚洲77美色| 亚洲精品免费网站| 色综合国产| 午夜a视频| 一级毛片在线播放免费观看| 丝袜国产一区| 亚洲精品免费网站| 国产精品女人呻吟在线观看| 国产福利观看| 亚洲成人在线网| 亚洲精品无码在线播放网站| 国产精品高清国产三级囯产AV| 精品福利国产| 热久久综合这里只有精品电影| 免费aa毛片|