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

    Chair-like Nin AlN3 with high-energy density

    2023-12-15 11:48:26ShiTaiGuo郭世泰ZhenZhenXu徐真真YanLeiGeng耿延雷QiRui芮琦DianChenDu杜殿臣JianFuLi李建福andXiaoLiWang王曉麗
    Chinese Physics B 2023年12期

    Shi-Tai Guo(郭世泰), Zhen-Zhen Xu(徐真真), Yan-Lei Geng(耿延雷), Qi Rui(芮琦),Dian-Chen Du(杜殿臣), Jian-Fu Li(李建福),?, and Xiao-Li Wang(王曉麗),?

    1School of Physics and Electronic Information,Yantai University,Yantai 264005,China

    2School of Physics and Electronic Engineering,Linyi University,Linyi 276005,China

    Keywords: high pressure,aluminum nitride,electronic property,high-energy-density materials

    1.Introduction

    Recently, all-nitrogen salts have drawn much attention in the high-energy-density materials (HEDMs) field due to their potential as ideal candidates in propellant or explosive applications.[1]Generally, the high energy of these compounds is primarily attributed to the substantial disparity in bond energy between N-N or N=N bonds in nonmolecular nitrogen (160/418 kJ/mol) and the N≡N triple bond in N2molecule (954 kJ/mol).In the past decade, several nitrogenrich metal compounds with high energy density have been reported, including Be-N,[2-6]Mg-N,[7,8]Al-N,[9-11]S-N,[12]Fe-N,[13-15]Cu-N,[16,17]Zn-N,[18-20]Se-N,[21]Hg-N,[22]Ta-N,[23]W-N,[24]Y-N,[25]Ne-N,[26]and others.Theoretical calculations and experimental synthesis have demonstrated the existence of diverse all-nitrogen ions or clusters ranging from N3to N13.[27,28]However,the successful preparation of such compounds remains rare.In 2001, Christe’s work demonstrated the stability ofand N+5Sb2F-11salts up to 70?C.[29]More recently, Lu’s group achieved a breakthrough in the field by experimentally synthesizing cyclo-N-5.[30,31]The achievement opens up new possibilities for further advancements in the field,particularly in the areas of explosives and propellants.As we all know, one notable example of a high-energy-density compound is cubic gauche nitrogen(cg-N),which is solely connected by N-N bonds and possesses a high energy density of 9.7 kJ/g-twice that of traditional energetic compounds like TNT and HMX.However,the synthesis of cg-N requires extremely high pressures and temperatures(2000 K,110 GPa)conditions.[32]Finding highenergy-density single-bond poly-nitrogen compounds that can remain stable at lower pressures is a challenging task that requires a specific solution approach.Fortunately, the discovery of t-N,[33]a dynamically and mechanically stable compound at ambient pressure with an estimated energy density of approximately 11.31 kJ/g, encourages further exploration of high-energy-density materials within poly-nitrogen compounds under ambient conditions.This successful theoretical prediction was made in 2018, prompting researchers to continue the search for other potential candidates with similar properties.More recently, it has been observed that cg-N can exhibit stability in an acidic condition or when the surface is saturated with less electronegative adsorbates.[34]

    Numerous planar or quasiplanar cyclo-N6anions have been predicted in Li, Mg, Ca, and Ba nitrides under high pressure.[7,35-37]However, they have not paid close attention to the microstructural characteristics of these cyclo-N6anions.It’s worth noting that an in-plane distortion may occur in the cyclo-N6sub-lattices when the neutralπ-aromatic system transforms into a charged anion.For instance, chair-like N6anions have been predicted in h-WN6,TeN6,and GdN6,[38-40]where strong chemical bonding and charge transfers are believed to contribute to the structural stability.The synthesis of WN6, as mentioned in Ref.[24], has led to an increased interest in the study of N6anions.However, it is important to note that the content of nitrogen in these compounds will be decreased, resulting in relatively lower energy density due to the larger molar mass of tungsten(W)and tellurium(Te).In a separate study conducted by Liuet al.,[41]they reported the existence of GaN5and GaN6compounds with high nitrogen contents.It is worth mentioning that gallium(Ga)and aluminum(Al) are elements belonging to the same group.Generally,aluminum has a small atomic radius, high valence electron density,and high electronegativity,which makes it capable of forming strong covalent bonds with nitrogen atoms.Meanwhile,the relatively lower molar mass of aluminum compared to Te, W, Gd, and Ga allows for higher nitrogen content in aluminum-based compounds.Moreover, Cuiet al.have predicted the existence of a notable monoclinic phase with polymeric nitrogen chains(P21/c-AlN3)for AlN3.[9]This phase is stable within a specific pressure range of 43-85 GPa,which indicates that aluminum and nitrogen can indeed form nitrogenrich compounds under high pressure.Based on these factors,it is reasonable to consider binary Al-N compounds as potential candidates for the search of N6anions.

    In the present study,theR-3m-AlN3phase has been successfully designed, featuring a chair-likeing with N-N single bonds.The phonon dispersion calculations and firstprinciples molecular dynamics (FPMD) simulations indicate that the compound exhibits good kinetics and thermodynamic stabilities.The calculated electronic properties and chemical bonding patterns ofR-3m-AlN3reveal that its high stability arises from a combination of Coulomb interactions and covalent bonds.These factors contribute to the compound’s structural integrity and overall stability.Furthermore,the study calculates the formation of enthalpy forR-3m-AlN3, suggesting that the compound can be synthesized by condensing AlN and N2under appropriate pressure.

    2.Computational details

    The structure prediction for the Al-N system was carried out using the particle swarm optimization(PSO)methodology implemented in the CALYPSO code.[42,43]This approach does not rely on any known structural information and has been successfully validated in various systems under high pressure,including binary and ternary systems.[44-50]For the structural optimization and electronic properties calculations, density functional theory (DFT) with the Perdew-Burke-Ernzerhof(PBE) exchange-correlation functional and generalized gradient approximation (GGA) was employed.[51]These calculations were performed utilizing the VASP code.[52]The allelectron projector augmented wave (PAW)[53]method was adopted to represent the ionic potentials, the valence electrons of Al and N atoms are 3s23p1and 2s22p3, respectively.The plane-wave basis set cutoff (700 eV) and appropriate Monkhorst-Packk-meshes[54](2π×0.03 ?A-1) were used to ensure all enthalpy calculations converged less than 1 meV/atom.In addition, phonon dispersion curves were computed using the supercell method implemented in the PHONOPY code[55]to determine the dynamical stability.FPMD simulations were performed in the NVT ensemble from 300 K to 700 K through the Nos′e-Hoover method[56,57]to assess the thermodynamic stability of the system at different temperatures.We calculated the elastic constants employing the strain-stress method[58]as adopted in the VASP code.Relative elastic property parameters were determined using the Voigt-Reuss-Hill approximations.[59]The chemical bonding mode of AlN3was executed by the solid state adaptive natural density partitioning(SSAdNDP)method.[60]Orbital calculations were carried out at the B3LYP/6-311+g(d) level of theory using the Gaussian 09 code.[61]Bader analysis[62]was utilized to determine the transfer of charge within the system.The crystal structure patterns were visualized through VESTA[63]software.

    3.Results and discussion

    3.1.Crystal structure

    We conducted a comprehensive search for all potential structures within the AlN3system across a pressure range from 0 to 100 GPa.Subsequently, we identified a crystal structure that features an N6ring and proceeded to optimize these structures.The total energy of structures containing N6rings was calculated under ambient conditions.The results indicated that the structure with theR-3mspace group exhibited the lowest energy among all the selected structures at 0 GPa, as shown in Fig.1(a).Furthermore,it was observed that structures with chair and boat N6rings generally had lower energies compared to structures with planar N6rings in AlN3.Simultaneously,the energy of structures with chair N6rings was generally lower than that with boat N6rings.This disparity arises since the bond angles between nitrogen atoms in chair N6rings tend to align more closely with the ideal angle (120?) in contrast to boat N6rings.This alignment reduction minimizes the electron cloud overlap and Coulomb repulsion, consequently elevating the energy level of the conformation.Additionally, a comparison was made between the energy ofP21/c-AlN3andR-3m-AlN3in the pressure range from 0 to 100 GPa(Fig.S1).The findings revealed that theP21/cstructure had lower energy compared to theR-3mstructure.However, recent highpressure experiments have encouraged us to continue the research onR-3m-AlN3.It is worth noting that the recently synthesized bp-N is not the lowest in energy,[64,65]and materials like silicon and ice can form different phases under varying pressure-temperature paths.[66,67]This suggests that the phases with higher energy can be formed by controlling the synthesis conditions.

    Table 1.Calculated structural parameter for R-3m-AlN3 structure at 0 GPa.

    Fig.1.Crystal structure of AlN3.(a)Energy of AlN3 with planar,chair,and boat N6 rings units.(b)and(c)The structure of R-3m-AlN3 phase.(d)Separate chaired N6 ring structure.(e)The structural unit of 4-fold bonded N.Light blue and gray spheres represent Al and N atoms,respectively.

    TheR-3m-AlN3phase(Figs.1(b)and 1(c))and the separate chaired N6ring structure(Fig.1(d))are illustrated.Under atmospheric conditions, the bond length and bond angle for the chaired N6ring structure are equal,measuring 1.486 ?A and 103.96?, respectively.In theR-3m-AlN3phase, six N atoms form an N6ring network, and each Al atom is bonded to six N atoms,forming an octahedral structure with a sixfold coordinated mode.Each N atom is driven by the sp3hybridization of two adjacent N atoms and Al atoms,forming four covalent bonds.The Al-N distances in the structure are 2.019 ?A and 2.149 ?A,respectively(Fig.1(e)).The detailed information on the structure is summarized in Table 1.

    3.2.Kinetic and thermal stability

    To assess the kinetics stability of theR-3m-AlN3structure at 0 GPa, phonon spectra were computed using the supercell method.In general, the presence of imaginary frequencies in the phonon spectra typically suggests the potential for vibration mode instability.Conversely, the absence of imaginary frequencies indicates a lack of vibration modes that could lead to structural instability.This also implies that the structure remains stable across a range of vibration modes.As shown in Figs.2(a)and 2(b),the absence of imaginary frequency in the over Brillouin zone demonstrates that the structure is kinetically stable.Based on the analysis of phonon density of states(PHDOS) for theR-3m-AlN3structure, it was observed that there are strongly coupled vibrations between the Al and N atoms,which contribute to the low-frequency vibration modes.Furthermore, the stretching mode of the N-N bond induces high-frequency vibrations.

    Fig.2.The dynamical and thermodynamic stabilities of the R-3m-AlN3 structure.(a) and (b) Phonon spectra and phonon DOS, respectively.(c) Fluctuations of the total energy of R-3m structure at the temperature of 300 K during FPMD simulations.The inset shows the structure pictures after 10 ps simulations at 300 K.

    Similarly, the thermal stability of the N6-6 anion in theR-3m-AlN3structure was assessed using a 2×2×2 supercell with 192 nitrogen atoms in the model.The evaluation was conducted through FPMD simulations,which were performed at a temperature of 300 K for 10 ps with a time step of 1 fs.The duration of 10 ps was chosen primarily because it allows for the observation of significant dynamic changes in the initial stage, such as structural changes and energy fluctuations,while also ensuring that the computational costs remain manageable.In Fig.2(c),the potential energy of the analog system is observed to fluctuate around 1330.2 eV.The inset shows the picture of the equilibrium structure at the end of 10 ps simulations, demonstrating that the N6-6 anion remains intact and is not broken, except for some thermal fluctuations.To further assess the thermal and kinetic stabilities, additional simulations were performed at higher temperatures.Figures S2 and S3 depict the temperature profiles at 500 K and 700 K,using a 1 fs time step.It is discovered that the N-N bond is broken at the 700 K temperature,leading to the disruption of the integrity of the N6ring.The plots of radial distribution function(RDF)for N-N separations are shown in Figs.S4 and S5.The first sharp peak in the RDF represents the nearest distance between the selected atom and its nearest neighboring atoms.It is shown that there is a single RDF peak at approximately 1.488 ?A corresponding to the N-N distance without significant change at 300 K and 500 K,which supports the conclusion that theR-3m-AlN3with N6-6 remains stable at ambient conditions and high temperatures.

    3.3.Mechanical stability

    It is widely known that the mechanical properties of any crystal are strongly influenced by its fundamental elastic constants(Cij).To verify the mechanical properties of theR-3m-AlN3structure at 0 GPa,the elastic constants were calculated using the strain-stress method.[58]TheR-3mphase requires the determination of six independent elastic constants,and the Born-Huang criterion[68]was employed to verify the mechanical stability of theR-3m-AlN3structure at 0 GPa.The mechanical stability criteria for the hexagonal phase are evaluated based on the following conditions:C44> 0,zC11>|C12|(C11+2C12)C33>2C213.Obviously,the elastic constants of the structure satisfy these criteria, indicating its mechanical stability.Through the elastic constants,the corresponding shear modulus (G), bulk modulus (B), Young’s modulus (E),and Poisson’s ratio(ν)can be deduced by the Voigte-Reusse-Hill approximation,[59]as shown in Table 2.

    Table 2.Calculated elastic constants Cij,bulk modulus B,shear modulus G,Young’s modulus E,and Poisson’s ratio v of R-3m-AlN3 structure at ambient pressure.

    It can be seen from Table 2, that theR-3m-AlN3phase exhibits a shear modulus of 133 GPa and a bulk modulus of 185 GPa.The Pugh criterion serves as an empirical approach used to predict the brittle or ductile behavior of crystalline materials.Importantly,it has been successfully applied to anticipate various material-related characteristics,and the outcomes of its predictions have gained recognition within the scientific community.According to Pugh’s criterion,[69]the ratio ofk=G/Bdetermines whether the materials behave in a brittle or ductile manner.Ifk>0.57, the material behaves in a brittle manner; otherwise, it is ductile.In the case ofR-3m-AlN3, the calculated value ofG/B(0.72) is larger than 0.57,similar to the diamond (1.10) and c-BN (0.86),[70]meaning the brittle materials.Additionally, the stiffness of a material can be characterized by Young’s modulus,[71]the high value(322 GPa)indicates that theR-3m-AlN3structure can be used as a potential hard material.Poisson’s ratio reflects the degree of covalent bonding in a material.[72]The value(0.21)is lower than 0.33, indicating a strong degree of covalent bonding.Taken together, these excellent mechanical parameters demonstrate the good mechanical stability of theR-3m-AlN3structure.

    3.4.Electronic properties

    The electronic band structure and projected density of states (PDOS) of theR-3m-AlN3phase at 0 GPa were calculated to explore its inner nature and formation mechanism,as shown in Fig.3.The calculations reveal that theR-3m-AlN3structure possesses semiconductor property with a large band gap of 4.14 eV.This suggests its potential for applications in electronic devices and optoelectronics.However,density functional calculations usually lead to a significant underestimation of the energy gap due to several primary factors,such as the approximation of exchange-correlation functionals,self-interaction errors,the size of the basis set,etc.Therefore,the actual band gap is expected to be greater than 4.14 eV.In addition,the electronic band structure shows a clear overlap between the Al 3p and N 2p orbitals,demonstrating the formation of a strong chemical bonding as depicted in Figs.3(b)and 3(c).To gain further insights into the behavior of the chemical bonds,Bader[62]charge analysis(shown in Table S1)was performed.The analysis revealed that there is a significant electron transfer (2.4e) between aluminum and nitrogen atoms.This information provides insights into the behavior of chemical bonds within theR-3m-AlN3structure.

    Fig.3.Electronic properties of the R-3m structure of AlN3.(a) Band structure.(b)and(c)Partial densities of states(PDOS)for Al and N,respectively.(d)and(e)The COHP between Al and N atoms and between N atoms,respectively.

    The crystal orbital Hamiltonian group (COHP) is a target of energy resolution to measure the overlap strength and crystal orbital overlap population, which helps to understand the contribution of atomic pairs to the structural stability determining the bonding information based on the COHP in the LOBSTER program.[73]Generally speaking, the positive and negative COHP values indicate the presence of antibonding and bonding states,respectively.The COHP plot of the Al-N pairs and N-N pairs are shown in Figs.3(d)and 3(e),respectively.It is obvious that bonding states are fully occupied,and the antibonding states are not occupied, which demonstrates the strong covalent bond feature between them.In addition,the ICOHP values quantify the strength of the antibonding interaction.In this case, the ICOHP values corresponding to Figs.3(d)and 3(e)are-4.49 and-9.18,respectively.These results reflect the strong covalent bonding characteristics between Al and N atoms in the Al-N pairs and the strong covalent bonding between neighboring N atoms in the N-N pairs.

    3.5.Chemical bonding pattern

    To illustrate the reason for the high stability of theR-3m-AlN3phase, the SSAdNDP method was used to analyze its chemical bonding patterns, relative results of the analysis are plotted in Fig.4.In the primitive cell,there are 36 2c-2eAl-Nσbonds(Fig.4(a))and 18 2c-2eN-Nσbonds(Fig.4(c)),the total number of electrons involved in these bonds is 108,which is in complete agreement with the sum of outer valence electrons in the AlN3unit cell.It is worth noting that one unit cell comprises six AlN3formulas,with Al and N elements contributing 3s23p1and 2s22p3, respectively.The reasonable occupation numbers of these bonds indicate that the computational results are reliable.Moreover, the bonding modes of Al-N and N-N in theR-3m-AlN3phase are illustrated in Figs.4(b)and 4(d), respectively.For theanion, each Al atom is bonded to six N atoms, forming a twisted ring structure.The distance between adjacent N atoms in the ring is measured to be 1.486 ?A,which is close to the typical N-N single bond distance of 1.45 ?A.[74]Furthermore,theR-3m-AlN3phase exhibits strong covalent bonding between N atoms and between Al and N atoms.This is evident due to the isosurface value exceeding 0.75,along with the presence of lone pair electrons, as revealed by the calculated electron localization functions(ELF)[75]shown in Figs.4(e)and S6.

    Fig.4.Schematic of the SSAdNDP chemical bonding pattern for AlN3 unit cell.(a) 36 2c-2e Al-N σ bonds, occupation number threshold:1.75|e|, the isosurface value=0.15.(c)18 2c-2e N-N σ bonds, occupation number threshold:1.75|e|,the isosurface value=0.2.(b)and(d)The bonding modes of Al-N and N-N,respectively.(e)The calculated ELF,the isosurface value=0.8.

    To demonstrate the bonding situation, Fig.5 displays a scheme of the molecular orbital diagram of the.The HOMO to HOMO-5 show that two pz-like lone pair electrons occur at each atom of the N6molecule.The HOMO-6 to HOMO-11 show additional py-like lone pair electrons on each N atom.HOMO-12 to HOMO-17 show strong s-pσbonds between N atoms.It is precise because the pzlone pair electrons ofenter the 3p orbital of the Al atom, leading to the formation of the coordinated covalent bond between the Al atom and the N atom.

    Fig.5.Occupied orbitals of the N6-6.The isosurface value is±0.05 and the blue and yellow colors denote negative and positive values,respectively.

    3.6.Stabilization mechanism

    Fig.6.Pressure-driven reaction mechanism between AlN with N2.(a)The relationship between atomic distance and external pressure.The blue line represents the sum (1.96 ?A) of Al-3p and N-2p covalent radii; the green line indicates the N-N distance for cg-N.(b) Calculated enthalpy difference (?H), internal energy (?E), and volume term (P?V) with pressure for the formation of AlN3.For the chemical reaction of AlN + N2 →AlN3, we defined ?H =HAlN3 -HAlN-HN2, ?E =EAlN3 -EAlN-EN2,P?V =PVAlN3-PVAlN-PVN2, where P and V represent the pressure and volume per formula unit,respectively.

    Subsequently, we investigated the formation mechanism of Al-N covalent bonds by analyzing the relationship between Al-N distances and the ever-increasing pressure, as depicted in Fig.6(a).Our calculations revealed that the distances between Al and N atoms gradually decrease as the pressure increases.Notably,the Al-N distance reaches the sum(1.96 ?A)of the covalent radii[76]of Al 3p and N 2p at a pressure of 20 GPa.Moreover, it is worth mentioning that throughout the pressure range studied, the Al-N distance remains larger than the N-N distance in the case of cg-N.This result provides strong evidence supporting the covalent bonding nature of the Al-N bonds.

    To illustrate the decomposition of AlN3into AlN and N2at high pressure,we plot the enthalpy difference(?H)between AlN3and AlN + N2.Additionally, we analyze the contributions to ?Hfrom the internal energy (?E) and volume term(P?V),as shown in Fig.6(b).It reveals that the decomposition of AlN3into AlN and N2occurs spontaneously at pressures higher than 46 GPa.Furthermore,the competition between the?EandP?Vterms is crucial in determining the actual chemical reaction.Although ?Eremains largely positive,theP?Vterm is negative due to the associated volume reduction, and this effect becomes more substantially reduced as the pressure increases.At 46 GPa, theP?Vterm dominates the competition, leading to a negative ?Hvalue and the formation of the stable AlN3compound.

    3.7.Explosive performance

    The energy density is a crucial parameter for evaluating the detonation properties of a material.Here, in the case of the predicted stable phase ofR-3m-AlN3,a significant amount of chemical energy is released when it dissociates into stable AlN solid and dinitrogen gas under ambient conditions.Considering the energy reduction of nitrogen gas concerning theαnitrogen phase(around 0.25 eV/atom),[80]the chemical energy densities (Ed) of denseR-3m-AlN3with a density of 3.45 g/cm3can be calculated to be approximately 5.04 kJ/g(Table 3).This calculation is based on the dissociation of AlN3(s)→AlN (s) + N2(g) and referencing to the energy of their most stable phase at 0 GPa.[81,82]Strikingly, the performance ofR-3m-AlN3in terms of energy density surpasses that of the typical explosive material TNT(4.3 kJ/g).[77]Furthermore,the volumetric energy densities(Ev)ofR-3m-AlN3are estimated to be 17.39 kJ/cm3,which is approximately 2.5 times that of TNT(~7.05 kJ/cm3)and 1.6 times that of HMX(~10.83 kJ/cm3).

    Table 3.Detonation properties of R-3m-AlN3 estimated by Kamlet-Jacobs empirical equations,compared to the experimental values of the known explosives of TNT and HMX.[77-79]

    3.8.Route to chair-like anion

    In the conducted research, an experimental synthesis of aluminum nitride (AlN) in theP63mcphase (wurtzite structure) has been successfully achieved and found to be stable within the pressure range of 0-20 GPa.[81]Here, to explore the feasibility of synthesizingR-3m-AlN3under high pressure,the formation enthalpy of AlN3relative to AlN and N2(Fig.7)has been calculated.Additionally, the phase transition process of AlN and N is shown in Figs.S7 and S8, respectively.The result shows that theR-3m-AlN3exhibits lower energy compared to AlN and N2under pressures exceeding 46 GPa.This suggests that the synthesis of N6-6 can be achieved by condensing AlN and N2at 46 GPa.Importantly,the structure remains mechanically stable under ambient conditions,which raises the possibility of recovery at ambient conditions.

    Fig.7.Enthalpy of formation of R-3m-AlN3 relative to AlN and N2.

    4.Conclusions

    Through swarm-intelligence structure searches,the presence of all-nitrogen anions N6-6 in the AlN3compound has been predicted.The calculation results of phonon dispersion and FPMD suggest that theR-3m-AlN3structure exhibits excellent kinetic and thermal stabilities at atmospheric pressure.Meanwhile,it demonstrates good mechanical stability at 0 GPa.The high stability ofR-3m-AlN3is attributed to the combined effects of Coulomb interactions and covalent bonds between Al3+and N6-6.Further analysis of the electronic properties reveals that the structure possesses a semiconductor character with a large band gap of 4.14 eV.In addition,the calculated formation enthalpies ofR-3m-AlN3indicate that it can be synthesized by condensing AlN and N2under appropriate pressure conditions.Upon decomposition, AlN3releases a substantial amount of energy, with a calculated energy release of 5.04 kJ/g.This indicates it possesses outstanding detonation properties, including high detonation velocity and pressure.For instance, the detonation velocity and pressure are approximately 1.9 and 5.3 times higher than that of TNT.These characteristics positionR-3m-AlN3as a promising candidate for high-energy-density material.The findings of this study are expected to facilitate the experimental synthesis of the AlN3phase and stimulate further research on all-nitrogen salts.

    Acknowledgements

    Project supported by the National Natural Science Foundation of China (Grant No.11974154), the Taishan Scholars Special Funding for Construction Projects, and the Natural Science Foundation of Shandong Province (Grant No.ZR2022MA004).

    久久久久久久国产电影| 18在线观看网站| 十分钟在线观看高清视频www| 国产精品成人在线| 亚洲精品自拍成人| av在线观看视频网站免费| 街头女战士在线观看网站| 中文字幕人妻熟女乱码| 精品久久久久久电影网| 我的亚洲天堂| 妹子高潮喷水视频| 日本黄色日本黄色录像| 日韩不卡一区二区三区视频在线| 男男h啪啪无遮挡| 美女视频免费永久观看网站| 黄色配什么色好看| 中文字幕人妻熟女乱码| 国产av国产精品国产| 人妻系列 视频| 欧美变态另类bdsm刘玥| 天天操日日干夜夜撸| 午夜av观看不卡| 国产日韩欧美视频二区| 人人妻人人添人人爽欧美一区卜| 色视频在线一区二区三区| 精品卡一卡二卡四卡免费| 国产精品99久久99久久久不卡 | 亚洲精品乱久久久久久| 亚洲综合色惰| 一级毛片 在线播放| www.熟女人妻精品国产| 久久青草综合色| 日韩成人av中文字幕在线观看| 免费观看无遮挡的男女| 亚洲情色 制服丝袜| 成人国语在线视频| 欧美bdsm另类| 亚洲精品乱久久久久久| 国产片内射在线| 免费观看性生交大片5| 少妇人妻精品综合一区二区| 亚洲,欧美,日韩| 亚洲精品自拍成人| 久久久久久久精品精品| www.自偷自拍.com| 免费女性裸体啪啪无遮挡网站| 国产成人91sexporn| 午夜日韩欧美国产| 免费少妇av软件| 寂寞人妻少妇视频99o| 中国三级夫妇交换| 亚洲欧美中文字幕日韩二区| 亚洲精华国产精华液的使用体验| 国产黄色免费在线视频| 午夜免费鲁丝| 精品99又大又爽又粗少妇毛片| 日韩,欧美,国产一区二区三区| 欧美人与性动交α欧美精品济南到 | 久久久久久伊人网av| 黄色毛片三级朝国网站| 丝瓜视频免费看黄片| 成人免费观看视频高清| 最近最新中文字幕免费大全7| 午夜老司机福利剧场| 精品亚洲乱码少妇综合久久| 欧美xxⅹ黑人| 成人18禁高潮啪啪吃奶动态图| 免费观看性生交大片5| 国产精品国产三级国产专区5o| 18+在线观看网站| 国产黄色免费在线视频| 国产一区二区 视频在线| 国产精品一二三区在线看| 熟女av电影| 国产精品一国产av| 午夜日韩欧美国产| 久久精品熟女亚洲av麻豆精品| 亚洲av欧美aⅴ国产| 男女边吃奶边做爰视频| 欧美成人精品欧美一级黄| www.av在线官网国产| 美女主播在线视频| 亚洲精品国产av成人精品| 国产成人精品一,二区| 免费女性裸体啪啪无遮挡网站| 国产男人的电影天堂91| 国产老妇伦熟女老妇高清| 成年动漫av网址| 国产av码专区亚洲av| 99香蕉大伊视频| 日韩制服骚丝袜av| 黑人猛操日本美女一级片| 青青草视频在线视频观看| 在线观看免费高清a一片| 美女福利国产在线| 欧美 亚洲 国产 日韩一| a 毛片基地| 老司机亚洲免费影院| 少妇 在线观看| 国产乱人偷精品视频| 日韩在线高清观看一区二区三区| 国产一区二区三区av在线| 亚洲精品久久午夜乱码| 久久精品国产亚洲av天美| 亚洲人成77777在线视频| 日韩不卡一区二区三区视频在线| 久久久久精品性色| kizo精华| 国产熟女午夜一区二区三区| 欧美人与性动交α欧美软件| 精品人妻在线不人妻| 亚洲精品自拍成人| 啦啦啦中文免费视频观看日本| 蜜桃国产av成人99| 久久人妻熟女aⅴ| 久久久久国产一级毛片高清牌| 久久久久精品性色| 精品酒店卫生间| 精品一区二区三区四区五区乱码 | 麻豆av在线久日| av天堂久久9| 丝袜在线中文字幕| 人人妻人人澡人人爽人人夜夜| 国产成人91sexporn| 在现免费观看毛片| 少妇的逼水好多| 一区二区三区精品91| 人人妻人人添人人爽欧美一区卜| 老女人水多毛片| 国产福利在线免费观看视频| 日本爱情动作片www.在线观看| 熟妇人妻不卡中文字幕| 日本免费在线观看一区| 久久综合国产亚洲精品| 天天躁夜夜躁狠狠躁躁| 电影成人av| 精品第一国产精品| videos熟女内射| 大话2 男鬼变身卡| 久久99热这里只频精品6学生| 飞空精品影院首页| 天堂中文最新版在线下载| 多毛熟女@视频| 激情视频va一区二区三区| 亚洲精品美女久久久久99蜜臀 | 91午夜精品亚洲一区二区三区| 久久ye,这里只有精品| 九色亚洲精品在线播放| 黄片小视频在线播放| 国产极品粉嫩免费观看在线| 亚洲欧美一区二区三区久久| 热99久久久久精品小说推荐| 国产精品无大码| 91精品伊人久久大香线蕉| 免费黄频网站在线观看国产| 视频在线观看一区二区三区| 大片免费播放器 马上看| 巨乳人妻的诱惑在线观看| freevideosex欧美| 好男人视频免费观看在线| 免费黄频网站在线观看国产| 两性夫妻黄色片| 色播在线永久视频| 久久亚洲国产成人精品v| 国产欧美日韩综合在线一区二区| 欧美 亚洲 国产 日韩一| 宅男免费午夜| 99久久中文字幕三级久久日本| 9热在线视频观看99| 尾随美女入室| 99久久中文字幕三级久久日本| 久久久久国产精品人妻一区二区| 欧美最新免费一区二区三区| √禁漫天堂资源中文www| 免费看av在线观看网站| 国产免费福利视频在线观看| 一二三四中文在线观看免费高清| 久久国产精品男人的天堂亚洲| 婷婷色麻豆天堂久久| 久久青草综合色| 一级爰片在线观看| 欧美日韩一级在线毛片| 美女国产高潮福利片在线看| 精品酒店卫生间| 大片电影免费在线观看免费| 久热久热在线精品观看| 国产成人精品无人区| 人人澡人人妻人| 涩涩av久久男人的天堂| 国产精品无大码| 婷婷色综合www| a级片在线免费高清观看视频| 视频区图区小说| 成年人午夜在线观看视频| 亚洲男人天堂网一区| 99热网站在线观看| 欧美亚洲日本最大视频资源| 人成视频在线观看免费观看| 久久 成人 亚洲| 日韩视频在线欧美| 91aial.com中文字幕在线观看| 一本久久精品| 亚洲精品美女久久久久99蜜臀 | 搡女人真爽免费视频火全软件| a级毛片黄视频| 看免费成人av毛片| 国产欧美日韩一区二区三区在线| 成人毛片a级毛片在线播放| 精品少妇久久久久久888优播| 涩涩av久久男人的天堂| 观看av在线不卡| 不卡av一区二区三区| 国产成人一区二区在线| av又黄又爽大尺度在线免费看| 亚洲精品在线美女| 高清黄色对白视频在线免费看| 最近中文字幕高清免费大全6| 天天躁狠狠躁夜夜躁狠狠躁| 日韩精品有码人妻一区| 一区二区三区激情视频| 交换朋友夫妻互换小说| 夜夜骑夜夜射夜夜干| 99re6热这里在线精品视频| 大码成人一级视频| 伊人久久大香线蕉亚洲五| 男的添女的下面高潮视频| 啦啦啦在线免费观看视频4| 日韩精品免费视频一区二区三区| 亚洲综合色惰| 在线观看www视频免费| 熟女电影av网| 国产精品一区二区在线不卡| 999精品在线视频| 亚洲成av片中文字幕在线观看 | av福利片在线| 亚洲精品国产色婷婷电影| 九色亚洲精品在线播放| 免费看不卡的av| 汤姆久久久久久久影院中文字幕| 丰满乱子伦码专区| 国产精品国产三级国产专区5o| 99re6热这里在线精品视频| 欧美日韩成人在线一区二区| 久久影院123| 久久精品久久久久久久性| 欧美老熟妇乱子伦牲交| 超色免费av| 久久综合国产亚洲精品| 最近2019中文字幕mv第一页| 五月开心婷婷网| 亚洲欧洲日产国产| 国产一区有黄有色的免费视频| 哪个播放器可以免费观看大片| 国产精品嫩草影院av在线观看| 十八禁网站网址无遮挡| av一本久久久久| 日日啪夜夜爽| 一区二区av电影网| 黄色 视频免费看| 满18在线观看网站| 免费在线观看视频国产中文字幕亚洲 | 最黄视频免费看| 人妻少妇偷人精品九色| 免费观看在线日韩| 一本久久精品| 最近中文字幕2019免费版| 午夜免费鲁丝| 高清视频免费观看一区二区| 国产视频首页在线观看| 精品一品国产午夜福利视频| 中文字幕色久视频| 午夜精品国产一区二区电影| 日产精品乱码卡一卡2卡三| 欧美 亚洲 国产 日韩一| 国产成人精品福利久久| 黄色怎么调成土黄色| h视频一区二区三区| 国产成人免费观看mmmm| 十八禁高潮呻吟视频| 欧美精品一区二区免费开放| 亚洲三区欧美一区| 欧美精品一区二区大全| 一本久久精品| 熟女电影av网| 国产在视频线精品| 国产av码专区亚洲av| 91在线精品国自产拍蜜月| 亚洲国产毛片av蜜桃av| 亚洲av成人精品一二三区| 丰满饥渴人妻一区二区三| 成人毛片60女人毛片免费| 国产免费一区二区三区四区乱码| 女人久久www免费人成看片| 在线观看国产h片| 国产xxxxx性猛交| 亚洲国产av新网站| 亚洲美女搞黄在线观看| 69精品国产乱码久久久| 国产无遮挡羞羞视频在线观看| 天天操日日干夜夜撸| 又黄又粗又硬又大视频| 午夜福利在线免费观看网站| 国产精品久久久久久av不卡| 伦理电影大哥的女人| 观看美女的网站| 少妇的逼水好多| 精品久久蜜臀av无| 亚洲第一av免费看| 成人国产av品久久久| 一个人免费看片子| 99热网站在线观看| 国产亚洲精品第一综合不卡| 国产一区二区在线观看av| 国产精品久久久久久久久免| 欧美人与性动交α欧美精品济南到 | 91成人精品电影| 最新中文字幕久久久久| 中文字幕亚洲精品专区| 男人操女人黄网站| 麻豆精品久久久久久蜜桃| 人成视频在线观看免费观看| 高清欧美精品videossex| 黄色一级大片看看| 免费观看无遮挡的男女| 亚洲中文av在线| 亚洲精品久久午夜乱码| 捣出白浆h1v1| 色视频在线一区二区三区| 国产又色又爽无遮挡免| 男女边吃奶边做爰视频| 一区福利在线观看| 久久久国产一区二区| 伊人久久国产一区二区| 少妇人妻久久综合中文| 日本91视频免费播放| 亚洲欧洲国产日韩| 熟妇人妻不卡中文字幕| 亚洲欧洲国产日韩| 十八禁高潮呻吟视频| 精品一区二区三区四区五区乱码 | 久久韩国三级中文字幕| 国产白丝娇喘喷水9色精品| 我的亚洲天堂| 欧美日本中文国产一区发布| 国精品久久久久久国模美| 91精品三级在线观看| 国产一区二区 视频在线| 免费看不卡的av| 欧美少妇被猛烈插入视频| 欧美97在线视频| 国产精品人妻久久久影院| 波多野结衣av一区二区av| 日韩制服骚丝袜av| 老女人水多毛片| 亚洲av.av天堂| 亚洲图色成人| 极品人妻少妇av视频| 18禁观看日本| 亚洲国产av新网站| 91成人精品电影| 十分钟在线观看高清视频www| 乱人伦中国视频| 精品久久久久久电影网| 国产日韩欧美亚洲二区| 十分钟在线观看高清视频www| 国产日韩欧美亚洲二区| 日韩三级伦理在线观看| 男人添女人高潮全过程视频| 欧美少妇被猛烈插入视频| 精品人妻在线不人妻| 欧美中文综合在线视频| 国产成人午夜福利电影在线观看| 一二三四在线观看免费中文在| 永久免费av网站大全| 在线 av 中文字幕| 成人毛片a级毛片在线播放| 男人爽女人下面视频在线观看| 国产成人精品无人区| 久久久久网色| 伊人亚洲综合成人网| 一级a爱视频在线免费观看| 久久免费观看电影| 国产精品一二三区在线看| 黄片无遮挡物在线观看| 亚洲国产成人一精品久久久| 国产成人91sexporn| 97在线人人人人妻| 日韩三级伦理在线观看| 一区二区日韩欧美中文字幕| 精品一区在线观看国产| 久久亚洲国产成人精品v| 久久综合国产亚洲精品| 精品人妻熟女毛片av久久网站| 精品少妇久久久久久888优播| 亚洲国产日韩一区二区| 好男人视频免费观看在线| 美女脱内裤让男人舔精品视频| 亚洲综合精品二区| 男的添女的下面高潮视频| 国产黄色视频一区二区在线观看| 麻豆av在线久日| 麻豆乱淫一区二区| 国产一区二区三区av在线| 国产精品久久久久久久久免| 国产精品秋霞免费鲁丝片| 国产精品一二三区在线看| 国产成人精品久久二区二区91 | 亚洲美女黄色视频免费看| 丰满少妇做爰视频| 另类亚洲欧美激情| 可以免费在线观看a视频的电影网站 | 中文字幕av电影在线播放| 一级爰片在线观看| 日本av免费视频播放| 亚洲一区中文字幕在线| 伦理电影大哥的女人| 久久久精品94久久精品| 国产在线一区二区三区精| 一级,二级,三级黄色视频| 久久精品国产综合久久久| 亚洲在久久综合| 午夜久久久在线观看| 亚洲精品中文字幕在线视频| 伊人亚洲综合成人网| 9191精品国产免费久久| 国产在线免费精品| 国产综合精华液| 尾随美女入室| 啦啦啦中文免费视频观看日本| 女性生殖器流出的白浆| 亚洲四区av| 男女下面插进去视频免费观看| 国产精品免费视频内射| 男人操女人黄网站| 精品亚洲成a人片在线观看| 最近最新中文字幕免费大全7| 国产又爽黄色视频| 青春草亚洲视频在线观看| 亚洲男人天堂网一区| 一二三四中文在线观看免费高清| 人人妻人人澡人人爽人人夜夜| 婷婷色综合www| 精品国产露脸久久av麻豆| 亚洲综合精品二区| 大话2 男鬼变身卡| 一级a爱视频在线免费观看| 欧美精品一区二区免费开放| 黄色 视频免费看| 久久精品久久久久久噜噜老黄| 新久久久久国产一级毛片| 自拍欧美九色日韩亚洲蝌蚪91| 久久精品熟女亚洲av麻豆精品| 性少妇av在线| 日韩电影二区| 精品第一国产精品| 国产成人精品在线电影| 女性被躁到高潮视频| 日韩av在线免费看完整版不卡| 久久久久久免费高清国产稀缺| 26uuu在线亚洲综合色| 欧美 日韩 精品 国产| 国产麻豆69| 亚洲国产精品国产精品| 一级片'在线观看视频| 国产精品国产av在线观看| 男女国产视频网站| 美女高潮到喷水免费观看| 午夜91福利影院| 制服丝袜香蕉在线| 欧美日韩精品网址| 欧美变态另类bdsm刘玥| 免费av中文字幕在线| 在线精品无人区一区二区三| 国产精品久久久久久av不卡| 国产成人精品无人区| 三级国产精品片| 中文字幕人妻熟女乱码| 国产一区有黄有色的免费视频| 国产激情久久老熟女| 久久av网站| 精品国产国语对白av| 永久网站在线| 另类精品久久| 桃花免费在线播放| 两个人免费观看高清视频| 在线观看美女被高潮喷水网站| 国产精品一区二区在线观看99| 老司机亚洲免费影院| 一区二区三区精品91| 国产深夜福利视频在线观看| 大片电影免费在线观看免费| 精品少妇黑人巨大在线播放| 一级毛片我不卡| a级毛片在线看网站| 国产精品久久久久久精品古装| 精品一区二区免费观看| 一区二区日韩欧美中文字幕| av在线播放精品| 黄频高清免费视频| 天堂俺去俺来也www色官网| 日韩制服骚丝袜av| 久久久精品免费免费高清| 我的亚洲天堂| 狠狠婷婷综合久久久久久88av| 一边摸一边做爽爽视频免费| 成人国产麻豆网| 国产免费又黄又爽又色| 精品第一国产精品| 精品国产露脸久久av麻豆| 97精品久久久久久久久久精品| 久久精品国产自在天天线| 巨乳人妻的诱惑在线观看| 亚洲情色 制服丝袜| 大陆偷拍与自拍| 精品少妇一区二区三区视频日本电影 | 国产成人精品无人区| 亚洲成人一二三区av| 亚洲精品乱久久久久久| 人人妻人人爽人人添夜夜欢视频| 一本久久精品| 国产精品偷伦视频观看了| 97在线人人人人妻| 制服人妻中文乱码| 美女脱内裤让男人舔精品视频| 久久人人爽人人片av| 亚洲国产av新网站| 少妇被粗大的猛进出69影院| 王馨瑶露胸无遮挡在线观看| a级毛片在线看网站| 九九爱精品视频在线观看| 国产激情久久老熟女| 九草在线视频观看| 水蜜桃什么品种好| 免费播放大片免费观看视频在线观看| 在线观看免费视频网站a站| 久久 成人 亚洲| 男的添女的下面高潮视频| 18禁动态无遮挡网站| xxx大片免费视频| 精品少妇黑人巨大在线播放| 国产精品熟女久久久久浪| 一区在线观看完整版| 亚洲三级黄色毛片| 少妇猛男粗大的猛烈进出视频| 电影成人av| 久久精品国产a三级三级三级| √禁漫天堂资源中文www| 国产免费现黄频在线看| 亚洲精品成人av观看孕妇| 亚洲av.av天堂| 国产成人一区二区在线| 欧美国产精品va在线观看不卡| 两性夫妻黄色片| 国产又爽黄色视频| 国产1区2区3区精品| 日韩人妻精品一区2区三区| 欧美成人午夜免费资源| 中文字幕最新亚洲高清| 国产熟女欧美一区二区| 精品国产超薄肉色丝袜足j| 桃花免费在线播放| 久久精品国产亚洲av涩爱| 美女国产高潮福利片在线看| 又大又黄又爽视频免费| 18+在线观看网站| 最近手机中文字幕大全| 黄片无遮挡物在线观看| 男人添女人高潮全过程视频| 国产成人a∨麻豆精品| 久久久亚洲精品成人影院| 十八禁高潮呻吟视频| 亚洲国产精品一区三区| 国产成人精品久久久久久| 纯流量卡能插随身wifi吗| 黑人猛操日本美女一级片| 精品国产一区二区三区久久久樱花| 性色av一级| 国产一区亚洲一区在线观看| 亚洲国产日韩一区二区| 一二三四在线观看免费中文在| 99热网站在线观看| 亚洲精品美女久久久久99蜜臀 | 国产欧美日韩综合在线一区二区| h视频一区二区三区| av国产久精品久网站免费入址| 亚洲精品日韩在线中文字幕| 天天躁狠狠躁夜夜躁狠狠躁| 亚洲,一卡二卡三卡| 亚洲一区二区三区欧美精品| 丁香六月天网| 久久精品久久精品一区二区三区| 久久免费观看电影| 国产野战对白在线观看| 观看美女的网站| 成人毛片a级毛片在线播放| 亚洲国产精品国产精品| 我的亚洲天堂| 国产有黄有色有爽视频| 国产精品久久久久久久久免| 欧美少妇被猛烈插入视频| 国产熟女欧美一区二区| 中文字幕亚洲精品专区| 制服人妻中文乱码| 亚洲国产精品一区二区三区在线| 纯流量卡能插随身wifi吗| 国产精品一区二区在线不卡| 寂寞人妻少妇视频99o| 久久久久久久大尺度免费视频| 又粗又硬又长又爽又黄的视频| 久久精品人人爽人人爽视色| 亚洲精品第二区| 亚洲av福利一区| 国产免费现黄频在线看| 在线观看www视频免费| 国产男女内射视频|