• <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
    日本精品一区二区三区蜜桃| 不卡一级毛片| 热99国产精品久久久久久7| 人妻丰满熟妇av一区二区三区| 精品日产1卡2卡| 性少妇av在线| 中文欧美无线码| 久久性视频一级片| 大型av网站在线播放| 国产精品自产拍在线观看55亚洲| www.999成人在线观看| 午夜a级毛片| 最好的美女福利视频网| 亚洲全国av大片| 12—13女人毛片做爰片一| 成人特级黄色片久久久久久久| 国产精品99久久99久久久不卡| 两性午夜刺激爽爽歪歪视频在线观看 | 成人三级黄色视频| 久久久久久免费高清国产稀缺| 亚洲欧美日韩另类电影网站| 欧美激情久久久久久爽电影 | 女人爽到高潮嗷嗷叫在线视频| 中文字幕av电影在线播放| 日日夜夜操网爽| 人人澡人人妻人| 看免费av毛片| 看免费av毛片| 久久青草综合色| 91精品国产国语对白视频| 91在线观看av| 一级毛片女人18水好多| 精品久久久久久久久久免费视频 | 精品少妇一区二区三区视频日本电影| 中文亚洲av片在线观看爽| 又紧又爽又黄一区二区| 黄色成人免费大全| 精品少妇一区二区三区视频日本电影| 午夜激情av网站| 9热在线视频观看99| 亚洲午夜精品一区,二区,三区| 成年版毛片免费区| 欧美人与性动交α欧美软件| 99国产精品一区二区三区| 一本综合久久免费| 91九色精品人成在线观看| 精品国产国语对白av| 久久久久久久久中文| 在线观看舔阴道视频| 国产亚洲精品综合一区在线观看 | 国产亚洲精品综合一区在线观看 | 无限看片的www在线观看| 色播在线永久视频| 91九色精品人成在线观看| 在线十欧美十亚洲十日本专区| 老熟妇乱子伦视频在线观看| 老司机在亚洲福利影院| 九色亚洲精品在线播放| 美女国产高潮福利片在线看| 欧美日韩av久久| 欧美日韩精品网址| 又紧又爽又黄一区二区| 亚洲中文字幕日韩| 亚洲精品国产区一区二| 波多野结衣av一区二区av| 国产精品影院久久| 亚洲全国av大片| 国产三级在线视频| 日本欧美视频一区| 亚洲成人精品中文字幕电影 | 久久中文字幕一级| 国产亚洲精品久久久久5区| 精品福利永久在线观看| 超色免费av| 欧美日韩视频精品一区| 久热爱精品视频在线9| 久热爱精品视频在线9| 悠悠久久av| 久久久国产欧美日韩av| 日本欧美视频一区| 日韩欧美一区二区三区在线观看| 亚洲,欧美精品.| 一二三四社区在线视频社区8| 亚洲性夜色夜夜综合| 日韩av在线大香蕉| 制服诱惑二区| 亚洲一卡2卡3卡4卡5卡精品中文| 无人区码免费观看不卡| 两人在一起打扑克的视频| √禁漫天堂资源中文www| xxx96com| 国产99久久九九免费精品| 国产成人精品久久二区二区免费| 高清欧美精品videossex| 水蜜桃什么品种好| 天堂俺去俺来也www色官网| 999久久久精品免费观看国产| 成人av一区二区三区在线看| 女人高潮潮喷娇喘18禁视频| 免费少妇av软件| 丰满的人妻完整版| 女性生殖器流出的白浆| 日韩欧美免费精品| 十八禁人妻一区二区| 久久香蕉激情| av片东京热男人的天堂| 亚洲欧美精品综合久久99| 国产av一区在线观看免费| 欧美最黄视频在线播放免费 | 国产av在哪里看| 精品国产一区二区三区四区第35| 国产精品久久久av美女十八| 亚洲,欧美精品.| 天天添夜夜摸| 99国产精品一区二区蜜桃av| 视频在线观看一区二区三区| 99热国产这里只有精品6| 亚洲欧美精品综合久久99| 黄片播放在线免费| 99久久99久久久精品蜜桃| 国产精品国产高清国产av| 激情视频va一区二区三区| 日韩欧美国产一区二区入口| avwww免费| 国产精品九九99| 一夜夜www| 中文字幕人妻熟女乱码| 久久久久国产一级毛片高清牌| 久久热在线av| 老司机深夜福利视频在线观看| 夜夜看夜夜爽夜夜摸 | 91精品三级在线观看| 亚洲av电影在线进入| 电影成人av| 国产成人av教育| 国产精品永久免费网站| 男男h啪啪无遮挡| 亚洲专区字幕在线| 一边摸一边做爽爽视频免费| 亚洲伊人色综图| 久久天堂一区二区三区四区| 在线天堂中文资源库| 中亚洲国语对白在线视频| 久9热在线精品视频| 国产单亲对白刺激| 成人18禁在线播放| 一区二区三区国产精品乱码| 曰老女人黄片| 91字幕亚洲| 国产精品久久久人人做人人爽| 淫秽高清视频在线观看| 亚洲精华国产精华精| 精品国产一区二区久久| 成人18禁在线播放| 国产精品野战在线观看 | 精品国产乱码久久久久久男人| 欧美成狂野欧美在线观看| 中出人妻视频一区二区| 国产成人免费无遮挡视频| 国产乱人伦免费视频| 黑人巨大精品欧美一区二区蜜桃| 在线观看一区二区三区| 长腿黑丝高跟| a级毛片在线看网站| 国产亚洲欧美98| 怎么达到女性高潮| 曰老女人黄片| 日韩精品青青久久久久久| 国产av一区二区精品久久| 午夜亚洲福利在线播放| 麻豆国产av国片精品| 亚洲欧洲精品一区二区精品久久久| 人人澡人人妻人| 免费日韩欧美在线观看| 妹子高潮喷水视频| 啦啦啦免费观看视频1| 亚洲av电影在线进入| 一区二区三区精品91| 成人亚洲精品av一区二区 | 亚洲视频免费观看视频| 一个人免费在线观看的高清视频| av网站免费在线观看视频| 国产精品秋霞免费鲁丝片| 欧美日韩一级在线毛片| 两人在一起打扑克的视频| 欧美最黄视频在线播放免费 | 一区福利在线观看| 国产三级黄色录像| 这个男人来自地球电影免费观看| 精品午夜福利视频在线观看一区| 黄色毛片三级朝国网站| 亚洲精品一卡2卡三卡4卡5卡| 久久欧美精品欧美久久欧美| 亚洲片人在线观看| 五月开心婷婷网| 18美女黄网站色大片免费观看| 国产97色在线日韩免费| 丰满迷人的少妇在线观看| 高清欧美精品videossex| 成人国产一区最新在线观看| 日韩欧美一区二区三区在线观看| 俄罗斯特黄特色一大片| 一级作爱视频免费观看| 91老司机精品| 国产精品久久久久成人av| 99久久人妻综合| 亚洲成av片中文字幕在线观看| 久久久久久久久久久久大奶| 欧美丝袜亚洲另类 | 麻豆成人av在线观看| 国产在线观看jvid| 精品久久久久久久毛片微露脸| www.www免费av| 狂野欧美激情性xxxx| 美女福利国产在线| 中文亚洲av片在线观看爽| 啦啦啦免费观看视频1| 一个人免费在线观看的高清视频| 欧美色视频一区免费| 欧洲精品卡2卡3卡4卡5卡区| 成人永久免费在线观看视频| 身体一侧抽搐| 国产高清激情床上av| 日本免费a在线| 免费在线观看亚洲国产| 我的亚洲天堂| 久久香蕉精品热| 亚洲国产毛片av蜜桃av| 美女高潮到喷水免费观看| 欧美色视频一区免费| 国产高清videossex| 一边摸一边抽搐一进一出视频| 亚洲成人免费av在线播放| 欧美日韩中文字幕国产精品一区二区三区 | 一级毛片精品| 身体一侧抽搐| 亚洲av美国av| 电影成人av| 欧美丝袜亚洲另类 | 日韩成人在线观看一区二区三区| 久久中文看片网| 一级,二级,三级黄色视频| 中文字幕色久视频| netflix在线观看网站| 国产成人欧美| 黄色女人牲交| 亚洲七黄色美女视频| 久久久国产成人精品二区 | 日本精品一区二区三区蜜桃| 1024视频免费在线观看| 精品久久久久久成人av| 激情在线观看视频在线高清| 欧美在线黄色| 国产亚洲精品久久久久5区| 美女 人体艺术 gogo| 一级毛片精品| 每晚都被弄得嗷嗷叫到高潮| 宅男免费午夜| 亚洲五月天丁香| 久久久国产一区二区| 欧美性长视频在线观看| 久久人人精品亚洲av| 国产一区二区三区在线臀色熟女 | 热re99久久精品国产66热6| 少妇被粗大的猛进出69影院| 电影成人av| 成人18禁高潮啪啪吃奶动态图| 水蜜桃什么品种好| 国产主播在线观看一区二区| 亚洲一码二码三码区别大吗| 女人被狂操c到高潮| 黄色毛片三级朝国网站| 黄色 视频免费看| 亚洲男人天堂网一区| 中文欧美无线码| 黄色怎么调成土黄色| 亚洲国产欧美日韩在线播放| 免费在线观看完整版高清| 一区二区三区精品91| 制服人妻中文乱码| 少妇 在线观看| 在线永久观看黄色视频| 嫁个100分男人电影在线观看| xxx96com| 在线观看免费视频日本深夜| 在线观看免费视频网站a站| 高清黄色对白视频在线免费看| 精品久久久久久电影网| 99国产精品99久久久久| 亚洲成人久久性| 女生性感内裤真人,穿戴方法视频| 午夜精品久久久久久毛片777| 50天的宝宝边吃奶边哭怎么回事| 欧美日韩黄片免| 欧美日本中文国产一区发布| 欧美色视频一区免费| 国产一区二区三区视频了| 国产主播在线观看一区二区| 亚洲一区二区三区欧美精品| 国产片内射在线| 中文亚洲av片在线观看爽| 深夜精品福利| 桃红色精品国产亚洲av| 精品第一国产精品| 久久精品91蜜桃| 亚洲熟妇中文字幕五十中出 | 亚洲熟女毛片儿| 久久精品国产清高在天天线| 久久九九热精品免费| 久久青草综合色| 欧美日韩亚洲高清精品| 国产一区二区在线av高清观看| 国产精品美女特级片免费视频播放器 | 亚洲美女黄片视频| 欧美日韩亚洲综合一区二区三区_| 免费不卡黄色视频| 一本综合久久免费| 无遮挡黄片免费观看| 少妇裸体淫交视频免费看高清 | tocl精华| 精品一区二区三区av网在线观看| 每晚都被弄得嗷嗷叫到高潮| 久久精品人人爽人人爽视色| av福利片在线| 看黄色毛片网站| 亚洲美女黄片视频| 欧美成人性av电影在线观看| 老司机午夜福利在线观看视频| 日韩精品中文字幕看吧| 男女做爰动态图高潮gif福利片 | 黄片小视频在线播放| 91国产中文字幕| 久久人人爽av亚洲精品天堂| 欧美日韩亚洲高清精品| 亚洲av日韩精品久久久久久密| 国产成人欧美在线观看| 欧美日韩精品网址| 黄片小视频在线播放| 一区在线观看完整版| 国产成人免费无遮挡视频| 韩国精品一区二区三区| 在线天堂中文资源库| 日本vs欧美在线观看视频| 两性午夜刺激爽爽歪歪视频在线观看 | 久久国产亚洲av麻豆专区| 91成人精品电影| 亚洲七黄色美女视频| 亚洲精品国产一区二区精华液| 精品少妇一区二区三区视频日本电影| 1024香蕉在线观看| 欧美成狂野欧美在线观看| 丁香六月欧美| 国产av在哪里看| av网站免费在线观看视频| 中文字幕最新亚洲高清| 亚洲成人免费电影在线观看| 亚洲精品美女久久久久99蜜臀| 视频区图区小说| www日本在线高清视频| 80岁老熟妇乱子伦牲交| 99riav亚洲国产免费| 麻豆一二三区av精品| 免费看a级黄色片| 黑丝袜美女国产一区| 亚洲一区高清亚洲精品| av电影中文网址| 18禁国产床啪视频网站| 99香蕉大伊视频| 看片在线看免费视频| 亚洲美女黄片视频| 黄频高清免费视频| 夫妻午夜视频| 色哟哟哟哟哟哟| 国产三级黄色录像| www国产在线视频色| 国产精品一区二区精品视频观看| 中文字幕人妻丝袜一区二区| 国产精品秋霞免费鲁丝片| 亚洲国产欧美一区二区综合| 成人三级做爰电影| 激情视频va一区二区三区| 亚洲,欧美精品.| 三级毛片av免费| 日韩免费高清中文字幕av| 国产又色又爽无遮挡免费看| 18禁裸乳无遮挡免费网站照片 | 一二三四在线观看免费中文在| 男女床上黄色一级片免费看| 女警被强在线播放| 久久久久久久久中文| 国产不卡一卡二| 国产精华一区二区三区| 国产熟女xx| 精品国产一区二区三区四区第35| 视频在线观看一区二区三区| 精品一区二区三区av网在线观看| 久久人人97超碰香蕉20202| 黄频高清免费视频| 国产人伦9x9x在线观看| 999精品在线视频| 免费观看人在逋| 精品一区二区三区av网在线观看| 精品国产一区二区三区四区第35| 成人手机av| 日韩免费高清中文字幕av| 黄片小视频在线播放| 极品人妻少妇av视频| 九色亚洲精品在线播放| 成年人免费黄色播放视频| 91麻豆av在线| 亚洲一区二区三区色噜噜 | 久久久精品国产亚洲av高清涩受| 成年版毛片免费区| 亚洲精品粉嫩美女一区| 欧美另类亚洲清纯唯美| 大型黄色视频在线免费观看| 久久精品91无色码中文字幕| 国产成人精品久久二区二区免费| 日韩中文字幕欧美一区二区| 亚洲av成人av| 中文字幕色久视频| 麻豆久久精品国产亚洲av | 99国产精品99久久久久| 日本 av在线| 在线观看免费午夜福利视频| 亚洲一区中文字幕在线| www.999成人在线观看| 婷婷丁香在线五月| 老熟妇仑乱视频hdxx| 最近最新免费中文字幕在线| 午夜日韩欧美国产| 国产精品国产av在线观看| aaaaa片日本免费| 亚洲中文字幕日韩| 精品免费久久久久久久清纯| 国产99久久九九免费精品| 日韩欧美一区视频在线观看| 18禁黄网站禁片午夜丰满| 欧美精品亚洲一区二区| 免费高清在线观看日韩| 亚洲欧美日韩高清在线视频| 天天影视国产精品| 99香蕉大伊视频| 少妇 在线观看| 香蕉丝袜av| 午夜免费鲁丝| 日韩 欧美 亚洲 中文字幕| 国产精品美女特级片免费视频播放器 | 国产黄a三级三级三级人| 熟女少妇亚洲综合色aaa.| 日韩欧美国产一区二区入口| 又大又爽又粗| 成人永久免费在线观看视频| 长腿黑丝高跟| 熟女少妇亚洲综合色aaa.| 乱人伦中国视频| 亚洲五月婷婷丁香| 不卡av一区二区三区| 久久精品人人爽人人爽视色| 国产精品一区二区在线不卡| 无遮挡黄片免费观看| 黄色片一级片一级黄色片| 夜夜看夜夜爽夜夜摸 | 在线av久久热| 妹子高潮喷水视频| av超薄肉色丝袜交足视频| 国产亚洲精品一区二区www| 长腿黑丝高跟| 中文字幕最新亚洲高清| 亚洲精品粉嫩美女一区| 欧美丝袜亚洲另类 | 色综合婷婷激情| 每晚都被弄得嗷嗷叫到高潮| 中文字幕人妻丝袜一区二区| 久久精品亚洲av国产电影网| 在线播放国产精品三级| 美国免费a级毛片| 在线观看一区二区三区激情| 别揉我奶头~嗯~啊~动态视频| 手机成人av网站| 日韩视频一区二区在线观看| 午夜免费鲁丝| 欧美午夜高清在线| 老司机福利观看| 新久久久久国产一级毛片| 久久午夜亚洲精品久久| 麻豆一二三区av精品| 欧美av亚洲av综合av国产av| 久9热在线精品视频| 9191精品国产免费久久| 国产真人三级小视频在线观看| 日本一区二区免费在线视频| 欧美精品亚洲一区二区| 日韩欧美免费精品| 一级毛片高清免费大全| 国产精品二区激情视频| 国产高清国产精品国产三级| 男人舔女人下体高潮全视频| 成人国语在线视频| 18禁观看日本| 国产欧美日韩一区二区精品| 国产视频一区二区在线看| 亚洲一码二码三码区别大吗| 好看av亚洲va欧美ⅴa在| 午夜福利免费观看在线| 久久国产亚洲av麻豆专区| 亚洲,欧美精品.| 亚洲精品国产色婷婷电影| 男女下面进入的视频免费午夜 | 搡老熟女国产l中国老女人| 啦啦啦 在线观看视频| 久久午夜综合久久蜜桃| av欧美777| 国产单亲对白刺激| 女人被躁到高潮嗷嗷叫费观| 青草久久国产| 国产亚洲精品久久久久久毛片| 一个人观看的视频www高清免费观看 | 岛国在线观看网站| 丝袜美腿诱惑在线| 黑人操中国人逼视频| 五月开心婷婷网| 国产精品电影一区二区三区| 99精品欧美一区二区三区四区| 在线播放国产精品三级| 亚洲国产看品久久| 可以免费在线观看a视频的电影网站| 婷婷精品国产亚洲av在线| 99久久国产精品久久久| 每晚都被弄得嗷嗷叫到高潮| 后天国语完整版免费观看| 欧美日韩黄片免| 十八禁人妻一区二区| 国产视频一区二区在线看| 视频在线观看一区二区三区| 9热在线视频观看99| 国产精品久久久久久人妻精品电影| 黑人操中国人逼视频| 日韩视频一区二区在线观看| 在线观看免费高清a一片| 久久香蕉国产精品| 欧美最黄视频在线播放免费 | tocl精华| 悠悠久久av| 别揉我奶头~嗯~啊~动态视频| 日韩中文字幕欧美一区二区| 老鸭窝网址在线观看| 欧美不卡视频在线免费观看 | 真人一进一出gif抽搐免费| 动漫黄色视频在线观看| 国产色视频综合| 精品一区二区三区视频在线观看免费 | 国产一区二区三区在线臀色熟女 | 国产亚洲欧美在线一区二区| 黄色怎么调成土黄色| 国产一区二区三区综合在线观看| 国产成人av激情在线播放| 亚洲男人天堂网一区| 久久热在线av| 岛国视频午夜一区免费看| 久久人人爽av亚洲精品天堂| 亚洲av美国av| 怎么达到女性高潮| 国产黄a三级三级三级人| 黄色a级毛片大全视频| 亚洲成av片中文字幕在线观看| 91老司机精品| 免费高清视频大片| 91精品国产国语对白视频| 美女高潮到喷水免费观看| 91精品国产国语对白视频| 琪琪午夜伦伦电影理论片6080| 亚洲av成人一区二区三| 成在线人永久免费视频| 一本大道久久a久久精品| 男女下面插进去视频免费观看| www日本在线高清视频| 交换朋友夫妻互换小说| 午夜激情av网站| 十八禁人妻一区二区| 黄色 视频免费看| 视频区图区小说| 最好的美女福利视频网| 中文字幕色久视频| 在线观看一区二区三区| 伦理电影免费视频| www.999成人在线观看| 一级毛片精品| 日韩视频一区二区在线观看| 又黄又粗又硬又大视频| 亚洲一码二码三码区别大吗| 亚洲色图av天堂| 国产成人影院久久av| 欧美亚洲日本最大视频资源| 18禁黄网站禁片午夜丰满| 色在线成人网| 熟女少妇亚洲综合色aaa.| 国产一卡二卡三卡精品| 免费观看精品视频网站| 女性生殖器流出的白浆| 搡老熟女国产l中国老女人| 天堂中文最新版在线下载| 国产精品98久久久久久宅男小说| 大香蕉久久成人网| 男人舔女人的私密视频| 久久精品国产99精品国产亚洲性色 | 精品久久久久久电影网| 欧美av亚洲av综合av国产av| 成人亚洲精品av一区二区 | 亚洲精品国产一区二区精华液| 久久中文看片网| 丁香六月欧美| 黄频高清免费视频| 国产乱人伦免费视频| 十八禁人妻一区二区|