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

快凝Ni70Ag30 納米顆粒中C15 單元的識別*

2022-09-14 10:09:08江明香田澤安2謝泉高廷紅梁永超陳茜
物理學報 2022年17期
關鍵詞:可視化結構分析

江明香 田澤安2)? 謝泉 高廷紅 梁永超 陳茜

1) (貴州大學大數據與信息工程學院,先進光電材料與技術研究所,貴陽 550025)

2) (湖南大學信息科學與工程學院,長沙 410082)

模擬計算已經成為材料科學的重要手段,從模擬計算輸出的原子坐標得到體系的結構特征是研究材料結構與性能相關性的前提.對于晶胞只含2—6 個原子的簡單(BCC,HCP 和FCC)晶體,數值分析方法只需要確定每個原子的局域特征,拓撲結構相同的原子相互連接即構成晶體區域.但要確定含有幾十上百個原子的晶胞,數值方法的計算量極大.數值分析與可視化相結合是解決此類問題的方法之一.本文采用分子動力學方法快凝得到Ni70Ag30 納米顆粒,發現納米顆粒含有FCC 晶體和大量結構復雜的拓撲密堆(TCP)結構.利用基于最大標準團簇的分析軟件提供的多種可視分析功能,結合晶體學相關知識,采用拓撲構型分析思路,確定了納米顆粒中的TCP 原子構成C15 相.本文使用的分析思路為將來開發復雜晶體結構數值識別軟件提供了算法邏輯.

1 引言

材料的性能由其微觀結構決定,結構復雜的晶體或非晶可能存在特殊物理和化學性質.研究發現金屬化合物中析出的拓撲密堆(topologically closepacked,TCP)相[1,2],對機械性能有很大影響.非晶合金因具有高強度、高斷裂韌性、高阻率和耐腐蝕性等性能也被各大領域廣泛應用.實驗上利用高分辨電子顯微鏡、X 射線單晶衍射等方法結合晶體的對稱特征,確定了A15 相[3-6]、H相[7,8]、μ相[9,10]等1500 多種復雜晶體相的結構單元.然而,由于實驗技術在探測材料內部結構細節上的限制,導致仍有大量實驗樣品的內部結構尚未確定;同樣,非晶結構因缺乏平移和旋轉對稱性,用傳統的衍射方法獲取結構信息時會出現結構信息相互重疊,內部結構細節難以獲取.計算模擬能夠突破實驗技術和成本的限制,是確定原子體系微觀結構重要且有效的手段,但對于結構復雜的材料,其特征結構(最小重復)單元涉及數十乃至上百個原子,數值分析和三維可視化相結合是重要的技術手段.

在原子體系微觀結構的數值分析中,確定每個原子的局域結構特征是分析中程和長程結構的基礎.常用的局域結構分析方法有中心對稱參數分析(centrosymmetry parameter,CSP)[11]、HA-Pair分析[12]、共有近鄰分析(common neighbor analysis,CNA)[13]和Voronoi 方法[14-16]等,這些方法對簡單和高度對稱的晶體結構十分有效,但難以分析非晶和復雜晶體結構,如拓撲密堆(topologically closepacked,TCP)晶體.謝卓成[17]最近開發了一種表征TCP 晶體結構的方法,結合CNA 和CSP 來表征Laves 晶體中原子中程序特征.然而,此方法不能識別復雜體系中的結構單元,且必須預先設定模板和小心設置截斷距離才能實施結構分析,導致分析結果存在參數依賴性.

最大標準團簇分析(the largest standard cluster analysis,LaSCA)[18,19]根據每個原子所處的局部特征自動確定截斷半徑,能準確識別其局部結構—最大標準團簇(the largest standard cluster,LaSC),并根據LaSC 類型對原子進行分類,最近還嚴格定義了TCP 團簇[20].TCP 團簇體現了金屬玻璃的基本特征[21],能完美解釋PDF 曲線第一峰和第二峰分裂的結構起源[22],能成功描述TCP 相中sigma 相[20]和A15 相[23]的結晶度,固態非晶體系內的TCP 原子的含量與其非晶形成能力正相關[21]與玻璃化轉變的結束溫度Tg隨壓強的變化高度一致[24].簡言之,LaSCA 能有效分析復雜體系中的短程序和中程序,結合可視化軟件可分析長程序結構.

本文采用分子動力學模擬Ni70Ag30納米液滴的快凝過程,用LaSCA 初步分析后發現,得到的納米顆粒中含有FCC 晶體和大量排列有序,但分布復雜的TCP 原子.利用可視化軟件的分析功能深入分析了納米顆粒中的TCP 原子的空間分布特征,確定了TCP 原子構成C15 晶體.本文使用的方法可用于分析其他的復雜晶體相,并為將來開發復雜晶體結構數值識別算法提供思路.

2 模擬條件與分析方法

2.1 分子動力學模擬

采用開源的大規模原子/分子并行模擬軟件(large-scale atomic/molecular massively parallel simulator,lammps)[25],使用周期性邊界條件模擬Ni70Ag30納米液滴的快速凝固過程.首先將5000個原子(1500 個Ag 原子和3500 個Ni 原子)放入設定好的立方盒子中,并在顆粒表面增加30 ?的真空層,確保納米顆粒的有限邊界特征同時不丟失原子,原子間的相互作用采用EAM 勢來描述[26],時間步長設為1.0 fs.其次將系統在初始溫度2200 K下等溫弛豫1 ns 達到平衡態.并在NVT 系綜下,以1010K/s 將系統從初始溫度2200 K 到最終溫度200 K 連續降溫,每隔1 K 記錄每個原子的空間位置和系統的其他物理量.

2.2 最大標準團簇分析法(LaSCA)

最大標準團簇分析(LaSCA)克服了分析結果對Rc的依賴性問題[18].當兩個原子的距離小于臨界值Rc時,互為近鄰.由一個原子(中心原子)及其近鄰原子組成的局部結構中(如圖1(a)),每個近鄰原子都可以與中心原子組成一個參考對(RP),如圖1(b)和圖1(c)中帶圈的原子對;然后由RP和它們的共有近鄰(common-near neighbor,CNN)組成中心原子近鄰子團簇(center neighbor subcluster,CNS),如圖1(b)和圖1(c)中的結構.CNS的拓撲結構(如圖1(d)和圖1(e))可采用CNS 指數法—Sijk來描述,其中“S”是Sub-cluster 的首字母,i為CNN 原子總數,j為CNN 成鍵數,k為所有j個鍵中最長鏈的長度.因此,圖1(b)和圖1(c)中的CNS 可標記為S555 和S666.

圖1 LaSC 結構特征示例 (a)以一個中心原子(4852)和它的16 個近鄰原子組成的Z16[12-555,4-666] LaSC;(b)S555由4852 和2367 組成的根原子對和5 個CNN (d)組成;(c)S666 由4852 和322 組成的根原子對及其6 個CNN (e)組成Fig.1.A sample for LaSC structural features: (a) a Z16[12-555,4-666] LaSC composed of a central atom (labelled 4852) and 16 neighbors;(b) a S555 composed of a bonded reference pair (labelled 4852 and 2367) and 5 CNNs (d),a CNS of S666 (c) and its topology represented by 6 CNNs (e).

顯然,隨著Rc的增加,每個原子的近鄰數目會增多,中心原子與近鄰原子之間的成鍵數目增加并變得復雜.當Rc超出某一臨界值時,會出現多鍵點(multi-bonded point,MBP)和公有近鄰子環(common neighbor sub-ring,CNSR),破壞CNS指數與其對應空間結構的一一對應關系.既沒有MBP 也沒有CNSR 的CNS 稱為標準CNS,如果一個團簇中的所有CNS 都是標準CNS,就稱此團簇為標準團簇(standard cluster,SC);以一個原子為中心可以找到多個SC,但足夠大的Rc將破壞SC的條件,因此每個原子周圍的最大標準團簇(LaSC)是唯一的,可以通過某種算法唯一地識別表征[18].

一個標準團簇可以用由CNS 的種類和數量構成的二元組的集合表示,如圖1(a)以4852 為中心原子和它的16 個CNN 組成的LaSC,由12 個S555和4 個S666 組成,記為[12-555,4-666].類似地,二十面體(ICO)、面心立方(FCC)、體心立方(BCC)則分別用[12-555],[12-421],[6-444,8-666]來表征.LaSCA 可標定以每個原子為中心的LaSC 類型,因此體系內的原子可根據LaSC 種類進行分類,如ICO LaSC 和FCC LaSC 中心原子分別稱之為ICO 原子和FCC 原子.

至少含有一個S555,且僅由S444(×n4),S555(×n5)和S666(×n6)三類CNS 組成的LaSC 稱為拓撲密堆團簇(TCP LaSC)[20].TCP LaSC 都是致密的局域結構,其內部只有四面體間隙,廣泛存在于TCP 晶體[27,28]和金屬玻璃中.由歐拉公式(V+F—E=2)很容易證明2×n4+n5=12.由于n5> 0,所以n4只能取6 個值(n4< 6);因此,每個TCP LaSC 都可簡單、清晰、唯一地由一個字母和一個整數組成的TCP 指數“Ln”表示.其中“L”用來表示S444 的個數,當n4=1—5 時,“L”分別為A—E,n4為0 時字母“L”為Z;“n”表示團簇中心原子數的配位.根據定義,圖1(a)中的[12-555,4-666]是一種TCP LaSC,可記為Z16;ICO 結構[12-555]為Z12.

以LaSC 為核心,本團隊自主開發了可視化分析軟件,本文后續分析中的數據和圖片均由該軟件獲得.

3 結果和討論

3.1 相變與結構

平均原子勢能隨溫度的變化(E-T)曲線能簡單直觀地反映相變的初步特征[29].如圖2(a)所示,Ni70Ag30納米顆粒快速凝固過程的E-T曲線上有兩個明顯的突變點,表明發生了一級相變.圖2(b)顯示,FCC,Z16 和Z12 原子數量百分比在相變期間急劇增加,200 K 時總量大于65%.

圖2 Ni70Ag30 納米液滴凝固過程中(a)平均原子勢能和(b)原子數量百分比隨溫度的變化,以及200 K 時(c)納米顆粒的結構分解,(d)FCC 晶區和(e)TCP 晶區.對(d)和(e)中的晶區按原子元素分類得到(f)和(g)Fig.2.(a)The average atomic potential energy and (b)atomic number percentage as a function of temperature for Ni70Ag30 nanoparticles.(c) The structural decomposition of nanoparticles,(d) FCC region and (e) TCP region at 200 K.The element-based color encoded versions (f) and (g) to the FCC and TCP regions,respectively.

圖2(c)中200 K 時的納米顆粒可以分為主要由Ni 原子構成的FCC 晶體區(圖2(d))和主要由Z12 和Z16 原子組成的TCP 晶體結構,如圖2(e).FCC 晶區具有五重對稱性,HCP 原子構成孿生共格界面,但是TCP 晶區的晶體類型需要深入分析才能確定.另外,從圖2(f)可以看到,FCC 晶體完全由Ni 原子構成,不含Ag 原子,而圖2(g)顯示TCP 晶體由Ni 和Ag 兩種原子構成.

3.2 TCP 晶體相結構分析

僅僅使用數值計算來確定復雜的拓撲密堆晶體的結構單元極其困難,借助具有分析功能的可視化軟件可大大降低工作的難度.本文使用基于最大標準團簇分析法的可視化軟件結合直覺思維,確定TCP 晶體的結構單元.首先研究此TCP 相主要構成元素(Z12 和Z16 原子)的空間分布特征.設置“同類直連”搜索方式,然后以位于TCP 區域中間位置編號為3659 的Z12 原子和編號為4852的Z16 原子為中心,連續兩次搜索“當前近鄰”;并將“原子著色”設置為“LaSC 類型”,將Z12 和Z16 原子分別設置為紅色和青色,得到圖3 所示的結果.

圖3 Z16 和Z12 原子的空間分布特征.Z16 原子(4852)為中心,包含(a)第一近鄰和(b)第二近鄰的結構;Z12 原子(3659)為中心,包含(c)第一近鄰和(d)第二近鄰的結構.圖(b)中兩個黃色三角形相互平行但取向相反Fig.3.The spatial distribution characteristics of Z16 and Z12 atoms.The structures composed of Z16 atoms,including (a) the nearest neighbors and (b) the second neighbors.The structures composed of Z12 atoms,including (c) the nearest neighbors and (d) the second neighbor.The two yellow triangles in panel (b) are parallel but opposite.

圖3(a)顯示,每個Z16 原子有4 個Z16 近鄰,構成典型的空間四面體結構,圖3(b)進一步明確Z16 原子占據閃鋅礦晶體的點陣位置[30];但是圖3(c)和圖3(d) 顯示Z12 原子在空間的分布沒有構成簡單晶體結構.因此Z16 原子的分布具有空間平移對稱性,可構成TCP 晶體骨架.

3.3 確定晶胞骨架

繼續圖3 的方式,只能得到具有中心的殼層性結構,這樣的結構單元難以直觀體現晶體材料的平移對稱性和旋轉對稱性.參照閃鋅礦晶胞的結構(如圖4(a)),通過可視化分析來確定這個TCP 晶體相中Z16 原子構成的骨架結構.

圖4(a)顯示,閃鋅礦晶體的晶胞包含的18 個原子可以劃分為5 層,結合數值分析可視化軟件的逐層搜索功能,從位于圖4(a)第三層的原子出發,在納米顆粒體系內逐步甄選符合空間分布特征的其他原子,得到TCP 晶體的晶胞.

在圖2(c)右下角的TCP 區域的中心位置找到一個原子(編號為933),把它定位于圖4(a)中第三層最左側的原子.使用“同類直連”近鄰搜索模式,查找第一近鄰得到圖4(b).參照圖4(a),去掉圖4(b)左側兩個原子2112 和4501;查找剩下的3 個原子的近鄰,可得到圖4(c),圖中已經出現了5 層原子.

在圖4(c)中,{933,4863,3852,4603,2753}處于晶胞的邊界位置(參照圖4(a),下同),且都是4882 和3403 的近鄰原子,可先去掉,查找剩余(圖4(c)虛線框中)4 個原子的近鄰得到圖4(d).

去掉不屬于本結構單元的4 個原子得到圖4(e),繼續去掉本結構單元最外端的5 個原子和{1662,2381}.查找剩下(圖4(e)虛線框中)4 個原子的近鄰即可得到圖4(f).圖中從上往下5 層各原子層依次是{4863,3852,2683},{4882,1523},{933,1662,2381,2141},{3403,3381}和{4603,2753,2423}.對照圖4(a)和圖4(f)可知,長方體8 個頂點中的4 個還未確定,且這4 個頂點與晶胞內其他原子不成鍵,因此不能通過簡單的近鄰查找來確定.

對于理想的閃鋅礦,在圖4(a)中每層原子都處于同一平面內,且五層原子所在的平面相互平行;在圖4(f)中所示的情況下只有第三層的4 個原子中任意三個不在同一直線上,可以用來確定平面的方向指數.然后按以下步驟,確定最后4 個頂點原子.

圖4 通過可視化軟件的限定條件搜索功能,確定由Z16 原子構成的TCP 晶體骨架結構的過程 (a)閃鋅礦晶體結構及其分層特征;(b)原子933 及其第一近 鄰;(c)原子4882 和3403 及其第一近 鄰;(d)原子4882,3403,1662 和2381 及其第 一近鄰;(e)從(d)中去掉不屬于本晶胞的原子后的結構;(f)原子4882,1523,3403 和3381 及其第一近鄰Fig.4.A step-by-step three dimensional (3D) visualization for the identification of a periodic structural unit composed of Z16 atoms through search function with specific conditions (only Z16 atoms are considered) provided by the software: (a) Crystal structure of sphalerite and its layering characteristics with a Z16 atom labelled 933 being leftmost in the third layer;(b) an atom labelled 933 and its nearest neighbors;(c) two atoms {4882,3403} and their nearest neighbors;(d) four atoms {4882,3403,1662,2381} and their nearest neighbors;(e) the structure after removing atoms (not in the unit) from (d);(f) five atoms {4882,1523,3403,3381}and their nearest neighbors.

利用軟件的“三點截面”功能,以原子{933,1662,2381}得到圖5(a)的截面,此時軟件自動保存當前截面的方向矢量.然后使用軟件的“指數截面”功能,指定截面的“通過顆粒”為第一層的中心原子3852,“截取厚度W”取0.5 (略大于一層原子厚度),可得到以3852 為中心的第一層平面內的所有原子如圖5(b)所示.如圖5 中虛線框所示,{4863,3852,2683}在一條直線上,與之垂直的直線上距離3852 最近的兩個原子是2813 和2711,則第一層原子為{4863,3852,2683,2813,2711}.同理,將平面通過原子2753,可以得到如圖5(c)的結果,得到第五層原子為{4603,2753,2423,322,1132}.至此,得到了TCP 相中以Z16 原子為骨架的晶體原胞結構,如圖6(a)所示.

圖5 尋找閃鋅礦結構單元中的4 個頂點原子(與內部原子不成鍵) (a) 第三層原子“933,1662,2381”所在平面;(b) 通過原子3852 的截面;(c) 通過原子2753 的截面.(b)和(c)中虛線框涉及的原子為圖4(a)中的第一和第五層Fig.5.Search for four vertex atoms in the sphalerite structural unit (not bonded to internal atoms): (a) The plane of the third layer atoms (933,1662,2381);(b),(c) the sections of 4852 and 2753 atoms respectively.The atoms involved in the dashed boxes in panels (b) and (c) are the first and fifth layers in Fig.4(a).

圖6 根據晶胞骨架兩兩正交的3 個均分面交集獲得晶胞完整結構的過程(a)晶胞骨架及3 個均分面.(b)在晶胞骨架基礎上,一次“近鄰搜索”得到的結構,截取厚度d∈(0.5d0,0.6d0),d0 為晶格常數;在圖(b)上依次取(a)中紅色虛線框截面得到(c),取橙色虛線框截面得到(d),取藍色虛線框得到(e).(f)基于原子元素分類得到的晶胞.(a),(e),(f)中的黑色實線為人工添加,以便觀察,不表示成鍵關系Fig.6.The process of obtaining the complete structure of the cell is based on the intersection of three equipartition planes of the biorthogonal cell skeleton: (a) Cell skeleton and three equipartition planes.(b) The structure obtained by one "nearest neighbor search" based on the cytoskeleton;section thickness d∈(0.5d0,0.6d0),where d0 is the lattice constant.In panel (b),atoms on the red,orange,and blue dashed boxes in panel (a) are taken for the section to obtain (c),(d) and (e).(f) Crystal cell based on classification of atomic elements.The solid black lines in panels (a),(e),and (f) are added manually for observation.

3.4 晶胞內的其他原子

基于圖6(a)的結構在“無類型限制”前提下進行一次“近鄰搜索”,可以得到圖6(b)所示的結構.這個結構除了包含晶胞內的所有原子外,還包含了一些屬于鄰近晶胞的原子.需要去掉這些多余的,得到單一晶胞的結構.

理想的閃鋅礦晶胞中有3 個兩兩垂直的均分面,分別是{3852,1662,2753,2381},{933,1662,2141,2381}和{3852,933,2753,2141},如圖6(a)中紅、橙、藍虛線框所示.晶胞內的原子到這3 個平面的距離均小于等于晶格常數的一半,而晶胞外的原子則大于晶格常數的一半.因此可通過求這3 個均分面的近鄰原子的交集得到單個晶胞內的所有原子,具體操作如下.

保持“截取厚度”為4.0(略大于晶格常數的一半)和“只考慮當前顯示粒子”,在圖6(b)的基礎上,用{3852,1662,2381}“三點截面”,可得到圖6(c).在圖6(c)基礎上,使用{933,1662,2141}“三點截面”,得到圖6(d).同樣,在圖6(d)基礎上繼續用{3852,933,2141}“三點截面”,得到圖6(e).

圖6(e)所示即為TCP 晶體的晶胞,由圖6(f)可知,Ag 原子占據Z16 位置,而Ni 原子占據Z12原子位置.其中Z16 原子構成閃鋅礦晶胞骨架,其間填充了4 個Z12 原子構成的四面體(稱為Z12四面體),而且這4 個Z12 四面體各出一個頂點構成一個四面體,此四面體的中心與整個晶胞的中心重合.整個晶胞內含有24 個原子: 頂點原子1 個(8×(1/8)=1),面心原子3 個(6×(1/2)=3),內部原子20 個(34—8—6=20).且Z12 與Z16 原子的比例為2∶1,根據這些特征可確定模擬得到的TCP 相屬于Laves 相[31,32]中的C15 相[33],且各原子的Wyckoff 占位如表1 所列[34,35].

表1 C15 相中各原子的Wyckoff 占位Table 1. Wyckoff positions of atoms in the C15 phase.

要完成本文所述的復雜晶體原胞結構的識別,軟件必須提供以下可視化分析功能: 1)標定每個原子的局域結構特征;2)在可視化中顯示每個原子的編號;3)可由用戶指定需要顯示的原子;4)在限定條件下查找近鄰原子;5)支持多種形式的截面分析功能;6)支持分析過程中的邏輯與操作.本文所述的識別過程對軟件使用者的晶體學專業水平要求高,一般從業者難以掌握,因此需要探索更容易使用的復雜晶體結構識別方法和開發相應的可視化分析功能.一種思路是發揮數值分析方法的優勢(擅長描述局域拓撲特征),人工分析每一種復雜晶體的拓撲特征,然后構造對應的拓撲規則,由程序根據拓撲規則來識別.另一種思路是使用AI 技術,構造所有已知的復雜晶體的標準結構,并進行數據增強處理,得到帶標簽的含有非標準結構的數據集,然后使用合適的神經網絡模型進行有監督學習,訓練好的模型就可以用來識別數值模擬系統中的復雜晶體.

4 總結

利用分子動力學方法,對Ni70Ag30納米顆粒的快速凝固過程進行模擬及分析.結果發現,快凝得到的納米顆粒包含不同的晶體區域,利用基于LaSCA 的三維可視化軟件結合專業知識確定了晶體的類型,得出以下結論.

1) Ni70Ag30納米液滴的凝固經歷了結晶相變,得到含有兩個晶相的納米顆粒—Z12 和Z16 組成的復雜晶體相,以及五重孿生FCC 晶相.

2) 利用基于LaSCA 的可視化軟件,確定復雜晶體相中Z16 原子占據閃鋅礦晶體的點陣位置,構成晶體骨架,而Z12 原子為間隙原子.這兩種原子構成C15 晶體相.

3) C15 晶體相中Ag 原子占據Z16 位置,而Ni原子占據Z12 原子位置,且原子數比例為Ni∶Ag=2∶1;而FCC 晶體僅由Ni 原子組成,不含Ag 原子.

猜你喜歡
可視化結構分析
基于CiteSpace的足三里穴研究可視化分析
基于Power BI的油田注水運行動態分析與可視化展示
云南化工(2021年8期)2021-12-21 06:37:54
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
隱蔽失效適航要求符合性驗證分析
基于CGAL和OpenGL的海底地形三維可視化
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
“融評”:黨媒評論的可視化創新
傳媒評論(2019年4期)2019-07-13 05:49:14
電力系統不平衡分析
電子制作(2018年18期)2018-11-14 01:48:24
電力系統及其自動化發展趨勢分析
論《日出》的結構
主站蜘蛛池模板: 亚瑟天堂久久一区二区影院| 91毛片网| 91国语视频| 亚洲无码91视频| 亚洲男人天堂2018| 少妇被粗大的猛烈进出免费视频| 免费三A级毛片视频| 国产97视频在线| 亚洲黄色激情网站| 亚洲国产精品无码AV| 福利在线免费视频| 国产黄色爱视频| 大陆精大陆国产国语精品1024| 老色鬼欧美精品| www.精品国产| 香蕉国产精品视频| 国产成人免费视频精品一区二区| 国产白浆在线观看| 国产美女在线免费观看| 欧美亚洲第一页| 全裸无码专区| 国产综合另类小说色区色噜噜 | 亚洲人成网址| 亚洲aaa视频| 在线免费亚洲无码视频| 国产精品污视频| 国产手机在线小视频免费观看| 欧美亚洲香蕉| 久久久久免费精品国产| 亚洲综合经典在线一区二区| 久久久久久久久18禁秘| 97精品伊人久久大香线蕉| 另类重口100页在线播放| 97视频在线观看免费视频| 日韩高清成人| 国产一区二区三区在线观看视频 | 日韩一区二区三免费高清| 欧美综合在线观看| 99久久这里只精品麻豆| vvvv98国产成人综合青青| 亚洲日本中文综合在线| yy6080理论大片一级久久| 成人福利一区二区视频在线| 国产精品第三页在线看| 999国产精品| 国产精品一区二区不卡的视频| 国产欧美日韩另类| 国产av无码日韩av无码网站| 亚洲精品波多野结衣| 中国特黄美女一级视频| 高清大学生毛片一级| 99国产精品国产高清一区二区| 欧美色99| 韩国自拍偷自拍亚洲精品| 少妇精品网站| 福利一区三区| 日韩在线永久免费播放| 亚洲an第二区国产精品| 亚洲天堂精品在线观看| 亚洲日韩高清在线亚洲专区| 欧美视频在线观看第一页| 久久综合激情网| 国产拍在线| 成人午夜久久| 日韩国产一区二区三区无码| 色网站在线视频| 第一页亚洲| 日日碰狠狠添天天爽| 久久性妇女精品免费| 九九视频免费看| 欧美亚洲日韩中文| 国产高潮流白浆视频| 久久精品只有这里有| 日韩av在线直播| 欧美日在线观看| 野花国产精品入口| 91色爱欧美精品www| 亚洲欧洲天堂色AV| 国产欧美日韩在线一区| 国产网站在线看| 日韩123欧美字幕| 成人国产小视频|