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

    Investigation into the thermal effect of the LIPS-200 ion thruster plume

    2022-08-01 11:33:34XinweiCHEN陳新偉BijiaoHE賀碧蛟ZuoGU顧左HaiGENG耿海NingGUO郭寧YongZHAO趙勇KaiSHI史楷KaiTIAN田愷TaoCHEN陳燾andYifanMA馬伊帆
    Plasma Science and Technology 2022年7期
    關(guān)鍵詞:趙勇

    Xinwei CHEN(陳新偉),Bijiao HE(賀碧蛟),Zuo GU(顧左),Hai GENG(耿海),Ning GUO(郭寧),Yong ZHAO(趙勇),Kai SHI(史楷),Kai TIAN(田愷),Tao CHEN(陳燾) and Yifan MA(馬伊帆)

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

    2 School of Astronautics,Beihang University,Beijing 100191,People’s Republic of China

    Abstract The distribution of the thermal effects of the ion thruster plume are essential for estimating the influence of the thruster plume,improving the layout of the spacecraft,and for the thermal shielding of critical sensitive components.In order to obtain the heat flow distribution in the plume of the LIPS-200 xenon ion thruster,an experimental study of the thermal effects of the plume has been conducted in this work,with a total heat flow sensor and a radiant heat flow sensor over an axial distance of 0.5–0.9 m and a thruster angle of 0°–60°.Combined with a Faraday probe and a retarding potential analyzer,the thermal accommodation coefficient of the sensor surface in the plume is available.The results of the experiment show that the xenon ion thruster plume heat flow is mainly concentrated within a range of 15°.The total and radial heat flow of the plume downstream of the thruster gradually decreases along the axial and radial directions,with the corresponding values of 11.78 kW m?2 and 0.3 kW m?2 for the axial 0.5 m position,respectively.At the same position,the radiation heat flow accounts for a very small part of the total heat flow,approximately 3%–5%.The thermal accommodation factor is 0.72–0.99 over the measured region.Furthermore,the PIC and DSMC methods based on the Maxwell thermal accommodation coefficient model(EX-PWS)show a maximum error of 28.6% between simulation and experiment for LIPS-200 ion thruster plume heat flow,which,on the one hand,provides an experimental basis for studying the interaction between the ion thruster and the spacecraft,and on the other hand provides optimization of the ion thruster plume simulation model.

    Keywords:LIPS-200 ion thruster,plume thermal effects,thermal accommodation factor,PIC-DSMC

    1.Introduction

    Ion thrusters[1–3]have the potential to dramatically enhance the payload,prolong the service life and reduce the launching budget of spacecraft owing to their high specific impulse,low thrust,prolonged lifespan,and repeatable start-up etc,and have been successfully employed in missions concerning attitude control,station keeping,drag compensation,and deep space exploration.The implementation of electric propulsion technology on spacecraft will alter the plasma environment in the vicinity of the spacecraft owing to its special plume environment,resulting in a series of effects such as charged spacecraft surfaces,electromagnetic interference,and subsequent changes in the thermal properties of thermally sensitive surfaces.The ion thruster thermal effect mainly refers to the heat exchange resulting from the collision of beam energetic ions against spacecraft surfaces,which places the spacecraft in a thermal environment,causing an increase in spacecraft surface temperature and material damage[4,5].This can cause serious degradation of the working quality of the spacecraft and,in severe cases,can lead to the overall inability of the spacecraft to function and the failure of the space mission.The diagnostic measurements of electric propulsion are currently primarily conducted using Faraday probes[6],emission probes[7],Langmuir probes[8],retarding potential analyzer s(RPA)[9],and other contact diagnostic devices placed directly in the plasma plume which,in turn,can provide information such as electron density,spatial potential,ion energy distribution,plasma potential,the current density in the plasma,etc.

    Figure 1.Sketch(a)and photograph(b)of the ion thruster plume thermal effect measurement system.

    The thermal effect of the electric thruster plume has also been studied by a minority of scholars.In 1998,Kinget al,at the University of Michigan,used a heat flow meter to measure the heat flow of the SPT-100 Hall thruster[10]in three-degree increments at the radii of 0.5 m and 1 m,respectively.The results show that the high-energy plume region is the main distribution area of the heat flow.Also,due to the interaction of the plume with the background particles,many active neutral particles are found during the exchange of low-energy ions.In 2008,Berisfordet al,at the University of Texas in Austin,carried out a diagnostic experimental study of the thermal effects of the spiral helicon plasma discharge process using an infrared camera[11].The experimental results show that the energy loss in the insulated tube increases continuously with increasing electric field strength.However,there is very little formal published literature on the experimental diagnosis of thermal effects of the ion thruster in the plasma.

    In this work,in order to investigate the specific distribution of thermal effects in the plume of the ion thruster,an experimental and numerical simulation study of the thermal effects in the plume field of a 200 mm xenon ion thruster has been carried out using a heat flow sensor diagnostic system.Combined with the ion current density determined by the faraday probe,and the ion energy of the plasma measured by the RPA probe,the thermal accommodation coefficient of the xenon plasma to the surface of the heat flow sensor is then derived from the data results that have been obtained for the total and radiated heat flow of the ion thruster plume.Furthermore,a simulation study of the thermal effect of the thruster plume was carried out with the EX-PWS software based on the Maxwell model of the thermal accommodation factor.The distribution of heat flow in the plume field of a xenon ion thruster is obtained,which,on the one hand,provides an experimental basis for studying the interaction between the ion thruster and the spacecraft,and on the other hand,provides optimization of the ion thruster plume simulation model.

    The remainder of this paper is organized as follows.Section 2 characterizes the experimental apparatus and methodology,section 3 presents an analysis of the experimental results and the simulations of heat flux and thermal accommodation coefficient,and section 4 yields the conclusions of the current study.

    2.Experimental setup

    This experiment for the thermal effects of a 200 mm xenon ion thruster(LIPS-200)plume was carried out in the vacuum plume experimental system(PES)at the Beihang University Shahe campus.The entire experimental system was composed of five parts:the vacuum plume experimental system(PES),the 200 mm xenon ion thruster(LIPS-200),the heat flow sensor(HFS),the Faraday probe(FP),and the retarding potential analyzer(RPA)probe.The sketch(a)and photograph(b)of the ion thruster plume thermal effect measurement system are shown in figure 1.

    As shown in figure 1,the HFS,FP,and RPA are placed on the movement mechanism support in a horizontal position.The center of rotation of the movement mechanism is the center of the ion thruster exit plane.The maximum travel of the movement mechanism is 1 m@360°.The total heat flow,radiation heat flow,ion current density,and ion energy distribution were available for different radii and angle points in the plume.

    2.1.Vacuum systems

    The PES of Beihang University principally constitutes a vacuum chamber,several pumps,a pumping system consisting of a molecular pump,a liquid nitrogen heat sink,an electrical control system,and a storage and supply system.The vacuum chamber[12,13]features a cylindrical horizontal structure with a maximum outer diameter of 5.5 m and a total length of 12.8 m composed of stainless steel.The effective internal diameter of the chamber after installation of the heat sink is 4.2 m,its length is 9.9 m and the volume is 250 m3.The vacuum chamber is capable of pumping(nitrogen)superior to 107l s?1and has an ultimate vacuum of 5×10?5Pa when the double layer heat sink is cooled with liquid nitrogen.The operating vacuum pressure of the system during the experiments was better than 9×10?4Pa.The photo of the vacuum system is shown in figure 2.

    2.2.LIPS-200 ion thruster

    The thruster used in this plume thermal effects test was a 200 mm diameter xenon ion thruster developed by the Lanzhou Institute of Physics,in which xenon was employed as the propellant,with the xenon purity reaching 99.9995%and the vacuum maintained at around 2.3×10?3Pa during operation.The thruster was successfully used onboard the SJ-9A[14]satellite.The performance parameters of the xenon thruster with a diameter of 20 cm are approximate:a nominal thrust of 40 mN,a specific impulse of >3000 s,and a power of 1000 W.The condition parameters of the thruster during this experiment were mainly:a discharge chamber mass flow rate of 1.24 sccm,a beam current of 814 mA,an anode mass flow rate of 11.24 sccm,a beam voltage of 1026 V.Photos of the LIPS-200 ion thruster are shown in figure 1(b).

    2.3.HFS

    The HFS system consists of a total and radiation heat flow meter,a signal conversion device,and a current signal acquisition system.The Schmidt–Boelter heat flow sensor was used to obtain the heat flow,which uses a thermopile based on a spatial temperature gradient to measure the temperature difference and determine the amount of total heat flow.The radiant heat flow is obtained by placing a sapphire glass sheet in front of the receiving front surface of the HFS.A photo of the HFS is shown in figure 3.When the ion thruster is in operation the HFS converts the received plume heat flow signal into a voltage signal(mV),the voltage signal is then transformed inside the chamber into a current signal(mA),this current signal is transmitted by means of a flange directly to the current acquisition system outside the chamber and finally to the computer terminal,where the total heat flow and radiation heat flow of the ion thruster plume are obtained.The relative positions of the heat flow sensor and thruster during the experiment are shown in figure 1(a).The range of the HFS is 0–10 kW m?2,with a maximum measurement of 12 kW m?2and an error of ±2% of the measured value.

    2.4.FP

    The planar FP[5,6,15,16]is a type of diagnostic apparatus used to identify the density of the ion current.Normally,it is composed of a planar collecting electrode(PCE)encompassed within a shield ring(SR).The collector surface of the FP used in this experiment is a molybdenum disk with a diameter of 10 mm.The material of the guard ring is stainless steel.The distance between the collector and the guard ring is 2 mm.The collector electrode and guard ring were supplied with bias voltages of ?30 V.The FP is placed on a vacuum movement mechanism.The electrical schematic and a photo of the FP are shown in figure 4.The relative positions of the FP and ion thruster during the experiment are shown in figure 1.

    The ion current density(j)is derived in equation(1).

    In equation(1),Vis the sampling resistor voltage,Ris the sampling resistor impedance,andArepresents the collection area.

    2.5.RPA

    The RPA[17,18]is principally employed to measure the energy distribution of ions in a plasma.The electrical schematic and a photo of the RPA probe are shown in figure 5.The relative positions of the RPA and ion thruster during the experiment are shown in figure 1.The RPA probe used in this experiment mainly consists of a housing,a four-layer grid,and a collection electrode,etc.The housing of the RPA consists of an inner diameter of 16 mm and a length of 22 mm stainless steel.The first electrode,the nearest plasma grid,is floating and is mainly used to reduce the disturbance caused by the plasma to the probe.The second grid has a voltage of?30 V,which aims to repel the electrons in the plasma and facilitate the entry of ions as much as possible.The third grid layer is scanned with a voltage of from 0 to 1200 V.This is mainly to screen the ions,with the aim of collecting the ion current bypassing only the higher energy ions.The fourth grid is applied with a voltage of ?30 V.The collection electrode consists of a 10 mm diameter molybdenum sheet.The ions from the plasma are filtered through the multi-layer grids,causing the higher energy ions to reach the collector surface of the RPA,where it is collected to obtain a current signal.The scanned voltage and the collected current constitute theI-Vcurve from which the ion energy in the plasma can be further obtained and its value is the magnitude of the firstorder derivative of the curve:

    Figure 3.Photo and principle diagram of the HFS.

    Figure 4.Electrical schematic(a)and photo of FP(b).

    Figure 5.Electrical schematic(a)and photo of the RPA probe(b).

    whereIandVdenote the current of the collection electrode and scanned voltage,respectively.mi,qeandedenote ion mass,species charge,and single-electron charge,respectively.Acis the area of the collection electrode andf(V)denotes a function related to the scanning voltage.

    Figure 6.Schematic diagram of the composition of the heat flow from the ion thruster plume to the target surface.

    3.Results and discussion

    3.1.Experimental results and discussions

    In conventional thermal effects measurements,heat flow sensors are often employed to detect both the energy transfer from the collision of plume particles with the sensor surface and the radiated energy transfer from the thruster and its plume to the sensor.The schematic diagram of the composition of the heat flow from the ion thruster plume to the target surface is shown in figure 6.

    The composition of the heat flow in the ion thruster plume can be expressed as:

    whereqtotalrepresents the total heat flow,which can be measured directly.Theqplumerepresents the heat flow generated by the action of the plasma plume particles on the surface.qradis the heat generated by the thruster and its plume radiation on the surface,which can be measured by adding sapphire glass to the head of the heat flow sensor.So,

    In the actual action,qplumecan be expressed as the difference between the incoming energy(qi)density and the energy of elastic reflection reflected(qr)energy density.

    To further investigate the incoming and reflected energy densities,we use the classical model of elastic collisions between gas molecules and walls—the Maxwell model.The Maxwell model assumes that the molecules are either fully energy-adapted diffuse reflection from a surface or mirror reflection with no energy exchange,as illustrated in figure 7.The thermal accommodation coefficient(α)is the fraction of diffuse molecules.Whenα=0 it is a complete mirror reflection and whenα=1 it is a complete diffuse reflection.

    Figure 7.Schematic diagram of Maxwell’s reflection model.

    According to the definition for the thermal accommodation coefficient(TAC)(α).

    Figure 8.Schematic layout of the heat flow measurement points.

    In summary,the thermal accommodation coefficient(α)in the Maxwell reflection model can be expressed as:

    Since the particle current parameterjof the incoming flow can be measured by the FP,Uis available by the RPA,qtotalandqradare determined by the total and radiant heat flow meters,the thermal adaptation coefficient of the xenon plasma interacting with the wall of the heat flow sensor is finally derived.The correlation between the energy flow density of the incoming flow and the thermal adaptation coefficient is established.

    3.1.1.HFS results.The experimental study of the heat flow distribution characteristics of the LIPS-200 thruster plume was carried out using a HFS.The principle of scanning from far-field to near-field,more points in the beam area than the non-beam area,and repeated measurements at crucial points were followed in the experiments.The measurement points of the probe are located at the horizontal plane consisting of the center of the ion thruster outlet plane.The specific distribution of measurement points:the distance from the thruster exit plane is 0.5 m,0.7 m,and 0.9 m respectively,and the range of 0°–60° angle with the thruster exit centerline.The center of the thruster exit plane is the center of the circle of the actual measurement point,the thruster exit centerline is 0°,and the counter-clockwise direction along the thruster exit centerline is the positive direction.The measurement points are arranged in 2.5°increments within the range of 0°–20°and angles in 5°increments within the range of 20°–60°.A schematic layout of the heat flow measurement points is shown in figure 8.

    Figure 9.Contour maps of total heat flow(a)and radiant heat flow(b)as a function of distance and angle.

    The magnitudes of the total and radiant heat flow densities in the plume field at different radii and angles are determined by the total HFS and the radiant heat flow meter respectively.The contour maps of total heat flow(a)and radiant heat flow(b)as a function of distance and angle are shown in figure 9.As can be seen from figure 9(a),the total heat flow at the centerline of the thruster at a radius of 0.5–0.9 m is about 6.8–12 kW m?2.At the identical angle and different radii,the smaller the radius,the larger the total heat flow;at the same radius,the total heat flow at the centerline of the thruster is the largest and decreases with the increasing angle,and then drops to the lowest at the angle of 30°,and does not change in the range of 30°–60°.Within the radius of 0.5–0.9 m and the angle of 0°–30°,the total heat flow at each radius accounts for most of the total heat flow at the respective radius,with the proportion being approximately 85%–95%.In the angle range of 30°–60°,the total heat flow at each radius ranges from roughly 0–1 kW m?2,accounting for approximately 5%–15% of the total heat flow.The analysis is attributed to the fact that the total heat flow obtained by the heat flow meter is mainly due to the ions impinging on the surface of the heat flow sensor,the kinetic energy of the ions is transformed into heat energy.Therefore,in the ion thruster beam region plasma beam density is great,the plasma velocity is high,the further away from the axis position,the plasma number density is relatively smaller as well as the ion velocity is also significantly decreased.Hence,the total heat flow density is greater in the range of 0°–30° and smaller in the range of 30°–60°.The ion action is basically in thermal motion and the heat flow density tends to be constant.

    As can be seen from figure 9(b),the distribution of radiation heat flow of the thruster is similar to that of the total heat flow.The radiation heat flow at the radius of 0.5–0.9 m of the thruster centerline is approximately in the range of 0.1–0.35 kW m?2.Under the same radius,the radiation heat flow at the centerline of the thruster is the maximum,with the increasing angle,the radiation heat flow decreases,the radiation heat flow decreases to the minimum at the angle of 30°.The radiation heat flow value essentially does not change in the range of 30°–60°.From the above analysis,it can be concluded that the total heat flow and radiation heat flow of the ion thruster plume accounts for a considerable part of the corresponding total in the range of 0°–30°,beyond which the total heat flow and radiation heat flow changes are essentially small.

    It is evident from figure 9 that at the same radius the radiant heat flow represents a minimal fraction of the total heat flow,around 3%–5%.At a radius of 0.5 m,the total and radiated heat flow is stable at around 30°,i.e.,at the junction of the beam and non-beam flows,which is generally consistent with the previous analysis.At a radius of 0.9 m,the radiated heat flow is stable in the region at an angle of about 15°.The possible reasons for this are the distance from the thruster,the smaller ion number density and thruster thermal radiation at this radius,and the larger sensor error.To better investigate and obtain the general law of the distribution of the thermal effect of the ion thruster plume,a nonlinear fitting of the total heat flow and the radiated heat flow in two variables,radius(r)and angle(θ)is carried out,and accordingly,the total(qtotal)and radiation(qrad)heat flow fitting equations are shown in equations(10)and(11),respectively.

    Figure 10.Diagram of raw FP data and fitted data for the position of the axis 500 mm downstream from the thruster.

    Figure 11.Contour maps of ion current density as a function of distance and angle.

    3.1.2.FP results.The raw FP data and fitted data for the position of the axis 500 mm downstream from the thruster are shown in figure 10.The value of the current density measured by the Faraday probe at 500 mm on the axis does not vary with time and is a constant value,which can be inferred from the ion thruster working very steadily during the measurement process.The raw data from the experiment were averaged to obtain the current at 500 mm and then this value was divided by the effective contact area of the FP to give a final current density of 1.2102 mA cm?2at 500 mm.Similarly,combined with the moving mechanism,the current density of the Faraday probe can be determined for a radius of 500–900 mm and an axial angle of 0°–15° to the thruster,as shown in figure 11.The LIPS-200 thruster beam divergence half-angle is less than 15°[19,20].From the heat flow probe data,it can be seen that the heat flow in the plume is mainly concentrated in the beam region(within 15°)and the non-beam region accounts for about 5%,so only the ion current density in the beam region of the thruster is presented here,as shown in figure 11.As can be seen from figure 11,the ion current density in the thruster plume exhibits a gradually decreasing characteristic along with the thruster axial and radial directions.For example,the ion current density decreases from 1.21 mA cm?2to 0.59 mA cm?2from 0.5 m to 0.9 m along the thruster axis.This data is used as data input for the calculation of the thermal accommodation factor.

    Figure 12.Schematic of the RPA probe I-V and differential curves located 500 mm from the thruster exit plane.

    3.1.3.RPA results.The schematic of the RPA probeI-Vand differential curves located 500 mm from the thruster exit plane is shown in figure 12.From the preceding analysis,it is noted that the peak of the differentiatedI-Vcurve corresponds to the ion energy of the RPA at that position of the plasma.It can be seen from figure 12 that the voltages in the energy distribution region are approximately concentrated in the range of 940–1012 eV,which corresponds to a peak value of approximately 968.5 eV,the magnitude of which is the ion potential of the plasma at this location.The corresponding ion energies of the plasma at axial positions of 600 mm,700 mm,800 mm,and 900 mm are 987.6 eV,976.8 eV,982.7 eV,and 988.4 eV.Since the acceleration voltage of the ion thruster is 1024 V,it can be seen that the energy of the ion thruster plasma in the axial position is still very stable and does not essentially depend on the axial distance,which in turn leads to the basic stability of the plasma plume energy in the beam region and indirectly also reflects the excellent stability of the ion thruster performance.

    3.1.4.TAC results.As can be seen from the previous information,the thermal accommodation coefficient(TAC)of the HFS surface in a xenon plasma plume can be derived from the combined results of the total heat flow meter,the radiation heat flow meter,the FP,and the RPA probe in the identical spatial location.The thermal adaptation coefficients of the xenon plasma to the surface of the HFS at various locations in the beam region are then derived,as shown in figure 13.

    Figure 13.Schematic of TAC as a function of angle for 0.5 m(a),0.7 m(b),and 0.9 m(c).

    It can be seen from figure 13(a)that the TAC of the xenon plasma interacting with the surface of the HFS at 500 mm does not vary with the angle,and the average value of the TAC in the range of 0°–15° in the beam region of the ion thruster is approximately 0.93.The minimum position occurs at an angle of 15°,which is considered to be at the boundary between the beam and non-beam regions.The plasma is considerably diluted and the amount of plasma and the velocity of motion is quite small,so there is a certain error compared to the core region.Similarly,it can be seen from figures 13(b)and(c)that the average TAC for both is 0.88 and 0.84 at locations with a radius of 0.7 m and 0.9 m,respectively.The TAC of the surface of the xenon plasma near the axial position is relatively high when acting with the surface of the heat flow meter,basically between 0.9 and 1.The TAC at the corresponding positions was between 0.8 and 0.9 when deviated from the axial direction.The possible reasons for this discrepancy are that the distance from the thruster(the density and velocity of the plasma cause errors in the measurement equipment),the vacuum level of the experimental process,the variation of the surface temperature of the heat flow meter during the heat flow measurement and the CEX plasma during the probe measurement all influence the measurement results.The main reason for this phenomenon is that,from the perspective of particle impact,most of the energy is converted into heat when the particles hit the surface of the heat flux meter in a vertical direction in the case of forwarding incidence,while the oblique particles have a certain outward velocity,which reduces the energy exchange.Therefore,as the angle increases,the adaptation factor is reduced.

    3.2.Numerical simulation of heat flows

    In this section,the thermal effects of the electric propulsion plume are simulated for the LIPS-200 ion thruster and compared with experimental results.The calculation of the thermal effects takes into account the effect of the thermometer in the flow field on the plume,in particular the fact that energetic particles hitting the surface of the thermometer are neutralized and lose energy to become slow neutral atoms,but still participate in the xenon atomic flow field.The LIPS-200 ion thrust plume thermal effect simulation was carried out with dedicated software(EX-PWS)for electric propulsion plume simulation developed by Beihang University.The plasma dynamics model and the collision dynamics model in this software are solved by the particle-in-cell(PIC)and direct simulation Monte Carlo(DSMC)approaches,respectively.The EX-PWS has been successfully used to simulate the plume flow field characteristic distribution and plume sputtering deposition effects for single ion thruster and multiple ion thrusters.

    Figure 14.Schematic diagram of the calculation area mesh.

    The calculation grid and boundary conditions are set as follows.The schematic diagram of the calculation area mesh is illustrated in figure 14.The computational domain is a cube with dimensions of 1.5 m×1.5 m×1.5 m.The simulated LIPS-200 ion thruster is modeled as a small cylinder.The coordinates of the center point of the thruster exit surface are set to(0,0,0).To study the flow field distribution over a range of 1.0 m including the thruster,the entire calculation area was set to 1.0 m in front of the thruster and 0.5 m behind the thruster.The number of computational grids is approximately 4 million,divided into four parallel partitions along thex-direction.The virtual heat flow meters(in practice,the boundary does not exist and the particles move directly through it)were placed at 0.5 m,0.7 m,and 0.9 m downstream of the thruster.By statistically capturing the flux of particles at the location of the heat flux meter,the energy density of the incoming flow is obtained and,using the Maxwell reflection model,the energy flow of the reflected particles is statistically calculated,and the difference between the two is obtained as the heat flow density.In this simulation,the background pressure is set to 9×10?4Pa,while all cube boundary surfaces and virtual heat flow meter are set to free boundary and Maxwell wall,considering that the gas molecules are then absorbed after impacting on the surface and removed from the calculation domain.The ratio of Xe2+to Xe1+in the plume is assumed to be 10% for the calculation.For the calculations,the time step was set to 1×10?7s,the particle weights were set to 1×109,the ambient and individual wall temperatures were set to 300 K.

    The EX-PWS simulates particle injection,particle movement,and particle collision physics processes of the LIPS-200 thruster as detailed in the published literatures[4,13,21,22].This article does not duplicate excessive descriptions,it just presents the calculation method for heat flow based on the Maxwell model in the EX-PWS.For incoming plasma,the heat flow is divided into two parts,namely the energy flow of the plasma translational motion and the energy flow of the plasma rotation plus vibrational internal energy(qro+vi).

    The energy of plasma translational motion(qtran)can be expressed as:

    Combined with thermal accommodation coefficients obtained from experiments,direct statistics of particle fluxes at specific locations in the LIPS-200 ion thruster plume are employed to calculate the heat flow density at different locations.

    The Xe+ion number density and velocity distribution are depicted in figure 15.It is apparent from figure 15 that the ion number and velocity of plume Xe+are symmetrically distributed along the axial direction of the thruster,with a characteristic maximum number(approximately 6.5×1015m?3)and velocity(approximately 36 km s?1)of Xe+in the vicinity of the center of the thruster outlet and then decreasing along the radial and axial directions.

    It can be seen from the experimental results that the thermal accommodation coefficient is not constant for different regions and that it varies with angle and distance,which brings some challenges to the simulation.The experimentally measured thermal accommodation coefficients ranged from 0.72 to 0.99(close to 1).To further weigh the accuracy and speed of the simulation,three cases of the mean(0.88)and the minimum(0.72)and maximum(1)values of the thermal accommodation coefficient were adopted in this simulation for comparison with the experimental results.

    Figure 15.Xe+ ion number density(a)and velocity distribution(b).

    Figure 16.Simulation and experimental comparison of radial heat flow at 0.5 m(a),0.7 m(b),and 0.9 m(c),respectively,from the thruster outlet under different thermal accommodation factors.

    The simulation and experimental comparison of radial heat flow at 0.5 m(a),0.7 m(b),and 0.9 m(c),respectively,from the thruster outlet under different thermal accommodation factors,is depicted in figure 16.It is evident in figure 16(a)that the experimental measurements in the plume beam region at angles from 0° to 15° are essentially smaller than the simulated values for heat flow with a thermal accommodation factor of 1 and higher than the simulated values for heat flow with a thermal accommodation factor of 0.88 and 0.72,respectively.Moreover,the heat flow error between the experiment and the numerical simulation decreases with the decreasing of the angle.The error between the experiment and the numerical simulation increases with the increasing angle,with a maximum error at the angle of 15°.The relative absolute maximum error is within 12.1%for a thermal accommodation factor of 1 and within 18.5% and 26%for the thermal accommodation factors of 0.88 and 0.72,respectively.The location of the larger error in the heat flow with a thermal adaptation factor of 1 is approximately 15°.The reason for this is that this position is the half-angle of the beam,i.e.,the junction of the beam and non-beam regions,where the plasma number density is thinner and both the simulation and experiment are difficult,hence the high relative error.The thermal accommodation coefficients of 0.88 and 0.72 have large errors in the heat flow due to the experimental mean and minimum values to calculate the heat flow,and the actual measured heat flow at different locations at 0.5 m from the thruster has a certain error,so there is a difference between the simulation and experimental heat flow after comparison.It can also be seen that the thermal accommodation coefficients of 1,0.88,and 0.72 for the simulated heat flow,respectively,fall more gently in the radial direction compared to the experimental measurements,suggesting that the accommodation coefficients used in the Maxwell model for particle-surface energy exchange should not be fixed values,rather it should vary with angle.This study is consistent with the experimentally-measured variation of the thermal accommodation coefficient with angle.Analogously,similar conclusions can be drawn from figures 16(b)and(c).The difference is that the error in the simulation and experimental comparison increases with increasing distance and angle for the identical thermal accommodation factor.For different accommodation factors,the error in heat flow between experiment and simulation is larger with smaller thermal accommodation.However,it can also be seen that the error between the simulated and experimental heat flow with a thermal accommodation factor of 0.72 is generally higher within a small radial angle at 0.9 m.This indicates that the simulation results for the low thermal accommodation factor used in the simulation differ significantly from the actual results.The position with the highest error between the heat flow simulation and the experiment is at an angle of 15°,with a thermal adaptation factor of 0.72,an error of 28.6%.The main reasons for this are the distance from the thruster at 0.9 m,the thinner plasma,and the large error between simulation and experiment.It is found experimentally that the thermal adaptation factor varies for different locations in the plume.The simulation results based on the Maxwell thermal accommodation coefficient model show that the simulation and experimental errors of the heat flow with an average thermal accommodation coefficient and a minimum thermal accommodation coefficient are small over a wide range of region for a thermal accommodation coefficient of 1.Therefore,in order to improve the efficiency of the simulation,we suggest that the thermal adaptation coefficient of the simulated ion thruster plume based on the thermal adaptation coefficient of the Maxwell simulation is taken to be 1.

    To further weigh the accuracy and speed of the simulation,three values of the thermal accommodation coefficient were adopted in this simulation for comparison with the experimental results.In comparison with the experiment,the simulation results can predict the heat flow distribution of the thruster plume.In order to better predict the heat flow distribution at different positions of the thruster,the measured radial distribution values are subsequently input into the model for calculation,and the error between the simulation and calculation of the heat flow will be further reduced.

    4.Conclusions

    An experimental study of the thermal effects of the ion thruster plume was implemented with a total heat flow sensor and a radiant heat flow sensor in the region of 0.5–0.9 m downstream of the LIPS-200 ion thruster,with an angle of?60° to 60°.The TAC of the heat flow sensor surface in the plume is derived using a FP,a RPA,and a HFS.The PIC and DSMC methods based on the Maxwell thermal accommodation factor model(EX-PWS)were also used to simulate the thermal effects of the ion thruster plume.The conclusions obtained in this study are as follows:

    (1)At the identical radius and different angles,the total and radiant heat flow is maximum at the axis of the thruster and decreases successively away from the axis.The closer to the thruster,the greater the total and radiant heat flow in the plume for different radii.At the identical radius and different angles,the radiant heat flow accounts for a minimal part of the total heat flow,approximately 3%–5%.

    (2)The thermal accommodation factor,determined from the FP,a RPA,and heat flow sensor,is 0.72–0.99 over the measured region and is not a constant value depending on the angle.

    (3)The PIC and DSMC methods based on the Maxwell thermal accommodation coefficient model(EX-PWS)show a maximum error of 28.6% between simulation and experiment for LIPS-200 ion thruster plume heat flow which,on the one hand,provides an experimental basis for studying the interaction between the ion thruster and the spacecraft,and on the other hand provides optimization of the ion thruster plume simulation model.

    Acknowledgments

    The authors are grateful to National Natural Science Foundation of China(No.12005087)and the Science and Technology Program of Gansu Province(Nos.2006ZCTF0054,HTKJ2019KL510003,and 20JR10RA478).

    猜你喜歡
    趙勇
    時(shí)來運(yùn)轉(zhuǎn)
    《電訊技術(shù)》2020-2023年度“十佳通訊員”名單揭曉
    A study of the influence of different grid structures on plasma characteristics in the discharge chamber of an ion thruster
    800萬“諜中諜”:贅婿詐死亢奮了一對在后黃雀
    Quantitative analysis and time-resolved characterization of simulated tokamak exhaust gas by laser-induced breakdown spectroscopy
    被毒狗針戳殺的負(fù)心男
    趙勇:抓住脫貧突破口 實(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*
    Investigation of non-cohesive levee breach by overtopping flow*
    亚洲 欧美一区二区三区| 国产一区二区三区综合在线观看| 欧美日韩亚洲国产一区二区在线观看 | 大片免费播放器 马上看| 亚洲伊人久久精品综合| 国产福利在线免费观看视频| av视频免费观看在线观看| 女人被躁到高潮嗷嗷叫费观| 男女高潮啪啪啪动态图| 国产成人影院久久av| 日韩欧美三级三区| 成人特级黄色片久久久久久久 | 国产成人av激情在线播放| 国产男女内射视频| 国产麻豆69| 日本一区二区免费在线视频| 国产欧美日韩一区二区精品| 欧美精品一区二区免费开放| 高清黄色对白视频在线免费看| 男人舔女人的私密视频| 久久久精品区二区三区| 精品卡一卡二卡四卡免费| 亚洲精品乱久久久久久| 欧美日韩视频精品一区| 母亲3免费完整高清在线观看| 性色av乱码一区二区三区2| 在线天堂中文资源库| av视频免费观看在线观看| 国产高清激情床上av| 国产成人系列免费观看| 免费少妇av软件| 巨乳人妻的诱惑在线观看| 丝袜喷水一区| 欧美av亚洲av综合av国产av| 精品亚洲成国产av| 乱人伦中国视频| 悠悠久久av| 建设人人有责人人尽责人人享有的| 亚洲成人免费电影在线观看| bbb黄色大片| 国产91精品成人一区二区三区 | 亚洲精品久久成人aⅴ小说| 少妇猛男粗大的猛烈进出视频| 老鸭窝网址在线观看| 免费在线观看完整版高清| 亚洲中文字幕日韩| 色播在线永久视频| 成人三级做爰电影| 老司机午夜福利在线观看视频 | 最新美女视频免费是黄的| 正在播放国产对白刺激| 欧美精品av麻豆av| av天堂久久9| 免费在线观看视频国产中文字幕亚洲| 法律面前人人平等表现在哪些方面| 久久久精品国产亚洲av高清涩受| 亚洲av国产av综合av卡| 性高湖久久久久久久久免费观看| 好男人电影高清在线观看| 久久久久网色| 人人妻人人添人人爽欧美一区卜| 久久久久久亚洲精品国产蜜桃av| 国产黄色免费在线视频| 少妇裸体淫交视频免费看高清 | 侵犯人妻中文字幕一二三四区| 久久 成人 亚洲| 国产人伦9x9x在线观看| aaaaa片日本免费| 一级,二级,三级黄色视频| 久久中文字幕一级| av有码第一页| 天天添夜夜摸| 黄色视频,在线免费观看| 狠狠狠狠99中文字幕| 母亲3免费完整高清在线观看| 久久天躁狠狠躁夜夜2o2o| 亚洲人成电影观看| 精品国产一区二区三区久久久樱花| 一进一出好大好爽视频| 在线观看一区二区三区激情| 国产亚洲精品第一综合不卡| 波多野结衣av一区二区av| 丝瓜视频免费看黄片| 飞空精品影院首页| 亚洲国产欧美日韩在线播放| 亚洲专区中文字幕在线| 天天躁狠狠躁夜夜躁狠狠躁| 在线观看免费视频网站a站| 在线永久观看黄色视频| 少妇被粗大的猛进出69影院| 人妻一区二区av| 啦啦啦 在线观看视频| 国产精品一区二区在线观看99| 肉色欧美久久久久久久蜜桃| 国产福利在线免费观看视频| 中文字幕人妻丝袜制服| 亚洲精品在线观看二区| 亚洲专区字幕在线| 后天国语完整版免费观看| 一级,二级,三级黄色视频| 精品人妻在线不人妻| 少妇被粗大的猛进出69影院| 涩涩av久久男人的天堂| 日本五十路高清| 无遮挡黄片免费观看| 欧美精品高潮呻吟av久久| av线在线观看网站| 久久久久久人人人人人| av视频免费观看在线观看| 婷婷成人精品国产| 自线自在国产av| 欧美性长视频在线观看| 九色亚洲精品在线播放| 在线观看免费视频日本深夜| 亚洲国产欧美日韩在线播放| 亚洲伊人色综图| 巨乳人妻的诱惑在线观看| av线在线观看网站| 日韩制服丝袜自拍偷拍| 国产欧美日韩一区二区三区在线| 又黄又粗又硬又大视频| 亚洲中文日韩欧美视频| 久久99热这里只频精品6学生| 老司机福利观看| 操出白浆在线播放| 十八禁人妻一区二区| 99riav亚洲国产免费| 日韩人妻精品一区2区三区| 满18在线观看网站| 亚洲av第一区精品v没综合| 国产伦人伦偷精品视频| 中文字幕人妻丝袜制服| videosex国产| 久久精品国产亚洲av香蕉五月 | 亚洲全国av大片| 欧美亚洲日本最大视频资源| 国产精品亚洲一级av第二区| 日韩中文字幕欧美一区二区| 亚洲全国av大片| 国产淫语在线视频| 亚洲精品美女久久久久99蜜臀| av有码第一页| 久久天躁狠狠躁夜夜2o2o| 极品人妻少妇av视频| 老司机福利观看| 天堂8中文在线网| 亚洲专区字幕在线| 欧美精品亚洲一区二区| 又黄又粗又硬又大视频| 精品久久久精品久久久| 在线观看舔阴道视频| 亚洲自偷自拍图片 自拍| a级毛片黄视频| 丰满饥渴人妻一区二区三| 国产深夜福利视频在线观看| 色在线成人网| 日韩 欧美 亚洲 中文字幕| 亚洲成人国产一区在线观看| 亚洲人成77777在线视频| 国产伦人伦偷精品视频| av网站免费在线观看视频| 叶爱在线成人免费视频播放| 欧美午夜高清在线| 国产日韩欧美亚洲二区| 亚洲国产欧美日韩在线播放| 亚洲欧美精品综合一区二区三区| 免费在线观看日本一区| 一级片'在线观看视频| 久9热在线精品视频| 亚洲欧美一区二区三区黑人| 一本一本久久a久久精品综合妖精| 国产日韩一区二区三区精品不卡| 捣出白浆h1v1| 欧美精品啪啪一区二区三区| 69av精品久久久久久 | 一个人免费看片子| 熟女少妇亚洲综合色aaa.| 他把我摸到了高潮在线观看 | 男女床上黄色一级片免费看| 欧美精品高潮呻吟av久久| 国产精品香港三级国产av潘金莲| 黄片小视频在线播放| 国产伦人伦偷精品视频| 久久九九热精品免费| 精品久久久久久电影网| 亚洲欧洲日产国产| 两个人看的免费小视频| 亚洲国产欧美网| 欧美日韩av久久| 精品一区二区三区视频在线观看免费 | 亚洲性夜色夜夜综合| 色在线成人网| 久热这里只有精品99| 日韩熟女老妇一区二区性免费视频| av超薄肉色丝袜交足视频| 亚洲全国av大片| 久久天堂一区二区三区四区| 日韩欧美免费精品| 男女免费视频国产| 亚洲第一青青草原| 国产精品av久久久久免费| 国产午夜精品久久久久久| 久久国产精品影院| 黑人巨大精品欧美一区二区mp4| a级毛片黄视频| 在线 av 中文字幕| 老司机在亚洲福利影院| 亚洲欧美精品综合一区二区三区| 中文亚洲av片在线观看爽 | 人妻久久中文字幕网| 黑人猛操日本美女一级片| √禁漫天堂资源中文www| 亚洲精品久久午夜乱码| 黑人操中国人逼视频| 老熟妇乱子伦视频在线观看| 变态另类成人亚洲欧美熟女 | 亚洲欧美一区二区三区久久| 九色亚洲精品在线播放| 国产亚洲精品第一综合不卡| 人妻一区二区av| 精品国产乱码久久久久久男人| 两性夫妻黄色片| 免费av中文字幕在线| 天天躁日日躁夜夜躁夜夜| 精品国产超薄肉色丝袜足j| 精品亚洲乱码少妇综合久久| 一区福利在线观看| 久久久欧美国产精品| 久久久水蜜桃国产精品网| 中文字幕另类日韩欧美亚洲嫩草| 亚洲精品自拍成人| 国产精品欧美亚洲77777| 成人永久免费在线观看视频 | 亚洲av日韩精品久久久久久密| 三上悠亚av全集在线观看| 亚洲精品国产区一区二| 国产高清激情床上av| 国产精品av久久久久免费| 精品国产乱码久久久久久小说| 色尼玛亚洲综合影院| 免费在线观看视频国产中文字幕亚洲| 丰满饥渴人妻一区二区三| 啦啦啦免费观看视频1| 国产一区二区三区综合在线观看| 午夜福利在线观看吧| 国产成人欧美在线观看 | 午夜激情久久久久久久| 久久国产亚洲av麻豆专区| 欧美久久黑人一区二区| 免费在线观看日本一区| 777米奇影视久久| 亚洲精品av麻豆狂野| 国产精品一区二区在线不卡| 国产精品欧美亚洲77777| 欧美另类亚洲清纯唯美| 亚洲精品乱久久久久久| 精品国产乱子伦一区二区三区| 黑人猛操日本美女一级片| 一级毛片精品| 国产免费视频播放在线视频| 亚洲三区欧美一区| 国产xxxxx性猛交| 午夜福利在线免费观看网站| 欧美乱妇无乱码| 国产精品欧美亚洲77777| 又黄又粗又硬又大视频| 精品人妻1区二区| 亚洲精品美女久久久久99蜜臀| 久久免费观看电影| 大码成人一级视频| netflix在线观看网站| 成人国产一区最新在线观看| 男女床上黄色一级片免费看| 亚洲中文av在线| 亚洲av电影在线进入| 亚洲av第一区精品v没综合| 成人手机av| 免费久久久久久久精品成人欧美视频| 日日摸夜夜添夜夜添小说| av不卡在线播放| 热re99久久精品国产66热6| 国产在线一区二区三区精| 一区福利在线观看| 国产在线免费精品| 欧美在线一区亚洲| 国产精品秋霞免费鲁丝片| 丰满人妻熟妇乱又伦精品不卡| 9191精品国产免费久久| 激情视频va一区二区三区| 亚洲国产看品久久| 女人久久www免费人成看片| 男女下面插进去视频免费观看| 日本一区二区免费在线视频| 欧美国产精品一级二级三级| 窝窝影院91人妻| 高清在线国产一区| 国产成人一区二区三区免费视频网站| 菩萨蛮人人尽说江南好唐韦庄| 91国产中文字幕| 免费在线观看黄色视频的| 视频区图区小说| 日本撒尿小便嘘嘘汇集6| 亚洲成国产人片在线观看| 国产成人欧美在线观看 | 亚洲成人免费av在线播放| 十八禁高潮呻吟视频| 久久久国产欧美日韩av| 自线自在国产av| 女人被躁到高潮嗷嗷叫费观| 中文亚洲av片在线观看爽 | 18禁黄网站禁片午夜丰满| e午夜精品久久久久久久| 老司机福利观看| 伦理电影免费视频| 757午夜福利合集在线观看| 十八禁高潮呻吟视频| 大陆偷拍与自拍| 久久精品国产99精品国产亚洲性色 | 日韩 欧美 亚洲 中文字幕| 男女免费视频国产| 免费久久久久久久精品成人欧美视频| 九色亚洲精品在线播放| 日韩精品免费视频一区二区三区| 亚洲一码二码三码区别大吗| 欧美日韩亚洲国产一区二区在线观看 | 12—13女人毛片做爰片一| 黑人巨大精品欧美一区二区蜜桃| 国产精品麻豆人妻色哟哟久久| 在线永久观看黄色视频| 中文字幕人妻丝袜一区二区| 69av精品久久久久久 | 下体分泌物呈黄色| 亚洲中文字幕日韩| 啦啦啦免费观看视频1| 免费在线观看完整版高清| 欧美黄色片欧美黄色片| 国产精品二区激情视频| 99香蕉大伊视频| 王馨瑶露胸无遮挡在线观看| 无遮挡黄片免费观看| 夜夜爽天天搞| 啦啦啦视频在线资源免费观看| 女人高潮潮喷娇喘18禁视频| 久久精品国产亚洲av高清一级| 国产片内射在线| 免费在线观看影片大全网站| 欧美日本中文国产一区发布| 国产精品美女特级片免费视频播放器 | 热99国产精品久久久久久7| 久热爱精品视频在线9| 少妇 在线观看| 久久久精品区二区三区| 欧美亚洲日本最大视频资源| 天天躁夜夜躁狠狠躁躁| 极品少妇高潮喷水抽搐| 国产xxxxx性猛交| 午夜福利一区二区在线看| 考比视频在线观看| 久久久精品免费免费高清| 18禁国产床啪视频网站| 欧美激情 高清一区二区三区| 亚洲视频免费观看视频| 国产伦理片在线播放av一区| 十八禁网站网址无遮挡| 视频区欧美日本亚洲| 久久人妻熟女aⅴ| 麻豆av在线久日| 波多野结衣一区麻豆| 国产成人一区二区三区免费视频网站| 亚洲成国产人片在线观看| 91麻豆av在线| 黑人猛操日本美女一级片| 久久久久精品国产欧美久久久| 国产成人av激情在线播放| 亚洲欧美精品综合一区二区三区| h视频一区二区三区| 男女免费视频国产| 老司机午夜福利在线观看视频 | 麻豆乱淫一区二区| 久久久久精品人妻al黑| 精品一区二区三区四区五区乱码| 中文字幕精品免费在线观看视频| avwww免费| 如日韩欧美国产精品一区二区三区| 丁香六月欧美| 亚洲国产毛片av蜜桃av| 视频在线观看一区二区三区| 精品国产乱子伦一区二区三区| 日韩欧美国产一区二区入口| 黑丝袜美女国产一区| 最黄视频免费看| 热re99久久国产66热| 欧美 日韩 精品 国产| 国产一区二区激情短视频| 久久这里只有精品19| 免费日韩欧美在线观看| 久久久精品免费免费高清| 国产av又大| 中文字幕高清在线视频| 欧美+亚洲+日韩+国产| 黄色视频,在线免费观看| 狂野欧美激情性xxxx| av免费在线观看网站| 欧美黄色片欧美黄色片| 午夜老司机福利片| a级毛片黄视频| 亚洲精品美女久久久久99蜜臀| 成年女人毛片免费观看观看9 | 欧美亚洲 丝袜 人妻 在线| 真人做人爱边吃奶动态| 丝袜喷水一区| 男女下面插进去视频免费观看| 国产精品 欧美亚洲| 热re99久久精品国产66热6| 午夜福利在线免费观看网站| 大码成人一级视频| 国产精品99久久99久久久不卡| 亚洲欧美色中文字幕在线| 天堂动漫精品| 热99re8久久精品国产| 中文字幕最新亚洲高清| 国产亚洲一区二区精品| 男女无遮挡免费网站观看| cao死你这个sao货| 久久久久久亚洲精品国产蜜桃av| 久久国产亚洲av麻豆专区| 日韩中文字幕欧美一区二区| 91麻豆精品激情在线观看国产 | 国产国语露脸激情在线看| 亚洲欧美日韩高清在线视频 | 狠狠婷婷综合久久久久久88av| 首页视频小说图片口味搜索| 国产精品久久久久成人av| av线在线观看网站| 香蕉丝袜av| 女人精品久久久久毛片| 老司机深夜福利视频在线观看| 国产在线观看jvid| 久久久久久久精品吃奶| 国产单亲对白刺激| 亚洲va日本ⅴa欧美va伊人久久| 国产一卡二卡三卡精品| 老司机午夜十八禁免费视频| 欧美人与性动交α欧美软件| 99久久人妻综合| 亚洲av第一区精品v没综合| 成人18禁高潮啪啪吃奶动态图| 999久久久国产精品视频| 99国产精品一区二区蜜桃av | 日韩人妻精品一区2区三区| 看免费av毛片| 在线观看一区二区三区激情| 十八禁人妻一区二区| 黄色毛片三级朝国网站| 欧美av亚洲av综合av国产av| 久久天躁狠狠躁夜夜2o2o| 91精品国产国语对白视频| 99久久人妻综合| 肉色欧美久久久久久久蜜桃| 女人精品久久久久毛片| 老司机亚洲免费影院| 国产又色又爽无遮挡免费看| 欧美乱妇无乱码| 两性午夜刺激爽爽歪歪视频在线观看 | 中文字幕高清在线视频| 免费在线观看完整版高清| 可以免费在线观看a视频的电影网站| 日本黄色视频三级网站网址 | 人人妻人人澡人人爽人人夜夜| 国产激情久久老熟女| 亚洲五月色婷婷综合| 国内毛片毛片毛片毛片毛片| 搡老乐熟女国产| av欧美777| 视频区图区小说| 黑人巨大精品欧美一区二区mp4| 99精国产麻豆久久婷婷| 中文字幕最新亚洲高清| 免费黄频网站在线观看国产| 久久午夜亚洲精品久久| 欧美激情久久久久久爽电影 | 国产亚洲欧美在线一区二区| 亚洲国产欧美一区二区综合| 中文字幕最新亚洲高清| 中文字幕制服av| 丝袜美腿诱惑在线| 久久99热这里只频精品6学生| 成人国产一区最新在线观看| 国产精品av久久久久免费| 免费久久久久久久精品成人欧美视频| 久久国产精品影院| 黄色视频,在线免费观看| 一级片'在线观看视频| 国产亚洲av高清不卡| 国产男女超爽视频在线观看| 大码成人一级视频| 女同久久另类99精品国产91| 久久人妻福利社区极品人妻图片| 最新美女视频免费是黄的| 亚洲第一av免费看| 亚洲人成电影免费在线| 亚洲精品一卡2卡三卡4卡5卡| 亚洲av日韩精品久久久久久密| 欧美久久黑人一区二区| av线在线观看网站| 成人黄色视频免费在线看| 国产精品99久久99久久久不卡| 俄罗斯特黄特色一大片| 黄色a级毛片大全视频| 国产91精品成人一区二区三区 | 亚洲成人免费电影在线观看| 亚洲av欧美aⅴ国产| 热re99久久国产66热| 免费在线观看视频国产中文字幕亚洲| 午夜成年电影在线免费观看| 国产亚洲精品第一综合不卡| 亚洲国产av新网站| 热re99久久精品国产66热6| 日韩大码丰满熟妇| 午夜两性在线视频| 日本av免费视频播放| 欧美+亚洲+日韩+国产| 国产精品秋霞免费鲁丝片| 自线自在国产av| 欧美成狂野欧美在线观看| 建设人人有责人人尽责人人享有的| 成年人黄色毛片网站| 老熟妇乱子伦视频在线观看| 国产一区二区三区在线臀色熟女 | 老司机靠b影院| 国产精品 国内视频| 成人永久免费在线观看视频 | 久久ye,这里只有精品| 亚洲国产毛片av蜜桃av| 欧美日韩一级在线毛片| 国产精品亚洲一级av第二区| 麻豆av在线久日| 天天躁日日躁夜夜躁夜夜| 欧美日本中文国产一区发布| 中文字幕最新亚洲高清| 久久亚洲精品不卡| 最新在线观看一区二区三区| 成人av一区二区三区在线看| 中国美女看黄片| 免费日韩欧美在线观看| 国产成人精品在线电影| 一级毛片女人18水好多| 亚洲精品美女久久久久99蜜臀| 久久午夜综合久久蜜桃| 午夜福利乱码中文字幕| 最新在线观看一区二区三区| 老熟妇仑乱视频hdxx| 香蕉国产在线看| 亚洲精品成人av观看孕妇| 国产免费福利视频在线观看| 中文字幕高清在线视频| 精品少妇黑人巨大在线播放| 成人影院久久| 国产高清激情床上av| 婷婷丁香在线五月| 精品亚洲成国产av| 国产成人一区二区三区免费视频网站| 国产高清视频在线播放一区| 视频在线观看一区二区三区| 欧美精品亚洲一区二区| 亚洲av国产av综合av卡| 久久中文字幕人妻熟女| 亚洲av成人不卡在线观看播放网| 别揉我奶头~嗯~啊~动态视频| 午夜两性在线视频| 国产精品久久久久久人妻精品电影 | 中文字幕精品免费在线观看视频| 嫩草影视91久久| 91国产中文字幕| 午夜福利影视在线免费观看| 女人被躁到高潮嗷嗷叫费观| 国产精品香港三级国产av潘金莲| 高清黄色对白视频在线免费看| 日日夜夜操网爽| 久久久欧美国产精品| 一区二区av电影网| 亚洲色图av天堂| 99精国产麻豆久久婷婷| 少妇猛男粗大的猛烈进出视频| 18禁国产床啪视频网站| 麻豆国产av国片精品| 国产不卡一卡二| 国产av又大| 国产精品98久久久久久宅男小说| 免费观看av网站的网址| 免费一级毛片在线播放高清视频 | 日韩视频一区二区在线观看| 黑人猛操日本美女一级片| 日韩免费av在线播放| 欧美精品人与动牲交sv欧美| 我的亚洲天堂| 99re在线观看精品视频| 久9热在线精品视频| 国产精品久久久久久精品电影小说| 亚洲久久久国产精品| 精品久久蜜臀av无| 蜜桃国产av成人99| 久久精品国产亚洲av高清一级| 黑人欧美特级aaaaaa片| 在线观看人妻少妇| 视频区图区小说| 91字幕亚洲| 18禁黄网站禁片午夜丰满| 999久久久国产精品视频| 人人妻人人澡人人看| 国产在线一区二区三区精| 一区二区av电影网| 涩涩av久久男人的天堂| 首页视频小说图片口味搜索| 国产一区有黄有色的免费视频|