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

        ?

        熔融Cu57團簇在急冷過程中弛豫和局域結構轉變的分子動力學研究*

        2010-09-08 06:05:06樊沁娜李蔚張林
        物理學報 2010年4期
        關鍵詞:局域熔融原子

        樊沁娜 李蔚 張林

        (東北大學理學院,沈陽110004)

        (2009年5月11日收到;2009年7月3日收到修改稿)

        熔融Cu57團簇在急冷過程中弛豫和局域結構轉變的分子動力學研究*

        樊沁娜 李蔚 張林?

        (東北大學理學院,沈陽110004)

        (2009年5月11日收到;2009年7月3日收到修改稿)

        采用基于嵌入原子方法的正則系綜分子動力學研究熔融Cu57團簇在急冷過程中的弛豫及其局域結構變化.通過對弛豫過程中均方位移、非相干中間散射函數(shù)和非Gauss參數(shù)三種函數(shù)和原子鍵對隨急冷溫度不同所發(fā)生變化的分析表明,在經過短時間的原子劇烈運動后,急冷溫度極大地影響著團簇內原子結構弛豫過程.急冷溫度較高時,原子在經歷短時間劇烈運動的β弛豫后,進入α弛豫區(qū)后以擴散運動為主,隨后原子運動表現(xiàn)為非擴散性的原子局域結構重排,團簇內沒有出現(xiàn)明顯的成核結構.隨著溫度的降低,原子局域結構的變化在經過短時間原子劇烈運動的β弛豫后,在α弛豫區(qū)原子運動表現(xiàn)為擴散性運動,并出現(xiàn)一定數(shù)量的不穩(wěn)定二十面體結構.當急冷溫度很低時,在進入α弛豫區(qū)后,團簇結構變化逐漸表現(xiàn)為非擴散性原子局域結構重排,形成相當數(shù)量的穩(wěn)定成核二十面體結構.

        團簇,分子動力學,計算機模擬,表面

        PACC:3640B,7115Q,6185,6800

        1. 引言

        貨幣金屬金、銀和銅都具有d10s1價電子結構,近些年來由這三種金屬所組成的團簇,因其具有迥異于塊體材料的豐富多樣的物理與化學性質日益受到研究者們的廣泛關注,并在生物及生物物理上具有廣泛的應用[1—5].在通過物理方法實驗制備金屬團簇的過程中,大多采用氣相或離化的金屬在不同方式下快速冷卻凝結得到,盡管實驗觀測到結構各異的團簇結構,但對于包含幾十個原子的小尺寸團簇在凝結過程中原子是如何運動及其所導致的結構變化目前尚知之甚少.了解這些小尺寸團簇在冷卻凝結過程中的弛豫過程、成核和局域結構的變化,對于預測和控制團簇尺寸進而實驗合成具有特定結構及性質的納米顆粒具有指導意義.

        為了研究金屬熔融液體在過冷時所發(fā)生的結晶轉變,1984年由Bengtzelius[6]等以及Leutheusser[7]在有關液態(tài)理論的基礎上提出了MCT(mode-coupling theory)理論,這種理論認為結晶轉變過程是一個動力學轉變,在發(fā)生凝結的液體中存在兩種弛豫過程,慢的α弛豫和快的β弛豫,α弛豫對應結構弛豫過程,該過程受溫度影響較大,β弛豫是二級弛豫過程,該過程受溫度影響較小.通過研究一些隨時間變化相關函數(shù)的變化,MCT理論提出當冷卻溫度接近某一溫度時,會在β弛豫末期也就是α弛豫前期出現(xiàn)一個過渡區(qū)間.應用該理論,研究者們已經對普通液體、金屬熔融液體和較大的金屬團簇在冷卻過程中的弛豫特性作了較多的探討[8—14].但是對于包含幾十個原子的小尺度自由表面金屬團簇,相關的研究則尚未見到.這些小團簇由于它們具有很高的比表面(即表面積與體積之比),使得它們隨溫度改變所發(fā)生的結構轉變過程與大尺寸團簇和塊體相比,展現(xiàn)出不同的特點.其中作為一種重要貨幣金屬團簇的銅團簇以其在納米催化劑、數(shù)據(jù)存儲介質、醫(yī)藥診斷以及生物傳感器等方面的潛在應用價值,備受研究者們的關注[15—20].

        本文采用基于嵌入原子方法(embedded atom method,EAM)的正則系綜分子動力學(molecular dynamics,MD),通過均方位移(mean squaredisplacement,MSD)、非相干中間散射函數(shù)(incoherent intermediate scattering function)和非Gauss參數(shù)(non-Gaussian parameter,NGP)等三種隨時間變化的函數(shù)研究熔融Cu57團簇在被急冷到不同溫度下的弛豫行為,通過鍵對分析(pair analysis,PA)方法給出了在不同急冷溫度下團簇局域結構的變化.

        2. 模擬方法

        原子間的相互作用勢采用Mei等[21]提出的EAM形式,應用該形式所給出勢的分子動力學計算已用于模擬包含幾十個原子小尺寸銅團簇在溫度變化過程中的結構轉變[15,17,19,20].在EAM形式中,體系的總能量Et為

        其中,ρi是其他原子在原子i處產生的電荷密度,F(xiàn)(ρi)是將原子i嵌入到電子密度為ρi處所具有的能量,rij是原子i和原子j之間的距離,(rij)是原子i和原子j之間的兩體勢.模擬采用正則系綜(NVT)分子動力學方法,計算中時間步長取為1.6×10-15s.計算的幾個函數(shù)如下,

        其中〈·〉表示對于統(tǒng)計時間步的平均值.對分布函數(shù)g(r)給出了在相同原子密度下原子隨機分布時,找到相距為r的一對原子的概率.V是所模擬體系的體積,N是體系中的原子數(shù).〈r2(t)〉是均方位移函數(shù)MSD,其中ri(t)和ri(0)是第i個原子分別在時刻t和0時的位移.Fs(q,t)表示非相干中間散射函數(shù),是von Hove函數(shù)自相關部分的Fourier變換,式中的波向量q選取為液態(tài)對分布函數(shù)的第一峰.α2(t)是非Gauss參數(shù)NGP,〈r2(t)〉是均方位移,〈r4(t)〉是平均四次方位移.非高斯參數(shù)被用于描述過冷液體的動力學非均勻特性,當液體是均勻分布并處于平衡態(tài)時,α2(t)接近于零.這里均方位移函數(shù)MSD,F(xiàn)s(q,t)和α2(t)都是隨時間變化描述團簇弛豫特性的函數(shù).

        研究中采用包含四個整數(shù)的鍵對分析技術研究急冷過程中原子局域結構的變化[22,23].第一個整數(shù)表示一對原子是否為設定截斷距離內的近鄰原子(或成鍵原子對),當它為1時表示成鍵.第二個整數(shù)表示某一成鍵原子對的共有近鄰原子數(shù),第三個整數(shù)表示這些共有近鄰原子間的成鍵數(shù)目.第四個為附加數(shù),表示某一特定的原子鍵對.例如,1551對應兩個成鍵原子共有五個近鄰原子,且這五個近鄰原子形成具有五個鍵的五邊形,該鍵對是二十面體結構的特征鍵對.1421和1422鍵對分別用于表征面心立方和密排六方的局域結構.

        模擬中首先構造一個20a0×20a0×20a0的包含32000個原子的Cu晶體(a0=0.3615nm為Cu的晶格常數(shù)).20a0×20a0×20a0為模擬中分子動力學元胞的體積,在這個Cu晶體塊體中截取出原子數(shù)目為57個原子的孤立自由表面銅團簇,將這個團簇直接升溫到1800K,得到處于熔融狀態(tài)的團簇.圖1顯示了這個團簇在1800K時的對分布函數(shù)曲線.由圖中可見,1800K時Cu57團簇的對分布函數(shù)曲線表現(xiàn)出典型的液態(tài)結構特征,即對分布函數(shù)僅明顯存在兩個已經展寬的峰.在圖1中右上側內插圖顯示了該團簇在1800K時處于熔融無序狀態(tài)的原子結構圖.

        圖1 Cu57團簇在1800K時的對分布函數(shù)內插圖為原子結構圖

        3. 結果分析

        在獲得處于熔融狀態(tài)的Cu57團簇后,將其從1800K急冷到1200,1100,1000,900,800,700,600,500,400,300和200K時,圖2給出了該熔融團簇被急冷到不同溫度時的原子平均能量.由圖中可見,Cu57團簇能量的降低可以分為三個階段:在急冷溫度高于1000K時,能量以較大的斜率近直線下降;當急冷溫度低于500K后,能量以低于高急冷溫度區(qū)間斜率的近直線形式降低;在500到1000K的急冷溫度區(qū)間內,能量的變化呈現(xiàn)為振蕩下降的形式.

        圖2 Cu57熔融團簇急冷到不同溫度的原子平均能量

        圖3(a),(b)和(c)分別顯示了熔融銅團簇在被急冷到1200,600和200K時均方位移(MSD)、非高斯參數(shù)(NGP)和非相干中間散射函數(shù)Fs(q,t)隨時間的變化.由圖3(a)可見,當急冷溫度為1200K時,MSD曲線表現(xiàn)為典型的液態(tài)行為.即在短時間0.1 ps的β弛豫區(qū)內,由于原子的集體劇烈運動,MSD以〈r2〉∝t2的形式迅速增加.隨著保溫時間的增加,原子的運動表現(xiàn)為擴散行為,進入α弛豫區(qū),這時〈r2〉∝t.當時間增加到50 ps后,MSD的增速進一步放緩,這表明此時團簇內結構變化以非擴散性原子局域結構重排為主.這里,擴散性原子運動是指原子掙脫了其周圍原子的束縛,此時該原子與周圍原子構成的局域結構也發(fā)生了改變,非擴散性原子局域重排是指雖然某一原子運動,但它仍處于其周圍原子的束縛中,即該原子與其周圍原子構成的局域結構沒有發(fā)生改變,從而在整體上看來是團簇內這些原子局域結構之間位置的重新排布.隨著急冷溫度的降低,MSD曲線出現(xiàn)了不同于較高溫度時的形式.當急冷溫度為600K時,在短時間0.1 ps內的β弛豫區(qū)內,MSD與時間的關系仍表現(xiàn)為〈r2〉∝t2,隨后增速變緩,在經過這個緩慢變化區(qū)到0.2 ps后進入α弛豫區(qū),MSD繼續(xù)增加并與時間具有〈r2〉∝t的關系.隨著急冷溫度進一步降低到200K,進入MSD增速變緩區(qū)的時間縮短到0.08 ps,并在經過一個曲線緩慢增加的階段后,在50 ps后進入一個平臺區(qū).如圖3(a)200K時的MSD曲線所示,這個平臺的高度要遠低于高溫時的MSD曲線.這說明當急冷溫度很低時,每個原子的熱運動降低,其移動范圍變小,在團簇內原子間出現(xiàn)了較多的局域有序結構.另外根據(jù)Einstein關系,〈r2〉=6Dt,即當時間足夠長時,自擴散系數(shù)D可以由MSD曲線的斜率計算得到,很明顯,MSD出現(xiàn)平臺意味著由于原子局域結構重排使得D趨于零.這樣根據(jù)MSD曲線,我們就能確定在什么時間上過冷熔融團簇內開始了成核的轉變過程.由圖中可見,當急冷溫度為600K時,過冷團簇內沒有發(fā)生穩(wěn)定的成核或者說團簇內局域有序結構還不能夠較長時間的存在,而當急冷溫度降低到200K時,過冷團簇在50 ps后出現(xiàn)了穩(wěn)定的成核,這時存在了較多的局域有序結構.這種轉變的進一步證據(jù)可由非相干中間散射函數(shù)Fs(q,t)和非高斯參數(shù)(NGP)顯示出來.

        圖3 熔融Cu57團簇被急冷到1200,600和200K時的均方位移、非相干中間散射函數(shù)和非Gauss參數(shù)隨時間的變化(a)均方位移MSD,(b)非相干中間散射函數(shù)Fs(q,t),(c)非Gauss參數(shù)α2(t)

        由圖3(b)的非相干中間散射函數(shù)Fs(q,t)可見,F(xiàn)s(q,t)曲線從Fs(q,0)=1開始,在短時間內即開始下降,隨著溫度的不同,該函數(shù)下降的行為也表現(xiàn)出差異.在急冷到較高溫度1200K時,F(xiàn)s(q,t)很快衰減為零,這表明此時團簇內的原子運動仍遵循通常液態(tài)的行為.隨著急冷溫度降低到600K時,團簇中的原子會出現(xiàn)結構重組,這時與MSD曲線的變化相對應,在經過一個快速下降的過程后,F(xiàn)s(q,t)的變化出現(xiàn)一個較緩的下降過程,隨后再次下降趨于零,這時非相干中間散射函數(shù)趨于零的時間明顯比較高急冷溫度時要長.當急冷溫度下降到200K時,F(xiàn)s(q,t)曲線在0.1 ps時即進入曲線以較緩坡度下降區(qū)間,并在50 ps后出現(xiàn)了一個平臺,這說明隨著急冷溫度的降低,團簇原子已不會以液態(tài)原子的形式運動,并發(fā)生結構重排,這時Fs(q,t)不會趨于零而是出現(xiàn)一個平臺.

        類似的情形也出現(xiàn)在隨時間變化的非Gauss參數(shù)變化上.如圖3(c)所示,在急冷溫度為1200K時,NGP沒有表現(xiàn)出明顯的峰.但當急冷溫度降低到600和200K時,NGP曲線的各峰出現(xiàn)了升高,這表明被冷卻系統(tǒng)結構變化的復雜程度增加,并且當急冷溫度為200K時,NGP曲線在溫度下降過程中出現(xiàn)了非零的平臺,平臺的出現(xiàn)說明此時團簇內原子已可以長時間處于局域有序的平衡位置.這里需要指出的是,這三個函數(shù)在急冷到不同溫度時隨時間的變化基本符合,即當急冷溫度較高時,團簇內原子所進行的β弛豫和α弛豫如同通常的簡單液體.隨著急冷溫度的降低,除了通常的β弛豫外,在α弛豫區(qū)內會出現(xiàn)一個短時間的非擴散性原子運動行為,但隨后原子們的運動又以擴散性的原子運動為主.當急冷溫度很低時,團簇內原子在經過集體劇烈運動的β弛豫后,α弛豫區(qū)內逐漸表現(xiàn)為非擴散性的原子局域結構重排.這樣上述三個函數(shù)的變化表明所研究團簇在急冷過程中會出現(xiàn)原子集體運動,即當一個被束縛的原子要移動加快或減慢時,它的近鄰原子也會如此,這必然導致所研究體系局域結構的變化.

        圖4 急冷溫度為1200,600和200K時原子鍵對相對數(shù)目隨時間的變化右側圖為其所對應的隨時間變化的原子結構圖. (a)1200K,(b)600K,(c)200K

        圖4(a),(b)和(c)分別給出了團簇內原子在被急冷到1200,600和200K溫度時的1421,1422和1551三種原子鍵對相對數(shù)目隨時間的演變,這里的原子鍵對相對數(shù)目由某一鍵對數(shù)目除以鍵對總數(shù).圖的右側為隨時間變化的原子結構圖.在模擬中的不同溫度處,還會出現(xiàn)其他表征菱形對稱結構的1201,1301和1311等鍵對以及表征液態(tài)和非晶態(tài)結構的1541和1431鍵對等.圖4之所以只選取1421,1422和1551三種鍵對是考慮到在實驗和模擬中發(fā)現(xiàn)大部分的銅團簇具有面心立方結構或二十面體結構,特別是包含五十幾個原子的小尺寸團簇具有基于二十面體的幾何形狀.由7個原子組成的二十面體基元具有1551鍵對特征,隨著由更多的原子組成稍大的二十面體,會出現(xiàn)較多數(shù)量的1422鍵對,并且當團簇尺寸更大時,所形成的二十面體還會有相當數(shù)量的1421鍵對.由圖4中可見,當急冷溫度為1200K時,所研究的各鍵對相對數(shù)目雖然在經過β弛豫后出現(xiàn)了少量的1421,1422和1551鍵對,但在隨后的α弛豫中僅表現(xiàn)出隨時間的振蕩變化,這說明團簇內原子在這個較高急冷溫度下仍表現(xiàn)為如右側原子結構圖所示的無序狀態(tài).隨著急冷溫度降低到600K,雖然團簇內1421和1422鍵對數(shù)目沒有發(fā)生大的改變,但值得注意的是在進入α弛豫區(qū)后,1551鍵對數(shù)目隨時間出現(xiàn)增加的趨勢,這說明隨著急冷溫度的降低,一定數(shù)量的二十面體局域結構開始出現(xiàn),這種趨勢隨著急冷溫度的進一步降低表現(xiàn)得更為明顯.如圖4(c)所示,當急冷溫度為200K時,在0.1 ps后,團簇內1551鍵對數(shù)目明顯增加,并且這種結構在50 ps后仍有相當數(shù)量的存在,這說明大部分的二十面體局域結構可以較長時間地存在.這種1551鍵對的出現(xiàn)也可以由圖右側100 ps時團簇原子排布出現(xiàn)環(huán)狀形狀得到驗證.這里還需指出的是,在急冷溫度200K為600K時,1421和1422鍵對都保持很低的值,這說明團簇內成核的結構主要為二十面體結構,并且這些成核的二十面體沒有發(fā)生明顯的生長.

        4. 結論

        應用基于嵌入原子方法的正則系綜分子動力學研究熔融Cu57團簇冷卻到不同溫度時隨時間變化函數(shù)及團簇局域結構的變化.模擬結果表明,盡管熔融Cu57團簇在被急冷到三個不同溫度時都出現(xiàn)相似的短時原子集體劇烈運動,但在隨后進入α弛豫后,團簇內原子運動及其局域結構隨急冷溫度不同出現(xiàn)了很大的差異.當急冷溫度較高時,α弛豫區(qū)的原子運動先是以原子擴散運動為主,然后出現(xiàn)非擴散性的原子局域結構重排.隨著急冷溫度的降低,在α弛豫前期出現(xiàn)非擴散性的原子局域結構調整,隨后原子局域結構的變化來自于原子的擴散性運動,這時出現(xiàn)一定數(shù)量的二十面體結構.當急冷溫度很低時,在α弛豫區(qū)逐漸表現(xiàn)為非擴散性原子局域結構重排,相當數(shù)量的成核二十面體結構形成,其中大部分可以保持其局域穩(wěn)定結構,這些成核的二十面體結構沒有發(fā)生明顯的生長.

        [1]Quan H J,Gong X G 2000 Chin.Phys.9 656

        [2]Li T X,Ji Y L,Yu S W,Wang G H 2000 Solid State Commun. 116 547

        [3]Xu S N,Zhang L,Zhang C B,Qi Y 2007 Acta Metall.Sin.43 379(in Chinese)[徐送寧、張林、張彩碚、祁陽2007金屬學報43 379]

        [4]Wu H,Desai S R,Wang L S 1996 Phys.Rev.Lett.77 2436

        [5]Liu H B,Ascencio J A,Alvarez M P,Yacaman M J 2001 Surf. Sci.491 88

        [7]Leutheusser E 1984 Phys.Rev.A 29 2765

        [8]Kob W,Andersen H C 1995 Phys.Rev.B 22 4

        [9]Kob W,Andersen H C 1995 Phys.Rev.E 51 5

        [10]Tokuyama M 2006 Physica A 23 62

        [11]Pang H,Jin Z H,Lu K 2003 Phys.Rev.B 67 094113

        [12]Kob W,Donati C,Plimpton S J,Poole P H,Glotzer S 1997 Phys.Rev.Lett.79 15

        [13]Donati C,Glotzer S C,Poole P H,Kob W,Plimpton S J 1999 Phys.Rev.E 60 3

        [14]Gleim T,Kob W,Binder K 1998 Phys.Rev.Lett.81 20

        [15]Zhang L,Zhang C B,Qi Y 2007 Chin.Phys.16 77

        [16]Yang Q W,Zhu R Z,Wen Y H 2005 Acta Phys.Sin.54 89(in Chinese)[楊全文、朱如曾、文玉華2005物理學報54 89]

        [17]Zhang L,Zhang C B,Qi Y 2008 Phys.Lett.A 372 2874

        [18]Zhang L,Zhang C B,Qi Y 2009 Physica B 404 205

        [19]Zhang L,Sun H X 2009 Solid State Commun.149 1722

        [20]Zhang L,Sun H X 2009 Chin.J.Chem.Phys.22 69

        [21]Mei J,Davenport J W,F(xiàn)ernado G W 1991 Phys.Rev.B 43 4653

        [22]Honeycutt J D,Andersen H C 1987 J.Phys.Chem.91 4950

        [23]Clarke A S,Jonsson H 1993 Phys.Rev.E 47 3975

        PACC:3640B,7115Q,6185,6800

        *Project supported by the State Key Development Program for Basic Research of China(Grant No.G2006CB605103).

        ?Corresponding author.E-mail:zhanglin@imp.neu.edu.cn

        Molecular dynamics study of relaxation and local structure changes in a rapidly quenched molten Cu57cluster*

        Fan Qin-Na Li Wei Zhang Lin?
        (College of Science,Northeastern University,Shenyang110004,China)
        (Received 11 May 2009;revised manuscript received 3 July 2009)

        Relaxation and local structure changes of a molten Cu57cluster during rapidly quenching have been studied by molecular dynamics simulation using embedded atom method.With decreasing quenching temperature,atom motion details are analyzed using three factors,including the mean square displacement,incoherent intermediate scattering function,and non-Gaussian parameter,while the local structure changes are identified by pair analysis.Simulation results reveal that after a drastic collective motion of atoms,the temperature greatly affects the relaxation processes of the cooled cluster.At a high quenching temperature,after atoms dramatically move in a β relaxation region,diffusion motion of the atoms plays a dominant roles followed by non-diffusion rearrangements of local atomic structures,and no nucleation occurs.When the temperature decreases,local structure changes of atoms occur as the initial dramatic motion,then through the diffusion of atoms in the α relaxation region,and some unstable icosahedral structures are observed.At a low quenching temperature,the structure changes in the α relaxation region result mainly from non-diffusion rearrangement of the atom positions,and a notable amount of icosahedral structures are formed.

        cluster,molecular dynamics,computer simulation,surface

        book=102,ebook=102

        *國家重點基礎研究發(fā)展計劃(批準號:G2006CB605103)資助的課題.

        ?通訊聯(lián)系人.E-mail:zhanglin@imp.neu.edu.cn

        猜你喜歡
        局域熔融原子
        原子究竟有多小?
        原子可以結合嗎?
        帶你認識原子
        局域積分散列最近鄰查找算法
        電子測試(2018年18期)2018-11-14 02:30:34
        PET成像的高分辨率快速局域重建算法的建立
        sPS/PBA-aPS共混物的結晶與熔融行為
        中國塑料(2015年7期)2015-10-14 01:02:40
        基于局域波法和LSSVM的短期負荷預測
        電測與儀表(2015年7期)2015-04-09 11:39:50
        基于非正交變換的局域波束空時自適應處理
        FINEX熔融還原煉鐵技術簡介
        新疆鋼鐵(2015年3期)2015-02-20 14:13:56
        惰性氣體熔融–紅外光譜法測定硅中的氧
        美女裸体无遮挡黄污网站| 亚洲综合欧美色五月俺也去| 男人的天堂av网站| 亚洲av日韩专区在线观看| 欧美成人在线A免费观看| 素人系列免费在线观看| 蜜桃久久综合一区二区| 永久免费毛片在线播放| 亚洲日本va中文字幕| 九九精品无码专区免费| 五月激情狠狠开心五月| 精品成人av人一区二区三区| 欧美奶涨边摸边做爰视频| 欧美人和黑人牲交网站上线| 人妻无码中文专区久久综合| 亚洲成av在线免费不卡| 极品尤物在线精品一区二区三区| 亚洲av无码一区二区一二区| 又粗又硬又黄又爽的免费视频| 国产视频最新| 女同另类一区二区三区| 色偷偷久久久精品亚洲| 成人网站免费看黄a站视频| 国产精品青草视频免费播放| 日韩人妻av不卡一区二区三区| 一区二区在线视频免费蜜桃| 天天躁夜夜躁狠狠躁2021a2| 国产精品久久婷婷六月丁香| 在线无码免费看黄网站| 中文字幕亚洲精品专区| 亚洲av午夜福利精品一区| 色翁荡息又大又硬又粗又视频图片 | 一本色道久久综合中文字幕| 熟女免费视频一区二区| √天堂中文官网在线| 日韩a毛片免费观看| 国产精品狼人久久久影院| 国产成人国产三级国产精品 | 黑色丝袜秘书夹住巨龙摩擦| 国产丝袜一区二区三区在线不卡| 国产三级三级精品久久|