• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    Novel structures and mechanical properties of Zr2N:Ab initio description under high pressures*

    2021-01-21 02:13:34MinruWen文敏儒XingXie謝興ZhixunXie謝植勛HuafengDong董華鋒XinZhang張欣FugenWu吳福根andChongYuWang王崇愚
    Chinese Physics B 2021年1期
    關(guān)鍵詞:張欣

    Minru Wen(文敏儒), Xing Xie(謝興), Zhixun Xie(謝植勛), Huafeng Dong(董華鋒),?, Xin Zhang(張欣),Fugen Wu(吳福根), and Chong-Yu Wang(王崇愚)

    1School of Physics and Optoelectronic Engineering,Guangdong University of Technology,Guangzhou 510006,China

    2School of Materials and Energy,Guangdong University of Technology,Guangzhou 510006,China

    3Department of Physics,Tsinghua University,Beijing 100084,China

    Keywords: phase transition,phonon dispersion,Zr2N,first-principles calculations

    1. Introduction

    As key materials for cutting coatings, super wearresistance, and thermal barrier coatings, group IVB transition metal nitrides (TMNs; TM=Ti, Zr, and Hf) have been widely used in mold-making industries,automobile manufacturing, and aerospace industries owing to their high melting point, excellent chemical stability, corrosion resistance, and metallic conductivity.[1–4]In the TMNs,the well-known stable phase is the stoichiometric B1-TMN compound with a rocksalt structure(space group Fmm). However,the TMNs may display a variety of compositions and phases featuring both cation and nitrogen enrichment.[5–8]Owing to the improvement of high-temperature and high-pressure experimental conditions, novel TMNs such as Th3P4-Zr3N4[9]and I4/mcm-TiN2[10]have been experimentally synthesized and exhibit excellent mechanical properties. With metalloid vacancies displaying long- or short-range order, or both, most TMNs(TMNx,for x=0–1)exist on the metal-rich side(i.e.,x<0.5)of the phase diagram.[5–7]However, compared to transition metal (TM) sub-nitrides, more research has focused on stoichiometric TMNs and nitrogen-rich TMNs. Consequently,an investigation of novel stable crystal structures for metal-rich TMNs under high pressure is essential for the development of practical applications.

    Numerous studies exist on the phase stability, electronic structure,and mechanical properties of Ti-N,[11–17]Zr-N,[11–16,18–20]and Hf-N[12–16,18,21]systems. Recent studies using the evolutionary algorithm/search algorithm and density functional theory(DFT)have investigated the stable structures and possible stoichiometries of the TMN[15–17,20,21]systems. Yu et al.[17]explored the stable structures and possible stoichiometries of the Ti-N system at the pressures of 0 GPa, 20 GPa, and 60 GPa and found the new ground structures Ti3N2, Ti4N3, and Ti6N5at atmospheric pressure and two novel stable phases Cmcm-Ti2N and I4/mcm-TiN2at high pressures (synthesized by Bhadram et al. using the laser-heated diamond anvil cell technique at 73 GPa and 2400 K[10]).Weinberger et al.[15]investigated candidate stable structures of the group IVB and group VB TMNs at ambient pressure and demonstrated that the stable phase of Ti2N and Zr2N is the tetragonal P42/mnm phase, while Hf2N is more energetically favorable in the orthorhombic Pnnm phase. Research by Yu et al.[20]revealed that Fmm-ZrN, P42/mnm-Zr2N,C2/m-Zr4N3,C2/m-Zr6N5,C2/m-Zr8N7,P-Zr15N16,P-Zr7N8, and C2/m-Zr4N5are ground-state structures at 0 GPa and low temperatures. Zhang et al.[21]studied the pressure-composition phase diagram for the Hf-N system under 0–60 GPa and revealed new candidate phases with a Hf:N ratio of 6:1,3:1,and 3:2 at zero pressure.Zirconium as well as Ti and Hf are group IVB TMs whose nitrides display similar phase transitions,chemical bonding,and elasticity. Unlike many studies for titanium and hafnium nitrides, the investigations of novel structural phases of high-pressure zirconium nitrides are limited. Specifically, no theoretical study has yet focused on new crystal structures of Zr-rich nitrides under high pressures.

    Considering this background, the present research focuses on novel crystal configurations of zirconium sub-nitride Zr2N at 0–200 GPa by combining evolutionary searches and first-principles calculations. By calculating the elastic constants and phonon spectra,the mechanical and dynamical stabilities for all candidate structures are evaluated. Additionally,the electronic and mechanical properties of the crystal phases are systematically analyzed to reveal the potential applications of the binary zirconium heminitrides.

    2. Computational methods

    To identify all potential novel Zr2N crystal structures,the evolutionary algorithm method implemented in the USPEX code[22–24]with variable cell number (up to 30 atoms)is performed at pressures of 0 GPa,30 GPa,60 GPa,90 GPa,120 GPa,150 GPa,and 200 GPa. We set the number of structures in the population to be 100,thus producing at random the first generation of structures. The 50 structures in the subsequent generations are produced by heredity(40%),soft mutation(20%),transmutation(20%),and randomness(20%).

    The structural relaxations (i.e., optimization of lattice shape, volume, and atomic positions) and energies of every structure along with all calculations of the properties of interest are executed in the VASP code[25]based on DFT.[26,27]The projector augmented-wave method[28]is adopted for the ionelectron interactions,and the generalized gradient approximation of Perdew et al.[29]is used for the exchange-correlation function. Electrons in the 4s24p64d25s2and 2s22p3orbitals are considered as valence electrons for Zr and N, respectively. The plane-wave kinetic energy cutoff is set to 600 eV.Monkhorst-Pack[30]k-meshes with a resolution of 2π×0.04 ?A are adopted to characterize the energy integration in the Brillouin zone. The criteria for the force on the atoms are set at 0.01 eV/?A for ion self-consistency, while the energy convergence of the electronic self-consistency is set to be 10-5eV per atom.

    Using a 2×2×2 supercell, the phonon dispersion relations are achieved in the framework of density functional perturbation theory[31,32]as implemented in the PHONOPY package[33]as well as the VASP code. The elastic constants are calculated by the strain–stress method,[34]and the elastic moduli (i.e., bulk modulus B, shear modulus G, and Young’s modulus E)as well as Poisson’s ratio ν,are obtained under the Voigt–Reuss–Hill approximation.[35–37]Employing Chen’s model,[38]the theoretical Vickers hardness values are estimated using

    where κ =G/B is the Pugh ratio. The directional dependence of Young’s modulus is also calculated to reveal the anisotropy of all potential Zr2N phases.

    3. Results and discussion

    Fig. 1. Schematic crystal structures of potential Zr2N for (a) P42/mnm, (b) Pnnm, (c) Cmcm, (d) P4/nmm, (e) I4/mcm, (f) P21/m,(g)Pm1,and(h)C2/m phases at 0 GPa. The schematic shows the N atoms(gray circles)in the centers of the gray polyhedrals and the Zr atoms(green circles).

    As displayed in Fig. 1, although the Cmcm, P4/nmm,and P21/m phases belong to different crystal systems and have different symmetries, they demonstrate similar characters wherein all N atoms are octahedrally coordinated by six Zr atoms and the NZr6octahedrons arrange along a uniform direction. In fact,the Cmcm,P4/nmm,and P21/m phases are versions of the rocksalt B1-TiN(space group Fmm,No.225)with ordered N-vacancies. Likewise,vacancy-ordered phases have also been reported for group IVB TMNs[17,21]and carbides.[39,40]In the Pm1 and C2/m phases, N atoms are also centered in the octahedron(composed of Zr atoms)interstices, though NZr6octahedrons arrange along different orientations. There are 12 atoms in a unit cell of the anti-rutile structure of the I4/mcm phase,of which the metal Zr and nonmetal N atoms occupy the Wyckoff 8h and 4a positions, respectively. Each N atom is coordinated by eight Zr atoms,forming a decahedron.This crystal structure naturally reminds us of the I4/mcm-TiN2phase belonging to the same space group,where the metal atoms and nitrogen atoms in I4/mcm-TiN2and I4/mcm-Zr2N phases occupy the opposite Wyckoff positions.However,unlike the short NN dumbbell bond length of I4/mcm-ZrN2,which is as short as 1.424 ?A at 1 atm,[16]the Zr–Zr bond length of I4/mcm-Zr2N is as long as 2.981 ?A at ambient pressure. It can therefore be predicted that the mechanical strength and hardness of I4/mcm-Zr2N will be significantly less than those of I4/mcm-ZrN2.

    To ascertain the thermodynamic stability of the predicted Zr2N configurations, we obtain the formation enthalpies per atom ΔH using

    where H(Zr2N),H(Zr),and H(N2)are the enthalpies per formula unit of Zr2N, P63/mcm-Zr, and α-N2, respectively, at a specific pressure. Our determined formation enthalpies of all Zr2N candidate phases at 0 GPa are summarized in Table 1.It can be seen that P42/mnm-Zr2N is the most stable phase of zirconium heminitride at atmospheric pressure with the formation enthalpy calculated to be-1.2914 eV per atom,which agrees well with the result of-1.249 eV/atom from Ref. [21].As a distorted version of ε-Zr2N(see Figs.1(a)and 1(b)),we note that the ΔH of the Pnnm phase is only slightly higher than that of P42/mnm-Zr2N.As shown in Table 1,the I4/mcm phase exhibits the highest formation enthalpy among the eight phases, indicating that I4/mcm Zr2N is the least energetically favorable phase for zirconium heminitride at 0 GPa.

    Table 1. Structural parameters(lattice constants and atomic coordinates)and formation enthalpies ΔH of the candidate phases of Zr2N at 0 GPa,along with other DFT-predicted values.

    To further investigate the phase transition of Zr2N at high pressure, the enthalpies of these candidate configurations are calculated as a function of pressure up to 200 GPa. Figure 2 plots the calculated enthalpy(relative to the Cmcm phase)per atom as a function of pressure for all candidate phases and the inset shows a detail of the enthalpy-pressure curves for the 3–17 GPa pressure range. As demonstrated in Fig.2,the diverse crossings between different structures in the pressure range of 0–200 GPa indicate the existence of various structural phase transformations for Zr2N under pressure. The P42/mnm-Zr2N transforms to Pnnm-Zr2N at approximately 6.5 GPa, whereupon Pnnm-Zr2N transforms to Cmcm-Zr2N at 16.3 GPa.This is different from the scenario of Ti2N and Hf2N. Specifically,the stable crystal structure of titanium heminitride at the ambient pressure (i.e., P42/mnm-Ti2N) transforms directly to Cmcm-Ti2N at 20.8 GPa.[17]In terms of hafnium heminitride,the orthorhombic Pnnm-Hf2N phase,which is the most stable configuration at atmospheric pressure, transforms to Cmcm-Hf2N at 12.5 GPa.[21]According to Fig.2,though the Cmcm-Zr2N is initially the most energetically favorable when the pressure exceeds 16.3 GPa, other phases become more energetically favorable as the pressure increases further. For example,the P4/nmm-Zr2N phase exhibits a lower enthalpy than Cmcm-Zr2N at 55.2 GPa, as does the P21/m-Zr2N phase at 65.0 GPa, where the P4/nmm-Zr2N phase possesses the lowest enthalpy from 55.2 GPa to 144.4 GPa. When the pressure further increases, the P4/nmm and P21/m phases both transform to the I4/mcm configuration at pressures of 122.8 GPa and 144.4 GPa, respectively. This signifies that the I4/mcm-Zr2N phase is the most energetically-favored structure in the 144.4–200 GPa pressure range. The enthalpies of the Pm1 and C2/m phases are higher than those of the Pnnm phase during the entire studied pressure range,but these phases possess a lower enthalpy than the P42/mnm phase when the pressure exceeds 34.5 GPa and 40.3 GPa,respectively.

    Fig. 2. Calculated enthalpy curves as a function of pressure for the different phases of Zr2N,with reference to the Cmcm phase. Inset: Detail of the enthalpy–pressure curves in the 3–17 GPa pressure range.

    Fig. 3. Phonon dispersion curves for (a) P42/mnm, (b) Pnnm, (c) Cmcm, (d) P4/nmm, (e) I4/mcm, (f) P21/m, (g) Pm1, and (h)C2/m phases at 0 GPa,respectively.

    Next,the dynamical and mechanical stabilities of the candidate Zr2N structures are investigated using the phonon spectra and elastic constants.The calculated phonon frequencies of the five novel structures(P4/nmm,I4/mcm,P21/m,Pm1,and C2/m)and the P42/mnm,Pnnm,and Cmcm phases of Zr2N at atmospheric pressure are plotted in Fig. 3. Clearly, there are no imaginary frequencies for all predicted phases at 0 GPa,signifying that these configurations are all dynamically stable.The calculated phonon dispersion curves of P42/mnm-Zr2N in this study are consistent with previous DFT results.[20]In addition to ambient pressure,we also study the phonon stability of all structures under high pressure. Soft phonons are found in ε-Zr2N,P4/nmm-Zr2N,and C2/m-Zr2N with an initial negative slope at the Γ point(related to crystal instability under a homogeneous deformation)for pressures of 20 GPa,200 GPa,and 120 GPa,respectively. This suggests the existence of mechanical instability in ε-Zr2N,P4/nmm-Zr2N,and C2/m-Zr2N at these corresponding pressures.[41]The Cmcm,I4/mcm,and P21/m phases,however,maintain their phonon stability above 200 GPa, while the Pnnm and Pm1 configurations hold dynamic stability up to 200 GPa and 160 GPa,respectively.

    The elastic constants were obtained using the strainstress method and are summarized in Table 2. The P42/mnm,P4/nmm,and I4/mcm structures belong to the tetragonal crystal system (Laue class of 4/mmm) and have six independent elastic constants. The mechanical stability criteria of the tetragonal crystal are given as[42–44]

    The Pnnm- and Cmcm-Zr2N are both orthorhombic crystal systems with nine independent elastic constants. The mechanical stability criteria for an orthorhombic class are[45]

    Table 2. Calculated elastic constants(Cij),elastic moduli(i.e.,bulk modulus BH,shear modulus GH,and Young’s modulus E),B/G ratio,Poisson’s ratio(v),and Vickers hardness(Hv)of the Zr2N compounds at zero pressure,along with other DFT calculated values. Ci j,BH,GH,E,and Hv are in GPa,B/G,and v are dimensionless.

    The constants C11,C22, and C33reflect the resistance of the crystal to normal strain in the[100],[010],and[001]directions,respectively. As mentioned above,the Cmcm,P4/nmm,and P21/m phases are versions of the cubic B1-TiN possessing ordered nitrogen vacancies.As shown in Table 2,P21/m-Zr2N possesses the lowest C11among the eight predicted structures.This is because the N-vacancies of the P21/m phase are distributed in the(100)plane and the least-dense atomic arrangement is in the [100] direction (see Fig. 1). Similarly, the Nvacancies of the P4/nmm phase are distributed in the (001)plane and the atomic compactness in the[001]direction is the worst (shown in Fig. 1). The P4/nmm-Zr2N phase possesses the lowest C33value among the eight candidate structures.Conversely, the Cmcm,C2/m, and P42/mnm phases respectively exhibit the highestC11,C22,andC33values.We note that our calculated elastic constants(except C11)of P42/mnm-Zr2N are consistent with the results of Yu et al.,[20]and our calculated C11of ε-Zr2N agrees well with the result of Ref. [48].

    The elastic moduli (i.e., bulk modulus BH, shear modulus GH,and Young’s modulus E),B/G ratio,Poisson’s ratio v,and Vickers hardness Hvof the Zr2N compounds at zero pressure are further studied by the methods mentioned in Secion 2 and are summarized in Table 2. Clearly, ε-Zr2N exhibits the largest bulk,shear,and Young’s moduli among the Zr2N compounds, which indicates that ε-Zr2N has the best resistance to volume change and shear shape deformation, and the best stiffness. Notably, the Pnnm phase, which is a distorted version of the ε-Zr2N phase,has similar but lower elastic moduli compared with those of the ε phase. The P4/nmm-Zr2N phase exhibits the smallest BHvalue,while I4/mcm-Zr2N possesses the smallest GHand E values.

    The B/G ratios (proposed by Pugh[49]) of all structures except the Cmcm phase are larger than 1.75, implying that the P42/mnm, Pnnm, P4/nmm, I4/mcm, P21/m, Pm1, and

    C2/m phases belong to ductile materials, where the I4/mcm phase is the most ductile among these configurations. In contrast,Cmcm-Zr2N possesses the smallest B/G ratio and Poisson’s ratio v,which is the most brittle phase among the Zr2N compounds.Intriguingly,owing to its brittleness,Cmcm-Zr2N exhibits the highest Vickers hardness (12.9 GPa), though its shear modulus is lower than that of ε-Zr2N. Contrarily, the I4/mcm phase has the lowest hardness as small as 2.1 GPa at ambient pressure. Note here that the I4/mcm phase has the lowest shear modulus,Young’s modulus,and hardness among the studied structures,where its abnormal and relatively higher bulk modulus primarily originates from its large mass density.

    Moreover, we investigate the directional dependence of Young’s modulus for all predicted Zr2N phases (see Fig. 4).Although the eight candidate configurations possess identical chemical composition ratios, their mechanical properties exhibit considerable differences. Young’s modulus of the Pm1 phase displays the most isotropic features among the Zr2N compounds because its three-dimensional Young’s modulus diagram shape is the closest to a sphere. Looking at the diagram shapes in Fig.4 and the degree of their deviation from spherical,the P21/m-Zr2N phase exhibits the most anisotropic features, followed by the Pnnm and I4/mcm phases. The anisotropy of the P21/m-Zr2N arises from the high C33and C22values and the low C23value,resulting in a large Young’s modulus in the y- and z-axis directions and a small E in the diagonal directions of the y–z plane. The anisotropy of the Pnnm and I4/mcm phases is mainly owing to their relatively low C11and C22values and high C12value,resulting in a small Young’s modulus in the x-and y-axis directions and a large E in the diagonal directions of the x–y plane.

    To understand the electronic properties and bonding features of these zirconium heminitrides, we further study the density of states (DOS) of each system. Figure 5 shows the s- and p-orbital partial DOS of N, the d-orbital partial DOS of Zr, and the total DOS for all Zr2N phases. Notably, all Zr2N compounds are metallic owing to the finite DOS at the Fermi level,which are dominated by the 4d orbitals of the zirconium atoms (see Fig. 5). Owing to the contribution of the 2s orbitals of the N atoms, all studied configurations exhibit a deep valence band below -12 eV. Similar to the scenario reported in the Tin+1Nn(n=1,2,3,and 5)structures,[17]interactions exist here between Zr and N atoms. According to Fig. 5, these interactions originate from the hybridization of the N-2p states and Zr-4d states in the energy region of-7 eV to -2.5 eV. We note that the interaction between Zr and N in the Cmcm phases exhibits greater hybridization than that in the P4/nmm and P21/m structures. This suggests that Cmcm-Zr2N will possess better mechanical strength and deformation resistance than the P4/nmm and P21/m phases, which is consistent with the calculated results of the elastic moduli in Table 2. Unlike the strong covalent N–N bonding nature in the I4/mcm-type TM pernitrides, weaker interactions are exhibited by the covalent-like Zr–N bond and metallic Zr–Zr bond in the I4/mcm-Zr2N phase. Consequently,compared with I4/mcm-type TM pernitrides,the I4/mcm-type TM heminitrides possess reduced mechanical strength and hardness.For example,compared with the values exhibited by I4/mcm-Zr2N, the I4/mcm-ZrN2structure exhibits larger calculated bulk modulus (1.7×), shear modulus (3.4×), Young’s modulus(3.1×),and hardness(12.3×).[16]

    Fig.4. Directional dependence of Young’s modulus(in GPa)of(a)P42/mnm,(b)Pnnm,(c)Cmcm,(d)P4/nmm,(e)I4/mcm,(f)P21/m,(g)Pm1,and(h)C2/m phases at 0 GPa.

    Fig.5. Projected density of states(DOS)of Zr2N compounds at ambient pressure,showing the partial DOS of the N 2s(red line)and N 2p(blue line)orbitals,the Zr 4d orbital(green line),and the total DOS(black line).

    4. Conclusion

    In summary, the pressure-dependent phase diagram of zirconium heminitride is investigated using ab initio evolutionary algorithm methods. Five novel high-pressure phases(i.e., P4/nmm, I4/mcm, P21/m, Pm1, and C2/m) of Zr2N are predicted in the present study. As the pressure increases from 0 GPa to 200 GPa,Zr2N moves sequentially through the phases as P42/mnm→Pnnm →Cmcm →P4/nmm→I4/mcm at the corresponding transition pressures of 6.5 GPa,16.3 GPa,55.2 GPa,and 144.4 GPa,respectively.The calculated phonon spectra reveal that all candidate Zr2N structures are dynamically stable at ambient pressure. According to the calculated elastic constants, all predicted Zr2N configurations are expected to be mechanically stable at atmospheric pressure.Although the eight candidate configurations possess identical chemical composition ratios, their mechanical properties exhibit considerable differences. Excepting the Cmcm phase,the Zr2N structures are ductile materials,among which the I4/mcm phase is the most ductile configuration. At atmospheric pressure, the Vickers hardness of the Zr2N compounds is ordered as Cmcm >P42/mnm >C2/m >P4/nmm >Pm1 >Pnnm >P21/m >I4/mcm. The calculated electronic structures demonstrate that the high mechanical strength and hardness of the Cmcm phase originate from the strong hybridization of the N-2p states and Zr-4d states, while the weaker interactions between the Zr and N atoms are likely responsible for the weaker mechanical behaviors of I4/mcm-Zr2N. These results provide important additional insight into the structureproperty relationships of TM heminitrides and the development of practical applications for Zr2N.

    Acknowledgment

    We thank Sara Maccagnano-Zacher, PhD, from Liwen Bianji, Edanz Editing China (www.liwenbianji.cn/ac), for editing the English text of a draft of this manuscript.

    猜你喜歡
    張欣
    《卯兔之年》
    《保護(hù)生態(tài)》
    《城》
    Boron at tera-Pascal pressures
    平面向量線(xiàn)性運(yùn)算的轉(zhuǎn)化思想的應(yīng)用
    隨筆四則
    作品(2020年4期)2020-05-11 06:21:45
    Dynamical stable-jump-stable-jump picture in a non-periodically driven quantum relativistic kicked rotor system?
    自作多情
    張欣現(xiàn)代重彩作品欣賞
    Estimating the clutch transmitting torque during HEV mode-switch based on the Kalman filter
    国产精品1区2区在线观看.| 亚洲不卡免费看| 亚洲国产日韩欧美精品在线观看| 久久国产乱子免费精品| 亚洲三级黄色毛片| 国产高清视频在线播放一区| 久久午夜福利片| 欧美不卡视频在线免费观看| 久9热在线精品视频| 日本精品一区二区三区蜜桃| 国内精品久久久久久久电影| 久久精品国产亚洲av香蕉五月| 亚洲中文日韩欧美视频| 最近最新免费中文字幕在线| 亚洲一区二区三区不卡视频| 亚洲自拍偷在线| 久久草成人影院| 高清日韩中文字幕在线| 91在线观看av| 性插视频无遮挡在线免费观看| 熟女电影av网| 日本撒尿小便嘘嘘汇集6| 日韩成人在线观看一区二区三区| 99热这里只有精品一区| 久久中文看片网| 国产成+人综合+亚洲专区| 性色avwww在线观看| 91字幕亚洲| 亚洲精品乱码久久久v下载方式| 国产免费一级a男人的天堂| 欧美色视频一区免费| 亚洲国产高清在线一区二区三| 国产精品女同一区二区软件 | 日韩欧美国产在线观看| 免费av毛片视频| 国产熟女xx| 中文字幕免费在线视频6| 人人妻人人看人人澡| 俺也久久电影网| 熟女电影av网| 亚洲一区高清亚洲精品| 欧美+亚洲+日韩+国产| 亚洲国产色片| eeuss影院久久| 日日摸夜夜添夜夜添av毛片 | 女人被狂操c到高潮| 亚洲内射少妇av| 动漫黄色视频在线观看| 1024手机看黄色片| 久久精品国产自在天天线| 少妇的逼水好多| 亚洲欧美清纯卡通| 婷婷精品国产亚洲av在线| 亚洲欧美日韩高清在线视频| 欧美成人一区二区免费高清观看| 亚洲专区中文字幕在线| 国内久久婷婷六月综合欲色啪| 久久精品影院6| 国内精品美女久久久久久| 亚洲第一欧美日韩一区二区三区| 久久热精品热| 久久国产精品人妻蜜桃| 成人高潮视频无遮挡免费网站| 18禁黄网站禁片免费观看直播| 老师上课跳d突然被开到最大视频 久久午夜综合久久蜜桃 | 三级毛片av免费| 国产欧美日韩一区二区三| 波多野结衣巨乳人妻| 非洲黑人性xxxx精品又粗又长| 亚洲性夜色夜夜综合| 又紧又爽又黄一区二区| 国产日本99.免费观看| 亚洲国产欧美人成| 老熟妇仑乱视频hdxx| 免费观看精品视频网站| 老鸭窝网址在线观看| 国产一级毛片七仙女欲春2| 青草久久国产| 熟女电影av网| 亚洲三级黄色毛片| 国产av不卡久久| 中文亚洲av片在线观看爽| 日韩有码中文字幕| 国产欧美日韩精品亚洲av| 亚洲成人精品中文字幕电影| 久久久精品大字幕| 美女 人体艺术 gogo| 女生性感内裤真人,穿戴方法视频| 波野结衣二区三区在线| 色综合婷婷激情| 亚洲aⅴ乱码一区二区在线播放| 国产白丝娇喘喷水9色精品| 日本三级黄在线观看| 久久6这里有精品| 国产三级中文精品| 黄色丝袜av网址大全| 深爱激情五月婷婷| 免费无遮挡裸体视频| 色综合站精品国产| 国产成人啪精品午夜网站| 午夜老司机福利剧场| 日韩欧美三级三区| 永久网站在线| 老熟妇乱子伦视频在线观看| 亚洲人成电影免费在线| 日韩中字成人| 一本综合久久免费| 国产麻豆成人av免费视频| 听说在线观看完整版免费高清| 我的老师免费观看完整版| 免费电影在线观看免费观看| 国产欧美日韩精品一区二区| 亚洲av一区综合| 两个人的视频大全免费| 天堂影院成人在线观看| 少妇裸体淫交视频免费看高清| 全区人妻精品视频| 中出人妻视频一区二区| 日本精品一区二区三区蜜桃| 欧美3d第一页| 成人av一区二区三区在线看| 最新在线观看一区二区三区| 88av欧美| av天堂在线播放| 真人做人爱边吃奶动态| 九九热线精品视视频播放| 亚洲欧美日韩无卡精品| 少妇的逼水好多| 午夜老司机福利剧场| 国产伦精品一区二区三区视频9| 亚洲无线观看免费| x7x7x7水蜜桃| 我要搜黄色片| 精品不卡国产一区二区三区| or卡值多少钱| av天堂在线播放| 亚洲人成网站高清观看| 最近在线观看免费完整版| 国产成年人精品一区二区| 99国产精品一区二区三区| 桃色一区二区三区在线观看| 精品久久国产蜜桃| av在线蜜桃| 国产精品人妻久久久久久| 免费人成视频x8x8入口观看| 日本免费一区二区三区高清不卡| 性插视频无遮挡在线免费观看| 日日摸夜夜添夜夜添小说| 给我免费播放毛片高清在线观看| 男人舔奶头视频| 久久精品国产亚洲av涩爱 | 能在线免费观看的黄片| 18禁在线播放成人免费| 免费高清视频大片| 黄色日韩在线| 午夜免费男女啪啪视频观看 | 亚洲一区二区三区不卡视频| 很黄的视频免费| or卡值多少钱| 国产精品久久久久久久电影| 亚洲国产色片| 97人妻精品一区二区三区麻豆| 久久人妻av系列| 偷拍熟女少妇极品色| 波多野结衣高清无吗| 亚洲成av人片在线播放无| 日本黄色片子视频| 性插视频无遮挡在线免费观看| 热99在线观看视频| 亚洲精品在线美女| 日韩欧美国产在线观看| 十八禁人妻一区二区| 欧美bdsm另类| 一区二区三区激情视频| 人妻丰满熟妇av一区二区三区| 国产精品美女特级片免费视频播放器| 亚洲狠狠婷婷综合久久图片| 国产精品影院久久| 99久久99久久久精品蜜桃| 欧美成人免费av一区二区三区| 午夜福利在线观看吧| 国产一区二区三区视频了| 国产精品久久久久久久久免 | 精品国内亚洲2022精品成人| 久久性视频一级片| 一区二区三区四区激情视频 | 国产亚洲精品综合一区在线观看| bbb黄色大片| 听说在线观看完整版免费高清| 亚洲av美国av| 婷婷精品国产亚洲av| 淫妇啪啪啪对白视频| 亚洲欧美精品综合久久99| 亚洲aⅴ乱码一区二区在线播放| 听说在线观看完整版免费高清| 欧美黄色片欧美黄色片| 亚洲欧美日韩无卡精品| 欧美色欧美亚洲另类二区| 亚洲欧美日韩卡通动漫| 俄罗斯特黄特色一大片| 欧美午夜高清在线| 日韩欧美 国产精品| 赤兔流量卡办理| aaaaa片日本免费| 亚洲aⅴ乱码一区二区在线播放| avwww免费| av福利片在线观看| 国产综合懂色| 精品人妻熟女av久视频| 一区二区三区四区激情视频 | 久久久成人免费电影| 天堂av国产一区二区熟女人妻| 成人亚洲精品av一区二区| 欧美xxxx黑人xx丫x性爽| 亚洲内射少妇av| 国产精品亚洲一级av第二区| 啦啦啦韩国在线观看视频| 99热这里只有是精品在线观看 | 欧洲精品卡2卡3卡4卡5卡区| 性色avwww在线观看| 国产精品爽爽va在线观看网站| 亚洲精品在线观看二区| 12—13女人毛片做爰片一| 永久网站在线| 亚洲第一区二区三区不卡| 久久久久免费精品人妻一区二区| 中文字幕久久专区| 欧美日韩亚洲国产一区二区在线观看| 国产成人av教育| 制服丝袜大香蕉在线| 亚洲内射少妇av| 久久这里只有精品中国| 亚洲一区二区三区色噜噜| 日韩国内少妇激情av| 精品久久久久久久人妻蜜臀av| 丰满人妻熟妇乱又伦精品不卡| 一级黄片播放器| 嫩草影院入口| a级毛片免费高清观看在线播放| 国产伦精品一区二区三区视频9| 日韩成人在线观看一区二区三区| 国产精品一区二区三区四区免费观看 | 18禁在线播放成人免费| 久久精品国产亚洲av香蕉五月| 国产高清激情床上av| 免费在线观看影片大全网站| 熟女电影av网| 日本黄大片高清| 午夜两性在线视频| 久久天躁狠狠躁夜夜2o2o| 成年女人看的毛片在线观看| 欧美绝顶高潮抽搐喷水| 人人妻人人看人人澡| 特大巨黑吊av在线直播| 一边摸一边抽搐一进一小说| 成年女人看的毛片在线观看| 久久久久久久久中文| 久久久久久久久中文| 久久草成人影院| 中文字幕精品亚洲无线码一区| 又爽又黄无遮挡网站| 18禁黄网站禁片午夜丰满| 亚洲国产精品成人综合色| 免费观看人在逋| 午夜福利免费观看在线| 亚洲精品色激情综合| www.999成人在线观看| 丰满人妻一区二区三区视频av| 我的女老师完整版在线观看| 亚洲自偷自拍三级| 欧美黄色片欧美黄色片| 国产成年人精品一区二区| 十八禁网站免费在线| 天堂动漫精品| 女生性感内裤真人,穿戴方法视频| 99热只有精品国产| av在线蜜桃| 亚洲最大成人中文| 亚洲无线观看免费| 精品国内亚洲2022精品成人| 一级作爱视频免费观看| 午夜亚洲福利在线播放| 一个人免费在线观看的高清视频| 五月玫瑰六月丁香| 久久精品国产清高在天天线| 免费看a级黄色片| 亚洲真实伦在线观看| 久久久久久久久中文| 夜夜躁狠狠躁天天躁| 亚洲avbb在线观看| 一级av片app| 亚洲欧美激情综合另类| 非洲黑人性xxxx精品又粗又长| 日韩有码中文字幕| 日韩亚洲欧美综合| 中国美女看黄片| 日韩av在线大香蕉| 国产精品久久久久久人妻精品电影| 欧美成人a在线观看| 免费搜索国产男女视频| 最近中文字幕高清免费大全6 | 很黄的视频免费| 一进一出抽搐动态| 丰满的人妻完整版| 男人舔奶头视频| 天堂动漫精品| 99在线人妻在线中文字幕| 色在线成人网| 久久精品人妻少妇| 国产三级在线视频| av天堂在线播放| 久久国产乱子免费精品| 黄色一级大片看看| 在线a可以看的网站| 尤物成人国产欧美一区二区三区| 亚洲第一电影网av| 老女人水多毛片| 淫妇啪啪啪对白视频| 国产成人a区在线观看| 国产aⅴ精品一区二区三区波| 直男gayav资源| 内地一区二区视频在线| 欧美3d第一页| 亚洲av电影不卡..在线观看| 成人性生交大片免费视频hd| 日本 欧美在线| 美女高潮喷水抽搐中文字幕| 精品一区二区免费观看| 日韩成人在线观看一区二区三区| 国产一区二区在线观看日韩| 18禁黄网站禁片免费观看直播| 欧美在线黄色| 日本精品一区二区三区蜜桃| 国产av不卡久久| 久久久精品欧美日韩精品| 亚洲av熟女| 午夜福利视频1000在线观看| 51国产日韩欧美| 51国产日韩欧美| 国产精品美女特级片免费视频播放器| 如何舔出高潮| 国产高清三级在线| 在线播放无遮挡| 久久人妻av系列| 脱女人内裤的视频| 在线播放无遮挡| 久久久久国产精品人妻aⅴ院| 免费无遮挡裸体视频| 国产精品久久久久久久电影| 熟女电影av网| 色综合婷婷激情| 久久久久久久久久成人| 国产乱人伦免费视频| 欧美日韩黄片免| 欧美又色又爽又黄视频| 小说图片视频综合网站| 国产免费一级a男人的天堂| 久久国产乱子免费精品| 黄色女人牲交| 亚洲18禁久久av| 极品教师在线免费播放| 免费搜索国产男女视频| 色精品久久人妻99蜜桃| 亚洲av第一区精品v没综合| 我要搜黄色片| 国产精品嫩草影院av在线观看 | 国产精品自产拍在线观看55亚洲| 嫩草影院新地址| 午夜免费激情av| 亚洲欧美激情综合另类| 级片在线观看| 国产伦一二天堂av在线观看| 伊人久久精品亚洲午夜| 动漫黄色视频在线观看| 亚洲乱码一区二区免费版| 婷婷色综合大香蕉| h日本视频在线播放| 成人av一区二区三区在线看| 一区福利在线观看| 在线免费观看的www视频| 有码 亚洲区| 搡老妇女老女人老熟妇| 男女床上黄色一级片免费看| 久久精品综合一区二区三区| 欧美成人性av电影在线观看| 成人无遮挡网站| 能在线免费观看的黄片| а√天堂www在线а√下载| 97热精品久久久久久| 亚洲人与动物交配视频| 欧美另类亚洲清纯唯美| 国内精品久久久久精免费| 中文字幕av在线有码专区| 久久精品国产亚洲av涩爱 | 欧美日本亚洲视频在线播放| 国产爱豆传媒在线观看| 精华霜和精华液先用哪个| 波多野结衣高清作品| 欧美性猛交╳xxx乱大交人| 高清毛片免费观看视频网站| 欧美日韩福利视频一区二区| 狠狠狠狠99中文字幕| 丰满乱子伦码专区| 一级毛片久久久久久久久女| 一级作爱视频免费观看| 中文字幕免费在线视频6| 蜜桃亚洲精品一区二区三区| 成年人黄色毛片网站| 很黄的视频免费| 亚洲人成网站在线播| 麻豆一二三区av精品| 精品不卡国产一区二区三区| 久久国产乱子伦精品免费另类| 国产高清视频在线观看网站| 性色av乱码一区二区三区2| 久久午夜福利片| 亚洲欧美日韩高清专用| 一级av片app| 欧美最黄视频在线播放免费| 伦理电影大哥的女人| 黄色日韩在线| 两性午夜刺激爽爽歪歪视频在线观看| 国产乱人伦免费视频| 精品人妻1区二区| 亚洲成a人片在线一区二区| 91麻豆av在线| 国产亚洲精品综合一区在线观看| 精品不卡国产一区二区三区| 国产三级在线视频| 国产私拍福利视频在线观看| 免费看美女性在线毛片视频| 亚洲av免费高清在线观看| 男人狂女人下面高潮的视频| 看免费av毛片| 国产乱人伦免费视频| 三级国产精品欧美在线观看| 在线观看午夜福利视频| 99精品久久久久人妻精品| 成人亚洲精品av一区二区| 亚洲午夜理论影院| 九九在线视频观看精品| 成人欧美大片| 动漫黄色视频在线观看| eeuss影院久久| 最近视频中文字幕2019在线8| 婷婷六月久久综合丁香| 精品国产三级普通话版| 美女xxoo啪啪120秒动态图 | 国产高清激情床上av| 色吧在线观看| 欧美一区二区亚洲| 青草久久国产| 中出人妻视频一区二区| 日日夜夜操网爽| 国产精品一区二区三区四区免费观看 | 日韩欧美在线乱码| 在线观看66精品国产| 久久国产精品影院| 日韩欧美国产在线观看| 十八禁国产超污无遮挡网站| 欧美激情国产日韩精品一区| 成人av一区二区三区在线看| 亚洲av.av天堂| 日本a在线网址| 一本一本综合久久| 午夜精品在线福利| www.熟女人妻精品国产| 内射极品少妇av片p| av天堂在线播放| 制服丝袜大香蕉在线| 色视频www国产| 国产在视频线在精品| 国产精品野战在线观看| 国产亚洲精品久久久久久毛片| 亚洲成人中文字幕在线播放| 男人的好看免费观看在线视频| 99久国产av精品| 人妻制服诱惑在线中文字幕| 久久人人爽人人爽人人片va | 黄色视频,在线免费观看| 在线观看午夜福利视频| 免费一级毛片在线播放高清视频| 色5月婷婷丁香| 亚洲,欧美,日韩| 亚洲真实伦在线观看| 俄罗斯特黄特色一大片| 中国美女看黄片| 麻豆av噜噜一区二区三区| 中文字幕高清在线视频| 国产免费一级a男人的天堂| 国内揄拍国产精品人妻在线| 国产伦人伦偷精品视频| 国产视频内射| 好看av亚洲va欧美ⅴa在| 欧美又色又爽又黄视频| 国产精品爽爽va在线观看网站| 日韩 亚洲 欧美在线| 男女那种视频在线观看| 国产激情偷乱视频一区二区| 亚洲在线自拍视频| 亚洲人成网站在线播放欧美日韩| 精品国产三级普通话版| 久久精品国产清高在天天线| 亚洲熟妇熟女久久| 1024手机看黄色片| 午夜久久久久精精品| 日韩精品中文字幕看吧| 一区二区三区高清视频在线| 亚洲无线在线观看| 热99re8久久精品国产| 看黄色毛片网站| 好男人电影高清在线观看| 国产av不卡久久| 真人一进一出gif抽搐免费| 成人一区二区视频在线观看| 99久久久亚洲精品蜜臀av| 国产精品,欧美在线| 真人一进一出gif抽搐免费| 亚洲在线自拍视频| 在线天堂最新版资源| 欧美黑人欧美精品刺激| 国产亚洲av嫩草精品影院| 国产色婷婷99| 亚洲专区中文字幕在线| 天堂av国产一区二区熟女人妻| 老师上课跳d突然被开到最大视频 久久午夜综合久久蜜桃 | 狂野欧美白嫩少妇大欣赏| 欧美日韩国产亚洲二区| 日韩成人在线观看一区二区三区| 99在线人妻在线中文字幕| 能在线免费观看的黄片| 免费av观看视频| 欧美bdsm另类| 十八禁人妻一区二区| 毛片女人毛片| 小蜜桃在线观看免费完整版高清| 又紧又爽又黄一区二区| 亚洲精品亚洲一区二区| 亚洲精品粉嫩美女一区| 欧美日韩福利视频一区二区| 在线播放无遮挡| 亚洲真实伦在线观看| 欧美成人性av电影在线观看| 国产精品久久久久久久久免 | 69人妻影院| 亚洲一区高清亚洲精品| 激情在线观看视频在线高清| 久久这里只有精品中国| 免费观看的影片在线观看| 丁香六月欧美| 国产午夜精品久久久久久一区二区三区 | 丝袜美腿在线中文| 校园春色视频在线观看| 国产亚洲精品久久久久久毛片| 一二三四社区在线视频社区8| 嫩草影院新地址| 色尼玛亚洲综合影院| 欧美三级亚洲精品| 露出奶头的视频| 成人亚洲精品av一区二区| 国产精品1区2区在线观看.| 波多野结衣高清无吗| 免费大片18禁| 精品人妻熟女av久视频| 中文字幕久久专区| 亚洲欧美日韩无卡精品| 中国美女看黄片| 亚洲不卡免费看| 国产免费一级a男人的天堂| 久久热精品热| 亚洲七黄色美女视频| 九色国产91popny在线| 我的老师免费观看完整版| 精品一区二区免费观看| 亚洲最大成人手机在线| 在线播放无遮挡| eeuss影院久久| 毛片女人毛片| 一本综合久久免费| 久久久久久久久大av| 啪啪无遮挡十八禁网站| 好男人在线观看高清免费视频| 亚洲午夜理论影院| 日韩免费av在线播放| 女生性感内裤真人,穿戴方法视频| 欧美另类亚洲清纯唯美| 国产精品国产高清国产av| 久久中文看片网| 能在线免费观看的黄片| 亚洲不卡免费看| 国产精品99久久久久久久久| www日本黄色视频网| 美女cb高潮喷水在线观看| 国产在视频线在精品| 日韩欧美国产一区二区入口| 中国美女看黄片| 国产熟女xx| 精品一区二区免费观看| 在现免费观看毛片| 免费看日本二区| 两个人的视频大全免费| 搡老妇女老女人老熟妇| 日韩亚洲欧美综合| 日本熟妇午夜| 88av欧美| 麻豆av噜噜一区二区三区| 中文字幕av在线有码专区| 亚洲av日韩精品久久久久久密| 亚洲男人的天堂狠狠| 亚洲专区国产一区二区| 一个人免费在线观看的高清视频| 欧美国产日韩亚洲一区| 国产在视频线在精品| 亚洲精品在线观看二区| 欧美黑人巨大hd|