• <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

    老司机影院毛片| 亚洲avbb在线观看| 免费在线观看视频国产中文字幕亚洲 | 亚洲少妇的诱惑av| 久久精品人人爽人人爽视色| 桃红色精品国产亚洲av| 黄色a级毛片大全视频| 午夜激情av网站| 国产在线观看jvid| 51午夜福利影视在线观看| 久9热在线精品视频| 久久热在线av| 黄网站色视频无遮挡免费观看| 午夜福利乱码中文字幕| 免费观看人在逋| 99国产综合亚洲精品| 午夜激情av网站| 色婷婷久久久亚洲欧美| 91九色精品人成在线观看| 正在播放国产对白刺激| 国产成人欧美在线观看 | √禁漫天堂资源中文www| 老司机福利观看| 桃花免费在线播放| 亚洲免费av在线视频| 亚洲天堂av无毛| 一级片免费观看大全| 久久精品成人免费网站| 丝袜人妻中文字幕| 99精品久久久久人妻精品| kizo精华| 亚洲精品第二区| 国产精品熟女久久久久浪| 免费在线观看日本一区| 黄色a级毛片大全视频| 亚洲国产欧美网| 欧美 亚洲 国产 日韩一| 免费在线观看日本一区| 亚洲一卡2卡3卡4卡5卡精品中文| 男女之事视频高清在线观看| 国产欧美日韩精品亚洲av| 久久久精品免费免费高清| 精品人妻一区二区三区麻豆| 国产亚洲精品第一综合不卡| 午夜精品国产一区二区电影| 999久久久精品免费观看国产| 手机成人av网站| 免费一级毛片在线播放高清视频 | 一本一本久久a久久精品综合妖精| 无遮挡黄片免费观看| 精品一区二区三卡| 久久中文看片网| 啦啦啦 在线观看视频| 国产精品一区二区在线观看99| 国产精品九九99| 国产精品影院久久| 久热爱精品视频在线9| 美女国产高潮福利片在线看| 一级片'在线观看视频| 久久免费观看电影| 久久精品国产亚洲av高清一级| 精品人妻熟女毛片av久久网站| 国产欧美亚洲国产| 成年美女黄网站色视频大全免费| 伊人亚洲综合成人网| 亚洲美女黄色视频免费看| 视频区图区小说| 中文精品一卡2卡3卡4更新| 亚洲精品国产色婷婷电影| 中文字幕av电影在线播放| 国产在视频线精品| 韩国高清视频一区二区三区| 十八禁高潮呻吟视频| 国产精品久久久人人做人人爽| 19禁男女啪啪无遮挡网站| 久久精品国产亚洲av高清一级| 精品欧美一区二区三区在线| 午夜两性在线视频| 十八禁网站网址无遮挡| 亚洲av男天堂| 亚洲va日本ⅴa欧美va伊人久久 | 精品国产一区二区久久| 免费在线观看视频国产中文字幕亚洲 | 亚洲精品日韩在线中文字幕| av国产精品久久久久影院| 午夜精品国产一区二区电影| 大陆偷拍与自拍| 国产97色在线日韩免费| 国产精品1区2区在线观看. | 成人av一区二区三区在线看 | 黄色视频在线播放观看不卡| 人人妻人人爽人人添夜夜欢视频| 久久久精品国产亚洲av高清涩受| av不卡在线播放| 巨乳人妻的诱惑在线观看| 满18在线观看网站| 欧美久久黑人一区二区| 狠狠精品人妻久久久久久综合| 日韩熟女老妇一区二区性免费视频| 日韩制服骚丝袜av| 亚洲自偷自拍图片 自拍| 美女中出高潮动态图| 国产精品久久久久久人妻精品电影 | 青春草亚洲视频在线观看| 日韩中文字幕欧美一区二区| 老司机影院毛片| 少妇的丰满在线观看| 色老头精品视频在线观看| 另类亚洲欧美激情| 一个人免费在线观看的高清视频 | 两人在一起打扑克的视频| 夜夜夜夜夜久久久久| 99久久精品国产亚洲精品| 亚洲色图 男人天堂 中文字幕| 国产色视频综合| 中亚洲国语对白在线视频| 人人妻人人澡人人看| 亚洲第一欧美日韩一区二区三区 | 免费少妇av软件| 日韩欧美一区二区三区在线观看 | 久久久久国产一级毛片高清牌| 精品一区二区三区av网在线观看 | 50天的宝宝边吃奶边哭怎么回事| 在线 av 中文字幕| 不卡一级毛片| 国产无遮挡羞羞视频在线观看| 亚洲欧美激情在线| 又紧又爽又黄一区二区| 国产成人精品在线电影| 欧美97在线视频| 蜜桃在线观看..| 男女床上黄色一级片免费看| 天天添夜夜摸| 十八禁网站免费在线| 成人免费观看视频高清| 久久免费观看电影| 在线观看免费午夜福利视频| 婷婷丁香在线五月| 国产99久久九九免费精品| 亚洲成人手机| 亚洲精品国产av成人精品| 婷婷色av中文字幕| 久久亚洲精品不卡| 天天添夜夜摸| 国产一区二区三区av在线| 亚洲欧洲日产国产| 99久久综合免费| 悠悠久久av| 热99久久久久精品小说推荐| 久久国产精品大桥未久av| 免费高清在线观看视频在线观看| 最新在线观看一区二区三区| 啦啦啦视频在线资源免费观看| 国产视频一区二区在线看| 久久久久精品人妻al黑| 亚洲精品国产一区二区精华液| 一区二区三区激情视频| 性色av乱码一区二区三区2| 伦理电影免费视频| 97在线人人人人妻| 久久久久久久久久久久大奶| 少妇猛男粗大的猛烈进出视频| 久久久国产欧美日韩av| 母亲3免费完整高清在线观看| 精品一区二区三区av网在线观看 | 99国产极品粉嫩在线观看| 69av精品久久久久久 | 飞空精品影院首页| 日韩欧美免费精品| 精品一区在线观看国产| av在线app专区| 国产精品 国内视频| av国产精品久久久久影院| 一个人免费看片子| 久久香蕉激情| 岛国在线观看网站| 午夜福利免费观看在线| 国产高清国产精品国产三级| 精品人妻一区二区三区麻豆| 国产av又大| 搡老熟女国产l中国老女人| 亚洲精品久久成人aⅴ小说| 亚洲色图 男人天堂 中文字幕| 国产在视频线精品| 国产伦理片在线播放av一区| 中文字幕色久视频| 亚洲国产成人一精品久久久| 久久免费观看电影| 亚洲精品国产av成人精品| 1024香蕉在线观看| 婷婷丁香在线五月| 99久久人妻综合| 伊人亚洲综合成人网| av国产精品久久久久影院| 巨乳人妻的诱惑在线观看| 日本av免费视频播放| 国产在视频线精品| 国产精品九九99| 成人三级做爰电影| 人人妻,人人澡人人爽秒播| 国产精品久久久久久人妻精品电影 | 免费一级毛片在线播放高清视频 | 男女边摸边吃奶| 日日爽夜夜爽网站| 1024香蕉在线观看| 人妻久久中文字幕网| 欧美 日韩 精品 国产| 成年美女黄网站色视频大全免费| 久久综合国产亚洲精品| 十八禁高潮呻吟视频| 日韩大码丰满熟妇| 不卡av一区二区三区| 精品少妇黑人巨大在线播放| 爱豆传媒免费全集在线观看| 中文字幕制服av| cao死你这个sao货| 啦啦啦视频在线资源免费观看| 一区福利在线观看| 一本—道久久a久久精品蜜桃钙片| 国产视频一区二区在线看| 老司机影院毛片| 国产精品一区二区在线不卡| 亚洲国产毛片av蜜桃av| 99精品久久久久人妻精品| 一边摸一边抽搐一进一出视频| 国产欧美日韩精品亚洲av| av网站在线播放免费| 三上悠亚av全集在线观看| 69av精品久久久久久 | 嫩草影视91久久| 麻豆av在线久日| 99久久99久久久精品蜜桃| 18禁国产床啪视频网站| 欧美日韩亚洲国产一区二区在线观看 | 国产一区二区在线观看av| 侵犯人妻中文字幕一二三四区| 色94色欧美一区二区| 黑人猛操日本美女一级片| 亚洲精品一区蜜桃| 亚洲人成电影免费在线| 亚洲三区欧美一区| av超薄肉色丝袜交足视频| 亚洲av成人一区二区三| 久久久水蜜桃国产精品网| 大片电影免费在线观看免费| 国产野战对白在线观看| 女警被强在线播放| 又紧又爽又黄一区二区| 纵有疾风起免费观看全集完整版| 精品少妇内射三级| 国产欧美日韩精品亚洲av| 一级片'在线观看视频| 伊人久久大香线蕉亚洲五| 免费看十八禁软件| 丝袜脚勾引网站| 亚洲av电影在线进入| 中国国产av一级| videos熟女内射| 国产av又大| av免费在线观看网站| 亚洲精品粉嫩美女一区| 正在播放国产对白刺激| 国产一区有黄有色的免费视频| av福利片在线| 美女视频免费永久观看网站| 青草久久国产| 一级黄色大片毛片| av网站免费在线观看视频| 亚洲av成人一区二区三| 性高湖久久久久久久久免费观看| 可以免费在线观看a视频的电影网站| 男人爽女人下面视频在线观看| 成年人免费黄色播放视频| av一本久久久久| 波多野结衣av一区二区av| 国产欧美日韩一区二区精品| 动漫黄色视频在线观看| 亚洲人成电影观看| 97精品久久久久久久久久精品| 精品国产乱码久久久久久男人| 国产99久久九九免费精品| 午夜福利乱码中文字幕| 2018国产大陆天天弄谢| 国产成人系列免费观看| 日本vs欧美在线观看视频| 亚洲三区欧美一区| 亚洲精品久久久久久婷婷小说| 欧美 亚洲 国产 日韩一| 国产亚洲欧美在线一区二区| 精品第一国产精品| 久久久欧美国产精品| 老司机午夜福利在线观看视频 | 亚洲国产欧美日韩在线播放| 国产黄频视频在线观看| 欧美成人午夜精品| 麻豆av在线久日| 不卡av一区二区三区| 国产亚洲欧美在线一区二区| 中文字幕av电影在线播放| 国产伦理片在线播放av一区| 91大片在线观看| av欧美777| 成人国产av品久久久| 午夜福利乱码中文字幕| 国产男女内射视频| 精品一区二区三卡| 久久久国产一区二区| 精品人妻一区二区三区麻豆| 日本撒尿小便嘘嘘汇集6| 飞空精品影院首页| 午夜精品久久久久久毛片777| 色老头精品视频在线观看| 国产精品av久久久久免费| 亚洲精品久久午夜乱码| 十八禁网站免费在线| 搡老乐熟女国产| av欧美777| 精品久久蜜臀av无| 国产成人系列免费观看| 日韩中文字幕视频在线看片| 巨乳人妻的诱惑在线观看| e午夜精品久久久久久久| 天天躁夜夜躁狠狠躁躁| 久9热在线精品视频| 天堂俺去俺来也www色官网| 男女边摸边吃奶| 手机成人av网站| 久久人妻福利社区极品人妻图片| 精品高清国产在线一区| 91成年电影在线观看| 日日摸夜夜添夜夜添小说| 日本精品一区二区三区蜜桃| 91老司机精品| 日本黄色日本黄色录像| 日日摸夜夜添夜夜添小说| 老汉色av国产亚洲站长工具| 久久青草综合色| 国产91精品成人一区二区三区 | 老司机在亚洲福利影院| 纵有疾风起免费观看全集完整版| 国产精品久久久久久精品电影小说| 中文字幕另类日韩欧美亚洲嫩草| 亚洲精品国产av成人精品| 天天添夜夜摸| 亚洲欧美清纯卡通| 丝袜美腿诱惑在线| 日韩免费高清中文字幕av| 性高湖久久久久久久久免费观看| 99久久综合免费| 视频区欧美日本亚洲| 看免费av毛片| 日韩,欧美,国产一区二区三区| 丰满迷人的少妇在线观看| 欧美午夜高清在线| 中国国产av一级| kizo精华| 我的亚洲天堂| 大香蕉久久网| 精品卡一卡二卡四卡免费| 久久精品人人爽人人爽视色| 欧美日韩国产mv在线观看视频| 精品久久久久久久毛片微露脸 | 国产91精品成人一区二区三区 | 动漫黄色视频在线观看| 免费观看a级毛片全部| 国产无遮挡羞羞视频在线观看| 两性夫妻黄色片| 性色av一级| 一区二区三区精品91| 麻豆乱淫一区二区| 中文欧美无线码| 建设人人有责人人尽责人人享有的| 亚洲人成电影观看| 欧美日韩黄片免| 久久精品亚洲av国产电影网| 成人av一区二区三区在线看 | 日本欧美视频一区| 色综合欧美亚洲国产小说| 亚洲欧美色中文字幕在线| 精品久久蜜臀av无| 国产精品香港三级国产av潘金莲| 国产精品 欧美亚洲| 久久精品aⅴ一区二区三区四区| 国产精品免费大片| 热99久久久久精品小说推荐| 久久人人爽av亚洲精品天堂| 亚洲精品美女久久av网站| 日日爽夜夜爽网站| 黄色视频,在线免费观看| 国产亚洲精品一区二区www | 成人亚洲精品一区在线观看| 91九色精品人成在线观看| 香蕉丝袜av| 99国产精品一区二区三区| 亚洲国产精品一区二区三区在线| 欧美日韩av久久| 久久天躁狠狠躁夜夜2o2o| 黑人巨大精品欧美一区二区mp4| 最新的欧美精品一区二区| 秋霞在线观看毛片| 麻豆国产av国片精品| 久久久久久久大尺度免费视频| 亚洲一卡2卡3卡4卡5卡精品中文| 99久久人妻综合| 国产成人免费观看mmmm| 青春草视频在线免费观看| √禁漫天堂资源中文www| 欧美午夜高清在线| 精品少妇久久久久久888优播| 黄色视频不卡| 天天躁日日躁夜夜躁夜夜| 免费女性裸体啪啪无遮挡网站| 欧美精品啪啪一区二区三区 | 免费在线观看视频国产中文字幕亚洲 | 中文字幕av电影在线播放| 日韩熟女老妇一区二区性免费视频| 精品卡一卡二卡四卡免费| 香蕉丝袜av| 欧美精品啪啪一区二区三区 | 欧美人与性动交α欧美软件| 交换朋友夫妻互换小说| 秋霞在线观看毛片| 精品熟女少妇八av免费久了| av在线app专区| 欧美一级毛片孕妇| 另类亚洲欧美激情| av国产精品久久久久影院| 90打野战视频偷拍视频| 亚洲欧美清纯卡通| 中亚洲国语对白在线视频| 日韩中文字幕欧美一区二区| 欧美日韩一级在线毛片| 日韩大码丰满熟妇| 老司机靠b影院| 18禁黄网站禁片午夜丰满| 女人高潮潮喷娇喘18禁视频| 精品久久蜜臀av无| 欧美午夜高清在线| 日日夜夜操网爽| 亚洲成人国产一区在线观看| www.精华液| 久久毛片免费看一区二区三区| 搡老乐熟女国产| 久久精品熟女亚洲av麻豆精品| av片东京热男人的天堂| 午夜两性在线视频| 国产国语露脸激情在线看| 涩涩av久久男人的天堂| 午夜久久久在线观看| 午夜福利,免费看| 国产1区2区3区精品| 国产有黄有色有爽视频| 一级a爱视频在线免费观看| 久久精品亚洲av国产电影网| 老熟女久久久| 国产亚洲欧美精品永久| 亚洲av美国av| 大片免费播放器 马上看| 日韩欧美免费精品| 亚洲欧美日韩另类电影网站| 国产熟女午夜一区二区三区| videos熟女内射| 免费久久久久久久精品成人欧美视频| 人人妻人人爽人人添夜夜欢视频| 王馨瑶露胸无遮挡在线观看| 欧美精品人与动牲交sv欧美| 欧美亚洲日本最大视频资源| 操美女的视频在线观看| 久久精品国产亚洲av香蕉五月 | 丰满人妻熟妇乱又伦精品不卡| 在线精品无人区一区二区三| 亚洲欧美一区二区三区久久| av超薄肉色丝袜交足视频| 亚洲精品中文字幕在线视频| e午夜精品久久久久久久| av欧美777| 人人妻人人添人人爽欧美一区卜| 成在线人永久免费视频| 国产精品久久久久久精品电影小说| 久久性视频一级片| 极品少妇高潮喷水抽搐| 精品一区在线观看国产| 国产精品.久久久| 亚洲va日本ⅴa欧美va伊人久久 | 777久久人妻少妇嫩草av网站| 亚洲五月婷婷丁香| 精品亚洲乱码少妇综合久久| 亚洲国产中文字幕在线视频| 国产在线视频一区二区| 不卡av一区二区三区| 亚洲国产欧美一区二区综合| 国产一级毛片在线| 国产精品亚洲av一区麻豆| 人妻 亚洲 视频| xxxhd国产人妻xxx| 成年美女黄网站色视频大全免费| 性色av一级| 日本撒尿小便嘘嘘汇集6| 久久久精品免费免费高清| 窝窝影院91人妻| 国产免费福利视频在线观看| 天堂俺去俺来也www色官网| 日韩制服骚丝袜av| 欧美变态另类bdsm刘玥| 少妇猛男粗大的猛烈进出视频| 高清黄色对白视频在线免费看| 国产无遮挡羞羞视频在线观看| e午夜精品久久久久久久| 在线永久观看黄色视频| 久久国产亚洲av麻豆专区| 国产成人精品久久二区二区91| 动漫黄色视频在线观看| 日韩大片免费观看网站| 久久中文看片网| 欧美日韩国产mv在线观看视频| 亚洲一区二区三区欧美精品| www.自偷自拍.com| √禁漫天堂资源中文www| 国产精品九九99| 在线 av 中文字幕| 美女主播在线视频| 日韩 欧美 亚洲 中文字幕| 高潮久久久久久久久久久不卡| 丝袜美腿诱惑在线| 天堂俺去俺来也www色官网| 亚洲av男天堂| 亚洲,欧美精品.| 咕卡用的链子| 欧美亚洲 丝袜 人妻 在线| 97在线人人人人妻| 亚洲av电影在线观看一区二区三区| 黑丝袜美女国产一区| 精品少妇一区二区三区视频日本电影| 亚洲伊人久久精品综合| 宅男免费午夜| 国产成+人综合+亚洲专区| 久久精品亚洲av国产电影网| 无遮挡黄片免费观看| 久久人人爽av亚洲精品天堂| 性色av乱码一区二区三区2| 超碰97精品在线观看| 日韩大片免费观看网站| 亚洲成av片中文字幕在线观看| 亚洲国产欧美网| 又大又爽又粗| 国产亚洲av片在线观看秒播厂| 9色porny在线观看| 中文字幕另类日韩欧美亚洲嫩草| 国产无遮挡羞羞视频在线观看| 成年女人毛片免费观看观看9 | 久久久久久人人人人人| www.精华液| 国产免费现黄频在线看| 精品少妇一区二区三区视频日本电影| 少妇人妻久久综合中文| 亚洲视频免费观看视频| 日韩免费高清中文字幕av| 日本av免费视频播放| 巨乳人妻的诱惑在线观看| 国产av国产精品国产| 交换朋友夫妻互换小说| 国产深夜福利视频在线观看| 伦理电影免费视频| 少妇裸体淫交视频免费看高清 | 无遮挡黄片免费观看| 可以免费在线观看a视频的电影网站| 搡老岳熟女国产| 免费av中文字幕在线| 在线永久观看黄色视频| 丰满迷人的少妇在线观看| 99国产精品免费福利视频| 久久人人爽av亚洲精品天堂| 国产在线视频一区二区| 国产免费视频播放在线视频| 国产精品自产拍在线观看55亚洲 | 成人亚洲精品一区在线观看| 国产成人a∨麻豆精品| 人人妻人人澡人人爽人人夜夜| 午夜福利视频精品| 少妇裸体淫交视频免费看高清 | 十分钟在线观看高清视频www| 亚洲男人天堂网一区| 国产在视频线精品| 别揉我奶头~嗯~啊~动态视频 | av又黄又爽大尺度在线免费看| 日韩电影二区| 亚洲欧洲精品一区二区精品久久久| 日本精品一区二区三区蜜桃| 精品国内亚洲2022精品成人 | 成人黄色视频免费在线看| 黄色 视频免费看| 午夜老司机福利片| 精品亚洲成a人片在线观看| 美女主播在线视频| 丝瓜视频免费看黄片| 两人在一起打扑克的视频| 一个人免费看片子| 久久久国产一区二区| 国产精品香港三级国产av潘金莲| 亚洲国产中文字幕在线视频| 一区二区日韩欧美中文字幕| 成人影院久久| www.999成人在线观看| 一级毛片精品| 国产精品自产拍在线观看55亚洲 | av欧美777| 国产精品欧美亚洲77777| 国产欧美日韩一区二区三 | 国产精品一区二区精品视频观看| 高清黄色对白视频在线免费看| 手机成人av网站| 国产黄色免费在线视频| 久久久久国产精品人妻一区二区| 纵有疾风起免费观看全集完整版| 久久国产精品影院|