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

    Phase behaviors of ionic liquids attributed to the dual ionic and organic nature

    2022-10-22 08:16:34ChenyuTang唐晨宇andYantingWang王延颋
    Communications in Theoretical Physics 2022年9期

    Chenyu Tang (唐晨宇) and Yanting Wang (王延颋)

    1 CAS Key Laboratory of Theoretical Physics,Institute of Theoretical Physics,Chinese Academy of Sciences,Beijing 100190,China

    2 School of Physical Sciences,University of Chinese Academy of Sciences,Beijing 100049,China

    Abstract Ionic liquids (ILs),also known as room-temperature molten salts,are solely composed of ions with melting points usually below 100 °C.Because of their low volatility and vast amounts of species,ILs can serve as ‘green solvents’ and ‘designer solvents’ to meet the requirements of various applications by fine-tuning their molecular structures.A good understanding of the phase behaviors of ILs is certainly fundamentally important in terms of their wide applications.This review intends to summarize the major conclusions so far drawn on phase behaviors of ILs by computational,theoretical,and experimental studies,illustrating the intrinsic relationship between their dual ionic and organic nature and the crystalline phases,nanoscale segregation liquid phase,IL crystal phases,as well as phase behaviors of their mixture with small organic molecules.

    Keywords: ionic liquids,phase behaviors,nanoscale segregation liquid,ionic liquid crystal

    1.Introduction

    Ionic liquids(ILs)are a type of salts with low melting points,often below 100 °C,meaning that they tend to remain in the liquid phase at room temperature and are believed to exhibit some unique features because of the strong electrostatic interactions among ions.Typical aprotic ILs are normally composed of small anions and bulky cations with a long alkyl side-chain and a charged head group,as shown in figure 1,which demonstrates the chemical structure of 1-butyl-3-methylimidazolium chloride,a typical imidazolium-based IL.Possessing both ionic and organic features,they are believed to have advantageous properties of both organic liquids and inorganic salts,such as good solvation ability and tunability,low melting temperature,good conductivity,wide electrochemical window,thermal and electrochemical stabilities,non-volatility,and non-flammability [1–3].They are thus regarded as ‘green’ and ‘engineer’ solvents that can be utilized under many industrial circumstances [4–10].

    Understanding the fundamental properties of ILs,particularly their phase behaviors,is essential to their applications.To investigate their phase behaviors,many computational and experimental methods have been employed to investigate phase behaviors of ILs.Molecular dynamics(MD)simulation has become an important means of studying structures and dynamics of ILs [11–14]where different modelling methods that employ various software packages including GROMACS,NAMD,LAMMPS,etc [15–17]have been developed.All-atom force fields are commonly used in addressing IL-related problems by means of MD simulation[18–20],and the applicability of some commonly used all-atom force fields,including the Amber force field[21],OPLS force field[22],CHARMM force field[23],etc,to IL systems,has been verified and the models have been constantly improved to be better applied to ILs [14,24–34].By considering the polarizable effect at the atomic level,a polarizable model has also been developed to better quantify the microscopic structures and dynamics of ILs [35–37].Another modelling strategy is to apply a coarse-grained (CG) model in MD simulation,which reduces the cost of computation and renders researchers with longer simulated times than all-atom models.One of the CG methods that is used in the context of ILs is the Multiscale Corse-Graining (MS-CG) method [38–40],which matches the instantaneous forces applied to atoms in the allatom MD simulation to determine the optimal empirical CG forces between CG sites (atomic groups).By contrast,the effective force coarse-graining (EF-CG) method directly calculates the effective averaged force between each pair of CG sites (atomic groups) to gain better transferability [41,42].Other MD methods applicable to study ILs include ab initio MD [25,43–45],MD with a polarized force field [46–53],and MD with other CG models [54–59].Apart from MD simulations,Monte Carlo (MC) simulations [60–62],electronic correlation method [63–66],and density functional theory (DFT) calculations [67–71]are also applied to studying ILs,coming up with many favorable results.

    Figure 1.Chemical structure of 1-butyl-3-methylimidazolium chloride.

    As typical complex liquids,ILs usually have abundant and distinctive phase behaviors beyond the description of simple liquid theories.Some unique phases,including the nanoscale segregation liquid (NSL) phase [72–76],ionic liquid crystal (ILC) phases [77–82],the ‘partially arrested’glassy phase[83],and the metastable crystal phase[84],were revealed by MD simulation and verified by experiment.Much attention has been aroused since the detailed knowledge of these phases can shed light on the understanding of not only ILs but also other amphiphilic complex liquids.It also provides guidance for industrial utilizations of ILs as novel solvents.In this review,we majorly focus on summarizing the phase behaviors of aprotic ILs with alkyl cationic side chains whose number of atomic groups is even.We will show that these behaviors tend to be affected highly by temperature,length of the cationic alkyl side-chains,charge distribution,and nature of cations and anions,indicating that these unique phase behaviors of ILs result from their dual ionic and organic nature [85].By reviewing these phase behaviors and the associated affecting factors and mechanisms,we aim to inspire further investigations in this direction and invite more researchers to delve into this field to develop a thorough understanding of ILs in the future.

    2.Dual ionic and organic nature of ILs

    One of the most significant features of ILs that broadens their industrial usage is that ILs inherit both the ionic nature of inorganic salts and the organic nature of organic solvents.They have the advantages of non-flammability,non-volatility,good stability,and conductivity compared to organic solvents,and low melting temperature compared to traditional inorganic salts.By varying cations and anions,millions of available ILs can be produced with various physical and chemical properties[6],which is highly beneficial to meeting certain requirements for designated applications.However,since selecting an appropriate combination of cation and anion via experiment is usually tedious or even unfeasible,it is critical to have a good knowledge of the dual ionic and organic nature of ILs to help achieve computer-aided systematic design of ILs.It has been well acknowledged that the competition between electrostatic and van der Waals (VDW)interactions is the key factor characterizing the ionic and organic nature[86,87],so analyzing these interactions in ILs by means of MD simulations with suitable force fields should be informative.

    In a recent work,Shi and Wang[85]performed a series of all-atom MD simulations for four representative ILs(1-butyl-3-methylimidazolium nitrate ([BMIM][NO3]),1-butyl-3-methylimidazolium tetrafluoroborate ([BMIM][BF4]),1-butyl-3-methylimidazolium hexafluorophosphate ([BMIM][PF6]),and 1-butyl-3-methylimidazolium bis(trifluoromethylsulfonyl)imide([BMIM][Tf2N])),and compared them with three molecular systems with different charge distributions: a typical molten inorganic salt(molten sodium chloride,NaCl),a strongly polar liquid (dimethyl sulfoxide,DMSO),and a weakly polar liquid(toluene).The dual ionic and organic nature of ILs can be depicted from the viewpoint of the cage energy landscape(CEL) by the obtained forces,vibrational force constants,intrinsic electric fields,cohesive energies,and cage energies.

    By introducing the concept of ion cage [88–96],which indicates how each ion is surrounded by several counter-ions in the first coordination shell,along with its counterpart of molecular cage [97–101]for molecular liquids,it is possible to depict the structures and dynamics of various types of liquids since the cage volume corresponds to density and the cage stability reflects dynamics.From cage structure,it is then possible to determine CEL by calculating the ensembleaveraged local energy landscape as a function of the dislocation of the central ion from the cage center.The curvature,slope,and depth of CEL can be determined as the force constantk,forceF,and activation energyEa,which is the average energy of a particle required to climb over the energy barrier and escape the cage calculated by using the harmonic approximation (figure 2(a)).

    Shi and Wang [85]employed the first moment of the vibrational density of state (VDOS) to qualitatively describe the average characteristic frequency of intermolecular vibrational modes,defined as [102]

    whereωis the frequency andI(ω) the VDOS.The total,VDW,and electrostatic forces,respectively,are compared in all the investigated liquids,and the liquid-phase cage energyUcageis defined as the average potential energy between an ion and a counter-ion in its ion cage to characterize the local ion–ion interaction in the liquid.

    Figure 2.(a)Cage energy landscape characterized by curvature,slope,and depth,corresponding to force constant,force,and activation energy experienced by molecules,respectively.(b) Schematic illustration of cage structures and cage energy landscape in inorganic salts,ionic liquids,and organic solvents.The cage energy landscape of inorganic salts is deep and steep,where as that of ionic liquids is still deep but much more gently.Organic solvents and ionic liquids have a similar slope and curvature near the minimum of the cage energy landscape,but the depths for the organic solvents are much lower.Reprinted with permission from [85].

    From these approaches,a conclusion has been drawn that similar molecular size,geometry,and component lead to comparable VDW forces in organic ions and organic molecules.They also weaken the electrostatic interactions in ILs because of charge delocalization and charge transfer.The cage energy of ILs induced by electrostatic interactions is drastically different from organic liquids.These findings indicate that the VDW interactions,which dominate the intermolecular forces and vibrational force constants,characterize the organic nature of RTILs,resulting in a similar geometry near the minimum of the CEL to organic liquids;whilst the cage energy,or the depth of the CEL,can characterize their ionic nature (figure 2(b)).Such a mechanism explains the similar and dissimilar characteristics between ILs and organic liquids,and it clarifies the blurry dual ionic and molecular nature of ILs where as the corresponding microscopic mechanism provides new insights into their phase behaviors.

    3.Nanoscale segregation liquid (NSL) phase

    A unique phase behavior in ILs that can be explained by the dual ionic and organic nature of ILs is the nanoscale segregation liquid(NSL) phase of ILs,which is different from either the simple liquid phase or the liquid crystal (LC) phase.The structure was first discovered by MD simulation and later verified by experiment [73,75,103,104].The discovery of the NSL phase may possibly boost new industrial applications of ILs,and the microscopic mechanism of NSL shall aid the engineering of ILs[75].The significance of this discovery is enhanced by the fact that the NSL exists in most IL systems with an amphiphilic cation,independent of a specific choice of the anion.

    In a series of MS-CG MD simulations on[CnMIm][NO3]with n = 4–12,Wang et al [73–75]observed a tail-aggregation phenomenon and later referred to it as the nanoscale segregation liquid (NSL) phase,as shown in figure 3.It has been found that,in a certain temperature range,while the whole liquid is macroscopically homogeneous,the tail groups of the cationic alkyl side chains form nanoscale nonpolar tail domains microscopically segregated with the continuous polar network formed by charged anions and cationic head groups.The existence of the NSL phase was later confirmed by experiments [103,105,106]using optical heterodynedetected Raman-induced Kerr effect spectroscopy(OHDRIKES),neutron diffraction,x-ray diffraction,electrospray ionization mass spectrometry (ESI-MS),and NMR measurements,as well as simulation works [107–109].

    To quantify the degree of tail aggregation in the NSL phase,the Gaussian-like heterogeneity order parameter(HOP)is defined as [75]

    to quantify the spatial heterogeneity for identical sites,whererijis the distance between sites i and j,andσ=L/N1/3with L being the side length of the cubic simulation box and N being the total number of identical sites.

    Both the radial distribution functions (RDFs) and HOPs of the simulated systems show that,compared with ILs with short alkyl side chains,which is basically in the simple liquid phase,the tail groups distribute quite heterogeneously in ILs with an intermediate alkyl side chain length.The charged domain,constituted by anions and cationic head groups,is dominated by the electrostatic interactions,where as the neutral cationic tail groups form a nonpolar tail domain with the VDW interactions among tail groups.It has also been discovered that such aggregation behavior is temperature-sensitive,which can be understood through analyzing the tail domain diffusion in ILs [75].Investigating the behavior of tail domains with increasing temperature has also revealed that the transition of ILs from NSL to simple liquid is characterized by the collective behavior of tail groups.Nevertheless,such a collective behavior is passively induced by the steric repulsion from the continuous polar network,as concluded by a further study on ILs with an external electric field applied that the repulsion from the polar part is the major cause of the aggregation of the nonpolar tail groups [76].

    4.ILC phases

    While ILs with an intermediate cationic alkyl side-chain length can form the unique NSL phase,it has been revealed by experiment that ILs with longer alkyl chains can easily form ILC phases,most of which are smectic [79,80,82,110–115].It is therefore critical to unveil the mechanism for ILs to transform from NSL to ILC phases when the length of alkyl side chains increases.To study this phase transition,Ji et al[81]and Li et al[116]adopted the EF-CG model for ILs[42]to perform MD simulations on[CnMIm][NO3]withn=6–22 [81]and 12–24 [116],respectively.Both simulations indicate that the phase transition from NSL to ILC exists when the side-chain length is increased.

    Figure 3.(a)Molecular structure and coarse-graining scheme for the simulated IL.(b)Snapshot illustrating the NSL phase.The white spheres represent the cationic terminal groups,the gold spheres represent the cationic head groups,and the red spheres represent the anions.The ellipses in blue indicate the approximate positions of the nonpolar tail domains.Reprinted with permission from [104].Copyright 2007 American Chemical Society.

    Ji et al [81]have simulated 512 ion pairs with an isotropic barostat and cubic simulation box.The HOP values reveal that the IL systems go through a structural transition when the side-chain length increases beyond 14,and the HOP of tail groups decreases drastically from C14to C16(figure 4(a)).It not only coincides with previous experimental studies [82]but also is verified through RDFs of the investigated systems.The orientation correlation function (OCF)for side chains has also been introduced,which is the ensemble-averaged correlation between the orientations of two side chains as a function of the distance between the center-of-masses (COMs) of cations:

    whereis the unit vector pointing from the head to the tail of the ith cation.

    As shown in figure 4(b),the maximum value of OCF increases with the side-chain length from 6 to 14,and stays close to 1 from 16 to 22,indicating that the side chains of neighboring cations change their relative spatial feature from aggregated to parallel-aligned.It can also be observed that,with increasing distance,the OCF values for C6to C14gradually decrease to zero and those for C16to C22gradually decrease to a finite value.This provides further evidence that the transformation is a phase transition from the NSL phase to the ILC-like phase since,in the NSL phase,side chains do not have long-range correlations,where as they have a long-range order in the ILC-like phase.It has been suggested that the increment of the VDW interactions might be the major cause for such a transition.

    Li et al’s work [116]further investigated this phase transition and interpreted it from the perspective of the percolation phase transition.Larger systems (4096 ion pairs)have been simulated with an anisotropic barostat allowing the simulation box size in three dimensions to change independently and an isotension-isothermal ensemble allowing the simulation box to change its shape.In the configurations equilibrated by simulating annealing procedures,two side chains are considered ‘connected’ when their COM distance is less than 0.72 nm and at the same time the twist angle between them is less than 30°.A set of ‘connected’ side chains is then defined as a cluster.In IL systems with an intermediate side-chain length,such as C12,the side chains have weak tendencies to form clusters (figure 5(a)).With increasing side-chain length,in the C16system,more and larger clusters formed locally because of the stronger tendency of parallel alignment of side chains,but the orientations of the clusters are still random (figure 5(b)).For C22,side chains are in parallel globally and the largest cluster is comparable to the size of the simulated box,and the secondlargest cluster is comparatively very small (figure 5(c)),implying the occurrence of a percolation phase transition[117,118].By analysing the size difference between the largest and second-largest clusters (figure 5(d)),it has been determined that the phase transition happens at C18where both the average cluster size and the correlation length reach their maxima of around 500 and 0.42,respectively,with large fluctuations (figures 5(e) and (f)).

    Figure 4.(a)HOP for cationic head groups,anions,and cationic tail groups.(b)OCF for side chains.(c)Molecular structure of the investigated IL system and the corresponding coarse-grained model.Reprinted with permission from [81].Copyright 2013 American Chemical Society.

    Figure 5.Percolation phase transition.(a)–(c)Several largest clusters in C12,C16 and C22 systems,respectively,with the largest cluster colored blue.(a)The largest cluster in C12 is very small and not well aligned.(b)The largest cluster in C16 is larger and better aligned in parallel,but still local with little orientation correlation between clusters.(c)The largest cluster in C22 almost fills in the whole simulation box,indicating that the majority of the side chains are globally aligned in parallel and well connected.(d) Normalized sizes of the largest and second-largest clusters for all systems.(e)Average cluster size versus side-chain length.After the percolation phase transition,the largest cluster is not counted in the calculation of the average cluster size.(f)Correlation length versus side-chain length.For a finite system,both the average cluster size and the correlation length reach their maxima at the phase transition point.The correlation length is directly related to the average cluster size.Reprinted with permission from [116].

    Figure 6.HOPs for some CG sites as a function of temperature.The error bars are almost invisible since they are smaller than the size of the markers.Cr represents the crystal phase,SmA the smectic A phase,and Iso the NSL phase.Reprinted with permission from[119].Copyright 2015 American Chemical Society.

    Figure 7.(a)HOP values for some CG sites as a function of time during the transition from the NSL phase into the smectic A phase,along with snapshots taken at T = 505 K after simulated for 0,48,and 72 ns.(b) HOP values for some CG sites as a function of time during the transition from the smectic A phase into the crystal phase,along with snapshots taken at T = 480 K after simulated for 0,5,and 20 ns.Reprinted with permission from [119].Copyright 2015 American Chemical Society.

    In the above simulations,although a CG model has been adopted,the simulated IL systems still suffer from limited equilibration time and finite-size effects.Therefore,the above results render future studies with even larger temporal and spatial simulation scales to refine the microscopic mechanism of the phase transition in ILs from the NSL state to the ILC state.

    5.Solid–solid and melting phase transitions

    Apart from varying side-chain lengths of the ILs,it has also been observed in computational and experimental studies that phase transitions can exist when varying temperatures or artificially tuning partial charges of ions [119,120].By MD simulations with the EF-CG model conducted on [C16MIm][NO3],Saielli et al [119]have discovered that,with increasing temperature,the phase of IL can transit firstly from crystalline to smectic A (SmA) ILC,and then from ILC to NSL,which is similar to experimental results conducted on[CnMIm][BF4][80,121].The HOP provides direct evidence for the existence of such transitions since the HOP of the tail groups (CG sites E) increases drastically between 470 K and 505 K,corresponding to the phase transition from crystal to ILC,and decreases after 550 K,from ILC to NSL (figure 6).The snapshots representing the development from ILC to crystal at 480 K and from NSL to ILC at 550 K,respectively,are illustrated in figure 7.

    Figure 8.Phase diagram of the IL system simulated with the EF-CG model.Cr represents the crystal phase,SmA the smectic A phase,and Iso the NSL phase.Reprinted with permission from [120].Copyright 2016 American Chemical Society.

    Similar phase behavior can also be observed while tuning partial charges [120].The phase transitions are observed and determined by calculating the orientational order parameter(OOP)and the translational order parameter(TOP),which are commonly used in describing the LC phase transitions.It is thus possible for the phase diagram to be determined as a function of temperature and partial charges (figure 8).This study suggests that the existence of the ILC phases is strongly increased by the total charge of the ions.When the partial charges on the CG sites are rescaled by a factor lower than 0.9,the ILC phase is absent and the crystal directly melts into the NSL phase.For the systems with larger partial charges,the thermal range for a stable ILC phase is significantly increased.The HOP of the NSL phase also suggests that the nanoscale segregation is largely affected by the partial charges: the IL systems tend to be homogenous with low partial charges while nanosegregate with larger partial charges.The IL systems melted from the ILC phase have a higher degree of nanosegregation,as measured by the HOP,than that melted directly from the crystal at the same temperature.These studies indicate that the increase of either the alkyl chain length or the total charge of the cation head group and anion intensifies the competition between hydrophobic and electrostatic interactions,which enlarges the existence ranges on the phase diagram for both the NSL phase and the ILC phase.

    Cao et al [84]have investigated the phase behavior of[CnMIm][NO3](n=4-12) during heating by manually constructing the initial crystal structures that constitute parallel polar layers composed of cationic head groups and anions as well as nonpolar regions composed of cationic alkyl side chains in between.The initial configuration is then heated by all-atom MD simulation with the AMBER force field.

    During the heating process,a solid–solid phase transition has been found below the melting point,manifested by the kinks on the caloric curves shown in figure 9.The jumps of the potential energies at the solid–solid phase transition points are much smaller than those at melting transition points.The C4system exhibits distinctiveness compared with others due to its weaker VDW interaction coming from the shorter side chains.The difference between such phase transitions of C4and C8is shown in figure 10.The TOPs and OOPs indicate that the conformations of alkyl chains of the C4system are not as ordered as those with longer side chains.

    Further increasing the simulation temperature indicates that the melting phase transition of the investigated systems consists of two steps except for C6,which has only one step.A metastable state exhibits during the melting transition when the crystalline solid phase is transformed into the NSL phase(for C4,C6,and C8)or SmA ILC phase(for C10and C12).The snapshots of the metastable states of C4and C8systems are shown in figure 9(b).The metastable state tends to be more stable along with a higher melting temperature when the alkyl chain becomes longer.These features of melting transitions can be comprehended by the competition between the free energy contributions of the effective VDW interaction in the nonpolar regions(or EF1 for abbreviation)and that of the effective electrostatic interaction in the polar layers(or EF2 for abbreviation).The existence of the unique solid–solid phase behavior of ILs can thus be attributed to the dual ionic and organic nature of ILs.

    As the side-chain length increases,EF1 becomes stronger and thus the system possesses ordered conformations during solid–solid phase transitions.During the melting process,the imbalance between EF1 and EF2 at the melting point leads to an uneven melting process of polar layers and nonpolar regions,which causes the metastable phase to appear.For the C4system,EF1 is weaker than EF2 at the melting point of 315 K,so the alkyl side chains lose the crystalline order before the polar layers do.For systems containing C8to C12,longer alkyl side chains lead to stronger EF1,and thus the increase of the melting temperatures,at which the alkyl side chains keep their orientations uniformly for a certain time whilst polar layers drift from their original lattice positions.For the C6system,the free energy contributions from the two interactions are well balanced,resulting in the absence of the metastable phase.

    The above phase transition in [CnMIm][NO3](n = 2,4,6,8,and 10) systems was later experimentally examined by Abe et al [122].In particular,they examined the phase transition of [C10mim][NO3]via simultaneous SWAXS and DSC measurements,which indicate the existence of a multistep phase transition process.

    6.Partially arrested glassy state

    γα≡ limt→∞〈 (ΔR(α))2〉,and the equation reads as [124]:

    The self-consistent generalized Langevin equation (SCGLE)theory has been used to investigate the dynamics of molten salts,and such implementation has predicted the existence of the so-called partially arrested state when the larger ions are still in the fluid state while the smaller counter-ions are arrested in the glassy state [123].The SCGLE theory puts forward a simple equation for the asymptotic value of the mean-squared displacement of speciesα,which is defined aswhereδαβ=0,α≠β;δαβ=1,α=β

    Figure 9.Caloric curves during heating in[CnMIm][NO3]systems.(a)From T = 240 K to 325 K.(b)From T = 325 K to 450 K.The solid–solid phase transition points are marked by black arrows,and the melting transition points are marked by red arrows.Reprinted with permission from [84].Copyright 2018 American Chemical Society.

    whereSis the matrix of partial structure factors,handcare the Ornstein-Zernike matrices of total and direct correlation functions,respectively,and {...}ααrepresents the diagonal elements of the matrices.The matrixis defined as andλ(k) aswherekc(α)=8.17/σαwithσαbeing the diameter of the investigated particle andkbeing the wave factor.The theoretical prediction for molten salts indicates that there are partially arrested states in the dynamic arrest line (F–G region in figure 11(a))where only one species of ions are arrested in the glassy state while the counter-ions are still in the liquid state.

    Ramírez-Gonzalez et al [83]have applied the SCGLE theory to investigating the [C2MIm][BF4]IL and found that the IL system is very likely to be in the partially arrested glassy state.Anions’ valuevaries from zero at higher temperatures than cations’ value,which suggests that the IL is in the F–G region with the anions arrested while the cations being liquid.As a consequence,although at short times,the mobility of BF4-is higher than that of C2MIm+due to the free flight in the ballistic regime,it is observed that the diffusion of the smaller BF4-becomes slower than the larger C2MIm+in the diffusive regime,as observed in the all-atom MD simulations with the Amber force field [21](figure 11(b)).By analyzing the RDFs between cations and anions,Ramírez-Gonzalez and co-workers have discovered that the anions behave as a typical Wigner glass whose mean distance of the constituents is very large.

    7.Liquid–liquid phase separation of IL/benzene mixture

    Aside from pure IL systems,the phase behaviors of IL mixtures have also attracted much attention.Specifically,since ILs are found to exhibit dual ionic and organic nature and amphiphilic features,the mixtures of ILs and organic molecules have been widely investigated [125–133].Unlike ILs directly solved in many small organic solvents,such as acetonitrile,dichloromethane,and chloroform [134],ILs mixed with benzene exhibit liquid–liquid phase separation.Holbrey et al [128]observed experimentally that aromatic molecules,such as benzene and toluene,are soluble in imidazolium-based ILs and that liquid–liquid phase separations occur in the mixtures of ILs and aromatic liquids.Later studies on the solubility of aromatic molecules in different ILs also indicate that electrostatic interactions between cations and anions can be of much decisiveness in understanding the phase separation [125,126,130–133].

    To understand the microscopic mechanisms of the above liquid–liquid phase separation,Li et al [135]performed both NMR experiments and all-atom MD simulations on a set of mixtures of benzene and viologen bistriflimide salts[Cm(bpCn)2][CF3(SO2)2N]with m and n being the numbers of carbon atoms on the two alkyl side chains.The viologen salts are crystal at room temperature but also exhibit ILC and NSL phases at higher temperatures,as common ILs do [135].

    Both NMR experiments and MD simulations confirm that,when mixed with benzene,the salt absorbs benzene molecules to form a sponge-like phase with benzene inserted in the nonpolar domains of the NSL-like structure compose of nonpolar alkyl chains surrounded by the continuous polar network formed by anions and charged cationic head groups(figure 12(b)).This mixture phase coexists with liquid benzene or crystalline viologen-based IL whichever is excessive.The upper boundaries separating the coexisting phases depend linearly on the cationic alkyl chain length of the IL(figure 12(a)) because a larger volume of nonpolar domains can absorb more benzene molecules.The lower boundaries also exhibit linear dependency on the side-chain length corresponding to the minimum amount of benzene required to liquidize the salt,which increases proportionally with the cationic side-chain length.

    Although many other investigations[125,126,128,131,133]have suggested that π–π or ion-π interactions between cations and benzene molecules may contribute to such kind of liquid–liquid phase separations,the MD simulations have indicated that they do not exist in the above benzene/viologen-based IL mixture systems,and the benzene molecules reside inside the nonpolar domains only due to their planar and nonpolar molecular feature.Further studies in this direction may investigate the influence of the alkyl chain length,the cation and anion types,and the temperature on the liquid–liquid phase separation behaviors of the mixtures composed of various ILs and small organic molecules and develop a general theoretical framework for this kind of phase behaviors.

    Figure 10.Snapshots of the structures before and after the solid–solid phase transitions.(a)C4 at 295 K before the phase transition.(b)C4 at 300 K after the first but before the second phase transition.(c) C4 at 305 K after the phase transition.(d) C8 at 250 K before the phase transition.(e)C8 at 250 K after the phase transition.All these snapshots are taken from the[100]direction.Reprinted with permission from[84].Copyright 2018 American Chemical Society.

    8.Other phase behaviors of ILs

    ILs mixed with water constitute an important topic of investigation.Many MD simulations on the mixtures of water with imidazolium-based ILs have focused on studying their structural characteristics,analyzing the RDFs and spatial distribution functions (SDFs),where as many others have focused majorly on depicting the mixture by analyzing hydrogen bonds or investigating the effect of alkyl chain length on their dynamics[136–146].Some investigators have also addressed the significance of unique structural features.Jiang et al [147]studied the influence of water molecules on the NSL structure of [C8MIm][NO3],and water molecules were found to be inserted inside the polar network of the IL.Abe and coworkers [148–151]further determined the existence of aggregations of water molecules in IL systems by discovering the percolation limit that Bernardes et al [152]predicted through MD simulations.They also brought up a hypothesis that the so-called ‘water pockets’ exist in the mixture.Bystrov et al [136]have confirmed the existence of nanoscale structures in IL/water mixtures.By using MD simulation and the experimental Pulsed Field Gradient-Stimulated Echo method,they have discovered that other than forming ‘water pockets’ or ‘nanodroplets’,water molecules tend to distribute within the hydrophilic region in the IL system.All of these studies indicate that the NSL structure of ILs plays an important role in the mixture of IL and water.

    Figure 11.(a) Arrest lines predicted by the SCGLE theory for the primitive model with a size asymmetry adequate for the simulated IL(1:3.5).The fluid region is labelled as F–F.The G–G region corresponds to fully arrested states.The partially arrested phenomenon occurs in the F–G region.The inset shows the values of the parameter 1/γα following the isobaric trajectory with a pressure of 1 atm indicated by the dotted line.Circles correspond to anions (smaller particles) and squares to cations (larger particles).(b) MSDs for [EMIM][BF4].Dashed lines represent cations and solid lines represent anions.The inset shows the theoretical calculation results of the MSDs over the isobaric trajectory(big circles in the left panel).The first two temperatures are located in the F–F region and the last one inside the F–G region.Solid and dashed lines correspond to anions and cations,respectively.Reprinted from [83]with permission from AIP Publishing.

    Figure 12.(a)Phase diagram of the mole ratio R = Nbenzene/Nviologen versus m + n,the total number of carbon groups in the cationic alkyl chains,measured by experiment.The blue symbols represent the coexistence line between liquid benzene and sponge-like liquid phase,the red symbols represent the coexistence line between solid viologen salt and sponge-like liquid phase,and the black symbols represent intermediate state points with only the sponge-like liquid phase present.(b) Schematic representation of the sponge-like phase.The typical nanoscale segregation normally observed in ionic liquids is reproduced in the sponge-like phase after benzene molecules are absorbed in the hydrophobic regions.For the sake of clarity,only a few alkyl chains are schematically drawn.Reprinted with permission from [135].Copyright 2020 American Chemical Society.

    Another set of mixing systems undergoing intense investigations is the binary IL mixtures.It has been widely accepted that,by varying cationic head groups,their sidechain length,and the nature of anions,the NSL structures of ILs can be fine-tuned[103,107,153–155].There is,however,an alternative option of mixing two or more types of ILs with different chemical structures [156].Cosby et al [157]have thoroughly studied how tuning the composition of binary IL mixtures can affect the mesoscale organization and dynamics of the mixing system with the example of[C8MIm][BF4]and[C2MIm][BF4]mixture.Through detailed x-ray scattering,neutron scattering,and MD simulation studies,Cosby has discovered that increasing the mole ratio of [C2MIm][BF4]can lead to the increase of static dielectric permittivity,εs,which is the consequence of a transition in the mesoscale morphologies of the binary mixing system.

    Unique phase separation behaviors are observed in many other mixing systems as well.In a computational study of a ternary electrolyte mixture of tetramethylene glycol dimethyl ether (tetraglymeor G4),[C2MIm][BF4],and Li salts (LiNO3and LiI),Fuladi et al [158]have discovered that resembling the phase separation of IL mixtures,the ternary mixture tends to separate into two domains with one majorly consisted of ionic species and the other consisted of only G4molecules.This phase separation has been analyzed by changing the volume fraction of G4and ILs,salt concentration,and temperature.It is believed that the phase separation process is driven primarily by entropy and is thus temperature sensitive,which slows the diffusive dynamics of Li+ions.This study will inspire future investigation and engineering of IL electrolytes that may improve the performance of Li batteries.

    9.Conclusions

    ILs are room-temperature organic salts that possess both ionic and organic characteristics,which leads to some distinctive phase behaviors both in the context of pure ILs and their mixture with other liquids.Previous theoretical,computational,and experimental studies on phase behaviors of ILs have suggested the existence of many phases in this amphiphilic liquid,whose phase behaviors can be fine-tuned by manipulating thermodynamic conditions and molecular structures,such as charge distribution,length of cationic alkyl side chains,and temperature.Mixing with other molecules can vary their phase behaviors as well.

    With an intermediate alkyl cationic side-chain length,ILs exhibit a unique NSL state which is macroscopically isotropic but microscopically phase-separated into a polar network and nonpolar domains.The ILC phase appears when the alkyl side chains are sufficiently long.For the crystalline solid of ILs,when increasing temperature,the solid–solid phase transition may occur and the melting transition from crystal to IL or ILC sometimes encounters a metastable state.At a certain thermal condition,an IL system may stay in the partially arrested glassy state.Liquid–liquid phase separation can occur when ILs mixed with some sorts of small organic molecules.

    By increasing alkyl side-chain length,ILs go over a percolation phase transition from the NSL phase to the ILC phase due to the increasing VDW interactions among alkyl side chains to transform the conformation of side-chain bundles from nanoscale aggregation to globally parallel alignment.Most of the phase behaviors of aprotic ILs can be qualitatively understood by the competition between the effective free energy contributions from the collective VDW interactions among alkyl side chains and the electrostatic interactions among charged anions and cationic head groups.

    It can be noted that the dual ionic and organic nature of ILs plays an important role in determining the unique phase behaviors of ILs.The studies mentioned in this review thus pave the way for future theoretical works investigating phase behaviors of ILs and their mixtures from this perspective.Through analyzing the free energy contributions from electrostatic and VDW interactions,analytical models may be viable to provide a unified mechanism to explain the existing computational and experimental results and provide researchers and engineers with proper guidance to finely tune the physical and chemical properties of ILs in favor of certain industrial applications.

    Acknowledgments

    This work was supported by the National Natural Science Foundation of China(Nos.11774357,22011530390,12047503)and the Chinese Academy of Sciences (No.QYZDJ-SSWSYS01).

    ORCID iDs

    免费看光身美女| 亚洲欧洲精品一区二区精品久久久| 欧美日韩福利视频一区二区| 91久久精品国产一区二区成人 | 国产一区二区三区视频了| 好看av亚洲va欧美ⅴa在| 一区二区三区国产精品乱码| 夜夜躁狠狠躁天天躁| 国产综合懂色| 久久久久国产一级毛片高清牌| 日韩有码中文字幕| 久久久久国内视频| 久久中文字幕一级| 十八禁网站免费在线| 巨乳人妻的诱惑在线观看| 亚洲 国产 在线| 欧美3d第一页| 长腿黑丝高跟| 国产极品精品免费视频能看的| 久9热在线精品视频| 久久99热这里只有精品18| 国产成+人综合+亚洲专区| 欧美成人性av电影在线观看| 又大又爽又粗| 国产激情久久老熟女| 香蕉久久夜色| 99久久久亚洲精品蜜臀av| 国产一区二区三区视频了| 久久精品亚洲精品国产色婷小说| 亚洲精品一卡2卡三卡4卡5卡| 91在线观看av| 免费在线观看视频国产中文字幕亚洲| 亚洲最大成人中文| 国产激情偷乱视频一区二区| 日韩欧美在线乱码| 色老头精品视频在线观看| 亚洲中文字幕日韩| 国产精品久久久久久亚洲av鲁大| 少妇丰满av| 亚洲av成人av| 欧美性猛交╳xxx乱大交人| 欧美黑人巨大hd| 99精品久久久久人妻精品| 免费看日本二区| 久久国产精品人妻蜜桃| 国产99白浆流出| 亚洲,欧美精品.| 性色avwww在线观看| 日韩欧美国产在线观看| 女同久久另类99精品国产91| 色综合亚洲欧美另类图片| 日本 欧美在线| 精品久久久久久成人av| 国产视频一区二区在线看| 色综合亚洲欧美另类图片| 午夜福利在线观看吧| 日韩欧美在线二视频| 国语自产精品视频在线第100页| 午夜福利欧美成人| 午夜福利在线观看吧| 亚洲第一电影网av| 午夜福利在线观看免费完整高清在 | 国产精品亚洲av一区麻豆| 国产伦一二天堂av在线观看| 巨乳人妻的诱惑在线观看| 亚洲av熟女| 久久久久国产一级毛片高清牌| 日本熟妇午夜| 欧美午夜高清在线| 亚洲国产色片| 国产午夜精品论理片| 亚洲美女黄片视频| 午夜亚洲福利在线播放| 日韩欧美国产在线观看| 亚洲中文字幕日韩| 国产一区二区在线av高清观看| 成人av在线播放网站| 国产激情欧美一区二区| 18禁观看日本| 国产精品女同一区二区软件 | 午夜a级毛片| 国产精品 国内视频| 99久久成人亚洲精品观看| 99视频精品全部免费 在线 | 国产午夜精品论理片| 成人三级黄色视频| 欧美xxxx黑人xx丫x性爽| 中文在线观看免费www的网站| 午夜两性在线视频| 琪琪午夜伦伦电影理论片6080| 一进一出抽搐动态| 岛国在线观看网站| 美女高潮喷水抽搐中文字幕| 毛片女人毛片| 免费看日本二区| 午夜福利视频1000在线观看| 国产成人福利小说| 少妇丰满av| 精品国产三级普通话版| av国产免费在线观看| 19禁男女啪啪无遮挡网站| 夜夜爽天天搞| 久久久久性生活片| 丰满人妻熟妇乱又伦精品不卡| 亚洲午夜理论影院| 亚洲电影在线观看av| 巨乳人妻的诱惑在线观看| 亚洲av日韩精品久久久久久密| 国产高清有码在线观看视频| 女警被强在线播放| 最近最新中文字幕大全电影3| 国产精品乱码一区二三区的特点| ponron亚洲| 精品国产美女av久久久久小说| 国产精品av久久久久免费| 亚洲国产色片| 国产 一区 欧美 日韩| 日韩国内少妇激情av| av天堂在线播放| 久久香蕉国产精品| 国模一区二区三区四区视频 | 日本黄色视频三级网站网址| 一级a爱片免费观看的视频| 19禁男女啪啪无遮挡网站| 又大又爽又粗| 午夜福利成人在线免费观看| 中文字幕高清在线视频| 日韩国内少妇激情av| 激情在线观看视频在线高清| 中文亚洲av片在线观看爽| 嫩草影院入口| 国产亚洲欧美98| 成人一区二区视频在线观看| 国产精品亚洲一级av第二区| 少妇人妻一区二区三区视频| 欧美日韩中文字幕国产精品一区二区三区| 中国美女看黄片| 国产v大片淫在线免费观看| 国产1区2区3区精品| 日本黄色视频三级网站网址| 国产高清激情床上av| 久久中文字幕一级| 美女 人体艺术 gogo| 欧美一级a爱片免费观看看| 亚洲欧美精品综合一区二区三区| 亚洲成人精品中文字幕电影| 久久久久久国产a免费观看| 久久久久九九精品影院| 精品久久久久久成人av| 国产免费男女视频| 嫩草影视91久久| 波多野结衣高清作品| 麻豆av在线久日| 亚洲avbb在线观看| 搞女人的毛片| 午夜福利在线在线| 亚洲熟妇中文字幕五十中出| 亚洲精品久久国产高清桃花| 日本黄色片子视频| 午夜福利成人在线免费观看| 午夜日韩欧美国产| 欧美日韩黄片免| 久久天堂一区二区三区四区| 亚洲欧美一区二区三区黑人| x7x7x7水蜜桃| 怎么达到女性高潮| 中文亚洲av片在线观看爽| 国产不卡一卡二| 精品国产乱码久久久久久男人| 亚洲人成伊人成综合网2020| 成人一区二区视频在线观看| 亚洲电影在线观看av| 国产主播在线观看一区二区| 麻豆成人午夜福利视频| 婷婷精品国产亚洲av| 岛国视频午夜一区免费看| 两人在一起打扑克的视频| 亚洲色图av天堂| www.999成人在线观看| e午夜精品久久久久久久| 亚洲av五月六月丁香网| 最近视频中文字幕2019在线8| 色综合欧美亚洲国产小说| 久久精品91蜜桃| 久久中文字幕人妻熟女| 久久久精品大字幕| 一a级毛片在线观看| 两性夫妻黄色片| 香蕉久久夜色| 变态另类丝袜制服| 午夜福利成人在线免费观看| 国产精品影院久久| 青草久久国产| 欧美不卡视频在线免费观看| 久久久久久大精品| 精品一区二区三区av网在线观看| 少妇人妻一区二区三区视频| 国产成人影院久久av| aaaaa片日本免费| www国产在线视频色| 亚洲av片天天在线观看| 国产午夜福利久久久久久| av在线蜜桃| 婷婷六月久久综合丁香| 精品无人区乱码1区二区| 色老头精品视频在线观看| 日韩三级视频一区二区三区| 成人亚洲精品av一区二区| 亚洲中文字幕一区二区三区有码在线看 | 人妻久久中文字幕网| 麻豆久久精品国产亚洲av| 极品教师在线免费播放| 国产av麻豆久久久久久久| 日韩 欧美 亚洲 中文字幕| 亚洲色图 男人天堂 中文字幕| 国产熟女xx| 久久精品国产清高在天天线| 99久国产av精品| 国产极品精品免费视频能看的| 午夜亚洲福利在线播放| 免费观看精品视频网站| 欧美日本亚洲视频在线播放| 久久香蕉国产精品| 这个男人来自地球电影免费观看| 免费高清视频大片| xxxwww97欧美| 中文字幕av在线有码专区| 亚洲国产精品999在线| 99国产精品一区二区蜜桃av| 国产成人啪精品午夜网站| 亚洲中文字幕一区二区三区有码在线看 | 日韩欧美免费精品| 桃色一区二区三区在线观看| 国产精品野战在线观看| 久久久久亚洲av毛片大全| 舔av片在线| 日韩欧美在线二视频| 国产乱人视频| 叶爱在线成人免费视频播放| 午夜成年电影在线免费观看| 亚洲国产高清在线一区二区三| 久久香蕉国产精品| 久久精品国产99精品国产亚洲性色| 白带黄色成豆腐渣| 亚洲片人在线观看| 日韩高清综合在线| 丁香六月欧美| 国产一区在线观看成人免费| 久久精品91无色码中文字幕| 夜夜躁狠狠躁天天躁| 久久人人精品亚洲av| 18禁裸乳无遮挡免费网站照片| 日韩大尺度精品在线看网址| 亚洲成人久久性| 精品久久久久久久毛片微露脸| 女人被狂操c到高潮| 黑人巨大精品欧美一区二区mp4| 亚洲av日韩精品久久久久久密| 精品午夜福利视频在线观看一区| 精品不卡国产一区二区三区| 五月玫瑰六月丁香| 日本一二三区视频观看| 亚洲 欧美一区二区三区| 不卡一级毛片| 亚洲中文av在线| www.熟女人妻精品国产| 男女午夜视频在线观看| 亚洲欧美日韩东京热| 国产三级在线视频| 亚洲色图 男人天堂 中文字幕| 白带黄色成豆腐渣| 国产精品精品国产色婷婷| 91老司机精品| 美女午夜性视频免费| 亚洲人成伊人成综合网2020| 国内精品美女久久久久久| 国产精品 国内视频| 久久精品亚洲精品国产色婷小说| 国产高清视频在线观看网站| 色av中文字幕| 此物有八面人人有两片| 午夜福利视频1000在线观看| 亚洲午夜精品一区,二区,三区| 国产精品九九99| 真实男女啪啪啪动态图| 午夜福利欧美成人| xxxwww97欧美| 午夜视频精品福利| www.www免费av| 欧美精品啪啪一区二区三区| 亚洲无线在线观看| 男女那种视频在线观看| 一区福利在线观看| xxxwww97欧美| 日本免费一区二区三区高清不卡| 十八禁人妻一区二区| www日本在线高清视频| 90打野战视频偷拍视频| 夜夜躁狠狠躁天天躁| 精品国产美女av久久久久小说| 亚洲欧美精品综合久久99| 国产三级中文精品| 国产亚洲精品久久久com| 免费av不卡在线播放| 国产伦人伦偷精品视频| 天堂av国产一区二区熟女人妻| 国产精品一及| 午夜福利成人在线免费观看| 无遮挡黄片免费观看| 在线免费观看不下载黄p国产 | 99久久国产精品久久久| 亚洲 欧美 日韩 在线 免费| 成人午夜高清在线视频| 久久精品亚洲精品国产色婷小说| 综合色av麻豆| 可以在线观看毛片的网站| 丰满人妻一区二区三区视频av | 午夜精品一区二区三区免费看| 男女做爰动态图高潮gif福利片| 久久精品国产99精品国产亚洲性色| 男人舔女人的私密视频| 麻豆国产av国片精品| 国产视频一区二区在线看| 一个人观看的视频www高清免费观看 | 久久国产乱子伦精品免费另类| 国产精品美女特级片免费视频播放器 | 亚洲精品中文字幕一二三四区| 少妇的逼水好多| 亚洲熟妇熟女久久| 久久久久久久久免费视频了| 婷婷六月久久综合丁香| 欧美另类亚洲清纯唯美| 精华霜和精华液先用哪个| 热99re8久久精品国产| 两个人的视频大全免费| 亚洲五月婷婷丁香| 亚洲精品一卡2卡三卡4卡5卡| 69av精品久久久久久| 亚洲国产日韩欧美精品在线观看 | 美女 人体艺术 gogo| 久久久久久人人人人人| 麻豆国产av国片精品| 国产99白浆流出| 美女cb高潮喷水在线观看 | 麻豆国产av国片精品| 欧美三级亚洲精品| 又大又爽又粗| 黑人巨大精品欧美一区二区mp4| 最近最新免费中文字幕在线| 免费看光身美女| 国产伦人伦偷精品视频| 成人永久免费在线观看视频| 亚洲狠狠婷婷综合久久图片| 99久久精品一区二区三区| 亚洲18禁久久av| 99热这里只有是精品50| avwww免费| 一卡2卡三卡四卡精品乱码亚洲| 国产成人系列免费观看| 精品久久蜜臀av无| 亚洲精品一区av在线观看| 一进一出抽搐动态| 成人国产一区最新在线观看| bbb黄色大片| 成人国产一区最新在线观看| 国产麻豆成人av免费视频| 色哟哟哟哟哟哟| 亚洲欧美日韩高清专用| 午夜福利在线观看免费完整高清在 | 精品乱码久久久久久99久播| 国产 一区 欧美 日韩| 18禁国产床啪视频网站| 国产aⅴ精品一区二区三区波| 国产精品久久久久久精品电影| 精品久久久久久久久久久久久| 午夜a级毛片| 欧美中文综合在线视频| 国内揄拍国产精品人妻在线| 日韩欧美在线二视频| 99久久精品热视频| 国产人伦9x9x在线观看| 热99在线观看视频| 亚洲国产中文字幕在线视频| 一区福利在线观看| 亚洲第一欧美日韩一区二区三区| h日本视频在线播放| 两个人的视频大全免费| 两性午夜刺激爽爽歪歪视频在线观看| 悠悠久久av| 久久精品综合一区二区三区| 欧美乱妇无乱码| 久久久久免费精品人妻一区二区| 久久精品国产亚洲av香蕉五月| 久久午夜亚洲精品久久| 搡老熟女国产l中国老女人| 99久久精品热视频| 午夜激情福利司机影院| 欧美不卡视频在线免费观看| 久久香蕉国产精品| 成人三级做爰电影| 久久久久久九九精品二区国产| 午夜亚洲福利在线播放| 制服人妻中文乱码| 无遮挡黄片免费观看| 美女cb高潮喷水在线观看 | 亚洲国产欧美人成| 99久久精品一区二区三区| 亚洲国产欧洲综合997久久,| 天堂网av新在线| 国产精品一区二区三区四区久久| 国产aⅴ精品一区二区三区波| 午夜免费观看网址| 亚洲av成人精品一区久久| 日韩欧美三级三区| 久久久久久久久免费视频了| 88av欧美| 国产成人av激情在线播放| 丰满人妻一区二区三区视频av | 成人三级做爰电影| 色噜噜av男人的天堂激情| 欧美日韩亚洲国产一区二区在线观看| 国产一区在线观看成人免费| 欧美中文综合在线视频| 九九久久精品国产亚洲av麻豆 | 每晚都被弄得嗷嗷叫到高潮| 一本一本综合久久| 我的老师免费观看完整版| 欧美乱妇无乱码| av在线蜜桃| av女优亚洲男人天堂 | 国产欧美日韩一区二区精品| 午夜福利在线观看免费完整高清在 | 亚洲精品乱码久久久v下载方式 | 亚洲av成人av| 校园春色视频在线观看| 午夜影院日韩av| 午夜激情福利司机影院| 日本成人三级电影网站| 在线观看美女被高潮喷水网站 | 亚洲九九香蕉| 亚洲专区字幕在线| 最近最新中文字幕大全免费视频| 色老头精品视频在线观看| 一个人观看的视频www高清免费观看 | 精华霜和精华液先用哪个| 日本三级黄在线观看| 在线免费观看的www视频| 国产伦精品一区二区三区四那| 亚洲狠狠婷婷综合久久图片| 757午夜福利合集在线观看| or卡值多少钱| 波多野结衣巨乳人妻| 俄罗斯特黄特色一大片| 国产精品98久久久久久宅男小说| 一本一本综合久久| 日本与韩国留学比较| 99热这里只有精品一区 | 人妻久久中文字幕网| 国产主播在线观看一区二区| 91老司机精品| 国产精品精品国产色婷婷| 亚洲午夜精品一区,二区,三区| 久久久久九九精品影院| 免费在线观看影片大全网站| 一进一出抽搐动态| 午夜福利欧美成人| 亚洲在线观看片| 亚洲狠狠婷婷综合久久图片| 欧美日韩精品网址| 亚洲国产欧美一区二区综合| 人妻夜夜爽99麻豆av| 亚洲人成伊人成综合网2020| 久久精品aⅴ一区二区三区四区| 国产成人av激情在线播放| or卡值多少钱| 婷婷精品国产亚洲av| 国产欧美日韩一区二区精品| 天天躁狠狠躁夜夜躁狠狠躁| 中文在线观看免费www的网站| 国产一区二区激情短视频| 免费在线观看视频国产中文字幕亚洲| 午夜免费激情av| 18禁观看日本| 校园春色视频在线观看| 国产av不卡久久| 又黄又粗又硬又大视频| 这个男人来自地球电影免费观看| 97碰自拍视频| 日本在线视频免费播放| 女生性感内裤真人,穿戴方法视频| 国产伦在线观看视频一区| 久久天躁狠狠躁夜夜2o2o| 麻豆久久精品国产亚洲av| 欧美黄色淫秽网站| 他把我摸到了高潮在线观看| 国产激情偷乱视频一区二区| 91久久精品国产一区二区成人 | 美女大奶头视频| 级片在线观看| 最近视频中文字幕2019在线8| 婷婷精品国产亚洲av| 99精品久久久久人妻精品| 91av网一区二区| 一本精品99久久精品77| 天天添夜夜摸| 国产精品亚洲一级av第二区| 亚洲欧美激情综合另类| 一进一出好大好爽视频| 9191精品国产免费久久| 久久久久国内视频| av中文乱码字幕在线| av黄色大香蕉| 久久久久久大精品| 久久国产精品影院| 成人av在线播放网站| 欧美日本亚洲视频在线播放| 级片在线观看| 女生性感内裤真人,穿戴方法视频| 欧美精品啪啪一区二区三区| 麻豆国产av国片精品| 美女免费视频网站| 免费观看精品视频网站| 成人午夜高清在线视频| 午夜视频精品福利| 国产精品久久电影中文字幕| 亚洲中文av在线| 国产乱人伦免费视频| 亚洲成人久久爱视频| 久久精品夜夜夜夜夜久久蜜豆| 中文字幕久久专区| 精品免费久久久久久久清纯| 我的老师免费观看完整版| 婷婷丁香在线五月| 国产成人福利小说| 成人国产综合亚洲| 国产成人aa在线观看| avwww免费| 国产1区2区3区精品| 久久亚洲真实| 久久欧美精品欧美久久欧美| 级片在线观看| 免费大片18禁| 999精品在线视频| 成熟少妇高潮喷水视频| 亚洲国产欧美人成| 三级男女做爰猛烈吃奶摸视频| 狂野欧美白嫩少妇大欣赏| 久久婷婷人人爽人人干人人爱| 国产成人精品久久二区二区免费| 亚洲精品久久国产高清桃花| 亚洲乱码一区二区免费版| 亚洲人成网站高清观看| 少妇裸体淫交视频免费看高清| 亚洲国产高清在线一区二区三| 欧美中文日本在线观看视频| 亚洲av成人av| 免费无遮挡裸体视频| 国产精品美女特级片免费视频播放器 | 免费在线观看成人毛片| 91麻豆精品激情在线观看国产| 国产午夜精品论理片| 欧美成狂野欧美在线观看| 啦啦啦免费观看视频1| 国产探花在线观看一区二区| 夜夜爽天天搞| 女人高潮潮喷娇喘18禁视频| 国产激情久久老熟女| 亚洲,欧美精品.| 国产精品电影一区二区三区| 听说在线观看完整版免费高清| 不卡av一区二区三区| а√天堂www在线а√下载| 国产乱人伦免费视频| 麻豆av在线久日| 成年女人看的毛片在线观看| 久久久国产成人免费| 国产av在哪里看| 人妻久久中文字幕网| 18禁国产床啪视频网站| 国产激情欧美一区二区| 俄罗斯特黄特色一大片| 国产毛片a区久久久久| 国产成人影院久久av| 老汉色∧v一级毛片| x7x7x7水蜜桃| 午夜福利欧美成人| 九色国产91popny在线| 亚洲中文av在线| 特级一级黄色大片| 国模一区二区三区四区视频 | 一级作爱视频免费观看| 精品人妻1区二区| 熟女少妇亚洲综合色aaa.| 18禁国产床啪视频网站| 成人欧美大片| 午夜亚洲福利在线播放| 久久香蕉精品热| 国产精品美女特级片免费视频播放器 | 窝窝影院91人妻| 亚洲国产欧美人成| 免费在线观看视频国产中文字幕亚洲| 久久亚洲真实| 麻豆久久精品国产亚洲av| 99久久国产精品久久久| 欧美在线黄色| 国产高清videossex| 国产乱人视频| 美女高潮的动态| 久久精品aⅴ一区二区三区四区| 亚洲av成人精品一区久久| 亚洲av电影在线进入| 欧美精品啪啪一区二区三区| 亚洲五月婷婷丁香| 国产精品久久久久久亚洲av鲁大| 欧美日韩亚洲国产一区二区在线观看|