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

碰撞能對H+CH+→C++H2反應立體動力學性質的影響?

2017-08-01 00:35:16唐曉平周燦華和小虎于東麒楊陽3
物理學報 2017年2期

唐曉平周燦華和小虎于東麒楊陽3)

1)(遼寧師范大學物理與電子技術學院,大連 116029)

2)(中國科學院大連化學物理研究所,分子反應動力學國家重點實驗室,大連 116023)

3)(中國科學院大連化學物理研究所,中國科學院化學激光重點實驗室,大連 116023)

碰撞能對H+CH+→C++H2反應立體動力學性質的影響?

唐曉平1)2)周燦華3)和小虎2)?于東麒1)?楊陽2)3)

1)(遼寧師范大學物理與電子技術學院,大連 116029)

2)(中國科學院大連化學物理研究所,分子反應動力學國家重點實驗室,大連 116023)

3)(中國科學院大連化學物理研究所,中國科學院化學激光重點實驗室,大連 116023)

(2016年8月23日收到;2016年10月31日收到修改稿)

準經典軌線方法,矢量相關,反應截面

1 引 言

近些年來由于在燃燒化學、煤氣化反應和天體化學方面扮演著重要的角色,含碳和氫的反應一直很受關注[1-3],C+和H2相碰撞可以形成CH+離子[4],CH+在星際空間中的含量很豐富[5],當與星際云中含量最豐富的H原子碰撞時發生抽取反應則生成C++H2,眾多的研究小組對體系的勢能面構建、反應的速率常數以及反應截面等方面進行了研究[5-15].2005年Stoecklin和Halvick[6]的研究小組基于多組態從頭算能量點擬合出了一個勢能面來研究基態的H+CH+反應.此勢能面曾被多個課題組所引用.Halvick等[7]基于這個勢能面運用準經典軌線法和相空間理論分析了H+CH+的反應,結果顯示在高能區時非彈性和交換兩個反應軌道中觀察到了動力學效應,非彈性反應截面隨碰撞能的增加而增大,抽取反應截面呈下降趨勢并且比理論預計結果下降得更快.當溫度為20-700 K時抽取反應的速率常數隨溫度的增加呈單調遞減趨勢,但是Plasil等[8]和Gerlich等[9]發現在溫度為60 K處速率常數有最大值,并且在低于60 K時速率常數的計算結果[6,7]與實驗值背道而馳.2011年,Warmbier和Schneider[10]發表了一篇關于基于從頭算能量點擬合出WS勢能面(以下簡記為WS勢能面)的論文,在此勢能面上運用準經典軌線法和非含時量子散射法計算了熱速率常數,發現抽取反應速率常數與實驗數據的結果在較低溫度50和100 K時呈現相同的下降趨勢,但是下降的幅度不同.基于WS勢能面[10],2014年Herráez-Aguilar等[11]使用準經典軌線法和高斯Binning算法研究了反應物的內在激發對C++H2反應的動力學性質,結果表明CH+在星際物質當中的含量異常豐富與反應物的振動激發和轉動激發有關.最近Bonfanti等[12]研究了最低的三個電子態勢能面以及在星際空間中的CH+的形成和破壞,最重要的結果是獲得一個可做全量子動力學計算的多功能絕熱勢能面,并幫助解決碳化學過程在星際介質中的一些問題.2015年Li等[13]構建了體系的勢能面(以下簡記為LYQ勢能面),并且運用含時量子方法對其動力學性質進行了計算,將反應截面隨碰撞能的變化曲線與2011年所得的結果[10]進行了對比,結果表明兩條曲線在20 meV之后基本處于重合狀態.同年Werfelli等[14]運用非含時量子散射法在新構建的勢能面上討論H+CH+→C++H2的反應在低溫條件下的速率常數,結果顯示溫度在50-800 K的范圍內理論值與實驗結果相符,而當低于50 K時實驗數值遠遠小于理論值.到目前為止,所有的研究都只關注在該反應體系的標量性質,并未發現前人任何關于其矢量性質的研究報道.為了更充分地揭示此體系的動力學性質,我們不應該局限于標題反應的標量屬性,也應該重視反應的矢量屬性,只有將標量性質和矢量性質結合起來才能更好地了解該反應的動力學圖像.準經典軌線是一種研究反應體系矢量性質的有效方法,盡管其計算代價較低,由于忽略了量子效應,對于反應概率、反應截面和反應速率常數等標量性質的計算結果不太準確,但是,對于計算矢量性質,例如速度分布和角動量分布等,則相對準確.因此,本文中我們基于LYQ勢能面[13]對反應立體動力學性質做了準經典軌線計算,以期能夠加深對該反應的動力學性質的認識.本文一共分四部分,第一節是引言,第二節主要對計算中涉及的立體動力學參數和準經典軌線的方法進行介紹,第三節是對結果的描述和討論,最后一節對整個文章的研究結果做了總結.

2 理論計算

2.1 矢量相關函數分布

我們采取質心坐標系描述反應物的相對速度k和產物相對速度k′的分布.坐標系的Z軸正方向平行于反應物相對速度矢量k的方向,Y軸垂直于含有反應物相對速度矢量k和產物相對速度矢量k′的X-Z平面(該平面為散射面).k和k′的夾角θt為散射角,產物轉動角動量j′的極角和方位角分別為θr和φr,見圖1.

描述產物分子k-j′兩矢量相關的分布函數P(θr)可以展開為Legendre[16]多項式

圖1 描述k,k′和j′分布的質心坐標系Fig.1.The center-of-mass coordinate system used to describe thek,k′andj′correlations.

描述k,k′和j′三矢量相關的函數P(φr)可以用Fourier級數[17]展開

其中,Ckq(θr,φr)是修正的球諧函數.

聯系k,k′和j′三矢量的角分布函數可以寫為[17]

其中,[k]=2k+1;ωt=θt,φt;ωr=θr,φr;σ表示積分截面;Ckq(θr,φr)是修正的球諧函數;是廣義極化微分反應截面.在很多光誘導的雙分子反應實驗中,人們對k=0的極化分量感興趣,我們只計算了兩個極化微分反應截面:(2π/σ)(dσ00/dωt) 和(2π/σ)(dσ20/dωt), 在計算中展開到k=7時就可以得到較好的收斂結果.

2.2 準經典軌線的計算

本文使用準經典軌線法基于LYQ勢能面,研究了碰撞能對H(2S)+CH+(X1Σ+)→反應的反應截面以及立體動力學性質的影響.當計算反應截面時,將碰撞能取為1,2,···,10,20,30,···,100,200,300,···,1000 meV共計28個能量點.在計算反應物分子立體動力學性質時碰撞能分別取為1,10,100,500,1000 meV共計五個能量點.總的軌線條數為50000.H原子和CH+離子質心間的初始距離取為15 ?,它附近的勢能約為10-10eV能量級,這個勢能值非常小,因此相互作用力微弱,對動力學過程的影響可以忽略.當把質心距離增大到50 ?時,發現勢能沒有變化,這也進一步說明了CH+離子質心間的初始距離為15 ?時已經可以認為是解離狀態了.積分步長為0.1 fs.反應的最大碰撞參數bmax的確定方法是:先選擇5000條軌線運行,初步確定bmax的范圍,然后再用50000條軌線進行運行,逐漸增加b使其反應的軌線條數不增加即可[18].

3 結果與討論

H(2S)+CH+(X1Σ+)→反應,是一個放熱反應.放熱量為0.496 eV,并且存在一個勢阱,勢阱深度超過4 eV.在表1中介紹了WS[10]和LYQ勢能面[13]的放熱量和分子解離能與相應實驗值之間的比較.通過對比發現,LYQ勢能面[13]不論是放熱量還是解離能都比WS勢能面[10]更接近于實驗值,表明LYQ勢能面[13]的精確度更高.

圖2表示的是反應截面隨碰撞能的變化情況,并與文獻中的含時量子計算結果[13]以及非含時量子散射方法計算的結果[10]做了比較.從圖上可以清晰看出,隨著碰撞能的增加三條曲線的反應截面都呈下降趨勢,當前的結果與Li等[13]計算的結果基本符合,在低于20 meV時兩條曲線稍有差距,原因可能是由于準經典軌線法自身的局限性,不能處理零點能效應,從而造成結果的偏差.三條曲線總體上符合得很好,這證明使用準經典軌線法對H+CH+反應進行的動力學研究是可靠的.

表1 CH體系的一些相對能量值,表中對比了WS勢能面[10],LYQ勢能面[13]和實驗值的偏差.De代表解離能,ΔDe代表兩個雙原子分子的解離能的差.所有能量的單位是eVTable 1.Relevant data of the energetic for the CHsystem,the differences between WS PES[10],LYQ PES[13]and experiment values are shown.Deis represents the dissociation energies,and ΔDerepresents the difference of the two dissociation energies.Energy is in eV.

表1 CH體系的一些相對能量值,表中對比了WS勢能面[10],LYQ勢能面[13]和實驗值的偏差.De代表解離能,ΔDe代表兩個雙原子分子的解離能的差.所有能量的單位是eVTable 1.Relevant data of the energetic for the CHsystem,the differences between WS PES[10],LYQ PES[13]and experiment values are shown.Deis represents the dissociation energies,and ΔDerepresents the difference of the two dissociation energies.Energy is in eV.

性質Property WS勢能面WS[10]/eV LYQ勢能面LYQ[13]/eV實驗值Exp./eV解離能的差 ΔDe-0.518-0.496-0.496解離能De(H2) 4.711 4.748 4.751[19]解離能De(CH+) 4.195 4.252 4.255[6]

圖2 (網刊彩色)H+CH+(v=0,j=0)→C++H2反應截面隨碰撞能的變化This work表示當前的工作,TD表示Li等的結果[13],TI表示基于WS勢能面的結果[10],E表示碰撞能Fig.2.(color online)Reactive cross sections of the H+CH+(v=0,j=0)→C++H2reaction as a function of the collision energy.TW gives the present job,TD gives the result of Liet al.[13],TI gives the result based on the WS potential energy surface[10],Erepresent the collision energy.

圖3表示的是不同的碰撞能下(E=1,10,100,500,1000 meV),H+CH+→C++H2反應的k和j′兩矢量相關的P(θr)的分布.從圖中我們可以看出,五個碰撞能對應的函數P(θr)在θr=90°處有一個明顯的峰值,并且關于θr=90°呈軸對稱分布,這表明矢量j′的取向分布傾向Y軸的方向.當碰撞能由1 meV增加到500 meV時P(θr)的峰變低,寬度基本沒有變化,這說明隨碰撞能的增加取向逐漸變弱.但是當碰撞能達到1000 meV時θr=90°處峰值增高,顯然碰撞能由500 meV到1000 meV時取向逐漸增強.可以觀察到:低能區時產物轉動角動量的取向隨碰撞能的增加呈減弱趨勢,高能區時隨碰撞能的增加取向增強.出現這種情況,可能是勢能面上存在較深的勢阱所致.

圖3 (網刊彩色)不同的碰撞能H+CH+(v=0,j=0)→C++H2反應的P(θr)分布,E表示碰撞能Fig.3.(color online)Angular distribution ofP(θr)atfive collision energies of the H+CH+(v=0,j=0)→C++H2reaction,Erepresent the collision energy.

圖4表示的是不同的碰撞能下(E=1,10,100,500,1000 meV),H+CH+→C++H2反應的k,k′和j′三矢量相關的P(φr)的分布.從圖上我們可以看出不同的碰撞能對應的P(φr)分布關于φr=180°不對稱,這說明產物分子的轉動角動量分布有強烈的極化.φr=90°和φr=270°時都有峰值,但是二者的峰值不同,在φr=90°處的峰值要大于在φr=270°的峰值,這表明了產物轉動角動量j′不僅沿Y軸有取向分布,還定向分布于Y軸的正方向.此外還可以觀察到,在φr=90°和φr=270°兩處的P(φr)分布峰值在碰撞能從1 meV增加到100 meV時,呈下降趨勢,而當碰撞能從100 meV增加到1000 meV時,峰值呈上升趨勢,同時還發現φr=90°處的峰值變化更大,這表明產物轉動角動量沿著Y軸的正定向分布在低能區隨碰撞能增加而減弱,在高能區則隨碰撞能增加而增強.這種分布情況可以使用三原子反應的排斥模型進行解釋,依據瞬時碰撞模型[20],產物轉動角動量j′可以表示成j′=Lsinβ2+jcosβ2+J1mB/mAB. 其中,L和j分別為反應物的軌道角動量和轉動角動量;和rCB分別是由B原子指向A原子和C原子的單位矢量,μBC是BC分子的約化質量,R為排斥能.在化學鍵斷開和重新形成的過程中,Lsinβ2+jcosβ2是對稱的,但是由于排斥能的影響使J1mB/mAB項不對稱,更加傾向于某個方向從而導致了產物轉動角動量矢量的定向效應[21].

圖4 (網刊彩色)不同的碰撞能下H+CH+(v=0,j=0)→C++H2反應的P(φr)分布,E表示碰撞能Fig.4.(color online)Angular distribution ofP(φr)atfive collision energies of the H+CH+(v=0,j=0)→C++H2reactions,Erepresent the collision energy.

為了更好地描述出反應的立體動力學信息,圖5給出反應產物的空間分布P(θr,φr)的立體圖形.圖5(a)-圖5(e)依次對應E=1,10,100,500,1000 meV下碰撞能的空間分布情況.為了方便比較,圖中的概率分布的顯示范圍全部調整為0.04-0.18.如圖所示:在θr=90°和φr=270°處P(θr,φr)有明顯的峰值,并且二者峰值大小不同,前者的峰值要高于后者.由圖5(a)-圖5(e)我們可以看出在低能區圖5(a)-圖5(c)時隨著碰撞能的增加峰值高度在逐漸變小,而在更高的碰撞能區域中峰值隨碰撞能增加而增高.反應產物的空間分布P(θr,φr) 與圖3 中的P(θr)和圖4中的P(φr)結果是一致的.

圖5 (網刊彩色)不同的碰撞能下H+CH+(v=0,j=0)→C++H2反應的空間分布P(θr,φr)(a)E=1 meV;(b)E=10 meV;(c)E=100 meV;(d)E=500 meV;(e)E=1000 meVFig.5.(color online)Spatial distribution ofP(θr,φr)at five collision energies of the H+CH+(v=0,j=0)→C++H2reaction:(a)E=1 meV;(b)E=10 meV;(c)E=100 meV;(d)E=500 meV;(e)E=1000 meV.

圖6 (網刊彩色)不同的碰撞能下H+CH+(v=0,j=0)→C++H2反應的兩個極化微分反應截面,E表示碰撞能Fig.6.(color online)Two polarization dependent differential cross-sections of the reaction H+CH+(v=0,j=0)→C++H2at five collision energy,Erepresent the collision energy.

圖6表示的是不同的碰撞能下(E=1,10,100,500,1000 meV),H+CH+→C++H2反應的極化微分反應截面的分布情況.從圖6(a)可以看出產物分子有明顯的前向和后向散射,(2π/σ)(dσ00/dωt)只與反應物相對速度矢量和產物相對速度矢量有關,隨著碰撞能的增加前向和后向散射程度依次增強,當碰撞能從500 meV增加到1000 meV時前向散射的強度減弱而后向散射依然增強.這反映了低能區與高能區下能量的變化對反應的產物角動量影響有很大的不同,在高能區所受到碰撞能的影響很小,而在低能區下,由于能量接近反應物跨過勢壘所需的能量,反應過程中發生的彈性碰撞概率增加,且平動能較小,因而出現了此種情況.圖6(b)中極化微分反應截面(2π/σ)(dσ20/dωt)不僅與θt有關,同時還是轉動取向因子的函數,其變化的趨勢與(2π/σ)(dσ00/dωt)的趨勢相反.(P2(cosθr))的期望值是負的,因此產物分子的轉動角動量j′在垂直于k的方向上有強烈的取向.觀察圖6(b)發現分布在散射角末端的產物其轉動角動量的取向要強于其他方向上的產物.

4 結 論

本文采用準經典軌線法,基于LYQ勢能面[13]研究了不同的碰撞能下反應截面以及立體動力學性質.計算發現,反應截面隨著碰撞能的增加逐漸降低,并且與文獻[10,13]計算的結果做了比較,結果顯示三組符合得非常好,從而驗證了在此勢能面上我們利用準經典軌線法求得的結果是可靠的.通過取不同的碰撞能值,該反應的兩矢量、三矢量分布都有明顯的變化,由此可見碰撞能對其影響很大.同樣極化微分反應截面也隨碰撞能的變化有所改變,(2π/σ)(dσ00/dωt)隨著碰撞能的增加前向散射增強,然而當E=1000 meV時前向散射隨著碰撞能的增加而減弱,后向散射依然保持增強的趨勢.(2π/σ)(dσ20/dωt)顯示的變化趨勢表明分布在散射角末端的產物其角動量分布具有更強的取向效應.綜上所述,此反應的立體動力學性質對碰撞能有很強的依賴性.

[1]Langer W 1978Astrophys.J.225 860

[2]Draine B T 1986Astrophys.J.310 408

[3]Zanchet A,Godard B,Bulut N,Roncero O,Halvick P,Cernicharo J 2013Astrophys.J.766 80

[4]Lique F,Werfelli G,Halvick P,Stoecklin T,Faure A,Wiesenfeld L,Dagdigian P J 2013J.Chem.Phys.138 204314

[5]Ervin K M,Armentrout P B 1986J.Chem.Phys.84 6738

[6]Stoecklin T,Halvick P 2005Phys.Chem.Chem.Phys.7 2446

[7]Halvick P,Stoecklin T,Larrégaray P,Bonnet L 2007Phys.Chem.Chem.Phys.9 582

[8]Plasil R,Mehner T,Dohnal P,Kotrik T,Glosik J,Gerlich D 2011Astrophys.J.737 60

[9]Gerlich D,Disch R,Scherbarth S 1987J.Chem.Phys.87 350

[10]Warmbier R,Schneider R 2011Phys.Chem.Chem.Phys.13 10285

[11]Herráez-Aguilar D,Jambrina P,Menéndez M,Aldegunde J,Warmbier R,Aoiz F 2014Phys.Chem.Chem.Phys.16 24800

[12]Bonfanti M,Tantardini G F,Martinazzo R 2014J.Chem.Phys.A118 6595

[13]Li Y Q,Zhang P Y,Han K L 2015J.Chem.Phys.142 124302

[14]Werfelli G,Halvick P,Honvault P,Kerkeni B,Stoecklin T 2015J.Chem.Phys.143 114304

[15]Grozdanov T,McCarroll R 2013Chem.Phys.Lett.575 23

[16]Chen M D,Han K L,Lou N Q 2003J.Chem.Phys.118 4463

[17]Aoiz F,Brouard M,Enriquez P 1996J.Chem.Phys.105 4964

[18]Wu V W K 2011Phys.Chem.Chem.Phys.13 9407

[19]Balakrishnan A,Smith V,Stoiche ffB 1992Phys.Rev.Lett.68 2149

[20]Han K L,He G Z,Lou N Q 1998Chin.J.Chem.Phys.11 525(in Chinese)[韓克利,何國忠,樓南泉1998化學物理學報11 525]

[21]Kong H,Liu X G,Xu W W,Zhang Q G 2009Acta Phys.-Chim.Sin.25 935(in Chinese)[孔浩,劉新國,許文武,張慶剛2009物理化學學報25 935]

PACS:34.50.Lf,31.15.xv,87.15.H- DOI:10.7498/aps.66.023401

Influence of collision energy on the stereodynamics of the H+CH+→C++H2reaction?

Tang Xiao-Ping1)2)Zhou Can-Hua3)He Xiao-Hu2)?Yu Dong-Qi1)?Yang Yang2)3)

1)(School of Physics and Electronic Technology,Liaoning Normal University,Dalian 116029,China)
2)(State Key Laboratory of Molecular Reaction Dynamics,Dalian Institute of Chemical Physics,Chinese Academy of Sciences.Dalian 116023,China)
3)(Key Laboratory of Chemical Lasers,Dalian Institute of Chemical Physics,Chinese Academy of Sciences,
Dalian 116023,China)

23 August 2016;revised manuscript

31 October 2016)

The reactive cross section and stereodynamics at selected collision energies for the H(2S)+CH+(X1Σ+)→reaction on a globally smoothabinitio potential surface of the 2A′state are calculated in detail by the quasi-classical trajectory(QCT)method.The calculated cross section decreases with the increase of the collision energy,which is found to be in overall good agreement with the previous time-dependent quantum results in the high collision energy regime(Ec>20 meV).The discrepancy between the QCT and previous quantum cross section below 20 meV can be attributed to the limitations of the classical trajectory method,because the QCT method cannot handle the effect of zero point energy.In general,QCT results show qualitative agreement with the quantum results,which confirmsthe validity of the QCT method.The research shows that the product rotational angular momentum vector is aligned and oriented.The alignment of the product rotational angular momentum vectorj′depends very sensitively on the collision energy.With the increase of the collision energy,the alignment effect recedesin the low collision energy region(1500 meV),while it is enhanced in the high collision energy region(500-1000 meV).Moreover,thek-k′-j′distributions tend to be asymmetric with respect to thek-k′scattering plane(or aboutφr=180?),with two peaks appearing atφr=90?andφr=270?,respectively.This indicates that the product rotational angular momentum is not only in theY-axis direction but also along the positiveY-axis direction.The peak intensity decreases with the collision energy increasing from 1 meV to 100 meV,while it increases with collision energy increasing from 100 meV to 1000 meV.Therefore theY-axis orientation effect turns weak with the enhancement of the collision energy in the low energy region,while it becomes strong in the high energy region.In addition,the polarization dependent differential cross sections(PDDCSs)(2π/σ)(dσ00/dωt)and(2π/σ)(dσ20/dωt)are calculated.PDDCS(2π/σ)(dσ00/dωt)results indicate that the products have almost symmetrically scattered forward and backward,and the intensity of the scattering increases with the increase of the collision energy.The PDDCS(2π/σ)(dσ20/dωt)shows that the alignment effect of the rotational angular momentum of the products is stronger at the terminal of the scattering angle than at the other directions.

quasi-classical trajectory method,vector correlation,reactive cross section

:34.50.Lf,31.15.xv,87.15.H-

10.7498/aps.66.023401

?國家自然科學基金(批準號:21403226,21503226)資助的課題.

?通信作者.E-mail:hxh@dicp.ac.cn

?通信作者.E-mail:useeu@163.com

*Project supported by the National Natural Science Foundation of China(Grant Nos.21403226,21503226).

?Corresponding author.E-mail:hxh@dicp.ac.cn

? Corresponding author.E-mail:useeu@163.com

主站蜘蛛池模板: 欧美亚洲日韩中文| 国产黄在线免费观看| 精品无码日韩国产不卡av| 久久亚洲美女精品国产精品| 日韩精品久久无码中文字幕色欲| 又猛又黄又爽无遮挡的视频网站| 亚洲国产午夜精华无码福利| 无码一区中文字幕| 国产激情影院| 9啪在线视频| 亚洲永久精品ww47国产| 4虎影视国产在线观看精品| 91麻豆精品视频| 天天综合网在线| 三上悠亚在线精品二区| 超碰aⅴ人人做人人爽欧美 | 婷婷五月在线| 日韩美一区二区| 精品国产美女福到在线不卡f| 欧美成人h精品网站| 激情五月婷婷综合网| 亚洲一区国色天香| 国产精品毛片在线直播完整版| 国产成人成人一区二区| 欧美人人干| 久久黄色免费电影| 久久综合久久鬼| 色婷婷综合在线| 凹凸精品免费精品视频| 精品偷拍一区二区| 日本三区视频| 欧美国产另类| 女人18一级毛片免费观看| 国产成人亚洲欧美激情| 国产成人久视频免费| 亚洲欧美不卡视频| 原味小视频在线www国产| 露脸一二三区国语对白| 尤物精品视频一区二区三区| 国产在线观看人成激情视频| 欧美日韩高清在线| 国产青榴视频在线观看网站| 日韩AV无码一区| 国产黑丝视频在线观看| 国产乱人伦精品一区二区| 国产成人精品高清不卡在线| 国产一区成人| 亚洲无码视频一区二区三区| 欧美人在线一区二区三区| 欧美一级99在线观看国产| 97影院午夜在线观看视频| 69视频国产| 亚洲成年人片| 国产综合色在线视频播放线视| 日韩大乳视频中文字幕| 欧美亚洲一区二区三区在线| 日韩在线欧美在线| 五月天在线网站| 国产亚洲精品yxsp| 91青青视频| 久久精品丝袜| 色妞永久免费视频| 9999在线视频| 一级黄色网站在线免费看| 国产精品播放| 亚洲日本中文字幕乱码中文| 日本免费高清一区| 久久国产精品麻豆系列| 99精品国产自在现线观看| 另类综合视频| 国产主播喷水| 99视频在线免费| 中文无码伦av中文字幕| 亚洲一区二区成人| 麻豆精品久久久久久久99蜜桃| 久久久国产精品免费视频| 国产精品网曝门免费视频| 久久久久国产一级毛片高清板| 四虎永久在线精品国产免费| 成人一级免费视频| 91精品国产麻豆国产自产在线| 亚洲成人精品在线|