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

    A study of the influence of different grid structures on plasma characteristics in the discharge chamber of an ion thruster

    2023-03-09 05:45:42MingmingSUN孫明明JianfeiLONG龍建飛WeilongGUO郭偉龍ChaoLIU劉超andYongZHAO趙勇
    Plasma Science and Technology 2023年1期
    關(guān)鍵詞:劉超趙勇

    Mingming SUN(孫明明),Jianfei LONG(龍建飛),Weilong GUO(郭偉龍),Chao LIU(劉超) and Yong ZHAO(趙勇)

    1 Science and Technology on Vacuum Technology and Physics Laboratory,Lanzhou Institute of Physics,Lanzhou 730000,People’s Republic of China

    2 Gravitational Wave Universe Taiji Laboratory,Hangzhou Institute for Advanced Study,UCAS,Hangzhou 310024,People’s Republic of China

    3 Key Laboratory of Gravitational Wave Precision Measurement of Zhejiang Province,UCAS,Hangzhou 310024,People’s Republic of China

    Abstract The grid structure has significant effects on the discharge characteristics of an ion thruster.The discharge performances of a 30 cm diameter ion thruster with flat,convex and concave grids are studied.The analysis results show that the discharge chamber with a convex grid has a larger‘magnetic-field free area’than the others,and the parallelism of the magnetic-field isopotential lines and anode is generally the same in the three models.Plasma densities of the three structures at the grid outlet are in the range of 3.1×1016–6.9×1017 m?3.Along the thruster axis direction,the electron temperature in the chamber with the convex and concave grids is in the range of 3.3–3.5 eV,while that with a flat grid is lower,in the range of 3.1–3.5 eV.In addition,the convex and the concave grids have better uniform distribution of electron temperature.Moreover,the collision frequency ratios show that the axial degree of ionization of the three models is the highest,and the flat grid has the highest discharge efficiency,followed by the convex grid and the concave grid is the least efficient.The test and simulation results of the 30 cm diameter ion thruster with the convex grid show that the measurement and calculation results are 3.67 A and 3.44 A,respectively,and the error above mainly comes from the ignorance of the doubly charged ions and parameter settings in the model.The comparison error between the simulation and measurement of beam current density is mainly caused by the actual thermal deformation of the grids during the discharge process,which leads to the change in electric potential distribution and variation of the focusing characteristics of the grids.Upon consideration of discharge performance and the thermal grid gap variation,it can be concluded that the flat and concave grids are more suitable for small-diameter ion thrusters,while the convex grid is a more reasonable choice for the higher-power and larger-diameter thrusters.

    Keywords:ion thruster,the grids,discharge characteristics

    1.Introduction

    As the core component of an ion thruster for the extraction and acceleration of ions,the structure of the grids(also called an ionic optical system)has a significant influence on the generation and transport of plasma in the discharge chamber and ionic optical system[1–3].The focusing of the ion beam is formed in the upstream and downstream regions of the grids[4,5].Meanwhile,due to the potential of the grids being the lowest of the entire discharge region,the electric field formed between the grid and the anode has a decisive effect on plasma characteristics and distribution in the chamber,and affects the thruster performance,such as thrust,specific impulse and efficiency.Figure 1 shows the structure of the discharge chamber of a 30 cm diameter ion thruster developed by Lanzhou Institute of Physics(LIP)and a schematic diagram of the ion beam extraction process.As shown in figure 1,the 30 cm diameter ion thruster adopts the triple grid as an ionic optical system[6,7].The discharge chamber structure is a typical combination of cylindrical and conical sections,and the screen grid and the accelerator grid load 1200 V and ?220 V voltages,respectively.Meanwhile,the decelerator grid,as the charge exchange(CEX)ions bombardment target only,does not apply a voltage but improves the service life of the accelerator grid.

    Figure 1.A schematic diagram of the ion beam extraction process.

    Figure 2.DCMs with a flat grid(a),convex grid(b)and concave grid(c)(unit:m).

    In 2005,Wirz and Katz investigated the plasma characteristics in the discharge chamber of the 30 cm diameter NSTAR ion thruster with a convex grid using a multicomponent hybrid 2D computational discharge chamber model(DCM)[8].The purpose of the study was to calculate plasma parameters for the long-life requirement of the thruster in orbit to carry out plasma etching simulation and life prediction for the core components,mainly the hollow cathode and the grids.The results show that under TH15 working conditions,the electron density is in the range of 1×1016–7×1017m?3when the anode flow,cathode flow,screen voltage and acceleration voltage are 23.4 sccm,3.73 sccm,1100 V and ?249 V,respectively.The electron temperature along the thruster axis is the highest,reaching 4.3 eV.The closer it is to the anode,the lower the electron temperature is,about 4 eV.The electron collision frequency ratio(ratio of electrons–neutrals and electrons–ions)is highest close to the anode,reaching 4.7–8.1,while the closer to the thruster axis,the higher the electron–ion collision frequency.The ion current density at the center of the grid is the highest,reaching 40 A m?2while,along the radial direction,the ion current density presents a parabolic decrease.The DCM results show that the ionization process can be more uniform throughout the discharge chamber by allowing the primary electrons to move away from the thruster axis:that is,the ionization efficiency can be improved by reducing the densities of the primary electrons and doubly charged ions on the thruster axis.In 2006,Goebel[9]developed ‘modified’NSTAR ion thruster designs that guide primary electrons away from the thruster center line to improve the performance.In addition,these designs had been validated by experiments in 2005[10],and the measurements of doubly charged ions and neutral density along the thruster axis are in good agreement with the simulations by the DCM.These previous research studies have efficiently simulated the plasma characteristics in the discharge chamber of the ion thruster,and proved the effectiveness of the simulation model in the optimization design of the thruster.However,it takes a long time to build a simulation model,and the cost is high;therefore,it is worth considering the use of mature commercial software to quickly obtain the parameters of discharge plasma to improve the optimization of a thruster or verify different structural designs.

    The main purpose of this work is to establish a simple and efficient discharge model with a certain accuracy using general commercial software.Then,the effects of different grid structures on plasma parameters are calculated and evaluated,and the results of the model are verified by test.According to the simulation results of the plasma parameters of different grid structures,theoretical support will be provided for the possible design improvement of the existing grid structure in the future.

    2.Discharge model and boundary settings

    Due to the high costs associated with processing grid components with different structures,a numerical simulation method is used to obtain the influence of the grid structure on the plasma characteristics of the discharge chamber.COMSOL multi-physics coupling software is used in the simulation process,mainly because it is a useful tool for carrying out quick and inexpensive parametric studies of the effects of geometry,but the convergence of calculation needs to be solved[11].

    As shown in figure 2,three DCMs with different grid structures are established:a flat grid,a convex grid and a concave grid.Since the study mainly focuses on the plasma characteristics in the discharge chamber,the accelerator grid does not appear in the three models.The main dimensions of the three discharge chambers are generally consistent,including the position of the hollow cathode,gas inlet and the magnet pole,and the obvious difference is the location of the screen grid.As shown in figure 2(a),the DCM is composed of column and cone structures,which are all made of tantalum material due to the small thermal deformation of titanium at high temperatures.The radius of the discharge chamber outlet,also the screen grid radius,is 0.143 m.Four annular magnets made of permanent magnetic material form a ringcusp field,which makes the discharge process more uniform[12,13].Middle magnet 2 is nearly mounted in the middle of the cone structure,and the downstream magnet is installed at the bottom of the discharge chamber.The upstream magnet is mounted near the outlet of the column structure and the location of middle magnet 1 is at the joint of the cone and column structures.The diameters of the cathode gas inlet and the anode gas inlet are 1 mm and 1.5 mm,respectively.

    All three DCMs are axially symmetric;the axial length of the model in figure 2(a)is 0.16 m,which is from the chamber outlet to the cathode outlet,and the radial length is 0.14 m,from the symmetry axis to the edge of the chamber.Figure 2(b)shows the DCM with the convex grid;the main difference is the grid shape compared with the model in figure 2(a).The arch height of the convex grid is 0.017 m,and the position and sizes of other parts are consistent with figure 2(a).Similarly,the arch height of the concave grid shown in figure 2(c)is ?0.017 m.

    Figure 3.Potential distribution of discharge models with a flat grid(a),convex grid(b)and concave grid(c)(unit:V).

    The boundary settings include the electric parameters,magnetic properties of materials,the gas supply and plasma discharge settings.The boundary settings have a significant effect on the discharge process,accuracy and convergence of the simulation results.Therefore,the working parameters of the thruster and the appropriate adjustment of the parameters according to the convergence of the simulation results should be considered.Table 1 gives the typical working conditions,such as the electric parameters,magnetic properties and gas supply settings in 5 kW operation conditions of the 30 cm diameter ion thruster.

    Table 1.Working conditions of a 30 cm diameter thruster.

    Table 2.Parameter settings of the DD module.

    As shown in table 1,the discharge model is based on the thruster power supply relationship,and the cathode voltage is floating on the screen grid voltage.Then,the screen grid voltage is 1200 V;therefore,the keeper voltage is 1209.5 V and the anode voltage is 1235.8 V.The electric field can be expressed as the gradient of the electric potential,E= ??V,and the negative sign comes from the convention that the electric field always points in the direction of ion motion.In addition,the electric field simulation is based on the assumption that plasma is quasi-steady and the self-consistent electric field of plasma is ignored(the self-consistent potential of the plasma is in the range of 0–2 V and barely affects the particle motion[14]).The static electric field(ES)module in COMSOL is used to solve the potential distribution within the discharge chamber.

    Although COMSOL is convenient discharge simulation software,the convergence of simulation calculation is a prominent problem during use,especially in the fluid simulation of the discharge chamber,with the size of the gas inlet,fluid property settings and pressure conditions,which all have a significant impact on the convergence.Therefore,according to the real size and the gas supply of the thruster,the modeling and boundary conditions are set first,and then the boundaries are slightly changed according to the convergence of the calculated results.For the cathode,the diameter of the orifice is 1 mm.The keeper is a tantalum round plate with a central aperture,and the thickness and the diameter of the aperture are 1 mm and 5 mm,respectively.Neutral xenon gas is injected into the inlet of the cathode tube through a gas insulator(a ceramic component connected to the inlet of the cathode to isolate the gas discharge process toward the supply tube).The temperature of the cathode orifice is heated to over 1500 °C during the discharge process,and then the thermal electrons emit from the orifice and create thermal electron currents[15].As shown in figure 1(a)and table 1,the diameter and mass flow of the cathode are set as 1 mm and 0.2 mg s?1,respectively.Based on the same method,the diameter and mass flow of the anode are set as 1.5 mm and 5.4 mg s?1,respectively.The SPF module in COMSOL is used to solve the pressure and velocity of the fluid within the discharge chamber.According to previous test results of the same type of cathode[16],before the discharge begins,the velocity of neutral gas in the discharge chamber is at a low level,and below 0.3 Mach number.Therefore,the fluid in the discharge zone is considered as incompressible and turbulence is ignored.The temperature of the fluid is set as 300 K to obtain the fluid parameters before gas discharge.All the walls of the discharge chamber are set as the no-slip boundary,which means the fluid velocity on the wall is 0.The boundaries of the inlets of the anode and the cathode are set as a mass flow condition with values of 5.4 mg s?1and 0.2 mg s?1,respectively,and the flow direction is perpendicular to the inlet section.The initial pressure of the discharge zone is set as 0.01 Pa,and the internal gas is assumed to be a steady flow.The screen grid shown in figure 1 is defined as the gas outlet,and the boundary of the gas outlet is set as a pressure condition with the value of 0.04 Pa according to previous test results[16].

    At present,the main plasma discharge simulation methods of ion thrusters include PIC-MCC(Particle-in-Cell-Monte-Carlo-Collision)and fluid simulation.The PIC-MCC simulation is relatively complex and has higher computational resource requirements and longer calculation time than the fluid simulation method.The fluid simulation method is relatively simple,but requires a large number of parameters to complete the calculations.The drift-diffusion(DD)module in COMSOL is a typical fluid simulation method based on the DD equation,which regards the motion process of charged plasma in the discharge chamber as a fluid flow process under the influence of an electric field and a magnetic field.Therefore,the DD module is adopted to obtain the characteristics of the plasma in the discharge chamber.

    With regard to the potential of the plasma sheath in the discharge chamber,according to the conclusion of Goebel[17],for a quasi-neutral ionization process,the potential of the sheath can be expressed as equation(1)by assuming that the ion density is equal to electron density,that isni=ne

    In equation(1),kis the Boltzmann constant,Teis electron temperature(1 eV is 11600 K),Mandmare the masses of the ion and electron,respectively,andAaandAare the areas of electron loss and the total area of ion loss,respectively.In this paper,AaandAare the areas of the anode and screen,respectively,by default,and the electron temperature is assumed to be 5 eV in the calculation(refer to the electron temperature of the NSTAR ion thruster when the discharge reaches a steady state[8]);thus,the sheath potential of the discharge model is about 36.5 V,which is essentially equal to the anode potential shown in table 1(the potential also floats on the screen grid).Therefore,the anode is taken as the sheath boundary in the simulation to simplify the calculation.

    The DD module includes the electron continuity equation,electron DD equation,electron flux equation and energy equation.The electron continuity equation is shown as:

    whereneis electron density,Γeis the electron flux in the cathode discharge zone,uis the electron velocity,Reis the electron generation rate and the unit is m?3s?1.The electrons are generated from elastic collision,excitation collision and ionization collision,and the expression of the reaction rate coefficient for xenon in different collision types has been given by Milleret al[18].The different reaction rate coefficients are denoted byk1,k2andk3,respectively.Meanwhile,all the different reaction rate coefficients have a similar expressed equation as 〈σν 〉(whereσ is the collision cross section,andν is the frequency of collisions between particles and is determined by the electron temperatureTe[19]).The elastic collision is ignored because there are almost no Maxwellian electrons being produced and,in most cases,the energy transfer is predominant.Then,the electronic excitation reaction rater2and ionization reaction rater3can be expressed ask2n0neandk3n0ne,respectively.Therefore,Reis expressed as:

    The electron flux Γecan be given by the DD equation,which is shown as:

    After all the boundary conditions of the discharge model are obtained,the solver setting is considered.The simulations in the ES and SPF modules are steady simulation.The simulation in the DD module is a transient simulation,and the time range is set from 10?8to 10?1s to ensure that the plasma properties are evolving with time until steady state conditions are reached.

    3.Simulation results and analysis

    The potential distributions in the discharge chamber of the three different grid structures are shown in figures 3(a)–(c).The movement direction of ions is perpendicular to the isopotential line,and the potential distribution of the convex and concave grids near the screen grid is opposite from the comparison results,which are shown in figures 3(b)and(c).Therefore,the ion beams extracted from convex and concave grids are divergent and focused,respectively.Because the focused ion beam is extracted from the grids and neutralized by the electrons emitted by the neutralizer,then,under the diffusion of fluid and the interaction between particles,the beam extraction process in figure 3(c)will therefore present the characteristics of focusing first and then diverging.The potential distribution shown in figure 3(a)determines that the beam extracted by the flat grid will present perpendicular to the grid outlet surface.

    Figures 4(a)–(c)show the magnetic-field distribution of the three different discharge chamber structures.As the results show in figure 4,there is no significant difference in the magnetic-field distributions of the three discharge chamber structures.The green line in figure 4 is the 50 Gs magneticfield isopotential,and the region enclosed by this line is called the ‘magnetic-field free area’[13].The magnetic field in this region has less of an effect on the plasma than the magnetic field near the anode,and the larger the volume of the ‘free area’,the better the flatness of the ion beam.Meanwhile,the closer and more parallel the edge of the ‘free area’ is to the anode,the better the plasma uniformity and the ion current density profile that can be obtained[21],and the influence of plasma on the anode can be reduced.It is obvious that the discharge chamber with the convex grid has the largest ‘free area’,but this does not mean that the larger the discharge chamber,and the larger the ‘free area’,the better the plasma uniformity.In addition,the combined effect of fluid diffusion,the magnetic pole position,and the continuity of the 50 Gs magnetic isopotential line should be considered.Therefore,the discharge chamber has the optimum ratio of length to diameter.In addition,due to the same position and boundary settings of the magnets,the parallelism of the magnetic-field isopotential lines and anode is essentially the same in the three models.

    Figure 4.Magnetic-field distribution of the discharge models with a flat grid(a),convex grid(b)and concave grid(c)(unit:Gs).

    Figure 5.Fluid velocity distribution of the discharge models with a flat grid(a),convex grid(b)and concave grid(c)(unit:m s?1).

    The simulation results of fluid velocity in the discharge chamber are shown in figures 5(a)–(c),and it can be concluded that fluid velocity is barely affected by the structure of the discharge chamber.The fluid velocity is the highest along the axis of the hollow cathode,and the value is in the range of 0.1–0.2 m s?1.The fluid ejecting from the anode gas inlet shows obvious diffusion characteristics,mainly because the size of the anode gas inlet is higher than that of the cathode gas inlet,and the gas mass flow is higher.In the previous boundary condition settings,the walls of the discharge chamber are all set as no-slip.In the calculation,it is further found that if the inertia term of the fluid,i.e.the friction inside the fluid,is not ignored,the calculation results can barely achieve convergence.Due to the low gas pressure in the discharge chamber,the friction inside the fluid is ignored.In addition,it is found that if the fluid outlet does not prohibit backflow,the convergence of the calculation results is difficult.According to previous calculation results[16],the backflow from the closer plume region(1–2 times the grid diameter from the grid surface)to the discharge chamber at room temperature only accounts for 10%of the total gas flow into the chamber.When the thruster is operating,the proportion of backflow will be lower.Therefore,to ensure the calculation process,the backflow is ignored.The neutral density calculated by the SPF module is in the range of 2.8×1017–7.2×1018m?3,and the neutral density in most areas is 3.2×1018m?3.

    Figures 6(a)–(c)show the plasma generation process in a transient state of the discharge chamber with the flat grid at 2×10?7,1×10?4and 0.1 s,respectively.As shown in figure 6(a),at 2×10?7s,which is the initial stage of the discharge process,the plasma density along the thruster axis(as shown in figure 1(a))from the cathode to the screen grid is within the range of 1.1×1017–4.1×1017m?3.The plasma density is the highest at the cathode outlet,and the trend of moving toward the magnetic pole is obvious.This is because the primary electrons emitted from the cathode drift along the magnetic-field line(E×Bmotion),and eventually deposit to the magnetic poles.As the ionization process progresses,the plasma density distribution results after discharging for 1×10?4s are shown in figure 6(b).At this time,the ionization has expanded to most areas of the discharge chamber,and the plasma density along the thruster axis is in the range of 1.6×1017–5.7×1017m?3.Figure 6(c)shows the plasma density distribution at 0.1 s;it can be concluded that the discharge process has reached stability,and the plasma distribution in the discharge chamber does not change with time.In addition,it is evident that both primary and secondary electrons are deposited to magnetic poles and form an anode current.The plasma density along the thruster axis is in the range of 1.7×1017–6.5×1017m?3.

    Figure 6.Distribution of plasma density at 2×10?7 s(a),1×10?4 s(b)and 0.1 s(c)(unit:m?3).

    Figure 7.Simulation results of plasma density(a)and electron temperature(b).

    Figure 7(a)gives the plasma density distribution along the radial direction of the screen grid after the discharge processes of the three different chamber structures reach stability.As shown in the figure,the plasma density of the three models is the highest in the center of the discharge chamber,and the density decreases as it gets closer to the edge.This is due to the combined influence of the potential and the neutral gas distribution,as shown in figures 3 and 5.

    Figure 7(b)shows the electron temperature distribution along the axial region of the three different discharge chambers in the stability discharge process.It can be concluded that the electron temperature of the chamber with the flat grid is the lowest,within the range of 3.1–3.5 eV,while the electron temperature of the convex and concave grids is not much different,which is essentially in the range of 3.3–3.5 eV along the axis.To better compare the electron temperature distributions of the different chamber structures,the distribution diagram is given in figures 8(a)–(c).The results show that the higher electron temperatures of the three structures are essentially in the axial region of the thruster,while the electron temperature decreases in the region close to the upstream magnet.This is due to the low energy of the electrons emitted from the cathode,generally in the range of 2–3 eV.When the electrons are accelerating,the magnetic-field intensity of the axis of the discharge chamber is lower,and the average free path of electrons is higher;therefore,the collision probability is lower.This makes the electrons on the axis have a longer acceleration time,and the electron temperature is higher.However,the region near the anode is the main ionization region;therefore,the electron temperature is lower.In addition,the difference in electron temperature near the grid outlet of the three models is mainly caused by the distribution of potential and plasma density,as shown in figures 3 and 7(a),respectively.As shown in figure 3,the closer it is to the grid outlet,the lower the electric field acceleration effect is.Meanwhile,the region near the grid outlet is also the ion beam extraction area and,with the extraction process,the energy to maintain the discharge becomes lower;therefore,the potential of the plasma sheath decreases,resulting in the reduction of electron temperature.

    Figure 8.Electron temperature distribution of the discharge models with a flat grid(a),convex grid(b)and concave grid(c)(unit:eV).

    Figure 9.The electron collision frequency ratios(ν e nνe i )of the discharge models with a flat grid(a),convex grid(b)and concave grid(c).

    Figure 10.The Faraday probe array(a),equivalence of the ion beam(b)and comparison between the calculated results and measured results of ion beam density(c).

    Finally,collisions between particles are an important part of the analysis,especially between electrons and neutral atoms,and between electrons and ions.The collision frequency of different particles can reflect the degree of ionization of different regions in the discharge chamber.The electron collision frequency ratio(νenνei)is shown in figures 9(a)–(c),and the ratio unit is dimensionless.The ratio is smaller,and the higher the collision frequency between electrons and ions is,the higher the density of ions in this region,leading to a higher degree of ionization.Conversely,the higher the ratio,the lower the degree of ionization.

    According to the calculation results of the collision frequency ratio,intermediate levels(0–2)of ionization exist throughout the three different chambers.However,the area of intermediate ionization levels of the flat grid model is apparently higher than those of the convex and concave grid models.This means that the discharge chamber with the flat grid has the highest discharge efficiency,followed by the convex grid,and the discharge chamber with the concave grid is the worst.

    4.Test results and comparison

    In the development of the 30 cm diameter ion thruster,due to the actual thermal deformation of the grids,the processing difficulties and cost of the large area of the flat and concave grids,LIP finally selected the convex grid for the 30 cm diameter ion thruster with reference to NSTAR and XIPS-25 ion thrusters.The flat and concave grids have been used in a 10 cm diameter Kaufman ion thruster developed by LIP.Plasma density and electron temperature in the discharge chamber of the thruster with the convex grid have not been measured at present.Therefore,the accuracy of the calculated results of the convex grid discharge model is verified by the relationship between the ion beam current and plasma density in the upstream region of the grids,which is shown in equation(6).If the comparison error is small,it is considered that the calculation results of the discharge models with the flat and concave grids are also accurate,and the three models can reflect the effects of different grid structures on plasma characteristics in the discharge chamber:

    whereIbis the ion beam,andAsandTsare the grid area and optical transparency(which is 67%),respectively.According to the calculation results of the electron temperature and plasma density in the upstream region of the grid shown in figure 7,the plasma density and the electron temperature from the grid edge to the center are in the range of 3.1×1016–6.3×1017m?3and 3.3–3.5 eV,respectively.The ion beam current calculated by equation(6)is within the range of 0.032–0.276 A,and the distribution along the radial length of the grid is shown in figure 10(c).Based on the plasma density distribution in figure 7(a),the total current on the screen grid is calculated to be 1.72 A and,considering that the simulation model in figure 2(b)is half of the real size of the discharge chamber,then the total beam current of the thruster with the convex grid is 3.44 A.The calculated value of the beam is obviously lower than the actual value of 3.67 A shown in table 1,because the beam current actually contains parts of doubly charged ions.However,in the simulation model,all the ions are set as singly charged and the existence of doubly charged ions is ignored.If the proportion of doubly charged ions in the total beam is assumed to be 10%,the actual singly charged ion current in the beam is 3.3 A,and the simulation results is essentially consistent with this value.In addition,the parameter settings of the model are one of the main reasons for the error.

    The performance of the 30 cm diameter ion thruster with the convex grid is tested on the specific vacuum facility of electric propulsion at LIP[22];the diameter of the facility is 3 m,the length is 5 m and the working vacuum degree of the facility reaches 2×10?5Pa.The ion beam density test is measured by a Faraday probe array,which is shown in figure 10(a).The Faraday probe array is mounted on a circular support frame with a diameter of 1 m,and the distance between the frame and the outlet of the thruster is 0.5 m.The support frame has four probe mounting rods of equal length across the center of the circular frame,and each mounting rod has 20 Faraday probes.The current received by the Faraday probes is transmitted through signal lines to a computer for processing.

    Since the Faraday probe array measures the received beam current,and the simulation results represent the beam value on the grid surface,the beam current density on the grid therefore needs to be equivalent to the probe array.As shown in figure 10(b),in the beam extraction process,assuming that the ion flux does not decay,then the current value does not change:mainly,the receiving area changes.The comparison between the calculation results via equivalence and measured beam density on the probe array is shown in figure 10(c).The measurements show that the beam density is the highest in the ring region in the middle of the grid,not in the center,while the beam density is lower the closer it is to the edge.There are some errors between the calculations and the measurements,as shown in figure 10(c),especially in the ring region from the middle of the grid to the edge.The analysis shows that the gap of the convex grids will be greatly reduced under the effect of thermal stress during thruster operation.Based on previous measurements[6],the gap will be decreased from 1 to 0.5 mm.Therefore,the beam density obtained by the test is measured after the thermal deformation of the ion optical system,while the focusing characteristics of the ion beam have changed significantly,which is the main reason for the existence of errors.In addition,the thermal deformation of the screen grid leads to the variation of the electric field isopotential line near the grids,which leads to the deviation between the distribution of plasma density of the upstream region of the grids and the real situation,and causes the errors between the calculations and measurements shown in figure 10(c).

    5.Conclusion

    The grid structure has obvious effects on the discharge characteristics of the ion thruster.According to the simulation results,the discharge chamber with the convex grid has a larger‘magnetic-field free area’and ionization region,and the distributions of collision frequency,electron temperature and plasma density are more uniform among the flat,convex and concave grids.In addition,according to the previous grid gap test results of the convex grid with 30 cm diameter,the variation of the grid hot gap is within an acceptable range(the minimum gap is greater than 0.25 mm).Therefore,the 30 cm diameter ion thruster finally chooses the convex grid as the optical system for ion acceleration and extraction.The comparison results of the ion beam show that there are some errors between the simulation results of 3.44 A and the actual value of 3.67 A because the existence of doubly charged ions is ignored and because of the parameter settings of the models.In addition,the errors between the simulation and measurement of the beam density are mainly due to the model,which does not consider the variation of the potential distribution near the grid and ion beam focusing caused by the grid thermal deformation during thruster operation.However,the analysis results show that,based on reasonable boundary condition settings and convergence of the calculation results,COMSOL can be used to obtain the discharge characteristics of the thruster and to evaluate its work performance.Meanwhile,the simulation results show that electron temperature distributions of the convex and the concave grids are more uniform.The flat grid has the highest discharge efficiency,followed by the convex grid,and the concave grid is the least efficient.This is because the convex grid gap is gradually reduced,while the concave grid gap is gradually increased during thruster operation,and the flat grid with a small area can obtain better beam flatness and its thermal deformation is small.Therefore,flat and concave grids are more suitable for small-diameter ion thrusters.At present,they have been applied in the 10 cm diameter Kaufman ion thruster developed by LIP.

    Acknowledgments

    This work is supported by National Natural Science Foundation of China(No.61901202),and by the Key Laboratory Funds for the Science and Technology on Vacuum Technology and Physics Laboratory,Lanzhou Institute of Physics(No.HTKJ2019KL510003).We are very grateful to the Science and Technology on Vacuum Technology and Physics Laboratory for providing the COMSOL software.

    猜你喜歡
    劉超趙勇
    時來運(yùn)轉(zhuǎn)
    老吳的拉面館
    Quantitative analysis and time-resolved characterization of simulated tokamak exhaust gas by laser-induced breakdown spectroscopy
    Theoretical study of novel B-C-O compounds with non-diamond isoelectronic
    Removal of GaN film over AlGaN with inductively coupled BCl3/Ar atomic layer etch
    征程萬里,初心如一
    雷鋒(2021年12期)2021-04-12 00:57:22
    Design and experiment of the centrifugal pump impellers with twisted inlet vice blades *
    趙勇:抓住脫貧突破口 實(shí)現(xiàn)老區(qū)又好又快發(fā)展
    Numerical Investigation on Drag Reduction Effect by Mass Injection from Porous Boundary Wall
    Ship hull optimization based on wave resistance using wavelet method*
    精品人妻1区二区| 99久久精品一区二区三区| 18美女黄网站色大片免费观看| 国产爱豆传媒在线观看| 欧美中文综合在线视频| 欧美成人一区二区免费高清观看 | 亚洲美女视频黄频| 亚洲在线观看片| 69av精品久久久久久| 久久久成人免费电影| 国产精品亚洲一级av第二区| 欧美黄色淫秽网站| 99热只有精品国产| 日本五十路高清| 色尼玛亚洲综合影院| 极品教师在线免费播放| 久久久成人免费电影| 动漫黄色视频在线观看| 99久久综合精品五月天人人| 国产精品女同一区二区软件 | 欧美激情久久久久久爽电影| 欧美丝袜亚洲另类 | www.999成人在线观看| 国产一区二区三区视频了| 夜夜躁狠狠躁天天躁| 欧美黑人巨大hd| 一区福利在线观看| 色老头精品视频在线观看| 国产aⅴ精品一区二区三区波| 国内毛片毛片毛片毛片毛片| ponron亚洲| 国产伦在线观看视频一区| 国产精品 欧美亚洲| 我要搜黄色片| 国产一级毛片七仙女欲春2| 精品一区二区三区av网在线观看| 18禁黄网站禁片午夜丰满| 国内精品一区二区在线观看| 亚洲成av人片免费观看| 午夜免费观看网址| 亚洲 欧美 日韩 在线 免费| 在线国产一区二区在线| 久久久国产成人免费| www.999成人在线观看| 久久久久亚洲av毛片大全| 久久国产精品影院| 99久久综合精品五月天人人| 黄频高清免费视频| 国产精品久久久人人做人人爽| 精品国产亚洲在线| 午夜成年电影在线免费观看| 琪琪午夜伦伦电影理论片6080| 亚洲国产精品999在线| 小说图片视频综合网站| 一卡2卡三卡四卡精品乱码亚洲| 久久久精品欧美日韩精品| 丁香欧美五月| 亚洲成人久久性| 亚洲在线观看片| 国产真实乱freesex| 久久久久久久精品吃奶| 国产欧美日韩一区二区三| 婷婷精品国产亚洲av| 亚洲欧洲精品一区二区精品久久久| 美女 人体艺术 gogo| а√天堂www在线а√下载| 人人妻,人人澡人人爽秒播| 色综合婷婷激情| 一级毛片精品| 男人和女人高潮做爰伦理| 中文字幕人妻丝袜一区二区| 国产高清视频在线播放一区| 老熟妇乱子伦视频在线观看| 成人三级做爰电影| 亚洲欧美日韩卡通动漫| 国产真人三级小视频在线观看| 中文亚洲av片在线观看爽| 免费在线观看日本一区| 国产美女午夜福利| 嫩草影院精品99| 国产真人三级小视频在线观看| 不卡av一区二区三区| 特级一级黄色大片| 国产97色在线日韩免费| 两性夫妻黄色片| 国产精品精品国产色婷婷| 亚洲真实伦在线观看| 99国产精品一区二区三区| 国产1区2区3区精品| 国产美女午夜福利| 黄色视频,在线免费观看| 国产精品一及| 亚洲美女黄片视频| 欧美日韩瑟瑟在线播放| 精品99又大又爽又粗少妇毛片 | 日本精品一区二区三区蜜桃| 亚洲成人中文字幕在线播放| 身体一侧抽搐| 男女那种视频在线观看| 欧美乱妇无乱码| 制服丝袜大香蕉在线| 亚洲欧美日韩无卡精品| 国产精品1区2区在线观看.| 不卡av一区二区三区| 精品久久久久久成人av| 欧美日韩福利视频一区二区| 国产亚洲精品久久久com| 亚洲av成人av| 波多野结衣高清无吗| 日韩精品青青久久久久久| 国产av麻豆久久久久久久| 午夜激情欧美在线| 女人被狂操c到高潮| 日韩av在线大香蕉| 亚洲av电影在线进入| 宅男免费午夜| 每晚都被弄得嗷嗷叫到高潮| 午夜影院日韩av| 成人高潮视频无遮挡免费网站| 天堂网av新在线| 中文字幕久久专区| 女警被强在线播放| 在线观看美女被高潮喷水网站 | 香蕉国产在线看| 99国产综合亚洲精品| 国产不卡一卡二| 精品国产三级普通话版| 亚洲精品在线美女| 欧美高清成人免费视频www| 两性午夜刺激爽爽歪歪视频在线观看| 亚洲第一电影网av| 亚洲成a人片在线一区二区| 国产成人一区二区三区免费视频网站| cao死你这个sao货| 欧美zozozo另类| 亚洲av美国av| 久久欧美精品欧美久久欧美| 午夜福利18| 欧美不卡视频在线免费观看| 两个人看的免费小视频| 国产精品电影一区二区三区| 桃色一区二区三区在线观看| 1024手机看黄色片| 国产亚洲av嫩草精品影院| 黄色视频,在线免费观看| 最近视频中文字幕2019在线8| av国产免费在线观看| 精品国产超薄肉色丝袜足j| 亚洲国产日韩欧美精品在线观看 | 久久这里只有精品19| 国产亚洲精品久久久com| 亚洲精华国产精华精| 国产精品日韩av在线免费观看| 免费高清视频大片| 亚洲片人在线观看| ponron亚洲| 国产一区在线观看成人免费| 国产一区二区三区视频了| 在线十欧美十亚洲十日本专区| 日韩国内少妇激情av| 高潮久久久久久久久久久不卡| 最新美女视频免费是黄的| 国产69精品久久久久777片 | 日本熟妇午夜| 法律面前人人平等表现在哪些方面| 日本三级黄在线观看| 一进一出抽搐动态| 久久久久国内视频| 欧美日韩亚洲国产一区二区在线观看| 国产成+人综合+亚洲专区| www日本在线高清视频| 悠悠久久av| 欧美黑人欧美精品刺激| 免费看十八禁软件| 色综合欧美亚洲国产小说| 国产熟女xx| 又粗又爽又猛毛片免费看| 老汉色av国产亚洲站长工具| 免费看十八禁软件| 丰满人妻熟妇乱又伦精品不卡| 在线观看舔阴道视频| 精品国产乱子伦一区二区三区| 免费无遮挡裸体视频| bbb黄色大片| 国产精品女同一区二区软件 | 老熟妇仑乱视频hdxx| 香蕉久久夜色| 一级毛片高清免费大全| 欧美中文综合在线视频| 成人一区二区视频在线观看| 欧美丝袜亚洲另类 | 亚洲成人中文字幕在线播放| 我的老师免费观看完整版| 村上凉子中文字幕在线| 精品国产超薄肉色丝袜足j| 日韩精品中文字幕看吧| 麻豆成人av在线观看| 亚洲国产欧美网| 一个人看视频在线观看www免费 | 日韩欧美国产在线观看| 看黄色毛片网站| 国产 一区 欧美 日韩| 在线国产一区二区在线| 久久久久久久久免费视频了| 一本精品99久久精品77| 一级毛片女人18水好多| av黄色大香蕉| 亚洲国产精品合色在线| 亚洲专区字幕在线| 国产视频一区二区在线看| 国产精品乱码一区二三区的特点| 亚洲无线观看免费| 欧美日韩精品网址| 午夜日韩欧美国产| 我的老师免费观看完整版| 又紧又爽又黄一区二区| 欧美日本亚洲视频在线播放| 大型黄色视频在线免费观看| 男人舔女人的私密视频| 亚洲精品在线观看二区| 日韩人妻高清精品专区| 国产精品免费一区二区三区在线| 美女午夜性视频免费| 国产精品免费一区二区三区在线| 亚洲片人在线观看| 黄色女人牲交| 亚洲av美国av| www.999成人在线观看| 亚洲性夜色夜夜综合| 亚洲国产欧洲综合997久久,| 母亲3免费完整高清在线观看| 精品国产乱码久久久久久男人| 久久性视频一级片| 18禁裸乳无遮挡免费网站照片| 桃色一区二区三区在线观看| 在线看三级毛片| 成人鲁丝片一二三区免费| 老司机福利观看| 成熟少妇高潮喷水视频| 91在线观看av| 一个人看视频在线观看www免费 | 国产精品国产高清国产av| 久久久久久人人人人人| 全区人妻精品视频| 色在线成人网| 国产成人系列免费观看| 黄色女人牲交| 精品乱码久久久久久99久播| 亚洲国产欧洲综合997久久,| 欧美性猛交╳xxx乱大交人| av在线天堂中文字幕| 高清毛片免费观看视频网站| 这个男人来自地球电影免费观看| 国产v大片淫在线免费观看| 欧美丝袜亚洲另类 | 精品99又大又爽又粗少妇毛片 | 亚洲欧洲精品一区二区精品久久久| 国产精品 欧美亚洲| 少妇熟女aⅴ在线视频| 757午夜福利合集在线观看| 悠悠久久av| 大型黄色视频在线免费观看| 免费高清视频大片| 99在线人妻在线中文字幕| 给我免费播放毛片高清在线观看| 久久久久九九精品影院| 日本一本二区三区精品| 国产探花在线观看一区二区| 老汉色∧v一级毛片| av在线蜜桃| 五月玫瑰六月丁香| 久久香蕉国产精品| 黄片大片在线免费观看| 好男人在线观看高清免费视频| 最新美女视频免费是黄的| 亚洲在线观看片| 1000部很黄的大片| 好男人电影高清在线观看| 国产欧美日韩精品亚洲av| 成年女人永久免费观看视频| 香蕉国产在线看| 女警被强在线播放| 午夜日韩欧美国产| 免费在线观看视频国产中文字幕亚洲| 精品人妻1区二区| 巨乳人妻的诱惑在线观看| av在线天堂中文字幕| av视频在线观看入口| 欧美中文综合在线视频| 精品国内亚洲2022精品成人| 岛国在线观看网站| 又紧又爽又黄一区二区| 精品国产美女av久久久久小说| 黄片小视频在线播放| 亚洲成av人片免费观看| 国产成人aa在线观看| 国产成人精品久久二区二区免费| 亚洲无线观看免费| 美女免费视频网站| 亚洲激情在线av| 亚洲五月天丁香| 精品一区二区三区视频在线 | 国产淫片久久久久久久久 | 国产精品国产高清国产av| 国产亚洲精品久久久久久毛片| 九九久久精品国产亚洲av麻豆 | 在线观看日韩欧美| 老司机午夜十八禁免费视频| 一本精品99久久精品77| 看片在线看免费视频| 91字幕亚洲| 国产真人三级小视频在线观看| 欧美成人一区二区免费高清观看 | 一级黄色大片毛片| 久久久久久九九精品二区国产| 国产高清激情床上av| 哪里可以看免费的av片| www.www免费av| 热99re8久久精品国产| 99热这里只有是精品50| 亚洲一区二区三区不卡视频| 一进一出抽搐gif免费好疼| 亚洲国产精品999在线| 哪里可以看免费的av片| 最近最新中文字幕大全免费视频| 在线观看免费午夜福利视频| 校园春色视频在线观看| 窝窝影院91人妻| 一级黄色大片毛片| 国产精品影院久久| 久久中文看片网| 黄色视频,在线免费观看| 国产激情久久老熟女| 国产精品电影一区二区三区| 伊人久久大香线蕉亚洲五| 国产不卡一卡二| 欧美xxxx黑人xx丫x性爽| 三级毛片av免费| 国产在线精品亚洲第一网站| 成人永久免费在线观看视频| 国产一区二区三区视频了| 亚洲成人中文字幕在线播放| 久久香蕉国产精品| 色精品久久人妻99蜜桃| 美女高潮的动态| 亚洲 欧美一区二区三区| 三级男女做爰猛烈吃奶摸视频| 男女床上黄色一级片免费看| 在线免费观看不下载黄p国产 | 久久草成人影院| 国产精品国产高清国产av| 日本免费一区二区三区高清不卡| 18禁美女被吸乳视频| 99在线视频只有这里精品首页| 久久久久久久久免费视频了| 成人高潮视频无遮挡免费网站| 黄色 视频免费看| 精品乱码久久久久久99久播| 男女之事视频高清在线观看| 热99re8久久精品国产| 亚洲天堂国产精品一区在线| 亚洲 欧美 日韩 在线 免费| av国产免费在线观看| 十八禁人妻一区二区| 女人被狂操c到高潮| 手机成人av网站| 日本与韩国留学比较| 麻豆国产97在线/欧美| 亚洲国产欧美一区二区综合| 两人在一起打扑克的视频| 99精品在免费线老司机午夜| 久久国产精品影院| 亚洲激情在线av| 欧美成人性av电影在线观看| 淫妇啪啪啪对白视频| 欧美激情久久久久久爽电影| 一a级毛片在线观看| 亚洲九九香蕉| 精品一区二区三区视频在线 | 国产精品永久免费网站| 99国产精品一区二区三区| 欧美又色又爽又黄视频| 一二三四社区在线视频社区8| 欧美黑人巨大hd| 精品久久久久久成人av| 国产精品亚洲美女久久久| 亚洲av成人不卡在线观看播放网| 国产91精品成人一区二区三区| 亚洲精品中文字幕一二三四区| 国产成人av教育| 两性午夜刺激爽爽歪歪视频在线观看| 日韩人妻高清精品专区| av视频在线观看入口| 午夜激情福利司机影院| 亚洲av成人精品一区久久| 国产99白浆流出| 国产私拍福利视频在线观看| 欧美色欧美亚洲另类二区| 日日夜夜操网爽| 久久久久亚洲av毛片大全| 久久久国产精品麻豆| 国产淫片久久久久久久久 | 欧美乱妇无乱码| 欧美大码av| 成年女人看的毛片在线观看| 久久精品国产亚洲av香蕉五月| 久久九九热精品免费| 中文字幕最新亚洲高清| 18禁观看日本| 久久精品综合一区二区三区| 后天国语完整版免费观看| 怎么达到女性高潮| 精品一区二区三区视频在线观看免费| 亚洲国产欧美人成| 久久午夜亚洲精品久久| 香蕉久久夜色| av天堂在线播放| 国产精品女同一区二区软件 | av女优亚洲男人天堂 | 久久亚洲精品不卡| 欧美激情在线99| 日韩大尺度精品在线看网址| 亚洲中文字幕日韩| 亚洲av电影在线进入| 久久草成人影院| 日本熟妇午夜| 国产一区二区激情短视频| 给我免费播放毛片高清在线观看| 啦啦啦观看免费观看视频高清| 99久久久亚洲精品蜜臀av| 国产精品一区二区三区四区久久| 操出白浆在线播放| 窝窝影院91人妻| 国产伦精品一区二区三区视频9 | 男女做爰动态图高潮gif福利片| 免费人成视频x8x8入口观看| 欧美最黄视频在线播放免费| 我的老师免费观看完整版| 老汉色av国产亚洲站长工具| 美女 人体艺术 gogo| 国产亚洲精品综合一区在线观看| 宅男免费午夜| 草草在线视频免费看| 日本一二三区视频观看| 真实男女啪啪啪动态图| 久久久国产欧美日韩av| 亚洲七黄色美女视频| 精品99又大又爽又粗少妇毛片 | www.自偷自拍.com| 99riav亚洲国产免费| 淫妇啪啪啪对白视频| 国产精品综合久久久久久久免费| 中亚洲国语对白在线视频| 日本在线视频免费播放| 国产精品久久久久久久电影 | 亚洲激情在线av| 国产欧美日韩一区二区精品| 日日夜夜操网爽| 国产精品一及| 精品欧美国产一区二区三| 国产欧美日韩精品一区二区| 色在线成人网| 日本免费a在线| 91老司机精品| 亚洲 国产 在线| 日本免费a在线| 最新美女视频免费是黄的| 1000部很黄的大片| 国产精品久久视频播放| 欧美一级a爱片免费观看看| 成人国产一区最新在线观看| 高清毛片免费观看视频网站| 免费在线观看日本一区| 欧美又色又爽又黄视频| 欧美成人性av电影在线观看| 在线播放国产精品三级| www日本黄色视频网| 男女视频在线观看网站免费| 国产精品一区二区免费欧美| 国产男靠女视频免费网站| 国产成人精品无人区| 人人妻人人看人人澡| 美女午夜性视频免费| 俺也久久电影网| 白带黄色成豆腐渣| 国产成年人精品一区二区| 亚洲人成伊人成综合网2020| 又紧又爽又黄一区二区| 国产一区二区在线av高清观看| 校园春色视频在线观看| 亚洲一区二区三区色噜噜| 国产精品自产拍在线观看55亚洲| 99久久成人亚洲精品观看| 日韩欧美 国产精品| 欧美zozozo另类| 99精品在免费线老司机午夜| 成人三级黄色视频| 非洲黑人性xxxx精品又粗又长| 天天一区二区日本电影三级| 色吧在线观看| 欧美中文日本在线观看视频| 三级男女做爰猛烈吃奶摸视频| 嫩草影院精品99| 亚洲天堂国产精品一区在线| 伊人久久大香线蕉亚洲五| 两性夫妻黄色片| 国产亚洲精品一区二区www| 日韩欧美一区二区三区在线观看| 老司机午夜福利在线观看视频| 国产亚洲精品久久久久久毛片| 九九热线精品视视频播放| 国产综合懂色| 午夜影院日韩av| 国产亚洲av嫩草精品影院| 不卡av一区二区三区| 亚洲熟妇中文字幕五十中出| 欧美中文综合在线视频| 久久精品人妻少妇| 日本撒尿小便嘘嘘汇集6| 麻豆久久精品国产亚洲av| 亚洲av第一区精品v没综合| 久久久久久久久免费视频了| 香蕉av资源在线| 欧美丝袜亚洲另类 | 欧美日韩精品网址| 久久中文字幕人妻熟女| 日韩欧美国产在线观看| 国产野战对白在线观看| 久9热在线精品视频| 成人国产综合亚洲| 综合色av麻豆| 男女那种视频在线观看| 国产成人系列免费观看| 在线免费观看不下载黄p国产 | 在线观看免费视频日本深夜| 最近视频中文字幕2019在线8| 国产精品野战在线观看| 男女做爰动态图高潮gif福利片| 久久这里只有精品19| 成人av一区二区三区在线看| 欧美在线一区亚洲| 亚洲av片天天在线观看| 亚洲,欧美精品.| 亚洲精品456在线播放app | 人人妻人人澡欧美一区二区| 给我免费播放毛片高清在线观看| 亚洲欧美日韩东京热| 男女之事视频高清在线观看| netflix在线观看网站| 美女高潮喷水抽搐中文字幕| 成人亚洲精品av一区二区| 男人舔奶头视频| 毛片女人毛片| 国产精品1区2区在线观看.| 欧美激情在线99| 欧美日本视频| 国产精品一区二区三区四区久久| 黑人欧美特级aaaaaa片| 美女免费视频网站| 中文字幕精品亚洲无线码一区| www.自偷自拍.com| 精品欧美国产一区二区三| 国产精品98久久久久久宅男小说| 18美女黄网站色大片免费观看| 免费在线观看视频国产中文字幕亚洲| 国产成人福利小说| 丁香六月欧美| 国产成年人精品一区二区| 日本三级黄在线观看| 中文字幕人成人乱码亚洲影| 国产97色在线日韩免费| 精品电影一区二区在线| 精品欧美国产一区二区三| 在线观看免费视频日本深夜| 99国产精品一区二区蜜桃av| 欧美绝顶高潮抽搐喷水| 偷拍熟女少妇极品色| 波多野结衣巨乳人妻| 国产精品一区二区三区四区久久| 精品国产乱子伦一区二区三区| 久久久成人免费电影| 午夜福利成人在线免费观看| 国产一区二区在线观看日韩 | 视频区欧美日本亚洲| 精品国产乱子伦一区二区三区| 国产91精品成人一区二区三区| 亚洲人成网站高清观看| 日韩欧美一区二区三区在线观看| 午夜福利成人在线免费观看| 男人舔女人下体高潮全视频| 人妻久久中文字幕网| 亚洲国产精品sss在线观看| 亚洲欧美激情综合另类| 久久人人精品亚洲av| 在线看三级毛片| 久久欧美精品欧美久久欧美| 首页视频小说图片口味搜索| 亚洲真实伦在线观看| 非洲黑人性xxxx精品又粗又长| 97碰自拍视频| 两人在一起打扑克的视频| 午夜福利在线在线| 亚洲人成伊人成综合网2020| 日韩 欧美 亚洲 中文字幕| 久久伊人香网站| 精品国产美女av久久久久小说| 一本精品99久久精品77| 一进一出好大好爽视频| 亚洲欧洲精品一区二区精品久久久| 五月玫瑰六月丁香| 亚洲av成人精品一区久久| 亚洲精品色激情综合| 91老司机精品| or卡值多少钱| 欧美日本视频|