• <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
    国产成人一区二区三区免费视频网站| 精品人妻1区二区| 成年版毛片免费区| 亚洲 欧美 日韩 在线 免费| 精品第一国产精品| 男人的好看免费观看在线视频 | 亚洲av熟女| 18禁美女被吸乳视频| 久久中文字幕一级| 久久香蕉国产精品| 国产高清有码在线观看视频 | 99久久国产精品久久久| 亚洲欧美日韩高清专用| 久久九九热精品免费| 18禁美女被吸乳视频| a级毛片在线看网站| 无人区码免费观看不卡| 午夜激情福利司机影院| 中出人妻视频一区二区| 亚洲男人的天堂狠狠| 美女扒开内裤让男人捅视频| 级片在线观看| 国产成人av激情在线播放| 国产亚洲精品综合一区在线观看 | 国产精品,欧美在线| 亚洲精品美女久久av网站| 99国产综合亚洲精品| 免费搜索国产男女视频| 1024手机看黄色片| 亚洲av电影不卡..在线观看| 丁香六月欧美| 天天一区二区日本电影三级| 国语自产精品视频在线第100页| 亚洲成人国产一区在线观看| 国产黄a三级三级三级人| videosex国产| 亚洲aⅴ乱码一区二区在线播放 | 搞女人的毛片| 亚洲色图av天堂| 欧美绝顶高潮抽搐喷水| 一本大道久久a久久精品| 色综合婷婷激情| 在线观看日韩欧美| 久久精品夜夜夜夜夜久久蜜豆 | 人成视频在线观看免费观看| 精品欧美国产一区二区三| xxxwww97欧美| 久久国产精品人妻蜜桃| 亚洲av成人一区二区三| 757午夜福利合集在线观看| 又黄又爽又免费观看的视频| 中文亚洲av片在线观看爽| 又紧又爽又黄一区二区| 国产精品98久久久久久宅男小说| tocl精华| 成人一区二区视频在线观看| 亚洲成人精品中文字幕电影| 亚洲国产高清在线一区二区三| 一区二区三区国产精品乱码| 欧美日本视频| 三级男女做爰猛烈吃奶摸视频| 大型黄色视频在线免费观看| 国产亚洲精品第一综合不卡| 久久99热这里只有精品18| 黄色视频不卡| 日本 欧美在线| bbb黄色大片| 热99re8久久精品国产| 精品国产乱子伦一区二区三区| 在线观看www视频免费| 日本精品一区二区三区蜜桃| av福利片在线| 日本熟妇午夜| 午夜免费观看网址| 日本一区二区免费在线视频| 精品高清国产在线一区| 啦啦啦免费观看视频1| 亚洲aⅴ乱码一区二区在线播放 | 丁香欧美五月| 99久久无色码亚洲精品果冻| 亚洲在线自拍视频| 一区二区三区高清视频在线| 香蕉丝袜av| 麻豆成人午夜福利视频| 黑人操中国人逼视频| 国产野战对白在线观看| 俄罗斯特黄特色一大片| 国产成人啪精品午夜网站| 国产精品1区2区在线观看.| 午夜a级毛片| 中文字幕久久专区| 999精品在线视频| 哪里可以看免费的av片| 在线播放国产精品三级| 午夜视频精品福利| 欧美日韩亚洲综合一区二区三区_| 97超级碰碰碰精品色视频在线观看| 精品久久久久久久人妻蜜臀av| 一进一出抽搐动态| 国内精品久久久久精免费| 亚洲中文字幕日韩| av有码第一页| 成年女人毛片免费观看观看9| 欧美高清成人免费视频www| 全区人妻精品视频| 哪里可以看免费的av片| 桃色一区二区三区在线观看| 国产精品免费视频内射| 国产精品,欧美在线| 亚洲五月婷婷丁香| 国产激情欧美一区二区| 国产精品一区二区精品视频观看| 欧美+亚洲+日韩+国产| АⅤ资源中文在线天堂| 中亚洲国语对白在线视频| 欧美极品一区二区三区四区| 国产精品一区二区精品视频观看| 哪里可以看免费的av片| 在线看三级毛片| 欧美绝顶高潮抽搐喷水| 一进一出好大好爽视频| bbb黄色大片| 黑人操中国人逼视频| 亚洲乱码一区二区免费版| 一个人免费在线观看电影 | 真人一进一出gif抽搐免费| 亚洲国产中文字幕在线视频| 777久久人妻少妇嫩草av网站| 欧美黑人精品巨大| 成人手机av| 中文字幕久久专区| 久久欧美精品欧美久久欧美| 日日干狠狠操夜夜爽| 成熟少妇高潮喷水视频| 亚洲 欧美 日韩 在线 免费| 黄色视频,在线免费观看| 国产视频内射| 激情在线观看视频在线高清| 国产精品久久久人人做人人爽| 男插女下体视频免费在线播放| 亚洲 欧美一区二区三区| 国产三级中文精品| 啪啪无遮挡十八禁网站| 怎么达到女性高潮| 日本五十路高清| 在线观看美女被高潮喷水网站 | 欧美一区二区精品小视频在线| 淫秽高清视频在线观看| 一卡2卡三卡四卡精品乱码亚洲| 宅男免费午夜| 国产成人av激情在线播放| 亚洲专区中文字幕在线| а√天堂www在线а√下载| 在线国产一区二区在线| 免费搜索国产男女视频| 国产1区2区3区精品| 伊人久久大香线蕉亚洲五| 嫩草影院精品99| 亚洲av电影不卡..在线观看| 我要搜黄色片| 91大片在线观看| 国产精品一区二区精品视频观看| 老汉色∧v一级毛片| 久久久久久亚洲精品国产蜜桃av| 国产精品乱码一区二三区的特点| 亚洲精品美女久久久久99蜜臀| av天堂在线播放| 久久久水蜜桃国产精品网| 亚洲精品久久国产高清桃花| 老司机午夜十八禁免费视频| 国产1区2区3区精品| 亚洲真实伦在线观看| 欧美黑人欧美精品刺激| 一本久久中文字幕| 亚洲精品国产一区二区精华液| 亚洲自偷自拍图片 自拍| 久久久久亚洲av毛片大全| 视频区欧美日本亚洲| 可以在线观看的亚洲视频| 欧美一级a爱片免费观看看 | 亚洲国产欧美人成| 十八禁人妻一区二区| 成人18禁在线播放| 熟妇人妻久久中文字幕3abv| 真人一进一出gif抽搐免费| 女警被强在线播放| 日韩欧美一区二区三区在线观看| 国产一区在线观看成人免费| 免费看美女性在线毛片视频| 很黄的视频免费| 一级毛片女人18水好多| 国产麻豆成人av免费视频| 亚洲人与动物交配视频| 久久人妻av系列| 亚洲天堂国产精品一区在线| 精品熟女少妇八av免费久了| 国产激情欧美一区二区| 成人高潮视频无遮挡免费网站| 亚洲 国产 在线| 国产三级黄色录像| 1024手机看黄色片| 亚洲一卡2卡3卡4卡5卡精品中文| 少妇的丰满在线观看| 日日爽夜夜爽网站| 免费av毛片视频| 欧美成人性av电影在线观看| 夜夜夜夜夜久久久久| 免费一级毛片在线播放高清视频| 999久久久精品免费观看国产| 国产精品久久久久久精品电影| 欧美zozozo另类| 亚洲av电影不卡..在线观看| 亚洲色图av天堂| 免费av毛片视频| www.www免费av| 欧美国产日韩亚洲一区| 国产成人影院久久av| 久久伊人香网站| 亚洲最大成人中文| 日本五十路高清| 国产不卡一卡二| 国产高清视频在线观看网站| 日日干狠狠操夜夜爽| 免费看日本二区| 99国产精品一区二区三区| 2021天堂中文幕一二区在线观| 亚洲熟妇熟女久久| 制服诱惑二区| 啦啦啦观看免费观看视频高清| 色哟哟哟哟哟哟| 中文字幕人妻丝袜一区二区| 欧美成狂野欧美在线观看| 亚洲欧洲精品一区二区精品久久久| 真人做人爱边吃奶动态| 国产高清视频在线播放一区| 免费高清视频大片| 麻豆国产av国片精品| 国产精品电影一区二区三区| 亚洲一区二区三区色噜噜| 深夜精品福利| 香蕉国产在线看| 可以在线观看毛片的网站| 好男人电影高清在线观看| 一本一本综合久久| 亚洲在线自拍视频| 99久久久亚洲精品蜜臀av| 美女黄网站色视频| 成人亚洲精品av一区二区| 国产精品综合久久久久久久免费| 免费看十八禁软件| 欧美日韩国产亚洲二区| 欧美性猛交黑人性爽| 日韩欧美国产在线观看| 国产高清激情床上av| 午夜激情福利司机影院| 亚洲成人久久性| 99国产精品一区二区蜜桃av| 香蕉久久夜色| 90打野战视频偷拍视频| 国产精品野战在线观看| 长腿黑丝高跟| 亚洲人成网站高清观看| 国产欧美日韩一区二区三| 又爽又黄无遮挡网站| 99久久久亚洲精品蜜臀av| 国产精品久久视频播放| 亚洲av五月六月丁香网| 51午夜福利影视在线观看| 亚洲精品av麻豆狂野| 久久这里只有精品中国| 国内精品久久久久久久电影| 好男人电影高清在线观看| 午夜免费激情av| 国产精品国产高清国产av| 天天添夜夜摸| 日韩欧美在线乱码| 国产精品一区二区三区四区久久| 日韩欧美在线二视频| 国产野战对白在线观看| xxxwww97欧美| 国产片内射在线| 欧美不卡视频在线免费观看 | 国产区一区二久久| 精品电影一区二区在线| 好男人在线观看高清免费视频| 欧美不卡视频在线免费观看 | av福利片在线观看| а√天堂www在线а√下载| 国产成人精品久久二区二区免费| 亚洲精品国产精品久久久不卡| 国产又黄又爽又无遮挡在线| 日韩 欧美 亚洲 中文字幕| 日本撒尿小便嘘嘘汇集6| 一区二区三区高清视频在线| 久久久国产成人免费| 亚洲人成网站在线播放欧美日韩| 一个人观看的视频www高清免费观看 | 亚洲中文av在线| 欧美大码av| 男女下面进入的视频免费午夜| 久久久久久亚洲精品国产蜜桃av| 2021天堂中文幕一二区在线观| 久99久视频精品免费| 午夜两性在线视频| 一边摸一边抽搐一进一小说| 午夜激情av网站| 91av网站免费观看| 国产aⅴ精品一区二区三区波| 999精品在线视频| 老鸭窝网址在线观看| 国产伦一二天堂av在线观看| 国产一区二区三区在线臀色熟女| 淫妇啪啪啪对白视频| 久久久久久九九精品二区国产 | 国产精品久久久久久久电影 | 成人三级黄色视频| 国产单亲对白刺激| 美女大奶头视频| 欧美色欧美亚洲另类二区| 久久久久久人人人人人| 精品欧美国产一区二区三| 一边摸一边做爽爽视频免费| 国产成人影院久久av| 777久久人妻少妇嫩草av网站| 一夜夜www| 中文资源天堂在线| 国产精品av久久久久免费| 国产成人精品久久二区二区免费| 淫秽高清视频在线观看| 国产亚洲精品综合一区在线观看 | 成人高潮视频无遮挡免费网站| av在线天堂中文字幕| 亚洲精品av麻豆狂野| 国产精品一区二区三区四区久久| 欧美一级a爱片免费观看看 | 一本大道久久a久久精品| 色综合站精品国产| 亚洲av片天天在线观看| 国产一区二区在线av高清观看| 久久人妻av系列| 丰满的人妻完整版| 最新美女视频免费是黄的| 国内揄拍国产精品人妻在线| 国产精品 欧美亚洲| 99re在线观看精品视频| 亚洲色图av天堂| 亚洲av第一区精品v没综合| 欧美日韩国产亚洲二区| 黄频高清免费视频| www日本在线高清视频| 欧美性猛交╳xxx乱大交人| 日日摸夜夜添夜夜添小说| 久久久水蜜桃国产精品网| 国内精品久久久久精免费| 国产成年人精品一区二区| 级片在线观看| 免费电影在线观看免费观看| xxx96com| 久久精品影院6| 国内精品久久久久精免费| 男女床上黄色一级片免费看| 久久人妻av系列| 国产精品,欧美在线| 一本精品99久久精品77| 宅男免费午夜| 精品国产乱子伦一区二区三区| 欧美黑人欧美精品刺激| 日日夜夜操网爽| www.999成人在线观看| 黄色a级毛片大全视频| 国产精品,欧美在线| 最近最新中文字幕大全电影3| 99久久精品热视频| 国产高清videossex| 免费av毛片视频| 正在播放国产对白刺激| 美女黄网站色视频| 国产探花在线观看一区二区| 久久久久久国产a免费观看| 免费一级毛片在线播放高清视频| 亚洲人成77777在线视频| 一区二区三区国产精品乱码| av在线天堂中文字幕| 国产亚洲精品一区二区www| 久久中文字幕一级| 日本熟妇午夜| 精品久久久久久久末码| 国产精品1区2区在线观看.| av超薄肉色丝袜交足视频| 黄色成人免费大全| 青草久久国产| 悠悠久久av| bbb黄色大片| 午夜福利在线在线| 久久精品国产99精品国产亚洲性色| 国产av在哪里看| aaaaa片日本免费| 午夜免费成人在线视频| 精品欧美一区二区三区在线| 久久精品国产亚洲av香蕉五月| 亚洲精品国产一区二区精华液| 亚洲片人在线观看| 亚洲欧洲精品一区二区精品久久久| 宅男免费午夜| 成人一区二区视频在线观看| 中文字幕熟女人妻在线| 欧美乱色亚洲激情| 一级片免费观看大全| 在线a可以看的网站| 日本一二三区视频观看| 不卡一级毛片| 中文字幕精品亚洲无线码一区| 欧美+亚洲+日韩+国产| 国产视频内射| 国产真人三级小视频在线观看| 欧美中文综合在线视频| 88av欧美| 午夜日韩欧美国产| 两个人视频免费观看高清| 国语自产精品视频在线第100页| 久久精品成人免费网站| 人妻夜夜爽99麻豆av| 久久久国产成人免费| 又粗又爽又猛毛片免费看| 两人在一起打扑克的视频| 久久久国产欧美日韩av| 色综合欧美亚洲国产小说| 久久99热这里只有精品18| 久久热在线av| 男男h啪啪无遮挡| 舔av片在线| 日韩国内少妇激情av| 男女视频在线观看网站免费 | 亚洲人成伊人成综合网2020| 久久人人精品亚洲av| 色av中文字幕| 久久99热这里只有精品18| 黑人欧美特级aaaaaa片| www.999成人在线观看| 中文字幕高清在线视频| 欧美成人性av电影在线观看| 熟女电影av网| 黄色片一级片一级黄色片| 村上凉子中文字幕在线| 成人亚洲精品av一区二区| 欧美午夜高清在线| 日日摸夜夜添夜夜添小说| 成在线人永久免费视频| 成年版毛片免费区| tocl精华| av天堂在线播放| 亚洲av电影不卡..在线观看| 美女高潮喷水抽搐中文字幕| 日韩欧美在线乱码| 人成视频在线观看免费观看| 999久久久精品免费观看国产| 国产麻豆成人av免费视频| 欧美av亚洲av综合av国产av| 日韩欧美国产在线观看| 亚洲国产欧美人成| 亚洲av电影不卡..在线观看| 国内少妇人妻偷人精品xxx网站 | 热99re8久久精品国产| 日本免费一区二区三区高清不卡| 最新在线观看一区二区三区| 久久99热这里只有精品18| 怎么达到女性高潮| 在线观看美女被高潮喷水网站 | 免费观看人在逋| av免费在线观看网站| 好男人在线观看高清免费视频| 国产精品久久久久久亚洲av鲁大| 欧美日韩亚洲国产一区二区在线观看| 88av欧美| 香蕉久久夜色| 欧美成人性av电影在线观看| 在线十欧美十亚洲十日本专区| 欧美性猛交╳xxx乱大交人| 国产精品一区二区三区四区免费观看 | 老司机在亚洲福利影院| 香蕉国产在线看| av天堂在线播放| 美女扒开内裤让男人捅视频| 一二三四社区在线视频社区8| 久久久久久久久中文| av视频在线观看入口| 欧美三级亚洲精品| 夜夜夜夜夜久久久久| 国产精品av视频在线免费观看| 高潮久久久久久久久久久不卡| 亚洲人成77777在线视频| 曰老女人黄片| 最近在线观看免费完整版| 国产野战对白在线观看| 中文在线观看免费www的网站 | 亚洲va日本ⅴa欧美va伊人久久| 免费看十八禁软件| 波多野结衣高清无吗| 国产激情久久老熟女| 十八禁网站免费在线| 国产精品亚洲美女久久久| 国模一区二区三区四区视频 | 九九热线精品视视频播放| 欧美激情久久久久久爽电影| 国产av又大| 精品久久久久久,| 制服人妻中文乱码| 国产精品自产拍在线观看55亚洲| 在线观看66精品国产| 亚洲人与动物交配视频| 午夜a级毛片| 成人一区二区视频在线观看| 岛国在线免费视频观看| 黄色 视频免费看| 久久久久久亚洲精品国产蜜桃av| 亚洲成人国产一区在线观看| 亚洲国产日韩欧美精品在线观看 | 日韩欧美免费精品| 国产精品1区2区在线观看.| 可以在线观看的亚洲视频| 国产精品香港三级国产av潘金莲| 国产探花在线观看一区二区| 岛国在线免费视频观看| 亚洲成人久久性| 特大巨黑吊av在线直播| 国产成年人精品一区二区| 国产成+人综合+亚洲专区| 午夜免费激情av| 最好的美女福利视频网| 久久人妻av系列| 欧美日韩亚洲综合一区二区三区_| 国产亚洲av高清不卡| 久久久久性生活片| 国产精品久久久人人做人人爽| 成人精品一区二区免费| 国产精品久久视频播放| 国产精品影院久久| 国内久久婷婷六月综合欲色啪| 久久久久久国产a免费观看| 亚洲av五月六月丁香网| 国产精品美女特级片免费视频播放器 | 午夜免费观看网址| 韩国av一区二区三区四区| 18禁裸乳无遮挡免费网站照片| av中文乱码字幕在线| 91大片在线观看| 久久人妻福利社区极品人妻图片| 日本黄色视频三级网站网址| e午夜精品久久久久久久| 欧美午夜高清在线| 国产午夜精品论理片| www.www免费av| 女人高潮潮喷娇喘18禁视频| 精品第一国产精品| 精品免费久久久久久久清纯| 我的老师免费观看完整版| 人妻久久中文字幕网| 久久香蕉国产精品| 黄色成人免费大全| 两人在一起打扑克的视频| 日韩欧美三级三区| 欧美绝顶高潮抽搐喷水| 国产熟女午夜一区二区三区| 亚洲成人中文字幕在线播放| 窝窝影院91人妻| av超薄肉色丝袜交足视频| 人人妻,人人澡人人爽秒播| 级片在线观看| 欧美日韩乱码在线| 日韩欧美 国产精品| 国产成人啪精品午夜网站| 露出奶头的视频| 日韩中文字幕欧美一区二区| 首页视频小说图片口味搜索| 亚洲国产欧洲综合997久久,| 天天一区二区日本电影三级| 日韩有码中文字幕| 久久天堂一区二区三区四区| 一本大道久久a久久精品| 国产熟女午夜一区二区三区| 亚洲美女视频黄频| 人人妻人人看人人澡| 午夜免费激情av| 巨乳人妻的诱惑在线观看| 日日爽夜夜爽网站| 免费看十八禁软件| 亚洲欧美日韩无卡精品| 国产精品电影一区二区三区| 精品久久久久久久末码| 亚洲乱码一区二区免费版| 啦啦啦观看免费观看视频高清| 啪啪无遮挡十八禁网站| 99国产极品粉嫩在线观看| www.精华液| 日韩三级视频一区二区三区| 亚洲成av人片在线播放无| 国产精品日韩av在线免费观看| 18美女黄网站色大片免费观看| 全区人妻精品视频| 久久香蕉激情| 日本五十路高清| 欧美另类亚洲清纯唯美| 免费在线观看完整版高清| 免费一级毛片在线播放高清视频| 黄色视频,在线免费观看| 国产精品一区二区三区四区免费观看 | 啪啪无遮挡十八禁网站| 精品国产超薄肉色丝袜足j| 精品一区二区三区四区五区乱码| 啦啦啦免费观看视频1| 久久人妻av系列| 露出奶头的视频| 免费无遮挡裸体视频| 男插女下体视频免费在线播放| 成人国产一区最新在线观看|