• <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在线精品国自产拍蜜月| 国产精品精品国产色婷婷| 欧美一区二区国产精品久久精品| 神马国产精品三级电影在线观看| 91久久精品电影网| 国产成人精品久久久久久| 精品一区二区三区人妻视频| 成人亚洲欧美一区二区av| 日韩,欧美,国产一区二区三区 | 精品午夜福利在线看| 国产精品一及| 日韩制服骚丝袜av| 久久99热这里只有精品18| 亚洲美女视频黄频| 一a级毛片在线观看| 午夜福利18| 国产一区二区三区av在线 | 日韩精品中文字幕看吧| 色综合站精品国产| 深夜a级毛片| 欧美日本亚洲视频在线播放| 人人妻,人人澡人人爽秒播| 男女啪啪激烈高潮av片| 国产女主播在线喷水免费视频网站 | 国产片特级美女逼逼视频| 好男人在线观看高清免费视频| 成人无遮挡网站| 亚洲av熟女| 日日摸夜夜添夜夜爱| av视频在线观看入口| 我要搜黄色片| 少妇裸体淫交视频免费看高清| av卡一久久| 丝袜美腿在线中文| 色播亚洲综合网| 色播亚洲综合网| 六月丁香七月| 亚洲婷婷狠狠爱综合网| 欧美成人一区二区免费高清观看| 亚洲高清免费不卡视频| 亚洲婷婷狠狠爱综合网| 18禁在线无遮挡免费观看视频 | 亚洲在线观看片| 午夜爱爱视频在线播放| 观看美女的网站| 日本精品一区二区三区蜜桃| 精品午夜福利在线看| 男女边吃奶边做爰视频| 久久久久久国产a免费观看| 舔av片在线| 搡老妇女老女人老熟妇| 夜夜爽天天搞| 一本久久中文字幕| 可以在线观看的亚洲视频| 亚洲成人久久爱视频| 18禁裸乳无遮挡免费网站照片| 久久韩国三级中文字幕| 最后的刺客免费高清国语| 99久久精品国产国产毛片| 在现免费观看毛片| 99热全是精品| 欧美成人精品欧美一级黄| 亚洲国产精品国产精品| 国产乱人偷精品视频| 成人美女网站在线观看视频| 麻豆一二三区av精品| 日日摸夜夜添夜夜添av毛片| 精品少妇黑人巨大在线播放 | 黄色日韩在线| 亚洲电影在线观看av| 小说图片视频综合网站| 中文亚洲av片在线观看爽| 欧美区成人在线视频| 色在线成人网| 综合色丁香网| 国产v大片淫在线免费观看| 淫妇啪啪啪对白视频| 尾随美女入室| 看片在线看免费视频| 在线观看66精品国产| 97人妻精品一区二区三区麻豆| 亚洲无线在线观看| 少妇人妻精品综合一区二区 | 欧美色视频一区免费| 日本三级黄在线观看| 久久精品国产亚洲av涩爱 | 国产老妇女一区| 亚洲欧美日韩无卡精品| 久久综合国产亚洲精品| 日本黄色片子视频| 91在线观看av| 亚洲美女黄片视频| 日韩欧美一区二区三区在线观看| 精品人妻视频免费看| 1000部很黄的大片| 午夜福利18| 91麻豆精品激情在线观看国产| 波多野结衣高清作品| 日韩强制内射视频| 综合色丁香网| 男人的好看免费观看在线视频| 久久精品夜色国产| 国产亚洲欧美98| 亚洲图色成人| 午夜视频国产福利| 日韩欧美在线乱码| 欧美另类亚洲清纯唯美| 最新在线观看一区二区三区| 一级黄片播放器| 极品教师在线视频| 亚洲无线观看免费| 一级a爱片免费观看的视频| 国产精品,欧美在线| 日本免费a在线| 丰满乱子伦码专区| 久久久久久大精品| 97超视频在线观看视频| 日韩强制内射视频| 精品久久久久久久末码| 国产激情偷乱视频一区二区| 精品人妻视频免费看| 男女视频在线观看网站免费| 欧美国产日韩亚洲一区| 日本免费一区二区三区高清不卡| 欧美性猛交╳xxx乱大交人| 麻豆精品久久久久久蜜桃| 免费人成在线观看视频色| 搞女人的毛片| 国产 一区精品| 久久综合国产亚洲精品| 亚洲av不卡在线观看| 国产午夜精品论理片| 男人舔奶头视频| 波多野结衣高清作品| 毛片女人毛片| 免费观看精品视频网站| 国产午夜福利久久久久久| 国产真实伦视频高清在线观看| 18禁在线播放成人免费| 美女cb高潮喷水在线观看| 日本撒尿小便嘘嘘汇集6| 国产精品美女特级片免费视频播放器| 欧美丝袜亚洲另类| 午夜精品国产一区二区电影 | 久久人人爽人人片av| 一卡2卡三卡四卡精品乱码亚洲| 午夜福利视频1000在线观看| 两个人视频免费观看高清| 一级毛片电影观看 | 干丝袜人妻中文字幕| 色尼玛亚洲综合影院| 青春草视频在线免费观看| 久久久久国产网址| 日本黄色视频三级网站网址| 男女啪啪激烈高潮av片| 婷婷色综合大香蕉| 五月玫瑰六月丁香| 91久久精品国产一区二区三区| 一级毛片久久久久久久久女| 国产精品女同一区二区软件| 欧美激情在线99| 给我免费播放毛片高清在线观看| 国产麻豆成人av免费视频| 女的被弄到高潮叫床怎么办| 精品人妻一区二区三区麻豆 | 美女被艹到高潮喷水动态| 亚洲国产精品合色在线| 插逼视频在线观看| 亚洲国产精品sss在线观看| 日韩制服骚丝袜av| 最新中文字幕久久久久| 91麻豆精品激情在线观看国产| 舔av片在线| 亚洲欧美清纯卡通| 成年免费大片在线观看| 国产大屁股一区二区在线视频| 中文字幕av在线有码专区| 两性午夜刺激爽爽歪歪视频在线观看| 日本成人三级电影网站| 亚洲av中文av极速乱| 久久婷婷人人爽人人干人人爱| 丰满乱子伦码专区| 国产成人影院久久av| 免费在线观看成人毛片| 午夜久久久久精精品| 日韩成人伦理影院| 免费av不卡在线播放| 色综合色国产| 永久网站在线| 18禁黄网站禁片免费观看直播| 深夜精品福利| 免费看美女性在线毛片视频| 99riav亚洲国产免费| 色视频www国产| 国产高清三级在线| 又黄又爽又免费观看的视频| 免费在线观看影片大全网站| 午夜a级毛片| 亚洲内射少妇av| 国产欧美日韩精品亚洲av| 十八禁网站免费在线| 男女下面进入的视频免费午夜| 免费高清视频大片| 我要搜黄色片| 国内久久婷婷六月综合欲色啪| 亚洲成人av在线免费| 久久中文看片网| 在现免费观看毛片| 亚洲性久久影院| 黄片wwwwww| 欧美成人a在线观看| 十八禁国产超污无遮挡网站| 大型黄色视频在线免费观看| 欧美激情在线99| 欧美最黄视频在线播放免费| 亚洲欧美日韩东京热| 99久久九九国产精品国产免费| 午夜影院日韩av| 欧美最黄视频在线播放免费| 亚洲色图av天堂| 97人妻精品一区二区三区麻豆| 午夜福利视频1000在线观看| 国产黄片美女视频| 免费大片18禁| 免费人成在线观看视频色| 午夜a级毛片| 最近的中文字幕免费完整| 国产高清视频在线播放一区| 亚洲精品色激情综合| 久久精品国产99精品国产亚洲性色| 亚洲欧美精品自产自拍| 亚洲婷婷狠狠爱综合网| 中文字幕久久专区| 白带黄色成豆腐渣| 久久精品国产亚洲av天美| 18+在线观看网站| av天堂在线播放| 我要搜黄色片| 99热这里只有精品一区| 精品不卡国产一区二区三区| 99热只有精品国产| 国产淫片久久久久久久久| 欧美最黄视频在线播放免费| 日韩av不卡免费在线播放| av在线天堂中文字幕| 性欧美人与动物交配| 黄色配什么色好看| 国产精品无大码| 床上黄色一级片| 一本一本综合久久| 午夜福利成人在线免费观看| 国产精品乱码一区二三区的特点| 五月伊人婷婷丁香| 十八禁网站免费在线| 成人av在线播放网站| 亚洲人成网站高清观看| 一进一出好大好爽视频| 99热这里只有精品一区| 日本色播在线视频| 久久久久性生活片| 男女啪啪激烈高潮av片| 简卡轻食公司| 免费在线观看影片大全网站| 成人无遮挡网站| 免费看光身美女| 97在线视频观看| 国产aⅴ精品一区二区三区波| 无遮挡黄片免费观看| 亚洲欧美日韩高清专用| 欧美日韩一区二区视频在线观看视频在线 | 国产精品综合久久久久久久免费| 99久国产av精品国产电影| 看十八女毛片水多多多| 欧美zozozo另类| 亚洲欧美日韩高清在线视频| 亚洲精品一区av在线观看| av在线天堂中文字幕| 一进一出抽搐动态| 色5月婷婷丁香| 91av网一区二区| 亚洲五月天丁香| 精品一区二区三区av网在线观看| 日韩欧美 国产精品| 丰满乱子伦码专区| 国产精品久久久久久久久免| h日本视频在线播放| 黑人高潮一二区| 一进一出抽搐动态| 熟妇人妻久久中文字幕3abv| 亚洲美女视频黄频| 老司机影院成人| a级一级毛片免费在线观看| 国产午夜精品久久久久久一区二区三区 | 免费大片18禁| 欧美高清性xxxxhd video| 欧美日韩综合久久久久久| 久久精品影院6| 少妇熟女欧美另类| 亚洲国产精品合色在线| 婷婷六月久久综合丁香| 国产成人影院久久av| 亚洲在线自拍视频| 国产三级中文精品| 日日摸夜夜添夜夜爱| 国产在线精品亚洲第一网站| 亚洲av中文av极速乱| 亚洲国产欧洲综合997久久,| 在线看三级毛片| 在线观看免费视频日本深夜| 国产精品三级大全| 国产精品美女特级片免费视频播放器| 日韩成人伦理影院| 国产伦精品一区二区三区视频9| 国内精品宾馆在线| 日韩在线高清观看一区二区三区| 少妇丰满av| 给我免费播放毛片高清在线观看| 国产亚洲精品综合一区在线观看| 五月伊人婷婷丁香| 国内精品宾馆在线| 国产亚洲精品av在线| 一a级毛片在线观看| 亚洲最大成人手机在线| 国产精品不卡视频一区二区| eeuss影院久久| 久久久国产成人免费| 一个人看视频在线观看www免费| 日韩制服骚丝袜av| 亚洲国产色片| 国产伦在线观看视频一区| 久久精品国产清高在天天线| 国产精品美女特级片免费视频播放器| 内射极品少妇av片p| 国产精品,欧美在线| 午夜爱爱视频在线播放| 精品少妇黑人巨大在线播放 | 免费观看精品视频网站| 亚洲成人精品中文字幕电影| 亚洲av免费高清在线观看| 国产精品不卡视频一区二区| 日韩制服骚丝袜av| 狂野欧美白嫩少妇大欣赏| 久久久久国产网址| 哪里可以看免费的av片| 国产单亲对白刺激| 中文亚洲av片在线观看爽| 国产一区二区三区在线臀色熟女| 日本三级黄在线观看| 神马国产精品三级电影在线观看| 国产成人a∨麻豆精品| 亚洲丝袜综合中文字幕| eeuss影院久久| 国产片特级美女逼逼视频| 日本熟妇午夜| 欧美精品国产亚洲| 晚上一个人看的免费电影| 免费观看在线日韩| 男女之事视频高清在线观看| 国产成人a区在线观看| 亚洲精品在线观看二区| 岛国在线免费视频观看| 免费观看精品视频网站| 男女边吃奶边做爰视频| 久久久欧美国产精品| 久久国产乱子免费精品| 校园春色视频在线观看| 亚洲精品成人久久久久久| 成人av在线播放网站| 成人高潮视频无遮挡免费网站| 成年女人永久免费观看视频| 亚洲在线观看片| 波多野结衣高清作品| 人妻丰满熟妇av一区二区三区| 深夜精品福利| 美女被艹到高潮喷水动态| 床上黄色一级片| 一级毛片aaaaaa免费看小| 天堂√8在线中文| 中文字幕久久专区| 欧美一区二区国产精品久久精品| 啦啦啦观看免费观看视频高清| 久久久久久久久大av| 国产精品永久免费网站| 99热网站在线观看| 亚洲在线观看片| 床上黄色一级片| 午夜视频国产福利| 免费观看人在逋| 日本熟妇午夜| 欧美成人一区二区免费高清观看| 成年女人看的毛片在线观看| 欧美色视频一区免费| 欧美一区二区精品小视频在线| 国产精品美女特级片免费视频播放器| 99国产极品粉嫩在线观看| 日本黄色片子视频| 男女下面进入的视频免费午夜| 亚洲精品乱码久久久v下载方式| 国产蜜桃级精品一区二区三区| 国产大屁股一区二区在线视频| 亚洲av免费在线观看| 中文字幕久久专区| 国产精品日韩av在线免费观看| 日韩欧美在线乱码| а√天堂www在线а√下载| 日韩成人伦理影院| 亚洲国产欧美人成| 老师上课跳d突然被开到最大视频| 蜜臀久久99精品久久宅男| 亚洲高清免费不卡视频| 国产午夜精品论理片| 国产成人福利小说| 国产在视频线在精品| 搡老岳熟女国产| 三级毛片av免费| 中文亚洲av片在线观看爽| 亚洲欧美日韩卡通动漫| 久久久久久久久中文| 美女黄网站色视频| 99热这里只有精品一区| 亚洲真实伦在线观看| 久久久色成人| 免费av不卡在线播放| 午夜亚洲福利在线播放| 精品久久久久久久人妻蜜臀av| 免费不卡的大黄色大毛片视频在线观看 | 秋霞在线观看毛片| 在线观看免费视频日本深夜| 日日干狠狠操夜夜爽| 国产精品不卡视频一区二区| 长腿黑丝高跟| 午夜久久久久精精品| 一a级毛片在线观看| 夜夜夜夜夜久久久久| 久久久久久九九精品二区国产| 国模一区二区三区四区视频| 久久精品国产亚洲av涩爱 | 中国国产av一级| 淫妇啪啪啪对白视频| 一a级毛片在线观看| 国产麻豆成人av免费视频| 久久精品91蜜桃| 美女内射精品一级片tv| 婷婷亚洲欧美| 免费搜索国产男女视频| 在线天堂最新版资源| 能在线免费观看的黄片| 成人三级黄色视频| 国产综合懂色| 亚洲国产欧美人成| 久久久久国产网址| 给我免费播放毛片高清在线观看| 神马国产精品三级电影在线观看| 成熟少妇高潮喷水视频| 亚洲av美国av| 欧美日本视频| 18禁在线无遮挡免费观看视频 | 女人被狂操c到高潮| 国内揄拍国产精品人妻在线| 国内精品久久久久精免费| АⅤ资源中文在线天堂| 天堂网av新在线| 亚洲美女视频黄频| 成年av动漫网址| 少妇丰满av| av在线老鸭窝| 最近视频中文字幕2019在线8| 久久精品国产亚洲av涩爱 | 男女那种视频在线观看| 最近视频中文字幕2019在线8| 婷婷亚洲欧美| 99在线视频只有这里精品首页| 国产男靠女视频免费网站| 内地一区二区视频在线| 麻豆国产97在线/欧美| av在线老鸭窝| a级毛片免费高清观看在线播放| 内射极品少妇av片p| 在线观看午夜福利视频| 日本免费一区二区三区高清不卡| 晚上一个人看的免费电影| 91麻豆精品激情在线观看国产| 国产精品野战在线观看| av在线播放精品| 欧美三级亚洲精品| 美女 人体艺术 gogo| 91久久精品电影网| 亚洲欧美中文字幕日韩二区| 夜夜看夜夜爽夜夜摸| 18禁裸乳无遮挡免费网站照片| 日本成人三级电影网站| 欧美成人一区二区免费高清观看| 白带黄色成豆腐渣| 免费搜索国产男女视频| 免费无遮挡裸体视频| 免费搜索国产男女视频| 人妻制服诱惑在线中文字幕| 一区二区三区四区激情视频 | h日本视频在线播放| 欧美激情在线99| 婷婷亚洲欧美| 69av精品久久久久久| 久久久久国产精品人妻aⅴ院| 综合色丁香网| 亚洲,欧美,日韩| 欧美性猛交╳xxx乱大交人| 久久久久久久久中文| 噜噜噜噜噜久久久久久91| 精品久久久久久久人妻蜜臀av| 国产精品日韩av在线免费观看| 亚洲丝袜综合中文字幕| 男女边吃奶边做爰视频| 国产伦精品一区二区三区四那| 国产成人福利小说| 国产精品野战在线观看| 久久久久久国产a免费观看| 久久久久久久久中文| 久久久精品94久久精品| 国产激情偷乱视频一区二区| АⅤ资源中文在线天堂| 久久久精品大字幕| 长腿黑丝高跟| 亚洲国产精品合色在线| 在线免费十八禁| 日韩中字成人| 蜜桃久久精品国产亚洲av| 少妇的逼好多水| 国产毛片a区久久久久| 在现免费观看毛片| av国产免费在线观看| 十八禁国产超污无遮挡网站| 国产精品日韩av在线免费观看| 国产色爽女视频免费观看| 久久这里只有精品中国| av黄色大香蕉| 黑人高潮一二区| 美女被艹到高潮喷水动态| 婷婷六月久久综合丁香| 国产高清三级在线| 久久人人爽人人爽人人片va| 小蜜桃在线观看免费完整版高清| 听说在线观看完整版免费高清| 亚洲成av人片在线播放无| 中文在线观看免费www的网站| 久久亚洲精品不卡| av在线老鸭窝| 国产黄a三级三级三级人| 91狼人影院| 亚洲七黄色美女视频| 欧美激情久久久久久爽电影| 毛片一级片免费看久久久久| .国产精品久久| 精品欧美国产一区二区三| 欧美+亚洲+日韩+国产| 久久精品国产亚洲网站| 国产av在哪里看| 国产精品久久久久久久久免| avwww免费| 国产久久久一区二区三区| 一进一出好大好爽视频| 成人综合一区亚洲| 国产一区二区激情短视频| 真人做人爱边吃奶动态| 色综合站精品国产| 国产午夜精品论理片| 一区二区三区四区激情视频 | 国产一级毛片七仙女欲春2| 亚洲欧美中文字幕日韩二区| 久久欧美精品欧美久久欧美| 免费看日本二区| 热99在线观看视频| 男女那种视频在线观看| 亚洲内射少妇av| 欧美高清成人免费视频www| 综合色丁香网| 亚洲经典国产精华液单| 高清毛片免费观看视频网站| 1000部很黄的大片| 韩国av在线不卡| 五月伊人婷婷丁香| 级片在线观看| 变态另类丝袜制服| 日韩在线高清观看一区二区三区| 一个人免费在线观看电影| 亚洲精品一卡2卡三卡4卡5卡| 插阴视频在线观看视频| 91久久精品电影网| aaaaa片日本免费| 亚洲精品成人久久久久久| 亚洲第一电影网av| АⅤ资源中文在线天堂| 国产精品人妻久久久久久| 岛国在线免费视频观看| 国产老妇女一区| 熟女人妻精品中文字幕| 国产精品福利在线免费观看| 亚洲美女搞黄在线观看 | 天堂网av新在线| 色综合站精品国产| 少妇熟女aⅴ在线视频| 18禁在线无遮挡免费观看视频 | 美女cb高潮喷水在线观看| 色综合色国产| av视频在线观看入口| 亚洲精品一卡2卡三卡4卡5卡| 免费观看的影片在线观看| 在线观看免费视频日本深夜| 性插视频无遮挡在线免费观看| 日本黄色视频三级网站网址| 三级毛片av免费| 亚洲乱码一区二区免费版| 精品99又大又爽又粗少妇毛片|