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

    Theoretical study of novel B-C-O compounds with non-diamond isoelectronic

    2022-02-24 08:59:02ChaoLiu劉超andPanYing應(yīng)盼
    Chinese Physics B 2022年2期
    關(guān)鍵詞:劉超

    Chao Liu(劉超) and Pan Ying(應(yīng)盼)

    1College of Rare Earths,Faculty of Materials Metallurgy and Chemistry,Jiangxi University of Science and Technology,Ganzhou 341000,China2State Key Laboratory of Metastable Materials Science and Technology,Yanshan University,Qinhuangdao 066004,China3Hebei Key Laboratory of Microstructural Material Physics,School of Science,Yanshan University,Qinhuangdao 066004,China

    Two novel non-isoelectronic with diamond(non-IED)B—C—O phases(tI16-B8C6O2 and mP16-B8C5O3)have been unmasked. The research of the phonon scattering spectra and the independent elastic constants under ambient pressure(AP)and high pressure(HP)proves the stability of these non-IED B—C—O phases. Respective to the common compounds,the research of the formation enthalpies and the relationship with pressure of all non-IED B—C—O phases suggests that HP technology performed in the diamond anvil cell(DAC)or large volume press(LVP)is an important technology for synthesis.Both tI16-B8C6O2 and tI12-B6C4O2 possess electrical conductivity.mP16-B8C5O3 is a small bandgap semiconductor with a 0.530 eV gap. For aP13-B6C2O5,mC20-B2CO2 and tI18-B4CO4 are all large gap semiconductors with gaps of 5.643 eV,6.113 eV, and 7.105 eV, respectively. The study on the relationship between band gap values and pressure of these six non-IED B—C—O phases states that tI16-B8C6O2 and tI12-B6C4O2 maintain electrical conductivity, mC20-B2CO2 and tI18-B4CO4 have good bandgap stability and are less affected by pressure. The stress-strain simulation reveals that the max strain and stress of 0.4 GPa and 141.9 GPa respectively, can be sustained by tI16-B8C6O2. Studies on their mechanical properties shows that they all possess elasticity moduli and hard character. And pressure has an obvious effect on their mechanical properties, therein toughness of tI12-B6C4O2, aP13-B6C2O5, mC20-B2CO2 and tI18-B4CO4 all increases,and hardness of mP16-B8C5O3 continue to strengthen during the compression. With abundant hardness characteristics and tunable band gaps,extensive attention will be focused on the scientific research of non-IED B—C—O compounds.

    Keywords: stability analysis,pressure,formation enthalpy,property

    1. Introduction

    Since non-stoichiometric ratio B—C—O ternary compounds have been synthesized via high pressure (HP) experiments in multianvil apparatus with high temperature,[1—3]due to the formation of rigid covalent bonds, the scientific research on design of B—C—O ternary system with high hardness has attracted great enthusiasm.[4—7]As a typical representative (isoelectronic with diamond, short as IED, also the simplest ternary B—C—O compounds), B2CO compounds have become the first to be of concern among all the B—C—O ternary system compounds. After pioneering research,Liet al.exposed two superhard B2CO structures with high hardness about 50 GPa (tP4-B2CO and tI16-B2CO).[4]In fact, more ternary compounds exist in B—C—O systems including B2CXO(X=2,3,4,...),are IED.As the increase of C content will result in the stacking of C in the structure, and forming shorter and stronger sp3C—C bonds. Hence, Liet al.believe that B2CXO (X=2,3,4,...) are extremely likely to have higher hardness than B2CO.[4]

    Encouraged by Li’s research, Zhanget al.focused the theoretical research on the B2CXO (X ≥2) system and foretold three new superhard phases:[5]B2C2O with space group(SG)I41/amdis denoted as tI20-B2C2O, B2C3O with SGdenoted as tI12-B2C3O,and tP8-B2C5O for B2C5O with.[5]Based on the theoretical calculation model of hardness proposed by Gaoet al.,[8]tI20-B2C2O,tI12-B2C3O and tP8-B2C5O all have superhard nature with Vickers hardness value significantly exceeding 40 GPa (57 GPa, 62 GPa,and 68 GPa, respectively).[5]The relationship between carbon ratio via mechanical properties of B2CXO was established,and high carbon content improves both the elastic moduli and hardness,[5]which is well confirmed by Li’s conjecture.[4]

    Afterward, the B2CXO (X ≥1) series of compounds as the IED B—C—O have received constant attention. Different from those B2CXO (X ≥1) phases with tetragonal structures reported by predecessors,Liuet al.firstly proposed one superhard B2CO phase (oP8-B2CO) with an orthorhombic system,[9]which opens a new structural system for the design of high hardness materials of B—C—O system. Whereafter the physical properties of oP8-,tI16-,and tP4-B2CO under HP have been investigated by Qiaoet al.[10]Inspired by the structural characteristic of newly predicted B—C—O, both tP4-B2CO and tI16-B2CO have diamond-like structures,[4]oP8-B2CO possesses a lonsdaleite-like structure,[9]even tI20-B2C2O, tI12-B2C3O and tP8-B2C5O all share the diamondlike structures and can be derived from different tP4-B2COsupercells with partial B and O replaced by C.[5]Two orthorhombic B2CO phases(oP16-B2CO and oC16-B2CO)with superhard property have been predicted by four-step manual construction,[11]and they share the similar structures to Cco-C8 and Bct-C4,respectively.

    Even with great success in structure design,the most stable structure of B2CO in the ground state remains a mystery.To this end, Yanet al.carried out in-depth research and finally identified a novel orthorhombic structure oI16-B2CO as the B2CO thermodynamic ground-state phase via an unbiased structure searching method.[12]Different from these previously proposed phases,oI16-B2CO is formed by covalent sp2—sp3B—C and B—O bonds,and has a considerably lower energy than that of all others.[12]Encouraged by the first sp2and sp3hybridization coexist B2CO structure,a new metastable phase for B2CO (tP16-B2CO) with sp2—sp3coexisted was explored by Chenet al.[13]The calculation of mechanical properties indicates its hard properties with high hardness of 23.19 GPa and anisotropy with the max stress along[001]far larger than that along[100].[13]

    However, there are few studies on B—C—O compounds as non-IED. In 2016, Wanget al.initiated the study of nonisoelectronic system among B—C—O compounds and predicted a new superhard phase B4CO4with a body-centered tetragonal structure;[6]here we denoted it as tI18-B4CO4. tI18-B4CO4is metastable under ambient pressure (AP), but thermodynamically stable at HP above 23 GPa.[6]And the thermodynamic stability,mechanical and electronic properties of tI18-B4CO4were systematically studied by Qiaoet al.[14]and Nuruzzamanet al.[15]Also two low-enthalpy metastable compounds(B6C2O5with primitive triclinic structure, B2CO2with (0.5,0.5, 0) centered monoclinic) have been discovered.[6]For unification and ease of description, here named these two as aP13-B6C2O5and mC20-B2CO2, respectively. They can stably exist to HP as 20 GPa, and the hardness of aP13-B6C2O5is 29.6 GPa, which is relatively softer and more easily machinable compared to mC20-B2CO2(41.7 GPa).[16]Two years later, only one non-IED (tI12-B6C4O2) has been reported.[17]And the research about the mechanical properties reveals that tI12-B6C4O2possesses hard nature with hardness of 21.9 GPa.[17]Since then, there has none report of B—C—O compounds with non-IED.

    Previous theoretical studies have made good progress in predicting high-performance structures,but have not provided theoretical guidance for experimental synthesis. And then due to stoichiometric ratios B—C—O compounds have not been synthesized experimentally,the research of B—C—O compounds is limited to difficulties. In this paper,we focus on the non-IED B—C—O compounds. Undergoing rigorous structural stability analysis including mechanical and dynamical stabilities, two novel phases were discovered. Together with these previously reported structures,the optimized reactants and HP conditions were obtained based on the formation enthalpies calculation,which will provide guidance for HP experimental synthesis.Based on density function theory, the mechanical and electronic properties of all non-IED B—C—O compounds have been systematically studied,and the regulatory effect of pressure on properties has also been further investigated.

    2. Computational methods and details

    The potential structures of non-IED B—C—O compounds have been generated through evolutionary simulation methods CALYPSO[18—20]and USPEX.[21—23]The following full relaxation geometry optimizations were performed in CASTEP models.[24]The gradient-corrected functionals adopted the Perdew—Burke—Ernzerhof(PBE)form.[25]The ultrasoft pseudopotentials were adopted with energy cutoff as 380 eV andKpoints were generated by separation as 8π×10?3nm?1.[26]During the geometry optimizations,[27]iterations were continued until the energy did not change beyond 5×10?3meV per atom;force tensor on atoms was less then 0.1 eV/nm;displacement on atoms was no more than 5×10?5nm and stress did not exceed 20 MPa. The phonon calculation was implemented via the finite displacement method[28]and the elastic constants were calculated with 9 patterns generated for each strain and max strain amplitude as 0.009.

    The unit cells were adopted within the following research,and the symmetry points and their path for the Brillouin zone are displayed in Table S1. For the body-centered tetragonal(Bct) system as tI12-B6C4O2and tI18-B4CO4, all belongs share the same symmetry points and the path for the Brillouin zone.

    3. Results and discussion

    3.1. Optimized structure

    Two newly discovered B—C—O structures have been proposed and their crystal structures schematic at AP are shown in Fig.1. For the first compound with space group(SG)is a Bct structure containing 16 atoms per unit cell,hence denoted as tI16-B8C6O2. tI16-B8C6O2possess the same symmetry points with tI12-B6C4O2and tI18-B4CO4,as displayed in Table S1. The lattice parametersaandcare 3.650and 7.899,respectively. All B atoms share the equivalent atomic Wyckoff positions 8i(0.243, 0.243, 0.378); and all O atoms occupy 2b(0.5, 0.5, 0), while for C atoms, C1 2a(0, 0, 0);C2 4d(0.5, 0, 0.25). All the atoms formed covalent bonds with heteroatom by sp3hybridization, however O is not directly bonded to C. That is to say, only B—C and B—O bondsformed in tI16-B8C6O2. There are three different kinds of B—C bonds: type I bonds formed by B and C1 atoms with length of 1.578and type II bonds are formed by B and C2 atoms and length of 1.642. The B—O bonds share the uniform length of 1.640. The theoretical densityρof tI16-B8C6O2at AP is 3.007 g/cm3.

    Fig.1.Structure graphs for novel non-isoelectronic B—C—O:(a)tI16-B8C6O2 unit cell and(b)mP16-B8C5O3 (2×2×2)supercell at AP.

    For the second compound, it belongs to the primitive monoclinic due to the SG asPm, and the unit cell contains 16 atoms, hence denoted as mP16-B8C5O3. The lattice parametersa,bandcare 9.977, 2.638and 4.404, respectively. And the angelβis about 89.16°. Because the symmetry is extremely low, no class of atomic positions are equivalent(Table S2),hence the formed bonds are all different lengths ranging from 1.423to 1.701.For mP16-B8C5O3,it is also a typical structure with sp2and sp3hybridization coexisting. Thereinto,all C atoms are sp3hybridization with 4 B atoms or 3 B atoms and 1 C atom, all O atoms are sp2binding with 3 B atoms, while B atoms can be divided into two categories: sp3hybridized with 2 C atoms and 2 O atoms or 3 C atoms and 1 O atom, sp2binding with 3 C atoms or 2 C atoms and 1 O atom or 2 O atoms and 1 C atom. In addition, due to the high carbon content, there is also C—C bonds in mP16-B8C5O3, which is quite different from the sp2and sp3hybridization coexisting structures as oI16-B2CO[12]and tP16-B2CO.[13]The calculated densityρis 2.791 g/cm3at AP.

    3.2. Stability analysis at ambient pressure

    The elastic stabilities of tI16-B8C6O2and mP16-B8C5O3at AP are evaluated by the calculation of the independent elastic constantsCijs,as presented in Table 1.

    For tI16-B8C6O2, it belongs to tetragonal system with Laue class 4/mmm,Born criteria are listed as[29]

    For mP16-B8C5O3with Laue class as 2/m,belonging to a monoclinic system, it has 13 independent elastic constants.The generic necessary and sufficient criterion that all eigenvalues of the elastic matrix are positive is easy to check with simple linear algebra routines.[29]

    The elastic matrices are symmetric about the main diagonal, and this symmetric nonzero data is represented by the apostrophe symbol “...”, and the blank cell represents data that is 0. There is no doubt that from the above elastic matrix in Table 1, theCi js satisfy the criteria above, indicating tI16-B8C6O2and mP16-B8C5O3are elastic stable.

    Here theCijs of other four previously reported non-IED B—C—O phases at AP have also been calculated,and the corresponding matrices are listed in Table S3. Both mC20-B2CO2and aP13-B6C2O5have low symmetry, the generic necessary and sufficient criterion that all eigenvalues of C be positive is easy to check with simple linear algebra routines.[29]For tI12-B6C4O2, which belongs to tetragonal system with 4/mmmLaue class, the necessary and sufficient conditions for elastic stability are shown as Eq. (1). For tI18-B4CO4, which possess the tetragonal crystal system and Laue class as 4/m, its criterion for elastic stability can also be evidenced by[29]

    All the independent elastic constants Cijs in the matrices have successfully verified the stability of their elastic, which agree well with the previous research,[6,14,15,17]therefore also proved the accuracy of our calculation methods and parameters.

    Table 1. The elastic matrices of two newly non-IED B—C—O phases at AP,the numerical unit is unified with GPa.

    Fig.2. Phonon dispersion spectra and phonon DOS of(a)tI16-B8C6O2 and(b)mP16-B8C5O3 at AP.

    Fig. 3. Calculated formation enthalpies ΔHf via the pressure of six non-IED B—C—O compounds, (a)—(c) represent aP13-B6C2O5, mC20-B2CO2 and tI18-B4CO4,respectively. (d)The schematic diagram of LVP and DAC.

    The imaginary frequency will lead to crystal distortiom,indicating the dynamical instability. Phonon dispersion spectrum and the related phonon density of states (DOS) of two newly predicted and four reported non-IED B—C—O phases at AP are researched and plotted in Fig. 2 and Fig. S1, respectively. There are no negative phonon modes in the entire Brillouin zone of their unit cell, suggesting they are dynamically stable.

    The stability of the newly predicted structures as tI16-B8C6O2and mP16-B8C5O3was further confirmed by molecular dynamics(MD)simulation with a thermostat to maintain a constant temperature of 273 K and with a barostat to maintain a constant AP.The total simulation time is 3 ps and time step as 6 fs. The simulation results are shown in Fig.S2 and there was no significant change in the structure during the whole simulation process. That is,tI16-B8C6O2and mP16-B8C5O3maintain their structural integrity and immutability.

    3.3. Potential synthesis reaction

    For providing guidance of experimental synthesis,it is of great significance and extremely urgent to deploy the potential synthesis of novel non-IED B—C—O phases. Based on the common reactants such as boron carbide B4C,boron trioxide B2O3, boron-rich oxide B6O, graphite C, andα-boronα-B,constructing different reaction path,then balancing the synthesis equations of these novel non-IED B—C—O phases,therefore the relationships between formation enthalpy ΔHfand pressure are revealed according to the following Eqs. (3)—(13).The ΔHf1and ΔHf2represent the formation enthalpy via path 1 and path 2,respectively. For better comparison and description, these six non-IED B—C—O phases are divided into three groups: tI16-B8C6O2, mP16-B8C5O3; tI12-B6C4O2, aP13-B6C2O5; mC20-B2CO2(treat as B4C2O4), tI18-B4CO4due to the same B ratio.

    Figure S3 summarizes the relationships between formation enthalpy (normalized on a per formula) and the pressure of six non-IED B—C—O compounds. The enthalpy of formation is less than zero, indicating the likelihood that the target non-IED B—C—O compounds will be synthesized by the selected reaction, and more negative more possible to be synthesized. Hence the critical pressure for formation enthalpy reaching zero is the key to determining the HP synthesis.

    As shown in Fig. S3, the critical pressures of these six non-IED B—C—O compounds have been listed. For tI16-B8C6O2and mP16-B8C5O3, the critical pressures decrease with the C content decrease and O content increase. By comparing the reaction path 1 and path 2,it was found that B6O as one reactant will benefit the decrease of pressure required to react. For tI12-B6C4O2and aP13-B6C2O5, the critical pressures also have similar relationships of C/O content. For the group of mC20-B2CO2and tI18-B4CO4, the conclusion that the critical pressures decrease with the C content decrease are also very well illustrated. Among all the six non-IED B—C—O phases,the max critical pressure is just slightly over 100 GPa.As the working pressure and temperature of the diamond anvil cell (DAC) can reach hundreds of GPa and thousands of degrees centigrade, and the high temperature will benefit to reduce the reaction barrier, hence synergistic high temperature and HP experiments in DAC maybe a potentially efficient synthesis technology for the target B—C—O product at a lower pressure than the theoretical critical pressure value at zero temperature.

    For aP13-B6C2O5, mC20-B2CO2and tI18-B4CO4, they all have critical pressures lower than 30 GPa. The relationship between formation enthalpy and pressure in the working pressure range of large volume press (LVP) has been studied in detail. As shown in Fig. 3, for aP13-B6C2O5, mC20-B2CO2and tI18-B4CO4, their critical pressures are 16.5/19.7 GPa,21.4/28.8 GPa, and 4.9 GPa, respectively. This shows that these three non-IED B—C—O structures can be obtained via HP experiments on LVP devices. Of course,special technologies need to be combined,such as rapid pressure relief,rapid cooling and so on,which will help to preserve the new phases.

    Fig.4. The relationship of pressure and volume per formula of six non-IED B—C—O phases at AP to 120 GPa.

    3.4. Stability analysis at high pressure

    It is well known that pressure is an important parameter to regulate the state of matter. In general, solid materialsundergo a continuous densification during HP compression.As exhibited in Fig. 4, for these six non-IED B—C—O phases we studied, during the pressure from AP to HP of 120 GPa,the volumes per molecular formula of these phases decrease monotonically. Detailed research has revealed that with the same B proportion,tI16-B8C6O2has the smaller volume drop(24.63/formula) and volume compressibility (23.4%) than mP16-B8C5O3(35.83/formula, 30.9%); tI12-B6C4O2also has a smaller volume decrease (18.73/formula) and volume compressibility (23.6%) than aP13-B6C2O5(23.33per formula, 26.1%); mC20-B2CO2also shares the advantage of resisting volume change(7.43/formula)and volume compressibility (22.6%) than tI18-B4CO4(15.53/formula,25.3%). For mP16-B8C5O3,the volume varies with the pressure in two intervals,with a rapid decline before 40 GPa and a continuous slow decline after 40 GPa.

    Fig.5. Calculated phonon and DOS of serval non-IED B—C—O phases at 100 GPa HP;(a)—(e)represent tI12-B6C4O2,aP13-B6C2O5,mC20-B2CO2,tI18-B4CO4 and mP16-B8C5O3,respectively.

    As the HP technology is a potential way to synthesize these non-IED B—C—O phases,it is necessary to analysis their stability at HP. The elastic matrix of six non-IED B—C—O phases at HP of 100 GPa are shown in Table S4. All the independent elastic constantCi js satisfy the criteria as Eqs. (1)and (2), or all eigenvalues of elastic matrix are positive, indicating tI16-B8C6O2, mP16-B8C5O3, tI12-B6C4O2, aP13-B6C2O5,mC20-B2CO2and tI18-B4CO4keep elastic mechanical stability at 100 GPa HP.

    Their phonon scattering at HP also have been investigated. Figure 5 shows the phonon scattering spectra and DOS of tI12-B6C4O2, aP13-B6C2O5, mC20-B2CO2, tI18-B4CO4and mP16-B8C5O3. No negative frequency is found in the entire Brillouin zone,indicating that these five non-IED B—C—O phases remain stable at 100 GPa. However, as the phonon scattering spectrum and DOS at different pressure exhibited in Fig. S4, for tI16-B8C6O2, when pressure surpasses 80 GPa,negative frequency appears in the phonon spectrum,resulting in HP instability of tI16-B8C6O2. While the pressure is within 80 GPa,tI16-B8C6O2is still stable(no virtual frequency).The independent elastic constant Ci js of tI16-B8C6O2at HP of 40 GPa and 80 GPa have been researched(Table S5),and the related research results show the elastic mechanical stability of tI16-B8C6O2under HP less than 80 GPa. Hence, combining the study of phonon spectrum and elastic constants, the stability of tI16-B8C6O2at HP un-beyond 80 GPa is thoroughly proved.

    In fact,the phonon vibration frequency is affected by the strength of chemical bonds. That is the max phonon vibration frequency is an indication of the strength of the chemical bond. Furthermore, zero-point energy is the energy that a quantum will vibrate under the zero point of absolute temperature. The max phonon frequency tells you how strong the bond is. The amplitude of zero-point vibration increases with increasing temperature, and the lighter the atom (that is, the smaller the element number), the more pronounced the zero point vibration. Based on the phonon scattering spectrum and its DOS, zero-point energy can be calculated. Hence, the relationship of max vibration frequency(zero-point energy)and pressure are further researched based on the phonon scattering spectrum and the related DOS at series pressures.

    As displayed in Fig. 6, both tI16-B8C6O2and mP16-B8C5O3have max vibration frequency larger than 30 THz,and their max vibration frequencies increases monotonically with increasing pressure. Therein, tI16-B8C6O2has a larger growth range (8.003 THz) and the growth rate (26.1%) than mP16-B8C5O3(6.649 THz, 16.9%). By the comparison of zero-point energy, it was found that even with similar molecular molar mass, the zero point energies of tI16-B8C6O2and mP16-B8C5O3are significantly different, which may be mainly contributed by the structural difference including the density, porosity, and so on. The zero-point energies of tI16-B8C6O2and mP16-B8C5O3continuously increase with the pressure increasing, and the increase range and growth rate of these two phases are 0.578 eV/31.3%and 0.437 eV/20.3%,respectively.

    Fig.6. The relationship of pressure and(a)max vibration frequency and(b)zero-point energy of tI16-B8C6O2 and mP16-B8C5O3 at AP to 80 GPa.

    3.5. Electrical property

    Based on the GGA-PBE functional, the electronic band structures of these six non-IED B—C—O phases along the symmetry points of the Brillouin zone at AP are quickly calculated and displayed in Fig. S5. For tI16-B8C6O2, mP16-B8C5O3and tI12-B6C4O2, there are energy bands that pass through the Fermi level, which indicates that they all possess electrical conductivity. As for the other three, there is no band across Fermi level, and the forbidden bands are formed by valence band max(VBM)and conduction band min(CBM). Hence, for aP13-B6C2O5, mC20-B2CO2, and tI18-B4CO4, they all have semiconductor properties with the gap of 3.927 eV,4.650 eV,and 5.434 eV,respectively.

    Although based on the GGA-PBE functional, the electronic band structure can be quickly calculated. However,based on the GGA-PBE functional, the calculated band gap trend to be underestimated. Hence, the more accurate but more time-consuming calculation based on the hybrid functionals HSE06[30,31]is adopted. As shown in Fig. 7, after the rigorous testing of HSE06, both tI16-B8C6O2and tI12-B6C4O2show good electrical conductivity due to the distinct bands cross the Fermi level. For aP13-B6C2O5,mC20-B2CO2and tI18-B4CO4,they all keep semiconductor properties. The max value points in VBM and the min value point in CBM of aP13-B6C2O5are located at high symmetry pointsFandQ, respectively. Hence, aP13-B6C2O5is an indirect bandgap semiconductor with 5.643 eV. For mC20-B2CO2, max VBM locates atZpoint and min CBM does not occupy high symmetry points. It is also an indirect bandgap semiconductor with a gap of 6.113 eV. For tI18-B4CO4, max VBM occupiesGpoints,min CBM locates atMpoint and the value of indirect band gap is as high as 7.105 eV. While, for mP16-B8C5O3,it is an indirect semiconductor with a small gap as 0.530 eV,which is quite different from the conclusion based on GGAPBE. The max VBM and min CBM are located atBandE,respectively.

    Fig.7.Calculated electronic band structures via HSE06 of non-IED B—C—O phases at AP;(a)—(f)represent tI16-B8C6O2,mP16-B8C5O3,tI12-B6C4O2,aP13-B6C2O5,mC20-B2CO2 and tI18-B4CO4,respectively. The horizontal red line,dark cyan curve and gold curve represent the Fermi energy level,VBM and CBM,respectively.

    The relationship between the band gap values and pressure of these six non-IED B—C—O phases was researched.Considering the high computational cost of HSE06 and the research purpose is to reveal the trend of band gap variation with pressure,GGAPBE calculation is adopted here. Both tI16-B8C6O2and tI12-B6C4O2maintain electrical conductivity over the entire pressure range which they are dynamics and elasticity are stable. For the others four phases,the relationship in the range of 0—100 GPa is plotted in Fig.8. For mP16-B8C5O3,it keeps a zero gap until 28 GPa,then the gap rapid increases,however after 36 GPa(gap as 1.113 eV),the band gap sluggishly increases to 1.206 eV with the pressure to 100 GPa. For aP13-B6C2O5,its gap increasessignificantly in the range of AP (3.927 eV) up to 6 GPa (4.227 eV), then the gap increases slowly to 4.330 eV at 100 GPa.However,for mC20-B2CO2and tI18-B4CO4,their gap values both have very slight change during the entail pressure range,the gap value of mC20-B2CO2first goes uphill to 4.695 eV at 10 GPa then downhill to 4.454 eV at 100 GPa,and the gap value of tI18-B4CO4slight increases to 5.496 eV at 50 GPa and then negligible decreases to 5.472 eV at 100 GPa.

    Fig.8. (a)Calculated band gaps as a function of pressure for mP16-B8C5O3,aP13-B6C2O5,mC20-B2CO2 and tI18-B4CO4 based on GGAPBE.Panels(b)and(c)correspond to the dotted areas in(a)with more smaller intervals samples.

    3.6. Mechanical property

    The material’s ideal strength,which is defined as the critical stress when a perfect crystal lost its elastic stability,equal to the upper limit for material strength. The atomistic mechanism for structural deformation and failure models can be clearly interpretated with the investigation of strain—stress relations and bond-breaking processes.[32]

    Figure 9 presents the simulated of tensile stress—strain relationships along the specific directions including[100],[110]and [001] for tI16-B8C6O2. It can be seen that tI16-B8C6O2has the max tensile stress of 141.9 GPa in the [001] direction with the strain at 0.4,which is larger than the max tensile stress (123.1 GPa at strain 0.4) in the [100] direction and the max tensile stress(92.1 GPa at strain 0.275)in the[110]directions. When the strain exceeds 0.4 in the[100]and[001]directions,the tensile tI16-B8C6O2will be destroyed,while for the[110]direction,the largest strain is only about 0.275. For tensile in the[001]direction,the stress and strain increase synchronously,however for[100]and[110]directions,there is a phenomenon of rising first and then falling in the low strain region;then stress and strain increase synchronously in the high strain region. As shown in Table S6,the B—O bonds continue to get longer until the bonds break in any direction we studied,the B—C bonds also get longer monotonous in[001]direction. However,for the[100]and[110]directions,the change of B—C bond is quite complex. In general, with the increase of strain, some B—C bonds continue to get longer, and some B—C bonds experience the process of combination of longer and shorter. During stretching, the bond length changes due to the presence of strain,and the population of the bonds also changes until they break and then the structure will collapse.The calculated strain-stress relations reveal the anisotropy of tI16-B8C6O2.

    Fig.9. The simulation of tensile stress—strain relationships along the specific direction for tI16-B8C6O2.

    As is well-known,bulk modulusB,shear modulusG,and Young modulusEare important parameters to reflect the mechanical properties about the resistance deformation of condensed materials. Generally,Bindicates the ability to resist volume deformation by loading pressure,Grepresents the ability to resist deformation upon shear stress.[33]BandGcanbe acquired by an independent elastic parameter, thenEcan be obtained byBandG. The value ofErepresents the rigidity of the condensed material; the greaterEis,the more difficult it is to deform.

    Here the mechanical properties, asB,G, andEof these six non-IED B—C—O phases and their relationship with pressure,are studied in detail. As shown in Fig.10,tI16-B8C6O2has larger elasticity moduli than mP13-B8C5O3at AP. The elasticity moduli of tI16-B8C6O2changes constantly and regularly,the increments ofB,G,andEare 256.6 GPa,139.0 GPa,and 353.9 GPa, respectively. According to data fitting, the modulusBof tI16-B8C6O2and pressurePsatisfies a linear relationship asB=289.54+3.171P. However, as the structure of mP13-B8C5O3changes under pressure, the moduli vary with the pressure also fall in two portions divided by the mutations during 30 GPa to 40 GPa. The moduli’ increments of mP13-B8C5O3are 324.8 GPa (B), 221.7 GPa (G),and 542.0 GPa(E). At 80 GPa HP,the moduli asEandGare an advantage rather than that of tI16-B8C6O2.

    Fig. 10. The relationship of pressure and mechanical properties including bulk modulus B, shear modulus G and Young’s modulus E of tI16-B8C6O2 and mP13-B8C5O3. The unit of modulus B,G,and E is unified as GPa.

    For the relationship between elasticity moduli and pressure of tI12-B6C4O2and aP13-B6C2O5are plotted in Fig. S6(a), and for mC20-B2CO2and tI18-B4CO4are presented in Fig. S6(b). For modulusEof tI12-B6C4O2go through first rapidly rise then slightly descend, and the max value is 653.5 GPa at pressure of 90 GPa.Except theEof tI12-B6C4O2,all the elasticity moduli of these four B—C—O phases increase monotonically. Compared with the reported bulk modulus of tI18-B4CO4in the range of 0—50 GPa,[15]the similarity of the data and the consistency of the variation indicate the reliability of our calculation. Among all the relationships between moduli and pressure,the linear relationships between moduliBof these four non-IED B—C—O phases and pressurePare more obvious, therein, tI12-B6C4O2:B= 283.90+3.245P;aP13-B6C2O5:B=236.42+3.386P;mC20-B2CO2:B=296.16+3.407P;tI18-B4CO4:B=247.85+3.190P.

    The hardness is extensively adopted to evaluate the mechanical properties of condensed matter. Here, the Vickers hardness(HV)of six non-IED B—C—O phases were calculated based on Eq. (14).[34]As reflect the brittleness or toughness of materials,Poisson’s ratioμis also an important mechanical property of materials and can be calculated by Eq. (15). We also studied in detail the effect of pressure on their hardness and Poisson’s ratioμ.

    As listed in Table S7,tI16-B8C6O2has high hardness as 21.774 GPa at AP,which is higher than that of mP16-B8C5O3(17.034 GPa). With the increase of pressure, the hardness of tI16-B8C6O2generally increases first and then decreases,and the max value ofHVis 31.808 GPa at 50 GPa. For mP16-B8C5O3,the hardness continues to increase with the pressure increasing. With the increase of pressure, the Poisson’s ratioμof tI16-B8C6O2first decreases and then increases, the min value is 0.223. For the Poisson’s ratioμof mP16-B8C5O3,except for 40 GPa, it also experienced a descent followed by an ascent during the pressure increasing. In the whole pressure range we studied, Poisson’s ratioμof tI16-B8C6O2and mP16-B8C5O3did not exceed the critical value of 0.333,indicating that they both are hard and brittle materials.

    For the other four non-IED B—C—O phases, the effects of pressure on Vickers hardnessHVand Poisson’s ratioμare shown in Fig. 11. With loading pressure,HVof mC20-B2CO2and tI18-B4CO4both decrease from 37.919 GPa to 23.809 GPa, and 36.226 GPa to 28.986 GPa, respectively.Based on Eq. (14) and the reported bulk and shear modulus data of mC20-B2CO2under pressure,[16]the effects of pressure on Vickers hardnessHVof mC20-B2CO2can be further developed, and the same downward trend confirms the accuracy of our calculations. The Poisson’s ratioμof mC20-B2CO2and tI18-B4CO4continues to increase.Compared with tI18-B4CO4,mC20-B2CO2has greater variation of Poisson’s ratioμand hardnessHV.

    For aP13-B6C2O5,except for the oddity point at 10 GPa,the general trend is either regularly down(HV)or up(μ)with the process of increasing pressure. Based on Eq. (14), it is found that aP13-B6C2O5decreases first and then rises within AP to 20 GPa,[16]this is consistent with the conclusion of our research as the same region,which also proves the complexity of the non-IED B—C—O system and the reliability of the calculation. For tI12-B6C4O2, with the pressure increasing, its hardness first increases from 23.121 GPa(AP)to 25.107 GPa(10 GPa) and then decreases to 16.018 GPa (100 GPa), its Poisson’s ratioμalmost remaining the same at first,and then increase to 0.323, which is closed to the critical value 0.333,indicating that with the increase of pressure,the brittleness of tI12-B6C4O2is weakened and the toughness is strengthened.Our study not only shows their properties at AP, but also explains the law of their mechanical properties effected by pressure.

    Fig.11. The relationship of pressure and(a)Vickers hardness HV, (b)Poisson’s ratio μ of tI12-B6C4O2, aP13-B6C2O5, mC20-B2CO2, and tI18-B4CO4.

    4. Conclusion

    Through the candidate structures generated by evolutionary simulation methods and structural stability analysis by first-principles, two novel non-IED B—C—O phases (tI16-B8C6O2and mP16-B8C5O3) have been unmasked. The research of the phonon scattering spectra and the independent elastic constants under AP and HP also proves the stability of these non-IED B—C—O phases. Respected to the common compounds containing carbon, boron and oxygen, the formation enthalpies and the relationship with pressure of all known(tI12-B6C4O2,aP13-B6C2O5,mC20-B2CO2and tI18-B4CO4) and two newly discovered non-IED B—C—O phases have been systematically studied, and the negative formation enthalpies at specific pressures indicate that HP technology performed in DAC is an important means of synthesizing all non-IED B—C—O phases. Especially, aP13-B6C2O5, mC20-B2CO2and tI18-B4CO4have the possibility of HP synthesis in LVP.

    Both tI16-B8C6O2and tI12-B6C4O2possess electrical conductivity.mP16-B8C5O3is an indirect semiconductor with a small bandgap as 0.530 eV. For the aP13-B6C2O5, mC20-B2CO2and tI18-B4CO4are all indirect semiconductors with the gap of 5.643 eV,6.113 eV,and 7.105 eV,respectively. The relationship of band gap values and the pressure of these six non-IED B—C—O phases are researched. Both tI16-B8C6O2and tI12-B6C4O2maintain electrical conductivity. For mC20-B2CO2and tI18-B4CO4,their gap values both have very slight change during the entail pressure range. mP16-B8C5O3and aP13-B6C2O5varies gaps in multiple stages. The stress-strain simulation found that the max stresses for tI16-B8C6O2are 141.9 GPa in the[001]direction with strain 0.4,123.1 GPa in the[100]direction at strain 0.4,and 92.1 GPa in[110]directions at strain 0.275. When the strain/strain exceed this, the tI16-B8C6O2will be destroyed. Studies on their mechanical properties show that they all possess elasticity moduli and hard character. And pressure has an obvious effect on their mechanical properties,therein toughness of tI12-B6C4O2,aP13-B6C2O5, mC20-B2CO2and tI18-B4CO4all increases and hardness of mP16-B8C5O3continue to strengthen during the compression. The enriched mechanical and electronic properties of non IED B—C—O compounds are very likely to attract extensive attention in scientific research and industrial applications.

    Acknowledgments

    Project supported by the National Natural Science Foundation of China (Grant No. 12064013), the Natural Science Foundation of Jiangxi Province, China (Grant No. 20202BAB214010), the Open Funds of the State Key Laboratory of Metastable Materials Science and Technology(Grant No.201906),Ganzhou Science and Technology Project(Grant No.202060), and the Program of Qingjiang Excellent Young Talents,Jiangxi University of Science and Technology.

    猜你喜歡
    劉超
    老吳的拉面館
    Removal of GaN film over AlGaN with inductively coupled BCl3/Ar atomic layer etch
    征程萬里,初心如一
    雷鋒(2021年12期)2021-04-12 00:57:22
    《站立式起跑》教學(xué)設(shè)計(jì)
    小數(shù)與小數(shù)點(diǎn)
    小數(shù)與小數(shù)點(diǎn)
    Design and experiment of the centrifugal pump impellers with twisted inlet vice blades *
    A simple method for estimating bed shear stress in smooth and vegetated compound channels*
    A randomized, controlled trial of Shutangluo fang in treating painful diabetic peripheral neuropathy
    Support Effects on Thiophene Hydrodesulfurization over Co-Mo-Ni/Al2O3and Co-Mo-Ni/TiO2-Al2O3Catalysts*
    www国产在线视频色| 极品教师在线免费播放| 国产一区二区在线av高清观看| 黄频高清免费视频| 琪琪午夜伦伦电影理论片6080| 国产高清激情床上av| 一级毛片精品| 电影成人av| 国产成人影院久久av| 大码成人一级视频| 欧美性长视频在线观看| 欧美在线黄色| 交换朋友夫妻互换小说| 一级片'在线观看视频| 亚洲欧美一区二区三区久久| 国产亚洲精品第一综合不卡| 国产一区二区三区在线臀色熟女 | 激情视频va一区二区三区| 在线观看一区二区三区| 国产成人精品久久二区二区91| 99久久精品国产亚洲精品| 人妻丰满熟妇av一区二区三区| 国产精品乱码一区二三区的特点 | 女人被躁到高潮嗷嗷叫费观| 亚洲国产精品合色在线| 无限看片的www在线观看| 成人永久免费在线观看视频| 日本三级黄在线观看| 欧美日韩亚洲国产一区二区在线观看| 伦理电影免费视频| 黄色怎么调成土黄色| 91在线观看av| 中文字幕av电影在线播放| 一级a爱视频在线免费观看| 国产在线观看jvid| 天天添夜夜摸| 性少妇av在线| 老司机深夜福利视频在线观看| 国产三级在线视频| 午夜亚洲福利在线播放| 成人国语在线视频| 激情在线观看视频在线高清| 国产av一区二区精品久久| 日本免费a在线| 美女 人体艺术 gogo| av国产精品久久久久影院| 妹子高潮喷水视频| 国产精品九九99| 国产熟女午夜一区二区三区| 男男h啪啪无遮挡| 日本免费一区二区三区高清不卡 | 黄片大片在线免费观看| 夜夜躁狠狠躁天天躁| 亚洲中文字幕日韩| 国产一区二区三区综合在线观看| 国产精品秋霞免费鲁丝片| 免费一级毛片在线播放高清视频 | 欧美激情极品国产一区二区三区| 三级毛片av免费| 久久国产亚洲av麻豆专区| 18禁观看日本| 视频区欧美日本亚洲| 日本一区二区免费在线视频| 午夜91福利影院| 国产亚洲精品第一综合不卡| av网站免费在线观看视频| 老司机福利观看| 亚洲午夜理论影院| 丁香欧美五月| 中国美女看黄片| 免费av中文字幕在线| 亚洲精品一卡2卡三卡4卡5卡| 一a级毛片在线观看| 久久久久久久精品吃奶| 欧美另类亚洲清纯唯美| av国产精品久久久久影院| 国产区一区二久久| 日韩有码中文字幕| 日本三级黄在线观看| 成人亚洲精品av一区二区 | 在线永久观看黄色视频| 19禁男女啪啪无遮挡网站| 18禁美女被吸乳视频| 国产精品成人在线| 成人国语在线视频| 少妇 在线观看| 久久天堂一区二区三区四区| 亚洲一卡2卡3卡4卡5卡精品中文| 在线观看舔阴道视频| 国产男靠女视频免费网站| 99精国产麻豆久久婷婷| 亚洲精品国产一区二区精华液| 色老头精品视频在线观看| av在线播放免费不卡| 夜夜爽天天搞| 国产精品一区二区免费欧美| 日韩精品青青久久久久久| 国产成人精品在线电影| 婷婷丁香在线五月| 久久精品91蜜桃| 久久久水蜜桃国产精品网| 国产极品粉嫩免费观看在线| 高清欧美精品videossex| 欧美av亚洲av综合av国产av| 亚洲国产欧美一区二区综合| 国产成人av激情在线播放| 精品国产亚洲在线| 精品久久久久久久毛片微露脸| 日本vs欧美在线观看视频| 好看av亚洲va欧美ⅴa在| 欧美av亚洲av综合av国产av| 中国美女看黄片| 国产人伦9x9x在线观看| 操出白浆在线播放| 激情视频va一区二区三区| 韩国精品一区二区三区| 国产成人精品在线电影| 99热国产这里只有精品6| 脱女人内裤的视频| 亚洲三区欧美一区| 久久天堂一区二区三区四区| 精品免费久久久久久久清纯| 精品熟女少妇八av免费久了| 国产aⅴ精品一区二区三区波| 欧美不卡视频在线免费观看 | 亚洲精品国产区一区二| 国产深夜福利视频在线观看| 少妇 在线观看| 女人被躁到高潮嗷嗷叫费观| 日韩欧美一区二区三区在线观看| 国产一区二区激情短视频| 满18在线观看网站| 一边摸一边抽搐一进一小说| 日韩大码丰满熟妇| 午夜精品久久久久久毛片777| 精品欧美一区二区三区在线| 我的亚洲天堂| 国产1区2区3区精品| 成年版毛片免费区| 日本免费一区二区三区高清不卡 | 欧美老熟妇乱子伦牲交| 90打野战视频偷拍视频| 中亚洲国语对白在线视频| 久久午夜亚洲精品久久| 免费久久久久久久精品成人欧美视频| 亚洲午夜精品一区,二区,三区| 1024香蕉在线观看| 精品久久久久久电影网| 又紧又爽又黄一区二区| 国产成人精品在线电影| 正在播放国产对白刺激| 国产男靠女视频免费网站| 亚洲av熟女| 国产精品一区二区在线不卡| 亚洲精品在线观看二区| 又大又爽又粗| 国产亚洲精品一区二区www| 久久狼人影院| 午夜老司机福利片| 久热爱精品视频在线9| 国产精品二区激情视频| 老司机福利观看| 丰满饥渴人妻一区二区三| 成人亚洲精品一区在线观看| 夫妻午夜视频| 日本五十路高清| cao死你这个sao货| 免费一级毛片在线播放高清视频 | 欧美激情高清一区二区三区| 国产日韩一区二区三区精品不卡| 91成人精品电影| 日韩人妻精品一区2区三区| 午夜免费鲁丝| 一级a爱片免费观看的视频| 国产精品久久久久久人妻精品电影| 伦理电影免费视频| 国产野战对白在线观看| 男人的好看免费观看在线视频 | 久久久久久久精品吃奶| 热re99久久国产66热| 天天影视国产精品| 黄色片一级片一级黄色片| 国产精品乱码一区二三区的特点 | 久9热在线精品视频| 国产成+人综合+亚洲专区| 精品人妻在线不人妻| 欧美日韩一级在线毛片| 亚洲专区国产一区二区| 天天添夜夜摸| 看免费av毛片| 久热这里只有精品99| 午夜精品久久久久久毛片777| 亚洲自偷自拍图片 自拍| 欧美午夜高清在线| 女人被狂操c到高潮| 女同久久另类99精品国产91| 中文字幕人妻丝袜制服| 淫妇啪啪啪对白视频| 久久久久亚洲av毛片大全| 啪啪无遮挡十八禁网站| 咕卡用的链子| 长腿黑丝高跟| 黄色丝袜av网址大全| 欧美一区二区精品小视频在线| 女性生殖器流出的白浆| 在线观看免费午夜福利视频| 免费看a级黄色片| 国产91精品成人一区二区三区| 免费高清视频大片| 一个人免费在线观看的高清视频| 淫秽高清视频在线观看| 国产97色在线日韩免费| 美女福利国产在线| 国产视频一区二区在线看| av视频免费观看在线观看| 大型av网站在线播放| 国产97色在线日韩免费| 韩国av一区二区三区四区| 午夜精品国产一区二区电影| 婷婷六月久久综合丁香| 国产伦人伦偷精品视频| 亚洲成国产人片在线观看| 9191精品国产免费久久| 一进一出抽搐gif免费好疼 | 欧美日韩一级在线毛片| 人人澡人人妻人| 动漫黄色视频在线观看| 琪琪午夜伦伦电影理论片6080| 免费久久久久久久精品成人欧美视频| 村上凉子中文字幕在线| 99在线视频只有这里精品首页| 岛国在线观看网站| 国产一区二区三区在线臀色熟女 | 国产三级黄色录像| 男女高潮啪啪啪动态图| 久久国产乱子伦精品免费另类| 91国产中文字幕| 香蕉久久夜色| av有码第一页| 久久人妻av系列| 国产亚洲精品第一综合不卡| 国产精品美女特级片免费视频播放器 | 午夜91福利影院| a级片在线免费高清观看视频| 性少妇av在线| 好男人电影高清在线观看| 最近最新免费中文字幕在线| 欧美乱色亚洲激情| 欧美日韩中文字幕国产精品一区二区三区 | 色婷婷av一区二区三区视频| 悠悠久久av| 免费搜索国产男女视频| 中文字幕精品免费在线观看视频| 波多野结衣一区麻豆| 女性被躁到高潮视频| 五月开心婷婷网| 91大片在线观看| 一级a爱片免费观看的视频| 欧美日韩福利视频一区二区| 女同久久另类99精品国产91| 啪啪无遮挡十八禁网站| 在线十欧美十亚洲十日本专区| 亚洲精品国产精品久久久不卡| 国产精品久久电影中文字幕| 亚洲专区中文字幕在线| 成年人黄色毛片网站| 曰老女人黄片| 性色av乱码一区二区三区2| 亚洲欧美日韩高清在线视频| 亚洲视频免费观看视频| 91麻豆av在线| 国产精品久久久av美女十八| 国产区一区二久久| 久久欧美精品欧美久久欧美| 国产欧美日韩综合在线一区二区| 久久婷婷成人综合色麻豆| av网站免费在线观看视频| 在线观看舔阴道视频| 长腿黑丝高跟| 亚洲欧美日韩高清在线视频| 99热只有精品国产| 欧美另类亚洲清纯唯美| 美女午夜性视频免费| 热99re8久久精品国产| 激情视频va一区二区三区| 免费人成视频x8x8入口观看| 99国产精品一区二区三区| 国产精品九九99| 两个人看的免费小视频| 久久亚洲精品不卡| 欧美一级毛片孕妇| 国产在线观看jvid| 亚洲一区高清亚洲精品| 新久久久久国产一级毛片| 久久久国产成人免费| 黄色成人免费大全| 高清黄色对白视频在线免费看| 热re99久久国产66热| 亚洲精品美女久久久久99蜜臀| 日本三级黄在线观看| 丁香欧美五月| 99热国产这里只有精品6| 老司机深夜福利视频在线观看| 国产1区2区3区精品| 人妻久久中文字幕网| 国产亚洲精品一区二区www| 波多野结衣av一区二区av| 亚洲av熟女| √禁漫天堂资源中文www| 亚洲精品一卡2卡三卡4卡5卡| 人人妻人人爽人人添夜夜欢视频| 一区在线观看完整版| 99国产精品99久久久久| 国产三级在线视频| 久久精品亚洲av国产电影网| 久久久久久久久中文| 日韩 欧美 亚洲 中文字幕| 制服诱惑二区| 亚洲午夜理论影院| 村上凉子中文字幕在线| av网站免费在线观看视频| 国产高清国产精品国产三级| 女同久久另类99精品国产91| 深夜精品福利| 很黄的视频免费| 一区二区三区国产精品乱码| 色在线成人网| 最近最新免费中文字幕在线| 怎么达到女性高潮| 久久人人精品亚洲av| cao死你这个sao货| 亚洲在线自拍视频| 日韩欧美国产一区二区入口| 啪啪无遮挡十八禁网站| av天堂久久9| 久久欧美精品欧美久久欧美| 日本vs欧美在线观看视频| 成人国产一区最新在线观看| 黄色片一级片一级黄色片| 午夜视频精品福利| 欧美日韩一级在线毛片| 欧美+亚洲+日韩+国产| 在线观看免费日韩欧美大片| 国产成人精品久久二区二区免费| 三上悠亚av全集在线观看| 欧美日韩亚洲综合一区二区三区_| 精品日产1卡2卡| 欧美成人午夜精品| 久久中文字幕人妻熟女| 天堂影院成人在线观看| 12—13女人毛片做爰片一| 满18在线观看网站| 黄色a级毛片大全视频| 身体一侧抽搐| 欧美成人性av电影在线观看| 久久青草综合色| 亚洲成人精品中文字幕电影 | 麻豆av在线久日| 免费一级毛片在线播放高清视频 | 久久亚洲精品不卡| 99香蕉大伊视频| 丝袜美足系列| 在线观看一区二区三区| bbb黄色大片| 欧美成人午夜精品| 男女之事视频高清在线观看| 99久久精品国产亚洲精品| 亚洲国产中文字幕在线视频| 免费看a级黄色片| 欧美成狂野欧美在线观看| 91九色精品人成在线观看| 亚洲视频免费观看视频| 麻豆国产av国片精品| 在线观看舔阴道视频| 久久久国产欧美日韩av| 性色av乱码一区二区三区2| 国产单亲对白刺激| 色老头精品视频在线观看| 国产黄a三级三级三级人| 美女扒开内裤让男人捅视频| 80岁老熟妇乱子伦牲交| 亚洲男人的天堂狠狠| 午夜视频精品福利| 午夜免费激情av| 国产熟女xx| 美女高潮喷水抽搐中文字幕| 国产97色在线日韩免费| ponron亚洲| 国产精品久久视频播放| 午夜福利影视在线免费观看| 十八禁网站免费在线| 精品欧美一区二区三区在线| 97碰自拍视频| 男人的好看免费观看在线视频 | 日韩免费av在线播放| 国产精品久久久久久人妻精品电影| 黄色毛片三级朝国网站| 亚洲精品一二三| 夜夜夜夜夜久久久久| 国产亚洲av高清不卡| 日韩欧美免费精品| 一区二区三区国产精品乱码| 欧美黑人欧美精品刺激| svipshipincom国产片| 午夜免费成人在线视频| 欧美色视频一区免费| 9热在线视频观看99| 亚洲av电影在线进入| 黑人欧美特级aaaaaa片| svipshipincom国产片| 免费看a级黄色片| 日本vs欧美在线观看视频| 身体一侧抽搐| 丝袜人妻中文字幕| 两性夫妻黄色片| 黄色丝袜av网址大全| 91字幕亚洲| 99香蕉大伊视频| 少妇 在线观看| 黄频高清免费视频| 91在线观看av| 19禁男女啪啪无遮挡网站| 人人澡人人妻人| 午夜福利欧美成人| 在线观看午夜福利视频| 亚洲专区国产一区二区| 最近最新中文字幕大全免费视频| 韩国av一区二区三区四区| 精品久久久久久久毛片微露脸| 丰满迷人的少妇在线观看| 日韩欧美一区视频在线观看| 亚洲欧美精品综合久久99| 久久久久久免费高清国产稀缺| 母亲3免费完整高清在线观看| 久久精品亚洲熟妇少妇任你| 国产精品爽爽va在线观看网站 | 精品久久久久久,| 日本撒尿小便嘘嘘汇集6| 黄色 视频免费看| 久9热在线精品视频| 丝袜在线中文字幕| 欧美成狂野欧美在线观看| 日本黄色视频三级网站网址| 午夜免费激情av| 成人特级黄色片久久久久久久| 午夜福利在线观看吧| 一二三四在线观看免费中文在| 精品福利永久在线观看| 色婷婷久久久亚洲欧美| 黑丝袜美女国产一区| 在线国产一区二区在线| 国产一区二区三区综合在线观看| 国产在线精品亚洲第一网站| x7x7x7水蜜桃| 亚洲精品成人av观看孕妇| 久久人妻熟女aⅴ| 在线国产一区二区在线| 国产精品一区二区三区四区久久 | 婷婷精品国产亚洲av在线| 国产熟女午夜一区二区三区| 一级a爱视频在线免费观看| 午夜亚洲福利在线播放| 久久精品国产综合久久久| 19禁男女啪啪无遮挡网站| 国产成人av激情在线播放| 美女高潮喷水抽搐中文字幕| 国产在线精品亚洲第一网站| 80岁老熟妇乱子伦牲交| 多毛熟女@视频| 国产欧美日韩一区二区三| 国产亚洲欧美精品永久| 日韩 欧美 亚洲 中文字幕| 精品久久蜜臀av无| 国产在线观看jvid| 中文字幕色久视频| 搡老岳熟女国产| 久久草成人影院| 在线视频色国产色| 成人国语在线视频| 久久久国产欧美日韩av| 久久久久久亚洲精品国产蜜桃av| 无遮挡黄片免费观看| 午夜精品久久久久久毛片777| 老熟妇仑乱视频hdxx| 18禁黄网站禁片午夜丰满| 精品电影一区二区在线| 中文字幕人妻熟女乱码| 久久婷婷成人综合色麻豆| 国产成人影院久久av| 亚洲成人免费电影在线观看| 国产精品久久久久久人妻精品电影| 岛国在线观看网站| 满18在线观看网站| 亚洲色图av天堂| 欧美成人性av电影在线观看| 亚洲一区中文字幕在线| 男女之事视频高清在线观看| 亚洲国产精品合色在线| 久久精品91无色码中文字幕| 欧美精品啪啪一区二区三区| 亚洲国产精品一区二区三区在线| 9热在线视频观看99| 一级片免费观看大全| 国产乱人伦免费视频| 亚洲一区二区三区欧美精品| 老司机亚洲免费影院| 18禁美女被吸乳视频| 一本综合久久免费| 在线观看一区二区三区激情| 亚洲五月婷婷丁香| 在线国产一区二区在线| 人人妻人人爽人人添夜夜欢视频| 在线av久久热| 妹子高潮喷水视频| 久久人妻福利社区极品人妻图片| 国产精品永久免费网站| 12—13女人毛片做爰片一| 交换朋友夫妻互换小说| 欧美日韩国产mv在线观看视频| avwww免费| 日韩有码中文字幕| 亚洲中文字幕日韩| 欧美成人性av电影在线观看| 在线观看免费高清a一片| 久久精品国产清高在天天线| 亚洲人成电影免费在线| 欧美不卡视频在线免费观看 | 亚洲精品一卡2卡三卡4卡5卡| av有码第一页| 成人18禁在线播放| 免费看a级黄色片| 国产色视频综合| 国产无遮挡羞羞视频在线观看| 波多野结衣一区麻豆| 亚洲专区中文字幕在线| 香蕉国产在线看| 妹子高潮喷水视频| 欧美色视频一区免费| 悠悠久久av| x7x7x7水蜜桃| 在线观看一区二区三区激情| √禁漫天堂资源中文www| 久久中文看片网| xxxhd国产人妻xxx| 国产成人精品在线电影| 韩国av一区二区三区四区| 香蕉丝袜av| a级毛片黄视频| 又大又爽又粗| 超碰97精品在线观看| 一区在线观看完整版| 日韩免费av在线播放| 久久这里只有精品19| 婷婷精品国产亚洲av在线| 成人国语在线视频| 美女午夜性视频免费| 日韩欧美在线二视频| 精品国产美女av久久久久小说| 两个人看的免费小视频| 亚洲精品美女久久av网站| 久久国产亚洲av麻豆专区| 嫩草影院精品99| 欧美 亚洲 国产 日韩一| 国产精品久久久久久人妻精品电影| 变态另类成人亚洲欧美熟女 | 欧美日韩视频精品一区| 日本免费一区二区三区高清不卡 | 夜夜躁狠狠躁天天躁| 午夜福利在线免费观看网站| 免费av毛片视频| 久久青草综合色| 日本黄色日本黄色录像| 国产成人精品久久二区二区免费| 天堂中文最新版在线下载| 少妇 在线观看| 免费在线观看黄色视频的| 黄色 视频免费看| 成人国语在线视频| 另类亚洲欧美激情| 国产99久久九九免费精品| 嫩草影院精品99| 国产精品98久久久久久宅男小说| 国产精品成人在线| 亚洲专区中文字幕在线| 两人在一起打扑克的视频| 村上凉子中文字幕在线| 一级片'在线观看视频| 在线观看免费视频网站a站| 免费久久久久久久精品成人欧美视频| 丰满的人妻完整版| 一级毛片女人18水好多| 成在线人永久免费视频| 激情在线观看视频在线高清| 十八禁人妻一区二区| 一进一出抽搐gif免费好疼 | 精品人妻1区二区| 精品午夜福利视频在线观看一区| 久久精品国产清高在天天线| 日日摸夜夜添夜夜添小说| 51午夜福利影视在线观看| 在线观看免费视频日本深夜| 一级作爱视频免费观看| 亚洲美女黄片视频| 最近最新中文字幕大全电影3 | 99精品在免费线老司机午夜| 欧美日韩亚洲综合一区二区三区_| 亚洲五月色婷婷综合| 日韩有码中文字幕| 在线观看免费视频网站a站| 国产99久久九九免费精品| 亚洲国产精品一区二区三区在线| 动漫黄色视频在线观看| 国产97色在线日韩免费| 另类亚洲欧美激情| 99国产精品99久久久久|