• <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
    日韩不卡一区二区三区视频在线| 日韩欧美一区视频在线观看| 欧美97在线视频| 巨乳人妻的诱惑在线观看| 丝袜喷水一区| 少妇被粗大的猛进出69影院| 看十八女毛片水多多多| 午夜福利在线观看免费完整高清在| 美女高潮到喷水免费观看| 亚洲国产日韩一区二区| 久久狼人影院| 天天影视国产精品| 香蕉精品网在线| 狂野欧美激情性bbbbbb| 亚洲av免费高清在线观看| 一本大道久久a久久精品| 最近手机中文字幕大全| 久久国内精品自在自线图片| 不卡av一区二区三区| 春色校园在线视频观看| 韩国av在线不卡| 日韩精品免费视频一区二区三区| 交换朋友夫妻互换小说| 夫妻午夜视频| 99热全是精品| 99re6热这里在线精品视频| 日韩av在线免费看完整版不卡| 亚洲成国产人片在线观看| 精品少妇内射三级| 自线自在国产av| 午夜福利,免费看| 久久久久久久久久久免费av| 老熟女久久久| 亚洲内射少妇av| 最近最新中文字幕免费大全7| 国产精品人妻久久久影院| 久久人人爽av亚洲精品天堂| 多毛熟女@视频| 少妇人妻 视频| 亚洲美女黄色视频免费看| 久久久亚洲精品成人影院| 国产精品亚洲av一区麻豆 | 男女边吃奶边做爰视频| 国产乱来视频区| 成年动漫av网址| 亚洲国产精品一区二区三区在线| 日本vs欧美在线观看视频| videossex国产| 亚洲美女搞黄在线观看| 丝袜美腿诱惑在线| 2018国产大陆天天弄谢| 免费黄色在线免费观看| 亚洲精华国产精华液的使用体验| 国产黄频视频在线观看| 日日摸夜夜添夜夜爱| 久久久久久人妻| 国产一区二区激情短视频 | 欧美精品高潮呻吟av久久| 欧美日韩视频高清一区二区三区二| 一二三四在线观看免费中文在| 成年女人在线观看亚洲视频| 2022亚洲国产成人精品| 免费观看性生交大片5| 亚洲av综合色区一区| 亚洲在久久综合| 国产精品一国产av| 毛片一级片免费看久久久久| 国产一区有黄有色的免费视频| 欧美激情高清一区二区三区 | 91在线精品国自产拍蜜月| 亚洲成人一二三区av| 高清视频免费观看一区二区| 黄片无遮挡物在线观看| 日日摸夜夜添夜夜爱| 国产精品蜜桃在线观看| videosex国产| 精品午夜福利在线看| 下体分泌物呈黄色| 一本一本久久a久久精品综合妖精 国产伦在线观看视频一区 | 嫩草影院入口| 国产成人精品一,二区| 亚洲av日韩在线播放| 大片电影免费在线观看免费| 麻豆精品久久久久久蜜桃| 国产精品不卡视频一区二区| 少妇被粗大猛烈的视频| 夜夜骑夜夜射夜夜干| 午夜影院在线不卡| 老鸭窝网址在线观看| 多毛熟女@视频| 一级黄片播放器| 伦精品一区二区三区| 亚洲av欧美aⅴ国产| 成人黄色视频免费在线看| 男女国产视频网站| 久久这里只有精品19| 日韩中字成人| 日韩熟女老妇一区二区性免费视频| 边亲边吃奶的免费视频| 日韩一本色道免费dvd| 婷婷色麻豆天堂久久| 汤姆久久久久久久影院中文字幕| 国产精品国产三级专区第一集| 观看av在线不卡| 亚洲av成人精品一二三区| 18禁裸乳无遮挡动漫免费视频| 亚洲av福利一区| 天堂8中文在线网| 欧美日韩国产mv在线观看视频| 午夜激情av网站| 最近的中文字幕免费完整| 一级爰片在线观看| 母亲3免费完整高清在线观看 | 极品少妇高潮喷水抽搐| 男女啪啪激烈高潮av片| 我的亚洲天堂| 亚洲国产精品国产精品| 免费高清在线观看日韩| 欧美另类一区| 熟妇人妻不卡中文字幕| 国产女主播在线喷水免费视频网站| 国产亚洲最大av| 久久久久久免费高清国产稀缺| 美女高潮到喷水免费观看| 亚洲av.av天堂| 在现免费观看毛片| 不卡视频在线观看欧美| 日韩精品有码人妻一区| 下体分泌物呈黄色| 欧美激情 高清一区二区三区| 哪个播放器可以免费观看大片| 91精品伊人久久大香线蕉| 涩涩av久久男人的天堂| 亚洲精品,欧美精品| 一二三四中文在线观看免费高清| 成年动漫av网址| 国产黄频视频在线观看| 国语对白做爰xxxⅹ性视频网站| 国产xxxxx性猛交| 一级片免费观看大全| 青春草亚洲视频在线观看| 国产精品久久久久久久久免| 亚洲国产精品一区二区三区在线| 亚洲 欧美一区二区三区| 伊人久久大香线蕉亚洲五| 你懂的网址亚洲精品在线观看| www.精华液| 亚洲av综合色区一区| 夫妻性生交免费视频一级片| 青春草亚洲视频在线观看| 久久久久久人妻| 男女免费视频国产| 久久99热这里只频精品6学生| 欧美成人午夜免费资源| 成人毛片a级毛片在线播放| 久久午夜综合久久蜜桃| 波多野结衣av一区二区av| 电影成人av| 2018国产大陆天天弄谢| 午夜久久久在线观看| 极品少妇高潮喷水抽搐| 久久精品熟女亚洲av麻豆精品| 亚洲三级黄色毛片| 亚洲成av片中文字幕在线观看 | 午夜福利乱码中文字幕| 秋霞伦理黄片| 伦精品一区二区三区| 永久网站在线| 亚洲美女黄色视频免费看| 激情五月婷婷亚洲| 啦啦啦在线免费观看视频4| 99热全是精品| 天天影视国产精品| 国产精品.久久久| 日韩av不卡免费在线播放| 99热国产这里只有精品6| 国产男人的电影天堂91| 午夜91福利影院| 晚上一个人看的免费电影| 欧美激情高清一区二区三区 | 久久人人爽人人片av| 免费大片黄手机在线观看| 十分钟在线观看高清视频www| 亚洲国产精品成人久久小说| 国产精品99久久99久久久不卡 | av在线观看视频网站免费| 美女福利国产在线| 亚洲精华国产精华液的使用体验| 国产精品人妻久久久影院| 国产无遮挡羞羞视频在线观看| 久久久国产精品麻豆| xxx大片免费视频| 精品国产超薄肉色丝袜足j| 久久女婷五月综合色啪小说| 国产精品女同一区二区软件| 国产免费又黄又爽又色| 人人妻人人澡人人看| 亚洲国产色片| 国产成人精品在线电影| 国产精品免费大片| 久久久久久久国产电影| 亚洲av中文av极速乱| 国产免费视频播放在线视频| 亚洲,欧美精品.| 亚洲国产av新网站| 精品国产超薄肉色丝袜足j| 久久国产精品大桥未久av| 丰满乱子伦码专区| 老熟女久久久| 久久人妻熟女aⅴ| 亚洲综合色网址| 人妻少妇偷人精品九色| 亚洲欧美色中文字幕在线| 视频在线观看一区二区三区| 婷婷色综合www| 尾随美女入室| 色婷婷av一区二区三区视频| 99热全是精品| 亚洲av欧美aⅴ国产| 美女福利国产在线| 国产精品二区激情视频| 伦理电影免费视频| 日本欧美视频一区| 国产麻豆69| 一级a爱视频在线免费观看| 久久久国产精品麻豆| 国产亚洲最大av| 精品亚洲成a人片在线观看| 91国产中文字幕| 春色校园在线视频观看| videosex国产| 成人18禁高潮啪啪吃奶动态图| 日日啪夜夜爽| 日韩中字成人| 性色av一级| 大话2 男鬼变身卡| 久久精品国产综合久久久| 久久av网站| 亚洲av福利一区| 国产欧美日韩一区二区三区在线| a级毛片在线看网站| 自拍欧美九色日韩亚洲蝌蚪91| 亚洲成人av在线免费| 一级片免费观看大全| 色视频在线一区二区三区| 欧美 日韩 精品 国产| 免费大片黄手机在线观看| 国产成人91sexporn| 成年女人毛片免费观看观看9 | 久久午夜福利片| 蜜桃在线观看..| 少妇的丰满在线观看| 国产精品蜜桃在线观看| 亚洲一码二码三码区别大吗| 日韩免费高清中文字幕av| 十八禁网站网址无遮挡| 多毛熟女@视频| 国产亚洲av片在线观看秒播厂| 成人黄色视频免费在线看| 97人妻天天添夜夜摸| 黑丝袜美女国产一区| 成人漫画全彩无遮挡| 欧美av亚洲av综合av国产av | 少妇猛男粗大的猛烈进出视频| 国产成人午夜福利电影在线观看| 一级片'在线观看视频| 免费黄网站久久成人精品| 一区二区三区精品91| 国产精品熟女久久久久浪| 大香蕉久久网| 天堂8中文在线网| 日韩一区二区视频免费看| 亚洲伊人久久精品综合| 亚洲国产成人一精品久久久| 亚洲天堂av无毛| 欧美精品亚洲一区二区| 亚洲av日韩在线播放| 大片免费播放器 马上看| 欧美激情 高清一区二区三区| 亚洲综合色网址| 伊人久久国产一区二区| 久久这里只有精品19| 另类精品久久| 1024视频免费在线观看| av女优亚洲男人天堂| 熟女电影av网| 一区二区三区乱码不卡18| 黄色配什么色好看| 国产综合精华液| av视频免费观看在线观看| 国产精品亚洲av一区麻豆 | 久久久久久久久久人人人人人人| av国产精品久久久久影院| 考比视频在线观看| 免费久久久久久久精品成人欧美视频| 欧美av亚洲av综合av国产av | 午夜福利影视在线免费观看| 99国产综合亚洲精品| 看免费av毛片| 欧美人与善性xxx| 国产亚洲最大av| 国产精品香港三级国产av潘金莲 | 日韩精品有码人妻一区| 久久久久久久国产电影| 街头女战士在线观看网站| a 毛片基地| 寂寞人妻少妇视频99o| 日韩欧美精品免费久久| 亚洲国产日韩一区二区| av免费在线看不卡| 午夜激情av网站| 国产深夜福利视频在线观看| 热re99久久精品国产66热6| 欧美黄色片欧美黄色片| 亚洲欧美成人综合另类久久久| 日韩免费高清中文字幕av| 久久婷婷青草| 欧美人与性动交α欧美精品济南到 | av.在线天堂| 日韩av免费高清视频| 免费黄色在线免费观看| 国产又色又爽无遮挡免| 大香蕉久久成人网| 日日摸夜夜添夜夜爱| 在线观看人妻少妇| 欧美精品一区二区大全| 亚洲欧洲精品一区二区精品久久久 | 欧美 日韩 精品 国产| 最近中文字幕2019免费版| 成人漫画全彩无遮挡| 另类亚洲欧美激情| 欧美精品av麻豆av| 午夜福利一区二区在线看| av网站在线播放免费| 午夜福利视频在线观看免费| 日韩大片免费观看网站| 亚洲欧美色中文字幕在线| 日日撸夜夜添| 涩涩av久久男人的天堂| 一本—道久久a久久精品蜜桃钙片| 日韩制服丝袜自拍偷拍| 搡女人真爽免费视频火全软件| 欧美成人精品欧美一级黄| 国产精品成人在线| 成人二区视频| 26uuu在线亚洲综合色| 国产精品蜜桃在线观看| 国产国语露脸激情在线看| 欧美日韩亚洲高清精品| 久久精品国产自在天天线| 国产精品久久久久久精品古装| 2022亚洲国产成人精品| 妹子高潮喷水视频| 日本色播在线视频| 18禁观看日本| 精品久久久久久电影网| 久久久a久久爽久久v久久| 日韩一卡2卡3卡4卡2021年| 狠狠精品人妻久久久久久综合| 美女主播在线视频| 伊人久久大香线蕉亚洲五| 久久久久精品久久久久真实原创| 黄色视频在线播放观看不卡| kizo精华| 男人舔女人的私密视频| 日本黄色日本黄色录像| 91aial.com中文字幕在线观看| 色播在线永久视频| 在线观看免费视频网站a站| 国产 一区精品| 男男h啪啪无遮挡| 午夜福利乱码中文字幕| 久久精品久久久久久噜噜老黄| 亚洲国产精品成人久久小说| 欧美激情极品国产一区二区三区| 午夜久久久在线观看| 丝瓜视频免费看黄片| 精品人妻偷拍中文字幕| 日日摸夜夜添夜夜爱| 精品人妻在线不人妻| 又大又黄又爽视频免费| 少妇熟女欧美另类| av.在线天堂| 高清不卡的av网站| 亚洲精华国产精华液的使用体验| 亚洲av.av天堂| 国产精品蜜桃在线观看| 日本av免费视频播放| 国产成人欧美| 亚洲av成人精品一二三区| 亚洲国产精品国产精品| xxx大片免费视频| 日韩电影二区| 中国国产av一级| 国产精品二区激情视频| 久久韩国三级中文字幕| 国产熟女欧美一区二区| 老司机影院成人| 男女下面插进去视频免费观看| av女优亚洲男人天堂| 男女啪啪激烈高潮av片| 婷婷成人精品国产| 亚洲熟女精品中文字幕| 久久精品国产a三级三级三级| 99热网站在线观看| 男人添女人高潮全过程视频| 久久午夜综合久久蜜桃| 成人午夜精彩视频在线观看| 亚洲美女黄色视频免费看| 日产精品乱码卡一卡2卡三| 性色avwww在线观看| 久久av网站| 一本—道久久a久久精品蜜桃钙片| 国产精品99久久99久久久不卡 | 国产免费又黄又爽又色| 婷婷色综合大香蕉| 一级a爱视频在线免费观看| 在线天堂最新版资源| 青青草视频在线视频观看| 日韩精品免费视频一区二区三区| 精品一区二区三区四区五区乱码 | 1024视频免费在线观看| 国产在视频线精品| 黄色一级大片看看| 精品第一国产精品| 免费少妇av软件| 国产极品粉嫩免费观看在线| 免费黄色在线免费观看| 成人毛片a级毛片在线播放| 免费黄色在线免费观看| 观看av在线不卡| 老熟女久久久| 欧美日韩亚洲国产一区二区在线观看 | 亚洲av电影在线观看一区二区三区| 亚洲久久久国产精品| 91精品国产国语对白视频| 国产av精品麻豆| 在线 av 中文字幕| 亚洲精品,欧美精品| 国产欧美日韩一区二区三区在线| 最近最新中文字幕免费大全7| 国产淫语在线视频| 天堂8中文在线网| 日韩熟女老妇一区二区性免费视频| 美国免费a级毛片| 人成视频在线观看免费观看| 亚洲精品在线美女| 日韩一卡2卡3卡4卡2021年| 少妇精品久久久久久久| av一本久久久久| 精品人妻熟女毛片av久久网站| 嫩草影院入口| 天天操日日干夜夜撸| 韩国av在线不卡| 少妇被粗大猛烈的视频| 一级爰片在线观看| 成人18禁高潮啪啪吃奶动态图| 国产在线视频一区二区| 亚洲欧美清纯卡通| 亚洲精品久久久久久婷婷小说| 久久久久久久大尺度免费视频| 亚洲美女视频黄频| 精品一区二区三卡| 欧美日韩视频高清一区二区三区二| 九草在线视频观看| 久久久久久久久免费视频了| 免费看不卡的av| 校园人妻丝袜中文字幕| 成人亚洲欧美一区二区av| 下体分泌物呈黄色| 国产熟女午夜一区二区三区| 涩涩av久久男人的天堂| 精品一区在线观看国产| 性高湖久久久久久久久免费观看| 女性被躁到高潮视频| 大香蕉久久成人网| 亚洲精品一区蜜桃| 欧美少妇被猛烈插入视频| 久久精品国产亚洲av涩爱| videossex国产| 欧美精品人与动牲交sv欧美| 欧美成人精品欧美一级黄| 99久久中文字幕三级久久日本| 男女无遮挡免费网站观看| 青青草视频在线视频观看| 国产男女超爽视频在线观看| 欧美日韩精品网址| 一区在线观看完整版| 精品人妻一区二区三区麻豆| av一本久久久久| 伦理电影大哥的女人| 亚洲四区av| 国产黄频视频在线观看| 十八禁高潮呻吟视频| 夫妻午夜视频| 日韩熟女老妇一区二区性免费视频| 老司机亚洲免费影院| 久久精品亚洲av国产电影网| 亚洲少妇的诱惑av| 成人影院久久| 成人亚洲欧美一区二区av| 日韩 亚洲 欧美在线| av福利片在线| 免费女性裸体啪啪无遮挡网站| 777久久人妻少妇嫩草av网站| 美女主播在线视频| 巨乳人妻的诱惑在线观看| 黑人欧美特级aaaaaa片| 熟女电影av网| 国产亚洲一区二区精品| 成人国语在线视频| 亚洲男人天堂网一区| 一区二区三区乱码不卡18| 亚洲av中文av极速乱| www.自偷自拍.com| 一级毛片电影观看| 亚洲成色77777| 日韩 亚洲 欧美在线| 久久久久精品久久久久真实原创| 午夜福利乱码中文字幕| 91精品国产国语对白视频| 纵有疾风起免费观看全集完整版| 国产一级毛片在线| 亚洲综合色网址| 国产在线视频一区二区| 成年人午夜在线观看视频| 午夜av观看不卡| 在线观看一区二区三区激情| 亚洲精品久久午夜乱码| 高清黄色对白视频在线免费看| 亚洲av综合色区一区| 可以免费在线观看a视频的电影网站 | 亚洲中文av在线| 欧美日韩av久久| 欧美日韩一级在线毛片| h视频一区二区三区| 久久97久久精品| 交换朋友夫妻互换小说| 免费黄网站久久成人精品| 多毛熟女@视频| 90打野战视频偷拍视频| 九色亚洲精品在线播放| 视频区图区小说| 午夜福利一区二区在线看| 大片免费播放器 马上看| 亚洲一区二区三区欧美精品| 女人高潮潮喷娇喘18禁视频| 美女xxoo啪啪120秒动态图| 国产av精品麻豆| 国产亚洲一区二区精品| 日韩中文字幕视频在线看片| 91精品三级在线观看| 五月天丁香电影| 亚洲天堂av无毛| 日韩电影二区| 热re99久久国产66热| 丝袜人妻中文字幕| 久久久a久久爽久久v久久| 人人妻人人添人人爽欧美一区卜| 高清黄色对白视频在线免费看| 中文字幕亚洲精品专区| 黑人猛操日本美女一级片| 夫妻性生交免费视频一级片| 建设人人有责人人尽责人人享有的| 亚洲av男天堂| 久久av网站| 老司机亚洲免费影院| 国产欧美日韩综合在线一区二区| 欧美变态另类bdsm刘玥| 精品久久久精品久久久| 如日韩欧美国产精品一区二区三区| 久久久a久久爽久久v久久| 日本爱情动作片www.在线观看| 天天操日日干夜夜撸| 久久精品国产亚洲av天美| 男女国产视频网站| 香蕉丝袜av| 晚上一个人看的免费电影| 免费av中文字幕在线| 黑丝袜美女国产一区| 超碰97精品在线观看| 在线观看三级黄色| videosex国产| 国产精品一二三区在线看| 国产精品人妻久久久影院| 亚洲中文av在线| 欧美精品亚洲一区二区| 亚洲熟女精品中文字幕| 99久久人妻综合| 视频在线观看一区二区三区| 日本猛色少妇xxxxx猛交久久| 国产午夜精品一二区理论片| 久久久久人妻精品一区果冻| 国产成人一区二区在线| 国产成人av激情在线播放| 亚洲中文av在线| 天天躁日日躁夜夜躁夜夜| 欧美日韩一区二区视频在线观看视频在线| 欧美国产精品va在线观看不卡| 日韩av在线免费看完整版不卡| 一二三四在线观看免费中文在| 91aial.com中文字幕在线观看| www.自偷自拍.com| 一本色道久久久久久精品综合| 国产av一区二区精品久久| 国产成人av激情在线播放| 欧美最新免费一区二区三区| 五月天丁香电影| 在线观看免费日韩欧美大片| 丝瓜视频免费看黄片| 大片电影免费在线观看免费| 尾随美女入室| 午夜日韩欧美国产| 欧美精品一区二区大全| 国产无遮挡羞羞视频在线观看|