亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        D + DBr反應(yīng)的態(tài)-態(tài)動力學(xué)研究

        2020-06-30 12:12:24袁方園朱子亮
        物理學(xué)報 2020年11期
        關(guān)鍵詞:勢阱波包勢能

        袁方園 朱子亮

        (濰坊科技學(xué)院通識學(xué)院, 壽光 262700)

        (2020 年3 月2日收到; 2020 年4 月2日收到修改稿)

        在0—2.0 eV的碰撞能范圍內(nèi)采用含時量子波包方法和二階分裂算符傳播子對初始量子態(tài)為(v0 = 0, j0 = 0)的D + DBr反應(yīng)進行了態(tài)分辨理論水平上的動力學(xué)計算. 計算了抽取反應(yīng)通道D + DBr → Br + D2和置換反應(yīng)通道 D′+DBr→D+D′Br 的反應(yīng)概率、積分截面、微分截面、產(chǎn)物的振動和轉(zhuǎn)動分布以及速率常數(shù)等動力學(xué)性質(zhì), 并與相應(yīng)的理論和實驗結(jié)果進行了比較. 結(jié)果表明: 計算的速率常數(shù)與實驗結(jié)果十分符合. 微分截面的結(jié)果表明: 對于抽取反應(yīng), 前向的剝離反應(yīng)機制在反應(yīng)過程中占據(jù)主導(dǎo)地位; 對于置換反應(yīng), 頭碰頭的反彈機制占據(jù)主導(dǎo)地位.

        1 引 言

        氫原子和鹵化氫的反應(yīng)在發(fā)展和理解基礎(chǔ)的化學(xué)反應(yīng)上有著重要意義, 因此其抽取和置換反應(yīng)在最近幾十年當(dāng)中受到了廣泛的關(guān)注. 對于H + HBr及其同位素反應(yīng), 大量的實驗工作集中在測量該體系的熱速率常數(shù)、同位素效應(yīng)、積分截面、以及產(chǎn)物的振動和轉(zhuǎn)動電子態(tài)分布等動力學(xué)信息[1-10].1989年, Aker等[7,8]采用連貫的反斯托克斯拉曼散射光譜測量了1.6 eV能量下H + HBr反應(yīng)抽取和置換反應(yīng)通道的積分截面和產(chǎn)物的振動和轉(zhuǎn)動態(tài)分布. 1990年, Umemoto等[6]采用一種脈沖輻射共振吸附技術(shù)在214—300 K的溫度范圍內(nèi)測量了H + HBr及其同位素反應(yīng)的速率常數(shù). 采用激光閃光光解HBr分子和連續(xù)波共振熒光檢測H原子方法, Seakins和Pilling[3]與Talukdar等[4]也分別測量了H + HBr反應(yīng)不同溫度范圍的速率常數(shù).

        在理論方面, H + HBr反應(yīng)的勢能面構(gòu)建也在最近的幾十年中有了極大的進展. 基于多參考組態(tài)相互作用方法(MRCI)計算的104個從頭算能量點, Lynch等[11]構(gòu)建了第一個extended London-Eyring-Polanyi-Sato (eLEPS)函數(shù)形式的從頭算勢能面. 然而基于該勢能面得到的速率常數(shù)往往低于實驗值并且100多個能量點也不足以描述復(fù)雜的勢能面地貌. 2003年, Kurosaki和Tayanagi[12]采用MRCI方法和AVTZ基組構(gòu)建了BrH2體系基態(tài)和兩個激發(fā)態(tài)的勢能面, 并且在基態(tài)的勢能面上考慮了自旋軌道耦合修正(MB1). 然而基于MB1勢能面的速率常數(shù)明顯大于實驗值. 不久之后, 他們又對MB1勢能面進行了修正, 降低了約0.3 kcal/mol的勢壘高度, 并將新的勢能面命名為MB2[13]. 之后又在MB2勢能面的基礎(chǔ)上修改了Br + H2漸進區(qū)T型結(jié)構(gòu)的勢阱, 并命名為MB3.Kurosaki和Tayanagi[12]的勢能面得到了廣泛的應(yīng)用. Quan等[14-16]采用MB1, MB2和MB3勢能面和含時量子波包方法進行了一系列的動力學(xué)計算并報道了反應(yīng)概率、積分截面和速率常數(shù)等物理量. 結(jié)果表明對于H + HBr反應(yīng), MB3勢能面明顯優(yōu)于MB1和MB2勢能面, 基于MB3勢能面的速率常數(shù)和實驗結(jié)果十分吻合. Fu和Zhang[17]和Fu等[18]也采用MB3勢能面和含時量子波包方法進行了動力學(xué)計算, 并報道了量子態(tài)分辨的積分和微分截面. 然而, Fu和Zhang[17]基于MB3勢能面得到的速率常數(shù)明顯大于實驗值. 這與Quan等[14-16]得到的結(jié)果有很大的差異.

        2011年, Jiang等[19]構(gòu)建了H + HBr反應(yīng)的一個新的勢能面. 采用MRCI方法進行了從頭算計算, 在計算中考慮了自旋軌道耦合修正并把基組外推到了完全基組極限. 此外, 基于新構(gòu)建的勢能面, Jiang等[19]采用含時量子波包方法進行了動力學(xué)計算. 結(jié)果表明基于該勢能面得到的速率常數(shù)與實驗結(jié)果十分吻合. 之后, Xie等[20]采用該勢能面進行了H + HBr反應(yīng)的態(tài)-態(tài)動力學(xué)計算, 報道了反應(yīng)概率、積分截面、微分截面、產(chǎn)物的振動和轉(zhuǎn)動態(tài)分布等動力學(xué)信息, 并與相應(yīng)的實驗結(jié)果進行了對比. 結(jié)果表明, 基于該勢能面得到的產(chǎn)物振動和轉(zhuǎn)動態(tài)分布與實驗值存在較大差距. 同年,Jiang等[21]又構(gòu)建了關(guān)于BrH2體系的透熱勢能面并計算了Br(2p3/2,2p1/2) + H2反應(yīng)的積分截面.2012年, Xie等[22]采用Jiang等[21]構(gòu)建的透熱勢能面在量子態(tài)分辨的理論水平上計算了Br(2p1/2) +H2反應(yīng)的猝滅過程. 2012年Zhang等[23]也構(gòu)建了BrH2體系的透熱勢能面并與Jiang等[21]的透熱勢能面進行了比較, 同時對D + DBr反應(yīng)進行了動力學(xué)計算, 報道了積分截面和速率常數(shù)等動力學(xué)信息. 2014年, Zhang等[24]采用他們構(gòu)建的透熱勢能面進行了H + HBr反應(yīng)體系的同位素計算. 將得到速率常數(shù)與實驗結(jié)果進行了對比, 發(fā)現(xiàn)理論結(jié)果與實驗結(jié)果有或高或低的偏差. 最近的勢能面是由Li等[25]在2019年構(gòu)建的, 在從頭算計算中采用MRCI方法、AVTZ和AVQZ基組計算了約63000個從頭算能量點, 并采用神經(jīng)網(wǎng)絡(luò)方法對勢能面進行了擬合. 勢能面的構(gòu)建過程中, 考慮了自旋軌道耦合修正并將基組外推到了完全基組極限. 此外, 基于新構(gòu)建的勢能面, Li等[25]采用含時量子波包方法進行了H + HBr反應(yīng)的動力學(xué)計算. 得到的抽取通道的產(chǎn)物振轉(zhuǎn)態(tài)分布以及速率常數(shù)都與實驗結(jié)果十分符合.

        綜上, 之前的理論研究大部分集中在H + HBr反應(yīng)上, 對于同位素的理論研究相對較少. 據(jù)我們所知, 到目前為止還沒有關(guān)于其同位素反應(yīng)在量子態(tài)分辨理論層次上的研究. 基于Li等的勢能面, 本文采用含時量子波包法和二階分裂算符傳播子對D + DBr反應(yīng)進行初態(tài)為(v0= 0, j0= 0)的動力學(xué)計算.

        2 理論方法

        含時量子波包方法在計算原子與雙原子的反應(yīng)散射以及包含多原子的散射計算中有著廣泛的應(yīng)用[26-30]. 這里, 僅簡單的介紹一下原子與雙原子分子反應(yīng)散射的動力學(xué)方法. 在空間坐標(biāo)系下, 采用反應(yīng)物坐標(biāo), 對于給定的總角動量J, D + DBr反應(yīng)體系的哈密頓量可以寫為

        在進行動力學(xué)計算之前, 具有確定初始振動和轉(zhuǎn)動量子態(tài)的初始波包需要提前構(gòu)建. 初始的含時波包是由體坐標(biāo)系下的平動、振動和轉(zhuǎn)動基矢構(gòu)成, 初始含時波包可以寫為

        其中, 下標(biāo)n為平動基矢,? 是體系的宇稱,(υ0j0K0)是體系的初始量子態(tài), M是總角動量在空間坐標(biāo)系下在z軸的投影.

        在傳播過程中, 二階分裂算符[31]進行初始波包的傳播. 假設(shè)t時刻的波函數(shù)為 ψ (R,r,t) , 那么經(jīng) 過時間步長 Δt, 波函數(shù)變?yōu)?ψ(R,r,t+Δt)

        在計算中, 采用快速傅里葉變換[32]計算徑向動能算符, 在離散變量表象[33]里計算勢能算符. 此外, 為了滿足邊界條件, 在邊界處設(shè)置一個吸收勢,其形式為

        最終, 通過反應(yīng)物坐標(biāo)變換方法(RCB)[34,35]來提取態(tài)-態(tài)散射信息. 利用態(tài)-態(tài)散射矩陣可以得到反應(yīng)概率、積分截面和微分截面等動力學(xué)信息.

        其 中 ψ (E) 是 通 過 ψ (t) 的 快 速 傅 里 葉 變 換 得 到,是 散 射 矩 陣,是 轉(zhuǎn) 動 矩陣,? 是散射角.

        指定初始量子態(tài)的速率常數(shù)是通過熱平均積分 截面的平動能得到, 其公式為

        其 中 kB是玻爾茲曼常數(shù).

        3 結(jié)果和討論

        本文的動力學(xué)計算是基于2019年Li等[25]報道的勢能面. 圖1給出了D + DBr反應(yīng)抽取反應(yīng)通道和置換反應(yīng)通道的最小能量路徑. 由圖1可知, 在抽取反應(yīng)通道入口處有一個微小的勢壘, 其高度約為1.54 kcal/mol. 該反應(yīng)是一個放熱反應(yīng),其放熱能約為18.77 kcal/mol, 并且在反應(yīng)路徑上沒有明顯的勢阱. 對于置換反應(yīng)通道, 在反應(yīng)路徑上有一個高度為12.14 kcal/mol的勢壘.

        在進行動力學(xué)計算之前需要對計算參數(shù)進行調(diào)試. 本文測試了不同參數(shù)J = 0的反應(yīng)概率, 并最終得到了最優(yōu)的收斂參數(shù). 計算中使用的參數(shù)如下:

        投影截面位置 R =12.0 ;

        總角動量J在z軸投影數(shù)目 Kmax=9 ;

        總的傳播時間 Ttot=13000 ;

        傳播時間步長 Tstep=10 .

        圖 1 勢能面的最小能量路徑 (a) 抽取反應(yīng)通道; (b) 置換反應(yīng)通道Fig. 1. Minimum energy path of the potential energy surface: (a) Abstraction channel; (b) exchange channel.

        3.1 反應(yīng)概率

        在0—2.0 eV的碰撞能范圍內(nèi), 總角動量J =0, 20, 40, 60, 80, 90的反應(yīng)概率見圖2. 通過Li等[25]的勢能面可知, 對于抽取反應(yīng)和置換反應(yīng)對應(yīng)的勢壘高度分別為1.54 kcal/mol (約0.067 eV)和12.14 kcal/mol (約0.525 eV). 由圖2可知, 對于抽取反應(yīng), 其閾值接近為0; 對于置換反應(yīng), 其閾值約為0.47 eV; 明顯低于勢能面上的勢壘高度, 這是由于無論是抽取反應(yīng)還是置換反應(yīng)都發(fā)生了隧穿效應(yīng), 從而降低了反應(yīng)閾值. 對于抽取反應(yīng), 所有J的反應(yīng)概率的外形都很類似. 都是先快速升高當(dāng)峰值接近0.1左右之后隨著碰撞能的升高反應(yīng)概率降低. 這是因為當(dāng)能量較低的時候只能通過隧穿效應(yīng)發(fā)生反應(yīng), 隨著能量升高, 隧穿效率增加反應(yīng)概率增高. 而當(dāng)碰撞能大于勢壘高度, 巨大的放熱能使得D原子的速度變快, DBr分子來不及與D原子形成化學(xué)鍵. 因此, 當(dāng)碰撞能升高時, 反應(yīng)概率反而降低, 這也是放熱反應(yīng)的典型特征. 對于置換反應(yīng), 反應(yīng)概率隨著碰撞能的升高而增高, 計算的最大J = 90, 由圖2可知, 在研究的碰撞能范圍內(nèi)所有的結(jié)果都是收斂的. 然而, 對于抽取反應(yīng)只有在碰撞能低于1.2 eV時, 得到積分截面和微分截面的結(jié)果才是收斂的.

        圖 2 若干總角動量J的反應(yīng)概率 (a) 抽取反應(yīng); (b) 置換反應(yīng)Fig. 2. The reaction probabilities of several selected total angular momentum J: (a) Abstraction reaction; (b) exchange reaction.

        3.2 積分截面

        抽取反應(yīng)的總積分截面和振動態(tài)分辨的積分截面分別列入圖3(a)和圖3(b). 為了與其他理論結(jié)果比較, Zhang等[24]的結(jié)果也列入到了圖3(a).如圖3(a)所示, 本文的結(jié)果與Zhang等[24]的結(jié)果十分符合尤其是當(dāng)碰撞能比較低的時候, 然而隨著碰撞能的增加, 本文的結(jié)果與Zhang等[24]的結(jié)果差距逐漸變大. 這可能是由于計算中使用了不同的勢能面. 本文采用的是2019年Li等[25]最新報道的絕熱勢能面, 該勢能面沒有考慮不同電子態(tài)之間的非絕熱效應(yīng). 而Zhang等[24]結(jié)果是基于他們自己構(gòu)建的非絕熱勢能面. 我們認(rèn)為高能部分的差距是由于非絕熱效應(yīng)引起的. 圖3(b)列入了積分截面的振動分辨, 發(fā)現(xiàn)了一個有趣的現(xiàn)象, 當(dāng)碰撞能低于0.45 eV時, 出現(xiàn)了振動粒子數(shù)反轉(zhuǎn)的現(xiàn)象.產(chǎn)物D2的分布主要集中在v' = 1上. 隨著碰撞能的升高, 反轉(zhuǎn)現(xiàn)象變得不明顯進而消失. 這種反轉(zhuǎn)分布符合Polanyi的規(guī)則[36,37]. 隨著碰撞能的升高,其他振動激發(fā)態(tài)的通道逐漸打開. 盡管如此, D2分子的振動態(tài)分布主要還是集中在基態(tài)和第一激發(fā)態(tài)上.

        圖 3 (a) 抽取反應(yīng)總積分截面和文獻[24]的結(jié)果比較;(b) 抽取反應(yīng)振動分辨的積分截面Fig. 3. (a) Comparison between present total integral cross section of abstraction reaction and the values obtained from Ref. [24]; (b) the vibrational state-resolved integral cross section of abstraction reaction.

        圖 4 在若干選擇的碰撞能下, 產(chǎn)物D2的振轉(zhuǎn)態(tài)分布Fig. 4. The ro-vibrational distribution of D2 at several selected collision energies.

        圖4 給出了抽取反應(yīng)在碰撞能為0.4, 0.6,0.8和1.0 eV時產(chǎn)物的振轉(zhuǎn)態(tài)分布. 由圖4可知,隨著碰撞能的升高, 振動電子態(tài)的數(shù)目和轉(zhuǎn)動量子數(shù)的數(shù)目都增加. 這反映了隨著碰撞能的升高會有更多的碰撞能轉(zhuǎn)化為內(nèi)能. 此外, 還可以看出越高的振動態(tài)具有越少的轉(zhuǎn)動量子數(shù), 這是因為當(dāng)總能固定的情況下, 內(nèi)部的轉(zhuǎn)動能量轉(zhuǎn)化為振動能量.我們還發(fā)現(xiàn)隨著碰撞能的升高, 出現(xiàn)了共振峰結(jié)構(gòu). 例如, 當(dāng)0.4 eV時, v' = 1是單峰結(jié)構(gòu), 0.6 eV時變?yōu)殡p峰, 0.8 eV時變?yōu)槿澹?1.0 eV時變?yōu)槎喾褰Y(jié)構(gòu). 在圖1, 抽取反應(yīng)的最小能量路徑上沒有發(fā)現(xiàn)明顯的勢阱. 本文認(rèn)為, 盡管最小能量路徑上沒有勢阱, 然而當(dāng)產(chǎn)物D2分子處于振動激發(fā)態(tài)的時候, 就會形成勢阱, 并且隨著振動態(tài)數(shù)目的變大形成的勢阱也會變深, 相應(yīng)的勢阱支持的束縛態(tài)和準(zhǔn)束縛態(tài)的壽命會變長并且不同電子態(tài)之間的波函數(shù)耦合也會變得強烈. 這也可以解釋為什么1.0 eV時基態(tài)的共振峰的數(shù)目多于其他能量.

        圖 5 置換反應(yīng)總的積分截面和振動分辨的積分截面Fig. 5. The total and vibrational state-resolved integral cross section of exchange reaction.

        圖5 給出置換反應(yīng)的總積分截面和振動分辨的積分截面. 它們的形狀與圖2中的反應(yīng)概率類似,都是隨著碰撞的增加單調(diào)遞增. 并且隨著碰撞的增加, 其他振動激發(fā)態(tài)的通道逐漸打開. 產(chǎn)物的振動態(tài)分布也是隨著振動激發(fā)態(tài)的增加逐漸遞減. 主要的振動態(tài)分布集中在基態(tài)和第一激發(fā)態(tài)上.

        圖6給出了置換反應(yīng)在碰撞為0.8, 1.2,1.6和2.0 eV時的產(chǎn)物分布. 與圖4類似, 隨著碰撞能的升高, 振動態(tài)的數(shù)目和轉(zhuǎn)動量子數(shù)的數(shù)目都增加, 表明更多的碰撞能轉(zhuǎn)化為內(nèi)能.在0.8和2.0 eV時的基態(tài)出現(xiàn)了有趣的雙峰結(jié)構(gòu).1.2和1.6 eV基態(tài)的峰偏向于較大轉(zhuǎn)動量子數(shù)的位置.

        3.3 微分截面

        圖7給出了不同碰撞能下抽取反應(yīng)和置換反應(yīng)的微分截面. 對于抽取反應(yīng), 在低碰撞能時主要有后向和側(cè)向散射信號, 這表明在低碰撞能時, 頭碰頭的反彈機制反應(yīng)機制占據(jù)主導(dǎo)地位. 隨著能量的升高, 后向散射變的不明顯, 前向散射變得越加明顯, 這表明反應(yīng)機制由頭碰頭的反彈機制轉(zhuǎn)變?yōu)閯冸x反應(yīng)機制. 側(cè)向的信號在Li等[25]的分析中是因為勢能面上兩個電子態(tài)交叉的勢阱引起的, 是不可靠的. 對于置換反應(yīng), 在低碰撞能時, 散射信號主要集中在后向散射. 隨著碰撞的增加, 后向散射信號變得明顯, 同時側(cè)向和前向的散射信號也變得明顯. 不難看出對于置換反應(yīng), 反彈反應(yīng)機制一直占據(jù)主導(dǎo)地位. 隨著碰撞能的增加, D的速率變得越來越大, 使其帶走質(zhì)量相對較重的Br原子成為可能, 因此出現(xiàn)了側(cè)向和前向的散射信號.

        3.4 速率常數(shù)

        抽取反應(yīng)和置換反應(yīng)在不同溫度范圍內(nèi)的速率常數(shù)分別示于圖8(a)和圖8(b)中. 為了與相應(yīng)的理論和實驗結(jié)果比較, 文獻[1,6,24]的結(jié)果也列入到了圖8(a)中. 對于抽取反應(yīng), 速率常數(shù)隨著溫度的升高線性增加. 此外, 本文結(jié)果與Zhang等[24]的非絕熱結(jié)果以及實驗值都十分吻合. 盡管在低溫時出現(xiàn)了一些微小的差距. 這可能是因為本文的結(jié)果是基于絕熱的勢能面, 另一方面與實驗結(jié)果比較需要包含轉(zhuǎn)動激發(fā)態(tài)的貢獻, 而本文只考慮了基態(tài)的結(jié)果. 圖8(b)給出了0—3000 K溫度范圍內(nèi)置換反應(yīng)的速率常數(shù). 如圖8(b)所示, 速率常數(shù)隨著溫度升高而增加.

        圖 6 在若干選擇的碰撞能下, 置換反應(yīng)的產(chǎn)物的振轉(zhuǎn)分布Fig. 6. The ro-vibrational distribution of products for exchange reaction at several selected collision energies.

        圖 7 (a) 在若干碰撞能下, 抽取反應(yīng)的微分截面; (b) 在若干碰撞能下, 置換反應(yīng)的微分截面Fig. 7. (a) The differential cross sections of abstraction reaction at several selected collision energies; (b) the differential cross sections of exchange reaction at several selected collision energies.

        圖 8 (a) 在200-1000 K溫度范圍內(nèi), 抽取反應(yīng)的速率常數(shù)以及文獻[1,6,24]的結(jié)果; (b) 在0-3000 K溫度范圍內(nèi), 置換反應(yīng)的速率常數(shù)Fig. 8. (a) The rate constant of abstraction reaction and the values obtained from Refs.[1,6,24] in the temperature range from 200 to 1000 K; (b) the rate constant of the exchange reaction in the temperature range from 0 to 3000 K.

        4 結(jié) 論

        基于Li等的勢能面采用含時量子波包方法和二階分裂算符傳播子對D + DBr反應(yīng)的抽取通道和置換通道進行了動力學(xué)計算. 在量子態(tài)分辨的理論層次上報道了反應(yīng)概率、積分截面、微分截面、產(chǎn)物振轉(zhuǎn)態(tài)分布以及速率常數(shù)等動力學(xué)信息. 積分截面結(jié)果與Zhang等[24]結(jié)果十分吻合, 盡管在高碰撞能部分有一些微小的差距. 速率常數(shù)以及實驗結(jié)果與Zhang等[24]的結(jié)果進行了比較. 結(jié)果表明,本文結(jié)果與實驗結(jié)果十分吻合. 微分截面信息表明, 對于抽取反應(yīng)在低碰撞能時, 反彈機制占據(jù)主導(dǎo)地位. 剝離反應(yīng)機制在高碰撞能時占據(jù)主導(dǎo)地位. 對于置換反應(yīng), 反彈機制一直占據(jù)主導(dǎo)地位,然而隨著碰撞能升高時, 逐漸出現(xiàn)前向和側(cè)向的散射信號.

        猜你喜歡
        勢阱波包勢能
        “動能和勢能”知識鞏固
        含有陡峭勢阱和凹凸非線性項的Kirchhoff型問題的多重正解
        作 品:景觀設(shè)計
        ——《勢能》
        文化縱橫(2022年3期)2022-09-07 11:43:18
        分?jǐn)?shù)階量子力學(xué)下的二維無限深方勢阱
        “動能和勢能”知識鞏固
        時空分?jǐn)?shù)階量子力學(xué)下的δ勢阱
        “動能和勢能”隨堂練
        對稱三勢阱玻色—愛因斯坦凝聚體的非線性效應(yīng)
        基于小波包Tsallis熵和RVM的模擬電路故障診斷
        基于小波包變換的電力系統(tǒng)諧波分析
        一本之道高清无码视频| 一区二区三区在线视频爽| 亚洲精品一区二区三区四区| 亚洲一区二区在线观看网址| 久久久精品午夜免费不卡| 国产太嫩了在线观看| 日韩精品内射视频免费观看| 亚洲av综合日韩| 麻豆国产人妻欲求不满| 亚洲AV永久无码精品一区二国| 最新国产av网址大全| 97成人精品在线视频| 国产精品久色婷婷不卡| 无码专区人妻系列日韩精品 | 亚洲一区二区三区中国| 色一情一乱一伦麻豆| 久久国产精品久久久久久| 国产熟妇搡bbbb搡bb七区| 日韩在线不卡免费视频| 亚洲素人日韩av中文字幕| 一本之道久久一区二区三区| 色一情一乱一伦麻豆| 超碰97人人做人人爱少妇| 中文字幕国产欧美| 中文字幕一区二区va| 亚洲女优中文字幕在线观看| 精品国偷自产在线视频九色| 在线观看国产成人av片| 香蕉亚洲欧洲在线一区| 国产免费一区二区三区三| 一区二区三区无码高清视频| 一区二区三区乱码在线 | 欧洲 | 99视频全部免费精品全部四虎| 国内精品视频成人一区二区| 日本一区二区不卡二区| 包皮上有一点一点白色的| 无码少妇一区二区浪潮av| 国产精品中文第一字幕| 久久综合激情的五月天| 18黑白丝水手服自慰喷水网站| 日日噜噜夜夜狠狠久久无码区|