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

    Growth mechanism and characteristics of electron drift instability in Hall thruster with different propellant types

    2024-01-25 07:28:52LongChen陳龍ZiChenKan闞子晨WeiFuGao高維富PingDuan段萍JunYuChen陳俊宇CongQiTan檀聰琦andZuoJunCui崔作君
    Chinese Physics B 2024年1期
    關(guān)鍵詞:陳龍高維

    Long Chen(陳龍), Zi-Chen Kan(闞子晨), Wei-Fu Gao(高維富), Ping Duan(段萍),Jun-Yu Chen(陳俊宇), Cong-Qi Tan(檀聰琦), and Zuo-Jun Cui(崔作君)

    School of Science,Dalian Maritime University,Dalian 116026,China

    Keywords: Hall thruster,electron drift instability,axial electron mobility,particle-in-cell simulation

    1.Introduction

    The Hall thruster is a well-developed annular electromagnetic propulsion device that is widely used because it meets the needs of space propulsion missions such as attitude control,orbit transfer, and inclination correction of space vehicles.[1–3]The Hall thruster forms a radially distributed magnetic field through coils arranged outside the inner wall and outer wall of the discharge channel, and the potential drop between the anode and cathode forms an axially distributed electric field.Under the orthogonal electromagnetic field,the electrons emitted from the cathode perform Hall drift motion along the azimuthal direction of the channel and collide with the propellant gas atoms injected into the channel from the anode to produce plasmas.The ions are then accelerated by the axial electric field to generate thrust.Since the ion mass is significantly greater than the electron mass and the ion Larmor radius exceeds the size of the discharge channel, ions are commonly regarded as non-magnetized.The electrons experience an azimuthal drift velocity on the order of 106m/s due to the action of the orthogonal electromagnetic field.This velocity difference between the electrons and ions in the azimuthal direction induces a significant plasma instability in the discharge channel,which is known as electron drift instability(EDI).Tsikataet al.[4]measured the oscillation of electron density in the discharge channel of a Hall thruster using a collective light scattering method.The experiment indicated that EDI exhibits high-frequency oscillation in the order of~MHz in both the azimuthal direction and radial directions.

    The evolution of high-frequency short-wave oscillations of EDI can be effectively observed by particle-in-cell (PIC)simulations.The majority of studies have used PIC simulation that contains azimuthal direction to simulate these oscillations.These models include one-dimensional (1D) azimuthal direction,[5–7]two-dimensional (2D) axial-azimuthal directions,[8–12]2D radial-azimuthal directions,[13–16]and three-dimensional (3D) full-channel[17,18]numerical simulations.Ducrocqet al.[19]investigated the development of EDI and suggested that when the electric drift velocity is smaller than the electron thermal velocity,the EDI fluctuations satisfy an ion-sound wave dispersion relation and the saturation of the EDI may be related to ion Landau damping.Laflueret al.[7,20]proposed the corresponding dispersion relation for EDI by establishing a 1D azimuthal direction PIC simulation.The study by Coche and Garrigues[21]found that reducing plasma density can effectively suppress EDI.Taccognaet al.[9]studied the development process of EDI by building two different 2D models of the channel and showed that EDI triggers off both radial and azimuthal plasma oscillations and that the nonlinear development in the phase exhibits ion-acoustic wave properties.Senguptaet al.[13]studied the nonlinear growth phase of EDI by full-scale 2D simulation in the radial direction and azimuthal direction.They showed that the wavelength of EDI increases with nonlinear growth phase and that the saturation of EDI is strongly correlated with the motion of the ions in the azimuthal direction.It is clear that EDI shows different characteristics in each growth stage due to different excitation mechanisms.Tavassoliet al.[22]showed that in numerical simulation studies of EDI,numerical noise often significantly affects the simulation results and obscures the actual growth mechanisms at each stage.Asadiet al.[23]showed that increasing the number of macro-particles in the model can eliminate the numerical noise to some extent.

    The EDI is considered to be one of the main reasons for anomalous axial electron transport.[24–30]Lazurenkoet al.[31]demonstrated experimentally that the EDI significantly increases the axial electron mobility of the thruster.Laflueret al.[5]established a modified model of the electron axial mobility under the influence of EDI through theoretical derivation,and the modified electron mobility is in the same order of magnitude as the actual value.Lafleur and Chabert[11]developed a 2D model in the axial-azimuthal direction and demonstrated that anomalous electron transport can be explained by an instability-enhanced friction force between electrons and ions.

    Xenon gas is widely chosen as the propellant in traditional Hall thrusters, owing to its low ionization energy and chemical stability.In recent years,other types of propellants,such as argon and krypton, have also been proposed as propellants for Hall thrusters to reduce costs.[32,32–35]The EDI is expected to exhibit different characteristics under different types of propellants,thereby affecting the discharge process of Hall thrusters.Asadiet al.[23]investigated the axial electron mobility variation under the effect of EDI for different propellants corresponding to ion weight.However, the study used the same collision cross-section, which may not fully reflect the true influences of different propellants on EDI characteristics.

    From the content mentioned above, it can be observed that the precise growth mechanisms of EDI at different stages(linear and nonlinear) remain undetermined, and the detailed characteristics of EDI in Hall thrusters with different propellant types have not been comprehensively studied.Additionally, the effects of electromagnetic field strength and propellant types on the axial electron transport caused by EDI have not been thoroughly investigated.Based on the typical size and operating parameters of the Hall thruster, a 1D physical model of the azimuthal discharge channel is established in this work.In this model, the number of particles is increased to reduce the influence of numerical noise on the simulation results, and the actual differential scattering cross-sections for different types of propellants are considered.The improved 1D azimuthal PIC model is employed to investigate the development process of EDI and the influences of different propellants on EDI and electron axial migration rate.

    The rest of the paper is organized as follows.The physical model is established and the basic equations are introduced in Section 2.In Section 3,the numerical simulation results are described and analyzed.Finally,the results are summarized in Section 4.

    2.Physical model and basic equations

    The discharge channel of the thruster is an annular space formed by an inner ceramic sleeve and outer ceramic sleeve.The wavelength of the EDI is within the millimeter range,which is significantly smaller than the dimension of the Hall thruster discharge channel.Hence, the influence of curvature in the azimuthal direction in this model is not considered.Figure 1 shows a typical 3D structure of the Hall thruster and a schematic of the simulation region for this study.The thruster discharge channel is expanded in the azimuthal direction,and the profile is taken in the channel axial direction to establish a 2D coordinate phase space and a 3D velocity phase space.The lengths of the axial boundary and azimuthal boundary of the respective discharge channels in the simulation area are represented byLzandLyrespectively.A radial magnetic fieldB0perpendicular to the axial-azimuthal plane is applied to the simulation area,and a uniform and constant axial electric fieldE0is applied in thezdirection.The electric field strength in the simulation area is set to 2×104V/m, and the magnetic field strength is assumed to be 200 Gs(1 Gs=10?4T).[37]

    The particle positions are tracked in the axial direction and azimuthal direction,and the axial electric field component is provided by the applied electric fieldE0,while the azimuthal electric field component is obtained by solving Poisson’s equation in the azimuthal direction.It is worth noting that the 1D model established in this study solves the Poisson equation only in the azimuthal direction.The particles in this simulation exist in a 2D positional coordinate space and a 3D velocity space.

    The electrons are magnetized and their collisions with xenon atoms are considered in the simulation,including single ionization, excitation, and elastic scattering reactions, which are processed by the Monte Carlo collision(MCC)algorithm.Only energy loss and scattering are considered in the treatment of ionization and excitation collisions,and neither newly generated electrons nor ions are added to the simulation area.Additionally,because of the significant difference in mass between ions and electrons,this model does not take into account the influence of the magnetic field on ions.The influence of the axial electric field on ion motion is limited to changing the effective transit time through the acceleration region.Due to the smaller axial length of the simulation region than the actual size of the Hall thruster discharge channel, the ions are given an additional constant axial velocity,vdito simulate the ion motion near the exit of the discharge channel.

    The boundary conditions of the simulation region are set as follows.The azimuthal angular boundary uses periodic boundary conditions, while the axial boundary uses the following conditions: ions passing through the exit plane of the axial region are removed from the simulation and replaced by new ions loaded at the entrance with an initial temperature of 0.1 eV; electrons lost at the entrance plane are removed and replaced by new electrons loaded at the exit plane with an initial temperature of 10 eV.The parameters are selected from the literature.[23,37,38]Consequently,the total number of“macro particles”and the overall plasma density are kept constant.When new particles are loaded, their azimuthal positions are randomly distributed within the range ofLy.Table 1 displays the parameters of the physical quantities used in the simulation region.The values of the magnetic field intensity and the electric field intensity are chosen based on the parameters near the exit of a classical Hall thruster discharge channel.[1,34,36,39]Simultaneously, considering program runtime, computational load, and the actual values within Hall thruster discharge channels,[36]in the model a simulation region is specified to have an axial length of 1 cm and an azimuthal length of 2.75 cm.

    In the model a fixed time step is set to 5×10?12s,and a fixed grid cell length to 5×10?5m.These set values ensure that the model continues to satisfy the Courant–Friedrichs–Lewy (CFL) condition even after electron temperature has rised.[40]

    Fig.1.Typical 3D structure of Hall thruster and schematic diagram of simulation region.

    Table 1.Physical parameters used in PIC simulation.

    Furthermore, to minimize the influence of numerical noise on the simulation results, the model increases the total number of“macro particles”compared with the conventional 1D model.Specifically,1000 particles are allocated to per cell.Table 2 displays the variations in oscillation frequency, particle energy,runtime duration,and coefficient of determination(R2)in the EDI’s quasi-steady phase with different numbers of particles per cell(Nppc)in the simulation program.Here,R2is calculated from the approximate fitting to the electron density oscillations at 6μs by using the sum of sin functions curve:

    whereyirepresents the actual value being fitted, ?yirefers to the fitted value, and ˉydenotes the mean of the actual values.The coefficient of determinationR2measures the quality of the fitting of the fitted values to the true values.When fitting electron density oscillates,R2can also reflect the smoothness of the oscillation curve of electron density.A smoother curve indicates less numerical noise in the program.

    From Table 2, it can be observed that whenNppc≤500,the oscillation frequency of EDI, and particle energy show noticeable and irregular variations due to the value ofNppcchanging.However,whenNppc>500,all parameters listed in Table 2 exhibit insignificant differences when particle counts change.It can also be seen from Table 2 that with the increase ofNppc, theR-squared value after fitting the electron density oscillations also increases.This implies that withNppc>500,it is sufficient to reduce the interference of numerical noise to the results of simulation program.Considering the program runtime,the program is set withNppc=1000.

    Table 2.Frequencies,particle energy,coefficients of determination,and runtimes in EDI’s quasi-steady phase with different numbers of particles per cell.

    3.Simulation results and discussion

    3.1.Driving force and development process of EDI

    The main features of EDI are high frequency reaching to~MHz,large-amplitude oscillations of the azimuthal electric field and the plasma density in the channel.The simulation uses a time step of 2.5×10?10s and observes oscillation in the azimuthal electric field and electron density at a frequency of 5 MHz.Figures 2(a) and 2(b) show the temporal evolution of azimuthal electric field and electron density, revealing that the electron density oscillation lags behind azimuthal electric field oscillation, which means that the plasma density oscillation is triggered off by the electric field.As EDI develops, the maximum amplitude of azimuthal electric field oscillation can reach to 8 times the magnitude of the applied axial electric field, while the electron density oscillation can reach to 2 times the plasma density.Figure 2(c) shows the oscillation of azimuthal electric field in the first 0.375μs.Figure 2(d)shows the temporal evolution of the amplitude of azimuthal electric field oscillation in 0–3 μs.Comparison of Fig.2(c) with Fig.2(d) shows that the oscillation of the azimuthal electric field initiates at 0.1 μs, and EDI undergoes a linear growth phase for the first 0.375 μs.The amplitude growth of azimuthal electric field oscillation exhibits an exponential growth characteristic during the linear growth period.Theγin Fig.2(d)shows the amplitude linear growth rate normalized by the ion plasma frequency.The azimuthal electric field distribution at each time point is subjected to a fast Fourier transform (FFT), and the amplitude corresponding to the maximum frequency is selected as the characteristic amplitude at the corresponding time, followed by a linear fit to calculate the amplitude growth rate.

    Fig.2.Temporal evolution of azimuthal electric field, electron density, and the temporal evolution of the amplitude of azimuthal electric field oscillations: (a) electric field in azimuthal direction, (b) electron density, (c) azimuthal electric field evolution during the linear growth phase,(d)amplitude of the azimuthal electric field.

    Fig.3.Temporal evolution of amplitude of azimuthal electric field oscillations and the corresponding frequency spectrum analysis: (a) amplitude evolution(I:linear growth phase;II:nonlinear growth phase;III:quasi-steady-state phase), (b) frequency spectrum at nonlinear growth and quasi-steady-state stages.

    Figure 3 shows the temporal evolution of the amplitude of azimuthal electric field oscillation and the corresponding frequency spectrum analysis, here, the amplitude is normalized by the axial electric field strength.As shown in Fig.3(a),the amplitude variation of EDI follows a linear growth relationship before 0.375 μs.At the timepoint of 0.375 μs,where the amplitude reaches 2.5 times the axial electric field,the linear growth relationship is disrupted, and the growth rate slows down.At the timepoint of 0.75 μs, the amplitude of the EDI experiences a sharp drop, followed by oscillation and slow growth.After the 2.75 μs, EDI enters a high-amplitude oscillation mode,with the amplitude reaching to 7 times the strength of the axial electric field.Three clear evolution phases are shown in this process of EDI: (I) linear growth phase (0.1 μs–0.375 μs), in which the amplitude of the azimuthal electric field oscillation increases linearly; (II)nonlinear growth phase (0.375 μs–2.75 μs): EDI enters the first nonlinear growth phase(II-1 shown in Fig.3(a)),in which the amplitude of EDI continues to increase, but no longer follows the linear growth rate; after 0.75 μs, EDI enters the second phase of nonlinear growth (II-2 shown in Fig.3(a)),where the oscillation amplitude of EDI increases slowly and exhibits pronounced oscillation;(III)quasi-steady-state phase(after 2.75 μs): EDI reaches saturation finally, where the oscillation frequency no longer changes and the amplitude exhibits irregular oscillation.The spectrum analysis of the nonlinear growth phase and quasi-steady-state period of the EDI are shown in Fig.3(b).As can be seen from this figure, the EDI has a major frequency of 5 MHz in the nonlinear growth phase,and the major frequency of EDI decreases to 4.9 MHz after entering the quasi-steady state.

    Figure 4 shows the temporal evolution of the average particle energy and its azimuthal energy within the simulation region.As shown in Fig.4(a),the average energy of ions rapidly increases, and reaching 95 eV at the 0.75μs.This is primarily due to the acceleration of ions in the axial direction by the applied electric field.Subsequently, the ion energy undergoes oscillation until the EDI reaches a quasi-steady state.The electron energy first increases briefly,then decreases,and eventually stabilizes with oscillations in a range from 30 eV to 45 eV.From Fig.4(b),it can be observed that the azimuthal ion energy starts to increase from 0.1 eV after the onset of EDI.It exhibits different growth patterns during different stages of EDI development until it stabilizes at 6 eV after the EDI enters a quasi-steady state.The variation of electron temperature in the azimuthal direction follows a similar trend to the overall energy change.This means that due to EDI,there occurs a significant ion heating mechanism in the azimuthal direction.It is worth noticing that the variation of azimuthal ion energy during the linear growth and the first nonlinear growth phase of EDI shows a similar trend to the amplitude variation shown in Fig.3(a),which exhibits a short and rapid initial increase, followed by a brief decrease and then a slow growth.These pieces of evidence indicate a close relationship between the azimuthal ion heating mechanism and the evolution of EDI.To investigate the ion heating mechanisms of EDI and their relationship with the development of instability,four time points,denoted ast1=0.1μs,t2=0.375μs,t3=0.75μs,andt4=2.75μs,to represent different stages of EDI development,are selected.The azimuthal particle density oscillations,velocity phase space evolution,and changes in velocity distribution fitting profiles at each of these time points are shown in Figs.5–7.

    Figure 5 shows the azimuthal distributions of electron density and ion density at different time points.It can be observed that the electron density distribution oscillates before the ion density distribution.As the particle energy increases,the amplitude of plasma density oscillation also increases.Att2=0.375 μs, the amplitude of the plasma density reaches twice that of the initial density and stabilizes.This is because electrons have a higher azimuthal velocity than ions.Because of this velocity difference, the electron density oscillates before the ion density, and the oscillation of ion density is triggered off through their interactions.The interaction between electrons and ions, which results in oscillation in ion density and energy loss for electrons in the azimuthal direction, is commonly referred to as electron–ion friction.By combining Figs.4(b) and 5, it can be observed that EDI induces a deceleration of electrons and an acceleration of ions through the increase of electron–ion friction in the azimuthal direction.When electron density oscillation and ion density oscillation are almost synchroinous,the ion azimuthal energy reaches its maximum value of 6 eV, and the amplitude of ion density is slightly higher than that of electron density.

    Fig.4.Temporal evolution of (a) average particle energy and (b) its azimuthal energy in simulation region.

    Fig.5.Ion density distributions and electron density distributions at different time points.

    Figure 6 shows the evolutions of the ion phase space plot in azimuthal direction at different time points.As observed in this figure,the ion phase space gradually forms a vortex structure.Att3=0.75μs,vortices start to emerge in the ion phase space, accompanied by the formation of a striped long-tail,while the amplitude of the instability oscillation increases to a larger amplitude.Att4=2.75 μs, the majority of ion velocities are limited within the range of±5000 m/s.Combining with Fig.3(a), it can be observed that EDI exhibits irregular amplitude oscillations att3=0.75μs.This is due to the fact that aftert3=0.75μs, the EDI wave of the phase velocity is close to the ion thermal velocity.Under the influence of nonlinear Landau damping, the EDI and ions are interacted with each other, entering the nonlinear development phase.After entering the nonlinear phase, resonance occurs between the EDI wave and ions with velocities close to the wave speed.When there is a larger number of ions with velocities lower than the wave velocity, the particles gain more energy from the wave than their contribution to it.This explains the difference in the growth trend of ion azimuthal energy between the nonlinear growth stage and the linear growth stage in Fig.4(b).This phenomenon hinders the development of EDI and further heats the ions in the azimuthal direction.This is consistent with the phenomenon observed in the experiment in Ref.[41].As the ion velocity increases,distinct vortex patterns appear in the ion phase space,indicating that the ions are trapped in potential wells.Starting from this stage,the trapped ions exhibit energy oscillation as shown in Fig.4(a).

    Figure 7 shows the fitted azimuthal ion velocity distribution profiles at different time points.The dotted lines in Fig.7 represent Maxwellian distributions of the velocity distribution function (VDF) for ions with the same temperature at the corresponding time.The dotted dash lines indicate the wave speed of EDI at the corresponding time points.From Fig.7, it can be observed that the ion azimuthal velocity follows the Maxwellian distribution when EDI starts to generate.However, as EDI grows, the ion azimuthal velocity gradually deviates from the Maxwellian distribution.The actual velocity distribution profile is broader than the Maxwellian distribution,with the proportion of high-speed ions increasing significantly.This broadening effect is attributed to ion acceleration in the azimuthal direction, caused by electron–ion friction which results in the ions velocity deviating from the typical Maxwellian distribution.As shown in Figs.7(b)–7(d),the velocity distribution profiles in the azimuthal direction gradually shift towards positive values with ion energy increasing.Moreover, the portion of the velocity distribution profile exceeds the value of the dotted dash line gradually approaches zero.This is due to the fact that after entering the nonlinear growth phase, EDI extracts energy from ions with velocities exceeding the wave speed,leading it to further grow.The interaction between waves and ions reduces the proportion of velocities exceeding the wave speed in the velocity distribution.When the portion of ions with velocities exceeding the wave speed approaches zero,ions become trapped in potential wells.At this stage,the ion heating mechanism of EDI is disrupted,and the ability to gain energy through wave-particle interaction is suppressed.Consequently,the EDI enters a quasisteady state where ions no longer interact with the wave, and ion energy tends to stabilize.

    Fig.6.Evolutions of ion phase space in azimuthal direction at different time points.

    Fig.7.Fitting of ion azimuthal velocity distributions at different time points.

    3.2.Effect of propellant type on EDI

    In this study, three different types of propellants,Ar,Kr,and Xe,are selected to investigate their influences on the EDI.The variations in the mass of the propellant ions, atomic ionization energy, and scattering cross-section are considered.Only the collisions between electrons and neutral particles,including elastic, excitation, and ionization collisions, are considered under different types of propellants.

    Fig.8.Spatiotemporal distribution of angle electric fields with different propellants: (a)Xe,(b)Kr,and(c)Ar.

    Figure 8 displays the spatial and temporal distributions of the azimuthal electric field with different propellants.It can be observed from the figure that changing the propellant type does not have a significant inhibitory effect on the generation of EDI,as the oscillations of the azimuthal electric field occur nearly simultaneously in different propellants.Table 3 shows the duration, initial amplitude, and linear growth rate of the oscillatory linear growth phase of EDI with different propellants,with the linear growth rate normalized by the ion plasma frequency of Xe+.From the table,it can be observed that as the ion mass decreases, the linear growth duration of EDI gradually decreases and the initial amplitude decreases,the linear growth rate increases.This is due to the lighter ions being more susceptible to the effect of electron–ion friction,resulting in an earlier onset of EDI and a faster growth in its amplitude.The oscillation amplitude,frequency,phase velocity, and ion temperature are shown in Figs.9 and 10 to study the influences of different propellants on the EDI.

    Table 3.Durations,initial amplitudes,and growth rates in EDI’s linear growth phase with different types of propellants.

    Figure 9 shows the temporal evolutions of the amplitude of azimuthal electric field oscillation and the variations in frequency and phase velocity with different propellants.In this figure,the amplitude of oscillation is normalized by the axial electric field,and the frequency and phase velocity of the oscillation are analyzed as the system reaches a quasi-steady state.As shown in the figure, reducing the ion mass leads to faster attainment of a quasi-steady state in EDI,and it can effectively suppress the amplitude of oscillation after EDI has reached a quasi-steady state.Notably, the oscillations under Kr propellant and Ar propellant enter a large oscillation mode faster than under the Xe propellant.The frequency and phase velocity in the quasi-steady state of the EDI gradually increase as the ion mass decreases.The steady-state frequency of EDI oscillations in Ar, compared with that in Xe, increases by 30.23%,and the phase velocity increases by 30.22%.This is because ions with larger mass are not easily accelerated by the axial electric field,which leads to longer residence times of ions in the simulation region and longer EDI oscillation period, thus reducing the frequency.The lower acceleration and longer residence time of ions with heavier mass in the region also contribute to a decrease in the phase velocity of EDI oscillations.

    Figure 10 shows the temporal evolutions of the average energy and azimuthal energy of particles for different propellants.From Fig.10 it is evidently seen that as the ion mass decreases,electron temperature deceases sightly.The total energy gained by ions from the axial electric field decreases and stabilizes more quickly.Moreover, the maximum energy of ions in the azimuthal direction increases, while the energy in the electron azimuthal direction gradually decreases.The ions are heated more rapidly in the azimuthal direction, resulting in a reduction in the amplitude of ion azimuthal energy oscillation during the second nonlinear growth phase of EDI and a faster transition to the quasi-steady state.However, after entering the quasi-steady phase, there is no significant change in the azimuthal energy of ions for different propellants.This indicates that the change in propellant types does not significantly affect the effectiveness of the ion heating mechanism of EDI.Additionally, with ion mass decreasing, ions become more susceptible to acceleration from the axial electric field and electron–ion friction.The former leads to the increase of axial acceleration,causing ions to exit from the simulation region more quickly.The latter results in a higher heating rate of ions in the azimuthal direction and facilitates faster interaction between ions and EDI.The faster interaction between ion and EDI wave enhances the possibility of ion capture,thereby suppressing the amplitude of ion energy oscillation and shortening the duration of the second nonlinear growth stage of EDI.Consequently,ions reach a steady state more rapidly.

    Fig.9.(a)Temporal evolutions of amplitude of azimuthal electric field oscillation and(b)variations of frequency and phase velocity of EDI with atomic mass for different propellants.

    Fig.10.Variations of(a)ion average energy,(b)electron average energy,(c)ion azimuthal energy,and(d)electron azimuthal energy with time for different propellants.

    3.3.Effect of EDI on axial electron mobility

    EDI may be one of the main reasons for abnormal axial electron mobility.Considering the electron migration triggered by electrical drift instability, the classical electron mobility is corrected by introducing the electron momentum conservation equation as follows:[6,7,38]

    wheremandqare the electron mass and charge,neandυdeare the electron density and drift velocity,respectively,andυmis the collision frequency of electrons with neutral atoms in the simulation region,Reiis the electron–ion friction force density,andΠerepresents the electron stress tensor.

    By neglecting the electron–ion friction term in Eq.(2),the classical electron mobility due to collisions can be obtained as follows:[7,24]

    whereωceis the electron cyclotron frequency, and the classical mobility of electrons in the simulation region is 0.0021 m2/(V·s)calculated from Eq.(3),which is in the same order of magnitude as the simulation results in Ref.[34].

    By taking into consideration the influence of EDI, the axial electron migration rate, equation (2) can be expressed as[7,24]

    When considering the influence of azimuthal electric field oscillation in EDI, the axial electron migration rate can also be expressed as follows:[7,24]

    Furthermore,in this study,the actual axial electron mobility in the simulation region is calculated by working out the average axial electron velocity in the simulation region:

    wherevjzis the the axial component of the electron velocity,Nis the number of electrons in the simulation region,andE0is the axial electric field set by the model.The ion heating mechanism is a complex process influenced by various factors, making it challenging to be analyzed quantitatively.By comparing Eq.(4)with Eq.(5),the electron–ion friction force density can be expressed as[7]

    In the study the parameterReiis used as a measure of the response of the EDI to the heating capacity of ions in the azimuthal direction.It is described in Ref.[34]as follows:

    whereCsis the ion sound velocity.In this simulation,the initial velocity of ions is very low, and uniformly distributed in the azimuthal direction and the radial direction.Thus,the expression ofReican be rewritten as[9]

    Therefore, combining Eqs.(4) and (10), it can be concluded that the axial electron migration rate can also be expressed under the condition satisfied by Eq.(10)as follows:[38]

    Figure 11(a) shows the axial electron mobility in the simulation region within the 2-μs interval for different propellants after the EDI has reached the quasi-steady state.The axial electron mobility is calculated from Eqs.(5) and(7).From the figure, it can be observed that the axial electron mobility for Xe propellant in the simulated region is 9.98 m2/(V·s), which is consistent with the result reported in Ref.[23].When considering the actual differential scattering cross-section, the axial electron mobility fluctuates within a range of 9.98 m2/(V·s)–11.50 m2/(V·s),showing no significant variation with change in the type of propellant.This indicates that changes in the type of propellant will not significantly affect the axial electron mobility at the exit of the thruster discharge channel.Figure 11(b)depicts the variation in axial electron mobility when neglecting change in the differential scattering cross-section, yielding similar results to those reported by Asadiet al.[23]which demonstrates that increasing the ion mass can suppress the axial electron mobility effectively.This is attributed to the fact that neglecting the chnges in differential scattering cross-sections can distort the electron energy loss caused by collisions in the model.By combining

    Fig.11.Variations of(a)axial electron mobility and(b)electron energy with atomic mass for different propellants.

    Fig.10(b)with Fig.12,it can be observed that Xe atoms predominantly experience elastic collisions with electrons in the simulation region.But Ar atom and Kr atom exhibit a higher frequency of excitation collisions with electrons, leading to more significant electron energy losses.Referring to Eqs.(8)–(10), it becomes evident that change in electron temperature directly influences variation inRie,subsequently leading to the variation of axial electron mobility.As depicted on the vertical axis of Fig.11,the electron temperature variation resulting from change in scattering cross-sections align with the corresponding trends in axial electron mobility under respective conditions.

    Fig.12.Variations of differential scattering cross-sections with electron temperature for different propellants.

    Fig.13.Variations of(a)axial electron mobility and(b)electron–ion friction force density with electric field strength for different propellants;variations of(c)axial electron mobility and(d)electron–ion friction force density with magnetic field for different propellants.

    Figure 13 shows the variations of axial electron mobility and the corresponding electron–ion friction force density with the applied electric and magnetic field strength for different propellants.It can be observed from Fig.11 that the axial electron mobility calculated from Eqs.(5) and (7) exhibit good consistency, therefore figure 13 only displays the axial electron mobility calculated from Eq.(7).As shown in Figs.13(a)–13(c), smaller axial electric field and larger magnetic field can effectively suppress the axial electron mobility.This aligns with the conclusion drawn from the experimental findings in Ref.[42].Changes in types of propellants with external electric strength and magnetic field strength do not have a significant influence on the axial electron mobility.According tovde~E/B, the increase in electric field strength and the decrease in magnetic field strength both result in a larger drift velocity of electrons in the azimuthal direction.Moreover, the smaller magnetic field strength also results in the lager electron cyclotron drift radius.This weakens the ability of the electromagnetic field to bind the electrons,allowing the electrons to quickly pass through the simulation region.As illustrated in Figs.13(b) and 13(d), the enhancement of the electron–ion friction force density further promotes the increase in axial electron mobility.

    4.Conclusions

    This study focuses on the physical processes of the electron drift instability in the near-exhaust region of Hall thruster channel.A 1D physical model is established, and the generation,development,and saturation process of EDI,as well as its influence on particle energy and density, are numerically studied by using particle-in-cell simulation.The influences of different propellants on the EDI and axial electron mobility are also discussed.The obtained results are given below.

    During the linear growth phase, EDI exhibits a significant ion heating mechanism in the azimuthal direction through electron–ion friction,resulting in ion acceleration.In the nonlinear growth phase,EDI wave speed is similar to the ion velocity and further grows to reach a quasi-steady state through the Landau damping effect.However,the ion heating mechanism of EDI becomes ineffective due to the ion wave trapping effect caused by Landau damping.

    The choice of propellant in Hall thruster can affect the frequency and amplitude of the EDI oscillation greatly.Increasing the mass of the propellant prolongs the duration of the linear growth phase of EDI and reduces the linear growth rate of the amplitude.Comparing with Xe as the propellant, the linear growth rate of the amplitude of EDI under Ar decreases by 38.34%.Once EDI reaches a quasi-steady state,the smaller ion mass leads to higher frequency and smaller amplitude of EDI oscillation.Comparing with Xe,the oscillation frequency and phase velocity of EDI under Ar increases by 30.23%and 30.22%,respectively.

    The EDI can lead to an abnormal increase in the axial electron mobility by increasing electron–ion friction force.Variations in propellant type,with the actual differential scattering cross-section considered, are not sufficient to cause a significant change in the axial electron mobility.Additionally,smaller electric field and higher magnetic field strengths can effectively suppress the axial electron mobility.It is consistent with relevant experimental research results.

    Acknowledgements

    Project supported by the National Natural Science Foundation of China (Grant Nos.11975062 and 11605021) and the Fundamental Research Funds for the Central Universities(Grant No.3132023192).

    猜你喜歡
    陳龍高維
    情書
    一種改進(jìn)的GP-CLIQUE自適應(yīng)高維子空間聚類算法
    茶,有點(diǎn)苦
    準(zhǔn)確審題正確列式精確驗(yàn)證
    教師·下(2017年10期)2017-12-10 12:35:13
    《機(jī)械工程測試技術(shù)》教學(xué)方法初探
    基于加權(quán)自學(xué)習(xí)散列的高維數(shù)據(jù)最近鄰查詢算法
    一般非齊次非線性擴(kuò)散方程的等價(jià)變換和高維不變子空間
    賀聰、胡軼群、張釗浩、陳龍作品
    高維Kramers系統(tǒng)離出點(diǎn)的分布問題
    肉雞腸毒綜合癥防治
    国产亚洲精品一区二区www | 欧美在线一区亚洲| 69精品国产乱码久久久| 免费av中文字幕在线| 99九九在线精品视频| 欧美日韩亚洲高清精品| a级片在线免费高清观看视频| 男女下面插进去视频免费观看| 夜夜夜夜夜久久久久| 精品一区二区三区av网在线观看 | 国产一区二区三区在线臀色熟女 | 1024视频免费在线观看| 亚洲avbb在线观看| 搡老熟女国产l中国老女人| 丝袜在线中文字幕| 午夜福利,免费看| 日韩人妻精品一区2区三区| 一区二区av电影网| 9191精品国产免费久久| 伊人久久大香线蕉亚洲五| 亚洲久久久国产精品| 这个男人来自地球电影免费观看| 国产精品 欧美亚洲| www日本在线高清视频| 亚洲少妇的诱惑av| 精品久久久久久久毛片微露脸| 麻豆av在线久日| 国产黄频视频在线观看| 丰满迷人的少妇在线观看| 国产成+人综合+亚洲专区| av线在线观看网站| 午夜福利免费观看在线| 菩萨蛮人人尽说江南好唐韦庄| 啦啦啦视频在线资源免费观看| 国产成人欧美| 亚洲精品美女久久av网站| 青草久久国产| 手机成人av网站| 国产在线观看jvid| 少妇裸体淫交视频免费看高清 | 精品国产乱码久久久久久男人| 99久久国产精品久久久| 精品高清国产在线一区| 三上悠亚av全集在线观看| 精品人妻熟女毛片av久久网站| 桃红色精品国产亚洲av| 国产精品电影一区二区三区 | 老汉色av国产亚洲站长工具| 亚洲色图 男人天堂 中文字幕| 久久热在线av| 999久久久国产精品视频| 亚洲精品美女久久av网站| 亚洲熟女精品中文字幕| 成人精品一区二区免费| 一个人免费看片子| 电影成人av| 国产激情久久老熟女| 欧美激情久久久久久爽电影 | 淫妇啪啪啪对白视频| 国产亚洲一区二区精品| 男女床上黄色一级片免费看| 18禁黄网站禁片午夜丰满| 国产av国产精品国产| 国产成人欧美| 五月开心婷婷网| 视频区图区小说| 91精品国产国语对白视频| 亚洲国产成人一精品久久久| 久久人妻福利社区极品人妻图片| 午夜激情久久久久久久| 国产亚洲av高清不卡| 在线永久观看黄色视频| 久久九九热精品免费| 国产精品一区二区免费欧美| 精品国产亚洲在线| www.999成人在线观看| 制服诱惑二区| 国产精品电影一区二区三区 | 亚洲性夜色夜夜综合| 日韩一卡2卡3卡4卡2021年| 女人爽到高潮嗷嗷叫在线视频| 日本五十路高清| 老司机影院毛片| av片东京热男人的天堂| 亚洲专区字幕在线| 精品视频人人做人人爽| 国产精品成人在线| 丝袜美足系列| 久久狼人影院| 日韩人妻精品一区2区三区| 黄色丝袜av网址大全| 在线观看免费视频日本深夜| 国产亚洲欧美在线一区二区| 欧美成人免费av一区二区三区 | 悠悠久久av| 啦啦啦 在线观看视频| 免费在线观看黄色视频的| 成人精品一区二区免费| 久久精品熟女亚洲av麻豆精品| 亚洲黑人精品在线| 黑人欧美特级aaaaaa片| 男女午夜视频在线观看| 亚洲黑人精品在线| 99精品在免费线老司机午夜| 免费少妇av软件| 久久久久久久精品吃奶| 高潮久久久久久久久久久不卡| 女人被躁到高潮嗷嗷叫费观| 日韩三级视频一区二区三区| 精品午夜福利视频在线观看一区 | 男女之事视频高清在线观看| 成人18禁高潮啪啪吃奶动态图| 丰满饥渴人妻一区二区三| 捣出白浆h1v1| 老司机午夜福利在线观看视频 | 黄频高清免费视频| 久久热在线av| 无遮挡黄片免费观看| 最近最新中文字幕大全免费视频| 日本av手机在线免费观看| 好男人电影高清在线观看| 亚洲精品粉嫩美女一区| 可以免费在线观看a视频的电影网站| 免费日韩欧美在线观看| 国产免费av片在线观看野外av| 久久九九热精品免费| 真人做人爱边吃奶动态| 国产日韩欧美在线精品| 久久精品成人免费网站| 国产精品电影一区二区三区 | 肉色欧美久久久久久久蜜桃| 国产一区二区 视频在线| 老熟妇仑乱视频hdxx| 亚洲性夜色夜夜综合| 人妻久久中文字幕网| 一级黄色大片毛片| 大陆偷拍与自拍| 色综合欧美亚洲国产小说| 九色亚洲精品在线播放| 国产一区二区激情短视频| 在线av久久热| 一二三四在线观看免费中文在| 亚洲久久久国产精品| 久久这里只有精品19| 五月天丁香电影| kizo精华| 久久久精品国产亚洲av高清涩受| 欧美国产精品一级二级三级| 日韩一区二区三区影片| 男女无遮挡免费网站观看| 亚洲性夜色夜夜综合| 美女扒开内裤让男人捅视频| 男人舔女人的私密视频| 亚洲国产欧美一区二区综合| 久久国产亚洲av麻豆专区| 国产日韩欧美亚洲二区| 黄片大片在线免费观看| 99国产精品99久久久久| 巨乳人妻的诱惑在线观看| 久久久久国内视频| 午夜福利影视在线免费观看| a级片在线免费高清观看视频| 99热国产这里只有精品6| 精品久久久精品久久久| 免费在线观看完整版高清| 久久久精品国产亚洲av高清涩受| 欧美激情久久久久久爽电影 | 中文亚洲av片在线观看爽 | 国产一区二区在线观看av| 美女国产高潮福利片在线看| 亚洲久久久国产精品| 国产成人免费观看mmmm| 午夜两性在线视频| 亚洲人成77777在线视频| 久久久国产一区二区| 激情视频va一区二区三区| 国产精品一区二区在线观看99| 人妻久久中文字幕网| 国产无遮挡羞羞视频在线观看| 国产一区二区三区在线臀色熟女 | 久久精品国产99精品国产亚洲性色 | 成年人黄色毛片网站| netflix在线观看网站| 亚洲欧美色中文字幕在线| 久久精品成人免费网站| 国产成+人综合+亚洲专区| 热99久久久久精品小说推荐| 国产精品美女特级片免费视频播放器 | 精品少妇黑人巨大在线播放| 老司机福利观看| 天堂俺去俺来也www色官网| 精品人妻熟女毛片av久久网站| 国产精品亚洲av一区麻豆| 精品少妇一区二区三区视频日本电影| 亚洲人成77777在线视频| 免费av中文字幕在线| 天天躁夜夜躁狠狠躁躁| 黄片大片在线免费观看| 色老头精品视频在线观看| 亚洲成av片中文字幕在线观看| av片东京热男人的天堂| 每晚都被弄得嗷嗷叫到高潮| 老鸭窝网址在线观看| av国产精品久久久久影院| 天天躁狠狠躁夜夜躁狠狠躁| 亚洲人成电影免费在线| 黑人巨大精品欧美一区二区mp4| videosex国产| 亚洲情色 制服丝袜| 亚洲少妇的诱惑av| 最黄视频免费看| 国产精品久久久久久精品电影小说| 日韩三级视频一区二区三区| 国产高清国产精品国产三级| 亚洲欧美激情在线| 99在线人妻在线中文字幕 | 人妻久久中文字幕网| 精品国产国语对白av| 久久精品国产99精品国产亚洲性色 | 啪啪无遮挡十八禁网站| 丰满迷人的少妇在线观看| 日韩免费av在线播放| 久久人妻福利社区极品人妻图片| 精品欧美一区二区三区在线| 精品少妇一区二区三区视频日本电影| 日韩欧美三级三区| 嫁个100分男人电影在线观看| 高潮久久久久久久久久久不卡| 韩国精品一区二区三区| 亚洲精品粉嫩美女一区| 亚洲人成伊人成综合网2020| 亚洲午夜精品一区,二区,三区| 免费在线观看完整版高清| 极品人妻少妇av视频| 我要看黄色一级片免费的| 午夜福利视频在线观看免费| 啪啪无遮挡十八禁网站| 国产人伦9x9x在线观看| 变态另类成人亚洲欧美熟女 | 一边摸一边抽搐一进一小说 | 91九色精品人成在线观看| 日韩成人在线观看一区二区三区| 久久精品国产亚洲av高清一级| 人人妻人人澡人人爽人人夜夜| 在线永久观看黄色视频| 91av网站免费观看| 国产有黄有色有爽视频| 日本撒尿小便嘘嘘汇集6| 久久久久精品国产欧美久久久| 99国产精品99久久久久| 99久久99久久久精品蜜桃| 亚洲 国产 在线| 国产日韩欧美在线精品| 18禁美女被吸乳视频| 老司机午夜福利在线观看视频 | 欧美乱码精品一区二区三区| 国产欧美日韩一区二区三区在线| 日韩欧美一区视频在线观看| 香蕉丝袜av| 99久久精品国产亚洲精品| 成人精品一区二区免费| 亚洲av成人不卡在线观看播放网| 免费高清在线观看日韩| 亚洲色图 男人天堂 中文字幕| 制服人妻中文乱码| 久久久精品区二区三区| 1024香蕉在线观看| 一级片免费观看大全| 国产精品免费大片| 国产精品久久久久成人av| 黑人巨大精品欧美一区二区蜜桃| 十八禁网站网址无遮挡| 一区二区三区精品91| 国产男女内射视频| 欧美乱码精品一区二区三区| 另类精品久久| 欧美激情 高清一区二区三区| 亚洲九九香蕉| a在线观看视频网站| 成人影院久久| 久久国产精品影院| 日韩欧美一区视频在线观看| 亚洲国产欧美一区二区综合| 欧美成人免费av一区二区三区 | 99国产极品粉嫩在线观看| 国产在线免费精品| 亚洲久久久国产精品| 国产熟女午夜一区二区三区| 精品国产乱子伦一区二区三区| 国产在线一区二区三区精| 久久中文字幕人妻熟女| 狂野欧美激情性xxxx| 一级a爱视频在线免费观看| 精品乱码久久久久久99久播| 麻豆国产av国片精品| 一区二区三区乱码不卡18| 精品一区二区三卡| 51午夜福利影视在线观看| 国产亚洲欧美精品永久| 精品福利观看| 亚洲熟妇熟女久久| 一级片'在线观看视频| 国产视频一区二区在线看| 国产在线视频一区二区| 国产欧美日韩综合在线一区二区| 亚洲欧美日韩高清在线视频 | 国产福利在线免费观看视频| 成人18禁高潮啪啪吃奶动态图| 午夜免费成人在线视频| 午夜福利在线免费观看网站| 最新在线观看一区二区三区| 9热在线视频观看99| 99在线人妻在线中文字幕 | 欧美国产精品va在线观看不卡| 男男h啪啪无遮挡| www.自偷自拍.com| avwww免费| 69av精品久久久久久 | 国产精品影院久久| 国产伦人伦偷精品视频| 夫妻午夜视频| 欧美 日韩 精品 国产| 满18在线观看网站| 中文字幕最新亚洲高清| 成人特级黄色片久久久久久久 | 午夜激情久久久久久久| 欧美亚洲日本最大视频资源| 欧美日韩亚洲综合一区二区三区_| av一本久久久久| 纯流量卡能插随身wifi吗| 久久中文字幕人妻熟女| 久久久精品94久久精品| 亚洲精品中文字幕在线视频| 丝袜美腿诱惑在线| 色尼玛亚洲综合影院| 久久久精品区二区三区| 91av网站免费观看| 91麻豆av在线| 亚洲精品国产区一区二| 国产精品av久久久久免费| svipshipincom国产片| 999精品在线视频| 久久人妻福利社区极品人妻图片| 成人影院久久| av天堂在线播放| 这个男人来自地球电影免费观看| 热re99久久国产66热| 999久久久国产精品视频| 岛国毛片在线播放| 成人国语在线视频| 无限看片的www在线观看| 我要看黄色一级片免费的| 国产日韩欧美亚洲二区| 日韩视频一区二区在线观看| 高清毛片免费观看视频网站 | 另类亚洲欧美激情| 午夜精品久久久久久毛片777| 亚洲午夜理论影院| 久久精品aⅴ一区二区三区四区| 美女主播在线视频| 岛国在线观看网站| 十八禁网站网址无遮挡| 99久久精品国产亚洲精品| 黑人巨大精品欧美一区二区mp4| 国产成人精品在线电影| 首页视频小说图片口味搜索| 免费一级毛片在线播放高清视频 | 久久婷婷成人综合色麻豆| 精品福利永久在线观看| 免费一级毛片在线播放高清视频 | 亚洲黑人精品在线| 日日爽夜夜爽网站| 一区二区日韩欧美中文字幕| 国产成人av教育| 一级黄色大片毛片| 99精国产麻豆久久婷婷| 高潮久久久久久久久久久不卡| 人妻 亚洲 视频| 色婷婷av一区二区三区视频| 少妇 在线观看| 岛国毛片在线播放| 91麻豆av在线| 日本黄色视频三级网站网址 | 久久天躁狠狠躁夜夜2o2o| 亚洲欧美一区二区三区久久| 国产精品国产av在线观看| 久久这里只有精品19| 丝袜在线中文字幕| 中亚洲国语对白在线视频| 午夜91福利影院| 亚洲男人天堂网一区| 精品午夜福利视频在线观看一区 | 中文字幕另类日韩欧美亚洲嫩草| 精品第一国产精品| 亚洲国产看品久久| 亚洲国产av新网站| 亚洲精品国产区一区二| 电影成人av| 夜夜骑夜夜射夜夜干| 女性生殖器流出的白浆| 丰满迷人的少妇在线观看| 波多野结衣av一区二区av| 69精品国产乱码久久久| 老熟妇乱子伦视频在线观看| 免费观看av网站的网址| 怎么达到女性高潮| 国产av一区二区精品久久| 99热网站在线观看| 国产视频一区二区在线看| 中文字幕人妻丝袜制服| 久久 成人 亚洲| 国产一卡二卡三卡精品| 欧美黄色淫秽网站| 欧美日韩黄片免| 纵有疾风起免费观看全集完整版| 中文字幕av电影在线播放| 国产91精品成人一区二区三区 | 无遮挡黄片免费观看| 亚洲色图 男人天堂 中文字幕| 欧美精品亚洲一区二区| 中文字幕人妻丝袜制服| 王馨瑶露胸无遮挡在线观看| 久久中文看片网| 国产片内射在线| 欧美在线黄色| 热99re8久久精品国产| 国产无遮挡羞羞视频在线观看| 啦啦啦中文免费视频观看日本| 国产高清国产精品国产三级| 亚洲少妇的诱惑av| 在线观看免费日韩欧美大片| 老司机在亚洲福利影院| 十八禁网站网址无遮挡| 熟女少妇亚洲综合色aaa.| 午夜精品国产一区二区电影| 最近最新免费中文字幕在线| 咕卡用的链子| 成年人免费黄色播放视频| 亚洲精品国产色婷婷电影| 免费看a级黄色片| 男人操女人黄网站| 国产精品98久久久久久宅男小说| 免费观看人在逋| 亚洲色图综合在线观看| 亚洲视频免费观看视频| 国产成人欧美在线观看 | 丝袜喷水一区| 国产极品粉嫩免费观看在线| 亚洲精品在线美女| 亚洲av国产av综合av卡| 一本大道久久a久久精品| 在线天堂中文资源库| 久久久久精品国产欧美久久久| 在线观看一区二区三区激情| 热99re8久久精品国产| 亚洲国产成人一精品久久久| 久久精品亚洲av国产电影网| 欧美精品亚洲一区二区| 久久人妻熟女aⅴ| 亚洲 欧美一区二区三区| 宅男免费午夜| 三级毛片av免费| 又黄又粗又硬又大视频| 亚洲av电影在线进入| 国产av国产精品国产| 1024香蕉在线观看| 美女扒开内裤让男人捅视频| 亚洲成人免费av在线播放| 老司机午夜十八禁免费视频| 国产成人欧美| 黑丝袜美女国产一区| svipshipincom国产片| 亚洲第一av免费看| 久久久久久亚洲精品国产蜜桃av| 欧美亚洲 丝袜 人妻 在线| 一本综合久久免费| 19禁男女啪啪无遮挡网站| 50天的宝宝边吃奶边哭怎么回事| 1024香蕉在线观看| 伦理电影免费视频| 久久久欧美国产精品| 亚洲国产中文字幕在线视频| 亚洲国产看品久久| 狂野欧美激情性xxxx| 国产精品久久久av美女十八| 国产成人免费观看mmmm| 丁香六月天网| 99久久国产精品久久久| 天天躁夜夜躁狠狠躁躁| 亚洲精品自拍成人| 欧美日韩国产mv在线观看视频| 国产免费现黄频在线看| 中文字幕高清在线视频| 久久精品人人爽人人爽视色| 国产成人精品久久二区二区免费| 欧美激情极品国产一区二区三区| 久久青草综合色| 啦啦啦 在线观看视频| 国产精品电影一区二区三区 | 亚洲九九香蕉| 久久精品亚洲精品国产色婷小说| 亚洲性夜色夜夜综合| 窝窝影院91人妻| 国产一区二区激情短视频| 国产免费视频播放在线视频| 亚洲中文av在线| 久久青草综合色| 19禁男女啪啪无遮挡网站| 国产成人免费无遮挡视频| 捣出白浆h1v1| 国产精品电影一区二区三区 | cao死你这个sao货| 亚洲精品一卡2卡三卡4卡5卡| 一边摸一边抽搐一进一小说 | 黑人巨大精品欧美一区二区蜜桃| 欧美日韩一级在线毛片| 自线自在国产av| 精品一区二区三区四区五区乱码| 99久久人妻综合| 在线av久久热| 精品一区二区三卡| 欧美中文综合在线视频| 考比视频在线观看| 久久午夜综合久久蜜桃| 人妻一区二区av| av国产精品久久久久影院| 久久久久精品人妻al黑| 亚洲色图综合在线观看| 在线天堂中文资源库| 一区二区三区国产精品乱码| 美女高潮到喷水免费观看| 精品卡一卡二卡四卡免费| 黑丝袜美女国产一区| av电影中文网址| 久久人妻熟女aⅴ| 在线 av 中文字幕| 国产免费福利视频在线观看| 欧美日韩av久久| 免费看十八禁软件| 精品国产超薄肉色丝袜足j| 欧美日韩一级在线毛片| 自线自在国产av| 桃花免费在线播放| 搡老岳熟女国产| 国产男女内射视频| 另类亚洲欧美激情| 亚洲情色 制服丝袜| 日韩三级视频一区二区三区| 啦啦啦视频在线资源免费观看| 露出奶头的视频| 成人永久免费在线观看视频 | 9热在线视频观看99| 男女边摸边吃奶| 亚洲精品在线美女| 日本一区二区免费在线视频| 午夜免费成人在线视频| 大码成人一级视频| 日韩欧美一区视频在线观看| 国产人伦9x9x在线观看| 欧美日本中文国产一区发布| 国产一卡二卡三卡精品| 黑人猛操日本美女一级片| 国产精品国产av在线观看| 50天的宝宝边吃奶边哭怎么回事| 岛国毛片在线播放| 国产在视频线精品| 亚洲一码二码三码区别大吗| 在线观看66精品国产| 欧美日韩视频精品一区| 国产在线一区二区三区精| 亚洲成国产人片在线观看| 久久精品aⅴ一区二区三区四区| 久久热在线av| 日韩熟女老妇一区二区性免费视频| 91九色精品人成在线观看| 69av精品久久久久久 | 十八禁人妻一区二区| 女人被躁到高潮嗷嗷叫费观| 99久久精品国产亚洲精品| 精品国产一区二区三区四区第35| 国产片内射在线| 亚洲久久久国产精品| 免费一级毛片在线播放高清视频 | 欧美另类亚洲清纯唯美| kizo精华| 91精品国产国语对白视频| 免费女性裸体啪啪无遮挡网站| 亚洲午夜精品一区,二区,三区| 又大又爽又粗| 久久久久久人人人人人| 国产一区二区在线观看av| 亚洲国产毛片av蜜桃av| 夫妻午夜视频| 欧美激情极品国产一区二区三区| 成人永久免费在线观看视频 | 下体分泌物呈黄色| 1024视频免费在线观看| 久久久久久人人人人人| 日韩大片免费观看网站| 亚洲精品久久成人aⅴ小说| 91大片在线观看| 亚洲精品粉嫩美女一区| 欧美人与性动交α欧美软件| 国产黄色免费在线视频| 午夜精品国产一区二区电影| 一本—道久久a久久精品蜜桃钙片| 国产免费av片在线观看野外av| www.熟女人妻精品国产| 国产在线视频一区二区| 免费一级毛片在线播放高清视频 | 亚洲av日韩精品久久久久久密| 成人亚洲精品一区在线观看| 精品少妇黑人巨大在线播放| √禁漫天堂资源中文www| 两人在一起打扑克的视频|