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

    Numerical simulation and experimental research on an inductively coupled RF plasma cathode

    2022-02-15 11:08:20ZongqiXU徐宗琦PingyangWANG王平陽ZhiweiHUA華志偉ShiyuanCONG從拾源andShengnanYU余盛楠
    Plasma Science and Technology 2022年1期
    關(guān)鍵詞:平陽

    Zongqi XU (徐宗琦), Pingyang WANG (王平陽), Zhiwei HUA (華志偉),Shiyuan CONG (從拾源) and Shengnan YU (余盛楠)

    1 Institute of Engineering Thermophysics, School of Mechanical Engineering, Shanghai Jiao Tong University, Shanghai 200240, People’s Republic of China

    2 Key Laboratory of Materials Modification by Laser, Ion and Electron Beams (Ministry of Education),School of Physics, Dalian University of Technology, Dalian 116024, People’s Republic of China

    Abstract In this study,numerical simulation and discharge current tests were conducted on an inductively coupled radio frequency (RF) plasma cathode.Numerical simulations and experimental measurements were performed to study the factors influencing the electron extraction characteristics, including the gas type, gas flow, input power and extracting voltage.The simulation results were approximately consistent with the experimental results.We experimentally found that the RF input power mainly determines the extracted electron current.An electron current greater than 1 A was acquired at 270 W (RF input power), 2.766 sccm(xenon gas).Our results prove that an inductively coupled RF plasma cathode can be reasonable and feasible, particularly for low power electric propulsion devices.

    Keywords: RF discharge, RF plasma cathode, electron extraction characteristics, RF power,electron current

    1.Introduction

    Cathodes are important structural components primarily used in electric propulsion applications to deliver an electron beam,particularly in Hall thrusters and ion thrusters.For example,cathodes restrict the lifetime of thrusters and determine whether satellites can complete space missions optimally.An object or a device that can emit electrons is called a neutralizer.Based on different electron emission principles,neutralizers can be divided into hot cathode, hollow cathode,and plasma cathodes, for example, radio frequency (RF)plasma cathode and electron cyclotron resonance (ECR) discharge cathode [1].

    A hot cathode generally refers to a pure metal that emits electrons at a high temperature, with the advantages of stable emission, resistance to high-pressure ion bombardment, and lack of reactivation required after atmospheric exposure.Tungsten is typically used as an electron emitter owing to its high melting point,low work function,high emission and low evaporation rate.However, it can not only be processed into filaments or bands, but can also chemically react with oxides at high temperatures.With the development of science and technology, hollow cathodes [2] with thermionic emitters made from low-work-function materials (e.g.LaB6) or composites (e.g.porous tungsten impregnated with a low-workfunction composite of barium oxide, calcium oxide, and aluminum oxide) are being extensively used in plasma propulsion owing to their high current densities and low cathodefall voltage [3].However, some issues associated with the hollow cathodes have seriously affected their utilization,such as difficulty in acquiring the insert material, considerable amount of time required to preheat before operation and cooling after operation, discharge oscillation due to temperature drop [4], low electron release efficiency of the emitter during operation and poor resistance to poisoning.Therefore,hollow cathodes have the characteristics of high production cost, limited service lifetime and delayed ignition process.In recent years,there have been significant efforts on developing other electron sources as practical alternatives, while plasma cathodes have attracted much attention.Different plasma methods such as RF (capacitively coupled plasma [5] or inductively coupled plasma (ICP)) [6], micro-wave ECR)[3, 7, 8], and helicon [9, 10] have been employed in plasma cathode devices.Among plasma cathodes, microwave and helicon neutralizers are equipped with a complicated electron extraction device and require an external magnetic field that consumes additional power.Conversely, RF excitation is considered a preferable plasma generation method for cathode applications owing to its characteristics of simple structure,low-power loss and high plasma density, particularly in the ICP discharge mode [11].

    Inductively coupled RF plasma cathode has been the focus of research for more than a decade both theoretically and experimentally.Watanabe et al [12] designed an ICP cathode for ion thrusters and presented a series of experimental studies on its performances.The emission electron current reached approximately 2.1 A with an RF power of 100 W, Xe mass flow rate of 2 sccm, and target voltage of 40 V.Subsequently,an emission electron current of 3.3 A was obtained in a combination test involving an anode-layer-type Hall thruster and an RF plasma cathode [13, 14].They also found a balance between the performances of plasma ignition and electron emission and determined the discharge vessel diameter [15].Scholze et al described the design of an inductive coupled RF plasma bridge neutralizer with a maximum electron current of 1.6 A at a Xe gas flow rate of 0.2 sccm and a power of 300 W [16].They developed a performance model to calculate the electrical parameters and provided the optimization ideas[17].Godyak et al[3]found that in the RF plasma discharge experiments, the maximum emission efficiency of the ICP was higher than that of microwave and helicon.Other studies concerning the RF ignition of hollow cathodes can also be found in their papers[18–20].Jahanbakhsh et al discussed the relationship between the electron extraction capacity and the gas flow rate,the size of the extraction aperture, and the absorbed RF power[10, 21].According to another theory, the anode spot has an impact on the electron extraction property [22, 23].A comparison of numerical and experimental results showed that the effects of target voltage, aperture diameter, gas flow rate and RF power were consistent [24].In addition, a ferromagnetic core was utilized for enhancing the power transfer efficiency of ICP by the Bogazici University Space Technologies Laboratory[25].Another corrosive gas,such as iodine vapor,was also considered as the working gas in the RF cathode [26, 27].

    However,there are no investigations on the ion collectors which have an important influence on the electron extraction.In this study, we designed a type of inductively coupled RF plasma cathode with a cylindrical ion collector and analyzed the influences of the gas flow rate,RF power and bias voltage on the current extraction through simulations and experiments.In section 2, we describe the characteristics of the inductively coupled RF plasma cathodes in detail.Section 3 provides the description of a geometric model and presents results and discussion of the numerical simulations.Section 4 describes the experimental setups and presents the experimental results and corresponding discussions.Finally, in section 5, we summarize our findings and compare the two results.

    2.Inductively coupled RF plasma cathode

    Figures 1 and 2 show the RF plasma cathode designed for this investigation.The plasma chamber for the discharge is made into a cylindrical structure, and its inner diameter and length are 40 mm and 76 mm,respectively.The wall of the discharge chamber is made of alumina, which has the properties of insulation and high temperature resistance.The orifice plate is made of graphite and there is a 2 mm diameter orifice at the center of the plate for electron emission.Considering the skin effect of high-frequency transmission current, a 5 mm diameter copper pipe is wrapped around the discharge chamber as the induction coil.The cylindrical collector made of graphite is placed on the axis of the discharge chamber as an ion collector, and it has six small holes whose diameter is 1 mm along the periphery for gas to flow into the chamber.

    Because the ignition of ICP obeys Paschen’s law, the orifice diameter is much smaller than the inner diameter of the chamber,and this helps sustain a sufficiently high pressure in the RF plasma cathode for the self-ignition of plasma.To maintain quasi-neutrality in the ICP discharge, an equivalent quantity of ions should be collected at the ion collector in the cathode when electrons are emitted from the orifice.Thus, a graphite collector, which is inserted into the chamber, is necessary for the steady operation of the RF plasma cathode.The shape of the ion collector has an important influence on the ion collection ability and further determines the values of the emitted electric current.

    3.Method, results, and discussion of numerical simulations

    3.1.Solution strategy

    In this study, a 2D geometric model is employed to investigate the RF plasma cathode with xenon gas.Three coupled processes are considered in the model: particle transport,neutral gas flow and electromagnetic induction.They are described by typical drift-diffusion approximation equations,a set of N–S equations and Maxwell equations, respectively.A series of assumptions is made as follows to simplify the analysis.

    (a) The induction currents are solved in the frequency domain.

    (b) The electron energy distribution function (EEDF) is assumed to obey the Maxwellian distribution.

    (c) The working gas is considered as an ideal gas.It is governed by the state equation and compressible.

    (d) Because of the effective energy exchange between ions and neutral particles, their temperatures are equal.

    Figure 1.Schematic of RF plasma cathode.

    Figure 2.Side view of RF plasma cathode and target.

    (e) The effect of a self-induced electromagnetic field is negligible

    (f) The electronically excited xenon species are lumped into a single species.

    The equations and boundary conditions are presented in the following subsections and solved numerically using the finite element method and the standard Galerkin method for space variation and back differentiation method for time variation.The commercial software package COMSOL Multiphysics is employed as the solver, which has been successfully demonstrated to be applicable for simulating gas discharge processes [28–30].Triangular and boundary layer meshes are employed in this model.To obtain a fast convergence of the calculation,the particle transport equations are discretized in a logarithmic form using the linear Lagrange shape function.To accelerate the model convergence, a transient solver is employed and the quasi-steady results are discussed [31].

    3.2.Governing equations

    Based on the aforementioned assumptions, the transport processes of the electrons and heavy particles are described by typical drift-diffusion approximation equations.The electron and electron energy transport equations are as follows[32]:

    whereneandnεare the electron number density and electron energy density, respectively (i.e.nε=ne,whereis the mean electron energy, and the electron temperature isTe= 2εˉ3); Γeand Γεare respectively the el ectron flux and electron energy flux;eμandμεare,respectively, the electron mobility and electron energy mobility;DeandDεare,respectively, the electron diffusion coefficient and electron energy diffusion coefficient.According to the assumption of the Maxwellian distribution of the EEDF, the mobility and diffusion coefficients satisfy Einstein’s relationship, i.e.με=5μe3,De=μeTe,Dε=με Tε.Therefore, only the electron mobilityeμis required to be solved by the Boltzmann solver (i.e.BOLSIG+) [33].The collision source termSerepresents the production or loss of electrons, and the electron collision reactions with the heavy particles are summarized in table 1.The energy source termSεrepresents the complete energy change in the entire collision reaction.ugis the flow velocity of the neutral gas, andEis the electric field intensity.SeandSεare described as follows:

    wherexjis the mole fraction of the species for reactionj;kjis the rate coefficient,calculated by the Boltzmann solver;Nnis the total number density of the neutral particle;jεΔ is the energy loss from the reactionj.

    The heavy particles also have a similar transport equation as follows [31]:

    where ρ denotes the mass density of the working gas;wiis the mass fraction of the ith species;Siis a collision term of the heavy particles derived from the plasma chemical reactions in table 1;iΓ is the flux of the ith species,i,mμandDi,mare,respectively, the average mobility and diffusion coefficients of the ith species;Ziis the charge number of the ith species.The subscript m refers to the average over all components of the mixing working gas, including ground atoms, excited atoms,and ions.Moreover,Mˉ refers to the mean molar massof the mixing working gas [31], expressed as follows:

    Table 1.Plasma chemistry processes used in the model.

    where N is the total number of species in the gas mixture;Miis the molar mass of the ith species.Among them, the state atom has the maximum proportion, and the mass fraction expression is

    The electric field intensity is calculated by combining the particle transport equations with Poisson’s equation [31]

    whereφ,ε0,ni,and e denote the plasma potential,permittivity of vacuum, ion number density, and elementary charge,respectively.

    3.3.Gas flow equations

    The flow of the neutral gas in the calculation domain is considered a laminar flow model, which is described by a set of N–S equations.The continuity and momentum balance equations are given as follows [37]:

    where p is the absolute pressure of the gas mixture;dμis the dynamic viscosity; I is the unit tensor.The pressure gradient force and viscous stress tensor are on the right-hand side of equation (13).The superscript T represents the matrix transposition.

    3.4.Electromagnetic induction equation

    For the nonmagnetized and nonpolarized plasma, the induction currents are computed in the frequency domain using the following equation:

    Here, jkis the imaginary part(i.e.j2

    k= ?1);ωis the angular frequency; A is the magnetic vector potential;0μis conductivity of vacuum.The plasma conductivity should be specified as a material property, typically from cold plasma approximation:

    where q is the electron charge;meis the electron mass;eνis the collision frequency.

    3.5.Geometric model description and boundary conditions

    Because the designed cathode is symmetrical about its central axis, the calculation domain can be set one-half of its cross section.As detailed in figure 3, domain d.1 is the anode target; domain d.2 is the RF cathode wall; domain d.3 is the collector; domain d.4 is the RF antenna; boundary b.1 is the gas inlet; boundary b.2 is the plasma outlet; boundary b.3 is the inner wall;b.4 is the symmetric axis along the direction of z-axis; b.5 is the collector wall; b.6 is the anode target wall.The length and radius of the discharge chamber along the axis are 80 mm and 20 mm, respectively.The thickness is 4 mm.The gap between the anode target and the RF cathode tip is 10 mm.The length and width of the calculation domain are 104 and 55 mm.The diameters of the gas inlet and plasma outlet are respectively 1 mm and 2 mm.

    Figure 3.Geometric model of RF plasma cathode.

    The fluxes of the electron, electron energy, and heavy particle have boundary conditions as follows [31]:where n is the normal direction of the boundaries; k is the Boltzmann constant;Teis the electron temperature;iεis the mean energy of the secondary electrons;iγis the secondary electron emission coefficient, which is determined by the incident ions and material surface properties.In our study,only the secondary electron emitted from the collector surface was considered because the ion bombardment dominantly occurs on the collector wall.Therefore, an empirical value= 0.25was set.Rs,idetailed in [38] is the rate of surface reaction on the wall in the first term on the right side of equation (18), and the second term is the ionic migration on account of the electric field acting toward the wall, which is satisfied if the expressionZ i(E·n) > 0.The collector is grounded, i.e.φ= 0,and the anode is connected to a power supply of which the detailed expression is as follows [31]:

    whereJi,Je,andJdare, respectively, the current densities of the ion,electron and displacement,respectively;Idis the total discharge current applied to the anode target boundary.The symbol Ω refers to the integral boundary on it,and the surface area of the anode target will be automatically calculated after a symmetric rotation.

    For the gas flow equations,the collector surface and inner wall are considered as no-slip condition.The gas inlet and gas outlet respectively have mass flow rate and fixed-pressure boundary conditions.For the electromagnetic induction equation,the normal direction of the calculation boundaries is parallel to that of the magnetic vector potential.The other boundaries in the calculated domain satisfy Maxwell’s equations.Table 2 presents all the boundary conditions.

    3.6.Characteristics of plasma parameters

    The RF power supply drives the antenna coil, external to the plasma and separated from it by the cathode wall.The RF current flowing in the coil launches an evanescent disturbance that decays over a distance of centimeters into the plasma[39].This induces RF current in the plasma and transfers energy to electrons.The plasma discharge is sustained after the RF electromagnetic is coupled with the electrons in the plasma.The energy transfer phenomenon is called ohmic heating or collisional heating [39].In low-pressure inductively coupled RF plasma additional electron heating mechanism exists that is collisionless heating, which is dominant due to the too low collision frequencies for efficient ohmic heating of electrons [40, 41].The heated electrons react with neutral particles by elastic, ionization, and excitation collisions.

    The RF current in the antenna also generates a weak induced magnetic field along the axial direction, which has a significant impact on the parameter distributions of the charged particles in the discharge chamber.Figure 4 shows the magnetic field profiles.

    The maximum magnetic flux density appears in the central region near the RF coil, as shown in figure 4, where most of the electrons are trapped and the neutral particles have a high ionization rate.

    The distributions of the various plasma parameters in the discharge chamber are acquired through numerical simulation,which is helpful for analyzing the performance of the RF plasma cathode.The main simulation conditions are as follows: xenon gas flow rate is 5 sccm; extracted voltage is 100 V;RF power is 100 W.Figure 5 shows the distribution of electron number density.The figure shows that the maximum electron number density occurs near the outlet of the discharge chamber and the induction coil,and the value is above 2.5 × 1018m?3.As mentioned before, this is because the electrons are bound in the central region of the discharge chamber by the induced magnetic field.Moreover, the electrons are also attracted by the anode target and emitted from the outlet of the discharge chamber.This indicates that a high plasma density can be acquired by ICP discharge; thus, sufficient electrons can be produced and extracted to form an electric current.

    Figure 6 shows the electric potential distribution.Clearly,the potential changes sharply near the collector.Since the collector is biased negatively, a large potential gradient is more helpful for ions to be absorbed.The electric potential is maintained at approximately 100 V between the discharge chamber and the anode target.The electrons can be easy to be extracted under a potential condition.

    3.7.Extracted current characteristics

    Figure 7 presents the response time of the extracted current between the anode target and the collector.As shown in the figure, the current becomes stable in a short time, and the value is approximately 0.3369 A.This result indicates that thedesigned structure of the RF plasma cathode meets the requirement of acquiring a stable extracted current.

    Table 2.Boundary conditions of model.The numbers (b.1–b.7) indicate the boundaries shown in figure 3.

    Figure 4.Distribution of the magnetic flux density in the calculation domain.

    Figure 5.Distribution of the electron number density in the calculation domain.

    Figure 6.Distribution of the electric potential in the calculation domain.

    Figure 7.Response time of extracted current.

    Figure 8.Experimental system of RF plasma cathode.

    4.Method, results, and discussion of experiments

    As described above, the operation of the RF plasma cathode comprises two parts: plasma generation and electron extraction.Experimental studies have been conducted to determine the relationships between the electron extraction and gas flow rate, RF power, and bias voltage.The mass flow controller used in the experiment is M3030 V/A (N2equivalent).The gas factor corresponding to Xe is 1.383.The experimental parameters are set as follows: the gas flow reading of the controller varied from 2.0 to 4.0 sccm at intervals of 0.5 sccm,and the RF input power varied from 50 to 300 W at intervals of 10 W.The bias voltage varied from 0 to 100 V at intervals of 10 V, and the electron current extracted by the target is measured.

    4.1.Experimental setup

    Figure 8 shows the experimental system.All the experiments are conducted in a vacuum chamber, whose diameter and length are 0.6 m and 1.2 m, respectively.The vacuum chamber reaches approximately 1.0 × 10?4Pa without any load using two turbo-molecular pumps and two rotary vane pumps.The RF cathode and the anode target are placed in the chamber, whereas an RF generator, an impedance matching box,and a DC power supply are attached from outside of the vacuum chamber.The anode target is located approximately 10 mm downstream of the orifice plate.The anode target and ion collector are biased positively and negatively by the DC power supply, respectively.The potential drop between the ion beam and the neutralizer in the actual operation is simulated by that between the ion collector and the anode target in the experiment.

    Figure 9.Image of RF plasma cathode at an PF power of 200 W,Xe mass flow rate of 2.5 sccm, and a bias voltage of 50 V.

    Figure 10.Current–voltage characteristics of RF cathode of an RF input power of 100 W.

    The neutral xenon, which is employed as the propellant,is fed into the discharge chamber via the mass flow controller for cathode ignition.The induction coil is connected to the RF generation through the impedance matching box,in which the matching network reduces the reflected power to maximize the input power to a load.The xenon gas is coupled by the RF electric field, which is generated in the circumferential direction inside the coil,and finally,plasma is produced.The skin depth is low due to the high frequency (13.56 MHz).

    When the plasma plume supplies electrons to the anode target, the extracted electron current rapidly appears in the circuit.Figure 9 shows the image of the device under operation.As shown in the figure, the extracted electron current has a maximum value near the outlet of the discharge chamber.

    4.2.I–V characteristics

    Figure 10 shows the current–voltage characteristics of the RF plasma cathode when the RF input power is 100 W.After the ignition of the ICP, a current of approximately 0.16 A is generated with no bias voltage.The extracted current increases quickly when the bias voltage varies from 0 to 50 V and is maintained nearly constant beyond 85 V.As the voltage is increased, the extracted electrons are accelerated to ionize the Xe atoms, and a plasma plume is formed betweenthe cathode and anode target.The electron current is built after the target absorbs electrons from the plasma plume.The upper limits of the extracted electron current are observed at different Xe mass flow rates.

    Table 3.Slopes of the fitted line in the bias voltage range of 30–85 V.

    Figure 10 reveals that the extracted current increases with the variation in the mass flow rate from 2.0 to 4.0 sccm (Xe,2.7660–5.5320 sccm) at a constant RF power.However, the variable quantity of the extracted current ΔI decreases at the same mass flow rate.As shown in figure 10, the whole extracted currents have an almost linear behavior in the bias voltage range of 30–85 V.The values of the slopes are acquired via a linear fit.Table 3 presents the results.In addition,as the bias voltage increases further above 85 V,the extracted current remains nearly a constant.

    As listed in table 3, the slopes are similar in the mass flow rate range of 3.4575–4.8405 sccm.This means that the growth rates of every extracted current are almost consistent.The minimum and maximum slopes correspond to the minimum and maximum mass flow rates of Xe, respectively.It is inferred that a small amount of gas is not conducive for increasing the increment of extracted current, which albeit increases rapidly after the gas flow rate reaches a certain value.

    4.3.RF power and mass flow rate

    Figure 11 shows the electron current extracted by the target as a function of the RF input power and mass flow rate (target voltage = 50 V constant).The cathode achieved approximately 0.18 A with 50 W of RF input power, and over 1 A with 260 W,2.0 sccm(Xe,2.7660 sccm).The stability of the cathode is tested for over an hour period.The result shows that after an initial automatic adjustment of matching network for some minutes,the electron current is constant without any change in the matching.

    Figure 11 reveals that at constant bias voltage, by increasing the RF power from 50 to 300 W, the extracted current increases.However, the variable quantity of the extracted current ΔI decreases at a bias voltage of 50 V.As shown in figure 11, the whole extracted currents have a largely linear behavior in the RF power range of 80–260 W.The value of slope is approximately 0.003 89 via a linear fit.In addition,as the RF power increases further above 260 W,the extracted current is saturated.

    5.Conclusion

    Numerical and experimental parametric studies were conducted to investigate the effects of the gas flow rate, RF power and bias voltage on the plasma generation and electron extraction performance of an RF plasma cathode device.

    Figure 11.Current–power characteristics of RF cathode at a bias voltage of 50 V.

    A comparison of figures 7 and 10 shows that the extracted current is approximately 0.34 A under the condition of 5 sccm Xe,100 V bias voltage and 100 W RF input power in the numerical simulation.In addition, in the experimental result, limited by the mass flow controller, the extracted currents are 0.43 A and 0.45 A at approximately 4.84 sccm Xe and 5.53 sccm Xe,respectively,with the other conditions kept the same.It can be considered approximately 0.44 A at 5 sccm in the experimental measurement.Therefore, compared with the experiment, the relative error in the simulation is approximately 22.7%.This is due to the low plasma density obtained in the simulation software during the low-pressure calculation.The presented experimental and simulated results are approximately consistent.

    As shown in figure 11, an extracted current of 1.03 A from the cathode is achieved using 2.766 sccm xenon,270 W of RF input power at 13.56 MHz, and 50 V DC bias on the anode target.Increasing gas mass flow rate was conducive for increasing the extracted current.For the bias voltage and RF power,similar results were obtained.Moreover,the extracted current seems to have a limit value.

    An RF plasma cathode is feasible as an electron source for electric propulsion devices.The optimization design of the cathode will be part of future work for further performance enhancement, as a replacement for the hollow cathodes.

    Acknowledgments

    This work is supported by Joint Fund for Equipment Preresearch and Aerospace Science and Technology (No.6141B061203).

    ORCID iDs

    猜你喜歡
    平陽
    An experimental study on the degradation of the C12A7 hollow cathode
    Early experimental investigation of the C12A7 hollow cathode fed on iodine
    虎園觀虎
    Masterpieces 佳作
    邊疆文學(2020年5期)2020-05-22 02:58:58
    泥火交融見初心 入窯淬煉現(xiàn)珍品
    科學之友(2020年4期)2020-05-19 15:13:02
    馳騁西北的平陽商幫
    西夏研究(2020年1期)2020-04-01 11:54:40
    雷平陽:行走在故鄉(xiāng)云南
    中國三峽(2017年9期)2017-12-19 13:27:43
    平陽木版雕印戲曲拂塵紙一幅一戲四圖年畫考釋
    中華戲曲(2016年2期)2016-01-22 08:18:55
    平陽霉素注射及手術(shù)切除綜合治療口腔頜面部海綿狀血管瘤
    平陽話“垕”的詞性、功能及來源
    一个人免费在线观看的高清视频| 香蕉国产在线看| 中文资源天堂在线| 国产精品,欧美在线| 色尼玛亚洲综合影院| 麻豆国产97在线/欧美| 国产v大片淫在线免费观看| 国产亚洲欧美在线一区二区| 巨乳人妻的诱惑在线观看| 亚洲午夜精品一区,二区,三区| 18禁黄网站禁片午夜丰满| 国产精品爽爽va在线观看网站| 久久精品国产亚洲av香蕉五月| 中文亚洲av片在线观看爽| 一级毛片高清免费大全| 99精品久久久久人妻精品| 亚洲五月婷婷丁香| 少妇裸体淫交视频免费看高清| 黑人巨大精品欧美一区二区mp4| 亚洲欧洲精品一区二区精品久久久| 亚洲av免费在线观看| 欧美激情久久久久久爽电影| 麻豆成人午夜福利视频| www.999成人在线观看| 国产精品1区2区在线观看.| 精品午夜福利视频在线观看一区| 噜噜噜噜噜久久久久久91| 国产又色又爽无遮挡免费看| 亚洲中文字幕一区二区三区有码在线看 | 亚洲熟妇中文字幕五十中出| 日韩欧美一区二区三区在线观看| 天天一区二区日本电影三级| 中文在线观看免费www的网站| 久久午夜亚洲精品久久| 国产成人精品无人区| 亚洲熟妇熟女久久| 夜夜看夜夜爽夜夜摸| 国产成+人综合+亚洲专区| 成人鲁丝片一二三区免费| 嫩草影院入口| 久久中文看片网| 国产一区二区三区视频了| 人人妻人人看人人澡| 国产久久久一区二区三区| 亚洲成a人片在线一区二区| 亚洲aⅴ乱码一区二区在线播放| 99在线人妻在线中文字幕| 日本与韩国留学比较| 国产av麻豆久久久久久久| 欧美成狂野欧美在线观看| 午夜福利免费观看在线| 香蕉久久夜色| 欧美日韩亚洲国产一区二区在线观看| 黄频高清免费视频| 亚洲av免费在线观看| 国产精品久久久av美女十八| 少妇熟女aⅴ在线视频| 在线看三级毛片| 叶爱在线成人免费视频播放| 岛国视频午夜一区免费看| 国语自产精品视频在线第100页| 午夜精品在线福利| 久久久水蜜桃国产精品网| 日韩欧美国产一区二区入口| 国产乱人伦免费视频| 中文字幕av在线有码专区| av片东京热男人的天堂| 久久久精品欧美日韩精品| 一a级毛片在线观看| 国产毛片a区久久久久| 亚洲在线自拍视频| 操出白浆在线播放| 国产成人影院久久av| 中文字幕熟女人妻在线| 日本黄色视频三级网站网址| 人妻久久中文字幕网| 又大又爽又粗| 免费搜索国产男女视频| 久久久国产成人免费| 亚洲国产色片| 精品欧美国产一区二区三| а√天堂www在线а√下载| 国产精品,欧美在线| 亚洲电影在线观看av| 在线a可以看的网站| 网址你懂的国产日韩在线| 国产亚洲精品久久久久久毛片| 大型黄色视频在线免费观看| 高清在线国产一区| 国产精品久久久久久人妻精品电影| 亚洲性夜色夜夜综合| 久久精品aⅴ一区二区三区四区| 亚洲av片天天在线观看| 国产精品永久免费网站| 香蕉国产在线看| 香蕉丝袜av| 99精品欧美一区二区三区四区| 久久欧美精品欧美久久欧美| av福利片在线观看| 久久热在线av| 亚洲五月婷婷丁香| 久久久久亚洲av毛片大全| 网址你懂的国产日韩在线| 91在线观看av| 99热精品在线国产| 国产精品一区二区免费欧美| 日日干狠狠操夜夜爽| 免费看a级黄色片| 熟女电影av网| 在线观看舔阴道视频| 一级毛片高清免费大全| 免费观看精品视频网站| 18禁黄网站禁片午夜丰满| 制服丝袜大香蕉在线| 亚洲av成人av| 国产三级在线视频| 色哟哟哟哟哟哟| 老熟妇乱子伦视频在线观看| 男人的好看免费观看在线视频| 国产午夜精品久久久久久| 亚洲美女黄片视频| 国产精品久久久久久久电影 | 又大又爽又粗| 国产主播在线观看一区二区| 国产成人系列免费观看| 欧美精品啪啪一区二区三区| 脱女人内裤的视频| 一本久久中文字幕| 成在线人永久免费视频| 欧美性猛交╳xxx乱大交人| 亚洲一区二区三区色噜噜| 日日夜夜操网爽| 国产不卡一卡二| 91麻豆精品激情在线观看国产| 深夜精品福利| 国产伦精品一区二区三区视频9 | 三级毛片av免费| 国产三级在线视频| 神马国产精品三级电影在线观看| 亚洲av成人不卡在线观看播放网| 五月伊人婷婷丁香| 日韩人妻高清精品专区| 波多野结衣巨乳人妻| 日韩大尺度精品在线看网址| 两个人看的免费小视频| 桃红色精品国产亚洲av| 免费大片18禁| 国产一区二区激情短视频| 国产免费av片在线观看野外av| 久久久久亚洲av毛片大全| 激情在线观看视频在线高清| 亚洲自偷自拍图片 自拍| 国产亚洲欧美在线一区二区| 久久精品国产亚洲av香蕉五月| 亚洲欧美日韩无卡精品| 两个人的视频大全免费| 亚洲一区二区三区色噜噜| 男女那种视频在线观看| 两个人视频免费观看高清| 19禁男女啪啪无遮挡网站| 亚洲乱码一区二区免费版| 国产精品亚洲一级av第二区| 国产成+人综合+亚洲专区| 国产1区2区3区精品| 婷婷丁香在线五月| 日韩高清综合在线| 午夜福利免费观看在线| 九色成人免费人妻av| 婷婷六月久久综合丁香| 少妇熟女aⅴ在线视频| 亚洲无线观看免费| 欧美日韩亚洲国产一区二区在线观看| 欧美xxxx黑人xx丫x性爽| 一区福利在线观看| 男女午夜视频在线观看| 国产乱人视频| 又黄又粗又硬又大视频| 亚洲 欧美一区二区三区| 特大巨黑吊av在线直播| 国产精品久久电影中文字幕| 女警被强在线播放| 亚洲aⅴ乱码一区二区在线播放| 性色av乱码一区二区三区2| 国产av不卡久久| 色综合亚洲欧美另类图片| 草草在线视频免费看| 999精品在线视频| www日本在线高清视频| 亚洲成人久久爱视频| 嫩草影院入口| 精品久久蜜臀av无| 久久草成人影院| 国产69精品久久久久777片 | 免费在线观看日本一区| 精品熟女少妇八av免费久了| 一个人免费在线观看电影 | 美女被艹到高潮喷水动态| 亚洲男人的天堂狠狠| 亚洲国产中文字幕在线视频| 高清在线国产一区| 久久伊人香网站| 日韩大尺度精品在线看网址| 免费在线观看日本一区| 精品一区二区三区av网在线观看| 国内精品久久久久精免费| 99热这里只有精品一区 | 国产成人啪精品午夜网站| 99热这里只有精品一区 | 看片在线看免费视频| 日本在线视频免费播放| 国产成人啪精品午夜网站| 99热这里只有是精品50| 九九久久精品国产亚洲av麻豆 | 18禁黄网站禁片免费观看直播| 国产成人欧美在线观看| 无限看片的www在线观看| 一本精品99久久精品77| 三级毛片av免费| 美女 人体艺术 gogo| 男插女下体视频免费在线播放| 国产精品久久视频播放| а√天堂www在线а√下载| 97超视频在线观看视频| 色噜噜av男人的天堂激情| 男女视频在线观看网站免费| 中文字幕久久专区| 欧美日本亚洲视频在线播放| 午夜免费激情av| 少妇的逼水好多| 观看美女的网站| 国产成人精品无人区| 午夜福利成人在线免费观看| 欧美一区二区精品小视频在线| 午夜亚洲福利在线播放| 精品午夜福利视频在线观看一区| 人妻夜夜爽99麻豆av| 亚洲狠狠婷婷综合久久图片| 高潮久久久久久久久久久不卡| 婷婷精品国产亚洲av在线| 99在线人妻在线中文字幕| 国产成人欧美在线观看| 又粗又爽又猛毛片免费看| 亚洲aⅴ乱码一区二区在线播放| 精品一区二区三区av网在线观看| 亚洲专区中文字幕在线| 夜夜爽天天搞| 亚洲国产日韩欧美精品在线观看 | 亚洲九九香蕉| 一级毛片女人18水好多| 国产高清激情床上av| 国产精品av久久久久免费| 色视频www国产| 欧美日韩精品网址| 嫩草影院精品99| 日韩欧美国产在线观看| 亚洲成人中文字幕在线播放| 久99久视频精品免费| 国语自产精品视频在线第100页| 长腿黑丝高跟| 成人18禁在线播放| 亚洲欧美日韩无卡精品| 亚洲成人久久爱视频| 曰老女人黄片| 99国产综合亚洲精品| 国产综合懂色| 亚洲国产看品久久| 首页视频小说图片口味搜索| 亚洲 国产 在线| 最近最新免费中文字幕在线| 无限看片的www在线观看| 一二三四社区在线视频社区8| 久久草成人影院| 身体一侧抽搐| 麻豆av在线久日| netflix在线观看网站| 国产欧美日韩一区二区三| e午夜精品久久久久久久| 国产精品久久久久久久电影 | 两人在一起打扑克的视频| 亚洲熟女毛片儿| 一夜夜www| 99精品在免费线老司机午夜| 在线免费观看的www视频| 性欧美人与动物交配| 欧美日韩黄片免| 99久久无色码亚洲精品果冻| 好看av亚洲va欧美ⅴa在| 日本精品一区二区三区蜜桃| 日本在线视频免费播放| 一级a爱片免费观看的视频| 午夜精品久久久久久毛片777| 老司机在亚洲福利影院| 我要搜黄色片| 中文资源天堂在线| 国产黄a三级三级三级人| 国内久久婷婷六月综合欲色啪| 人人妻,人人澡人人爽秒播| 十八禁人妻一区二区| 国产午夜精品论理片| 亚洲av成人av| 国产又黄又爽又无遮挡在线| 亚洲片人在线观看| 老司机午夜福利在线观看视频| 18禁观看日本| 好男人在线观看高清免费视频| 日本熟妇午夜| 国产伦人伦偷精品视频| 亚洲av成人一区二区三| 天堂动漫精品| 最新美女视频免费是黄的| 搡老岳熟女国产| 精品一区二区三区四区五区乱码| 999久久久国产精品视频| 久久久久国内视频| 搡老熟女国产l中国老女人| 国产一区二区在线av高清观看| 高清在线国产一区| 男女下面进入的视频免费午夜| 一本一本综合久久| 亚洲 欧美 日韩 在线 免费| 午夜福利在线观看吧| 国产伦一二天堂av在线观看| 亚洲中文av在线| 久久午夜综合久久蜜桃| 男女做爰动态图高潮gif福利片| 夜夜夜夜夜久久久久| 一进一出抽搐gif免费好疼| 成人无遮挡网站| 亚洲无线观看免费| 亚洲无线在线观看| 日韩大尺度精品在线看网址| 深夜精品福利| 国产 一区 欧美 日韩| 搡老岳熟女国产| 国产视频内射| 国产av在哪里看| 久9热在线精品视频| 国产精品98久久久久久宅男小说| 最近最新中文字幕大全免费视频| 欧美乱色亚洲激情| 观看美女的网站| netflix在线观看网站| 麻豆国产97在线/欧美| 国产又黄又爽又无遮挡在线| 久久中文看片网| 亚洲男人的天堂狠狠| 亚洲自拍偷在线| 午夜激情欧美在线| 日韩人妻高清精品专区| 欧美日韩福利视频一区二区| 亚洲熟女毛片儿| 九九热线精品视视频播放| 观看免费一级毛片| 丰满的人妻完整版| 美女 人体艺术 gogo| 日本成人三级电影网站| 精品久久久久久久末码| 黄色日韩在线| 啪啪无遮挡十八禁网站| 曰老女人黄片| 中亚洲国语对白在线视频| www日本在线高清视频| 叶爱在线成人免费视频播放| 搡老熟女国产l中国老女人| 国产一区二区三区在线臀色熟女| 免费观看精品视频网站| 九九久久精品国产亚洲av麻豆 | 2021天堂中文幕一二区在线观| 久99久视频精品免费| 国产高潮美女av| 午夜精品久久久久久毛片777| 身体一侧抽搐| 久99久视频精品免费| 国产高清视频在线播放一区| 欧美日韩一级在线毛片| 老熟妇仑乱视频hdxx| 在线观看日韩欧美| 日韩免费av在线播放| 99在线视频只有这里精品首页| 日韩中文字幕欧美一区二区| 亚洲五月婷婷丁香| 精品国产亚洲在线| 别揉我奶头~嗯~啊~动态视频| 不卡一级毛片| 久9热在线精品视频| 悠悠久久av| 1024手机看黄色片| 12—13女人毛片做爰片一| 亚洲国产色片| 18禁观看日本| 欧美色欧美亚洲另类二区| 国产av在哪里看| 免费观看人在逋| 欧美高清成人免费视频www| 久久久久久大精品| 日韩欧美在线乱码| 我要搜黄色片| 老熟妇乱子伦视频在线观看| 亚洲五月天丁香| h日本视频在线播放| 中亚洲国语对白在线视频| 免费av不卡在线播放| 夜夜爽天天搞| 亚洲中文字幕日韩| 国产欧美日韩一区二区精品| 日本免费一区二区三区高清不卡| 亚洲人与动物交配视频| 婷婷精品国产亚洲av在线| 国产97色在线日韩免费| 色噜噜av男人的天堂激情| 国产激情久久老熟女| 男女下面进入的视频免费午夜| 亚洲国产中文字幕在线视频| 免费av毛片视频| 午夜久久久久精精品| 亚洲av美国av| 国产精品电影一区二区三区| 舔av片在线| 国产精品99久久久久久久久| 欧美成人性av电影在线观看| 欧美日韩精品网址| 美女cb高潮喷水在线观看 | 欧美丝袜亚洲另类 | 国产日本99.免费观看| 村上凉子中文字幕在线| 色综合欧美亚洲国产小说| 国产真人三级小视频在线观看| 成人三级做爰电影| 久久99热这里只有精品18| 国产av一区在线观看免费| 国产精品久久久人人做人人爽| 国产aⅴ精品一区二区三区波| 免费电影在线观看免费观看| 人妻丰满熟妇av一区二区三区| 国产精品香港三级国产av潘金莲| 亚洲人成伊人成综合网2020| 亚洲专区中文字幕在线| 免费在线观看亚洲国产| 少妇人妻一区二区三区视频| 熟女少妇亚洲综合色aaa.| 亚洲成人中文字幕在线播放| 国产日本99.免费观看| 免费在线观看影片大全网站| 可以在线观看的亚洲视频| 在线观看免费视频日本深夜| 亚洲av熟女| 超碰成人久久| 嫩草影院入口| 国产午夜福利久久久久久| 精品熟女少妇八av免费久了| 人人妻,人人澡人人爽秒播| 久久精品国产99精品国产亚洲性色| 亚洲国产精品sss在线观看| 成人亚洲精品av一区二区| 午夜a级毛片| 久久久久精品国产欧美久久久| 色吧在线观看| 国产激情欧美一区二区| 国产淫片久久久久久久久| 麻豆一二三区av精品| 最近最新中文字幕免费大全7| 午夜精品国产一区二区电影 | 97在线视频观看| 十八禁国产超污无遮挡网站| 亚洲欧美成人综合另类久久久 | 亚洲精品日韩av片在线观看| 欧美xxxx性猛交bbbb| 亚洲精品aⅴ在线观看| 97热精品久久久久久| 日韩,欧美,国产一区二区三区 | 男女啪啪激烈高潮av片| 亚洲精品456在线播放app| 亚洲欧洲国产日韩| 国产一区有黄有色的免费视频 | 毛片一级片免费看久久久久| 人妻夜夜爽99麻豆av| 国产视频内射| 菩萨蛮人人尽说江南好唐韦庄 | 久久99热这里只有精品18| 成人综合一区亚洲| or卡值多少钱| 人人妻人人澡人人爽人人夜夜 | 色视频www国产| 婷婷六月久久综合丁香| 欧美97在线视频| 欧美一区二区亚洲| 欧美最新免费一区二区三区| 国产老妇伦熟女老妇高清| 久久欧美精品欧美久久欧美| 中文字幕亚洲精品专区| 国产精品一区二区三区四区久久| 亚洲欧美成人精品一区二区| 91久久精品国产一区二区成人| 全区人妻精品视频| 国产精品国产三级国产专区5o | 99热这里只有是精品50| 18禁在线播放成人免费| 国产高清不卡午夜福利| 热99re8久久精品国产| 国产一区二区三区av在线| 直男gayav资源| 高清av免费在线| 麻豆成人av视频| 亚洲国产最新在线播放| 国产女主播在线喷水免费视频网站 | 99久久无色码亚洲精品果冻| 中文天堂在线官网| 成人亚洲精品av一区二区| 亚洲18禁久久av| 91精品一卡2卡3卡4卡| 亚洲aⅴ乱码一区二区在线播放| 国产成人精品久久久久久| 国产69精品久久久久777片| 日韩制服骚丝袜av| 免费一级毛片在线播放高清视频| 最近中文字幕高清免费大全6| 精品国产一区二区三区久久久樱花 | 成人一区二区视频在线观看| 综合色丁香网| 精品国产露脸久久av麻豆 | 亚洲精品乱码久久久v下载方式| 国产精品一区二区三区四区久久| 国产伦在线观看视频一区| 欧美xxxx黑人xx丫x性爽| 欧美一区二区精品小视频在线| 晚上一个人看的免费电影| 亚洲av福利一区| 一个人看视频在线观看www免费| 美女cb高潮喷水在线观看| 久久久午夜欧美精品| 内地一区二区视频在线| 又粗又硬又长又爽又黄的视频| 国产伦一二天堂av在线观看| 美女cb高潮喷水在线观看| av播播在线观看一区| 99久久人妻综合| 91久久精品国产一区二区三区| 国产精品1区2区在线观看.| 亚洲欧美中文字幕日韩二区| 免费观看a级毛片全部| 国产免费又黄又爽又色| 免费观看性生交大片5| 全区人妻精品视频| 欧美3d第一页| 免费不卡的大黄色大毛片视频在线观看 | 国产毛片a区久久久久| 色尼玛亚洲综合影院| 九色成人免费人妻av| 精品酒店卫生间| 久久久色成人| 亚洲成人中文字幕在线播放| 啦啦啦韩国在线观看视频| 超碰av人人做人人爽久久| 午夜福利在线观看免费完整高清在| 人妻系列 视频| 国产精品久久久久久精品电影| 国产成人午夜福利电影在线观看| 嫩草影院入口| 汤姆久久久久久久影院中文字幕 | 人人妻人人看人人澡| 欧美成人一区二区免费高清观看| 亚洲精品,欧美精品| 国产高清三级在线| 免费观看人在逋| 国产精品一区二区性色av| 午夜精品国产一区二区电影 | 欧美又色又爽又黄视频| 国产综合懂色| 日本爱情动作片www.在线观看| 亚洲精品aⅴ在线观看| 亚洲人成网站在线播| 国产精品一及| 亚洲成av人片在线播放无| 日本一二三区视频观看| 亚洲精品成人久久久久久| 成人综合一区亚洲| 毛片一级片免费看久久久久| 亚洲精品自拍成人| 亚洲国产精品国产精品| 免费看美女性在线毛片视频| av.在线天堂| 特大巨黑吊av在线直播| 床上黄色一级片| 日韩,欧美,国产一区二区三区 | 欧美bdsm另类| 99久久精品一区二区三区| 午夜福利在线观看吧| 纵有疾风起免费观看全集完整版 | 99久国产av精品| 亚洲熟妇中文字幕五十中出| 成年版毛片免费区| 国产乱来视频区| 国产男人的电影天堂91| 黄色配什么色好看| 99热全是精品| 色吧在线观看| 久久久久国产网址| 欧美xxxx性猛交bbbb| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 久久人人爽人人片av| 亚州av有码| 日韩,欧美,国产一区二区三区 | 国内揄拍国产精品人妻在线| 国产亚洲精品av在线| 国产精品嫩草影院av在线观看| 看非洲黑人一级黄片| 亚洲成av人片在线播放无| 好男人在线观看高清免费视频| 国产精品av视频在线免费观看| 亚洲欧美成人精品一区二区| 老司机影院成人| 亚洲欧洲国产日韩| 麻豆成人av视频| 国产美女午夜福利| av免费在线看不卡|