999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

丁酸甲酯單分子解離的非諧振效應

2016-12-29 05:42:36宋立國余憶玄林圣賢
物理化學學報 2016年11期
關鍵詞:效應

丁 楊 宋立國 余憶玄 姚 麗,* 林圣賢

(1大連海事大學輪機工程學院,遼寧大連116026;2臺灣交通大學應用化學系,臺灣新竹10764)

丁酸甲酯單分子解離的非諧振效應

丁 楊1宋立國1余憶玄1姚 麗1,*林圣賢2

(1大連海事大學輪機工程學院,遼寧大連116026;2臺灣交通大學應用化學系,臺灣新竹10764)

使用MP2/6-311++G(2d,2p)方法和基組,計算了丁酸甲酯單分子解離反應體系詳細的勢能面。應用RRKM理論,計算了在1000-5000 K的溫度范圍內(nèi)的正則系綜的速率常數(shù)。與此同時,在微正則系綜下,我們計算了溫度為1000-5000 K對應的能量從451.92到1519.52 kJ·mol-1的速率常數(shù)。計算結(jié)果表明反應通道2、4和5的非諧振效應比較明顯。因此對于丁酸甲酯單分子解離反應體系來說其非諧振效應是不能忽視的。

非諧振效應;單分子解離反應;RRKM理論;速率常數(shù)

1 Introduction

The utilization of petroleum-based fuels has been associated in recent years with various social and environmental issues including energy and national securities and air pollution by net CO2emissions which has been linked to climate change1.The extensive consumption of these fuels has motivated researchers to evaluate alternative solutions.In this sense,biofuels are alternative fuels suitable for being used in the current transport sector infrastructure,since they have similar physical properties than conventional petroleum-derived fossil fuels2.Biodiesel,consisting of long chain alkyl(e.g.,methyl,ethyl,propyl)esters,has emerged as a viable alternative to petroleum-based fuels.

Direct studies of typical biodiesel fuels(i.e.,methyl esters of fatty acids)are currently beyond our capabilities since the laboratory experiments would have to be carried out with complex and largely involatile mixtures.Methyl butanoate(MB),whose formula is CH3CH2CH2C(=O)OCH3,has been widely used as a biodiesel surrogate since it essentially possesses the chemical structure of long-chain alkyl esters including the methyl ester termination and a shorter,but similar,alkyl chain.Thus it is convenient to develop detailed reaction mechanisms for MB at a manageable size3-14.

In this follow-up paper,the rate constants for the decomposition reaction of MB as well as the transition state were calculated according to Rice-Ramsperger-Kassel-Marcus(RRKM)theory. Additionally the anharmonic effect was discussed herein.The Morse oscillators(MOs)were employed in the calculation for convenience.The RRKM theory15was previously applied to calculate the rate constants and the microcanonical and canonical cases showing similar results than those reported herein.

2 Calculation methodology and computational details

2.1 Ab initio calculations

The MP2 functional in conjunction with the 6-311++G(2d,2p) basis set were used to explore the geometry optimization of the reactants and transition states(TS)for the most important MB breakdown pathways.The vibrational frequency calculations were used to identify all of the stationary points as either minima(i.e., zero imaginary frequency)or transition states(i.e.,one imaginary frequency).Intrinsic reaction coordinate16(IRC)calculations at the same level were traced to confirm that the TS corresponded with the minima along the reaction pathways.Vibrational harmonic and anharmonic frequencies were also calculated at the same level. With the aim to obtain more accurate and reliable data,the singlepoint energies(SPEs)were recalculated by employing the CCSD (T)method with the 6-311++G(2d,2p)basis set.All the electronic structure calculations were carried out using the Gaussian 09 suite of programs17.

2.2 Microcanonical case

According to the famous equation of RRKM theory18,the microcanonical unimolecular reaction rate k(E)for a reaction at a given energy E and with an activation energy E≠can be expressed as follow:

where σ is the reaction degeneracy(herein σ=1),h is Plank′s constant,ρ(E)represents the total density of the states of the reactant at energy E,and W≠(E-E≠)stands for the total number of states for the transition state with an excess energy lower or equal to E-E≠,E is the total energy and E≠represents the activation energy.

The W(E)and ρ(E)function of can be generally defined by their respective definitions19-21:

H(E-Ei)denotes Heaviside function while Eirepresents the energy levels.

The Laplace transformation was used in the calculation of W(E) and ρ(E):

where β=1/kT,k is Boltzmann′s constant,T is the system temperature in K,and Q(β)is the partition function of the system. Therefore,in case Q(β)is acquired,with the use of above equations,W(E)and ρ(E)can be determined by employing the inverse Laplace transformation.

2.3 Canonical case

For a canonical system,the rate constant k(T)for the decomposition reaction can be calculated by the well-known equation of the transition state theory(TST)19,21-23.

where Q(T)and Q≠(T)represent the partition functions of the reactant and the activated complex,respectively.Thus:

where N stands for the number of the vibrational modes of the reactant.For each mode,qi≠(T)and qi(T)represent the vibrational partition functions of the activated complex and the reactant molecule for the ith single mode,respectively.

The aforementioned discussion demonstrates that the partition function has a significant weight in the calculation of k(E)and k(T). To calculate the partition function,the MO was taken as a simple form.And the energy of the ith vibrational mode can be calculated as follows:

where niand ωiare the vibrational quantum number and the frequency of the ith vibrational mode,respectively.χistands for the MO parameter and can be expressed as:

where Direpresents the well depth of the MO.In this study,χifor various molecules were obtained from the anharmonic frequencies which were calculated by the Gaussian 09 program.According to the RRKM theory,all the vibrational modes were treated as anharmonic MO.In the calculation of the density of states ρ(E),the harmonic and anharmonic degrees of freedom(DOF)of the reactant were calculated as 45(3N-6,N=17).For calculating the total number of states W(E),44 DOF were obtained excluding the harmonic and anharmonic imaginary frequencies for the transition states.The harmonic frequencies and χiwere chosen as effective dissociation energy parameters for the Morse potential in the calculations for each vibrational mode.

3 Results and discussion

The anharmonic and harmonic rate constants for six reaction channels were calculated for MB:

Fig.1 depicts with an illustration of the potential energy surfaces for the decomposition reaction(labeled as pathways 1-6)computed at the MP2/6-311++G(2d,2p)level of theory.The barrier heights were also reported with the CCSDT/6-311++G(2d,2p) method.The geometric and energetic parameters of the reactant and transition states for the unimolecular dissociation of the MB are summarized in Table 1.The optimized geometries of the reactants,transition states and products were optimized with a different method and levels as compared to Ref.24.The energy barrier of the decomposition reaction was higher as compared to previous work24.

As clearly seen from Fig.1,that CH3CH2CH2OOCH3reacts along six pathways which could be divided into three species:(i) C(8)―H(9)scission:H(9)from C(8)transfers to O(2);(ii)C(1)―O(3)scission,H(10)from C(1)transfers to O(3)or H(5)from C(4)transfers to C(1);(iii)C(8)―C(11)scission while H(15)from C(14)transfers to C(8)or O(2)or C(8)―C(1)scission and H(12) from C(11)transfers to C(1).

3.1C―H scission

To calculate the above mentioned energy as a function of the temperature,we applied the method by means of the relation between the total energy of a microcanonical system and the temperature of a canonical system18.

The energy in the microcanonical system can be calculated by equation(18),and it is listed in Table 2.

The MB reaction passes though TS2 and produces C2H5CH=C(OH)OCH3.The anharmonic and harmonic rate constants for the canonical system are summarized in Table 2 at the temperatures ranging from 1000 to 5000 K,with the energies being lower than the calculated activation energy,(i.e.317.67 kJ·mol-1).Thus,the rate constant in a microcanonical system at higher energy have to be calculated.Table 3 summarizes the harmonic and anharmonic rate constants of the microcanonical system as a function of the corresponding energy.

Fig.1 Potential energy surface(PES)schematic of the CH3CH2CH2C(O)OCH3dissociation reactions

Table 1 Parameters used in the rate constant calculations obtained from the MP2/6-311++G(2d,2p)calculations

From Table 2 and Fig.2,it is clear that both the harmonic and anharmonic rate constants increased with the temperature increasing.The rate constants for the reaction(Table 2)were plotted in Fig.2.The harmonic(from 2.35×10-4to 1.06×1010s-1)and the anharmonic rate constants(from 5.46×10-4to 1.34×1011s-1) were found to change with the temperature(from 1000 to 5000 K, respectively).The gap between anharnonic and harmonic rate constants changed with the increasing temperatures.When the temperature is 1000 K,the anharmonic rate constant(5.46×10-4s-1)is 2.20 times more than the harmonic one(2.35×10-4s-1),and the anharmonic rate constant(1.34×1011s-1)is 12.64 times thanharmonic one(1.06×1010s-1)at 5000 K.The harmonic and the anharmonic rate constants increased with the total energy for the microcanonical system(Table 3 and Fig.2).The harmonics rate constants(from 0.51×102to 0.56×1010s-1)and the anharmonic rate constants(from 1.46×102to 6.36×1010s-1)changed with the increasing energy(from 452.05 to 1519.56 kJ·mol-1).The gap between the anharnonic and harmonic rate constants changed with the total energy(anharmonic/harmonic rate constant ratios are 2.86 and 11.36 at 452.05 and 1519.56 kJ·mol-1,respectively).The ratio difference of canonical rate constant was compared with that of the microcanonical case.An increment of temperature or total energy resulted in significant anharmonic effects in both canonical and microcanonical systems.And there was a similar conclusion in Ref.25.

Table 2 Rate constants of the TS2 pathway at different temperatures for the canonical system

Table 3 Rate constants of the TS2 pathway at different energies for the microcanonical system

3.2C―O scission and C―H scission

The C―O scission can react with H from methyl group to O or C to form methanol and ketene or form formaldehyde and butyraldehyde.Two pathway reactions were denoted as TS3 and TS4.Similar to the TS2,the energy in the microcanonical system can be calculated through equation(18).The anharmonic and harmonic rate constants for the canonical system are summarized in Table 4 at temperatures ranging from 1000 to 5000 K.The energies are lower than the calculated activation energy(305.64 and 318.51 kJ·mol-1).Thus,the rate constants in a microcanonical system at higher energy values have to be calculated.Table 5 illustrates the harmonic and anharmonic rate constants of the microcanonical system with the corresponding energy.

Fig.2 Microcanonical and canonical rate constants for TS2

Similar to the TS2,it is clear that both the harmonic and anharmonic rate constants increase with the temperature for the reaction of TS3 and TS4(Table 4 and Figs.3-4).Both of the harmonic and anharmonic rate constants of TS3 and TS4 increasedsharply while increasing temperature from 1000 to 5000 K.The gap between anharnonic and harmonic rate constants were found to change with the temperatures(anharmonic/harmonic rate constant ratio=2.38,2.05 at 1000 K;5.86,30.67 at 5000 K,for TS3 and TS4,respectively).The harmonic and the anharmonic rate constants increased with the total energy(from 452.05 to 1519.56 kJ·mol-1)for the microcanonical system(Table 5 and Figs.3 and 4).The gap between the anharnonic and harmonic rate constants changed with the total energy(anharmonic/harmonic rate constant ratio=3.12,2.53 at 452.05 kJ·mol-1;6.46,26.63 at 1519.56 kJ·mol-1,for TS3 and TS4,respectively).The ratio differences of in the anharmonic and the harmonic rate constants in the canonical and microcanonical showed nearly similar results, Thus,with an increasing temperature or total energy,the anharmonic effect were not very pronounced in both canonical and microcanonical systems for TS3,more intense in both canonical and microcanonical systems for TS4.

Table 4 Rate constants of TS3 and TS4 pathways at different temperatures for the canonical system

Table 5 Rate constants of TS3 and TS4 pathways at different energies for the microcanonical system

Fig.3 Microcanonical and canonical rate constants for TS3

3.3C―C scission and C―H scission

The C―C scission can react with H from the methyl group to O or C to form ethene and methyl acetate,or propene and methyl, or ethene and 1-methoxy-ethenol.The three pathway reactions can be denoted as TS1,TS5 and TS6.The anharmonic and harmonic rate constants for the canonical system were summarized in Table 6 at the temperatures ranging from 1000 to 5000 K.The energies obtained were lower than the calculated activation energy(287.63, 415.50 and 441.08 kJ·mol-1).Thus,the rate constants at higher energy have to be calculated for the microcanonical system.

Fig.4 Microcanonical and canonical rate constants for TS4

Similar to theTS2,the reactions ofTS1,TS5 andTS6 proceeded such that both the harmonic and anharmonic rate constants increased with the temperature(Table 6 and Figs.5-7).The harmonics and anharmonic rate constants of TS1,TS5 and TS6 increased sharply while increasing temperature from 1000 to 5000 K.The gap between the anharnonic and harmonic rate constants changed with the temperature(anharmonic/harmonic rate constantratio=2.13,2.92,0.88 at 1000 K;6.49,32.00,1.30 at 5000 K,for TS1,TS5 and TS6,respectively).The harmonic and the anharmonic rate constants increased with the total energy from 452.05 to 1519.56 kJ·mol-1for the microcanonical system(Table 7 and Figs.5-7).The gap between the anharnonic and harmonic rate constants changed with the total energy(anharmonic/harmonic rate constant ratio=2.80,1.69,1.62 at 452.05 kJ·mol-1;7.49, 26.20,1.36 at 1519.56 kJ·mol-1,for TS1,TS5 and TS6,respectively).The ratio difference between the anharmonic and the harmonic rate constants in the canonical system was lower than that in the microcanonical system.Thus,the anharmonic effect was more intense in both canonical and microcanonical systems for TS5,and weaker in both canonical and microcanonical systems for TS1 and TS6.

Table 6 Rate constants of TS1,TS5 and TS6 pathways at different temperatures for the canonical system

Fig.5 Microcanonical and canonical rate constants for TS1

Additionally,we employed the following formula to obtain the tunneling probabilities for the unimolecular dissociation of MB26:

where

here,ωbis the magnitude of the imaginary frequency,V0is the barrier height relative to the reactant,and V1is the barrier height relative to the products.

Fig.6 Microcanonical and canonical rate constants for TS5

Fig.7 Microcanonical and canonical rate constants for TS6

Table 7 Rate constants of TS1,TS5 and TS6 pathways at different energies for the microcanonical system

Table 8 Tunneling probabilities for the CH3CH2CH2C(O)OCH3dissociation reactions

The tunneling probabilities increased with the total energy (Table 8).The values for the decomposition of MB in the barrier were obtained by MP2 methods.Note that the tunneling effect for the decomposition of MB was very small,which can be neglected in this work.

4 Conclusions

The anharmonic and harmonic rate constants of the unimolecular decomposition reaction of MB were calculated by using the RRKM theory.The reaction took place along with three kinds of pathways,including six reaction channels:(i)C―H scission: H from C transfers to O;(ii)C―O scission,H from C transfers to O or C;(iii)C―C scission and while H from C transfers to C or O.The rate constants of the reaction were evaluated by using the MP2/6-311++G(2d,2p)and CCSD(T)/6-311++G(2d,2p)methods in the temperature range of 1000-5000 K.TS1 showed the lowest energy barrier,and this reaction could thus be achieved.TS6 showed the highest energy barrier such that this reaction could not be reached.The anharmonic effect was represented in the Figs.2-7 and Tables 2-7.The difference between the harmonic and anharmonic rate constants increased with both the temperature and energy level.The anharmonic of the title reaction was also examined.Thus,the anharmonic rate constants were higher than the harmonic ones in both the microcanonical and the canonical systems,especially at high total energies and temperatures.The anharmonic effect played an important role in the unimolecular dissociation,such that the anharmonic effect could not be neglected.For the different models and vibrational states,the total number of states and density of states were counted,which affected the dissociation rate constant.For the first kind of reaction, it was a process of isomerization.With an increasing temperature or total energy,the anharmonic effect were more intense in bothcanonical and microcanonical systems for TS2.For the final kind of reaction,when one of the product was ethene for TS1 and TS6, with an increasing temperature or total energy,the anharmonic effect was not very pronounced in both canonical and microcanonical systems.These computational studies would be useful in providing a further insight into the chemical kinetics of MB,and further experimental studies are expected to be carried out with this reaction.

(1) Solomon,S.;Qin,D.;Manning,M.;Chen,Z.;Marquis,M.; Averyt,K.B.;Tignor M.;Miller,H.L.“IPCC,Climate Change 2007:The Physical Science Basis.Contribution of Working Group I to the FourthAssessment Report of the Intergovernmental Panel on Climate Change”;Cambridge University Press:Cambridge and New York,2007;p 996.

(2) Robert,L.H.;Roger,B.;Robert,W.AIChE J.2006,52(1),2. doi:10.1002/aic.10747

(3) Gail,S.;Thomson,M.J.;Sarathy,S.M.;Syed,S.A.;Dagaut, P.;Diévart,P.;Marchese,A.J.;Dryer,F.L.Proc.Combust.Inst 2007,31(1),305.doi:10.1016/j.proci.2006.08.051

(4) Metcalfe,W.K.;Dooley,S.;Curran,H.J.;Simmie,J.M.;El-Nahas,A.M.;Navarro,M.V.J.Phys.Chem.A 2007,111(19), 4001.doi:10.1021/jp067582c

(5) Huynh,L.K.;Violi,A.J.Org Chem.2008,73(1),94. doi:10.1021/jo701824n

(6) Huynh,L.K.;Lin,K.C.;Violi,A.J.Phys.Chem.A 2008,112 (51),13470.doi:10.1021/jp804358r

(7) Westbrook,C.K.;Pitz,W.J.;Curran,H.J.J.Phys.Chem.A 2006,110(21),6912.doi:10.1021/jp056362g

(8) Herbinet,O.;Pitz,W.J.;Westbrook,C.K.Combustion and Flame 2008,154(4),507.doi:10.1016/j. combustflame.2008.03.003

(9) Hill,J.;Nelson,E.;Tilman,D.;Polasky,S.;Tiffany,D.Proc. Natl.Acad.Sci.2006,103(30),11206.doi:10.1073/ pnas.0604600103

(10) Farooqa,A.;Davidson,D.F.;Hanson,R.K.;Huynh,L.K.; Violi,A.Proc.Combust.Inst.2009,32,247.doi:10.1016/j. proci.2008.06.084

(11) Dooley,S.;Curran,H.J.;Simmie,J.M.Combustion and Flame 2008,153(1-2),2.doi:10.1016/j.combustflame.2008.01.005

(12) Fisher,E.M.;Pits,W.J.;Curran,H.J.;Westbrook,C.K.Proc. Combust.Inst.2000,28,1579.

(13) Ali,M.A.;Violi,A.J.Org.Chem.2013,78(12),5898. doi:10.1021/jo400569d

(14) Huynh,L.K.;Violi,A.J.Org.Chem.2008,73(1),94. doi:10.1021/jo701824n

(15) (a)Yao,L.;Mebel,A.M.;Lu,H.F.;Neusser,H.J.;Lin,S.H. J.Phys.Chem.A 2007,111(29),6722.doi:10.1021/jp069012i (b)Yao,L.;Liu,Y.L.;Lin,S.H.Mod.Phys.Lett.B 2008,22 (31),3043.doi:10.1142/S0217984908017552 (c)Yao,L.;Lin,S.H.Sci.China.Ser.B 2008,51(12),1146. doi:10.1007/s11426-008-0125-1 (d)Yao,L.;He,R.X.;Mebel,A.M.;Lin,S.H.Chem.Phys. Lett.2009,470(4-6),210.doi:10.1016/j.cplett.2009.01.074 (e)Shao,Y.;Yao,L.;Lin,S.H.Chem.Phys.Lett.2009,478 (4-6),277.doi:10.1016/j.cplett.2009.07.051 (f)Yao,L.;Mebel,A.M.;Lin,S.H.J.Phys.Chem.A 2009,113 (52),14664.doi:10.1021/jp9044379 (g)Shao,Y.;Yao,L.;Mao,Y.C.;Zhong,J.J.Chem.Phys.Lett. 2010,501(1-3),134.doi:10.1016/j.cplett.2010.10.041 (h)Gu,L.Z.;Yao,L.;Shao,Y.;Yung,K.;Zhong,J.J.J.Theor. Comput.Chem.2010,9(1),813.doi:10.1142/ S0219633610006006 (i)Gu,L.Z.;Yao,L.;Shao,Y.;Liu,W.;Gao,H.Mol.Phys. 2011,109(16),1983.doi:10.1080/00268976.2011.602648 (j)Li,Q.;Xia,W.W.;Yao,L.;Shao,Y.Can.J.Chem.2012,90 (10),186.doi:10.1139/v11-137 (k)Li,Q.;Yao,L.;Shao,Y.CheM 2012,2(12),1. doi:10.5618/chem.2012.v2.n1.1 (l)Li,Q.;Yao,L.;Shao,Y.;Yang,K.J.Chin.Chem.Soc.2014, 61(3),309.doi:10.1002/jccs.201300277

(16) Gonzalez,C.;Schlegel,H.B.J.Chem.Phys.1989,90(4). 2154.doi:10.1063/1.456010

(17) Frisch,M.J.;Trucks,G.W.;Schlegel,H.B.;et al.Gaussian 09, Revision C.02;Gaussian,Inc.:Wallingford,CT,2009.

(18) Steinfeld,J.I.;Francisco,J.S.;Hase,W.L.Chemical Kinetics and Dynamic;Prentice-Hall:Englewood Cliffs,NJ,1989.

(19) (a)Forst,W.;Prasil,Z.J.Chem.Phys.1970,53(12),3065. doi:10.1063/1.1674450 (b)Forst,W.Chem.Rev.1971,71(4),339.doi:10.1021/ cr60272a001 (c)Forst,W.Theory of Unimolecular Reactions;Academic Press:New York,1973.

(20) Hoare,M.R.;Ruijgrok,T.W.J.Chem.Phys.1970,52(1),113. doi:10.1063/1.1672655

(21) Eyring,H.;Lin,S.H.;Lin,S.M.Basic Chemical Kinetics; AWiley-interscience Publication:New York,1980.

(22) Baer,T.;Hase,W.L.Unimolecular Reaction Dynamic:Theory and Experiment;Oxford University Press:New York,1996.

(23) Gilbert,R.G.;Smith,S.C.Theory of Unimolecular and Recombination Reactions;Blackwell:Oxford,1990.

(24) El-Nahas,A.M.;Navarro,M.V.;Simmie,J.M.;Bosselli,J.W.; Curran,H.J.;Dooley,S.;Metcalfe,W.J.Phys.Chem.A 2007, 111(19),3727.doi:10.1021/jp067413s

(25) Zhang,L.W.;Yao,L.;Li,Q.;Wang,G.Q.;Lin,S.H. Molecular Physics 2014,112(21),2853.doi:10.1080/ 00268976.2014.915066

(26) Miller,W.H.J.Am.Chem.Soc.1979,101(23),6810. doi:10.1021/ja00517a004

Anharmonic Effect of the Decomposition Reaction of Methyl Butanoate

DING Yang1SONG Li-Guo1YU Yi-Xuan1YAO Li1,*LIN Sheng-Hsien2
(1Marine Engineering College,Dalian Maritime University,Dalian 116026,P.R.China;2Department of Applied Chemistry,National Chiao-Tung University,Hsin-chu 10764,Taiwan,P.R.China)

In this paper,we have used the MP2/6-311++G(2d,2p)method to conduct a detailed investigation of the potential energy surface for the unimolecular dissociation reaction of methyl butanoate(MB).We have also used the Rice-Ramsperger-Kassel-Marcus(RRKM)theory to calculate the rate constants of the canonical and microcanonical systems at temperatures and total energies ranging from 1000 to 5000 K and 451.92 to 1519.52 kJ·mol-1,respectively.The results indicated that there was an obvious anharmonic effect for the TS2, TS4 and TS5 pathways,and that this effect was too pronounced to be neglected for the unimolecular dissociation reactions of MB.

Anharmonic effect;Unimolecular decomposition reaction;RRKM theory;Rate constant

O643

10.3866/PKU.WHXB201607212

Received:March 16,2016;Revised:July 20,2016;Published online:July 21,2016.

*Corresponding author.Email:yaoli@dlmu.edu.cn;Tel:+86-13130432506.

The project was supported by the Major Research Plan of the National Natural Science Foundation of China(91441132)and Fundamental Research Funds for the Central Universities,China(3132016127,3132016326).

國家自然科學基金(91441132)和中央高校基本科研業(yè)務費專項資金(3132016127,3132016326)資助項目

猜你喜歡
效應
鈾對大型溞的急性毒性效應
懶馬效應
場景效應
雨一直下,“列車效應”在發(fā)威
科學大眾(2020年17期)2020-10-27 02:49:10
決不能讓傷害法官成破窗效應
紅土地(2018年11期)2018-12-19 05:10:56
死海效應
應變效應及其應用
福建醫(yī)改的示范效應
福建醫(yī)改的示范效應
偶像效應
主站蜘蛛池模板: 美女高潮全身流白浆福利区| 日本成人精品视频| 91青青草视频| 国产精品无码影视久久久久久久| 欧美精品高清| 激情乱人伦| 免费在线a视频| a欧美在线| 国产精品片在线观看手机版| 国产日韩精品一区在线不卡 | 日韩毛片基地| 亚洲AⅤ无码日韩AV无码网站| 99在线视频免费| 九色在线观看视频| 久一在线视频| 国产色网站| 中国一级毛片免费观看| 国产精品青青| 国产人人乐人人爱| 日韩少妇激情一区二区| 日本人又色又爽的视频| 欧美a在线看| 伊人久久久久久久久久| 亚洲男人在线| 精品无码日韩国产不卡av| 91小视频版在线观看www| 欧美精品在线视频观看 | 国产99热| 97国产在线播放| 99久久精品视香蕉蕉| 欧美国产菊爆免费观看 | a级毛片免费看| 国产成人综合网在线观看| 亚洲视频在线观看免费视频| 日韩av手机在线| 亚洲国产日韩欧美在线| 国产一级二级在线观看| 女高中生自慰污污网站| 亚洲国产精品日韩av专区| 欧美成人综合在线| 亚洲制服中文字幕一区二区| 国产欧美日韩精品综合在线| 国模在线视频一区二区三区| 这里只有精品在线| 午夜高清国产拍精品| 精品国产自在现线看久久| 国产午夜无码片在线观看网站 | 亚洲成人在线网| 国产精品林美惠子在线观看| 国产裸舞福利在线视频合集| 国产在线视频二区| 国产xx在线观看| 毛片在线看网站| 亚洲av无码人妻| 国产视频入口| 欧美日本在线播放| 大香网伊人久久综合网2020| 国产拍揄自揄精品视频网站| 97青草最新免费精品视频| 情侣午夜国产在线一区无码| 中文字幕在线观看日本| 国产成人h在线观看网站站| 国产精品真实对白精彩久久| 日韩欧美国产精品| 亚洲天堂网在线视频| 人妻精品久久无码区| 91精品人妻互换| 人妻无码一区二区视频| 狠狠躁天天躁夜夜躁婷婷| 国产无吗一区二区三区在线欢| 亚洲男人的天堂网| 91在线日韩在线播放| 久久人搡人人玩人妻精品| 在线观看亚洲人成网站| 国产精品欧美日本韩免费一区二区三区不卡 | 国产国语一级毛片| 亚洲第一黄片大全| 77777亚洲午夜久久多人| 亚洲人在线| 国产自在线拍| 国产视频一二三区| 国产成人综合久久精品尤物|