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

    Solvation Structure and Dynamics of Mg(TFSI)2Aqueous Electrolyte

    2022-04-15 11:49:32ZhouYuTaylorJuranXinyiLiuKeeSungHanHuiWangKarlMuellerLinMaKangXuTaoLiLarryCurtissandLeiCheng
    Energy & Environmental Materials 2022年1期
    關(guān)鍵詞:機采棉模范帶頭傾情

    Zhou Yu ,Taylor R.Juran,Xinyi Liu,Kee Sung Han,Hui Wang,Karl T.Mueller,Lin Ma ,Kang Xu ,Tao Li*,Larry A.Curtiss*,and Lei Cheng*

    Using ab initio molecular dynamics(AIMD)simulations,classical molecular dynamics(CMD)simulations,small-angle X-ray scattering(SAXS),and pulsed- field gradient nuclear magnetic resonance(PFG-NMR),the solvation structure and ion dynamics of magnesium bis(trifluoromethanesulfonyl)imide(Mg(TFSI)2)aqueous electrolyte at 1,2,and 3 m concentrations are investigated.From AIMD and CMD simulations,the first solvation shell of an Mg2+ion is found to be composed of six water molecules in an octahedral configuration and the solvation shell is rather rigid.The TFSI-ions prefer to stay in the second solvation shell and beyond.Meanwhile,the comparable diffusion coefficients of positive and negative ions in Mg(TFSI)2aqueous electrolytes have been observed,which is mainly due to the formation of the stable[Mg(H2O)6]2+complex,and,as a result,the increased effective Mg ion size.Finally,the calculated correlated transference numbers are lower than the uncorrelated ones even at the low concentration of 2 and 3 m,suggesting the enhanced correlations between ions in the multivalent electrolytes.This work provides a molecular-level understanding of how the solvation structure and multivalency of the ion affect the dynamics and transport properties of the multivalent electrolyte,providing insight for rational designs of electrolytes for improved ion transport properties.

    Keywords

    Mg(TFSI)2aqueous electrolyte,molecular dynamics simulation,pulsed- field gradient nuclear magnetic resonance,small-angle X-ray scattering,ion dynamics

    1.Introduction

    Multivalent-ion batteries play an important role in the quest for energy-dense alternatives to Liion batteries(LIBs).Since the 1970s,lithiumbased battery chemistries have been sought after as the ultimate solution to energy storage needs,with the eventual fruition of LIBs as the dominant component of current energy storage technologies.However,as the demand for an energy storage device of high energy,high safety,low-cost,and risk-free supply chain becomes ever imminent,battery chemistries based on earth-abundant elements rather than lithium emerge as candidates,targeting applications in personal electronics,vehicles,and grid storage.In the early 2000s,bivalent Mg was reported as such a promising candidate to substitute for Li.[1]The Mg battery concept brings benefits such as its divalent nature,which provides two electrons per redox center with a much higher theoretical volumetric capacity(3832 mAh cm-3)compared with monovalent ions such as Li+(2062 mAh cm-3)and Na+(1129 mAh cm-3).[2]However,the development of a functional Mg battery has encountered challenges in every component.Specifically,to Mg electrolytes,slow ion transport,limited electrochemical stability,and concomitant safety are key performance metrics that need to be improved.[3,4]

    Aqueous electrolytes have drawn increasing attention over the past few years,due to the contribution of water in electrolytes to activate the electrode and help the Mg2+insertion process[5]and the significantly improved electrochemical stability brought by the“water-in-salt”concept.[6]Chen et al.[7]proposed an aqueous Mg-ion battery composed of a Prussian blue-type nickel hexacyanoferrate cathode,polyimide anode,and 1M(molarity,molsoluteper Lsolution)MgSO4aqueous electrolyte.In addition to achieving the maximum cell voltage of 1.5 V and the comparable operating voltage to the unstable non-aqueous electrolyte,the battery cell maintains 60% capacity after 5000 cycles,with almost 100% Coulombic efficiency.Furthermore,a~40 Wh kg-1energy density was achieved at a 1 A g-1current.[7]Wang et al.[8]proposed a highly concentrated aqueous Mg-ion battery composed of Li3V2(PO4)3cathode,polypyromellitic dianhydride anode,and 4 m (molality,molsoluteper kgsolvent)magnesium bis(trifluoromethanesulfonyl)imide(Mg(TFSI)2)aqueous electrolyte.The stable electrochemical window is increased to 2.0 V.Additionally,92% initial capacity can be retained after 6000 cycles.Moreover,an unprecedented high power density of 6400 W kg-1was achieved,which is~30 times higher than the non-aqueous Mg-ion batteries in the literature.[8]

    The cation solvation structure and dynamics in liquid electrolytes dictate critically important properties essential to battery performance.For example,the ion-pairing can impact the ion transport and charge transfer,[9]while the solvation structure often predicts the interfacial chemistry and,consequently,the electrochemical stability.[6]The use of aqueous Mg electrolytes is still in its infancy and only limited efforts were reported.Buchner et al.[10]found the effective hydration number,including the first and second solvation shells around the Mg2+ion,to decrease as the Mg(SO4)2concentration is increased.Wahab et al.[11]proposed that ion pairs do not form in the 1MMg(NO3)2aqueous solution.Recently,Bucur et al.[4]proved that the hydration of Mg ion might benefit the intercalation process in oxides.

    While these elegant works explored the fundamental physicochemical properties of Mg aqueous electrolytes,many key questions to the design of these electrolytes,such as their structural and dynamical properties,remain unanswered.In this work,we aimed to study the relationship between the solvation and dynamics in divalent(i.e.,Mg2+ion)aqueous solutions and compare these properties with monovalent(i.e.,Li+ion)solution.We utilize ab initio and classical molecular dynamics(MD)simulations,small-angle X-ray scattering technique(SAXS),and pulsed- field gradient nuclear magnetic resonance measurements(PFG-NMR)to investigate the multivalent and concentration effects on the solvation structure and dynamics of the Mg(TFSI)2aqueous electrolyte.We found that a rigid octahedral hydration shell around the Mg2+ion leads to a significant cage effect at the ps timescale.Meanwhile,the long lifetime hydration shell around Mg2+ion and similar ion size of[Mg(H2O)6]2+and TFSI-ion could be the decisive factors in the comparable diffusion coefficients of cations and anions in Mg(TFSI)2aqueous electrolyte for various concentrations.Furthermore,we found the transport properties of electrolytes at high concentrations deviate from the Einstein relation.

    2.Results and Discussion

    2.1.Solvation Structure

    While a recent study has demonstrated that nearly all ions maintain the solvent-separated single-ion state in the 1 m LiTFSI aqueous electrolyte,[12]clear evidence of the solvation state of the ion in the divalent Mg(TFSI)2aqueous electrolyte remains absent in the literature.The initial configuration plays an important role in the AIMD simulations due to the limited simulation timescale used(i.e.,~100 ps).A carefully set up initial configuration helps the system reach a meaningful configuration faster.We performed AIMD simulations with two distinct initial salt configurations:solvent-separated single-ion and contacting-ion pair,respectively,to elucidate the ion behavior.The solvent-separated single Mg2+ion system maintains its configuration for the entire simulation(~50 ps).Conversely,in the simulation starting with the contacting ion pair configuration,we observe the Mg2+ion dissociating from the TFSI-ion(Figure 1a)within~12 ps of simulation time.These results suggest that ions prefer the solvent-separated single-ion state in the 1 m Mg(TFSI)2aqueous electrolyte.Therefore,the following analysis of the structure and dynamics was conducted for the solvent-separated single-ion structure for both the Li+and Mg2+ion systems.

    Figure 1.AIMD simulation results of 1 m LiTFSI and 1 m Mg(TFSI)2aqueous electrolytes.a)Time evolution of the distance between Mg2+ion and the coordinated oxygen from TFSI-ion.b)Radial distribution function g(r)and coordination number N(r)of water molecules around Li+and Mg2+ions.c)Histogram of the angle formed between oxygen in water,cation,and oxygen in the water.The probability ratio of the two marked peaks for the Mg case is 4:1.d)Potential of mean force of water molecules as a function of their distance from Li+and Mg2+ions.The black arrow between the dashed lines indicates the energy barrier for a water molecule to escape from the first hydration shell around Li+ion.

    To characterize the multivalent effects on the solvation structure around the cation,we compute the radial distribution function(RDF)and the coordination number of water molecules surrounding Li+and Mg2+ions(Figure 1b).In the Li-1S system,the first hydration shell is composed of 4 water molecules at a distance of~0.20 nm from the cation.The water molecules at ~0.40 nm from the Li+ion compose the second hydration shell,which corresponds to a much broader RDF peak than the first one in Figure 1b.In the Mg-1S system,the first hydration shell is composed of 6 water molecules at a slightly further distance of 0.21 nm.A distinct second hydration shell represented by another sharp peak is located at~0.43 nm radially from the Mg2+ion.The analysis on the histogram of the angles formed by oxygen in water,cation,and oxygen in the water of the first hydration shell(Figure 1c)shows that the mean angle is 109.0± 10.8°in the Li-1S system,which indicates that the four water molecules in the first hydration shell form a tetrahedral coordination environment.In the Mg-1S system,there are two peaks at angles of 90.0 ± 6.3°and 170.7± 5.0°,respectively,with the probability ratio of 4:1.This angle distribution suggests an octahedral coordination configuration around the cation in the first hydration shell.These observations are consistent with previous studies.[13]Furthermore,the Mg-1S system possesses a smaller full width at half maximum(FWHM)and a reduced standard deviation of angle distribution compared with the Li-1S system.This suggests that the oscillation of water molecules in the first hydration shell surrounding Mg2+ion is rather confined.In other words,the hydration shell around the Mg2+ion is more rigid compared with that surrounding the Li+ion.

    It is also essential to understand the thermodynamics of the exchange of water molecules around the cations,which can affect the intercalation process and dynamic properties.[14]We also computed the potential of mean force(PMF),which reflects the potential energy change in water molecules as a function of its distance from the cation.It is computed using PMF( r)=-kBTln[g( r)],where kBis Boltzmann’s constant,T is the absolute temperature,and g(r)is the RDF from the targeted cation to the water molecules.The water molecule in the first hydration shell around the Li+ion needs to overcome an energy barrier of~17 kJ mol-1(the black arrow in Figure 1d)to migrate to a neighboring hydration shell.No exchanges of water molecules occur between the two hydration shells at ~0.20 and ~0.43 nm surrounding the Mg2+ion over the course of the simulation period,which implies that a water molecule within the first hydration shell surrounding an Mg2+ion must overcome a larger energy barrier to migrate to a neighboring shell than that required for a hydrated Li+ion system.

    付江錄是十三師火箭農(nóng)場的個體工商戶。他先后成立了哈密江盛有限責任公司、哈密鼎舜有限責任公司,在此期間,他還從事過個體運輸和機采棉等工作。無論在何種崗位,他都是干一行愛一行,一心撲在工作上,兢兢業(yè)業(yè),勤勤懇懇,一絲不茍,各項工作想在前、干在前,充分起到了模范帶頭作用。他用自己200多萬元的資金幫扶了近40人脫貧致富,其中有8名火箭農(nóng)場的少數(shù)民族兄弟在他的傾情支持下,從一無所有踏上了小康之路。

    Furthermore,to characterize the charge-delocalization effects between the cation and the hydration shell,we performed Bader’s charge analysis[15-18]on one hundred frames that were extracted from a 10 ps trajectory of each simulation with a time interval of 100 fs.We found the Li and Mg ions had 0.90 and 1.75 positive charges,respectively.Meanwhile,the first hydration shell of the Li and Mg ions had 0.11 and 0.19 negative charges,respectively.Considering the distance between cation and O in water shown in Figure 1b,we can estimate the electrostatic interaction potential between the ion and hydration shell(Ecation-water)following Ecation-water=ke(qcationqwater/r),where keis Coulomb’s constant,qcationand qwaterare the atomic charges of the cation and hydration shell,respectively,and r is the distance between the cation and hydration shell.The ratio of EMg-waterto ELi-wateris~3.The strong electrostatic interaction between Mg2+ion and first hydration shell contributes to the high solvation energy of Mg2+ion is shown in Figure S2.Meanwhile,this observation is in line with our observation on the PMF in Figure 1d.That is,water molecules in the first hydration shell around an Mg2+ion need to overcome a much larger energy penalty in order to escape from this shell than that required for a hydrated Li+ion system.

    The dynamics of ions in high concentrations are sluggish and AIMD simulations are computationally expensive.Therefore,we use the CMD simulations to study the concentration effects on the structure and dynamics of Mg(TFSI)2aqueous electrolytes.The force field used in our CMD simulations is firstly validated against the structural property results from our AIMD simulations and experiments.From Figure 2a,we see that the RDF and the coordination number of oxygen atoms from water molecules around the Mg2+ion in Mg-1 system,calculated with CMD simulation,are consistent with those calculated with AIMD simulation in Mg-1S system.There is a slight difference of 0.01 nm on the location of the first peak in the RDF from Mg2+ion to the oxygen in water obtained from the two methods.Furthermore,the structure factor shifts calculated with CMD simulation have similar trends to the experimental SAXS data across the three different concentrations,as shown in Figure 2b1,b2.These results show that the OPLS-aa force fields in CMD capture the main structural features of Mg(TFSI)2aqueous electrolytes.Further decoupling of the structure factor by the interactions between different components of the electrolyte can be found in Figure S3.

    Figure 2.a)Radial distribution function g(r)and coordination number N(r)of oxygens in water around Mg2+ion in 1 m Mg(TFSI)2aqueous electrolyte calculated using AIMD and CMD.Structure factors of the 1,2,and 3 m aqueous electrolytes b1)are measured using SAXS and b2)calculated using CMD simulation.

    Using CMD,we analyzed the solvation environments of Mg2+ions,including water molecules and TFSI-ion,at various concentrations(Figure 3).We see two distinct hydration shells within a 0.5 nm radius of the Mg2+ion.The first hydration shell is composed of 6 water molecules and is maintained at all three concentrations considered.The radial density and coordination number of water molecules in the second hydration shell decrease as the salt concentration is increased.This is due to TFSI-ion appearing at~0.5 nm and replacing water molecules in the second hydration shell at increased concentrations(Figure 3b).The radial distribution and coordination number of oxygen atoms in TFSI-ion increase as the concentration of salt increases.We find two prominent peaks at 0.44 and 0.68 nm,respectively,in the radial density distribution of oxygen atom in TFSI-ion.We note that the two peaks correspond to the two oxygen atoms binding to the same sulfur in TFSI-ion(Figure 3b).

    Figure 3.Radial density distribution and coordination number of a)oxygens in water and b)oxygens in TFSI-ion around Mg2+ion in 1,2,and 3 m aqueous electrolyte calculated using CMD simulation.The two vertical black dashed lines in Figure 3 guide the locations of the two peaks.The inset TFSI-ion shows the two oxygen atoms binding with the same sulfur attributing to the two peaks.Blue,yellow,brown,green,red,and white balls denote the N,S,C,F,O,and H atoms.

    2.2.Dynamics

    There are two theoretically equivalent methods to quantify the self-diffusivity of MD simulation systems,which are known as the Green-Kubo method and the Einstein method.In practice,the Einstein method is often preferred due to the precision and reproducibility of the post-simulation data analysis.[19]The general procedure of the Einstein method is to calculate the slope of the diffusive regimes in the mean-square displacements(MSDs)versus time plot.Unfortunately,within the AIMD simulation timescale of less than~100 ps,the dynamics of many systems do not reach the diffusive regimes that are meaningful for the self-diffusivity calculations.[20]We use a statistical parameter introduced by Burov et al.[21]to distinguish the dynamics of Li+and Mg2+ion at 1 m concentration within the relatively short AIMD trajectories we have.Specifically,three successive coordinates of the targeting ion,with a time interval Δ, are used to define two displacement vectors:V1(t)=X(t+Δ)-X(t);V2(t)=X(t+2Δ)-X(t+Δ),where X(t)is the coordinate of a targeted ion at time t.A small angle between V1(t)and V2(t)would indicate that the ion moves along almost the same direction as the previous time step,and the movement is barely restricted by the solvation shell.We can consider such trajectories as forward ballistic motion.A large angle between V1(t)and V2(t)would indicate the ion movement changed its direction,because the forward motion is blocked by other atoms occupying the position.This means a cage effect of the solvation shell is prominent.

    We compute the ensemble average of the angle between these two displacement vectors as a function of the time interval for the two systems,and the results are shown in Figure 4a-d.Specifically,the angles between displacement vectors at different time t from the simulation trajectory were calculated and then used to obtain the average at a given time interval Δ.A similar angle distribution for Li+and Mg2+ion has been observed within 1 ps(Figure 4a,c).Initially,a high frequency on the distribution occurs around the small angles,which suggests that the cation maintains forward ballistic motion.The histograms of the angle probability distribution for the time interval of 10 fs are shown in Figure 4e.We can see that the ion does not undergo the Brownian motion(random walk of particles),which would correspond to an equal distribution at all angles,as indicated by the horizontal dotted line.As the angle distribution increases from 0 to ~180 degrees in ~0.4 ps,as guided by the white dashed line in Figure 4a,c,we depict the transition from forward ballistic to motion that is affected by the solvation cage,which limits the forward movement of the ion.The angle distribution of Li+and Mg2+ions stabilizes around 110 and 140 degrees in 8 ps,respectively(Figure 4f).On average,the angle distribution in the Mg-1S system is higher than the Li-1S system,suggesting that the rigid octahedral hydration shell surrounding Mg2+ion leads to a stronger cage effect,which hinders the ionic motion.

    Figure 4.Histograms of the displacement vector angle probability distribution as a function of the time interval Δ for the a)0-1 ps and b)1-10 ps in the Li-1S system,as well as c)0-1 ps and d)1-10 ps in the Mg-1S system from AIMD simulations.Histograms of angle probability distribution with a time interval of e)10 fs and f)8 ps.BM stands for the Brownian motion.White dashed lines in a,c)indicate the transition trend from forward ballistic motion to the cage effect dominating motion.White dotted lines in a-d)represent a time interval of 10 fs and 8 ps.The horizontal dotted line in e,f)represents the distribution of Brownian motion for comparison.

    CMD was used to explore the concentration effects on the dynamics of Mg(TFSI)2aqueous electrolyte.The lifetime of the water and TFSI-ion in the first two hydration shells was characterized by calculating the residence correlation function for water molecules within the first and second hydration shells,the oxygen atoms in TFSI-ion,and the TFSI-ion in the second hydration shell.The residence correlation function is defined as <c( 0)c( t)>,where c( t)is an indicator of the targeted particle’s category(hydration shell).The c( t)is defined as 1.0 if the targeted particle in the corresponding shell at time t=0 resides continuously in this shell by time t,and 0 if the atom or molecule is dissociated from the ion according to the cutoff distances defined by the valleys in the RDF plot.For the definition of the residence time of the TFSI-ion,as long as one of the four oxygen atoms from the ion remains in the solvation shell,we consider the ion resides.In Figure 5,we plot out the time evolution of the statistical average residence correlation function for all cations in the simulation box with a time interval of 1 ps.A rapid decay in correlation function indicates that the targeted particle has a stronger tendency to escape from the corresponding hydration shell.From Figure 5a-c,we see the residence correlation functions for water and TFSI-ion in the first two hydration shells are similar for the three concentrations.Water molecules in the first hydration shell do not leave the Mg2+ion throughout the timescale used for our simulations(~30 ns).This observation is consistent with a previous NMR study,which reports an exchange timescale for water on the order of microseconds for Mg2+in the first hydration shell.[22]Furthermore,we see that when compared with water molecules,the coordinated oxygen atoms in the TFSI-ion from the second hydration shell escapes from Mg2+ion more readily.However,the dissociation of one oxygen is usually replaced by another oxygen atom from the same TFSI-ion.As a result,the TFSI-ion in the second hydration shell“resides”for a longer time than the water molecule in the first and this is evidenced by the slower decay of TFSI-ion on the residence correlation function.If we define the relaxation time as the time it takes for <c( 0)c( t)> to decrease to 5% ,we see that the relaxation time increases as the concentration increases(Figure 5d).Interestingly,the concentration affects the relaxation time of TFSI-ion more pronouncedly than that of the water molecules,as evidenced by the steeper slope for the relaxation time of TFSI-ion as a function of concentration.

    Figure 5.Residence correlation function for water in the first(r<0.30 nm in Figure 3a)and second(0.30 nm<r<0.50 nm in Figure 3a)hydration shells,the TFSI-ion and the oxygen atom in TFSI-ion in the second(0.30 nm<r<0.50 nm in Figure 3b)hydration shell in a)1,b)2,and c)3 m concentration Mg(TFSI)2aqueous solutions.d)Relaxation time defined as the time when <c( 0)c( t)> =0.05 as a function of concentration.

    Then,we calculated the diffusion coefficients of Mg2+,TFSI-ions,and water molecules as a function of salt concentration from the CMD simulations using the Einstein method.From Figure 6a,we can see the diffusion coefficients of ions and water decrease with the increase in salt concentration.The diffusion coefficients of Mg2+and TFSI-are comparable within the error bars,and water diffusion is about 3-4 times faster than for both Mg2+and TFSI-for various concentrations.All these critical features have also been captured in the PFG-NMR experiments(Figure 6b).The error in the diffusion coefficient of Mg2+ions is slightly larger than that of TFSI-anions and H2O molecules due to the lower signal intensity,that is,the signal-to-noise ratio in the25Mg PFG-NMR(inset of Figure S1).Here,we note that the diffusion coefficients calculated in MD simulation are 3-4 times larger than those measured in experiments.This is mainly because the TIP3P force field and the small scaling factor on the atomic charges in the force field of ions(i.e.,0.8)can significantly overestimate the mobility of water molecules and ions.The diffusion coefficients were also calculated in the MD simulations with a scaling factor of 0.9 and 1.0 on the atomic charges in the force fields of ions(see Figure S4).The key features,such as the high diffusion coefficient of water molecules and comparable diffusion coefficient of cation and anion,remain.

    The similar diffusion coefficients of the anion and cation in the aqueous Mg(TFSI)2electrolytes differ from that in the LiTFSI aqueous electrolyte,in which the diffusion coefficient of Li+is 50-130% larger than that of the TFSI-ion at 1-21 m concentrations using CMD simulations and PFG-NMR measurements.[12,23]This feature can be explained through the electrostatic interaction and ion size.The strong electrostatic field exerted by the divalent cation leads to a strong correlation between cation and anion.Therefore,the diffusion coefficient differences should be less than those in the monovalent electrolyte with similar concentration.Considering that the rigid octahedral hydration shell has a relatively long lifetime,the Mg2+ions in the electrolyte exist as[Mg(H2O)6]2+complex cations,which increases the “effective”ion size.The size of the first hydration shell around Mg2+ion(i.e.,0.4-0.5 nm)is similar to the distance between carbon atoms in one TFSI-ion(inset of Figure 3b),which also contributes to the comparable dynamics of ions.

    Figure 6.Diffusion coefficient of Mg2+,TFSI-,water,and the ratio of Mg2+and TFSI-as a function of concentration calculated in a)simulation and b)measured through PFG-NMR.Diffusion coefficients of TFSI-and water were calculated based on the center of mass of TFSI-ion and the oxygen atom in a water molecule in a CMD simulation.c)Uncorrelated and correlated transference numbers(t+uncand t+c)calculated in the CMD simulations and the transference number measured through PFG-NMR.d)The diffusion coefficient and viscosity relationship.Three points represent the experimental diffusivity and viscosity,and the black curve was calculated based on the Stokes-Einstein relations by taking the experimental diffusivity and viscosity of 1 m electrolyte as the reference.

    Furthermore,we measured the viscosities of the electrolytes and explored diffusion coefficient-viscosity relations.First,we estimated the diffusion coefficient of the electrolytes by the approximation D=xcDc+xaDa+xwDw,where xc,xa,xw,Dc,Da,and Dwrepresent the mole fractions and experimental diffusion coefficients of cation,anion,and solvent,respectively.[28]If the electrolyte is in the dilute regime and the correlation between ions can be ignored,the ions can be viewed as independent in their respective motions,and the Stokes-Einstein relation holds as ηiDi=ηjDj,where ηi,ηj,Di,and Djrepresent the viscosity and diffusion coefficient of the two systems i and j,respectively.By taking the experimental diffusion coefficient and viscosity of 1 m electrolyte as the reference,we can calculate the relationships between the viscosity and diffusion coefficient shown as the black curve in Figure 6d.However,the experimental viscosity deviates from the Stokes-Einstein relation as the concentration increases.For example,the Stokes-Einstein relation underestimates the viscosity by~2.3% and~7.3% at 2 and 3 m,respectively,due to ignoring the ionic correlations.These observations on the correlated transference number and viscosity once again remind us to exercise caution when applying these transport laws,which were originally derived for dilute/ideal electrolytes,to the highly concentrated regimes.

    3.Conclusion

    Using AIMD and CMD simulations,SAXS,and PFG-NMR,we investigated the solvation structure and dynamics of the cation in the divalent Mg(TFSI)2aqueous electrolyte and in the LiTFSI aqueous electrolyte for comparison.With AIMD,we found the octahedral hydration shell composed of 6 water molecules around Mg2+ion is more rigid than the tetrahedral shell composed of 4 water molecules around Li+ion.The strong electrostatic interaction between Mg2+ion and hydration shell is an important component of the high solvation energy and results in slow water exchange between hydration shells.Using experimentally verified CMD,we found TFSI-ions can only exist in the second solvation shell or beyond.The dynamics of 1,2,and 3 m Mg(TFSI)2aqueous electrolytes were further studied using CMD simulations.We found that the two oxygen atoms from the same TFSI-ion can interexchange;the TFSI-ion itself does not escape the second solvation shell readily.Unlike the LiTFSI electrolyte in which the cation dynamic is faster than the anion,the diffusion coefficients of the Mg cation and the anion are comparable.This is due to the formation of the stable[Mg(H2O)6]2+complex,resulting in the increased effective size that slows down the cation.As expected,the diffusion coefficient of water is significantly higher than the diffusion coefficients of both ions.Remarkably,the calculated correlated transference numberis much smaller than the uncorrelatedeven at the low concentrations of 2 and 3 m.This highlights the highly nonideal and correlated nature of the multivalent salts.

    This work has contributed to the understanding of the solvation structure and the resulting dynamic behaviors in the multivalent electrolytes.We showed that the strong electrostatic interaction between a multivalent ion and its negatively charged solvation shell and the counterion leads to the enlarged effective size and strong correlation of the Mg2+ion,ultimately slowing down the ion diffusion and decreasing correlated transference number.A further design goal for fast ion-conducting multivalent electrolytes should focus on achieving a moderate solvent-ion interaction and,consequently,a more flexible solvation shell through careful selections of solvents or diluents.

    4.Computational and Experimental Methods

    Ab Initio/Classical MD Simulations:Ab initio molecular dynamics(AIMD)simulations were performed on 1 m LiTFSI and 1 m Mg(TFSI)2aqueous electrolyte models,composed of 1 salt molecule and 55 water molecules.Classical MD(CMD)simulations were performed on 1,2,and 3 m Mg(TFSI)2aqueous electrolyte solutions with 90,180,and 270 salt molecules in the simulation models,respectively.Initial configurations were first set up using Packmol code.[29]The simulation box is periodic in all three directions.Additional details of the AIMD and CMD system setups are summarized in Table 1.

    Table 1.Setup of the AIMD and CMD systems studied in this work.

    The AIMD simulations were carried out using the Born-Oppenheimer approximation with the VASP code.[30-32]The projector augmented wave(PAW) method[33]was used to compute the interatomic forces with the Perdew--Burke-Ernzerhof(PBE)generalized-gradient approximation[34]for the exchangecorrelation energy.The plane-wave energy cutoff was set as 500 eV,and the Brillouin zone was sampled at the Γ-point.All AIMD simulations were performed with the NVT ensemble.The Nosé-Hoover thermostat[35]was used to stabilize the temperature of the simulations at 298 K.A time step of 0.5 fs was utilized for all simulations.Finally,simulations were run for a total time of 60 ps.The initial 20 ps trajectory was used to equilibrate the system.The standard deviation of the system energy during the last 5 ps of the equilibration run is less than 1 meV atom-1ps-1,which is small enough that the system can be considered to reach equilibration.[36]The 40 ps trajectory that follows was then considered as the “production”run,and the results were used for analysis.

    The CMD simulations were performed using GROMACS 5.1.4 with a 2 fs time step.[37]The bonded and non-bonded parameters of the force fields for Li+,Mg2+,and TFSI-ion were obtained from previous studies,[38,39]which were developed and optimized within the framework of an OPLS-aa force field.[40]The atomic charges in the force field of ions were scaled by 0.8 to compensate for the charge transfer,polarization,and binding effects.[41]The TIP3P model was employed for the water molecules.[42]Simulations were initially equilibrated under an NPT ensemble for 20 ns.Then,the production simulations were executed with an NVT ensemble for 250 ns.The temperature of the system was maintained at 298 K using the velocity-rescaling thermostat.[43]Non-electrostatic interactions were computed by direct summation with a cutoff length of 1.2 nm.Electrostatic interactions were computed using the particle mesh Ewald(PME)method.[44]The real space cutoff and Fourier spacing were 1.2 and 0.12 nm,respectively.The box size and periodic boundary condition of the CMD simulations have been evaluated by comparing the results of simulations with a variety of cell sizes.

    SAXS Measurements:Mg(TFSI)2was purchased from Sigma-Aldrich.Aqueous electrolytes were prepared by molality(molsoluteper kgsolvent).SAXS experiments were performed at the Advanced Photon Source(APS)12ID-B and C station of Argonne National Laboratory.The 2D SAXS data were collected on a PILATUS 2 M area detector(Dectris Ltd.)with a 2 meter distance from the sample to the detector and the incident energy of 13.3 keV.The two-dimensional scattering images were radially averaged over all orientations to produce plots of scattered intensity I(q)versus scattering vector q,where q=4π sin θ/λ.The scattering vector was calibrated using silver behenate.The samples were loaded into 1.5 mm diameter quartz capillary tubes for the SAXS measurements.

    PFG-NMR Measurements:PFG-NMR measurements were performed on a 600 MHz NMR spectrometer(Agilent,USA)with a 5 mm liquid NMR probe(Doty Scientific,USA),which can generate a z-gradient up to~31 T m-1.Diffusion coefficients of Mg2+cations,TFSI-anions,and H2O molecules were determined from25Mg,19F and,1H PFG-NMR,respectively,at 25°C.All diffusion coefficients were measured using the bipolar gradient stimulated echo sequence(Dbppste,vendor-supplied sequence,VNMRJ,Agilent).The PFG-echo profiles were recorded as a function of gradient strength with 16 equal steps and fitted with the Stejskal-Tanner equation:[45]S(g)=S(0)exp[-D(γgδ)2(Δ - δ/3)],where S(g)and S(0)are the echo heights at the gradient strengths of g and 0,D is the diffusion coefficient,γ is the gyromagnetic ratios of25Mg,19F or1H,Δ is the diffusion delay,which is the time interval between the two pairs of bipolar pulse gradients,and δ is the gradient length.The time interval Δ was 15 ms for25Mg PFGNMR and 40 ms for both1H and19F PFG-NMR.The gradient length δ was 2 ms for all measurements.The maximum gradient strength required for25Mg PFGNMR was~19 T m-1(see Figure S1),which is stronger approximately 16 times than that of1H/19F PFG-NMR(~1.2 T m-1)due to the fact that γ(1H)/γ(25Mg)≈16.4.

    Viscosity Measurements:Viscosity was measured using an Ostwald viscometer(Sibata Scientific Technology,Japan)with a capillary diameter of 0.75 mm.The temperature of the electrolyte in the viscometer was controlled by a circulating water bath(Sibata Scientific Technology,Japan)at 25°C.

    Acknowledgements

    This research was supported by the Joint Center for Energy Storage Research(JCESR),a U.S.Department of Energy,Energy Innovation Hub.The submitted manuscript has been created by UChicago Argonne,LLC,Operator of Argonne National Laboratory(“Argonne”).Argonne,a U.S.Department of Energy Office of Science laboratory,is operated under Contract no.DE-AC02-06CH11357.This research used resources of the Advanced Photon Source,a U.S.Department of Energy(DOE)Office of Science User Facility,operated for the DOE Office of Science by Argonne National Laboratory under Contract No.DE-AC02-06CH11357.The PFG-NMR measurements were performed at the Environmental Molecular Sciences Laboratory(EMSL),a national scientific user facility,sponsored by the DOE’s Office of Biological and Environmental Research and located at Pacific Northwest National Laboratory(PNNL).We acknowledge Dr.Yuyue Zhao for the density measurement of Mg(TFSI)2aqueous electrolytes.We gratefully acknowledge the computing resources provided on Bebop,a high-performance computing cluster,operated by the Laboratory Computing Resource Center at Argonne National Laboratory.

    Conflict of Interest

    The authors declare no conflict of interest.

    Supporting Information

    Supporting Information is available from the Wiley Online Library or from the author.

    猜你喜歡
    機采棉模范帶頭傾情
    如何發(fā)揮黨員在基層工作中的模范帶頭作用
    新疆北疆植棉區(qū)機采棉品種篩選試驗
    棉花科學(2017年1期)2017-03-10 20:38:43
    增強企業(yè)黨員先鋒模范作用的思考與探討
    機采棉打模機和運模車的故障維修
    短季棉品種華棉3109在荊門的種植表現(xiàn)
    為君傾情為君妖
    傾情大學生健康成長
    中國火炬(2014年8期)2014-07-24 14:30:18
    唉,我真后悔呀!
    新疆實施機采棉之研究
    心系學生 傾情關(guān)愛
    中國火炬(2010年2期)2010-07-24 14:36:02
    91成人精品电影| 欧美另类一区| 18禁黄网站禁片午夜丰满| 国产精品99久久99久久久不卡| 久久精品国产亚洲av高清一级| 中文字幕人妻丝袜一区二区| 12—13女人毛片做爰片一| www日本在线高清视频| av天堂久久9| 777米奇影视久久| 日本撒尿小便嘘嘘汇集6| 99精品久久久久人妻精品| 久久精品成人免费网站| 最近中文字幕2019免费版| 捣出白浆h1v1| 久久中文看片网| 精品一品国产午夜福利视频| 精品亚洲成国产av| 亚洲精品一区蜜桃| 男人舔女人的私密视频| 亚洲精品中文字幕在线视频| 精品少妇一区二区三区视频日本电影| a级毛片黄视频| 王馨瑶露胸无遮挡在线观看| 久久国产亚洲av麻豆专区| 国产成人免费无遮挡视频| 黑人欧美特级aaaaaa片| 欧美精品av麻豆av| 一级毛片电影观看| 中文字幕色久视频| 自拍欧美九色日韩亚洲蝌蚪91| 精品国产乱码久久久久久小说| 69精品国产乱码久久久| 午夜精品久久久久久毛片777| 精品熟女少妇八av免费久了| 亚洲美女黄色视频免费看| 亚洲一卡2卡3卡4卡5卡精品中文| a在线观看视频网站| 亚洲精品国产精品久久久不卡| 国产精品一区二区在线观看99| 亚洲av日韩精品久久久久久密| 国产一区二区三区综合在线观看| 99精品欧美一区二区三区四区| 久久av网站| a级片在线免费高清观看视频| av欧美777| 国产真人三级小视频在线观看| 亚洲欧美色中文字幕在线| 亚洲欧美色中文字幕在线| 欧美av亚洲av综合av国产av| 一个人免费看片子| 国产真人三级小视频在线观看| 波多野结衣av一区二区av| 国产欧美日韩一区二区三区在线| 飞空精品影院首页| 午夜精品国产一区二区电影| 亚洲国产欧美日韩在线播放| 秋霞在线观看毛片| 国产精品久久久久久精品古装| 国产视频一区二区在线看| 秋霞在线观看毛片| 一级a爱视频在线免费观看| 两性午夜刺激爽爽歪歪视频在线观看 | 我的亚洲天堂| 久久久精品94久久精品| av天堂久久9| 日日爽夜夜爽网站| 午夜福利视频精品| 人妻人人澡人人爽人人| 天天躁日日躁夜夜躁夜夜| 热re99久久国产66热| 国产成人免费无遮挡视频| 91精品国产国语对白视频| netflix在线观看网站| 欧美乱码精品一区二区三区| 国产成人精品久久二区二区免费| 少妇猛男粗大的猛烈进出视频| 精品久久久精品久久久| 99国产精品一区二区三区| 精品福利观看| 午夜两性在线视频| 男人舔女人的私密视频| 在线观看一区二区三区激情| 精品亚洲乱码少妇综合久久| 国产成人系列免费观看| 久久精品aⅴ一区二区三区四区| 12—13女人毛片做爰片一| 亚洲专区字幕在线| 两性午夜刺激爽爽歪歪视频在线观看 | 日韩 欧美 亚洲 中文字幕| videos熟女内射| 国产日韩欧美在线精品| 青春草视频在线免费观看| 丁香六月天网| 欧美黑人精品巨大| 欧美97在线视频| kizo精华| 国产成人精品在线电影| 精品国产一区二区三区久久久樱花| 国产精品自产拍在线观看55亚洲 | 十八禁高潮呻吟视频| 国产深夜福利视频在线观看| 欧美av亚洲av综合av国产av| 肉色欧美久久久久久久蜜桃| 精品一品国产午夜福利视频| 国产主播在线观看一区二区| 亚洲第一青青草原| 一本一本久久a久久精品综合妖精| 成年人免费黄色播放视频| 女性生殖器流出的白浆| 99精国产麻豆久久婷婷| 精品福利永久在线观看| 免费看十八禁软件| 一区在线观看完整版| e午夜精品久久久久久久| 多毛熟女@视频| 天堂8中文在线网| 午夜91福利影院| 免费日韩欧美在线观看| 一边摸一边抽搐一进一出视频| 男人爽女人下面视频在线观看| 99热全是精品| 另类精品久久| 91成人精品电影| 国产精品香港三级国产av潘金莲| 久久av网站| 美女脱内裤让男人舔精品视频| 69av精品久久久久久 | 交换朋友夫妻互换小说| 黄片大片在线免费观看| 一区二区三区精品91| 欧美日韩一级在线毛片| 大码成人一级视频| 制服人妻中文乱码| 亚洲成人手机| 下体分泌物呈黄色| 亚洲国产成人一精品久久久| 国产真人三级小视频在线观看| 久久久欧美国产精品| 久久久久精品人妻al黑| 男女免费视频国产| 免费高清在线观看视频在线观看| 别揉我奶头~嗯~啊~动态视频 | 秋霞在线观看毛片| 丝瓜视频免费看黄片| av欧美777| www.999成人在线观看| 肉色欧美久久久久久久蜜桃| 最新的欧美精品一区二区| 久久精品久久久久久噜噜老黄| av又黄又爽大尺度在线免费看| 国产亚洲欧美精品永久| xxxhd国产人妻xxx| 欧美性长视频在线观看| 亚洲欧美激情在线| 自线自在国产av| 少妇精品久久久久久久| 午夜福利一区二区在线看| 大香蕉久久网| 狠狠婷婷综合久久久久久88av| 亚洲人成电影观看| 极品人妻少妇av视频| 久久国产精品人妻蜜桃| 午夜福利影视在线免费观看| 老熟妇仑乱视频hdxx| 国产成人系列免费观看| 欧美日本中文国产一区发布| 亚洲一码二码三码区别大吗| 日本一区二区免费在线视频| 黑人猛操日本美女一级片| 一区二区三区四区激情视频| 下体分泌物呈黄色| 老熟女久久久| 亚洲第一青青草原| 满18在线观看网站| 麻豆国产av国片精品| 中文字幕色久视频| 久久青草综合色| 欧美日韩精品网址| 嫩草影视91久久| 久久久精品免费免费高清| 一个人免费在线观看的高清视频 | 欧美少妇被猛烈插入视频| 久久久国产精品麻豆| 欧美日韩亚洲综合一区二区三区_| 久久精品久久久久久噜噜老黄| 欧美精品啪啪一区二区三区 | 婷婷丁香在线五月| 欧美日韩亚洲高清精品| 1024视频免费在线观看| 一区二区三区乱码不卡18| 在线永久观看黄色视频| 别揉我奶头~嗯~啊~动态视频 | 中文欧美无线码| 精品久久久精品久久久| 日本黄色日本黄色录像| 国产男人的电影天堂91| 18禁观看日本| 亚洲精品日韩在线中文字幕| 中文字幕精品免费在线观看视频| 亚洲国产精品999| e午夜精品久久久久久久| 亚洲欧美日韩另类电影网站| 中文字幕另类日韩欧美亚洲嫩草| 人人澡人人妻人| 新久久久久国产一级毛片| 丝袜人妻中文字幕| 五月天丁香电影| 亚洲久久久国产精品| 老司机午夜福利在线观看视频 | 999久久久国产精品视频| 日本a在线网址| 人人妻人人添人人爽欧美一区卜| 男女无遮挡免费网站观看| 久久久国产成人免费| 在线十欧美十亚洲十日本专区| 久久中文看片网| 成人国产av品久久久| 精品人妻一区二区三区麻豆| 久9热在线精品视频| www日本在线高清视频| 捣出白浆h1v1| 丰满人妻熟妇乱又伦精品不卡| 老熟女久久久| 国产欧美日韩一区二区精品| 亚洲精品一卡2卡三卡4卡5卡 | 久久综合国产亚洲精品| 精品视频人人做人人爽| 日本撒尿小便嘘嘘汇集6| 亚洲熟女毛片儿| 国产亚洲精品久久久久5区| 国产成人av激情在线播放| 午夜影院在线不卡| 亚洲伊人久久精品综合| 午夜精品国产一区二区电影| 两性夫妻黄色片| 免费在线观看视频国产中文字幕亚洲 | 亚洲第一欧美日韩一区二区三区 | 搡老岳熟女国产| 久久久国产成人免费| 下体分泌物呈黄色| 一个人免费在线观看的高清视频 | 日韩精品免费视频一区二区三区| 桃花免费在线播放| 伊人久久大香线蕉亚洲五| 亚洲av国产av综合av卡| 丝袜脚勾引网站| 午夜福利免费观看在线| 精品国产乱子伦一区二区三区 | 国产99久久九九免费精品| 亚洲国产看品久久| 亚洲成国产人片在线观看| 婷婷色av中文字幕| 日韩大码丰满熟妇| 国产精品1区2区在线观看. | 在线观看免费视频网站a站| 黄色片一级片一级黄色片| 亚洲精品成人av观看孕妇| 精品国产乱码久久久久久男人| xxxhd国产人妻xxx| 人妻一区二区av| 一级片'在线观看视频| 久久人人97超碰香蕉20202| 国产黄色免费在线视频| 97人妻天天添夜夜摸| 成在线人永久免费视频| 一级毛片精品| 91精品三级在线观看| 国产精品欧美亚洲77777| 极品少妇高潮喷水抽搐| 一区二区av电影网| 成人av一区二区三区在线看 | 男人舔女人的私密视频| 精品福利永久在线观看| 夜夜骑夜夜射夜夜干| 国产精品九九99| 欧美大码av| 久久国产精品男人的天堂亚洲| 三级毛片av免费| 久久久久久免费高清国产稀缺| 亚洲av日韩在线播放| 真人做人爱边吃奶动态| 精品人妻熟女毛片av久久网站| 一级片'在线观看视频| 人人妻,人人澡人人爽秒播| 中文欧美无线码| 九色亚洲精品在线播放| 男女免费视频国产| 欧美性长视频在线观看| 少妇裸体淫交视频免费看高清 | 亚洲精品美女久久av网站| 精品一区在线观看国产| 欧美变态另类bdsm刘玥| 国产男女内射视频| 久久国产精品男人的天堂亚洲| 国产麻豆69| 日韩中文字幕欧美一区二区| 18在线观看网站| 男人爽女人下面视频在线观看| 精品国产一区二区三区久久久樱花| 日日夜夜操网爽| 国产伦理片在线播放av一区| 国产高清国产精品国产三级| 亚洲视频免费观看视频| 久久久久精品国产欧美久久久 | 日韩制服骚丝袜av| 黄色 视频免费看| a在线观看视频网站| 老司机在亚洲福利影院| 美女大奶头黄色视频| 免费在线观看视频国产中文字幕亚洲 | 搡老乐熟女国产| 男女无遮挡免费网站观看| 首页视频小说图片口味搜索| 波多野结衣av一区二区av| 久久这里只有精品19| 悠悠久久av| 国产黄频视频在线观看| 性色av一级| 这个男人来自地球电影免费观看| 亚洲精品成人av观看孕妇| 纵有疾风起免费观看全集完整版| 曰老女人黄片| 啪啪无遮挡十八禁网站| 欧美精品亚洲一区二区| 亚洲av欧美aⅴ国产| 热99国产精品久久久久久7| 天天添夜夜摸| 美国免费a级毛片| 女性生殖器流出的白浆| 精品少妇内射三级| 午夜福利免费观看在线| 久久精品国产综合久久久| 欧美日韩中文字幕国产精品一区二区三区 | 91麻豆av在线| 美女午夜性视频免费| 欧美人与性动交α欧美精品济南到| 亚洲熟女毛片儿| 最近最新免费中文字幕在线| 国产精品99久久99久久久不卡| 永久免费av网站大全| 中文字幕精品免费在线观看视频| 日本a在线网址| 亚洲精品国产av成人精品| 国产精品一区二区在线观看99| 中文字幕av电影在线播放| 99re6热这里在线精品视频| 青春草亚洲视频在线观看| 国产又色又爽无遮挡免| 国产人伦9x9x在线观看| 国产xxxxx性猛交| 各种免费的搞黄视频| 日本五十路高清| 亚洲国产av新网站| 丁香六月天网| 国产男人的电影天堂91| 国产主播在线观看一区二区| 国产高清videossex| h视频一区二区三区| 熟女少妇亚洲综合色aaa.| 国产xxxxx性猛交| 99热网站在线观看| videos熟女内射| 两人在一起打扑克的视频| 精品亚洲乱码少妇综合久久| 首页视频小说图片口味搜索| 精品一区二区三卡| av网站在线播放免费| 国产又爽黄色视频| 亚洲欧洲日产国产| 亚洲精品一二三| av天堂久久9| 久久中文看片网| 丰满迷人的少妇在线观看| 精品乱码久久久久久99久播| 亚洲欧美日韩高清在线视频 | 操出白浆在线播放| 国产在线一区二区三区精| 日韩制服骚丝袜av| 亚洲欧美色中文字幕在线| 精品欧美一区二区三区在线| 成年动漫av网址| 多毛熟女@视频| 国产精品国产av在线观看| 伦理电影免费视频| 亚洲久久久国产精品| 亚洲精品久久午夜乱码| 五月开心婷婷网| 丝袜美足系列| 性色av乱码一区二区三区2| 老司机深夜福利视频在线观看 | 啦啦啦啦在线视频资源| 免费观看a级毛片全部| 欧美日韩亚洲高清精品| 精品久久久久久电影网| 亚洲第一av免费看| svipshipincom国产片| 少妇被粗大的猛进出69影院| 青青草视频在线视频观看| 男女床上黄色一级片免费看| 在线观看免费高清a一片| 色94色欧美一区二区| 丰满少妇做爰视频| 大型av网站在线播放| a 毛片基地| 亚洲精品自拍成人| 999精品在线视频| 精品国产乱码久久久久久小说| 丝瓜视频免费看黄片| 99久久精品国产亚洲精品| 国产亚洲精品第一综合不卡| 精品欧美一区二区三区在线| 久久毛片免费看一区二区三区| 精品第一国产精品| 女性生殖器流出的白浆| 高清欧美精品videossex| 男女下面插进去视频免费观看| 欧美精品亚洲一区二区| 中文欧美无线码| 中文字幕制服av| 亚洲伊人久久精品综合| 美女中出高潮动态图| av一本久久久久| 欧美老熟妇乱子伦牲交| 国精品久久久久久国模美| 一个人免费在线观看的高清视频 | 久久精品人人爽人人爽视色| 日韩免费高清中文字幕av| 国产免费福利视频在线观看| 又紧又爽又黄一区二区| 麻豆国产av国片精品| 丝瓜视频免费看黄片| 亚洲成国产人片在线观看| 久久精品熟女亚洲av麻豆精品| 伦理电影免费视频| 啦啦啦 在线观看视频| 9191精品国产免费久久| 亚洲一码二码三码区别大吗| 国产视频一区二区在线看| 一区在线观看完整版| 亚洲精品一区蜜桃| 国产男女内射视频| 男女床上黄色一级片免费看| 精品国产国语对白av| 天堂中文最新版在线下载| 日韩 亚洲 欧美在线| 久久天躁狠狠躁夜夜2o2o| 久久久久久免费高清国产稀缺| 99热国产这里只有精品6| 韩国高清视频一区二区三区| 99热国产这里只有精品6| 丝袜喷水一区| 国产亚洲午夜精品一区二区久久| 欧美黄色片欧美黄色片| 久久国产亚洲av麻豆专区| 少妇裸体淫交视频免费看高清 | 人妻一区二区av| 欧美大码av| 国产一卡二卡三卡精品| 国产色视频综合| 亚洲欧美成人综合另类久久久| 99热全是精品| 麻豆av在线久日| 午夜激情av网站| 肉色欧美久久久久久久蜜桃| 少妇猛男粗大的猛烈进出视频| 国产精品二区激情视频| 免费av中文字幕在线| 青春草视频在线免费观看| 亚洲av男天堂| 12—13女人毛片做爰片一| 黄片大片在线免费观看| 老司机亚洲免费影院| 精品人妻在线不人妻| 午夜老司机福利片| 久热爱精品视频在线9| 人人妻人人澡人人爽人人夜夜| 日本vs欧美在线观看视频| 中文字幕制服av| 日本av免费视频播放| 美女国产高潮福利片在线看| 国产一区二区三区在线臀色熟女 | 久久 成人 亚洲| 青春草视频在线免费观看| 自拍欧美九色日韩亚洲蝌蚪91| 午夜福利视频在线观看免费| 免费不卡黄色视频| 国产一级毛片在线| 91国产中文字幕| 99热全是精品| 18在线观看网站| 久久人人爽人人片av| av福利片在线| 欧美日韩黄片免| 久久九九热精品免费| 建设人人有责人人尽责人人享有的| 1024视频免费在线观看| 水蜜桃什么品种好| 亚洲av欧美aⅴ国产| 亚洲国产精品成人久久小说| 永久免费av网站大全| 成年人午夜在线观看视频| 免费观看av网站的网址| 丰满少妇做爰视频| 国内毛片毛片毛片毛片毛片| 又黄又粗又硬又大视频| 欧美日韩成人在线一区二区| 亚洲少妇的诱惑av| 亚洲欧洲精品一区二区精品久久久| 国产亚洲欧美在线一区二区| 亚洲一区中文字幕在线| 男人操女人黄网站| 国产日韩欧美在线精品| 人人妻人人添人人爽欧美一区卜| 肉色欧美久久久久久久蜜桃| 亚洲精品在线美女| 国产精品亚洲av一区麻豆| 亚洲第一青青草原| 精品久久久精品久久久| 欧美日韩成人在线一区二区| av在线app专区| 精品一区在线观看国产| 午夜福利在线免费观看网站| 少妇 在线观看| 国产精品秋霞免费鲁丝片| 久久久精品免费免费高清| 真人做人爱边吃奶动态| 电影成人av| 久久久久久久大尺度免费视频| 色综合欧美亚洲国产小说| 深夜精品福利| 免费在线观看视频国产中文字幕亚洲 | 精品高清国产在线一区| 最近最新免费中文字幕在线| 国产精品影院久久| 久久午夜综合久久蜜桃| 日韩欧美免费精品| 一区二区三区激情视频| 爱豆传媒免费全集在线观看| 深夜精品福利| 在线 av 中文字幕| 久久香蕉激情| av欧美777| 91精品伊人久久大香线蕉| 久久天躁狠狠躁夜夜2o2o| 久久精品国产a三级三级三级| 久久久久精品国产欧美久久久 | av国产精品久久久久影院| 亚洲精品中文字幕一二三四区 | 性高湖久久久久久久久免费观看| 亚洲精品一卡2卡三卡4卡5卡 | 欧美日韩亚洲高清精品| 国产成人啪精品午夜网站| 久久久国产精品麻豆| 亚洲欧美日韩另类电影网站| 99热网站在线观看| 深夜精品福利| 久久久久网色| 午夜老司机福利片| kizo精华| 超碰成人久久| 日韩有码中文字幕| 热99久久久久精品小说推荐| 亚洲精品国产av成人精品| 制服诱惑二区| h视频一区二区三区| 丝袜美足系列| 91国产中文字幕| 久久国产精品男人的天堂亚洲| 窝窝影院91人妻| 电影成人av| 十八禁高潮呻吟视频| 国产精品一区二区免费欧美 | 欧美日韩黄片免| 大片电影免费在线观看免费| 国产免费现黄频在线看| 国产精品免费视频内射| 桃花免费在线播放| 久久久久精品人妻al黑| 精品福利观看| 女人精品久久久久毛片| 国产在线视频一区二区| 国产精品免费视频内射| 成人国产av品久久久| 国产成人精品无人区| 成年av动漫网址| 亚洲精品国产色婷婷电影| 99国产综合亚洲精品| 免费在线观看影片大全网站| 国产男女内射视频| 天天添夜夜摸| 婷婷成人精品国产| 美女扒开内裤让男人捅视频| 成人亚洲精品一区在线观看| 久久久精品区二区三区| 亚洲精品国产一区二区精华液| 女人久久www免费人成看片| 精品少妇一区二区三区视频日本电影| 18禁裸乳无遮挡动漫免费视频| 亚洲精品粉嫩美女一区| 一二三四在线观看免费中文在| 欧美黄色淫秽网站| 岛国毛片在线播放| 亚洲熟女毛片儿| 大片电影免费在线观看免费| 午夜福利在线免费观看网站| 不卡一级毛片| 丰满人妻熟妇乱又伦精品不卡| 肉色欧美久久久久久久蜜桃| 久久久久久久国产电影| 亚洲欧洲精品一区二区精品久久久| 久久人妻熟女aⅴ| 自线自在国产av| 淫妇啪啪啪对白视频 | 亚洲va日本ⅴa欧美va伊人久久 |