• <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)的分布問題
    肉雞腸毒綜合癥防治
    热99在线观看视频| 日韩大尺度精品在线看网址| videossex国产| 免费高清视频大片| 日韩人妻高清精品专区| 一个人看视频在线观看www免费| 精品欧美国产一区二区三| 久久久久久久久久黄片| 成人精品一区二区免费| 99九九线精品视频在线观看视频| 搡女人真爽免费视频火全软件 | 日本欧美国产在线视频| 欧美性猛交╳xxx乱大交人| 亚洲欧美日韩无卡精品| 亚洲人成网站在线播放欧美日韩| 久久国产精品人妻蜜桃| aaaaa片日本免费| 亚洲精品久久国产高清桃花| 国产一区二区三区在线臀色熟女| 美女免费视频网站| 日韩中文字幕欧美一区二区| 日韩在线高清观看一区二区三区 | 老司机午夜福利在线观看视频| 男女视频在线观看网站免费| 在线免费观看的www视频| 91午夜精品亚洲一区二区三区 | 亚洲av成人av| 人妻夜夜爽99麻豆av| 婷婷亚洲欧美| 日韩一区二区视频免费看| 亚洲成人久久爱视频| 69av精品久久久久久| 男女边吃奶边做爰视频| 国产精品av视频在线免费观看| 亚洲专区中文字幕在线| 国语自产精品视频在线第100页| 男女那种视频在线观看| 亚洲avbb在线观看| 欧美xxxx性猛交bbbb| 成人无遮挡网站| 丰满乱子伦码专区| 日日撸夜夜添| 国产爱豆传媒在线观看| 白带黄色成豆腐渣| 亚洲中文日韩欧美视频| 又爽又黄a免费视频| 国产视频内射| 免费看光身美女| 干丝袜人妻中文字幕| 成年人黄色毛片网站| 色播亚洲综合网| 亚洲熟妇熟女久久| 99热只有精品国产| 日韩欧美精品v在线| 国产午夜福利久久久久久| .国产精品久久| 国产探花极品一区二区| 国产亚洲精品av在线| 精品久久久噜噜| 国产精品永久免费网站| 国产日本99.免费观看| 禁无遮挡网站| 天堂√8在线中文| 国产精品女同一区二区软件 | 麻豆国产97在线/欧美| 精品一区二区免费观看| 老女人水多毛片| 午夜爱爱视频在线播放| 精品久久久久久久末码| 嫩草影院新地址| 欧美性猛交╳xxx乱大交人| 亚洲国产欧洲综合997久久,| 久久久成人免费电影| 亚洲美女视频黄频| 亚洲国产精品sss在线观看| 2021天堂中文幕一二区在线观| 亚洲真实伦在线观看| 成人特级黄色片久久久久久久| 日韩精品中文字幕看吧| 欧美激情国产日韩精品一区| 国产伦在线观看视频一区| 国产高清三级在线| 亚洲精品乱码久久久v下载方式| 麻豆av噜噜一区二区三区| 两人在一起打扑克的视频| 午夜精品久久久久久毛片777| 嫁个100分男人电影在线观看| 一区二区三区四区激情视频 | av黄色大香蕉| 国产真实乱freesex| 亚洲aⅴ乱码一区二区在线播放| 我要看日韩黄色一级片| 亚洲自偷自拍三级| 国产精品国产三级国产av玫瑰| 最新中文字幕久久久久| 日韩,欧美,国产一区二区三区 | 深夜精品福利| 亚洲欧美日韩无卡精品| 人人妻,人人澡人人爽秒播| 女同久久另类99精品国产91| 国产免费男女视频| 亚洲,欧美,日韩| 内射极品少妇av片p| 国产激情偷乱视频一区二区| 免费人成在线观看视频色| 日本一本二区三区精品| 午夜久久久久精精品| 亚洲欧美精品综合久久99| 桃色一区二区三区在线观看| 亚洲成人久久爱视频| 亚洲av免费高清在线观看| 国产亚洲精品久久久com| 亚洲精华国产精华精| 亚洲av第一区精品v没综合| 久久6这里有精品| 丝袜美腿在线中文| 99国产精品一区二区蜜桃av| 一个人观看的视频www高清免费观看| 日本爱情动作片www.在线观看 | 久久精品久久久久久噜噜老黄 | 欧美区成人在线视频| 老司机深夜福利视频在线观看| 国产亚洲精品综合一区在线观看| 免费人成在线观看视频色| 大又大粗又爽又黄少妇毛片口| 国产高清三级在线| 亚洲人成网站高清观看| 观看免费一级毛片| 老师上课跳d突然被开到最大视频| 又爽又黄a免费视频| 国产又黄又爽又无遮挡在线| 18禁黄网站禁片免费观看直播| 亚洲美女黄片视频| 少妇丰满av| 久久精品国产亚洲av香蕉五月| 婷婷色综合大香蕉| 国产精品无大码| 中文字幕av成人在线电影| 一区二区三区高清视频在线| 精品一区二区三区人妻视频| 国产精品98久久久久久宅男小说| 日本成人三级电影网站| 全区人妻精品视频| 三级国产精品欧美在线观看| 网址你懂的国产日韩在线| 欧美日韩综合久久久久久 | 久久热精品热| av在线观看视频网站免费| 欧美又色又爽又黄视频| 久久热精品热| 午夜福利在线在线| 久久久午夜欧美精品| 一级a爱片免费观看的视频| 国内揄拍国产精品人妻在线| 久久精品国产亚洲av香蕉五月| 日日撸夜夜添| 国产精品免费一区二区三区在线| 久久久久久久久久黄片| 国产成人福利小说| 91久久精品国产一区二区三区| 久久精品综合一区二区三区| 乱人视频在线观看| 一个人免费在线观看电影| 日日干狠狠操夜夜爽| 久久午夜福利片| 美女cb高潮喷水在线观看| a在线观看视频网站| 99久久精品国产国产毛片| 狂野欧美白嫩少妇大欣赏| 三级毛片av免费| 非洲黑人性xxxx精品又粗又长| 国产单亲对白刺激| 精品日产1卡2卡| 老师上课跳d突然被开到最大视频| 久久人妻av系列| 国产伦精品一区二区三区四那| 99久久精品一区二区三区| 91在线观看av| 欧美精品啪啪一区二区三区| 亚洲人成网站在线播放欧美日韩| 亚洲精品日韩av片在线观看| 国产伦精品一区二区三区视频9| 久久精品夜夜夜夜夜久久蜜豆| 搡女人真爽免费视频火全软件 | 亚洲午夜理论影院| 亚洲无线观看免费| 亚洲av二区三区四区| av福利片在线观看| 女同久久另类99精品国产91| 国产女主播在线喷水免费视频网站 | 精品久久久久久久久久久久久| 欧洲精品卡2卡3卡4卡5卡区| 91麻豆av在线| 天堂网av新在线| 久久久久久久久久成人| 成人高潮视频无遮挡免费网站| 麻豆成人av在线观看| 中文字幕av成人在线电影| 97人妻精品一区二区三区麻豆| 久久婷婷人人爽人人干人人爱| 美女高潮喷水抽搐中文字幕| 久久久久久大精品| 男女视频在线观看网站免费| 好男人在线观看高清免费视频| 亚洲成a人片在线一区二区| 精品午夜福利视频在线观看一区| 级片在线观看| 亚洲精华国产精华液的使用体验 | 精品久久国产蜜桃| 午夜亚洲福利在线播放| 午夜老司机福利剧场| 国产中年淑女户外野战色| 免费黄网站久久成人精品| 国产精品一区二区免费欧美| 欧美人与善性xxx| 女同久久另类99精品国产91| 婷婷精品国产亚洲av在线| 久久久成人免费电影| 亚洲无线在线观看| 亚洲欧美日韩卡通动漫| 亚洲四区av| 美女高潮的动态| 蜜桃久久精品国产亚洲av| 国产精品野战在线观看| 成人永久免费在线观看视频| h日本视频在线播放| 久久热精品热| 色综合亚洲欧美另类图片| 人人妻人人澡欧美一区二区| 一进一出好大好爽视频| 免费观看人在逋| 久久精品国产鲁丝片午夜精品 | 中文字幕免费在线视频6| 日本免费一区二区三区高清不卡| 色精品久久人妻99蜜桃| 人妻制服诱惑在线中文字幕| 日韩强制内射视频| 成人欧美大片| 久久99热这里只有精品18| 国产麻豆成人av免费视频| 又紧又爽又黄一区二区| 国内毛片毛片毛片毛片毛片| 亚洲人成网站高清观看| 一a级毛片在线观看| 精品人妻熟女av久视频| 国产一级毛片七仙女欲春2| 天美传媒精品一区二区| 国产伦精品一区二区三区四那| 草草在线视频免费看| 国产成人av教育| 熟女电影av网| 最近最新免费中文字幕在线| 在线播放无遮挡| 五月玫瑰六月丁香| 国产aⅴ精品一区二区三区波| 美女大奶头视频| 亚洲成人久久爱视频| 精品不卡国产一区二区三区| 男女视频在线观看网站免费| 黄色配什么色好看| 极品教师在线免费播放| 午夜免费男女啪啪视频观看 | 亚洲欧美清纯卡通| 变态另类成人亚洲欧美熟女| 少妇猛男粗大的猛烈进出视频 | 亚洲三级黄色毛片| 亚洲av免费在线观看| 国产v大片淫在线免费观看| 精品人妻偷拍中文字幕| 最近在线观看免费完整版| 亚洲av五月六月丁香网| 黄色视频,在线免费观看| 国内精品久久久久精免费| 少妇丰满av| 在线看三级毛片| av在线蜜桃| 中亚洲国语对白在线视频| 亚洲,欧美,日韩| 久久久久久国产a免费观看| 级片在线观看| 午夜福利欧美成人| 精品午夜福利视频在线观看一区| 深夜a级毛片| 美女 人体艺术 gogo| 国产亚洲精品久久久久久毛片| 欧美成人性av电影在线观看| 中文资源天堂在线| 综合色av麻豆| 国产欧美日韩精品一区二区| 精品午夜福利在线看| 蜜桃久久精品国产亚洲av| 国产精品无大码| 三级毛片av免费| 身体一侧抽搐| 精品欧美国产一区二区三| av福利片在线观看| 精品久久久久久久末码| 亚洲经典国产精华液单| 亚洲av熟女| 亚洲性久久影院| 九色成人免费人妻av| 国产综合懂色| 免费av毛片视频| 91精品国产九色| 午夜福利在线在线| 在线看三级毛片| 久久热精品热| 国产麻豆成人av免费视频| 国内精品美女久久久久久| 在线观看午夜福利视频| 小蜜桃在线观看免费完整版高清| 国产黄a三级三级三级人| 国产色爽女视频免费观看| 日韩精品有码人妻一区| 国产精品久久久久久久久免| 亚洲成人精品中文字幕电影| 桃色一区二区三区在线观看| 国产精品福利在线免费观看| 欧美3d第一页| 亚洲美女黄片视频| 久久久国产成人精品二区| 国产高潮美女av| 男人和女人高潮做爰伦理| 国产精品亚洲一级av第二区| 久久午夜福利片| 白带黄色成豆腐渣| 亚洲精品影视一区二区三区av| 国产探花极品一区二区| 一区二区三区四区激情视频 | 小蜜桃在线观看免费完整版高清| 欧美性猛交黑人性爽| 日日夜夜操网爽| 在线免费观看的www视频| 亚洲美女搞黄在线观看 | 性插视频无遮挡在线免费观看| 老司机深夜福利视频在线观看| 国产精品爽爽va在线观看网站| 99热这里只有是精品50| 亚洲精品色激情综合| 亚洲国产精品久久男人天堂| 精品一区二区免费观看| 国产 一区 欧美 日韩| 他把我摸到了高潮在线观看| 国产综合懂色| 精品午夜福利在线看| 日本一本二区三区精品| 欧美精品国产亚洲| 日日摸夜夜添夜夜添小说| 99在线人妻在线中文字幕| 国产色婷婷99| 九色成人免费人妻av| 国产探花在线观看一区二区| 九色成人免费人妻av| 精品国内亚洲2022精品成人| 欧美日韩中文字幕国产精品一区二区三区| 国产伦精品一区二区三区四那| 床上黄色一级片| 精品久久久噜噜| 天堂动漫精品| 国产激情偷乱视频一区二区| 噜噜噜噜噜久久久久久91| 午夜视频国产福利| a级毛片a级免费在线| 国产激情偷乱视频一区二区| 嫩草影院精品99| 欧美日本视频| 少妇裸体淫交视频免费看高清| 中国美白少妇内射xxxbb| 欧美性猛交╳xxx乱大交人| 亚洲无线观看免费| 岛国在线免费视频观看| 白带黄色成豆腐渣| 亚洲精品久久国产高清桃花| 国产极品精品免费视频能看的| 亚洲欧美清纯卡通| 99久久久亚洲精品蜜臀av| 成人av在线播放网站| 久久久久久九九精品二区国产| 好男人在线观看高清免费视频| 国内精品美女久久久久久| 最近在线观看免费完整版| 成人av一区二区三区在线看| 18禁黄网站禁片午夜丰满| 此物有八面人人有两片| 99热这里只有是精品在线观看| 2021天堂中文幕一二区在线观| 最新在线观看一区二区三区| 美女免费视频网站| 直男gayav资源| 久久精品国产清高在天天线| 人妻丰满熟妇av一区二区三区| 欧美黑人欧美精品刺激| a级一级毛片免费在线观看| 亚洲精品国产成人久久av| 91在线观看av| 女人十人毛片免费观看3o分钟| 老司机深夜福利视频在线观看| av视频在线观看入口| 国产白丝娇喘喷水9色精品| 久久亚洲真实| 国产视频一区二区在线看| 午夜老司机福利剧场| 亚洲电影在线观看av| 中亚洲国语对白在线视频| 特级一级黄色大片| 亚洲国产欧美人成| 日本 av在线| 能在线免费观看的黄片| 国产精品国产高清国产av| 亚洲av日韩精品久久久久久密| 久久久久久久久久久丰满 | 日韩中字成人| 在现免费观看毛片| 亚洲最大成人中文| 哪里可以看免费的av片| 国产亚洲精品综合一区在线观看| 亚洲七黄色美女视频| h日本视频在线播放| 亚洲aⅴ乱码一区二区在线播放| 国产精品1区2区在线观看.| 桃红色精品国产亚洲av| 日本免费一区二区三区高清不卡| 干丝袜人妻中文字幕| 色噜噜av男人的天堂激情| 黄色丝袜av网址大全| 88av欧美| 免费av观看视频| 一夜夜www| 亚洲av二区三区四区| 国产成人影院久久av| 国产一级毛片七仙女欲春2| 99精品久久久久人妻精品| 日韩人妻高清精品专区| 久久99热6这里只有精品| av国产免费在线观看| 校园春色视频在线观看| 美女高潮喷水抽搐中文字幕| 看黄色毛片网站| 午夜亚洲福利在线播放| 国产老妇女一区| 99在线视频只有这里精品首页| 日韩精品青青久久久久久| 麻豆久久精品国产亚洲av| 日韩亚洲欧美综合| 亚洲精品一区av在线观看| 人妻久久中文字幕网| 久久99热这里只有精品18| 欧美一级a爱片免费观看看| 女生性感内裤真人,穿戴方法视频| 日日摸夜夜添夜夜添av毛片 | 日本黄色视频三级网站网址| 热99re8久久精品国产| 成人三级黄色视频| 99久久成人亚洲精品观看| 免费高清视频大片| 91狼人影院| 亚洲五月天丁香| 麻豆久久精品国产亚洲av| xxxwww97欧美| 九九久久精品国产亚洲av麻豆| 久久久精品大字幕| 18禁在线播放成人免费| 一个人看的www免费观看视频| 在现免费观看毛片| 色在线成人网| 国产综合懂色| 亚洲狠狠婷婷综合久久图片| 日本熟妇午夜| 欧美黑人巨大hd| 在线观看66精品国产| 亚洲av免费高清在线观看| 老熟妇仑乱视频hdxx| 国产成人aa在线观看| 国产午夜精品久久久久久一区二区三区 | 亚洲va在线va天堂va国产| 22中文网久久字幕| 成熟少妇高潮喷水视频| 国产亚洲91精品色在线| 久久久久免费精品人妻一区二区| 国产精品女同一区二区软件 | 18禁黄网站禁片午夜丰满| 少妇人妻一区二区三区视频| 观看免费一级毛片| 国产精品自产拍在线观看55亚洲| 欧美不卡视频在线免费观看| 久久午夜亚洲精品久久| 亚洲精品一区av在线观看| 别揉我奶头 嗯啊视频| 变态另类丝袜制服| 看免费成人av毛片| 1024手机看黄色片| 亚洲成人精品中文字幕电影| 老女人水多毛片| 露出奶头的视频| 成年版毛片免费区| 精品一区二区三区视频在线| 色在线成人网| av在线老鸭窝| 听说在线观看完整版免费高清| 亚洲av二区三区四区| 欧美色视频一区免费| 日日啪夜夜撸| 亚洲国产精品久久男人天堂| 久久精品综合一区二区三区| 精品福利观看| 亚洲人与动物交配视频| 热99在线观看视频| 午夜视频国产福利| 亚洲五月天丁香| 国产高潮美女av| 婷婷丁香在线五月| 99在线人妻在线中文字幕| 内射极品少妇av片p| 午夜福利视频1000在线观看| 亚洲美女黄片视频| 国产高清三级在线| 日本免费a在线| 免费观看在线日韩| 一级av片app| 亚洲18禁久久av| 中文字幕av在线有码专区| 国产精品久久久久久av不卡| 日本-黄色视频高清免费观看| 很黄的视频免费| 一进一出好大好爽视频| 亚洲国产日韩欧美精品在线观看| 午夜福利在线观看吧| 日韩欧美国产一区二区入口| 精品99又大又爽又粗少妇毛片 | 日韩,欧美,国产一区二区三区 | www日本黄色视频网| 成熟少妇高潮喷水视频| 精品久久国产蜜桃| 亚洲av二区三区四区| 三级毛片av免费| 国产免费av片在线观看野外av| av专区在线播放| av视频在线观看入口| 国产高清视频在线播放一区| 亚洲久久久久久中文字幕| 有码 亚洲区| 白带黄色成豆腐渣| 热99re8久久精品国产| 国产蜜桃级精品一区二区三区| 国产激情偷乱视频一区二区| 最近最新免费中文字幕在线| 午夜福利在线在线| 观看免费一级毛片| 一夜夜www| 99久久久亚洲精品蜜臀av| 亚洲avbb在线观看| 丰满乱子伦码专区| 九九在线视频观看精品| 免费看av在线观看网站| 成人三级黄色视频| 一本精品99久久精品77| 日本免费a在线| 免费在线观看影片大全网站| 亚洲avbb在线观看| 淫秽高清视频在线观看| 日韩中文字幕欧美一区二区| 成年女人毛片免费观看观看9| 天堂√8在线中文| 免费无遮挡裸体视频| 一级黄片播放器| 亚洲一级一片aⅴ在线观看| 一个人看视频在线观看www免费| 精品一区二区三区视频在线| 热99在线观看视频| 中国美白少妇内射xxxbb| 国内精品宾馆在线| avwww免费| 免费av观看视频| 久久国内精品自在自线图片| 精品久久久久久,| 啦啦啦韩国在线观看视频| 毛片一级片免费看久久久久 | 午夜精品在线福利| av国产免费在线观看| 99久久精品一区二区三区| 久久久久国内视频| 少妇人妻一区二区三区视频| 成人欧美大片| 成人无遮挡网站| 欧美日本亚洲视频在线播放| 日韩大尺度精品在线看网址| 欧美成人a在线观看| 国产三级在线视频| 大又大粗又爽又黄少妇毛片口| 久久久久久久精品吃奶| 人人妻,人人澡人人爽秒播| 日本一本二区三区精品| 亚洲国产精品成人综合色| 嫩草影视91久久| 久久人妻av系列| 99国产精品一区二区蜜桃av| 日日啪夜夜撸| 国产色爽女视频免费观看| 久久久久久久精品吃奶| 成人国产一区最新在线观看| 国产色爽女视频免费观看| 亚洲aⅴ乱码一区二区在线播放| 久久久午夜欧美精品| 又黄又爽又刺激的免费视频.| 久久久久久久精品吃奶| 国产av在哪里看| 欧美最黄视频在线播放免费| 欧美黑人巨大hd| 精品久久久久久,| 免费一级毛片在线播放高清视频| 国模一区二区三区四区视频| 91av网一区二区| 午夜日韩欧美国产| 高清毛片免费观看视频网站| 久久国产乱子免费精品|