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

    Numerical simulation of carbon arc discharge for graphene synthesis without catalyst

    2021-07-07 02:40:50KaiYE葉凱DaZHANG張達(dá)KaiwenZHANG張楷文HaoyuWANG王灝宇FengLIANG梁風(fēng)WenhuiMA馬文會(huì)BinYANG楊斌andYongnianDAI戴永年
    Plasma Science and Technology 2021年7期
    關(guān)鍵詞:楊斌

    Kai YE(葉凱),Da ZHANG(張達(dá)),Kaiwen ZHANG(張楷文),Haoyu WANG(王灝宇),Feng LIANG(梁風(fēng)),2,Wenhui MA(馬文會(huì)),Bin YANG(楊斌) and Yongnian DAI(戴永年)

    1 Faculty of Metallurgical and Energy Engineering,Kunming University of Science and Technology,Kunming 650093,People’s Republic of China

    2 State Key Laboratory of Complex Nonferrous Metal Resources Clean Utilization,Kunming University of Science and Technology,Kunming 650093,People’s Republic of China

    Abstract In this study,graphene sheets are prepared under a hydrogen atmosphere without a catalyst,and the growth mechanism of graphene by direct current arc discharge is investigated experimentally and numerically.The size and layer numbers of graphene sheets increase with the arc current.Distributions of temperature,velocity,and mass fraction of carbon are obtained through numerical simulations.A high current corresponds to a high saturation temperature,evaporation rate,and mass density of carbon clusters.When the carbon vapor is saturated,the saturation temperatures are 3274.9,3313.9,and 3363.6 K,and the mass densities are 6.4×1022,8.42×1022,and 1.23×1023 m?3 under currents of 150,200,and 250 A,respectively.A hydrogen-induced marginal growth model is used to explain the growth mechanism.Under a high current,the condensation coefficient and van der Waals force increase owing to the higher saturation temperature and mass density of carbon clusters,which is consistent with experimental results.

    Keywords:graphene,arc discharge,numerical simulation,growth mechanism,saturation temperature,mass density

    1.Introduction

    Direct current(DC)arc discharge is an efficient and environment-friendly method for the controllable preparation of graphene[1–3].The products yielded contain high purity and few layers with a high preparation efficiency[4,5].Currently,DC arc discharge for the synthesis of graphene includes catalyst-based and catalyst-free methods.However,the large-scale fabrication of graphene sheets using the DC arc discharge method[6,7],which is controlled by parameters such as the atmosphere,pressure,current,and catalyst[3,5,8–10],remains challenging.To obtain the best synthesis conditions,the growth mechanism of graphene must be further investigated.

    Because the DC arc discharge method involves hightemperature and complex electromagnetic interactions[11,12],the growth process of graphene is difficult to explain using only limited experimental results.Numerical simulation is effective for understanding the physical characteristics of arc plasma and has been adopted to explain the growth mechanism of carbon nanotubes,thereby enabling the growth mechanism of graphene to be elucidated[13,14].For instance,Tewariet al[15]presented a growth model of a graphene–carbon nanotube hybrid under catalyst conditions via numerical simulations.The growth of the graphene–carbon nanotube hybrid occurred after the fragmentation of catalyst particles.It was discovered that the size of the products was affected by the ionization degree,gas type,and input power.High concentrations and ionization of hydrocarbons were more conducive to the growth of the graphene–carbon nanotube hybrid.Kundrapuet al[16]obtained the distribution of temperature and species density by simulating a carbon arc.The growth region of the carbon nanotube was determined based on catalyst clusters,which corresponded to a temperature range between the solidification point and the condensation point of nickel.Based on the above-mentioned numerical simulation studies involving the use of catalysts,it was discovered that the growth of graphene was governed by the nucleation of catalyst particles,which were easy to control.Furthermore,the growth temperature of graphene was reduced owing to catalytic DC arc discharge[16],thereby resulting in a low cost for large-scale synthesis.

    However,in the growth of graphene,an obstacle to the development of catalytic DC arc discharge is that the final product comprises graphene and the catalyst.These elements must be separated via post-transfer or magnetic separation[2].Cracks on the surface and interface,contamination,and wrinkles are present in the separation process[17],resulting in the poor quality of graphene.As for catalyst-free DC arc discharge in a reactive atmosphere,the introduction of impurities and damage to the graphene lamellar structure is prevented.This technique is more promising for the largescale preparation of high-quality graphene,although the growth mechanism is more complex and remains unclear.Studies regarding the growth mechanism of graphene without a catalyst via numerical simulation have not been reported hitherto.However,some inspirations can be obtained from relevant studies pertaining to carbon nanotubes and carbon nanoparticles.Mansouret al[18]simulated carbon arcs for the synthesis of single-walled carbon nanotubes.The energy equations for the electrons and heavy species were solved separately to obtain their temperatures.To obtain the number densities of multiple species,individual continuity equations,electron impact reactions,and chemical reactions among heavy species were coupled and solved.It was discovered that the change in the number density distributions was driven by the forms of anode ablation,which significantly affected the growth process.Yatomet al[19]investigated the growth region of carbon nanoparticles in an atmospheric pressure carbon arc.Combining with the two-dimensional computational fluid dynamic simulations and the thermodynamic modeling,the equation for the average diameter of carbon nanoparticles was established.It was proportional to the residence time in the growth region and the initial densities of particles at the starting location of condensation.In the above researches,the method for simulating the mass densities of carbon clusters has been reported,but the growth mechanism could not be explained using a specific growth model.Therefore,it is indispensable to investigate the growth mechanism of graphene without a catalyst and establish a growth model,which will contribute significantly to the development of graphene synthesized via DC arc discharge.

    In this study,graphene sheets were synthesized via DC arc discharge without a catalyst.The effect of the current on the size and number of layers of graphene is discussed.Moreover,the distributions of temperature,velocity,and mass density of carbon clusters were calculated via numerical simulations.Based on our proposed hydrogen-induced marginal growth model[20],the effects of saturation temperature and mass density of carbon clusters on the growth process of graphene are explained.

    2.Experiment

    Graphene sheets were synthesized using the DC arc discharge method.Figure 1 shows a schematic diagram of a DC arc discharge furnace.The two electrodes used were 99.999%pure graphite rods measuring 10 mm in diameter and 10 cm in length.The anode was placed on top of the water-cooled copper substrate,and the cathode was located above it.The cathode was sharpened,and the vertical length of the tip was 7 mm.Before discharge,the chamber was pumped into a vacuum.Subsequently,H2(99.99% purity)with a constant pressure of 70 kPa was injected as the buffer gas.The currents were set to 150,200,and 250 A to investigate their effects on the growth mechanism of graphene.By mechanically adjusting the height of the cathode,the distance between the two electrodes was maintained at 1 mm to ensure the stability of the carbon arc.After 10 min of evaporation,the products were obtained from the inner walls of the chamber.The morphology and nanostructure of the graphene sheets were analyzed using transmission electron microscopy(TEM,JEM-2100,JEOL,Japan)with an accelerating voltage of 200 kV.

    Figure 1.Schematic diagram of DC arc discharge furnace for preparing graphene sheets.

    3.Numerical simulation

    Numerical simulations were performed based on the synthesis conditions using COMSOL Multiphysics software.First,a free mesh was used to partition the mesh.Owing to the large current density gradient of the electrode boundaries and the velocity gradient of the chamber boundaries,fine meshes were required near these boundaries.Subsequently,the Joule heat and Lorentz force were set to the loads of the temperature and velocity fields,respectively.The electromagnetic,temperature,and velocity fields were synchronously coupled using the above-mentioned settings.To simplify the simulation calculation of plasma,the following conditions were assumed[18,21–27].(a)The arc plasma is under local thermodynamic equilibrium(LTE).Based on the same temperature for the electrons and heavy particles,the scale of the arc plasma is much larger than the Debye length.Therefore,the arc plasma is in a quasi-electric neutral state.(b)The laminar flow is in a steady-state and exhibits two-dimensional axial symmetry owing to the constant evaporation rate achieved by maintaining the distance between the two electrodes.(c)The plasma is optically thin,therefore,the radiation trapping in the plasma is negligible.(d)Gravity is negligible,and the flow direction of the arc plasma is only affected by the Lorentz force.(e)The effect of the electrode sheath is not considered.For the DC arc discharge,the gas dynamics of the arc discharge is primarily based on the assumption of LTE.

    3.1.Magnetohydrodynamic equations

    Temperature and velocity fields were obtained by solving the magnetohydrodynamic equations.The mass conservation,momentum conservation,and energy conservation equations are shown in equations(1)–(3)[16,23],respectively.Maxwell’s equations[28]are shown in equations(4)–(6)

    where ρ andare the density and velocity vector,respectively

    where μ is the viscosity coefficient;pis the pressure;andare the vectors of the current density and magnetic induction,respectively.The last term on the right-hand side of the equation represents the Lorentz force

    In equation(3),His the enthalpy;Cpis the specific heat at constant pressure;kBis the Boltzmann constant;kand σ are the thermal and electrical conductivities,respectively;Qrrepresents the radiation loss per unit volume.The last term on the left-hand side of the equation represents the Joule heat.To solve the terms of the Lorentz force and Joule heat,Maxwell’s equations are expressed as

    where Φ and μ0are the electric potential and vacuum permeability,respectively.

    A two-dimensional axisymmetric domain with a length and width of 70 and 61 mm,respectively,was used as the computational region of the plasma.The boundary conditions are illustrated in figure 2.In the electromagnetic field,the current density was applied to the anode surface,and the electric potential of the cathode tip was maintained at a constant value.All other boundaries were set as electric and magnetic insulation.In reality,the cathode tip is the closest to the anode and exhibits the largest potential gradient at its boundary layer,resulting in the maximum cathode current density.The cathode might ionize at high temperatures.This is consistent with the experimental observation that the cathode tip will be flattened slightly after discharge.However,anode evaporation is primarily considered in the preparation of graphene via arc discharge,and the products are obtained from only the inner walls instead of the cathode deposits[29].Based on this point,the ionization effect is only considered at the cathode tip,and the temperature of the cathode tip was set to 4000 K.The temperature of the cathode surface was a linear function from 4000 to 1000 K.Using the Langmuir model,a heat flux boundary was applied to the anode surface,where the anode sheath drop was disregarded and set to zero[12]:

    Figure 2.Axisymmetric domain of plasma region with boundary conditions.

    whereRAis the radius of the anode;Teis the electron temperature on the anode surface in eV;Φwis the work function of graphite;ΔH,Cp,andT0are the latent heat of evaporation,specific heat of carbon,and initial anode temperature,respectively.TA,which is the average temperature on the anode surface,is higher than the melting point of graphite,whereas τvaprepresents the mass flux of evaporation in kg m?2s?1.TAwas maintained at a constant value during the steady-state evaporation process of the anode.The temperatures of the other electrode and surrounding boundaries were set to 1000 and 600 K,respectively[30].

    For the boundaries of the velocity field,the anode surface was the inlet,and the surrounding boundaries were outlets with a constant pressure of 70 kPa.The other boundaries of the electrode and chamber were defined as walls.The inlet boundary condition is expressed by the mass flux of evaporation as follows[31]:

    Equation(8)is based on the Hertz–Knudsen relation,wheremCis the mass of a carbon atom,andpCis the pressure of carbon vapor in mmHg underTA,calculated using the Antoine equation as follows:

    3.2.Carbon vapor mass conservation equation

    Snyderet al[32]investigated the diffusion of hydrogen species in free-combustion arc discharge.It was shown that the formation temperature of hydrogen atoms by the ionization of H2molecules was 3000 K,and only hydrogen atoms existed in the temperature range of 4000–10 000 K.According to Byszewski[33],carbon atoms,carbon ions,C2molecules,and C3molecules dominated the arc plasma region.At temperatures above 4000 K,the carbon atom content was the highest among these carbon clusters.Therefore,hydrogen and carbon atoms were used to calculate the diffusion coefficients.

    The concentration of carbon was calculated using the equation of carbon vapor mass conservation as follows[34]:

    In equation(10),the three terms on the left-hand side of the equation represent the time,diffusion,and convection terms,respectively,whereYCandJCrepresent the mass fraction and mass flux of carbon,respectively.The carbon vapor source term(mass per unit volume and time)can be written as[34]

    In the equation above,ΔV/ΔSrepresents the width of the control body adjacent to the surface,and τdepis the deposition flux of carbon expressed as[35]

    The combined diffusion coefficient method was used to calculate the mass flux of carbon vapor.Among its simplified forms,a good approximation was indicated when the ‘combined ordinary diffusion coefficient mass fraction gradient’approach was used[34].In this method,only the ordinary diffusion coefficient was calculated.The pressure,electric field,and temperature diffusion coefficients were not considered.Therefore,the mass flux of carbon vapor can be written as[34,36]

    wherenis the number density,mHthe mass of a hydrogen atom,andxHthe mole fraction of hydrogen atoms.is the ordinary diffusion coefficient between carbon and hydrogen atoms;it can be expressed by the simplified form of Fick’s law,i.e.

    In equation(14),MC,MH,VC,andVHrepresent the molar masses and molecular diffusion volumes of carbon and hydrogen atoms,respectively.In this study,VC=1.59×10?5m3mol?1andVH=2.31×10?6m3mol?1.

    For the boundary condition,the anode surface was maintained at a constant mass fraction as follows[37]:

    wherepatmis the atmospheric pressure,andMH2is the molar mass of H2molecules.

    4.Results and discussion

    4.1.TEM images

    Figure 3 shows the TEM images of graphene sheets synthesized via DC arc discharge at 150,200,and 250 A.The products exhibited a two-dimensional lamellar structure.The size of the graphene sheets was calculated from the statistical results by counting 200 graphene sheets measuring 50–200 nm from 20 TEM images.The average sizes of the graphene sheets under 150,200,and 250 A were 93.5,119.6,and 138.0 nm,respectively.The rich-wrinkled morphology was located at the edges of the graphene sheets,which indicated that the stacking process occurred at these positions.Furthermore,the number of layers of graphene sheets can be estimated from the graphene edges in high resolution transmission electron microscopy images,and they were primarily 4–7,4–10,and 5–12 layers under 150,200,and 250 A,respectively.The results based on different currents indicate that a high current corresponds to graphene sheets of a larger size and higher layer number that tends to stack at the edges.

    Figure 3.TEM images of graphene sheets and HRTEM images of graphene edges under different currents:(a),(b)150 A,(c),(d)200 A,and(e),(f)250 A.

    4.2.Simulation results

    Figure 4 shows the distributions of temperature under 150,200,and 250 A,respectively.As shown in figure 4(a),the highest temperature of approximately 9250 K was recorded at the center of the two electrodes.Subsequently,the temperature decreased in the direction away from the electrodes owing to the decrease in the Joule heat with the current density.Moreover,the temperature distribution on the anode surface was not uniform,and the values ofTAwere 4119.6,4203.1,and 4299.7 K under 150,200,and 250 A,respectively.The corresponding mass fluxes of evaporation were 5.81,8.22,and 12.08 kg m?2s?1,respectively,as calculated using equation(8).Figure 5 shows the profiles of the velocity field at currents of 150,200,and 250 A,respectively.As shown in figure 5(a),the carbon vapor evaporated from the anode surface with a velocity of 90.2 m s?1and accelerated to 153.3 m s?1near the cathode.The velocity decreased gradually along the trajectory owing to expansion and reduced to zero in the walls of the chamber.To calculate the mass fraction,equation(10)is written in the form of the custom ‘partial differential equations’ module in COMSOL Multiphysics and then solved.Figure 6 shows the mass fraction profiles of carbon at currents of 150,200,and 250 A.The highest mass fractions of carbon were 66%,76%,and 83% on the anode surfaces,respectively.At a high current,the anode evaporation enhanced owing to the higher anode temperature,which corresponds to a higher initial condition of the mass fraction.In addition,the diffusion coefficient between carbon and hydrogen atoms was proportional to temperature,i.e.it was larger near the electrodes.Therefore,the mass fraction decreased rapidly at first in the direction away from the electrodes and then decreased slowly because of the lower diffusion coefficient.

    Figure 4.(a)–(c)Temperature profiles under 150,200,and 250 A,respectively.

    Figure 5.(a)–(c)Velocity profiles under 150,200,and 250 A,respectively.

    Figure 6.(a)–(c)Mass fraction distributions of carbon under 150,200,and 250 A,respectively.

    4.3.Growth region of graphene

    During the growth of graphene under hydrogen plasma,the enclosing of graphene sheets is prevented by the reaction between the hydrogen species and carbon clusters[8].Therefore,the growth of graphene begins at the saturation temperature of the carbon vapor and ends at the formation temperature of hydrogen atoms.S=pC/pSexpresses the supersaturation ratio of carbon,wherepCis the pressure of carbon in the carbon–hydrogen mixture,andpSis the saturation pressure of carbon vapor[38].The value ofpCwas provided by the simulation results of the carbon mass fraction,andpSwas obtained from equation(9).WhenS=1,the lamellar structure of graphene began to form through the condensation of carbon clusters.The saturation temperatures of carbon vapor were 3274.9,3313.9,and 3363.6 K at currents of 150,200,and 250 A,respectively.To provide important theoretical bases for explaining the growth mechanism of graphene,distributions of mass density of carbon clusters in the growth region under different arc currents are presented in figure 7.Because of the different streamline directions,the temperature in the horizontal direction from the center between two electrodes was selected to express the temperature range in the growth region.Figure 7 demonstrates the relationship among temperature,mass density of carbon clusters,and coordinates of streamlines on the horizontal axis;as shown in the figure,the mass density of the carbon clusters depended only on the temperature.Moreover,the temperature of the growth region ranged from the saturation temperature of carbon vapor to the termination temperature of 3000 K,and the range was wider for high currents.At the saturation of carbon vapor,the mass densities were 6.4×1022,8.42×1022,and 1.23×1023m?3under 150,200,and 250 A,respectively.It was clear that the mass density of the carbon clusters increased with the current on the same isothermal lines.

    Figure 7.Distributions of carbon clusters mass density in growth region under different arc currents.

    5.Growth mechanism of graphene

    Figure 8 shows the growth mechanism of graphene via DC arc discharge without a catalyst.The growth mechanism of graphene was investigated based on the hydrogen-induced marginal growth model[20].First,different types of carbon clusters(C+,C,C2,C3…)were formed by the evaporation of anode graphite.Subsequently,graphene islands with C–H bonds were generated at the edges owing to the reaction between the carbon clusters and hydrogen atoms.Because of the lower stability of C–H bonds in the high-temperature region,they were broken and replaced by C–C bonds.The size of graphene increased gradually during the cycling process of hydrogen transfer and carbon addition.In the longitudinal stacking process,the van der Waals force was the intermolecular force between the graphene islands.As the size of graphene increased,the van der Waals force was enhanced by the greater relative molecular mass of the graphene islands.According to the minimum energy principle,graphene islands tend to combine mutually[39].Therefore,graphene sheets with more stable structures were formed after the stacking process.

    Figure 8.Growth mechanism of graphene via DC arc discharge.

    The size and number of layers of the graphene sheets were dependent on the number density of the carbon clusters,which can be affected by the change in current.For the transverse condensation of graphene,the condensation coefficient can be expressed as[40]

    wherensis the number density of carbon clusters at saturation;v1andm1are the volume and mass of the monomer,respectively.It was discovered that the condensation coefficient was proportional to the saturation temperature of carbon vapor.For currents of 150,200,and 250 A,the condensation coefficients were 2.86×108,3.79×108,and 4.64×108,respectively.The transverse condensation process of graphene was enhanced by a higher condensation coefficient under a high current,thereby resulting in larger graphene sheets.This is consistent with the experimental results.Furthermore,the enhancement of the longitudinal stacking process under a high current is related to the increase in the van der Waals force.However,the van der Waals force is inversely proportional to the distance between the two molecules[41].Under a constant total volume,the distance between different graphene islands is closer for a larger carbon cluster number density.This corresponds to an increase in the van der Waals force.In addition,the van der Waals force is related to molecular deformation,which is proportional to the relative molecular mass[41].With an increase in current,the relative molecular mass is larger under a higher saturation temperature,thereby resulting in a greater van der Waals force and a larger number of graphene sheets.Therefore,a relationship between the carbon cluster number density and saturation temperature calculated via numerical simulation as well as the number of graphene layers was established.

    6.Conclusion

    Graphene sheets measuring 200–400 nm were synthesized via DC arc discharge without a catalyst.The formation of products with larger sizes and higher layer numbers was attributed to the enhanced evaporation process of anode graphite at high arc currents.Moreover,the saturation temperature and mass density of carbon clusters were obtained through numerical simulation.The hydrogeninduced marginal growth model was explained intuitively based on these parameters.For the high arc current,the greater saturation temperature and mass density of carbon clusters yielded a larger condensation coefficient and van der Waals force.Transverse condensation and longitudinal stacking were promoted,which is consistent with the larger size and higher layer numbers of graphene sheets demonstrated experimentally.

    Acknowledgments

    This work is financially supported by National Natural Science Foundation of China(No.11765010),the National Key Research and Development Program of China(No.2019YFC1907900),the Applied Basic Research Programs of Yunnan Provincial Science and Technology Department(No.202001AW070004),the Freely Exploring Fund for Academicians in Yunnan Province(No.2018HA006),the Key Laboratory of Resource Chemistry,Ministry of Education(No.KLRC_ME2001).

    猜你喜歡
    楊斌
    索要協(xié)調(diào)費(fèi)的金融“蛀蟲”
    公民與法治(2025年2期)2025-03-07 00:00:00
    半條毛毯見證80年革命友情
    大江南北(2024年11期)2024-12-31 00:00:00
    作品賞析
    雨花·藝術(shù) 楊斌作品
    雨花(2022年8期)2022-08-08 01:46:28
    信息化履職
    今日重慶(2020年1期)2020-09-10 07:22:44
    精神的漫游——楊斌工筆畫探微
    時(shí)間之象——楊斌自述
    虛懷以待
    寶藏(2018年10期)2018-10-18 02:29:24
    二婚丈夫中了千萬巨獎(jiǎng),福禍雙至打開婚姻結(jié)局
    追逐陽光的創(chuàng)業(yè)者楊斌昌
    大社會(huì)(2016年5期)2016-05-04 03:42:01
    少妇的逼水好多| av女优亚洲男人天堂| 男人爽女人下面视频在线观看| 下体分泌物呈黄色| 国产人妻一区二区三区在| 亚洲图色成人| 在线a可以看的网站| 国产成人精品福利久久| 日韩不卡一区二区三区视频在线| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 丰满人妻一区二区三区视频av| 免费不卡的大黄色大毛片视频在线观看| 午夜福利高清视频| 高清视频免费观看一区二区| av线在线观看网站| 欧美3d第一页| 精品亚洲乱码少妇综合久久| 久久久久精品性色| 99久国产av精品国产电影| 色播亚洲综合网| 好男人视频免费观看在线| 国产成人福利小说| 国内少妇人妻偷人精品xxx网站| 我的老师免费观看完整版| 午夜日本视频在线| 亚洲精品,欧美精品| 日韩不卡一区二区三区视频在线| 男人舔奶头视频| 欧美日韩精品成人综合77777| av又黄又爽大尺度在线免费看| 久久影院123| 久久人人爽av亚洲精品天堂 | 国模一区二区三区四区视频| 真实男女啪啪啪动态图| 高清午夜精品一区二区三区| av又黄又爽大尺度在线免费看| 禁无遮挡网站| 一级毛片aaaaaa免费看小| 久久亚洲国产成人精品v| 久久热精品热| av一本久久久久| 亚洲欧美一区二区三区国产| 日韩一区二区视频免费看| 1000部很黄的大片| 舔av片在线| 少妇 在线观看| a级一级毛片免费在线观看| av在线蜜桃| 啦啦啦在线观看免费高清www| 午夜福利视频1000在线观看| 欧美97在线视频| 成年女人看的毛片在线观看| 一级a做视频免费观看| 久久久久网色| 日韩欧美一区视频在线观看 | 看十八女毛片水多多多| 一级毛片 在线播放| 精品久久久久久久末码| 日本一本二区三区精品| 日韩欧美 国产精品| 身体一侧抽搐| 777米奇影视久久| 日本欧美国产在线视频| 欧美变态另类bdsm刘玥| 91aial.com中文字幕在线观看| 一级毛片电影观看| 三级国产精品片| kizo精华| 97人妻精品一区二区三区麻豆| 日本爱情动作片www.在线观看| 国产精品久久久久久精品古装| 搞女人的毛片| 亚洲av.av天堂| 日韩欧美精品免费久久| 真实男女啪啪啪动态图| 婷婷色麻豆天堂久久| 欧美另类一区| 亚洲精品中文字幕在线视频 | 成人亚洲精品av一区二区| 免费大片黄手机在线观看| 久久精品久久久久久噜噜老黄| 色播亚洲综合网| 亚洲精品一区蜜桃| 国产精品一区二区性色av| 久久久国产一区二区| 男人爽女人下面视频在线观看| 久久久精品免费免费高清| 精品人妻偷拍中文字幕| 国产精品爽爽va在线观看网站| 国产乱来视频区| 日日摸夜夜添夜夜添av毛片| 日韩精品有码人妻一区| 国产精品熟女久久久久浪| 午夜精品一区二区三区免费看| av在线蜜桃| 国国产精品蜜臀av免费| 干丝袜人妻中文字幕| 搞女人的毛片| 各种免费的搞黄视频| 国产老妇伦熟女老妇高清| 精品视频人人做人人爽| 国产精品无大码| 日日摸夜夜添夜夜爱| 日本三级黄在线观看| 国产有黄有色有爽视频| 久久ye,这里只有精品| 超碰97精品在线观看| 天天一区二区日本电影三级| 亚洲av男天堂| 国产人妻一区二区三区在| 97超碰精品成人国产| 91aial.com中文字幕在线观看| 久久6这里有精品| 尤物成人国产欧美一区二区三区| 久久久久久九九精品二区国产| 久久久精品欧美日韩精品| 一级av片app| 少妇高潮的动态图| 精品午夜福利在线看| 欧美xxxx性猛交bbbb| 高清视频免费观看一区二区| 99久久中文字幕三级久久日本| 久久久色成人| 国产伦在线观看视频一区| 国产一区二区在线观看日韩| 丰满少妇做爰视频| 赤兔流量卡办理| 久久99蜜桃精品久久| 69av精品久久久久久| 极品少妇高潮喷水抽搐| 午夜日本视频在线| 少妇猛男粗大的猛烈进出视频 | 欧美另类一区| 99久久人妻综合| 亚洲欧美精品自产自拍| 一级a做视频免费观看| 久久99热6这里只有精品| 国产精品一区二区三区四区免费观看| 国产精品一及| 免费在线观看成人毛片| .国产精品久久| 国产精品人妻久久久影院| 少妇人妻精品综合一区二区| 久久久久久久久久久免费av| 日本猛色少妇xxxxx猛交久久| 久久精品国产亚洲网站| 国产男女超爽视频在线观看| 国产视频首页在线观看| 久久久久久久久久久丰满| 制服丝袜香蕉在线| 久久影院123| 中文字幕免费在线视频6| 日本免费在线观看一区| 久久久精品欧美日韩精品| 午夜精品国产一区二区电影 | 毛片一级片免费看久久久久| 天堂网av新在线| 永久免费av网站大全| 欧美性感艳星| 如何舔出高潮| 国产爱豆传媒在线观看| 国产黄片视频在线免费观看| 国产色婷婷99| 国产精品成人在线| 丰满少妇做爰视频| a级一级毛片免费在线观看| 亚洲熟女精品中文字幕| 国产淫片久久久久久久久| 日韩在线高清观看一区二区三区| 色吧在线观看| 嘟嘟电影网在线观看| 大陆偷拍与自拍| 大片电影免费在线观看免费| 极品教师在线视频| 肉色欧美久久久久久久蜜桃 | 成年版毛片免费区| 亚洲内射少妇av| 三级国产精品片| 久久女婷五月综合色啪小说 | 真实男女啪啪啪动态图| 五月伊人婷婷丁香| 最近2019中文字幕mv第一页| 国产毛片a区久久久久| 亚洲精品乱久久久久久| 国产一级毛片在线| 中文天堂在线官网| 性插视频无遮挡在线免费观看| 午夜精品一区二区三区免费看| 亚洲自拍偷在线| 能在线免费看毛片的网站| 亚洲人成网站在线观看播放| 久久久久久久久大av| 亚洲欧美一区二区三区国产| 全区人妻精品视频| 久久99热这里只频精品6学生| 又爽又黄a免费视频| 国产乱人视频| 一区二区三区乱码不卡18| 80岁老熟妇乱子伦牲交| 久久久久久久久大av| 韩国av在线不卡| 久久久久久久久久人人人人人人| 18禁在线播放成人免费| 成人一区二区视频在线观看| 丝袜脚勾引网站| 日韩国内少妇激情av| 日日撸夜夜添| 国产一级毛片在线| 亚洲精品日韩在线中文字幕| 日本一二三区视频观看| 亚洲欧美日韩东京热| 久久人人爽人人爽人人片va| 特大巨黑吊av在线直播| 亚洲欧美一区二区三区黑人 | 日本猛色少妇xxxxx猛交久久| 特大巨黑吊av在线直播| 蜜臀久久99精品久久宅男| 一级毛片我不卡| 国产精品国产av在线观看| 久久99热这里只有精品18| 黄色视频在线播放观看不卡| 色网站视频免费| 综合色av麻豆| 国产色婷婷99| 七月丁香在线播放| 国产精品蜜桃在线观看| 亚洲精品日韩在线中文字幕| 亚洲精品日本国产第一区| 不卡视频在线观看欧美| a级毛片免费高清观看在线播放| 99久久人妻综合| 2022亚洲国产成人精品| 老司机影院成人| 黄色一级大片看看| 搡老乐熟女国产| 国产色爽女视频免费观看| 少妇人妻一区二区三区视频| 狂野欧美激情性bbbbbb| 日本黄大片高清| 啦啦啦中文免费视频观看日本| 男人添女人高潮全过程视频| 最近最新中文字幕免费大全7| 久久99热6这里只有精品| 天堂中文最新版在线下载 | 小蜜桃在线观看免费完整版高清| 各种免费的搞黄视频| 极品教师在线视频| 国产精品偷伦视频观看了| 2021少妇久久久久久久久久久| 搞女人的毛片| 久久久精品欧美日韩精品| 国产高清三级在线| 日本色播在线视频| 日韩一区二区三区影片| 国产免费福利视频在线观看| 亚洲aⅴ乱码一区二区在线播放| 波多野结衣巨乳人妻| 色吧在线观看| 久久精品国产鲁丝片午夜精品| 只有这里有精品99| 欧美xxxx性猛交bbbb| 国产精品久久久久久久电影| 国产精品久久久久久av不卡| 国产毛片a区久久久久| 狂野欧美激情性xxxx在线观看| 免费黄频网站在线观看国产| 免费播放大片免费观看视频在线观看| 蜜桃亚洲精品一区二区三区| av在线观看视频网站免费| 久久久色成人| 夜夜看夜夜爽夜夜摸| 啦啦啦中文免费视频观看日本| 波多野结衣巨乳人妻| 各种免费的搞黄视频| 国产爽快片一区二区三区| 蜜臀久久99精品久久宅男| 久久久色成人| 成年女人看的毛片在线观看| 亚洲综合精品二区| av福利片在线观看| 精品少妇久久久久久888优播| 国产在线男女| 国产国拍精品亚洲av在线观看| 尤物成人国产欧美一区二区三区| 亚洲av成人精品一二三区| 成人鲁丝片一二三区免费| 22中文网久久字幕| 99久国产av精品国产电影| h日本视频在线播放| 成人漫画全彩无遮挡| 国产欧美日韩精品一区二区| 国产男人的电影天堂91| 在线观看一区二区三区激情| 国产成人免费观看mmmm| 日本三级黄在线观看| 亚洲天堂av无毛| 麻豆乱淫一区二区| 男女边摸边吃奶| 99久久精品国产国产毛片| 亚洲经典国产精华液单| 国产成人免费无遮挡视频| 亚洲欧美成人精品一区二区| 夫妻午夜视频| 真实男女啪啪啪动态图| 精品久久久久久电影网| 成年av动漫网址| 国产精品福利在线免费观看| 青春草国产在线视频| 国产午夜福利久久久久久| 国产免费一区二区三区四区乱码| 国产白丝娇喘喷水9色精品| 激情 狠狠 欧美| 成人欧美大片| 中文精品一卡2卡3卡4更新| 熟妇人妻不卡中文字幕| 搡女人真爽免费视频火全软件| 婷婷色av中文字幕| 免费看a级黄色片| 国产精品99久久99久久久不卡 | 狂野欧美激情性bbbbbb| av在线蜜桃| 美女cb高潮喷水在线观看| 一个人看的www免费观看视频| 久久综合国产亚洲精品| 日韩成人伦理影院| 深爱激情五月婷婷| 欧美老熟妇乱子伦牲交| 偷拍熟女少妇极品色| 国产精品国产av在线观看| 国产黄a三级三级三级人| 免费黄频网站在线观看国产| 国产精品国产三级国产av玫瑰| 国产久久久一区二区三区| 看非洲黑人一级黄片| 啦啦啦啦在线视频资源| 精品人妻视频免费看| 韩国高清视频一区二区三区| 黄片无遮挡物在线观看| 国产一级毛片在线| 少妇被粗大猛烈的视频| 欧美高清性xxxxhd video| 看非洲黑人一级黄片| 免费看日本二区| av在线观看视频网站免费| 亚洲欧美日韩另类电影网站 | 久久精品久久久久久噜噜老黄| 国产一级毛片在线| 少妇被粗大猛烈的视频| 亚洲不卡免费看| 日韩中字成人| 成人亚洲精品av一区二区| 五月伊人婷婷丁香| 欧美成人午夜免费资源| 日韩不卡一区二区三区视频在线| 久久久久久久久久久丰满| 国产精品av视频在线免费观看| 18禁动态无遮挡网站| 国产精品无大码| 久久精品久久久久久久性| 亚洲精品,欧美精品| 国产乱来视频区| 蜜臀久久99精品久久宅男| 麻豆乱淫一区二区| 精品人妻视频免费看| 男人狂女人下面高潮的视频| 欧美日韩精品成人综合77777| 欧美成人精品欧美一级黄| 亚洲天堂av无毛| 美女视频免费永久观看网站| 黄片wwwwww| 成人亚洲精品av一区二区| 黄片wwwwww| 久久久久久久国产电影| 欧美日韩精品成人综合77777| 欧美日韩亚洲高清精品| 久久久久精品久久久久真实原创| 99久久人妻综合| 久久精品国产亚洲av涩爱| 一级二级三级毛片免费看| 色5月婷婷丁香| 久久影院123| 色5月婷婷丁香| 少妇裸体淫交视频免费看高清| av在线亚洲专区| 国产综合精华液| 不卡视频在线观看欧美| 丝袜美腿在线中文| 在线a可以看的网站| 精品久久久精品久久久| 肉色欧美久久久久久久蜜桃 | av播播在线观看一区| 超碰97精品在线观看| 日本色播在线视频| 别揉我奶头 嗯啊视频| 国产精品av视频在线免费观看| 国产乱人偷精品视频| 七月丁香在线播放| 国产探花在线观看一区二区| 大片电影免费在线观看免费| kizo精华| 国产av国产精品国产| 伦理电影大哥的女人| 午夜免费鲁丝| 99久久九九国产精品国产免费| 大片免费播放器 马上看| 人人妻人人澡人人爽人人夜夜| 好男人视频免费观看在线| 精品国产露脸久久av麻豆| 99久国产av精品国产电影| 国产高潮美女av| 在线免费观看不下载黄p国产| 亚洲一区二区三区欧美精品 | 男女下面进入的视频免费午夜| 免费黄网站久久成人精品| 亚洲美女视频黄频| 国产成人a区在线观看| 国产成人福利小说| 91狼人影院| 久久97久久精品| 夫妻性生交免费视频一级片| 国产精品嫩草影院av在线观看| 两个人的视频大全免费| 久久久久性生活片| 哪个播放器可以免费观看大片| 中文字幕亚洲精品专区| 少妇人妻久久综合中文| 亚洲av成人精品一区久久| 一级片'在线观看视频| 美女被艹到高潮喷水动态| 在线观看国产h片| 精品少妇黑人巨大在线播放| 久久午夜福利片| 高清午夜精品一区二区三区| 日韩 亚洲 欧美在线| 街头女战士在线观看网站| 成人亚洲精品av一区二区| 国产成人91sexporn| 又黄又爽又刺激的免费视频.| 少妇的逼好多水| 97超碰精品成人国产| 91在线精品国自产拍蜜月| 网址你懂的国产日韩在线| 丰满少妇做爰视频| 国产爱豆传媒在线观看| 一级av片app| 九九久久精品国产亚洲av麻豆| 国产亚洲av片在线观看秒播厂| 免费电影在线观看免费观看| 欧美潮喷喷水| 久久久久久久精品精品| 国产精品一区二区在线观看99| 大香蕉97超碰在线| 91午夜精品亚洲一区二区三区| 啦啦啦啦在线视频资源| www.色视频.com| 国产探花在线观看一区二区| 超碰av人人做人人爽久久| 看十八女毛片水多多多| 亚洲精品日本国产第一区| 国产极品天堂在线| 中文字幕制服av| 国产成人91sexporn| 联通29元200g的流量卡| 成人漫画全彩无遮挡| 成年版毛片免费区| 日本wwww免费看| 新久久久久国产一级毛片| 国产男女内射视频| 伦精品一区二区三区| av卡一久久| 国产精品一区二区性色av| a级毛色黄片| 久久精品久久久久久噜噜老黄| 国产高清三级在线| 久久影院123| 亚洲色图av天堂| 禁无遮挡网站| 国产亚洲5aaaaa淫片| 欧美极品一区二区三区四区| 中文精品一卡2卡3卡4更新| 欧美日本视频| 国产高清有码在线观看视频| 免费观看av网站的网址| 97精品久久久久久久久久精品| 五月伊人婷婷丁香| 天堂网av新在线| 成人漫画全彩无遮挡| 一个人看的www免费观看视频| 五月天丁香电影| 蜜臀久久99精品久久宅男| 中文字幕亚洲精品专区| 五月玫瑰六月丁香| 国产成人一区二区在线| 最后的刺客免费高清国语| 久久久久久九九精品二区国产| 国产精品伦人一区二区| 国产精品国产三级国产av玫瑰| 免费av观看视频| 男男h啪啪无遮挡| 国产高潮美女av| 嘟嘟电影网在线观看| 中国三级夫妇交换| 免费av观看视频| av卡一久久| 亚洲av日韩在线播放| 青春草国产在线视频| 菩萨蛮人人尽说江南好唐韦庄| 大又大粗又爽又黄少妇毛片口| 99九九线精品视频在线观看视频| 九九在线视频观看精品| 日韩一区二区三区影片| 99热网站在线观看| 国产在线男女| 毛片一级片免费看久久久久| 日韩人妻高清精品专区| 午夜精品国产一区二区电影 | 亚洲美女搞黄在线观看| 2018国产大陆天天弄谢| 赤兔流量卡办理| 久久ye,这里只有精品| 亚洲欧美成人精品一区二区| 天天躁日日操中文字幕| 欧美xxⅹ黑人| 午夜爱爱视频在线播放| 国产极品天堂在线| 在线亚洲精品国产二区图片欧美 | 精品少妇黑人巨大在线播放| 2021少妇久久久久久久久久久| 国产精品成人在线| 久久久a久久爽久久v久久| 成人漫画全彩无遮挡| 插逼视频在线观看| 国产视频内射| 爱豆传媒免费全集在线观看| 在线精品无人区一区二区三 | 欧美一级a爱片免费观看看| 欧美性感艳星| 亚洲第一区二区三区不卡| 久久99热6这里只有精品| 免费电影在线观看免费观看| 伦精品一区二区三区| 国产亚洲精品久久久com| 亚洲欧美日韩无卡精品| 国产色婷婷99| 91在线精品国自产拍蜜月| 我的老师免费观看完整版| 欧美日韩亚洲高清精品| 有码 亚洲区| 国产黄片美女视频| 狂野欧美激情性bbbbbb| 亚洲国产精品专区欧美| 国产精品国产av在线观看| av国产久精品久网站免费入址| 亚洲激情五月婷婷啪啪| 又爽又黄a免费视频| 一级二级三级毛片免费看| 肉色欧美久久久久久久蜜桃 | 各种免费的搞黄视频| 久久精品国产亚洲av涩爱| 五月天丁香电影| 久久精品国产自在天天线| 秋霞伦理黄片| 王馨瑶露胸无遮挡在线观看| 亚洲真实伦在线观看| 观看免费一级毛片| 久久99热这里只频精品6学生| 国产高潮美女av| 亚洲不卡免费看| 久久久精品欧美日韩精品| 肉色欧美久久久久久久蜜桃 | 亚洲自偷自拍三级| 国产欧美亚洲国产| 欧美xxxx性猛交bbbb| 亚洲av.av天堂| 99re6热这里在线精品视频| 成年人午夜在线观看视频| 3wmmmm亚洲av在线观看| 久久久精品免费免费高清| 成人无遮挡网站| 只有这里有精品99| 午夜精品一区二区三区免费看| 日日摸夜夜添夜夜爱| 亚洲一区二区三区欧美精品 | 中国国产av一级| 国产午夜福利久久久久久| 国产欧美日韩一区二区三区在线 | 三级经典国产精品| 麻豆乱淫一区二区| 精品午夜福利在线看| 欧美成人精品欧美一级黄| 国产成人福利小说| 国产男女超爽视频在线观看| 国产69精品久久久久777片| 又粗又硬又长又爽又黄的视频| 黄色配什么色好看| 国产黄色视频一区二区在线观看| 国产精品一及| 成人高潮视频无遮挡免费网站| 国产 精品1| 国产精品一及| 人妻制服诱惑在线中文字幕| 国产中年淑女户外野战色| 国产精品一及| 国产欧美亚洲国产| 99热6这里只有精品| 我的女老师完整版在线观看| 久久精品国产亚洲av涩爱| 亚洲成人久久爱视频| 免费电影在线观看免费观看| 免费大片黄手机在线观看| 少妇人妻精品综合一区二区| 爱豆传媒免费全集在线观看| 国产精品.久久久| 韩国高清视频一区二区三区| 国产成人福利小说|