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

提高天然氣燃燒機理對點火延遲時間預測的準確性研究

2025-05-23 00:00:00張軍毛功平齊虹
車用發動機 2025年2期

天然氣是一種重要的內燃機替代清潔低碳燃料。預計到2040年,天然氣的同等能源價格仍將低于石油燃料[1]。從環境和成本角度來綜合考慮,天然氣將是未來一段時間的主流燃料。因此,實現天然氣燃料的高效清潔燃燒,仍然是當前的研究熱點。

將天然氣的燃燒反應機理與CFD軟件進行耦合計算,是深入分析天然氣發動機燃燒過程的最有效手段之一[2]。近幾十年來,國內外學者提出了多種版本的天然氣機理。M.FRENKLACH等[3]提出的GRI-Mech3.0機理包含53種組分和325個基元反應,是一種使用非常廣泛的天然氣機理。文獻[4]提出了包含58種組分和270個基元反應的UniversityofCalifornia SanDiego 2016(UCSD2016)機理,主要適用于層流火焰速度預測、高溫點火和爆炸等工況。為了擴展天然氣機理的適用工況范圍,北京航空航天大學的C.W.ZHOU等[5在AramcoMech2.0機理的基礎上添加了低溫CO-C4含氧碳氫化合物的子機理,開發了包含581種物質和3037個基元反應的AramcoMech3.0機理。該機理能較準確地預測CO-C4燃料在低、中溫條件下的燃燒特性。沈陽航空航天大學的宋晨星基于Aramco2.0反應機理,通過路徑敏感性分析、生成速率分析與反應路徑分析,形成了包含21種組分、150個反應的天然氣 的初始簡化反應機理,基于解耦法,耦合初始簡化反應機理中的C1-C3反應機理、 詳細反應機理和 簡化反應機理,構建了包含40種組分、189個反應的天然氣的簡化反應機理。南京師范大學的徐傳義等8采用耦合誤差傳播和敏感性分析的直接關系圖(DRGEP-SA)法,對GRI-Mech3.0機理進行了簡化,構建出了含有23種組分和110個基元反應的簡化機理,采用充分攪拌反應器模型和一維層流預混反應器模型,分別對甲烷燃燒的溫度和組分摩爾分數分布進行了對比分析。研究結果表明,利用簡化機理分析得到的預測數據與詳細機理吻合良好。然而,由于化學反應路徑的復雜性、化學反應速率的不確定性等原因,在低溫 (800~1200K) 燃燒條件下,天然氣機理的預測性能有著不同程度的降低。

點火延遲時間是燃料最重要燃燒特征參數之一,對發動機的燃燒和排放特性有重要的影響。能否準確地預測點火延遲時間,是衡量燃料機理預測性能高低的重要標準。為此,國內外學者以點火延遲時間的準確預測為主要目標,開展了提高天然氣機理預測性能的研究。在GRI-Mech1.2機理的基礎上,E.L.PETERSEN等9增加了6種物質和13個基元反應,新機理包含38種物種和190個基元反應,預測的點火延遲時間與試驗值的誤差在15% 以內。S.HEYNE等[10]以GRI-Mech3.0詳細機理為基礎,采用多目標優化的進化算法,修正了敏感性系數較高的6個關鍵基元反應的速率常數,提高了對點火延遲時間的預測精度。J.HUANG等[11]發現GRI-Mech3.O機理預測的點火延遲時間與試驗值的誤差隨溫度的降低而增大,因此,修改了 的化學反應速率常數,并添加了 兩個基元反應,新機理在溫度為 時預測點火延遲時間的誤差僅為 5% 。長沙理工大學的黃章俊等[12]采用帶誤差傳播的直接關系圖法、全物種敏感性分析和人工神經網絡(ANN)聯合方法,以點火延遲時間和CO摩爾分數為優化目標,對USCMech2.0機理進行了簡化和優化,優化機理預測點火延遲時間、OH摩爾分數峰值和CO摩爾分數峰值的相對誤差均在 10% 以下。

可見,國內外學者通過添加基元反應、完善反應過程、修正化學反應速率常數等方法,提高了天然氣機理的預測性能。

在阿倫尼烏斯公式中,活化能 )為指數項,對機理的預測性能有重要影響。現有機理中基元反應的活化能數據主要來源于經驗值、試驗和理論計算,溫度不同時,活化能的取值差異較大,具有較大的不確定性[13]。此外,現有的天然氣機理大多無法體現低溫條件下活化能的變化[11]。已有文獻中也鮮有通過活化能的修正來提高機理預測性能的相關報道。隨著量子化學計算、超級計算機等技術的飛速發展,活化能的計算也變得日益成熟。因此,本研究擬采用五種量子化學計算方法(B3LYP、PM7、M06-2X、CBS-QB3、W1RO),在溫度為 800~ 1200K的條件下,計算影響點火延遲時間的關鍵基元反應的活化能,對機理中的活化能數據進行修正,比較分析修正前后機理預測的點火延遲時間,并與試驗值進行對比,以提高天然氣燃燒機理在低溫條件下預測點火延遲時間的準確性。

1天然氣詳細燃燒機理的對比

國內外學者構建了多種天然氣的詳細燃燒機理,其中,使用較廣泛的主要有三種:1)包含58種物質和270個基元反應的UniversityofCaliforniaSanDiego2016機理,記為UCSD;2)包含53種物種和325個基元反應的GRI3.0機理,記為GRI;3)包含581種物質和3037個基元反應的AramcoMech3.0機理,記為Aramco。為了篩選出合適的研究對象,在溫度 T 為 800~1400K 、壓力 p為1MPa 、當量比 為1.0和2.0的條件下,將上述3種詳細機理與CHEMKIN軟件中的閉式均相反應器模型進行耦合計算,將點火延遲時間的計算值與文獻中的試驗值[14]進行對比。

天然氣的成分比較復雜,一般由甲烷( )、乙烷 和丙烷 等短鏈的碳氫化合物組成,其含量隨產地的不同而各異[15]。為方便研究,本研究在計算過程中假設天然氣由體積分數為90% 的 F 9% 的 和 1% 的 組成。

圖1示出三種天然氣詳細機理計算的點火延遲時間與試驗值的對比。由圖1可以看出,三種機理預測的點火延遲時間均隨著溫度的升高而縮短,主要是因為溫度較高時,燃料更容易著火。在兩種不同當量比工況下,溫度高于 時,Aramco的點火延遲時間較其他兩種機理略長;當溫度為 900~ 時,Aramco與GRI的點火延遲時間很接近;當溫度低于 900K 時,三種機理的點火延遲時間均有最大誤差。本研究重點關注 的溫度范圍,因此,對該溫度區間的點火延遲時間進行定量分析。采用式(1)計算點火延遲時間與試驗值的相對誤差。

式中:e為點火延遲時間計算值與試驗值的相對誤差; 為詳細機理對點火延遲時間的計算值; 為點火延遲時間的試驗值。

圖1三種天然氣詳細機理預測的點火延遲時間與試驗值的對比

圖2示出點火延遲時間與試驗值相對誤差的對比。由圖2可以看出,兩種工況下,與試驗值相比較,在高溫 左右,GRI對點火延遲時間的預測誤差較小,分別為 4.6% 和 23.9% ;在其余溫度條件下,Aramco對點火延遲時間的預測誤差較小,一般在 15%~40% ,誤差較大的情況出現在低溫970K左右,最大誤差接近 40% 。由此可見,在大多數工況條件下,Aramco機理預測點火延遲時間的誤差均小于GRI和UCSD機理。因此,本研究選擇Aramco為詳細機理,重點針對低溫時誤差較大的情況進行改進。

圖2三種天然氣詳細機理預測的點火延遲時間與試驗值的相對誤差的對比

2天然氣詳細燃燒機理的改進

為了提高Aramco機理預測點火延遲時間的準確性,采用量子化學計算方法,計算機理中關鍵基元反應在 溫度下的活化能,對原機理中的活化能數據進行修正,從而對機理進行改進,主要分為以下5個步驟:

1)計算不同溫度條件下各基元反應對點火延遲時間的敏感性系數,找出敏感性系數最大的基元反應,將其作為關鍵基元反應;2)利用Gaussian-09W軟件[16],計算不同溫度下關鍵基元反應的反應物和產物的焓值;3)基于反應物和產物的焓值,計算反應物和產物的焓變,將其近似為關鍵基元反應的活化能;4)利用不同溫度下活化能的平均值,替代原機理中關鍵基元反應活化能數據,得到優化機理;5)將優化機理與閉式均相反應器模型進行耦合,計算點火延遲時間,與詳細機理的計算值和試驗值進行對比分析,檢驗優化機理的預測性能。

2.1 敏感性分析

圖3示出在不同工況下主要基元反應對點火延遲時間的敏感性系數。可以看出,在溫度為 950K 時,兩種壓力和當量比條件下,與其他基元反應相比,基元反應 (M代表反應的中間體)的敏感性系數最大,具有最高的反應活性。這主要是因為 的溫度下會快速分解,產生兩個高反應性羥基(—OH)自由基[17]。根據著火理論,—OH是一種重要的活化分子,是導致鏈鎖反應的重要中間物質,一OH濃度的增加為著火創造了必要的條件[18]。該反應被認為是在 850~1200K 時控制碳氫化合物點火的主要鏈支化反應,是火花點火發動機爆震、柴油機著火和均質混合壓燃發動機(HCCI)運行的核心動力學特征[19]

圖3天然氣點火延遲時間的敏感性分析

2.2 活化能的計算

屬于無過渡態的吸熱反應,因此,采用反應焓近似活化能的方法進行計算[20],如式(2)所示。

式中: 為活化能; 為焓差; 分別為產物(一OH)和反應物 )的焓值。

為提高計算結果的準確性,本研究采用 ,CBS-QB3[23],B3LYP/6-31G 種不同的理論方法,對 和 -OH 的幾何結構進行優化和頻率計算。這些方法覆蓋了從半經驗方法到高精度復合方法的多種理論水平,以確保對比全面性。對于每種方法,首先對反應物( )和產物(一OH)進行幾何優化,確保振動頻率均為實數(表明處于穩定態);然后進行頻率計算,獲取體系的焓值校正項。多種理論方法的選擇旨在平衡計算精度和資源需求。其中,PM7由于計算速度快可用于初步分析,但存在誤差較大問題;M06-2X在描述分子間相互作用方面具有優勢;B3LYP和CBS-QB3則在計算小分子體系的熱力學性質時表現優異,尤其是CBS-QB3,被廣泛應用于高精度的焓值和活化能預測。利用式(2),分別計算了不同溫度下基元反應 的活化能,計算結果如表1所示。

表1利用不同方法計算的基元反應 在不同溫度下的活化能

從表1可以看出,PM7和M06-2X/6-31G(d)的計算結果偏差較大,這是由于PM7是一種半經驗方法,摻雜了較多的經驗公式[26],而M06-2X/6-31G(d)雖然在計算諧波振動頻率方面的精度較高,但在計算焓值方面存在不足[27]。因此,本研究重點考慮B3LYP/6-31G(d)、CBS-QB3和W1RO方法的計算結果,并對不同溫度下求得的活化能取平均值,以拓寬機理的適用范圍。

2.3優化機理的預測性能分析

為了驗證優化機理的準確性,采用CHEMKIN中的閉式均相反應器(closedhomogeneousbatchreactor,CHBR)模型與詳細機理、優化機理進行耦合計算,將不同工況下計算的點火延遲時間與試驗數據進行對比分析。

圖4示出詳細機理和優化機理點火延遲時間計算值與試驗數據的對比。由圖4可以看出,4種工況下,當溫度為 時,與三種新機理相比,詳細機理的預測值最接近試驗值;當溫度為900~1025K 時,三種優化機理的預測值較詳細機理更接近試驗值,其中,基于B3LYP方法得到的優化機理的預測值與試驗值最為接近。其原因可能是:利用B3LYP方法求解的活化能比其他兩種方法都要小(如表1所示),而活化能是決定化學反應速率的內因,活化能越小,則化學反應速率越大,使點火延遲時間縮短[28],從而能減少試驗過程的散熱損失。

圖4不同機理預測的點火延遲時間與試驗值的對比

天然氣的著火延遲時間在 10~100ms ,是非常短暫的過程,本研究的仿真計算采用了閉式均相反應器模型,將反應過程視為絕熱反應,未考慮散熱損失[10]。因此,活化能較小時,數值計算的結果也較準確,表現為較小的活化能對應的點火延遲時間計算值和試驗值更為接近。當溫度為 .當量比為2.0時的誤差比當量比為1.0時的誤差更小。這是由于在富燃情況下,高濃度的燃料提供了更多的燃料分子和一OH自由基,這更加易于著火,點火延遲時間的計算值也較小,較低的傳熱損失也間接降低了計算值與試驗值的誤差[16]。

值得注意的是,B3LYP是一種基于密度泛函的方法,近期提出的相關理論還包括了對色散效應的處理和對長程相互作用的處理,在一定程度上提高了計算精度[29],而W1RO和CBS-QB3都是后自洽場方法,理論上W1RO和CBS-QB3計算活化能的精度應該更高[30]。然而,在本研究的實際應用過程中,利用B3LYP方法計算得到的活化能對機理進行改進后,機理的預測性能更優。這說明B3LYP方法可能更適用于對天然氣低溫燃燒機理活化能數據的修正。在后續研究中,還可以將這種方法用于對其他燃料燃燒機理的修正,在多種工況下檢驗該方法的合理性。

為了定量評價誤差,根據式(1),計算了3種優化機理預測結果的相對誤差。圖5示出3種優化機理預測的點火延遲時間與試驗值相對誤差的對比。由圖5可以看出,四種工況下,與其他兩種方法相比,B3LYP的相對誤差分別為 4%,1.8%,7.3% 18.4% ,均小于W1RO和CBS-QB3的誤差,也小于圖2中詳細機理的相對誤差 33.6% 30% . 42.1% ,62.5% 。可見,采用B3LYP方法對機理進行改進,能有效提高機理在低溫條件下預測點火延遲時間的準確性。

圖5三種新機理預測的點火延遲時間的對比

3結論

a)在壓力為 1MPa ,溫度為970K左右,當量比為1.0和2.0工況條件下,采用Aramco機理預測的點火延遲時間與試驗值的平均誤差為 33.6% ,小于UCSD機理和GRI機理計算的平均誤差,能較準確地預測低溫條件下的點火延遲時間;

b)與W1RO方法和CBS-QB3方法的計算結果相對比,利用B3LYP/6-31G(d)方法計算的活化能對 進行改進后,4種工況下優化機理預測的點火延遲時間與試驗值的誤差分別為 4%,1.8%,7.3%,18.4% ,且小于原機理的誤差 33.6%,30%,42.1%,62.5%

參考文獻:

[1] FETISOVV,GONOPOLSKYAM,DAVARDOOST H,etal.Regulationand impact of VOC and emissions on low-carbon energy systems resilient to climate change:A case study on an environmental issue in the oil and gas industry[J].Energy science amp; engineering, 2023,11(4):1516-1535.

[2] LAMIONIR,BRONZONI C,FOLLI M,etal.Impact of -enriched natural gas on pollutant emissions from domestic condensing boilers:Numerical simulations of thecombustion chamber[J].International journal of hydrogen energy,2023,48(51):19686-19699.

[3]FRENKLACH M,WANG H,GOLDENBERG M,et al.GRI-Mech:an optimized detailed chemical reaction mechanism for methane combustion[C].Chicago:Gas Research Institute Chicago,1995.

[4]加州大學圣地亞哥分校燃燒研究小組.圣地亞哥機理 [EB/OL].(2016-12-14).https://web.eng.ucsd.edu/ mae/groups/combustion/mechanism.html.

[5]ZHOU C W,LI Y,BURKE U,et al.An experimental and chemical kinetic modeling study of 1,3-butadiene combustion:Ignition delay time and laminar flame speed measurements[J].Combustion and flame,2018, 197:423-438.

[6]METCALFEWK,BURKESM,AHMEDSS,et al.A hierarchical and comparative kinetic modeling study of C1-C2 hydrocarbon and oxygenated fuels[J].International journal of chemical kinetics,20l3,45(10):638- 675.

[7]宋晨星.天然氣燃燒反應簡化機理研究[D].沈陽:沈陽 航空航天大學,2023.

[8]徐傳義,邵亞麗,盧平,等.基于DRGEPSA方法的甲烷 燃燒簡化機理的構建及驗證[J].上海電力大學學報, 2022,38(5):415-420.

[9]PETERSENEL,DAVIDSONDF,HANSONRK. Kinetics modeling of shock-induced ignition in low-dilution mixtures at high pressures and intermediate temperatures[J].Combustion and flame,1999, 117(1/2):272-290.

[10]HEYNE S,ROUBAUD A,RIBAUCOUR M,et al. Development of a natural gas reaction mechanism for engine simulations based on rapid compression machine experiments using a multi-objective optimisation strategy[J].Fuel,2008,87(13/14):3046-3054.

[11]HUANG J,HILL PG,BUSHE WK,et al.Shocktube study of methane ignition under engine-relevant conditions:experiments and modeling[J].Combustion and flame,2004,136(1/2):25-42.

[12]黃章俊,徐通,何洪浩,等.基于人工神經網絡的甲烷 富氧燃燒機理優化[J].動力工程學報,2024,44(4): 520-527.

[13]美國國家標準與技術研究院.化學動力學數據庫[EB/ OL].(2024-10-1).https://kinetics.nist.gov/kinetics.

[14]DROST S,AZNAR MS,SCHIEBL R,et al.Reduced reaction mechanism for natural gas combustion in novel power cycles[J].Combustion and flame,2021, 223:486-494.

[15]SPADACCINI L J,COLKET ⅢI M B.Ignition delay characteristics of methane fuels[J].Progress in energyand combustion science,1994,20(5):431-460.

[16]TOMBERG A.Gaussian 09W tutorial[M]//An intro duction to computational chemistry using and avogadro software.[S.l.]:[s.n.],2013:1-36.

[17]HONG Z,FAROOQ A,BARBOUR EA,et al.Hydrogen peroxide decomposition rate:a shock tube study using tunable laser absorption of near .The journal of physical chemistry A, 2009,113(46);12919-12925.

[18]SHANG Y,SHIJ,NING H,et al.Significance of reaction in two-stage ignition of nitromethane[J].Fuel,2019,256:115956.

[19]WESTBROOKC K.Chemical kinetics of hydrocarbon ignition in practical combustion systems[J].Proceedings of the combustion institute,20oo,28(2):1563- 1577.

[20]MORA-DIEZ N,ALVAREZ-IDABOYJR,BOYD R J.A quantum chemical and TST study of the OH hydrogen-abstraction reaction from substituted alde hydes:FCHO and ClCHO[J].The journal of physical chemistry A,2001,105(39):9034-9039.

[21] STEWAR TJJP.Optimization of parameters for semiempirical methods V:Modification of NDDO approximations and application to 70 elements[J].Journal of molecular modeling,2007,13:1173-1213.

[22]XIA Y,ZHANG H,ZHU X,et al.Fluorescent probes for detecting glutathione:bio-imaging and two reaction mechanisms[J].Dyes and pigments,2o19,163: 441-446.

[23]TAHAN A,SHIROUDI A.Oxidation reaction mechanism and kinetics between OH radicals and alkylsubstituted aliphatic thiols:OH-addition pathways [J].Progress in reaction kinetics and mechanism, 2019,44(2):157-174.

[24]BECKE A D.Density-functional thermochemistry III: the role of exact exchange[J].The journal of chemical physics,1993,98(7):5648-5652.

[25]ZHAO Y,TRUHLARDG.The MO6 suite of density functionals for main group thermochemistry,thermochemical kinetics,noncovalent interactions,excited states,and transition elements:two new functionals and systematic testing of four Mo6-class functionals and 12 other functionals[J].Theoretical chemistry accounts,2008,120(1):215-241.

[26]MINENKOVY,SHARAPAD I,CAVALLO L.Application of semiempirical methods to transition metal complexes:Fast results but hard-to-predict accuracy[J].Journal of chemical theory and computation, 2018,14(7):3428-3439.

[27] LYNCHBJ,ZHAOY,TRUHLARDG.Effectiveness of diffuse basis functions for calculating relative energiesbydensity functional theory[J].The journal of physical chemistry A,2003,107(9):1384-1388.

[28] ZHUW T,QIUXP.Study ontherelationship betweenactivation energy,heat of reaction and enthalpy of activation in chemical kinetics[J].Journal of tsinghua university:science and technology,1999,39(6): 23-24.

[29] AUSTINA,PETERSSONGA,FRISCHMJ,etal. A density functional with spherical atom dispersion terms[J].Journal of chemical theory and computation,2012,8(12):4989-5007.

[30] SONNENBERGJL,WONGKF,VOTHGA,etal. Distributed Gaussian valence bond surface derived from ab initio calculations[J].Journal of chemical theoryand computation,2009,5(4):949-961.

Prediction Accuracy Improvement of Natural Gas Combustion Mechanisms forIgnition Delay Time

ZHANG Jun1,MAO Gongping2,QI Hong3 (1.School of Automotive Engineering,Nantong Institute of Technology,Nantong 226o02,China; 2.School of Automotive and Traffic Engineering,Jiangsu University,Zhenjiang 2l2ol3,China; 3.School of Mechanical Engineering,Nantong Institute of Technology,Nantong226o02,China)

Abstract:The methodofquantum chemistry wasused tocalculate theactivation energyof key elementary reactions in the naturalgasmechanismundertheconditionoflow temperaturecombustion,soastomodifytheactivationenergydatainthe originalmechanism,andthereforeimprovetheacuracyofnaturalgas mechanismtopredicttheignitiondelaytime.Underdifferentworking conditions(temperatureof90o-1200K,pressuresof0.5,1,1.5MPa,equivalentratiosof0.5,1.0,1.5), threewidely-used detailed mechanisms(San Diego 2016,GRI3.0,Aramco Mech 3.O)werecomparedandanalyzed.The Aramco Mech3.Omechanismwhichcould predictthe ignitiondelaytimeaccuratelywasdeterminedas theresearchobject,and thelargestsensitivitycoeficientofreactionwasidentifiedasthekeyfundamentalreactionbasedonanalysisofthesensitivityof ignitiondelaytime.FivequantitativecalculationmethodsB3LYP,PM7,Mo6-2X,CBS-QB3andW1ROwereusedtocalculate theactivationenergyoffundamentalreaction.TheactivationenergiescalculatedbyB3LYP,CBSQB3andWROwereselect ed,andtheactivationenergiesofkeyelementaryreactions intheoriginalmechanismwereupdated,andthreediferent new mechanismswerehenceobtained.Theignitiondelaytimeobtainedbycouplingcalculationofthethreemechanismsanddetailed mechanisms withtheclosed homogeneous reactormodel werecomparedwith the experimental values.Theresultsshow that the average activation energy calculated by B3LYP method is 212.34kJ/mol when the temperature is .Usingthe datatocorrecttheactivationenergyintheoriginalmechanism,thecorrctednewmechanismcanpredicttheignitiondelaytime more accurately,and is more suitable for low-temperature combustion conditions.

Key Words:natural gas;quantum chemical calculation;activation energy;chemical reaction mechanism

[編輯:袁曉燕]

主站蜘蛛池模板: 国产成人狂喷潮在线观看2345| 青青草国产在线视频| 成人免费一级片| 91视频精品| 午夜人性色福利无码视频在线观看| 国产福利免费观看| 2021精品国产自在现线看| 免费AV在线播放观看18禁强制| 麻豆国产在线观看一区二区| 国产呦视频免费视频在线观看| 日韩无码视频播放| 老司国产精品视频91| 国产免费久久精品99re丫丫一| 欧美成人免费午夜全| 91精品最新国内在线播放| 欧美日韩激情在线| 亚洲男人天堂2018| 久久久久青草线综合超碰| 国产97视频在线| 台湾AV国片精品女同性| 成人毛片在线播放| 99精品影院| 欧美日韩在线亚洲国产人| 国产最新无码专区在线| 九九精品在线观看| 亚洲AV无码不卡无码| AV不卡无码免费一区二区三区| 亚洲天堂日韩av电影| 午夜精品久久久久久久无码软件| 久久人午夜亚洲精品无码区| 2020国产在线视精品在| 四虎永久在线精品影院| 亚洲男人天堂久久| 91免费片| 国内精品视频区在线2021| 国产精品丝袜视频| 精品久久777| 亚洲成av人无码综合在线观看| 亚洲视频一区在线| 国产精品久久久久久久伊一| 亚洲AV色香蕉一区二区| 久夜色精品国产噜噜| 凹凸精品免费精品视频| 国产区福利小视频在线观看尤物| 国产91小视频在线观看| 青青草国产一区二区三区| 国产成人一区| 欧美日韩成人| 一区二区三区国产精品视频| av天堂最新版在线| 五月天丁香婷婷综合久久| 老熟妇喷水一区二区三区| 色综合天天娱乐综合网| 欧美成人综合在线| 国产成人永久免费视频| а∨天堂一区中文字幕| 无码综合天天久久综合网| 亚洲 欧美 偷自乱 图片| 凹凸国产分类在线观看| 亚洲天堂免费观看| 538国产在线| 国产精品原创不卡在线| 99热这里都是国产精品| 亚洲浓毛av| 国产亚洲欧美另类一区二区| 亚洲天堂免费在线视频| 亚洲经典在线中文字幕| 免费无码在线观看| 1769国产精品视频免费观看| 伊人久久综在合线亚洲91| 亚洲精品无码不卡在线播放| 亚洲欧美日韩另类在线一| 亚洲成人网在线播放| 久久精品娱乐亚洲领先| 欧美激情视频一区二区三区免费| 浮力影院国产第一页| 人妻精品久久久无码区色视| 日本黄色a视频| 久久成人国产精品免费软件| 国产精品九九视频| 欧美 国产 人人视频| 99久久亚洲精品影院|