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

二沖程點燃式航空煤油發動機爆震預測模型研究

2023-06-08 09:22:12楊浩鵬楊海青魏民祥趙卓文
關鍵詞:發動機模型

楊浩鵬,楊海青,魏民祥,吳 昊,趙卓文,吳 昭

(1.南京航空航天大學 能源與動力學院, 南京 210016;2.南京長空科技有限公司, 南京 211800)

0 引言

航空活塞發動機以汽油作為燃料已經有很長一段歷史,從研究開發到生產制造也發展得非常成熟。但由于汽油飽和蒸氣壓高、閃點低、揮發性強的特點,使其在儲存、運輸和使用過程中存在相當大的安全隱患[1]。而航空煤油飽和蒸氣壓低、不易揮發且閃點較高(45~51 ℃)[2],滿足軍用燃料一體化的要求,因此航空活塞發動機燃用航空煤油已經成為國防所需。點燃式航空煤油發動機可以在現有航空汽油機的基礎上進行改造,研發成本低,并且該類型的發動機能夠繼承航空汽油機高功重比的優點。

發動機爆震是汽油機改燒煤油過程中不可忽視的技術障礙。相比于汽油,航空煤油的辛烷值(23~51)、自燃溫度(240~260 ℃)更低[3],并且航空煤油的火焰傳播速度較慢。這使得同樣為預混燃燒方式的點燃式航空煤油發動機更容易出現爆震現象。因此,研究以煤油為燃料的航空活塞發動機的爆震問題具有工程實際意義。目前存在多種理論解釋產生爆震的原因,但學者們廣泛接受的是自燃理論,即在火花塞點火之后,末端混合氣在火焰前鋒面到達之前自燃的現象。

對于發動機爆震的研究普遍采用試驗法和數值模擬的方法。胡春明等[4]對某雙火花塞點火的煤油發動機進行試驗,發現通過增大異步點火相位差可以降低爆震時的瞬間放熱率,能夠抑制爆震燃燒,但發動機的動力性能有所下降。Hess等[5]認為末端混合氣的自燃與燃燒開始前的預反應有關。他們提出了一種新的標準,將自燃開始之前未燃燒混合物的預反應考慮在內。通過燃用不同辛烷值的汽油機試驗進行驗證,發現預測的爆震起始曲軸轉角精度更高。李志銳等[6]對一臺單缸四沖程汽油機進行爆震試驗,建立了整機的一維仿真模型。通過對比兩類爆震模型的仿真結果,認為現有的爆震模型只能在部分工況下獲得較為準確的結果。陳龍華等[7]針對一臺增壓汽油機建立了同時考慮過量空氣系數和廢氣再循環的爆震預測模型。通過試驗數據確定了爆震模型的系數,仿真結果表明該爆震模型的預測性能較好,對不同發動機有較好的適應性。

綜上,由于發動機爆震會帶來很大的不可控風險,爆震試驗只能在有限的工況下進行,近年來采用數值模擬的方法預測爆震現象成為學者們研究的熱點。爆震預測模型在汽油機一維仿真模擬過程中的應用已發展得相當成熟[5-9],但是在航空煤油活塞發動機上的應用還并不完善,針對煤油燃料模型參數的設置略顯粗糙。此外,在這類模型的應用過程中,燃料的燃燒特性(如層流火焰傳播速度等)會對爆震預測精度以及一維仿真模擬的準確度產生較大影響,但相關的研究情況卻少有報道。

本文首先介紹了爆震預測模型的建立過程,并分析了模型計算所需的數據,然后基于一款二沖程點燃式航空煤油發動機的爆震試驗及相關爆震評價指標,根據已公開的RP-3航空煤油燃燒特性實驗數據建立煤油的燃燒模型,以一維仿真軟件GT-Power為平臺進行模型搭建及驗證,利用Matlab軟件編寫程序評估2種爆震預測模型的預測性能。

1 爆震預測模型

末端混合氣自燃引起的壓力振蕩被作為爆震發生的標志,為了預測爆震發生的時刻,學者們建立了眾多數學模型[8-11]。其中基于化學反應動力學模型的爆震預測計算量大,即使是簡化的動力學模型也因其龐大的計算量很難應用于工程實際。普遍使用的方法是監測燃料的著火延遲時間,可燃混合氣的著火延遲時間τ采用Arrhenius形式的經驗公式表示為壓力和溫度的函數,如式(1)所示。

(1)

式中:τ為著火延遲時間;P和T分別為壓力和溫度;A、n和B為經驗常數。

為了確定爆震是否發生,Livengood等[12]在1957年首次提出了著火延遲積分判斷法。這也是最初的爆震預測模型,如式(2)所示,當積分值達到1時,爆震發生時刻為tknock。

(2)

式中:τ為Arrhenius形式的經驗公式;tIVC為進氣門關閉時刻(對應二沖程發動機為掃氣口關閉時刻);tknock為爆震發生時刻。

Douaud等[13]將此方法應用于汽油發動機爆震預測的研究,通過發動機爆震實驗數據,開發并驗證了著火延遲時間的標準表達式,如式(3)所示。

(3)

式中:FON為燃料辛烷值;P為氣缸內壓力;T為末端混合氣溫度。

為了拓寬著火延遲積分預測爆震的應用范圍,Swarts等[14]考慮了過量空氣系數對著火延遲時間的影響,對標準表達式進行了修正。Hoepke等[15]在標準表達式中引入了EGR率,重新改進了著火延遲時間的表達式,使其能應用于帶有廢氣再循環系統的發動機。

然而這些工作大多是針對以汽油為燃料的發動機開展的,為了研究航空煤油的著火延遲特性,近年來國內外學者[16-19]通過化學激波管實驗獲得了航空煤油不同溫度、壓力條件下的著火延遲時間。劉靖等[17]在實驗數據的基礎上,利用多元線性回歸擬合得到了RP-3航空煤油及其多組分替代燃料著火延遲時間關系式,發現兩者的著火延遲特性基本吻合,其中RP-3航空煤油的單一Arrhenius形式的著火延遲時間如式(4)所示。

(4)

式中:φ為過量空氣系數。

已有研究表明[20],單一Arrhenius形式的著火延遲時間不能反映真實燃料復雜的燃燒放熱,Douaud模型只是反映了燃料在中溫狀態下的著火延遲響應,并不能很好地預測低溫或高溫狀態的自燃反應。考慮到燃料在燃燒放熱過程中的復雜性[21],為了反映真實燃料具有的兩級放熱和負溫度系數現象,建立了3-Arrhenius形式的著火延遲時間表達式,以此來模擬末端混合氣自燃時的低溫、中溫和高溫狀態,如式(5)和式(6)所示。

(5)

(6)

式中:τ1+τ2表示兩段低溫狀態的著火延遲時間;τ3表示高溫狀態的著火延遲時間。式(6)中,i=1,2或3;Ai、ni和Bi為經驗常數,對于這些經驗常數,Swarts[14]提出并驗證了正庚烷、異辛烷和甲苯這3種燃料的Ai、ni和Bi,分別在表1中列出。

表1 不同燃料的著火延遲模型參數

由于航空煤油化學組分復雜,并且存在負溫度系數現象,本文基于式(7)提出了一種由兩組分替代燃料計算的航空煤油混合著火延遲時間模型。

(7)

式中:j為燃料組分的編號;vj為燃料組分的體積分數;ω為經驗混合指數,取1.14[20]。

根據文獻[22]的研究,由89%正庚烷和11%甲苯組成的兩組分替代燃料能夠近似模擬航空煤油的化學動力學反應機理,并且通過射流攪拌器實驗驗證了該機理的可靠性。綜上所述,根據式(5)—(7)可計算航空煤油的混合著火延遲時間。圖1為航空煤油混合著火延遲時間表達式的計算流程,最終得到了航空煤油3-Arrhenius形式的混合著火延遲時間τk。

圖1 航空煤油混合著火延遲時間計算流程

可以看出缸內壓力P和末端混合氣溫度T對于爆震預測必不可少,后文在實驗和仿真的基礎上分別評價了單一爆震預測模型和基于3-Arrhenius公式的混合爆震預測模型對煤油發動機爆震起始角的預測性能,從而得到點燃式航空煤油活塞發動機最佳的爆震預測模型。

2 試驗系統及爆震評價指標

2.1 發動機試驗系統搭建

試驗所用發動機為德國LIMBACH公司生產的兩缸、水平對置、二沖程、風冷發動機,其具體參數見表2。試驗設備主要包括CWAC型交流電力測功機、Kistler 6113A火花塞式壓力傳感器、DH5960動態信號采集儀和ALM-ADV型空燃比分析儀等。試驗所用燃料為RP-3航空煤油,噴射方式為進氣道電噴,試驗過程中噴油壓力維持在0.3 MPa,進氣溫度為22~26 ℃。缸內壓力采樣頻率為200 kHz,保證氣缸壓力數據采集分辨率為 0.1°CA,每個工況點至少采集300個循環。

表2 發動機技術參數

本次發動機試驗包含穩態工況和爆震工況,分別記錄了發動機穩定運行時和發生爆震時的功率、燃油消耗率和氣缸壓力等數據,其中穩態工況下的試驗數據將用于發動機一維性能仿真模型的標定。發動機穩態工況的轉速為4 800~6 600 r/min、節氣門開度為48%~100%,共計27個工況點。由于爆震對發動機的危害極大,且在此發動機的試驗過程中發現轉速為4 800 r/min、節氣門開度為50%時,監測的氣缸壓力信號會出現輕微的振蕩,據此判斷發動機在該工況下易發生爆震。因此,為減小爆震對發動機的危害,爆震工況僅選擇轉速為4 800 r/min、節氣門開度為50%,試驗過程中只掃描點火提前角,使其從31°CA BTDC以1°CA的幅度不斷增加至36°CA BTDC。采集了發動機每個循環的缸內壓力變化。圖2為發動機試驗臺架。

圖2 發動機試驗臺架

2.2 試驗數據處理

發動機缸內燃燒的不穩定性導致運行過程中每個循環的燃燒壓力、燃燒放熱率等指標存在差異。有研究表明二沖程火花點火發動機燃用航空煤油后其指示平均有效壓力的循環變動率比燃用汽油時更高[23]。因此,為了獲得模型標定所需的有效數據,提取了每個工況點下連續300個壓力循環,并求得平均缸內壓力。圖3為發動機轉速在6 200 r/min、92%節氣門開度下的平均氣缸壓力處理結果,通過對每度曲軸轉角對應的不同循環缸壓求取平均值,獲得了該工況下的平均缸內壓力曲線。

爆震測量的方法有多種,而直接可靠的信號源只有爆震發生時燃燒室內的壓力振蕩,爆震時振蕩的典型頻率為5~10 kHz[24]。為了盡可能過濾正常燃燒時的低頻缸壓振蕩以及高頻的燃燒噪聲,保留有效的爆震信號,本文使用帶通濾波處理氣缸壓力信號來獲取爆震評價指標,帶通濾波的頻率范圍為3~25 kHz,可以覆蓋到爆震時壓力振蕩的典型頻率。取濾波后的最大壓力振蕩幅值(maximum amplitude of pressure oscillation,MAPO)作為爆震強度的評價指標,以此來區分爆震循環和非爆震循環[6-7]。圖4為發生爆震時缸內壓力的變化及帶通濾波的處理結果。根據試驗數據,采用文獻[6]中所述方法,對臨界爆震循環進行統計分析,最終確定爆震循環MAPO的閾值為0.08 MPa。圖5為點火提前角從31°CA BTDC增大到36°CA BTDC時900個循環的MAPO統計結果,縱坐標為不同MAPO出現的頻率,篩選獲得了405個循環的爆震缸壓。

圖4 爆震工況下缸內壓力變化及壓力振蕩幅度

圖5 不同MAPO的統計結果

獲得準確爆震起始角(knock onset,KO)是驗證爆震預測模型的前提。根據爆震機理,發動機發生爆震燃燒時缸內壓力會出現突變,因此采用缸內壓力曲線上第一個瞬間凸起對應的曲軸轉角作為爆震起始角。

利用Matlab讀取發生爆震時的缸壓數據,通過ginput函數獲得了氣缸壓力突變時刻所對應的曲軸轉角。圖6為逐漸增大點火提前角時連續200個循環的爆震起始角的統計結果。隨著點火提前角的增大,缸內最大爆發壓力所對應的曲軸轉角提前,導致爆震起始時刻逐漸提早,由于發動機燃用煤油后循環變動率相比于燃用汽油時更高[23],使得測得的爆震起始時刻存在較大波動。

圖6 連續200個循環的爆震起始角

3 末端混合氣溫度獲取

決定爆震發生時刻的著火延遲期與末端混合氣的壓力和溫度有關,氣缸內壓力隨曲軸轉角的變化可以通過實驗獲得,但是末端混合氣的溫度很難通過實驗獲得,普遍采用的方法是基于末端混合氣絕熱壓縮求解得到[7]。定義進氣口關閉時刻的缸內溫度為初始溫度,在絕熱壓縮條件下,可由式(8)計算得出末端混合氣溫度。

(8)

式中:T2為末端混合氣溫度;T1為初始溫度;p1為進氣口關閉時刻的氣缸壓力;p2為隨曲軸轉角變化的氣缸壓力;γ為混合氣的比熱比;其中p1和p2均由實驗獲得。

由于初始溫度受到缸體換熱以及燃油混合的影響,導致其與實測的進氣溫度并不相等,因此,為了獲得較為準確的初始溫度,需要建立該發動機的一維性能仿真模型。根據實驗發動機的工作原理和結構尺寸,利用GT-Power對發動機的進排氣系統、曲柄連桿等結構進行了建模,發動機的一維仿真模型如圖7所示。

圖7 發動機一維仿真模型

一維仿真模型中,傳熱和燃燒模型會顯著影響初始溫度的計算,本文選用普遍用于發動機數值模擬中的WoschniGT傳熱模型[25],該模型具有良好的計算精度。燃燒模型選用火花點火湍流燃燒模型(SI turbulent flame combustion model),然而軟件中并沒有關于航空煤油的模型參數。因此本文基于文獻[16]中RP-3航空煤油燃燒試驗數據計算獲得了煤油的層流火焰衰減速度、溫度指數以及壓力指數等信息,應用于GT-Power燃燒模型參數的設置,通過實驗數據對整個模型進行了校核,最終獲得較為準確的初始溫度,并利用式(8)計算得出末端混合氣溫度。

3.1 燃燒模型參數的獲取及校核

層流火焰傳播模型的參數需要通過計算獲得,不同燃料的層流火焰傳播速度也不盡相同,因此模型參數與燃料的種類有很大的關系。球形擴散火焰是火焰傳播模型的基本假設,為了得到層流火焰傳播速度隨壓力、溫度變化的明確表達式,常用的方法是將實測的層流燃燒速度擬合為壓力和溫度的指數關系[26-27],如式(9)所示。

(9)

式中:Sl為層流火焰傳播速度;S0為參考壓力溫度下的層流火焰傳播速度;T為混合氣溫度;P為壓力;T0為參考溫度;P0為參考壓力;α為溫度指數;β為壓力指數。

煤油的燃燒實驗表明式(9)中的層流火焰傳播速度由過量空氣系數φ、混合氣溫度T以及壓力P決定,方程系數α和β根據試驗測得的RP-3航空煤油燃燒數據計算得到,圖8為文獻[16]中不同溫度壓力下層流火焰傳播速度與過量空氣系數之間的關系。

圖8 不同溫度壓力下層流火焰傳播速度隨過量空氣系數的變化

據此可計算出溫度指數α和壓力指數β隨過量空氣系數φ的變化關系,結果如表3所示。

表3 溫度指數和壓力指數

GT-Power的燃燒模型對式(9)進行了部分修正,考慮了實際工程中可能會應用的廢氣再循環技術,同時將S0替換為與過量空氣系數相關的函數,其層流火焰傳播速度SL表示為

(10)

式中:Bm為參考溫度和參考壓力下的最大層流火焰傳播速度;φm為最大層流火焰傳播速度對應的當量比;Bφ為層流火焰衰減速度;φ為過量空氣系數;fEGR為廢氣再循環率的函數。由于暫不考慮廢氣再循環技術,此時fEGR為1。根據系數定義和試驗數據可知φm為1.2,Bm為0.761 m/s,為方便計算,Bφ取不同壓力溫度下層流火焰衰減速度的平均值,即-0.701 m/s。

實際的燃燒過程非常復雜,有規律可循的層流火焰傳播只是其一部分,在火焰傳播中湍流的影響也不可忽視[28]。湍流與詳細化學反應機理之間存在強烈的非線性相互作用,且湍流運動具有無序性和隨機性,目前尚未有普遍適用的湍流火焰傳播模型。

為方便建立不同發動機的湍流火焰傳播模型,需要在GT-Power軟件內設置3個湍流火焰傳播的模型參數,分別為火核生長速度因子FKG(flame kernel growth multiplier),湍流火焰速度因子TFS(turbulent flame speed multiplier),泰勒微尺度縮放因子TLS(Taylor length scale multiplier)。3個參數的默認值均為1,對于不同發動機存在一組最合適的參數值,因此需要根據發動機穩態工況下的試驗數據優化這3個參數。本文利用GT-Power軟件內置的參數優化板塊DOE(design of experiments)對這3個模型參數進行優化,根據GT-Power使用手冊,選擇參數優化范圍為0.5~3.0,選擇拉丁超立方(Latin hypercube sampling)的采樣方法,分別對這3個參數在優化范圍內取值,得到了1 000種模型參數的組合方式。隨后對這1 000個樣本進行試算,以發動機穩態工況下的功率、燃油消耗率和缸內壓力作為參數優化依據,篩選出仿真值與實驗值誤差最小的樣本。圖9為轉速在4 800 r/min、節氣門開度為82%時針對功率的參數優化結果。對篩選出的213個樣本依次進行針對燃油消耗率和缸內壓力的參數優化,得到一組此工況下發動機湍流火焰傳播的模型參數。

使用上述方法對該發動機的其他工況點進行模型參數優化,最終獲得發動機全工況下完整的一維仿真模型。圖10為轉速在4 800、5 200、5 600和6 200 r/min時不同節氣門開度的實驗功率和仿真功率對比,兩者誤差在5%以內。圖11為轉速在4 800 r/min、不同節氣門開度下實驗的燃油消耗率與仿真的燃油消耗率對比,兩者誤差在8%之內。圖12為轉速在4 800 r/min、節氣門開度為50%時,300個循環的實驗平均缸壓與仿真缸壓的對比,兩者誤差在6%之內。

圖10 不同節氣門開度下的實驗與仿真功率對比

圖11 實驗與仿真燃油消耗率對比

圖12 實驗與仿真缸壓對比

3.2 計算末端混合氣溫度

通過對上節燃燒模型參數的計算和校核,得到了該發動機完整的一維仿真模型,可以較為準確地計算進氣結束時刻的缸內溫度。為了計算末端混合氣溫度,還需確定混合氣的比熱比。根據文獻[7]可知,混合氣的溫度、組分以及各組分所占的比例均會影響混合氣的比熱比,隨著過量空氣系數的減小,混合氣的比熱比會隨之減小,而隨著溫度的升高混合氣的比熱比也有減小趨勢,當發動機活塞上下運動時,氣缸內的溫度和組分會有較大的改變。因此采用定比熱比計算末端混合氣溫度會產生很大的計算誤差,而末端混合氣溫度會直接影響到著火延遲積分的計算結果,從而影響爆震預測的精度。

本文采用校核后的一維仿真模型計算出比熱比隨曲軸轉角的變化關系,通過式(8)計算每個循環的末端混合氣溫度。圖13為4 800 r/min、節氣門開度為50%時,某一循環采用變比熱比計算末端混合氣溫度的計算結果。

圖13 某一循環末端混合氣的計算結果

4 爆震預測模型的性能評估

為了定量分析單一Arrhenius形式的爆震預測模型與本文提出的基于3-Arrhenius公式的混合爆震預測模型對爆震起始時刻預測的準確性,使用Matlab輸入溫度、壓力和2種模型的表達式等數據,找出了著火延遲積分恰好等于1的時刻,并將其換算為對應的曲軸轉角記為預測的爆震起始角(predicted knock onset,PKO),然后對比PKO與試驗的爆震起始角(experimental knock onset,EKO)之間的誤差即可評估爆震預測模型的性能。

圖14為單一爆震預測模型和混合爆震預測模型對同一輕微爆震循環的預測結果。當著火延遲積分等于1時,對應的曲軸轉角即為預測的爆震起始角。由圖14(a)可以看出,單一爆震預測模型對著火初期的低壓區很敏感,著火延遲積分在發生爆震之前就快速增加,不能識別爆震發生時缸內壓力的突然升高,導致預測的爆震起始角略小于實際的爆震起始角。由14(b)可以看出,混合爆震預測模型能夠識別著火初期的低壓區。當缸內壓力因著火而快速升高時,著火延遲積分也隨之快速增大,從而能夠對爆震發生時缸內壓力的突然升高做出反應,預測的爆震發生時刻更加接近實際的爆震發生時刻。

圖14 不同爆震預測模型對同一輕微爆震循環的預測結果

圖15為單一爆震預測模型和混合爆震預測模型對同一嚴重爆震循環的預測結果。

嚴重爆震時缸內壓力會突然升高同時伴隨著高頻的振蕩,而單一爆震預測模型不能很好地對缸內壓力的突然升高做出反應,導致預測的爆震起始角滯后于實際的爆震起始角。混合爆震預測模型在這方面表現出了更好的預測性能,當壓力因爆震而突然升高時,著火延遲積分也隨之快速增大,最終預測的爆震起始角更加接近實際的爆震起始角。

基于3-Arrhenius公式的混合爆震預測模型對于輕微爆震和嚴重爆震都具有良好的預測性能,混合爆震預測模型能夠區分正常燃燒與爆震燃燒的氣缸壓力變化。圖16為單一爆震預測模型和混合爆震預測模型對所有爆震循環的預測結果,2個模型預測的爆震起始角與實際的爆震起始角之間的標準誤差分別為3.23°CA和1.19°CA。

圖16 不同爆震預測模型對所有爆震循環的預測結果

5 結論

1) 根據正庚烷和甲苯的著火延遲模型參數,提出了基于3-Arrhenius公式的航空煤油兩組分混合著火延遲時間計算式。

2) 采用已公開的航空煤油燃燒數據,在GT-Power軟件平臺上搭建了煤油的燃燒模型,并根據發動機實驗數據對該燃燒模型進行了驗證。

3) 混合爆震預測模型的PKO與EKO之間的標準誤差為1.19°CA,而單一爆震預測模型的PKO與EKO之間的標準誤差為3.23°CA,混合爆震預測模型更加符合真實燃料的燃燒過程,表現出了更好的預測性能。

猜你喜歡
發動機模型
一半模型
元征X-431實測:奔馳發動機編程
2015款寶馬525Li行駛中發動機熄火
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
新一代MTU2000發動機系列
發動機的怠速停止技術i-stop
新型1.5L-Eco-Boost發動機
主站蜘蛛池模板: 麻豆国产在线不卡一区二区| 免费福利视频网站| 91色在线观看| 伊人久久青草青青综合| 色噜噜综合网| 人人艹人人爽| 成人年鲁鲁在线观看视频| 久久a级片| 国产91麻豆视频| 丰满少妇αⅴ无码区| 日韩福利视频导航| 99热在线只有精品| 狠狠操夜夜爽| 国产白浆在线| 国产成人精品优优av| 亚洲色偷偷偷鲁综合| 国产主播福利在线观看| 午夜久久影院| 欧美成人午夜视频免看| 天天综合色网| 国产精品页| 国产精品人成在线播放| 久久人妻xunleige无码| 国产精品无码久久久久AV| 欧美成人亚洲综合精品欧美激情| 日本黄色a视频| 国产乱论视频| 欧美日韩另类在线| 国产波多野结衣中文在线播放 | 一级毛片基地| 少妇极品熟妇人妻专区视频| 国产喷水视频| 中文字幕有乳无码| 青青青国产视频| 国产一级在线观看www色| 国内精品伊人久久久久7777人| 中文一级毛片| 亚洲综合天堂网| 天堂成人在线视频| 性激烈欧美三级在线播放| 成人看片欧美一区二区| 国产精品视频白浆免费视频| 丁香婷婷激情综合激情| 波多野结衣久久高清免费| 亚洲综合色吧| 亚洲欧美一区在线| 国产精品无码AV片在线观看播放| 日韩 欧美 小说 综合网 另类 | 亚洲一区二区精品无码久久久| a级毛片在线免费| 欧美不卡在线视频| 久久国产精品波多野结衣| 精品亚洲国产成人AV| 久久国产乱子| 91无码人妻精品一区二区蜜桃| 天堂成人在线| 青青久视频| 波多野吉衣一区二区三区av| 特级欧美视频aaaaaa| 欧美成人亚洲综合精品欧美激情| 国模私拍一区二区| 国产亚洲精品无码专| 亚洲手机在线| 欧美有码在线| 国产小视频在线高清播放| 精品国产电影久久九九| 亚洲AV无码久久天堂| AV在线麻免费观看网站| 国产亚洲精品91| 亚洲国产午夜精华无码福利| 熟女视频91| 国产第三区| 中文字幕久久波多野结衣| 无码'专区第一页| 成人国产精品视频频| 亚洲欧美另类日本| 国产一区二区三区免费| 美女被狂躁www在线观看| 9久久伊人精品综合| a级毛片免费在线观看| 亚洲品质国产精品无码| 国产原创自拍不卡第一页|