Study on the characteristics of non-Maxwellian magnetized sheath in Hall thruster acceleration region

2022-08-01 11:34LongCHEN陈龙YuhaoAN安宇豪ShaojuanSUN孙少娟PingDUAN段萍BoruiJIANG姜博瑞YehuiYANG杨叶慧andZuojunCUI崔作君
Plasma Science and Technology 2022年7期
关键词:陈龙

Long CHEN(陈龙),Yuhao AN(安宇豪),Shaojuan SUN(孙少娟),Ping DUAN(段萍),Borui JIANG(姜博瑞),Yehui YANG(杨叶慧)and Zuojun CUI(崔作君)

College of Science,Dalian Maritime University,Dalian 116026,People’s Republic of China

Abstract The secondary electron emission(SEE)and inclined magnetic field are typical features at the channel wall of the Hall thruster acceleration region(AR),and the characteristics of the magnetized sheath have a significant effect on the radial potential distribution,ion radial acceleration and wall erosion.In this work,the magnetohydrodynamics model is used to study the characteristics of the magnetized sheath with SEE in the AR of Hall thruster.The electrons are assumed to obey non-extensive distribution,the ions and secondary electrons are magnetized.Based on the Sagdeev potential,the modified Bohm criterion is derived,and the influences of the non-extensive parameter and magnetic field on the AR sheath structure and parameters are discussed.Results show that,with the decrease of the parameter q,the high-energy electron leads to an increase of the potential drop in the sheath,and the sheath thickness expands accordingly,the kinetic energy rises when ions reach the wall,which can aggravate the wall erosion.Increasing the magnetic field inclination angle in the AR of the Hall thruster,the Lorenz force along the x direction acting as a resistance decelerating ions becomes larger which can reduce the wall erosion,while the strength of magnetic field in the AR has little effect on Bohm criterion and wall potential.The propellant type also has a certain effect on the values of wall potential,secondary electron number density and sheath thickness.

Keywords:Hall thruster,magnetized sheath,non-extensive distribution,secondary electron emission

1.Introduction

Hall thruster is one of the most effective space propulsion devices due to its high specific impulse and excellent precision,which is widely used in satellite attitude control,position maintenance and orbit adjustment[1,2].The main structure of Hall thruster includes buffer cavity,anode,discharge channel,electromagnetic coil,magnetic circuit structure and hollow cathode outside the discharge channel,the discharge channel is formed by two coaxial insulated cylinders.The electric field generated in the channel is usually along the central axis which can accelerate ions to produce thrust and electromagnetic coils are located outside the discharge channel to generate a quasi-radial magnetic field[3–5].The orthogonal electromagnetic fields force electrons to perform E× B drift motion along the circumferential direction and collide with the propellants to generate plasma.Since the electron velocity is far greater than the ion’s,the wall contacting with the plasma always accumulates negative charges and forms the sheath,which is a non-electrically neutral region and its potential is negative to the plasma.The potential drop forms an electric field perpendicular to the wall,attracting ions to the wall while repelling electrons.A stable sheath is formed when the flow of electrons and ions reaches equilibrium[6,7].The electrons with high kinetic energy can induce SEE which can significantly change wall floating potential and sheath features[8].The SEE yield varies according to the wall material,for example,boron nitride(BN)is usually used as the wall material of Hall thruster and its SEE coefficient is usually in the range of 0.2–1[9].Besides,the magnetic field in the discharge channel constrains the electrical conductivity across magnetic field,so the sheath characteristics can be greatly altered[10].As the wall of AR is close to the thruster magnetic pole,the magnetic field strength in the sheath region around AR is relatively high(up to 0.15 T)and the magnetization effect of the ions cannot be ignored.

The kinetic simulations studied by Sydorenko[11]and Smirnov[12]have shown that the high-energy tail of the electron velocity distribution function(EVDF)in the thruster is almost exhausted due to the collision between the electrons and the wall;meanwhile,experimental studies[13,14]have proved that the electron energy distribution function in the thruster obviously deviates from the Maxwellian distribution due to the strong interaction between the electrons and the wall.Maxwellian distribution approximation has its limitations for systems such as long-range Coulomb interactions and non-thermodynamic equilibrium plasmas.However,many of the fluid simulations and mixed simulation on Hall thruster sheath still adopted the Maxwellian distribution of electrons as approximation[3,15–17].In 1988,Tsallis[18]introduced a statistical method to extend the Boltzmann entropy concept to nonextensive generalized entropy,and its entropy form is

wherekBis Boltzmann’s constant,N is the sum of micro states,piis the probability of a certain micro state,and q is the nonextensive degree of the system,which works as an indicator to reflect the system’s degree of deviation away from the Maxwellian equilibrium.The most special property of q entropy is its non-additive(namely non-extensive),which means that the sum entropy of each independent system is different from the entropy of the entire system.Therefore,the total entropy of the two independent subsystems satisfies the following relationship[19]

where I and J are two independent subsystems.Equation(2)shows that when q=1,the entropy of the composite system is equal to the sum entropy of the two independent subsystems,which means thatqentropy degenerates to the standard extended Boltzmann–Gibbs entropy;whenq< 1,the generalized entropy of the composite system is greater than the sum entropy of each component,which is called super-extension distribution;when q> 1,the generalized entropy of the composite system is less than the sum entropy of each component,which is called sub-extension distribution.

In recent years,the non-extensive distribution of electrons has been applied to study the fundamental problems of plasma sheath under different physical conditions[19–30].Dhawanet al[21]considered the non-extensive distribution of electron,ion temperature and ion neutral collision in the simulation.The ion enters the sheath region at a velocity lower than the typical Bohm velocity,which would change with the non-extensive parameter,and finally affect the sheath thickness.Zouet al[22]established a fluid model of collisional magnetized plasma sheath with nonextensive electrons and derived the modified Bohm criterion.The study has found that Bohm velocity decreased with the increase of parameter q,indicating that the physical structure of the collisional magnetized sheath with Maxwellian distribution electrons is significantly different from that with non-extensive electrons.Safaet al[23]studied the magnetized sheath structure with nonextensive distribution and Ghaniet al[24]studied the dusty plasma sheath structure with super-extensive electron distribution and dust particle SEE.They have found that the q value can affect the width of plasma sheath.Basnetet al[25]studied the characteristics of the magnetized plasma sheath in the presence of both Boltzmann and non-extensive distribution electrons using the fluid model.The results have shown that the parameter q and density ratio have significant effects on the parameters of the magnetized plasma sheath.

In conclusion,the non-Maxwellian electron distribution has a significant effect on the sheath characteristics,which should also be considered in the sheath simulation of Hall thruster.Furthermore,the magnetic field strength in the sheath around AR is about 0.15 T,which means the magnetized effects of ions and electrons cannot be ignored.In this work,assuming that electrons obey non-extensive distribution,the effects of non-extensive parameter and magnetic field on Bohm criterion,wall potential,spatial charge density,secondary electron density distribution and ion kinetic energy are studied by numerical simulation.In section 2,the magnetohydrodynamics equations used in this work which consider the magnetization effect,the self-consistent Bohm criterion and the wall floating potential are described.In section 3,the numerical results and discussions are presented.Finally,conclusions are given in section 4.

2.Physical model and basic equations

A magnetized plasma sheath model is established near the wall of the Hall thruster AR,as shown in figure 1.Here,x,y,z represent the radial,circumferential and axial directions of the Hall thruster channel,respectively.For simplification,the physical parameters of the sheath are assumed to be consistent along the circumferential and axial directions.Thus,the spatial coordinate can be reduced into one-dimension(radial direction).The sheath edge is at the position x=0,the wall is at x=xw,and the area between x and xwis the magnetized sheath region.The velocities of magnetized ions and magnetized electrons are three-dimensional in velocity space due to the existence of inclined magnetic field.The applied constant magnetic field is located in the(x,z)plane,and the angle between the magnetic field and the x-axis isθ,then the form of the inclined magnetic feild is:B=B(c os θex+sin θez),where B is the magnitude of magnetic field.At the interface(x=0)between the bulk plasma and the sheath,the electrostatic potential is set toφ=0.Besides ions and electrons,the sheath also contains magnetized secondary electrons generated by high-energy primary electrons impacting the wall.

Using the fluid model and assuming the electrons in the plasma sheath obey the non-extensive distribution of Tsallis statistical theory,the one-dimensional velocity distribution function of the electrons is[26]

whereCqis a normalization constant,and its form is

whereme,veandTeare mass,velocity and temperature of electrons in the sheath region,respectively.eis the elementary charge,kBis the Boltzmann constant,φis the electric potential.Γ is the standard gamma function,ne0is the electron number density at the sheath edge.Whenq=1,equation(3)reduces to the Maxwellian–Boltzmann velocity distribution function;whenq<-1,the non-extensive velocity distribution function cannot be normalized.Whenq> 1,the velocity distribution function is truncated,its form is[23,27]

According to the non-extensive EVDF,and integrating equation(3)in space,the electron density with non-extensive distribution can be obtained as follows:

whereneis the electron number density in the sheath region.

The magnetized ions are described by the fluid equations namely,the continuity equation and momentum transport equation,which are shown,in the steady state,in equations(7)and(8):

whereni,viandmiare number density,velocity vector and mass of ions,respectively.

The secondary electrons are generated by high-energy electrons impacting the wall and treated as fluid.Thus,the magnetized secondary electrons in the sheath satisfy the continuity and momentum equations,as follows:

wherens,msandTsare number density,mass and temperature of secondary electrons,vsxis the velocity of secondary electrons in thex-direction,andvsis the three-dimensional velocity vector of the secondary electrons.

Poisson’s equation that relates the charge density to sheath electrostatic potential is given as:

wheree0is the vacuum permittivity.

At the wall surface,the secondary electrons flow satisfies the following form:

wherejsis the secondary electron flux ejected from the wall andjeis the electron flux reaching the wall,γis the SEE coefficient.When the sheath is stable,the ion and electron current densities at the insulation wall are equivalent,that is:

In equation(13)

wherejiis the current density of ions,andφwis the floating potential at the wall.At the sheath edge,charged particle number density satisfies the plasma quasi-neutrality condition:ni0=ns0+ne0.

whereusxwis the initial velocity of the secondary electrons emitted from the wall in thex-direction,Φwis the normalized wall potential.Equations(27),(28),(30)and(31)show the relationship between the ion Bohm velocityuix0,the wall potential Φwand the secondary electron number densityδat the sheath edge,which form a self-consistent sheath model.By numerically solving these four equations,the evolving trends of ion Bohm velocity,wall potential and sheath edge secondary electron number density with non-extensive parameter and magnetic field inclination angle are studied.Taking these derived parameters as the initial conditions,the equations(17)–(26)can be solved numerically.Since it is a set of stiff ordinary differential equations,the‘ode15s’solver from MATLAB software is used to treat this problem,the specific numerical method used in the solver is the numerical differentiation formulas(NDFs)[31].After that,the physical parameters of sheath around AR wall,such as sheath potential and secondary electron number density distribution can be obtained.

3.Numerical results and discussion

The physical parameters used in this work are set as follows:Φ(ξ=0)=0,E0=0.01,n0=1 × 1017m-3,Te=10 eV,ts=0.01,γ=0.01 –0.8,B=0.04 –0.1 T,θ=0° – 80°,the axial velocity of ionsuiz0in the AR of Hall thruster is set to Mach number between 7 and 24(the corresponding ion velocity is about --

15000 45000 m s1for xenon propellant andT=

e 10 eV).Krypton(Kr)and xenon(Xe)are selected as the propellant gas to study the effect of ion mass on discharge.

3.1.The effect of non-extensive parameter on the AR sheath

Figure 2(a)shows the EVDF profile for different non-extensive parameterqaccording to equation(3).Whenq=1,the non-extensive distribution function reduces to the Maxwellian distribution.Whenq<1,the function profile becomes wider and the peak value is smaller,so the number of high-energy electrons in the system is relatively large.On the contrary,whenq>1,the EVDF shows a clear truncation in the high energy tail,this means in large-qcase,high velocity electrons extinct and the number of high-energy electrons in the system is relatively small,this feature could cause narrower sheath thickness due to smaller electron flux to the channel wall.Besides,the influence ofqon low velocity electron proportion is nearly linear.Figures 2(b)–(d)show the effects of nonextensive parameterqon the wall potential,secondary electron number density and Bohm velocity with different propellant gas,respectively.Figure 2(b)shows that the wall potential becomes larger with increasing the value of parameterq.When the value of the parameterqis relatively large,according to figure 2(a),there are less high-velocity electrons in the sheath,the electron flux towards the wall decreases relatively.Less negative charges accumulate in the wall,thus the wall potential is high.It also can be seen in figure 2(b)that the wall potential is lower when the propellant gas is Xe.The potential drop in the sheath region becomes larger for increasing mass number of ions,which also causes high energy ion bombardment on the wall.Figure 2(c)shows that the number density of secondary electrons reaching the sheath edge increases with the increase of the parameterq.This is due to the higher electron density at the wall whenqis relatively large,which induces high yield of SEE.In addition,as the propellant gas is Xe,the number density of secondary electrons reaching the sheath edge is less.Compared with Kr,when Xe is selected as the propellant gas,the wall potential is relatively lower,resulting in fewer electron flux reaching the wall,therefore the SEE also decreases in the condition of sameγ,secondary electron density also diminishes at the sheath edge.

Figure 1.Schematic diagram of the discharge channel and the magnetized plasma sheath at channel wall in AR of Hall thruster.

Figure 2.Sheath parameters versus parameter q.(a)Non-extensive electron velocity distribution,(b)wall potential,(c)secondary electrons density at the sheath edge,(d)Bohm velocity.

Figure 3.Sheath parameters versus parameter q.(a)Potential,(b)space charge density,(c)secondary electron density,(d)ion kinetic energy.

Figure 2(d)shows that Bohm velocity is lower at high values of parameterq.When the value ofqis large,the average velocity of electrons is relatively slow,resulting in the reduction of electron flux to the wall.Consequently,the ion flux at the wall also reduces to preserve the stability of the sheath.According to the continuity equation(7),the Bohm velocity of ions should also decrease.By comparing the two lines in figure 2(d),it can be seen that the ion mass number has little effect on the Bohm velocity.

Figure 4.Sheath parameters versus magnetic field inclination angle.(a)Wall potential,(b)secondary electron density at the sheath edge,(c)Bohm velocity.

Figure 3 shows the effect of non-extensive parameterqon several physical quantities in the sheath region when the propellant gas is Xe,two cases of super-extensive distribution(q=0.8 and 0.9)and two cases of sub-extensive distribution(q=1.1 and 1.2)are selected in the simulation.It is clearly seen in figure 3(a)that the potential in the sheath region falls more rapidly for large parameterq,which also causes larger electric field in the sheath.On the other hand,by decreasing the value ofq,the wall potential is reduced due to the accumulation of high-energy electrons on the wall,and the sheath potential drop becomes larger.

Figure 5.Sheath parameters versus magnetic field inclination angle.(a)Potential,(b)sheath thickness,(c)secondary electron density,(d)ion kinetic energy in x direction.

It is shown in figure 3(b)that the space charge density profile has a peak,the value of which declines as parameterqdecreases,and the position of which is closer to the wall with smallerq.According to figure 2(a),when the value ofqis relatively small,there are more high-energy electrons which can penetrate deep into the sheath region,and reduce the net space charge density.Figure 3(c)shows that the increase of the parameterqcan enhance the secondary electron number density in the sheath,which is consistent with figure 2(c).Figure 3(d)shows that the growth trend of ion kinetic energyEk(the dimensionless unit of energy iskBTe)transforms with the increase ofq.Whenqtakes a high value,the ion energy entering the sheath is relatively small due to the influence ofqon Bohm velocity.However,at a certain position in the sheath,the line sequence of ion energy is reversed owing to the different growth rates as shown in the enlarged figure window in figure 3(d).Moreover,when the value ofqis small,the sheath potential drop is large,which causes ions get more energy through the sheath.

3.2.The effect of magnetic field on the AR sheath

The magnetic field strength around the wall at thruster AR is relatively large,and its influence on the ions and secondary electrons cannot be ignored.Especially in the type of magnetic shielding Hall thruster,which has attracted extensive attention recently,magnetic lines gather at the outlet wall to form a high potential preventing ions from damaging the wall,in such sheath area,the magnetic field strength grows high,and the magnetization on ions and secondary electrons have significant effects on the physical properties of the sheath.

Figure 4 shows the effects of magnetic field inclination angle on wall potential,ion Bohm velocity and secondary electron number density at the sheath edge.Figure 4(a)shows the variation trend of wall potential with magnetic field angle under differentq.Whenq>1,the wall potential hardly changes with the increase of angle.Whenq<1,the increase of the inclination angle can reduce the wall potential.This is because the magnetic field angle affects the Lorentz force of ions in thex-direction,which hampers the movement of ions towards the wall.It is shown in figure 4(b)that the secondary electron number density at the sheath edge decreases with the increase of magnetic field inclination angle.As the inclination angle takes a high value,the ion density at the wall is less.In order to ensure the stability of the sheath at the wall,the electron density also reduces accordingly.Thus,the secondary electron density decreases under the sameγconditions.Moreover,when the inclination angle is small,qhas a great influence on the secondary electron density at the sheath edge.When the inclination angle is large,the influence ofqon the secondary electron density becomes weak.Figure 4(c)clearly shows that the value of Bohm velocity is lower at higher values of inclination angle.That is,the kinetic energy of ions entering the sheath is smaller.Furthermore,when the magnetic field angle takes a high value,Bohm velocity changes little with different values ofq.

Figure 5 shows the effect of different magnetic field inclination angle on the potential,space charge density,secondary electron density and positive ion kinetic energy in the sheath region.It is seen in figure 5(a)that as the inclination angle decreases,the space potential falls more rapidly and the sheath thickness becomes thinner.The effect of magnetic field inclination angle on wall potential is consistent with that shown in figure 4(a).In addition,the result shows that the inclination angle in the AR has a significant effect on the sheath thickness.Figure 5(b)shows the effect of angle on sheath thickness at different values ofq.It is found that the sheath thickness increases with the increase of inclination angle.Figure 5(c)shows that the secondary electron number density in the sheath is relatively less at high values of the inclination angle,which is due to the short sheath length and higher electron flux to the wall.As shown in figure 5(d),with the increase of the inclination angle,the kinetic energy of the ions in thex-direction decreases,and the kinetic energy reaching the wall is also smaller.The increase of the inclination angle makes the Lorentz force on the ion in thex-direction stronger,which reduces the ion velocity,and enhances the sheath length.Thus,the component of the magnetic field parallel to the wall plays the dominant role in the influence on the sheath characteristics.

Figure 6 shows the effect of magnetic field on sheath potential and ion kinetic energy.It is seen in figure 6(a)that as the magnitude of magnetic field enhances,the sheath wall potential and the sheath thickness increase.This is mainly because when the magnitude of magnetic field is small,it has little effect on the sheath potential.Figure 6(b)shows that the kinetic energy of ionsEkis small at high value of the magnetic field.When the magnitude of magnetic field increases,the Lorentz force on ions becomes larger,resulting in the slowdown of ions inx-direction,therefore the kinetic energy of ions reduces.Figure 6(c)shows that the secondary electron number density of the wall is relatively less at high magnitude of the magnetic field,which is mainly because the enhancement of electron magnetization leads to the reduction of electron flux to the wall.In addition,according to the equations(27),(28)and(30),the magnitude of magnetic field has no effect on the Bohm criterion.

3.3.Wall erosion analysis of acceleration region

In the radial direction of AR in Hall thruster,potential drop created by the sheath area plays a critical role in the service life of the Hall thruster.Ions elevate the kinetic energy in the sheath and pre-sheath,then bombard and erode the wall.Figure 7 shows the effects of non-extensive parameterqand magnetic field inclination angleθon the ion kinetic energy and the incident angle of ion impact on the wall when the propellant gas is Xe.It is seen in figure 7(a)that as the parameterqand the magnetic field inclination angle increases,the initial kinetic energy of ionsEk0decreases.By increasing the magnetic field inclination angle,the influence of parameterqbecomes weaker.It is also shown in figure 7(b)that the kinetic energy of the ions reaching the wallEkwincreases with the increase of the magnetic field inclination angle,and decreases with the increase of the parameterq.When the parameterqtakes a smaller value,the kinetic energy of the ions at the wall changes significantly with the magnetic field angle.Figure 7(c)shows that as the magnetic field inclination angle and parameterqincrease,the incident angle between the ion and the wall normal increases gradually,which indicates that the trajectory of ions is more parallel to the wall when parameterqorθincreases.

Figure 6.Sheath parameters versus magnetic field strength.(a)Potential,(b)ion kinetic energy,(c)secondary electron density.

Figure 7.Energy and the wall erosion versus parameter q and magnetic field inclination angle.(a)Ion initial kinetic energy,(b)ion kinetic energy at the wall,(c)the incident angle of the ions,(d)wall erosion rate.

Figure 8.Sheath thickness versus the axial velocity of ions.(a) B=0.05 T,(b) B=0.1 T.

The erosion rate at the wall of Hall thruster is directly proportional to the ion current density perpendicular to the wall and the material sputtering rate which depends on the ion incident energy and incident angle.Thus,the empirical formula for erosion rate can be expressed ase=ji⊥Y(K,θ)[32],whereji⊥is the ion current density normal to the wall andY(K,θ)is the sputtering yield of the wall material.θ is the average incident angle of ions,Kis the total impact energy of the ions at the wall.The sputtering yield can be further expressed asY(θ,K)=fθ(θ)fK(K),a function determined by incident angle and kinetic energy of ions.By applying the coefficients provided by[32],the effects of non-extensive parameterqand magnetic field inclination angle on erosion rate can be derived,as shown in figure 7(d).It can be seen that as non-extensive parameterqand magnetic field inclination angle increase,the erosion rate decreases gradually.When the magnetic field inclination angle is large and perpendicular to the axial direction of the discharge channel of the Hall thruster,the erosion rate is small,indicating that magnetic shielding configuration can reduce the erosion to a certain extent.

Figure 8 shows the effect of ion axial velocity on AR sheath thickness.The thickness of the sheath gradually expands as the ion axial velocity increases.As shown in figure 8(a),when the magnetic field is about 500 Gs,the sheath thickness has a linear relationship with the ion axial velocity.Once the ion velocity grows rapidly in the AR,the sheath thickness increases accordingly.This phenomenon can be interpreted as that when ion velocity grows,the component of Lorentz force on the ions along the radial direction becomes larger and prevents the ions from reaching the wall,thus the sheath becomes wider in order to maintain the sheath stability.As shown in figure 8(b),when the strength of magnetic field is larger,the sheath thickness changes more obviously.In addition,the atomic mass of the propellant also has a certain influence on the sheath thickness.With the decrease of propellant mass,the change of sheath thickness is clearer,especially when the strength of magnetic field is large,the sheath width with Kr propellant changes nonlinearly.

4.Conclusion

In this work,a 1D3V fluid model is established to study the characteristics of magnetized sheath with non-Maxwellian distribution elections in the AR of Hall thruster discharging channel.The modified Bohm criterion is deduced according to Sagdeev potential theory,which is self-consistent with the wall potential and SEE yield.It is found that the non-extensive distribution of electrons has a great impact on the potential,net charge,secondary electron number density,Bohm velocity and wall erosion rate.

With the increase of non-extensive parameterq,the number of high-energy electrons in the system is relatively small,the wall potential enhances,the secondary electron number density at the sheath edge increases,the value of the Bohm criterion decreases,and the number of ions entering the sheath enhances.When the magnetic field inclination angle in the sheath of AR is relatively large,the wall potential is low,the secondary electron number density at the sheath edge is less,the value of Bohm criterion is small.However,the sheath magnetic field strength in the AR of Hall thruster has no significant effect on Bohm criterion and wall potential.Moreover,the kind of propellant gas also has a certain effect on the values of wall potential and sheath edge secondary electron number density.

With the decrease of the parameterq,the high-energy electrons lead to an increase of the potential drop in the sheath region,the sheath thickness expands accordingly,the kinetic energy rises when the ions reach the wall,which can aggravate the wall erosion.When the magnetic field inclination angle decreases in the AR of the Hall thruster,the Lorenz force along the x-direction acting as a resistance on ions becomes smaller,the potential falls more rapidly,and the secondary electron number density shows an increasing trend.Under these circumstances,the incident angle between the ion and the wall normal is smaller and the wall erosion is more serious.As the parameterqand magnetic field angle increase,the wall erosion rate of the Hall thruster decreases.

Acknowledgments

The work is supported by National Natural Science Foundation of China(Nos.11975062,11605021,11975088)and the China Postdoctoral Science Foundation(No.2017M621120).

猜你喜欢
陈龙
情书
When weird weather strikes 当怪天气来临时
陈龙:“史上最帅武松”的幸福生活
茶,有点苦
章龄之 “选男友就选胡歌,要结婚就选陈龙”
准确审题正确列式精确验证
《机械工程测试技术》教学方法初探
女老板倒追“扫帚哥”,爱你就把你照亮
“餐馆女王”情定环卫工,我就在你左右
陈龙 慢节奏萌叔