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

        ?

        基于晶格化的動(dòng)力學(xué)Monte Carlo方法研究鈦酞菁配合物氧化還原反應(yīng)催化能力

        2019-05-08 00:44:04蘇,哲,揚(yáng),
        關(guān)鍵詞:氫氧催化活性比值

        周 蘇, 胡 哲, 高 揚(yáng), 姜 縝

        (同濟(jì)大學(xué) a. 汽車學(xué)院; b. 中德學(xué)院, 上海 201804)

        一般情況下,比較催化劑的性能,可以通過實(shí)驗(yàn)獲得以下兩個(gè)參數(shù):① 單位時(shí)間內(nèi)單位活性中心所轉(zhuǎn)化的反應(yīng)物分子數(shù)(Turnover Frequency,TOF),即單位時(shí)間內(nèi)在催化劑單位面積每個(gè)活性位置的轉(zhuǎn)化數(shù);② 活性位置密度(Site Density,SD),即催化劑單位面積上的活性位置.但是,對于未知或者具有新型結(jié)構(gòu)的催化劑,很難通過實(shí)驗(yàn)得到TOF值或SD值.例如,在研究一個(gè)可能的催化反應(yīng)機(jī)理時(shí),整個(gè)反應(yīng)過程持續(xù)時(shí)間很短,即使超高速攝像機(jī)也只能獲得若干張宏觀變化的圖像.若想了解其微觀層級(jí)的變化并進(jìn)一步推測其催化性能,并與其他催化劑進(jìn)行比較,需借助其他方法或工具.

        隨著計(jì)算機(jī)技術(shù)的發(fā)展,計(jì)算化學(xué)已經(jīng)能夠在沒有實(shí)驗(yàn)情況下給出定性或定量結(jié)論,為尋找或設(shè)計(jì)潛在的新型催化劑材料提供方向性的信息[1].例如,以鈦酞菁配合物(TiPc)為研究對象,周蘇等[2]基于密度泛函理論構(gòu)建了TiPc氧化還原反應(yīng)機(jī)理模型,通過計(jì)算其Gibbs自由能等指標(biāo),分析預(yù)測了TiPc的催化能力.動(dòng)力學(xué)Monte Carlo方法(KMC)也是一種應(yīng)用廣泛的計(jì)算化學(xué)工具[3].例如,湯星舟[4]結(jié)合動(dòng)態(tài)Monte Carlo法與分子動(dòng)力學(xué)研究了C面藍(lán)寶石襯底的氮化現(xiàn)象,分析總結(jié)了表面氮化的相關(guān)機(jī)理;Franziska Hess等[5]利用KMC研究了催化劑RuO2對于CO和HCl氧化催化反應(yīng)的作用.

        本文以TiPc為研究對象,基于晶格化的KMC并借助Kmos軟件,通過求解Markov鏈主方程,計(jì)算TiPc的TOF值并以其作為催化活性的參考評價(jià)指標(biāo),研究反應(yīng)溫度和氫氧氣體分壓比對TiPc催化活性的影響;通過比較未知TiPc和已知Pt兩者的TOF值,預(yù)測未知催化劑TiPc的催化活性;另一方面,從計(jì)算化學(xué)工具的角度,以TiPc作為對象,旨在驗(yàn)證KMC用于研究未知催化劑催化活性的可行性.

        1 模型與方法

        1.1 計(jì)算軟件

        圖1 Kmos軟件GUI界面(氫氧反應(yīng)圖)Fig.1 GUI interface of Kmos software (hydrogen and oxygen reaction)

        Kmos是基于Python的一款開源軟件,其模型建立是由用戶通過Python編寫py文件作為輸入,然后由Kmos模擬給出結(jié)果.為了展示催化劑表面的微觀變化,該軟件也支持UI界面,圖1所示為氫氣和氧氣在催化劑表面某一反應(yīng)時(shí)刻的截圖,圖中深色點(diǎn)代表氧原子,淺色點(diǎn)代表氫原子,白色點(diǎn)代表空位,基底圓形代表催化劑.

        1.2 基于晶格化的KMC及應(yīng)用的假設(shè)條件

        基于晶格化的KMC將整個(gè)催化劑薄膜簡化為10×10個(gè)晶胞所組成的平面,同時(shí)以單個(gè)晶胞作為研究對象,以基本反應(yīng)的動(dòng)力學(xué)參數(shù)作為輸入,在晶胞表面采用KMC求解Markov鏈主方程[5-6,7,9]及微觀動(dòng)力學(xué)模型;根據(jù)給定的反應(yīng)條件(如溫度和所有反應(yīng)物的分壓),可以計(jì)算或獲取催化劑的TOF值、晶胞表面的反應(yīng)物或/和生成物成分、反應(yīng)鏈中的步驟、反應(yīng)機(jī)理等信息[10-11].針對所研究的對象,應(yīng)用晶格化的KMC時(shí),通常給出如下4個(gè)方面的假設(shè).

        1.2.1特殊位置的吸附和晶格排布 假設(shè)所有反應(yīng)物和中間物都吸附在特殊位置,即晶體表面的活性位置.在實(shí)際的KMC模擬時(shí),只考慮這一個(gè)晶格,整個(gè)晶體表面可以看成由這個(gè)晶格重復(fù)排列組合而成.例如,低指數(shù)單晶表面中一個(gè)晶格是由多個(gè)相同的晶胞組成.圖2所示為一個(gè)面心立方排列(100)的催化劑模型,它有2個(gè)活性位置[5].

        圖2 催化劑表面CO氧化反應(yīng)(包含吸附、擴(kuò)散、反應(yīng)和解吸附過程)模型示意圖Fig.2 Catalyst surface diagram of CO oxidation reaction on catalyst surface (including adsorption, diffusion, reaction and desorption process)

        1.2.2Rare-event動(dòng)力學(xué)和Markov鏈主方程 假設(shè)時(shí)間演化是基于Rare-event動(dòng)力學(xué)的.利用時(shí)間尺度分離,根據(jù)微觀動(dòng)力學(xué)理論,假設(shè)事件的發(fā)生是相互獨(dú)立的,可以基于Markov近似來推演反應(yīng)過程.因此,事件(從一種組態(tài)到另一種組態(tài))的發(fā)生可以通過Markov鏈主方程描述:

        (1)

        1.2.3基本反應(yīng)的確定 假設(shè)前文提到的一個(gè)晶格中組態(tài)變化只考慮1~10個(gè)位置.因此,可以根據(jù)局部產(chǎn)物Ea的晶格狀態(tài)和變化后的局部產(chǎn)物Pa的晶格狀態(tài)以及相應(yīng)的伴隨速率Ka來確定一個(gè)基本反應(yīng),即

        (2)

        1.2.4狀態(tài)轉(zhuǎn)移矩陣的大小和結(jié)構(gòu) 式(1)中w可分解成一系列基本反應(yīng):

        (3)

        (4)

        式中:ka為基本反應(yīng)a的反應(yīng)速率.

        以上假設(shè)用于簡化Markov鏈主方程的求解.

        1.3 模型

        1.3.1催化反應(yīng)速率計(jì)算 KMC模擬的主要特點(diǎn)之一就是可以計(jì)算各種催化劑的反應(yīng)速率,包括TOF值.如果基本反應(yīng)只有一個(gè)最終生成物,那么TOF的計(jì)算只考慮這一種生成物.

        對于基本反應(yīng)a,在任意時(shí)刻t,

        TOFa=

        (5)

        1.3.2參數(shù)計(jì)算 基于密度泛函理論(DFT),利用Amsterdan Density Functional(ADF)等軟件,根據(jù)構(gòu)造的TiPc在氧化還原過程中形成的各粒子能量變化,可計(jì)算TiPc的氧吸附能:

        EbO2=EML-O2-(EO2+EML)

        (6)

        式中:EML-O2為氧分子和吸附粒子的總能量;EO2為孤立的氧分子的總能量;EML為吸附粒子的能量.此式也可用于計(jì)算Pt的氧吸附能.

        同理,計(jì)算TiPc的氫吸附能:

        EbH2=EML-H2-(EH2+EML)

        (7)

        式中:EML-H2為氫分子和吸附粒子的總能量;EH2為孤立的氫分子的總能量;EML為吸附粒子的能量.此式也用于計(jì)算Pt的氫吸附能.

        2 結(jié)果與討論

        2.1 TiPc的TOF值

        在模擬TiPc表面氫氣和氧氣反應(yīng)時(shí),反應(yīng)溫度T和TiPc表面的氫氣分壓與氧氣分壓的比值pH2/pO2是兩個(gè)重要的參數(shù).在Kmos軟件運(yùn)行環(huán)境下,取T在470~570 K(間隔5 K)之間的21個(gè)值和pH2/pO2在 0.1~10之間的40個(gè)值,可以計(jì)算得到相應(yīng)的TOF值.以溫度為470、520和570 K的情況為例,分析反應(yīng)溫度T和氫氧氣體分壓比值pH2/pO2對TiPc催化劑TOF值的影響.

        2.1.1TOF值隨反應(yīng)溫度T的變化 由化學(xué)反應(yīng)基本原理可知,反應(yīng)溫度T是影響催化劑活性(TOF值)的重要因素.溫度較低時(shí),催化劑活性較小,對于反應(yīng)促進(jìn)效果較差;隨著溫度上升,催化效果逐漸增強(qiáng).但是,大多數(shù)催化劑都有其催化活性適宜溫度范圍.溫度過高,會(huì)使催化劑發(fā)生燒結(jié),從而破壞其催化活性.為了分析反應(yīng)溫度對TOF值的影響,圖3所示為3種溫度條件下不同氫氧氣體分壓比值對應(yīng)的TOF值(實(shí)線部分).與溫度條件520 K及570 K相比,溫度470 K下的TOF值過低,所以圖3的虛線部分也給出了lg(TOF)值,以便于對比分析.

        圖3 不同溫度下pH2/pO2對應(yīng)的TiPc催化劑TOF及l(fā)g(TOF)值Fig.3 TOF and lg(TOF) values of TiPc under different pH2/pO2 at 470, 520, 570 K

        由圖3可明顯看出:① 隨著反應(yīng)溫度上升,同一氣體分壓比值對應(yīng)的TOF值也依次增高,這符合催化劑基本特性.但反應(yīng)溫度超過催化劑適宜反應(yīng)溫度區(qū)間后,催化劑催化性能下降.Choi等[14]及陳維民等[15]通過各自的研究發(fā)現(xiàn),催化劑會(huì)隨著溫度上升而發(fā)生聚集,超過一定程度會(huì)降低催化劑催化活性;② 與溫度區(qū)間520~570 K相比,在溫度區(qū)間470~520 K內(nèi)TOF值的增長幅度較小,這說明TiPc的活性雖然與溫度正相關(guān),但并不是線性相關(guān);③ lg(TOF)值變化隨溫度變化近似呈線性關(guān)系.

        2.1.2TOF值隨氫氧氣體分壓比值pH2/pO2的變化 除反應(yīng)溫度T外,氫氧氣體分壓比值pH2/pO2同樣會(huì)影響催化劑活性(即TOF值).為了研究相關(guān)規(guī)律,將3種溫度條件下TOF值(圖3實(shí)線部分)在橫坐標(biāo) 0.1~0.33 區(qū)間內(nèi)的曲線放大得到圖4.

        圖4 不同溫度下pH2/pO2在 0.1~0.33 區(qū)間時(shí)對應(yīng)的TiPc催化劑TOF值Fig.4 TOF values of TiPc within pH2/pO2 0.1-0.33 at 470, 520, 570 K

        圖5 溫度470至570 K條件下不同氧氫氣體分壓比值對應(yīng)的TiPc催化劑lg(TOF)值Fig.5 lg(TOF) values of TiPc under different O2 and H2 partial pressure ratios between 470 K and 570 K

        結(jié)合圖3和4的信息可知,3種不同溫度條件下TOF值的最大值出現(xiàn)在氫氧氣體分壓比值0.289~0.326區(qū)間內(nèi).在氫氧氣體分壓比值0.1~0.2之間,TOF值單調(diào)緩慢上升;在氫氧氣體分壓比值0.2~2區(qū)間內(nèi),TOF值先急速上升至最高點(diǎn),然后單調(diào)急速下降;從氫氧氣體分壓比2開始,TOF值單調(diào)平穩(wěn)下降至接近0值.催化反應(yīng)主要經(jīng)歷以下3個(gè)步驟[12-13]:① 反應(yīng)物被催化劑吸附在表面;② 反應(yīng)物在催化劑表面進(jìn)行反應(yīng);③ 反應(yīng)生成物從催化劑表面脫附.其中,步驟①對應(yīng)前文所計(jì)算吸附能,步驟③保證反應(yīng)生成物能夠及時(shí)富集及保證催化劑再生.根據(jù) 1.2 中的假設(shè),所有的反應(yīng)物和中間物都吸附在表層催化活性位置.因此,根據(jù)化學(xué)平衡移動(dòng)原理,當(dāng)氫氧分壓比從零開始升高時(shí),催化劑表面反應(yīng)物氫氣的濃度增加,反應(yīng)傾向于向正反應(yīng)方向移動(dòng),對于催化劑表面的活化中心具有很好的活化作用.但是,氫氧分壓比進(jìn)一步增大,使得氫氣的覆蓋度過大,會(huì)出現(xiàn)“氫中毒”現(xiàn)象,反而使催化劑活性下降,抑制反應(yīng)的進(jìn)行.本文所用方法不能從定量上解釋TOF值急速變化的原因,但是,從上述的分析中可以推測,TiPc催化劑的活性區(qū)域在氫氧分壓比值 0.2 ~2之間,此區(qū)間內(nèi)氫氣和氧氣在催化劑表面以較為合適正向反應(yīng)的比例共存.

        2.1.3采用2D Map圖分析TOF 為了能夠同時(shí)研究反應(yīng)氫氧氣體分壓比及反應(yīng)溫度對TiPc催化活性或反應(yīng)速率的影響,利用線性插值,將上述21組數(shù)據(jù)繪制成2D Map圖;為了分析方便,取氫氧氣體分壓比值的倒數(shù)即氧氫氣體分壓比值為橫坐標(biāo),縱坐標(biāo)為反應(yīng)溫度(見圖5);二維圖中區(qū)域的顏色深淺則對應(yīng)TOF值(反應(yīng)速率),顏色越淺代表TOF值越大.從TiPc的計(jì)算數(shù)據(jù)及上文的結(jié)論可知,在氫氧氣體分壓比值 0.1~0.2 區(qū)間內(nèi),TiPc的TOF值非常小,接近于零.為了增強(qiáng)這種區(qū)別,取TOF值的常用對數(shù).圖5所示為反應(yīng)溫度 470~570 K 下不同氧氫氣體氣體分壓比值對應(yīng)的 lg(TOF) 值.

        圖5中暗部區(qū)域?qū)?yīng)的橫坐標(biāo)范圍為氧氫氣體分壓比值 3.8~10,對應(yīng)于氫氧氣體分壓比值 0.1~0.26,基本符合 2.1.1 和 2.1.2 得出的結(jié)論,該區(qū)域內(nèi)TiPc反應(yīng)活性遠(yuǎn)低于左側(cè)亮部區(qū)域,且亮部區(qū)域和暗部區(qū)域突變位置即為TOF值陡降陡升位置.

        2.2 Pt的TOF值

        采用上述的相同條件,在Kmos軟件運(yùn)行環(huán)境下計(jì)算得到相應(yīng)的Pt催化劑TOF值,用于分析反應(yīng)溫度T及氫氧氣體分壓比值pH2/pO2對TOF值的影響,從而進(jìn)一步比較TiPc與Pt兩者的催化活性.

        2.2.1TOF值隨反應(yīng)溫度T的變化 由Pt催化劑的計(jì)算數(shù)據(jù)可知,Pt催化劑TOF值在氫氧氣體分壓比值 0.8~10 區(qū)間內(nèi)幾乎不變,為了便于分析,只選用氫氧氣體分壓比值 0~0.8 區(qū)間內(nèi)的計(jì)算數(shù)據(jù).圖6所示為3種溫度條件下不同氫氧氣體分壓比值對應(yīng)的TOF值(實(shí)線部分),且與TiPc催化劑情況相同,為突出470 K與520 K溫度下TOF值曲線間的差別,圖6的虛線部分也給出了 lg(TOF)值.

        圖6 不同溫度下pH2/pO2對應(yīng)的Pt催化劑TOF及l(fā)g(TOF)值Fig.6 TOF and lg(TOF) values of Pt under pH2/pO2 at 470, 520, 570 K

        由圖6可看出,不同溫度條件下,Pt催化劑表現(xiàn)出與TiPc催化劑相似的特性:① 隨著反應(yīng)溫度上升,同一氣體分壓比值對應(yīng)的TOF值也依次增高;② 雖然TOF值與溫度不是線性正相關(guān),但是 lg(TOF) 值隨溫度變化近似呈線性關(guān)系.

        2.2.2TOF值隨氫氧氣體分壓比值pH2/pO2的變化 為了研究氫氧氣體分壓比值pH2/pO2對催化劑活性(即TOF值)的影響關(guān)系,將3種溫度條件下TOF值(圖6實(shí)線部分)在橫坐標(biāo) 0.1~0.26 區(qū)間內(nèi)的曲線放大得到圖7.

        由圖7可以看出,3種不同溫度條件下TOF值的最大值均出現(xiàn)在氫氧氣體分壓比值0.229處,此時(shí)Pt催化劑的催化活性最高;與TiPc催化劑相比,Pt催化劑TOF值的最大值明顯小幾個(gè)數(shù)量級(jí),所以也可以預(yù)測TiPc催化劑的活性高于Pt催化劑.

        圖7 不同溫度下pH2/pO2 在 0.1~0.26 區(qū)間時(shí)對應(yīng)的Pt催化劑TOF值Fig.7 TOF values of Pt within pH2/pO2 0.1-0.26 at 470, 520, 570 K

        2.2.3采用2D Map圖分析TOF 將Pt催化劑21組原始數(shù)據(jù)利用線性插值繪制成2D Map圖.由圖6實(shí)線部分可知,在氫氧氣體分壓比值 0.1~0.3 區(qū)間內(nèi),Pt催化劑TOF值遠(yuǎn)大于其他區(qū)域.為了增強(qiáng)這種區(qū)別,取TOF值的常用對數(shù).圖8所示為反應(yīng)溫度470~570 K下不同氧氫氣體分壓比值對應(yīng)的lg(TOF)值.

        圖8中亮部區(qū)域?qū)?yīng)的橫坐標(biāo)范圍為氧氫氣體分壓比值 3.9~10,對應(yīng)氫氧氣體分壓比值 0.1~0.256,基本符合本節(jié) 2.2.1 和 2.2.2 所得結(jié)論,該區(qū)域內(nèi)Pt催化劑反應(yīng)活性遠(yuǎn)高于左側(cè)暗部區(qū)域,且亮部區(qū)域和暗部區(qū)域突變位置即為TOF值陡降陡升位置.

        圖8 470~570 K條件下不同氧氫氣體分壓比值對應(yīng)的Pt催化劑lg(TOF)值Fig.8 lg(TOF) values of Pt under different O2 and H2 partial pressure ratios between 470 K and 570 K

        2.3 分析與討論

        TiPc催化劑TOF值與反應(yīng)溫度T及氫氧氣體分壓比值pH2/pO2存在一定復(fù)雜關(guān)系.TiPc催化劑TOF值會(huì)隨溫度變化而發(fā)生改變,隨著反應(yīng)溫度由470 K增長到570 K,TOF值同樣呈現(xiàn)出遞增的趨勢;氫氧氣體分壓比值同樣會(huì)影響催化劑活性(TOF值),TiPc催化劑TOF值的最大值出現(xiàn)在氫氧氣體分壓比值 0.289~0.326 區(qū)間內(nèi),氫氧氣體分壓比值在 0.1 至該區(qū)間,TOF值首先接近于零,之后單調(diào)遞增;氫氧氣體分壓比值在該區(qū)間至10,TOF值單調(diào)遞減直至接近0值.此外,TiPc催化劑TOF值的最大值出現(xiàn)在反應(yīng)溫度570 K和氫氧氣體分壓比值 0.326 處,其值為 2.034×107.

        與TiPc催化劑不同,Pt催化劑TOF值在氫氧氣體分壓比值 0.1~0.229 之間單調(diào)遞增,在氫氧氣體分壓比值 0.229~0.26 之間單調(diào)遞減,可以推測氫氧氣體分壓比值 0.1~0.26 區(qū)間是Pt催化劑的活性區(qū)域;當(dāng)氫氧氣體分壓比值處在 0.26~10區(qū)間內(nèi),TOF值接近于零且單調(diào)遞減.Pt催化劑TOF值的最大值出現(xiàn)在反應(yīng)溫度570 K和氫氧氣體分壓比值 0.229 處,其值約為82.

        與Pt催化劑相比,TiPc催化劑TOF值的最大值以及TiPc催化劑的2D Map圖中亮部區(qū)域?qū)?yīng)氫氧氣體分壓比值區(qū)間都明顯更大,可以推測TiPc催化劑的活性區(qū)域更廣且其催化活性更高.

        3 結(jié)語

        基于晶格化的動(dòng)力學(xué)Monte Carlo方法,采用TOF值作為參考評價(jià)指標(biāo),評價(jià)TiPc催化劑在反應(yīng)中的催化活性,模擬氣-固界面氫、氧分子在催化劑表面的反應(yīng).討論了TiPc在不同反應(yīng)溫度及不同氫氧氣體分壓比值的反應(yīng)條件下的催化活性.TiPc催化劑TOF值隨反應(yīng)溫度上升而遞增,同時(shí)也隨氫氧氣體分壓比值上升而首先遞增之后遞減.TiPc催化劑TOF值的最大值出現(xiàn)在反應(yīng)溫度為570 K和氫氧氣體分壓比值 0.326 處.采用與TiPc催化劑相同的反應(yīng)條件,在Kmos軟件運(yùn)行環(huán)境下計(jì)算得到相應(yīng)的Pt催化劑TOF值,并與TiPc催化劑TOF值進(jìn)行比較,得到了與文獻(xiàn)[2]相同的預(yù)測結(jié)果,即TiPc催化劑催化活性強(qiáng)于Pt催化劑.

        因計(jì)算方法和計(jì)算能力的局限性,本文僅研究了TiPc在相對較簡單的氣-固界面上反應(yīng)的催化性能,沒有涉及質(zhì)子交換膜燃料電池中催化劑在固-液界面的真實(shí)反應(yīng)環(huán)境.

        猜你喜歡
        氫氧催化活性比值
        氫氧燃料電池演示實(shí)驗(yàn)的改進(jìn)
        自制液壓儲(chǔ)氣式氫氧燃料電池
        連鑄坯氫氧切割應(yīng)用與碳排放研究
        比值遙感蝕變信息提取及閾值確定(插圖)
        河北遙感(2017年2期)2017-08-07 14:49:00
        不同應(yīng)變率比值計(jì)算方法在甲狀腺惡性腫瘤診斷中的應(yīng)用
        稀土La摻雜的Ti/nanoTiO2膜電極的制備及電催化活性
        環(huán)化聚丙烯腈/TiO2納米復(fù)合材料的制備及可見光催化活性
        大推力氫氧火箭發(fā)動(dòng)機(jī)試驗(yàn)噪聲治理技術(shù)概述
        Fe3+摻雜三維分級(jí)納米Bi2WO6的合成及其光催化活性增強(qiáng)機(jī)理
        LaCoO3催化劑的制備及其在甲烷催化燃燒反應(yīng)中的催化活性
        日本人妻少妇精品视频专区| 天堂av在线一区二区| 一级做a爱视频在线播放| 在线一区二区三区视频观看| 亚洲国产精品免费一区| 亚洲av手机在线观看| av网站不卡的av在线| 国产一级一区二区三区在线播放| 亚洲av老熟女一区二区三区 | 亚洲国产精品中文字幕日韩| 国色天香精品亚洲精品| 精品日本免费观看一区二区三区| 国产自拍三级黄片视频| 色综合久久中文综合网亚洲| 久久99精品久久久大学生| 久久婷婷五月综合色欧美| 日本理伦片午夜理伦片| 国产精品天堂avav在线| 久久精品无码一区二区三区不卡 | 亚洲精品白浆高清久久| 用力草我小逼视频在线播放| 国产女主播大秀在线观看| 不卡一区二区视频日本| 婷婷色综合视频在线观看| 亚洲欧美牲交| 国产精品伦一区二区三级视频 | 一二三区无线乱码中文在线| 国产国产精品人在线视| 日韩精品视频一区二区三区| 亚洲av片一区二区三区| 亚洲av纯肉无码精品动漫| 成人国产精品一区二区网站| 丁香六月久久| 久久精品国产亚洲av夜夜| 久久免费看黄a级毛片| 亚洲成a v人片在线观看| 亚洲精品国精品久久99热一| 亚洲综合偷自成人网第页色| www久久久888| av网站在线观看二区| 国产一区二区三区十八区|