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

    Numerical simulation of the effects of protrusion on DC arc anode attachment

    2021-10-31 08:14:50ChongNIU牛沖XianMENG孟顯HejiHUANG黃河激TaoZHU朱濤SurongSUN孫素蓉andHaixingWANG王海興
    Plasma Science and Technology 2021年10期
    關(guān)鍵詞:黃河

    Chong NIU (牛沖), Xian MENG (孟顯), Heji HUANG (黃河激),Tao ZHU (朱濤), Surong SUN (孫素蓉),3,?and Haixing WANG (王海興),3,?

    1 School of Astronautics, Beihang University, Beijing 100191, People’s Republic of China

    2 Institute of Mechanics, Chinese Academy of Sciences, Beijing 100190, People’s Republic of China

    3 Ningbo Institute of Technology, Beihang University, Ningbo 315800, People’s Republic of China

    Abstract The attachment of the DC arc on the anode is usually affected by surface morphology such as protrusions due to ablation or melting deformation.A three-dimensional thermodynamic and chemical non-equilibrium model is used to numerically simulate the effect of artificially assumed surface protrusions on the arc anode attachment.The numerical simulation results show that the arc deflects toward the protrusions on the anode and attaches to them in a constricted mode,resulting in an increase in the temperature of the arc attachment region.The analysis shows that the presence of protrusion on the anode surface changes the electric field distribution,intensifies the degree of thermodynamic and chemical non-equilibrium in its vicinity,further influences the chemical kinetic process of the plasma around it, which is the main reason for the deflection of the arc toward the protrusions and the arc anode attachment in a constricted mode.In order to verify the numerical simulation results,verification experiments are also performed using similar size scale anode protrusion, and the results showed that the presence of protrusion can indeed cause the deflection of the arc and even cause the ablation of the protrusion.

    Keywords: DC arc, arc anode attachment, anode protrusion, chemical non-equilibrium

    1.Introduction

    The anodes of the DC arc devices used in the fields of electric propulsion and material processing are often subjected to high current and heat flow fluxes.In the case of extremely high heat flux or insufficient cooling, the anode suffers noticeable damage at the attachment region, especially for the constricted arc attachment mode.The temperature of the arc anode attachment region continues to rise, which may cause deformation, melting, and even ablation of the anode surface to form some structures with different sizes, such as protrusions or craters.These protrusions or craters have important effects on the arc anode attachment,causing the arc parameter distribution lose symmetry and even stability.Although there have been intensive investigations on the anode phenomena of DC arc which are related to the anode-plasma interaction,such as the physics of anode spots [1–5] and arc attachment mode [6–9], most studies are based on the assumption of a smooth anode surface.Our understanding on the effects of anode morphology on the arc attachment mode is still inadequate, mainly due to the complicated interaction mechanisms of the arc plasma and anode involved.

    Previous studies have shown that surface morphology of the electrodes, such as bumps, burrs, roughness, etc, has substantial influence on the mode of electrode arc attachment,while most of the studies focused to the interaction between the cathode and the arc [10–12].In the field of vacuum arc research, the effect of the protrusion on the cathode surface has attracted wide attention, which is because the protrusion plays an important role in the arc spot formation and evolution [13, 14].Actually, the effect of protrusions on the DC arc attachment on the anode surface is also important which is closely related to the ablation of anode.Since the protrusions are inside the anode boundary layer which is characterized by large gradients of the plasma parameters,so it is very difficult to analyze and study their effect on arc attachment.Therefore,the development of numerical simulation approach is important to investigate this process.As demonstrated in previous studies [15–19],at the arc fringes and inside the arc anode boundary layer, the plasma deviates significantly from the thermodynamic and chemical equilibrium [20–30].Therefore thermodynamic and chemical non-equilibrium should be important features of the region near the anode protrusions and should be considered in the physical model.A reasonable chemical kinetic model should be used to simulate the effect of protrusions on arc anode attachment[31–35].

    The main purpose of this paper is to provide a preliminary understanding of the influence of protrusions on arc anode attachment.An artificially assumed millimeter-level protrusion on the anode surface is adopted to study the influence of the protrusion on the arc attachment.In the study,the anode solid zone is also coupled to computational domain,allowing to obtain the temperature distribution of anode protrusion.Numerical simulation results of the effect of protrusions on arc anode attachment are presented, followed by preliminary experimental validation.

    2.Modeling approach

    2.1.Computation domain

    Figure 1 shows the schematic diagram of the wall-stabilized transferred arc device used in the calculation.The length and radius of the water-cooled constrictor tube are 40 mm and 5 mm, respectively, which are the same sizes with the experimental setup in[36].The distance between the outlet of the water-cooled constrictor tube and the anode surface is 12 mm,the radius of the anode is 15 mm,and the thickness of the anode is 10 mm.The cathode is assumed to be a round rod with a radius of 2 mm, with a top cone angle of 60°, and a small platform with a radius of 0.5 mm.In this study, we artificially set a trapezoidal protrusion on the anode surface.The radii of the bottom and the top of the protrusion are 0.5 mm and 0.2 mm,and the height is 1 mm.Considering that the protrusion of the anode surface may cause the arc to lose axi-symmetry,three-dimensional simulation is adopted in this study.Since the arc and the device are symmetric about theyz-plane (x=0 mm), ABCDEFGHIJ in the figure 1 is set as computational domain, which is the half of the arc device in this simulation.In this simulation,the governing equations of arc plasma and anode region are solved together.

    Figure 1.Schematic diagram of wall-stabilized transferred arc setup and computational domain used in calculation.

    Figure 2.Comparison between the calculated electron temperature (lines) and the experimental data (symbols).

    Figure 3.Heavy-species temperature distributions for the cases without(a)and with protrusion at different radial positions.(b)y=1 mm,(c)y=2 mm, (d) y=3 mm.

    Figure 4.Temperature distributions for the cases without and with protrusion at the radial position y=2 mm.

    Figure 5.Electron number density distributions for the cases without(a)and with protrusion at different radial positions.(b)y=1 mm,(c)y=2 mm, (d) y=3 mm.

    2.2.Main assumptions and governing equations

    The following assumptions are used in this simulation.(1)The flow is considered as steady and laminar.(2) The arc anode is operated without ablation, even if it is considered that the artificially set protrusion may be heated to a higher temperature.(3) Arc plasma is electrically neutral, and the radiation loss is evaluated by net emission coefficient with absorption radiusR=0 mm.(4) The thermodynamic nonequilibrium chemical non-equilibrium are considered, and a five-species (ground state atoms, excited state atoms, atomic ions,molecular ions and electrons)is include in the chemical kinetic model.

    Based on the above assumptions, the required governing equations are as follows.

    Gas species equation

    Here,ρsandare the mass density and diffusion flux of speciess.In this simulation,the diffusion fluxis calculated by the self-consistent effective binary diffusion approximation developed by Ramshaw [37–39].Sc,sis the chemical reaction source term which is determined by different kinetic processes shown in table 1.The diffusion flux→are obtained from equation (2) [40, 41].Here,pis the gas pressure.Msis the molecular weight andDsis effective diffusion coefficient.This the heavy-species temperature solved by heavy-species energy equation.Rgis the universal gas constant andysis the mass fraction of speciess.is the volumetric driving force defined asHere,qsis the charge per unit mass of speciessandis the ambipolar diffusion electric field.→is the driving force shown in equation (3)whereis the partial pressure of speciess,andβsjis the thermal diffusion coefficient between speciessandj

    Table 1.The reactions taken into account in the numerical model.

    Momentum conservation equation

    Current continuity and Maxwell’s equation

    whereμ0is the the permeability.

    The electrical conductivity is given by equation (7).Here,ecis the elementary charge.neandmeare number density and mass of electron respectively.vmis the sum of collision frequencies between electrons and heavy-species including the atoms and ions

    Electron energy equation

    whereTeandkBare respectively the electron temperature and Boltzmann constant.is the electron energy flux density,which can be expressed as follows

    wherekeis the electron thermal conductivity,heis the specific enthalpy of electron and→is the diffusion flux of electron which can be calculated by relation[35].In the right-hand of equation (8),represents the Joule heating andis the electric field which is calculated by relationφ= ?,here,φis the electric potential.Q,elQin,eare respectively the electron elastic collision term and the reaction heat due to electron impact reactions,which can be calculated by

    whereqe‐srepresents the cross sections for collisions between electrons and heavy-species, which are obtained from [49].Δε iis the energy loss due to inelastic collisions which is related to the corresponding kinetic process.vi,srepresents the stoichiometric coefficients corresponding to speciessin reactioni.In equation(8),Qradis net radiation losses obtained from [50].

    Table 2.The boundary conditions of the numerical model.

    Heavy-species energy equation

    Here,nhandkhare respectively the number density and the thermal conductivity of heavy-species.hsand represent the species specific enthalpy.Qin,his the reaction heat due to heavy-species impact reactions.

    In this paper, five-species are considered in the model.The excited state atoms (Ar*), molecular ions (Ar2+), and atomic ions (Ar+) are obtained from species equation.The electron density is determined by the charge neutrality condition,The density of ground state atom(Ar) is calculated by state equation.

    Since in this simulation, the copper is chosen as anode material.The anode is coupled to plasma region and its energy equation is as follows

    Here,Tsis the temperature of the solid anode.ρa,Cp,a,kaandσaare the density,specific heat,thermal conductivity and electrical conductivity in copper anode, the data are obtained from [51].

    Thermodynamic and transport properties involved in the previous governing equations are important for obtaining reasonable arc parameters and distributions, especially in the region with large parameter gradients.In this paper, these properties are calculated based on the temperature and species number density of the grid nodes during each iteration.The calculation methods of thermodynamic and transport properties of plasma under non-equilibrium conditions can be found in [52–54].The transport properties are calculated based on the Chapman–Enskog method,in which the collision between the particles is evaluated by collision integral.The collision integral can be obtained by integrating the interaction potential between particles.The collision integrals between neutral particles are calculated by HFDTCS2 potential [55],and the interaction potential of ions and neutral particles is obtained from [56].The Coulomb potential is used for the calculation of collision integral between charged particles,and the collision integrals of electrons and neutral particles are obtained from [57].

    2.3.Boundary conditions

    The boundary conditions used in the calculation are shown in table 2.The gas enters from the annular channel between the inner wall of the constrictor and the outer wall of the cathode with flow rate of 0.35 slpm.The gas flows out from the gap between the bottom wall of the constrictor and the anode surface and the pressure is set as 1 atm at outlet.Non-slip flow conditions are set at the walls of water-cooled constrictor and the anode surface.

    For heavy-species temperature,the gas temperature is set as 300 K at the inlet, outlet and inner and bottom walls of constrictor.Along the cathode surface, the temperature is assumed to gradually change from 300 K at the rod part to 3500 K at the cathode flat tip.The anode bottom is watercooled,and its temperature is set to be 300 K,while the upper surface of anode is coupled to the arc plasma and its temperature can be obtained self-consistently.For the electron temperature, it is set to be 300 K at the inlet.At other boundaries, zero electron temperature gradient is adopted.

    For anode,several important heat transfer processes,such as electron condensation heat, electron enthalpy transport,conduction heat and radiation loss are taken into account[23,51,58,59].The conduction heat transfer between the arc plasma and anode is solved by coupling.Therefore, the boundary conditions of anode heat transfer in this model can be set as follow:

    Here, the first term is the electron enthalpy transport and the electron condensation heat, and the second term is the anode radiation loss.Aa,εCuandβare respectively the work function, surface emissivity and Stefan–Boltzmann constant,which are taken from [51].

    At the cathode surface, the current density is assumed to be a Gaussian distribution as in [60], the peak value ofjmax=3.0×107A m?2is chosen forI=50 A.At the anode,the electric potential on the bottom surface of anode is connected to ground.Electric insulation is adopted along the other boundaries.For the magnetic vector potential equation, the Dirichlet boundaries are set at inlet, outlet and the walls of constrictor tube, and zero gradients are adopted at other boundaries.For species equation, zero gradients are set at all boundaries.

    It is expected that the presence of protrusion would lead to the ablation of the copper at the tip of protrusion,and then the copper vapor would enter the plasma zone.Therefore it is necessary to properly consider the effects of copper vapor on the plasma properties near the protrusion, especially the electrical conductivity.A more reasonable method is to develop a reasonable chemical non-equilibrium kinetic model that includes the mixing of metallic copper vapor and its ionization process with argon plasma.However,this is a very challenging task.In this simulation, we use a method similar to previous literatures[61,62]to improve the conductivity of the boundary layer.In this simulation, an artificially high electrical conductivity of a mixture of 50%Ar and 50%Cu is used to increase the conductivity of the boundary layer within the thickness of 0.2 mm around the protrusion.The electrical conductivity of the Ar–Cu mixture is calculated based on the Chapman–Enskog method, and the collision integrals required for the calculation are obtained from the [63].

    A non-uniform mesh with 30 (x-direction) × 56 (y-direction) × 94 (z-direction) is used in the simulation, and a refined grid is adopted near the electrodes, especially around the protrusion region of the anode surface.The parameter distributions are the same compared with those calculated by finer mesh so that the numerical results do not depend on grid accuracy.For decrease of computational effort,therefore,the mesh with 30(xdirection) × 56 (y-direction) × 94 (z-direction) is used.The complete set of coupled governing equations listed above is solved based on the COMSOL software platform [64].

    3.Results and discussion

    3.1.Comparison and verification of numerical simulation results and experimental measurement results

    To verify the reliability of code in this simulation, we compared the calculated results with the experimental results reported in [36] based on the same conditions.The anode region of arc with a superimposed cold flow was investigated in[36]for the wall-stabilized arc device shown in figure 1.In[36], three-dimensional electron temperature distributions are obtained with a laser Thomson scattering system,for the case with arc current of 100 A and the argon flow rate of 5 slpm,and cross flow rate of 12 slpm.

    From the comparison of the calculated electron temperature and the experimental measurement results given in figure 2,it can be seen that the calculated temperatures are in good agreement with the experimental results.However, in the upstream region of the lateral gas flow,which is the region with a larger parameter gradient,the numerical results deviate from the experimental data.The temperature gradients of experimental data are steep, while the gradients at the same position are underestimated by the numerical model.The predicted results of the same numerical approach are also compared with the experimental results presented in [65], as shown in [40].Generally, the results of numerical simulation predictions are acceptable in this simulation.

    3.2.Effects of the existence of protrusion on arc temperature and electron number density distributions

    In this section, the effects of protrusion on temperature and electron number density are discussed.The arc device is operated at current of 50 A and argon flow rate of 0.35 slpm.It can be seen from figure 3 that in the case without protrusion,the temperature distribution of the arc column is axi-symmetric.However,in the case with the protrusion,the arc anode attachment shifts to the position of the protrusion.For the cases where the protrusion is located 1 mm and 2 mm away from the device axis, the protrusion makes the arc attachment more constricted,and a high temperature zone with a temperature of 10 000 K is formed in the arc column above the protrusion.In the case where the protrusion is located 3 mm away from the device axis, the constriction and attraction effects of the protrusion on the arc is weakened so that the maximum temperature around the protrusion decreases.The protrusions of the anode are intensively heated by arc plasma.The highest temperature at the tip of the protrusions in three cases has already exceeded the melting point of copper, 1357 K, which increases the likelihood of the protrusion being ablated.

    Figure 4 presents the distributions of electron temperature and heavy-species temperature along they-direction 1.1 mm above the anode,i.e.0.1 mm above the protrusion.In the case without protrusion, the deviation of heavy-species temperature from electron temperature is small at the device axis, while the degree of thermodynamic non-equilibrium increases with the increase of distance from device axis.For the case with protrusion located aty=2 mm, the electron temperature around protrusion increases to 12 000 K due to the arc constriction, while the heavy-species temperature decreases to 5500 K.This indicates that the presence of anode protrusion strengthens the degree of thermodynamic nonequilibrium in this region.

    Figure 6.Distributions of (a) electric field strength, (b) current density and (c) electrical conductivity for the cases without and with protrusion at different radial positions of y=1 mm, y=2 mm, and y=3 mm.

    Figure 7.Ionization non-equilibrium degree distributions for the cases without and with protrusion at the radial position y=2 mm.

    The comparisons of electron number density for different cases are shown in figure 5.The presence of protrusion leads to the electron number density near the protrusion being comparable to that in the central area of the outlet of the water-cooled constrictor tube.As the distance between the protrusion and the device axis decreases, the influence of the protrusion on the electron number density distribution around it increases.

    3.3.Analysis of the effects of protrusion on the plasma nonequilibrium process near the anode

    As presented above, the presence of protrusion causes the changes of parameters distribution,further leading to changes in the non-equilibrium plasma processes.Therefore, it is necessary to further analyze the influence of the presence of protrusion on the electric field strength, current density and electrical conductivity.The electric field, current density and electrical conductivity distributions along they-direction 1.1 mm above the anode, i.e.0.1 mm above the protrusion,are given in figure 6.It can be seen that the presence of protrusion greatly increases the electric field distribution.The current density and electrical conductivity also increase.Moreover, the closer the position of the protrusion is to the device axis, the more significant the increase in electric field strength, current density and electrical conductivity.

    The chemical non-equilibrium process is essential to accurately and reasonably describe the characteristics of the plasma near the arc fringe and the electrode regions.The ionization non-equilibrium degree is calculated by the following equation [66]

    Here,kiionandkirecrepresent the rate coefficient of ionization reaction and recombination reaction respectively.The ionization non-equilibrium degreeθis between ?1 and 1.A value ofθequal to 1 indicates that only ionization reactions occur in plasma, whereas a value ofθequal to ?1 indicates that only recombination reactions occur in plasma.Moreover, a value ofθequal to zero means that ionization-recombination equilibrium.

    Figures 7 and 8 show the distributions of ionization nonequilibrium degree and ionization and recombination processes along they-direction 1.1 mm above the anode, i.e.0.1 mm above the protrusion respectively.In the case without protrusion,the degree is approximately 0.1 at the device axis,which indicates that the ionization and recombination reactions are close to equilibrium.However, the ionization non-equilibrium degree approaches ?1 as the distance from device axis increases.This means that the ions and electrons are recombined to atoms in this region.

    Figure 8.Ionization and recombination processes for the cases (a) without and (b) with protrusion at the radial position y=2 mm.

    Figure 9.The main heating processes along the anode surface for the cases with and without protrusion at the radial position y=2 mm.

    Figure 10.Experimental photo of a free burning argon arc with protrusion on anode.

    It can be seen from the figure 8 that the main ionization and recombination mechanisms are respectively the ionization of the excited state of argon atom by electron impact and the three-body recombination of electrons and ions at the device axis.The chemical reaction rate of three-body recombination is slightly lower than that of argon atom excited state ionization.However, recombination reactions of molecular ions and electrons become the dominant reactions at the arc fringe.For the case with protrusion, ionization non-equilibrium degree increases to 0.5 at the protrusion, which is consistent with the increase of electron number density shown in figure 5.It can be seen that the increase of ionization nonequilibrium degree is caused by the considerable increase of rate of excited state atom ionization reaction.Moreover, the rates of recombination reactions have also increased.In addition to the three-body recombination reaction of electrons and atom ions,the recombination reactions of molecular ions and electrons also become important.

    As shown in figures 3(b)–(d), the protrusion region on the anode surface is intensely heated by the plasma, making the tip temperature of protrusion much higher than other areas on the anode surface, which may further lead to anode material melting and ablation.Therefore,it is necessary to further analyze the main heating process that affects the temperature distribution of the protrusions.Figure 9 shows the distribution of heat flux on the anode surface with or without protrusion.In the attachment region, the electron condensation heat is dominant, and its value is higher than the electron enthalpy transport and heat conduction to the anode.For the case with protrusion,the peak heat flux density is almost twice that of the case without protrusion.The presence of protrusion penetrates into the high temperature plasma environment like pin ribs.The water cooling on the lower side of the anode is not enough to reduce the temperature of the top of protrusion,so that the heat at the top of the protrusion accumulates, and the temperature even exceeds the melting point of the metal material.

    Figure 11.Spectral lines of copper vapor in free burning argon arc.

    3.4.Preliminary experimental verification and discussion on the influence of protrusion on arc anode attachment

    As mentioned above in this study, an artificially assumed millimeter-level protrusion on the anode surface is adopted to study the influence of the protrusion on the arc attachment.In order to verify the results of the numerical simulation,we also set up a protrusion of which the size is the same with that in numerical simulation, and conducted experimental observations on a free burning arc device as shown in figure 10.In the experiment, the free burning arc is operated at 50 A current and 15 mm electrodes gap.The presence of protrusion on the anode surface does cause the arc to deflect to the protrusion.It is worth noting that the high temperature of the surface of the protrusion extending into the arc leads to ablation, and the green metallic copper vapor can be observed in the arc fringe.The anode ablation and the presence of copper vapor can be verified by the spectral lines shown in figure 11.

    Preliminary experiments verify the rationality of the numerical simulation in this paper, that is, the presence of protrusion on the anode surface will cause the arc deflection and the change of the arc attachment state.However,it should be pointed out that the current numerical simulation has certain limitations.The main problem is that the process of the metal evaporation in the protrusion area is not considered in a consistent manner.The ablation of protruding metallic copper will lead to two effects on the arc.First, it would change the chemical kinetic process,that is,the chemical nonequilibrium model used in the numerical simulation should take the chemical kinetic process related to copper into account, and second, it should further consider the impact of copper vapor on the plasma transport process and properties.Even with these shortcomings, considering the difficulty and complexity of the above-mentioned problems, our current work is still helpful to improve the understanding of the influence of protrusion on the arc anode attachment.

    4.Conclusion

    In this simulation, artificially assumed protrusion on the anode surface is used to study the influence of protrusion on arc attachment.Thermodynamic and chemical nonequilibrium models are used considering the importance of the non-equilibrium processes in the boundary layer and the arc fringe.Numerical simulation results show that the presence of protrusion causes the arc column deflect to the protrusion and increase of temperatures near the protrusion because of arc constriction.Compared to the case without protrusion,the presence of protrusion strengthens the degree of thermodynamic non-equilibrium around it.Moreover,protrusion extends into the high temperature plasma zone like the pin ribs, so that the protrusion is intensively heated, making the temperature much higher than other anode areas.Therefore, the anode ablation is more likely to occur at the protrusion.

    The presence of protrusion significantly increases the electric field and current density around it.The closer the position of the protrusion is to the device axis,the more significant the increase in electric field strength and current density.The chemical reaction rates also change significantly near the protrusion.For the case with protrusion, the ionization non-equilibrium degree increases resulting in the production of electrons and ions around the protrusion.Moreover, both the electron–ion three-body recombination reaction and molecular ions recombination reactions become dominant recombination processes.

    Preliminary experiments have qualitatively verified the rationality of the numerical simulation in this paper.In the experiment, the arc deflection caused by the presence of protrusion and the ablation caused by the high temperature in the protrusion tip were observed.Although the current nonequilibrium model used in numerical simulation has some limitations, it can still provide a reference for in-depth understanding of related processes.

    Acknowledgments

    This work was supported by National Natural Science Foundation of China (Nos.11735004 and 12005010).

    猜你喜歡
    黃河
    黃河娃
    多彩黃河
    金橋(2020年11期)2020-12-14 07:52:46
    黃河寧,天下平
    金橋(2020年11期)2020-12-14 07:52:42
    黃河黃河,我愛你
    黃河之聲(2020年15期)2020-10-16 01:03:44
    『黃河』
    當代陜西(2019年21期)2019-12-09 08:36:12
    黃河知道我愛誰
    當代音樂(2019年2期)2019-06-11 21:17:05
    去看最大黃河象
    黃河遐想
    黃河放歌
    青年歌聲(2017年5期)2017-03-15 01:21:40
    我家住在黃河邊
    民族音樂(2016年6期)2016-08-28 20:07:16
    一进一出好大好爽视频| 亚洲视频免费观看视频| www.www免费av| 欧美黑人精品巨大| 免费无遮挡裸体视频| 高清在线国产一区| 两性夫妻黄色片| 国产精品久久久久久精品电影 | 黑丝袜美女国产一区| 久久久久久人人人人人| 50天的宝宝边吃奶边哭怎么回事| cao死你这个sao货| 国产一区二区三区综合在线观看| 少妇粗大呻吟视频| 亚洲欧美日韩高清在线视频| 两个人视频免费观看高清| 国产精品香港三级国产av潘金莲| 成人精品一区二区免费| 性少妇av在线| 色播亚洲综合网| 亚洲激情在线av| 99久久99久久久精品蜜桃| 欧美日韩亚洲国产一区二区在线观看| 一区二区三区精品91| 啦啦啦 在线观看视频| 久久久久亚洲av毛片大全| 久久午夜综合久久蜜桃| 午夜福利在线观看吧| 黄色女人牲交| 久久久久久人人人人人| 天天一区二区日本电影三级 | 两性午夜刺激爽爽歪歪视频在线观看 | 涩涩av久久男人的天堂| 一级a爱视频在线免费观看| 久久午夜综合久久蜜桃| 97超级碰碰碰精品色视频在线观看| 99国产综合亚洲精品| 无遮挡黄片免费观看| 女同久久另类99精品国产91| 一级毛片女人18水好多| 精品国产乱子伦一区二区三区| 男男h啪啪无遮挡| or卡值多少钱| 美女午夜性视频免费| 日韩精品免费视频一区二区三区| 成人三级黄色视频| 精品久久久精品久久久| 日本三级黄在线观看| 亚洲中文日韩欧美视频| 99国产精品99久久久久| 久久精品国产亚洲av高清一级| 国产免费av片在线观看野外av| 波多野结衣高清无吗| 少妇熟女aⅴ在线视频| 看黄色毛片网站| 最新在线观看一区二区三区| 午夜精品久久久久久毛片777| 欧美激情高清一区二区三区| 最好的美女福利视频网| 久久精品91蜜桃| 亚洲在线自拍视频| 性少妇av在线| 亚洲成av人片免费观看| 最近最新中文字幕大全免费视频| videosex国产| 国产精品一区二区三区四区久久 | www.精华液| 老司机午夜十八禁免费视频| 高清黄色对白视频在线免费看| 91精品三级在线观看| 日本 av在线| 动漫黄色视频在线观看| 国产精品电影一区二区三区| 欧美日韩亚洲国产一区二区在线观看| 69精品国产乱码久久久| 日本免费a在线| 啦啦啦免费观看视频1| 久久人人精品亚洲av| 一级a爱视频在线免费观看| 免费在线观看视频国产中文字幕亚洲| 久久久久久久久免费视频了| 国产精品永久免费网站| 法律面前人人平等表现在哪些方面| 午夜成年电影在线免费观看| 国产亚洲欧美98| 丝袜人妻中文字幕| 国产真人三级小视频在线观看| 亚洲精品av麻豆狂野| 大香蕉久久成人网| 国产精品免费一区二区三区在线| 精品福利观看| 国产精品永久免费网站| 51午夜福利影视在线观看| 午夜视频精品福利| 国产日韩一区二区三区精品不卡| 精品免费久久久久久久清纯| 日本免费一区二区三区高清不卡 | 精品人妻1区二区| 久久国产亚洲av麻豆专区| 可以在线观看的亚洲视频| 亚洲国产精品久久男人天堂| 午夜福利视频1000在线观看 | 制服诱惑二区| 久久精品91无色码中文字幕| 在线观看免费日韩欧美大片| 日本一区二区免费在线视频| 18禁黄网站禁片午夜丰满| 精品熟女少妇八av免费久了| 午夜福利视频1000在线观看 | 91九色精品人成在线观看| 男男h啪啪无遮挡| 国产99久久九九免费精品| xxx96com| √禁漫天堂资源中文www| 欧美成人免费av一区二区三区| 如日韩欧美国产精品一区二区三区| 久久人妻熟女aⅴ| 香蕉国产在线看| 自拍欧美九色日韩亚洲蝌蚪91| 正在播放国产对白刺激| 18禁黄网站禁片午夜丰满| 在线观看免费午夜福利视频| 别揉我奶头~嗯~啊~动态视频| 国产亚洲av嫩草精品影院| 亚洲aⅴ乱码一区二区在线播放 | 欧美精品亚洲一区二区| 亚洲午夜理论影院| 看片在线看免费视频| 亚洲 欧美一区二区三区| 人妻丰满熟妇av一区二区三区| 国产激情久久老熟女| 国产极品粉嫩免费观看在线| 高清毛片免费观看视频网站| 亚洲精华国产精华精| 午夜福利在线观看吧| 色播亚洲综合网| 久久中文字幕一级| 一a级毛片在线观看| 男女下面进入的视频免费午夜 | 很黄的视频免费| 免费人成视频x8x8入口观看| 精品国产乱码久久久久久男人| 国内毛片毛片毛片毛片毛片| 制服诱惑二区| 久久草成人影院| 精品久久久久久久人妻蜜臀av | 精品久久久久久久人妻蜜臀av | 亚洲专区字幕在线| 亚洲在线自拍视频| www.精华液| 亚洲欧美激情综合另类| 大香蕉久久成人网| 视频在线观看一区二区三区| 国产av精品麻豆| 19禁男女啪啪无遮挡网站| 亚洲,欧美精品.| 亚洲第一av免费看| 午夜免费成人在线视频| 在线观看免费视频日本深夜| 天天一区二区日本电影三级 | 丰满的人妻完整版| 日韩精品免费视频一区二区三区| 乱人伦中国视频| 亚洲全国av大片| 亚洲一区中文字幕在线| 91九色精品人成在线观看| 欧美最黄视频在线播放免费| 法律面前人人平等表现在哪些方面| 97超级碰碰碰精品色视频在线观看| 国产亚洲精品久久久久5区| 天堂√8在线中文| 亚洲久久久国产精品| 欧美黑人精品巨大| 男女午夜视频在线观看| 亚洲国产精品sss在线观看| 自线自在国产av| 久久久久国产精品人妻aⅴ院| 一级毛片精品| 日韩大尺度精品在线看网址 | 亚洲成人免费电影在线观看| 欧美激情久久久久久爽电影 | 国产高清激情床上av| 亚洲国产精品999在线| 此物有八面人人有两片| 精品国产亚洲在线| 亚洲电影在线观看av| 夜夜看夜夜爽夜夜摸| 黄网站色视频无遮挡免费观看| 又紧又爽又黄一区二区| 亚洲最大成人中文| 国产精品野战在线观看| 欧美久久黑人一区二区| 黑人欧美特级aaaaaa片| 午夜成年电影在线免费观看| 成人三级做爰电影| av片东京热男人的天堂| 桃红色精品国产亚洲av| 亚洲五月天丁香| 丝袜美腿诱惑在线| 在线天堂中文资源库| АⅤ资源中文在线天堂| 国产精品一区二区三区四区久久 | 亚洲欧美激情综合另类| 国产高清有码在线观看视频 | 一级,二级,三级黄色视频| 丝袜美腿诱惑在线| 欧美日韩亚洲综合一区二区三区_| 后天国语完整版免费观看| 免费在线观看日本一区| 女人被躁到高潮嗷嗷叫费观| 国产私拍福利视频在线观看| 成年人黄色毛片网站| 中文字幕av电影在线播放| 天天躁夜夜躁狠狠躁躁| av在线天堂中文字幕| 无人区码免费观看不卡| 亚洲成a人片在线一区二区| 国产精品亚洲美女久久久| 成人国产综合亚洲| 两个人免费观看高清视频| 成人亚洲精品一区在线观看| 一级,二级,三级黄色视频| 最新美女视频免费是黄的| 亚洲电影在线观看av| 黄色女人牲交| 国产av精品麻豆| 香蕉国产在线看| 91av网站免费观看| 后天国语完整版免费观看| 成年女人毛片免费观看观看9| 欧洲精品卡2卡3卡4卡5卡区| 国产成人免费无遮挡视频| 亚洲自偷自拍图片 自拍| 少妇 在线观看| 波多野结衣高清无吗| 1024香蕉在线观看| 丝袜美腿诱惑在线| 久久九九热精品免费| 中出人妻视频一区二区| av天堂久久9| 亚洲午夜理论影院| 中文亚洲av片在线观看爽| 人人妻,人人澡人人爽秒播| 少妇被粗大的猛进出69影院| 成在线人永久免费视频| 99riav亚洲国产免费| 国内精品久久久久久久电影| 日韩欧美国产在线观看| 老司机深夜福利视频在线观看| 在线观看www视频免费| 欧美在线黄色| 琪琪午夜伦伦电影理论片6080| 欧美日韩乱码在线| 亚洲一卡2卡3卡4卡5卡精品中文| 精品日产1卡2卡| 精品国产国语对白av| 夜夜看夜夜爽夜夜摸| 日韩精品免费视频一区二区三区| 日韩三级视频一区二区三区| a级毛片在线看网站| 夜夜夜夜夜久久久久| 免费观看精品视频网站| 岛国在线观看网站| 亚洲 国产 在线| 免费av毛片视频| 制服丝袜大香蕉在线| 欧美丝袜亚洲另类 | 久久狼人影院| 午夜亚洲福利在线播放| 亚洲第一欧美日韩一区二区三区| 9色porny在线观看| 亚洲色图综合在线观看| 熟妇人妻久久中文字幕3abv| 又紧又爽又黄一区二区| 国产精品免费视频内射| 在线观看www视频免费| 国产欧美日韩一区二区三区在线| 久久人人97超碰香蕉20202| 国产欧美日韩一区二区三区在线| 成人av一区二区三区在线看| 男人舔女人的私密视频| 97碰自拍视频| 精品人妻1区二区| 嫩草影视91久久| 亚洲成人精品中文字幕电影| 免费在线观看影片大全网站| 18禁观看日本| 涩涩av久久男人的天堂| 午夜福利18| 久久狼人影院| 在线免费观看的www视频| 久久中文看片网| 久久香蕉国产精品| 日本vs欧美在线观看视频| 看黄色毛片网站| 欧美日本亚洲视频在线播放| 亚洲人成77777在线视频| 国产亚洲欧美98| 熟妇人妻久久中文字幕3abv| 日本 av在线| 一级,二级,三级黄色视频| 国产日韩一区二区三区精品不卡| 国产一区二区三区视频了| 欧美激情 高清一区二区三区| a在线观看视频网站| 色综合欧美亚洲国产小说| 国产熟女xx| 欧美+亚洲+日韩+国产| 久久久水蜜桃国产精品网| 男人的好看免费观看在线视频 | 精品国产超薄肉色丝袜足j| 男人舔女人下体高潮全视频| 满18在线观看网站| 女人被狂操c到高潮| 国产成+人综合+亚洲专区| av中文乱码字幕在线| avwww免费| 男女之事视频高清在线观看| 非洲黑人性xxxx精品又粗又长| 99精品在免费线老司机午夜| 精品国产超薄肉色丝袜足j| 一进一出好大好爽视频| 亚洲国产欧美一区二区综合| 午夜福利成人在线免费观看| 欧美 亚洲 国产 日韩一| 久久精品91蜜桃| 伊人久久大香线蕉亚洲五| 亚洲熟妇熟女久久| 亚洲五月婷婷丁香| 日韩欧美在线二视频| 亚洲最大成人中文| 成人三级黄色视频| 日本免费a在线| 免费不卡黄色视频| 久久青草综合色| 久久婷婷成人综合色麻豆| 天堂√8在线中文| 国产欧美日韩一区二区精品| 我的亚洲天堂| 亚洲一区二区三区色噜噜| 1024香蕉在线观看| 高清毛片免费观看视频网站| 国产精品一区二区在线不卡| 男男h啪啪无遮挡| 亚洲五月色婷婷综合| 好男人电影高清在线观看| 国产午夜精品久久久久久| 69av精品久久久久久| 国产熟女午夜一区二区三区| 久久人人爽av亚洲精品天堂| 精品高清国产在线一区| 999久久久国产精品视频| 成人亚洲精品一区在线观看| 男人的好看免费观看在线视频 | 美女午夜性视频免费| 国产成人精品无人区| 午夜福利成人在线免费观看| 欧美不卡视频在线免费观看 | 自拍欧美九色日韩亚洲蝌蚪91| 超碰成人久久| 亚洲伊人色综图| 女性生殖器流出的白浆| 男人操女人黄网站| 精品久久久久久久久久免费视频| 亚洲专区国产一区二区| 国产私拍福利视频在线观看| 成人手机av| 精品久久久久久成人av| 日韩国内少妇激情av| 精品国产一区二区三区四区第35| 十八禁网站免费在线| 日韩欧美国产一区二区入口| 欧美日本中文国产一区发布| 亚洲黑人精品在线| 可以免费在线观看a视频的电影网站| 免费在线观看视频国产中文字幕亚洲| 无遮挡黄片免费观看| 两个人视频免费观看高清| 狂野欧美激情性xxxx| 两性午夜刺激爽爽歪歪视频在线观看 | 午夜两性在线视频| 天堂动漫精品| 欧美日韩福利视频一区二区| 91精品国产国语对白视频| 国产精品影院久久| 国产精品精品国产色婷婷| 国产亚洲欧美精品永久| 午夜激情av网站| 国产精品免费一区二区三区在线| 国产亚洲av嫩草精品影院| 久久国产亚洲av麻豆专区| 精品久久久久久,| 亚洲黑人精品在线| 欧美大码av| 激情视频va一区二区三区| 十八禁网站免费在线| 国产成人一区二区三区免费视频网站| 韩国av一区二区三区四区| 神马国产精品三级电影在线观看 | netflix在线观看网站| 村上凉子中文字幕在线| 欧美成人性av电影在线观看| 久久伊人香网站| 欧美成人午夜精品| 美女高潮喷水抽搐中文字幕| 麻豆国产av国片精品| 亚洲av片天天在线观看| 人人妻人人澡欧美一区二区 | 色播在线永久视频| 亚洲成人久久性| 成人亚洲精品一区在线观看| 精品日产1卡2卡| 亚洲电影在线观看av| 精品久久久久久久毛片微露脸| 精品欧美一区二区三区在线| 久热爱精品视频在线9| 欧美激情极品国产一区二区三区| 免费不卡黄色视频| 国产精品98久久久久久宅男小说| 丰满的人妻完整版| 久久人人97超碰香蕉20202| 亚洲精品中文字幕一二三四区| 亚洲中文日韩欧美视频| 国产精品免费一区二区三区在线| 午夜老司机福利片| 一区二区三区激情视频| 少妇裸体淫交视频免费看高清 | 日韩中文字幕欧美一区二区| 91av网站免费观看| 久久青草综合色| 日韩欧美三级三区| 校园春色视频在线观看| 色哟哟哟哟哟哟| 国产一区二区激情短视频| www国产在线视频色| 成人特级黄色片久久久久久久| 亚洲中文av在线| 成人三级做爰电影| 动漫黄色视频在线观看| 精品久久久久久久毛片微露脸| 午夜久久久在线观看| 亚洲中文字幕一区二区三区有码在线看 | 老司机在亚洲福利影院| 99精品欧美一区二区三区四区| 久热爱精品视频在线9| 国产xxxxx性猛交| 青草久久国产| 99久久99久久久精品蜜桃| av欧美777| 在线观看免费午夜福利视频| 曰老女人黄片| 91精品三级在线观看| 女人被狂操c到高潮| 身体一侧抽搐| 51午夜福利影视在线观看| 50天的宝宝边吃奶边哭怎么回事| 天天躁狠狠躁夜夜躁狠狠躁| 久久午夜亚洲精品久久| 一本大道久久a久久精品| 免费av毛片视频| 成人国语在线视频| 成人特级黄色片久久久久久久| 国产成人精品在线电影| 色av中文字幕| 999久久久国产精品视频| av网站免费在线观看视频| 日日摸夜夜添夜夜添小说| 午夜a级毛片| 欧美成狂野欧美在线观看| 欧美乱码精品一区二区三区| 人人妻,人人澡人人爽秒播| 国产三级黄色录像| 午夜精品在线福利| 人人澡人人妻人| 亚洲成人免费电影在线观看| 少妇被粗大的猛进出69影院| 免费搜索国产男女视频| 亚洲第一电影网av| 黄色成人免费大全| 久久精品aⅴ一区二区三区四区| 两个人视频免费观看高清| 老司机午夜福利在线观看视频| 成人永久免费在线观看视频| 久久久国产欧美日韩av| 黄色视频,在线免费观看| 免费不卡黄色视频| 一本综合久久免费| 美女大奶头视频| 国产高清有码在线观看视频 | 国产成人影院久久av| 国产亚洲精品一区二区www| 亚洲,欧美精品.| 精品日产1卡2卡| 97碰自拍视频| 午夜精品在线福利| 老司机在亚洲福利影院| 亚洲性夜色夜夜综合| 看免费av毛片| 男女做爰动态图高潮gif福利片 | 亚洲国产日韩欧美精品在线观看 | 久久 成人 亚洲| 国产精品久久久久久精品电影 | 欧美中文日本在线观看视频| 亚洲少妇的诱惑av| 青草久久国产| 国产精品免费视频内射| 午夜亚洲福利在线播放| 国产激情欧美一区二区| 亚洲男人天堂网一区| 欧洲精品卡2卡3卡4卡5卡区| 国产男靠女视频免费网站| 日本 欧美在线| 曰老女人黄片| 啪啪无遮挡十八禁网站| 99国产精品99久久久久| 美女高潮到喷水免费观看| 久久影院123| 国产xxxxx性猛交| 桃色一区二区三区在线观看| 欧美亚洲日本最大视频资源| 国产精品久久视频播放| 久久精品国产综合久久久| 黄色视频不卡| 国产精品一区二区在线不卡| 日韩 欧美 亚洲 中文字幕| 亚洲成a人片在线一区二区| 午夜福利高清视频| 一区二区三区激情视频| 女性生殖器流出的白浆| 国内精品久久久久精免费| 久久国产精品男人的天堂亚洲| 亚洲欧美激情综合另类| 男人舔女人的私密视频| 欧美乱妇无乱码| 黄网站色视频无遮挡免费观看| 久久亚洲真实| 色老头精品视频在线观看| 大陆偷拍与自拍| 黑人操中国人逼视频| 国产亚洲欧美在线一区二区| 久久热在线av| 一进一出好大好爽视频| 精品日产1卡2卡| 免费少妇av软件| 久久久国产成人精品二区| 老熟妇仑乱视频hdxx| 亚洲av第一区精品v没综合| 久久午夜综合久久蜜桃| 国产高清视频在线播放一区| 精品免费久久久久久久清纯| 久久久久精品国产欧美久久久| 亚洲成人久久性| 午夜免费激情av| 免费一级毛片在线播放高清视频 | 黄色视频,在线免费观看| 婷婷丁香在线五月| av免费在线观看网站| 人人妻人人爽人人添夜夜欢视频| 欧美日韩福利视频一区二区| 欧美 亚洲 国产 日韩一| 国产91精品成人一区二区三区| 无人区码免费观看不卡| 国产精品乱码一区二三区的特点 | 亚洲精品久久国产高清桃花| 中文字幕久久专区| 美女免费视频网站| 制服诱惑二区| 女人爽到高潮嗷嗷叫在线视频| 久久狼人影院| 精品熟女少妇八av免费久了| 在线十欧美十亚洲十日本专区| 90打野战视频偷拍视频| 国产人伦9x9x在线观看| 国产真人三级小视频在线观看| 18禁黄网站禁片午夜丰满| 亚洲中文字幕一区二区三区有码在线看 | √禁漫天堂资源中文www| 国产成年人精品一区二区| 日日爽夜夜爽网站| 色综合站精品国产| 亚洲一区二区三区不卡视频| 一夜夜www| 黄色成人免费大全| 一区二区三区精品91| 身体一侧抽搐| 很黄的视频免费| 男女午夜视频在线观看| 免费不卡黄色视频| 97人妻精品一区二区三区麻豆 | 十八禁网站免费在线| 国产精品一区二区免费欧美| 女性生殖器流出的白浆| 精品国产国语对白av| 亚洲一码二码三码区别大吗| 久久国产精品影院| 69精品国产乱码久久久| 午夜免费成人在线视频| 一卡2卡三卡四卡精品乱码亚洲| 久久狼人影院| 亚洲国产精品久久男人天堂| 岛国视频午夜一区免费看| 高潮久久久久久久久久久不卡| 在线观看一区二区三区| 国产精品永久免费网站| 日韩欧美三级三区| 人妻丰满熟妇av一区二区三区| 女生性感内裤真人,穿戴方法视频| 美国免费a级毛片| 欧美在线一区亚洲| 黄网站色视频无遮挡免费观看| 99精品久久久久人妻精品| 日韩精品青青久久久久久| 首页视频小说图片口味搜索| 国产xxxxx性猛交| 成人亚洲精品av一区二区|