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

    Numerical simulation and experimental validation of multiphysics field coupling mechanisms for a high power ICP wind tunnel?

    2021-06-26 03:29:44MingHaoYu喻明浩ZheWang王哲ZeYangQiu邱澤洋BoLv呂博andBoRuiZheng鄭博睿
    Chinese Physics B 2021年6期

    Ming-Hao Yu(喻明浩) Zhe Wang(王哲) Ze-Yang Qiu(邱澤洋)Bo Lv(呂博) and Bo-Rui Zheng(鄭博睿)

    1Faculty of Mechanical and Precision Instrument Engineering,Xi’an University of Technology,Xi’an 710048,China 2School of Automation and Information Engineering,Xi’an University of Technology,Xi’an 710048,China

    Keywords: inductively coupled plasma,multiphysics field,coupling mechanism,simulation and experiment

    1. Introduction

    Reconnaissance satellites,manned spacecrafts,and space shuttles should return to the ground in time after they have completed their orbital missions. However, because of their high speed flight, they encounter severe aerodynamic heating problems when they re-enter the dense atmosphere. Accordingly,a perfect thermal protection system is required to ensure the safety of the spacecraft and astronauts. To improve the thermal protection system of re-entry vehicles,an inductively coupled plasma (ICP) wind tunnel is built on the ground for investigating the thermal protection materials for these vehicles.

    The ICP wind tunnel can provide flawless long-term stable operation of a high enthalpy plasma jet and accurately simulate the flight environment of re-entry aircraft.It is,therefore,a piece of ideal equipment for investigating the plasma sheath of hypersonic aircraft. The ICP wind tunnel is mainly composed of an intake part, a high frequency discharge part, and a test part. The structural diagram of an ICP wind tunnel is shown in Fig.1. The high frequency discharge part is mainly comprised of the high frequency power supply,induction coil,and quartz tube.The electrified induction coil is wound around the high temperature resistant quartz tube. After the working gas enters the quartz tube from the air inlet,it is rapidly heated in the induction coil area to form a high temperature high enthalpy induction coupling plasma. The plasma produced in the quartz tube (also called the plasma torch) usually continues to flow forward into the large test chamber. In this chamber,the ICP can be used to test the aerodynamic performance of the aircraft,develop thermal protection materials[1]and investigate communication blackouts and electromagnetic wave transmissions.[2]

    The ICP has been widely utilized in aerospace, materials, chemical industry, and other fields since the 1960s and 1970s. Presently,with rapid industrial development,many researchers have performed experiments and numerical simulations on ICPs.[1–5]In an experimental research,Moonet al.[4]measured the plasma density and coil voltage in ICP by using floating and high voltage probes,respectively,and studied the influence of coil turns on the ICP.They found that the plasma density increased with the number of coil turns because of the enhancement of inductive coupling. Moreover, the coil voltage was affected by both the plasma density and the coil.Their research can considerably aid the design of future ICP sources.Kimet al.[5]performed experiments on a CO2ICP that had a pressure of 100 mTorr (1 Torr=1.33322×102Pa) and power of 13.56 MHz. Electron energy distribution functions were measured using a Langmuir probe, and other plasma parameters, such as electron number density and electron temperature,were also obtained. According to the electrostatic probe technology, Panet al.[6]developed an electrostatic probe diagnosis system suitable for the test environment of a high frequency plasma wind tunnel. This system was employed to diagnose the flow field parameters of the ICP wind tunnel. Fujitaet al.[7]measured the temperature parameters of plasma in an ICP wind tunnel by means of a radiation spectrum.

    Fig.1. Structural diagram of inductively coupled plasma wind tunnel.

    Regarding the numerical studies of the ICPs,Ivanov and Zverev[8]numerically studied the influence of hydrogen on the parameters of plasma in the ICP torch with considering chemical reactions of Ar with H2. They found that as the hydrogen percentage increases over 10% at a constant total Ar+H2mass flow rate, the electrical conductivity of plasma increases slightly, but the enthalpy of plasma increases significantly due to high value of the dissociation energy of hydrogen molecules. Degrezet al.[9]also numerically studied the air ICP under a local thermodynamic equilibrium and a chemical non-equilibrium. They provided the calculation results of the streamline, temperature, and chemical concentration field in the air ICP torch. Alavi and Mostaghimi[10]numerically and experimentally studied a novel ICP torch which possesses a conical nozzle. The inductive coil were wound around the conical nozzle rather than around the conventional quartz tube. They found that the new torch has higher power density, better plasma stability, and better resistance against extinguishing factors than the traditional torch. Namet al.[11]conducted a numerical analysis of the power distribution and heat flow characteristics of Ar–N2ICP with nitrogen content of 0 mol%–50 mol%at a plasma power level of 50 kW.They obtained relevant data, such as heat loss and load resistance.They also employed numerical simulation and equivalent circuit analysis to optimize the design and analysis of a high power ICP torch system with a vacuum tube oscillator. Lu and Feng[12]investigated the effects of the operating parameters on the ICP properties by a non-equilibrium magnetohydrodynamic model. Their numerical model can aid in determining the optimal conditions for producing the ideal plasma.Sunet al.[13,14]studied ICP plasma characteristics by using two-dimensional self-consistent fluid model. It was found that the spatial distribution of plasma characteristics was controlled by changing the bias power of single frequency bias or high frequency bias power ratio of double frequency offset. The plasma density and ion species ratio were adjusted by changing bias voltage and pulse frequency. Gaoet al.[15]found that different RF bias voltage has different effect on plasma characteristics: in the low RF voltage range, the electron number density decreases and the effective electron temperature increases; in the high RF voltage range the electron number density increases and the effective electronic temperature decreases.

    In addition, with the rapid development of commercial softwares such as ANSYS, COMSOL,etc., many research groups prefer to simulate the ICP using these as simulation tools. For example, Zhuet al.[16]used the ANSYS FLUENT software to simulate a two-dimensional ICP torch and derived the electromagnetic field distribution,temperature distribution,and velocity distribution in the torch. Chenet al.[17]employed the ANSYS FLUENT software to calculate the temperature of pure argon thermal plasma and the flow field distribution in the plasma generator device. They found that the temperature distribution and flow field distribution can be controlled by changing the corresponding device structure and working frequency. Zhao[18]established a fluid model of argon ICP by using the COMSOL software.It was observed that the temperature of the argon ICP changed non-monotonically with power. To explain this phenomenon, a new equation for the average electron energy was proposed. Leiet al.[19]employed the COMSOL software to simulate the ICP discharge and studied the influences of different coil designs and gas compositions on the discharge characteristics of a twodimensional ICP. The best discharge characteristic is considered to be afforded by a 20-mm coil gap, with the first coil positioned 130 mm away from the top of the device, and the best discharge mode is considered to be argon discharge mode.This numerical model is a typical ICP discharge,and its results are important for further experimental studies.

    Although the above-mentioned research groups have conducted in-depth research on the numerical simulation of ICP,the majority of their simulation results have not been validated experimentally. Moreover,most of researches focused on low power(1 kW–50 kW)ICP wind tunnels. There is a lack of detailed studies to reveal the flow characteristics and the multiphysics coupling mechanism of the high power(such as a 1.2-MW class) ICP wind tunnel. For example, on the one hand,for a low power ICP wind tunnel,[20–23]its coil current,working pressure, plasma density are usually small;[20–22]hence,its electromagnetic field is small,[22,23]and the coupling effect between the electromagnetic field and the flow field is also relatively weak. One the other hand, for a high power ICP wind tunnel,[24,25]its working parameters such as the input power,operating pressure and the discharge frequency are quite different from those for a low power ICP wind tunnel.Its plasma density, electric-field intensity, plasma temperature, and heat flux are all higher than those of the low power one.[26,27]Thus,for the high power ICP wind tunnel,the coupling effects among the electromagnetic, flow, chemical, and thermodynamic fields are much stronger and more complex.Therefore,further researches on the high power ICP wind tunnel should be conducted to make clear such complex coupling effects. While so far there has been few such work available.

    In the present study,a 1.2-MW high power ICP wind tunnel is studied as a research object. Furthermore,the coupling mechanisms of the electromagnetic field,flow field,chemical field,and thermal field of the wind tunnel are explored and discussed by combining numerical simulation with experimental validation. The electron number density, temperature, pressure,electric field intensity,Lorentz force,and other important parameters are discussed and investigated through numerical simulation. The plasma parameters under different working conditions are compared with each other to reveal the distribution characteristics of the flow field, electromagnetic field,chemical field, and thermal field of the high power wind tunnel. The accuracy of the numerical simulation is verified by comparing the simulation results with experimental data.[24]

    In the following sections, the computational grid, operating conditions, control equations, numerical solution methods,and boundary conditions of the high power ICP wind tunnel are presented. Thereafter,by analyzing the distribution of flow field parameters(e.g.,temperature,electron density,and pressure) and electromagnetic field parameters (e.g., electric field intensity and Lorentz force) of the ICP wind tunnel under two different power values, the interaction and influence mechanism between the flow field and electromagnetic field in the high power ICP wind tunnel are determined. Finally,the accuracy of the research results is verified by comparing the experimental results with the numerical results.

    2. Research methods

    The research object of this study is an ICP wind tunnel with a power of up to 1.2 MW. The plasma torch is made of a quartz tube with an inner diameter and a wall thickness of

    80 mm and 3 mm,respectively;six electrified coils are wound around this tube. As shown in Fig.2, the calculation domain of the electromagnetic field is?100 mm≤x ≤570 mm and 0 mm≤y ≤450 mm,with each interval consisting of 161×97 grid nodes;the ICP torch consists of 111×41 grid nodes. The whole flow field grids including the vacuum chamber are composed of 329×41 nodes as shown in Fig.3.

    To calculate the electromagnetic field near the coil and air velocity at the air inlet accurately,the calculation grid is encrypted in both the area of the induction coil and the area of inlet.In this study,it is assumed that the plasma is steady,electrically quasineutral and two-dimensionally axisymmetric. The fluid micro-element is in a thermochemical non-equilibrium state. In the simulation calculation, the Lorentz force and induced magnetic field generated by the plasma are considered.

    Fig.2. Computational mesh of electromagnetic-field and ICP torch.

    Fig.3. Computational mesh of the whole flow field.

    Usually, the plasma discharge is ignited by an appropriate argon injection and thereafter,gradually replaced by an air mixture (16 g/s standard dust-free atmospheric air) until stable working conditions are reached. To measure the temperature of plasma particles,spectrum analyzers are located at 195,320,and 445 mm away from the plasma torch outlet. For the simulation study,two typical working conditions: cases 1 and 2 for the 1.2-MW ICP wind tunnel are selected. The specific working conditions are listed in Table 1.The influences of different working conditions on the plasma temperature,velocity,pressure,electron density,electric field intensity,and other parameters are compared to demonstrate the flow field and electromagnetic field distribution characteristics of the high power ICP wind tunnel. These comparisons are based on the multifield coupling solutions of the flow field,electromagnetic field,thermodynamic field,and chemical field.

    Table 1. Working conditions of computational cases.

    2.1. Governing equations

    The vector form of the system of non-equilibrium magnetohydrodynamic equations is as follows:

    The gas temperature consists of the translational temperature(Ttr),rotational temperature(Trot),vibration temperature(Tvib),and electronic temperature(Te). The conservation vectorQ, flux vectorsFandFv, and source term vectorWin Eq.(1)are expressed as follows:

    whereδi,j,FLi,Sjoule,Sint, ?h0s, andΘvib,sdenote the Kronecker function,Lorentz force,Joule heating efficiency,internal energy exchange rate,component formation enthalpy,and vibration characteristic temperature, respectively; subscriptsrepresents the chemical species;Mis the molecular number;˙ωsis the mass generation rate of substancesand is the source term of the mass conservation equation of the particle at line 4 inW,Xs,andDsare the mole fraction and diffusion coefficient ofs,respectively;ui(oruj)andxjrepresent the velocity and tensor in the rectangular coordinate system, respectively;τijandqjare the stress tensor and the heat flux vector,respectively,and expressed as

    The total internal energy(E)of the gas is defined as follows:

    The translational energy(Etr),rotational energy(Erot),vibration energy(Evib),and electronic internal energy(Ee)of gas particle are calculated by the corresponding equations given in the work of Yuet al.[22]

    The gas viscosityμ, and translationalλtr, vibrationalλvib, rotational thermal conductivityλrotin the above equations are all computed from Yos’ formulas.[28,29]The diffusion coefficientDsis calculated from the formula of Curtiss and Hirschfelder.[30]The conductivity of plasmaσand the electron thermal conductivityλeare evaluated from the thirdorder Sonine polynomials of air.[31–33]

    The Joule heating rateSjoule, and the axial and radial Lorenz force,i.e.,FLxandFLycan be calculated from the following equations,respectively:

    Here,Eθ=ER+ iEI.ERandEIare the real and imaginary part of the electric field intensity, respectively.They are obtained by solving the magnetic vector-potential equations[22]at each time step.

    In addition,in Eq.(2),the rotational energy exchange rate(Sint,rot),vibration energy exchange rate(Sint,vib),and electron energy exchange rate(Sint,e)are calculated from the following equations:

    Here,νb,rsandνf,rsare the stoichiometric coefficient of componentsbefore and after therchemical reaction, respectively;Msrepresents the molar mass of substances;sandjare the subscripts representing the component numbers of chemical species. The mass production rate of speciess,i.e.,˙ωsacts as a mass source term in the mass conservation equation of chemical species in the fourth row of the source termWin Eq. (2). It participates in iterative computations of the non-equilibrium magnetohydrodynamic equations mentioned above,and updates solutions at each time step.

    2.2. Iterative computational methods

    To iteratively solve the governing equations mentioned above, the flow field equations are discretized by the finitevolume method. The viscous term is calculated by the secondorder center difference method. To compute the inviscid fluxes of the flow field equations at the cell interfaces, the simple low dissipation advection upwind[36](SLAU)splitting method, which is a parameter-free and simple compressible numerical flux function of the AUSM family method,is used.The monotonic upstream-centered scheme for conservation law(MUSCL).[37,38]is adopted to avoid computational oscillations and spurious solutions in the region of high gradients.In the presently-used MUSCL, the physical properties at the cell interface are determined by linearly interpolating properties between neighboring cells. Thus, the accuracy of the inviscid numerical flux is kept to the second order. To overcome the numerical stiffness for a chemical nonequilibrium plasma,the point-implicit method[39]is utilized in this study. In addition, the point-implicit lower and upper symmetric Gauss–Seidel(LU-SGS)method[40]is used to update the solutions of the flow field results at each time step.

    The finite-difference method is adopted to discretize the magnetic vector-potential equations to figure out the distributions of the electric field intensitiesERandEI. The secondorder, central-difference method is used to evaluate their numerical flux. The sub-relaxation iterative method is used to solve the discretized electromagnetic-field equations stably.The relaxation factor is set to be 0.5. When the relative error of the magnetic-vector potential betweennandn+1 steps is less than 10?5,the calculation is considered to be converged.

    In this study,the local time stepping scheme is adopted to update the flow field results. When the residual of the mass,the momentum and the energy conservative equation are less than 10?5, 10?5, and 10?3, respectively, then the simulation is thought to be converged. In the present study, because the chemical nonequilibrium process of the chemical reaction of air is considered and modeled in the simulation,the residual of the energy conservative equation can hardly reach 10?5.Thus,when the residual of the energy conservative equation reaches 10?3and the flow parameters such as the temperature,velocity and pressure are no more changed after 1.2 million iterative computational steps,then we think the simulation is converged. In addition, all the equations and models mentioned above are programmed by the computer language Fortran 90 by ourselves. The computational mesh of the ICP wind tunnel is divided by Gridgen V15. All simulations are carried out by using the Tianhe-2 Supercomputer.

    2.3. Boundary conditions

    The working gas is injected into the quartz tube; the mass flow rate and initial gas temperature arem=16 g/s andT0=300 K, respectively. At the exit of the ICP wind tunnel laboratory,the working pressure is set to be a fixed value,and other flow properties are evaluated through the linear interpolation of adjacent internal grid points. For the central axis, the axisymmetric boundary condition (Qi,j=Qi,j+1) is adopted, whereQdenotes the conservation variable. For the magnetic vector-potential equations,the outer boundary of the electromagnetic field calculation area is set to bex=?10,x=570, andy=450 mm. The electric field intensitiesERandEIare set to be zero at these outer boundaries:Ei,0=Ei,1.It is assumed that there is no slippage nor catalytic effect nor pressure gradient on the wall of the plasma torch. The wall temperature is determined using the heat conduction equation,?Tw/?n=qjmax, whereqjmaxis the heat flux density at the inner wall. For the wall in the test room,a constant wall temperature boundary condition is adopted,i.e.,Tw=300 K.

    3. Results and discussion

    3.1. Numerical results

    In this subsection, we will show the simulated results of flow,electromagnetic,and chemical fields,respectively. As is well known,when the input power is different,the flow properties such as the pressure, temperature, velocity also change and may demonstrate different physical results. So in order to make clear this issue for the high power ICPWT, we take case 1 and case 2 listed in Table 1 as baseline cases to show the flow field differences under the different input power conditions. The simulated results of these two cases will be discussed comparatively below.

    Fig. 4. Distribution of axial Lorentz force (upper half) and radial Lorentz force(lower half)for case 1.

    Figures 4 and 5 show the distributions of the axial and radial Lorentz forces in the plasma torches for case 1 and case 2,respectively. First,it can be observed that the maximum axial and radial Lorentz force position are under the third coil and fourth coil,respectively.The axial Lorentz forces in both cases are positive and negative, whereas the radial Lorentz forces are both always negative. The radial Lorentz force direction is from the wall to the central axis. Second, the maximum values of the radial Lorentz force are 2.23 and 3.05 times greater than the maximum values of the axial Lorentz force for case 1 and case 2, respectively. This proves that the radial momentum transfer of the plasma is stronger than the axial momentum transfer. Finally,because the greater the input power,the greater the electric field intensity and Lorentz force are.Hence, the Lorentz force for case 2 is greater than that for case 1.

    Fig. 5. Distribution of axial Lorentz force (upper half) and radial Lorentz force(lower half)for case 2.

    Figures 6 and 7 show the distribution diagrams of the electronic number density and imaginary part of the electricfield intensity in the plasma torch for case 1 and case 2,respectively. The electronic number density and electric field intensity considerably vary in the plasma torch,and their maximum value positions are under the third and fourth coil,respectively.The upper half of the figure shows that the electron number density is maximum at the central axis and minimum at the wall. This is because the radial Lorentz force is negative,and its direction points from the wall to the central axis. The electrons are gathered at the central axis under the action of the radial Lorentz force. The maximum electron number densities for cases 1 and 2 are 1.01×1022m?3and 9.05×1021m?3,respectively. The lower half for each of the figures shows that the maximum negative electric-field intensity is on the plasma torch wall and gradually decreases away from the wall. The maximum electric-field intensity in case 2 is greater than that in case 1. This is because a negative electric field is generated by the current of the external induction coil; it is maximum near the wall and minimum near the central axis. The positive electric field generated by the plasma is the largest near the central axis;it can offset part of the negative electric field.The electric field intensity in case 2 is greater than that in case 1,because the input power of the coil in case 2 is greater than that in case 1.

    Fig.6. Distribution of electron number density(upper half)and imaginary part of electric-field intensity(lower half)for case 1.

    Fig.7. Distributionof electron numberdensity(upper half)and imaginary partof electric-field(lower half)for case2.

    Figures 8 and 9 show the distribution of the Joule heating rate and radial Lorentz force in the plasma torch for case 1 and case 2,respectively. The electromagnetic field and flow field of the wind tunnel are coupled and connected through the Joule heating rate and Lorentz force. First,it is observed that the distribution shape of the Joule heating rate is similar to that of radial Lorentz force,and the maximum value positions of Joule heating rate are also below the third and fourth coil.This proves that the position of the Joule heating phenomenon is mainly governed by the radial Lorentz force. Second, according to the formula, the Joule heating rate is directly proportional to the electric field intensity. Figures 6 and 7 show that the maximum electric field intensity is close to the wall,

    whereas the maximum Joule heating rate is not on the wall but approximately 40 mm away from the wall. This is because a cooling device is used for the wall of the plasma torch,thereby reducing the Joule heating rate. The power in case 2 is greater than that in case 1. According to the proportional relationship between power and electric field intensity,the greater the electric power,the greater the electric field intensity and the Joule heating rate are. The maximum Joule heating rate in cases 1 and 2 are 7.0×107W/m3and 8.2×107W/m3,respectively.

    Fig. 8. Distribution of Joule heating rate (upper half) and radial Lorentz force(lower half)for case 1.

    Fig. 9. Distribution of Joule heating rate (upper half) and radial Lorentz force(lower half)for case 2.

    Figures 10 and 11 show the streamline and velocity distribution in the entire flow field for case 1 and case 2,respectively;it can be observed that the maximum velocity in case 1 and case 2 are 110 m/s and 129 m/s, respectively. The maximum plasma flow velocity is found to be at the torch outlet.The maximum of flow velocity is on the central axis, and its minimum is on the test chamber wall. In the plasma torch, a strong vortex can be observed from the gas inlet to the second coils in Fig. 11. Such one or two recirculation flows are frequently found in the ICP simulations.[22,26,41]It has been reported that this vortex results mainly from the effect of the Lorentz force and the negative pressure gradient[22]generated in the coil region (see Figs. 12 and 14 below). On the one hand,the radial and axial component of the Lorentz force can affect corresponding momentum transfer of charged species significantly for a plasma flow.The strong radial Lorentz force strengthens the pinch effect,which can narrow the radial electrical channel of electrons and squeeze the plasma flow in the radial direction in the coil region.[22]On the other hand,both the negative axial Lorentz force and the negative pressure gradient beneath the inductive coils act as resistive forces to prevent the plasma flow from moving forward smoothly. As a result, some of the plasma elements move backwards and form a vortex near the center axis in the coil region.

    Fig. 10. Distribution of streamlines (upper half) and velocity (lower half)for case 1.

    Fig. 11. Distribution of streamlines (upper half) and velocity (lower half)for case 2.

    In addition, near the wall of the test chamber, there are several notable vortexes. Their formation and intensity are affected mainly by the pressure characteristics of the spacious vacuum chamber. First,due to the high speed plasma flowing along the center axis from the torch outlet to the test chamber outlet, it will results in pressure differences between the chamber wall and the outer boundary of the plasma column in the spacious test chamber. The local pressure difference in a spacious chamber easily results in recirculation flow in the chamber. While such recirculation flows near the chamber wall hardly affect the flow properties of the plasma column near the center axis due to the far distance.

    To check the correctness of the velocity simulated in our study,we consult several relevant references.[12,25,42]Approximate values of the simulated velocity are found for the high power ICPs.[25,42]For example, Vasil’evskii and Kolesnikov[25]simulated an air ICP flow under similar working conditions: the nominal input powerP= 90 kW,f=0.44 MHz,p=10 kPa,m=2.8 g/s. Their simulated maximum velocity at the torch exit is about 150 m/s.[25]Although their input power is lower thanP=160 kW in our case 1,their working pressure and mass flow rate are also lower thanp=15 kPa,m=16 g/s in case 1,respectively.Lu and Feng[12]has reported that the plasma velocity decreases as the pressure and the mass flow rate increase. Thus, for case 1 the simulated maximum velocity (110 m/s) has the same order of magnitude as the one (150 m/s) reported by Vasil’evskii and Kolesnikov.[25]Thus,the maximum plasma velocities obtained in our simulations are reasonable.

    Figures 12–14 show the distributions of the static pressure for case 1 and case 2. The maximum pressure positions are both below the third and fourth coil. The maximum pressures for cases 1 and 2 are 15027.9 Pa and 15076.7 Pa,respectively; thus, the maximum pressure for case 2 is greater than that for case 1. From Fig.14 it follows that the pressure variation trends of these two cases are similar along the center axis,Along the positive direction of thexaxis, the pressure value under the induction coil first increases and then decreases until the pressure value tend to be an identical constant of approximately 1.5×104Pa in the test chamber for these two computational cases. However, for these two cases, their maximum pressure drops in the whole flow field are 27.9 Pa and 76.7 Pa,respectively. They are both small.

    Fig.12. Distribution of static pressure in the whole flow field for case 1.

    Fig.13. Distribution of static pressure in the whole flow field for case 2.

    Fig.14. Comparison of simulated static pressure along center axis between case 1 and case 2.

    Fig.15.Radial distribution of mole fractions for air species at the coil center x=225 mm for case 1.

    Figure 15 shows the radial distributions of mole fractions for air chemical species at the coil centerx=225 mm for case 1. As can be seen from the figure,the atomic N and O are the dominant chemical species beforey ≤40 mm. It indicates that in the vicinity of the central axis,after a series of dissociation and ionization reactions, nitrogen molecules decompose into atomic nitrogen N and ionic nitrogen N+, oxygen molecules also decompose into oxygen atomic O and ionic oxygen O+.The mole fractions of electrons e?and ionic N+are about 0.1,and almost the same. That is to say, the negative charge of electrons just offsets the positive charge of ionic atomic nitrogen,so the simulated plasma flow is electrically quasi-neutral.With the increase of the distance from the central axis, the mole fractions of species N+2,O+,N+,e?,NO+decrease fast aftery ≥15 mm.Molecular nitrogen becomes the main chemical species aftery ≥40 mm.

    3.2. Experimental validation of simulated results

    To validate the correctness of the numerical methods and the simulated results presented in this study, the quantitative comparisons between the simulated and the measured temperatures are performed. Four simulations listed in Table 1 are carried out according to the corresponding experimental study in the 1.2-MW ICP wind tunnel for comparison.

    First, following the spectroscopic experimental study of Cipulloet al.,[24]the air inductive plasma jet produced by the 1.2-MW ICPWT is investigated. Intrusive and nonintrusive measurement technique are combined to characterize the plasma flow at a broad range of pressure and input power. The high speed camera imaging technique is used to study the flow features of the plasma,and a frequency-spatial-domain elaboration procedure is adopted to analyze the measured data. In the experiment, the high speed camera imaging is performed by using a complementary metal–oxide–semiconductor camera to ensure the complete visualization of the plasma flow field. Three spectrometers are used to observe the plasma at three different distances from the torch inlet (x1=665 mm,x2=790 mm,x3=915 mm)along the central axis. The temperatures at these three positions are estimated by fitting the measured spectrum with the theoretical spectra.

    Second, the simulated and measured temperatures are compared at these three positions for cases 1 and 3 (P=160 kW), and cases 2 and 4 (P=200 kW), respectively, in Figs.16 and 17. In these two figures,the curves and the symbols denote the numerical and experimental results, respectively. Additionally,because the error range for cases 1 and 3 is not obtained in the experimental study,[24]none of the error bars of the experimental points are depicted in Fig.16.

    Third, as shown in Fig. 16, it can be observed that the temperature value sharply increases to a maximum and then slightly decreases in the plasma torch. Due to the strong external electromagnetic field generated by the electrified coil,the working gas entering from the inlet is heated rapidly by the electrified coil. So the gas temperature rises rapidly in the coil region. The simulated maximum temperatures for cases 1 and 3 are 10514 K and 10693 K, respectively. In addition,the measured (red circle and blue delta) and the simulated(curves) temperatures in Fig. 16 decrease with the increase of the torch outlet distance due to the interaction between the plasma and its surrounding cold gas. In the vacuum chamber,the simulated temperatures accord well with the temperatures measured at the three positions. For example, at the positionx3=915 mm, the experimental and numerical temperatures for case 1 are 6108 K and 5941 K, respectively. The relative error between them is 2.7%.

    Fig. 16. Comparison between simulated and the measured temperatures along center axis for cases 1 and 3(P=160 kW).

    Additionally,as can be seen in Fig.16,the simulated temperatures show a slightly increasing tendency near the chamber outlet and in a region of 800 mm≤x ≤880 mm for case 1 (p=15 kPa) and case 3 (p=10 kPa), respectively. This abrupt temperature fluctuation might be caused by uncertainties of chemical reactions of air in the simulation,because similar abrupt temperature increment and fluctuation were also found in the arc-heated air plasma simulation.[43]However,the fluctuation of the temperature simulated in our study is about 200 K. The relative error of the simulated temperature caused by this fluctuation is within 3.5%. So the effect of this temperature fluctuation might be ignorable. To make clear the definite reasons for this temperature fluctuation,this issue will be studied further in the future.

    Fig. 17. Comparison between simulated and measured temperatures along center axis for cases 2 and 4(P=200 kW).

    Finally,as shown in Fig.17,generally the simulated and measured temperatures also accord well with each other in the test chamber. Although a few deviations can be seen at the positionx3=915 mm for the 10-kPa case(blue curve and blue symbol delta), the maximum relative error is 6.8%, because the measured and the simulated temperatures are 6465 K and

    6024 K respectively.In summary,the maximum relative errors between the simulated and measured results for cases 1–4 are all within 6.8%, and thus indicating that the numerical methods and results presented in this study are reliable and suitable.

    4. Conclusions

    In this study, the basic flow and electromagnetic properties of a 1.2-MW high power ICP wind tunnel,such as its distribution of the temperature,pressure,electronic density,electric field intensity, Lorentz force, Joule heating rate,etc., are obtained by solving the multiphysics coupling fluid dynamic equations. The simulated temperatures are compared with the experimentally measured ones,and they accord well with each other. So the numerical methods used in this study are reliable and suitable.

    Through the present study, it is also found that the maximum values of the temperature,pressure,electric field intensity, Lorentz force, Joule heating rate, and speed are all generated near the center of the induction coil. The temperature increases to the maximum value in the plasma torch,and then decreases with the increase in the distance from the torch outlet. The pressure value is kept basically the same in the test chamber,but slightly changes in the plasma torch. The change trend of the pressure is similar to that of the translational temperature. In addition, it is also observed that the maximum value of the radial Lorentz force is practically twice or thrice bigger than that that of the axial the Lorentz force. As a result,the momentum transfer in the radial direction is stronger.Moreover,the movement of electrons is also strongly affected by the radial Lorentz force pointing to the central axis. The electron number density and positive electric field intensity generated by the electron are largest at the central axis. Consequently,the negative electric-field intensity generated by the coil current is partially offset by the positive electric field intensity generated by electrons. At the entrance of the torch and near the wall of the test chamber,it is found that the eddy flows are easily produced.

    Acknowledgment

    All simulations were performed by using the Tianhe-2 Supercomputer at the National Supercomputer Center in Guangzhou,China.

    嫩草影院入口| 建设人人有责人人尽责人人享有的 | 国产精品久久久久久久电影| 亚洲综合精品二区| 亚洲成人精品中文字幕电影| .国产精品久久| 久久久久免费精品人妻一区二区| 青春草视频在线免费观看| 国产精品永久免费网站| 女的被弄到高潮叫床怎么办| 好男人在线观看高清免费视频| 欧美性猛交╳xxx乱大交人| 视频中文字幕在线观看| 国产精品三级大全| 国产黄a三级三级三级人| 亚洲成人中文字幕在线播放| 婷婷色麻豆天堂久久 | 18禁在线无遮挡免费观看视频| 日本av手机在线免费观看| 狂野欧美激情性xxxx在线观看| or卡值多少钱| 丰满人妻一区二区三区视频av| 在线天堂最新版资源| 国产精品1区2区在线观看.| 伊人久久精品亚洲午夜| 日韩成人伦理影院| 天堂av国产一区二区熟女人妻| 国产色爽女视频免费观看| 国产精品伦人一区二区| 国产乱来视频区| 午夜激情福利司机影院| 中文字幕精品亚洲无线码一区| 老司机影院成人| 最近手机中文字幕大全| 免费观看的影片在线观看| 免费观看人在逋| 综合色av麻豆| 亚洲欧洲日产国产| 欧美一区二区国产精品久久精品| 五月玫瑰六月丁香| 亚洲国产精品专区欧美| 99久久精品热视频| 久久99热这里只有精品18| 色噜噜av男人的天堂激情| 国产真实伦视频高清在线观看| 日韩av在线免费看完整版不卡| 日韩欧美 国产精品| 国产爱豆传媒在线观看| 亚洲国产高清在线一区二区三| www.av在线官网国产| 国产白丝娇喘喷水9色精品| 亚洲在久久综合| 国产精品99久久久久久久久| 欧美潮喷喷水| 成人美女网站在线观看视频| 亚洲av成人精品一区久久| 十八禁国产超污无遮挡网站| 一级爰片在线观看| 99久久人妻综合| 内地一区二区视频在线| 桃色一区二区三区在线观看| 国产麻豆成人av免费视频| 久久久午夜欧美精品| 在线免费观看的www视频| 乱人视频在线观看| 国产高清视频在线观看网站| 99热这里只有是精品50| videos熟女内射| 天堂影院成人在线观看| 国产高清三级在线| 亚洲在线自拍视频| 可以在线观看毛片的网站| 一级毛片aaaaaa免费看小| 十八禁国产超污无遮挡网站| 18禁动态无遮挡网站| 中文字幕制服av| 特大巨黑吊av在线直播| 99热6这里只有精品| 国产淫片久久久久久久久| 一级二级三级毛片免费看| 成人漫画全彩无遮挡| 一个人免费在线观看电影| 久久久久久国产a免费观看| 最近视频中文字幕2019在线8| 久久久久久大精品| 国产精品一区二区三区四区久久| 国产一区亚洲一区在线观看| 变态另类丝袜制服| 精品久久国产蜜桃| 国产欧美另类精品又又久久亚洲欧美| 变态另类丝袜制服| 亚洲最大成人av| 成人综合一区亚洲| 日韩av在线大香蕉| 国产精华一区二区三区| 男人和女人高潮做爰伦理| 成人av在线播放网站| 在线免费十八禁| 亚洲综合精品二区| 久久久久九九精品影院| 男人和女人高潮做爰伦理| 免费观看的影片在线观看| 麻豆成人午夜福利视频| 国产黄色视频一区二区在线观看 | 国语自产精品视频在线第100页| 欧美+日韩+精品| 人人妻人人澡人人爽人人夜夜| 九九在线视频观看精品| 久久人妻熟女aⅴ| 交换朋友夫妻互换小说| 国产精品偷伦视频观看了| 国产一区二区在线观看日韩| xxx大片免费视频| 亚洲,一卡二卡三卡| 成人漫画全彩无遮挡| 午夜免费鲁丝| 少妇猛男粗大的猛烈进出视频| 亚洲精品国产av成人精品| 人妻少妇偷人精品九色| 国产无遮挡羞羞视频在线观看| 亚洲美女搞黄在线观看| 久久狼人影院| 欧美日本中文国产一区发布| 精品国产露脸久久av麻豆| 最近最新中文字幕大全免费视频 | 这个男人来自地球电影免费观看 | 亚洲精品456在线播放app| 久久精品国产a三级三级三级| 久久婷婷青草| 伦精品一区二区三区| 免费高清在线观看日韩| a级毛片在线看网站| 91国产中文字幕| 亚洲国产精品专区欧美| 男女边吃奶边做爰视频| 91久久精品国产一区二区三区| 9色porny在线观看| 嫩草影院入口| av电影中文网址| 日韩一区二区三区影片| 啦啦啦啦在线视频资源| 国产乱人偷精品视频| 成年动漫av网址| 亚洲国产日韩一区二区| 女性生殖器流出的白浆| 国产精品久久久久久久电影| 国产成人精品久久久久久| 欧美精品人与动牲交sv欧美| 亚洲国产精品一区二区三区在线| 国产精品国产三级专区第一集| av播播在线观看一区| 亚洲av成人精品一二三区| av片东京热男人的天堂| 丰满乱子伦码专区| 午夜福利视频精品| 下体分泌物呈黄色| 欧美 日韩 精品 国产| 伊人久久国产一区二区| 欧美变态另类bdsm刘玥| 亚洲成色77777| 日产精品乱码卡一卡2卡三| 亚洲国产av影院在线观看| 久久久久视频综合| 在线观看国产h片| 亚洲,一卡二卡三卡| 亚洲精品第二区| 午夜精品国产一区二区电影| 国产不卡av网站在线观看| 精品熟女少妇av免费看| 午夜日本视频在线| 精品一区二区免费观看| 久久久久人妻精品一区果冻| 日韩欧美精品免费久久| 曰老女人黄片| av线在线观看网站| 男人添女人高潮全过程视频| 不卡视频在线观看欧美| 在线天堂最新版资源| 久久99精品国语久久久| 久久久久视频综合| 久久97久久精品| 国产精品久久久久久精品电影小说| 18禁观看日本| 成年女人在线观看亚洲视频| 欧美人与善性xxx| 久久久国产精品麻豆| a级毛片在线看网站| 精品国产一区二区久久| 精品第一国产精品| 亚洲欧美日韩卡通动漫| 免费黄色在线免费观看| 亚洲精品一区蜜桃| 尾随美女入室| 国产免费福利视频在线观看| 国产精品国产三级专区第一集| 亚洲伊人久久精品综合| 久久ye,这里只有精品| 夜夜骑夜夜射夜夜干| 91在线精品国自产拍蜜月| 久久韩国三级中文字幕| 最近的中文字幕免费完整| 久久精品国产亚洲av天美| 赤兔流量卡办理| 午夜影院在线不卡| 深夜精品福利| 欧美xxⅹ黑人| 国产精品久久久久久久电影| 夫妻性生交免费视频一级片| 欧美最新免费一区二区三区| 黄色一级大片看看| 视频在线观看一区二区三区| 亚洲美女视频黄频| 啦啦啦在线观看免费高清www| 久久青草综合色| 女性被躁到高潮视频| 成人二区视频| 国产色爽女视频免费观看| 亚洲,欧美精品.| 精品国产一区二区三区久久久樱花| 亚洲精品久久久久久婷婷小说| 777米奇影视久久| 久久ye,这里只有精品| 美女内射精品一级片tv| 国产精品嫩草影院av在线观看| 飞空精品影院首页| 亚洲欧美色中文字幕在线| 久久久久精品人妻al黑| av网站免费在线观看视频| 亚洲国产av影院在线观看| 欧美少妇被猛烈插入视频| 美女内射精品一级片tv| 亚洲精品成人av观看孕妇| 亚洲国产成人一精品久久久| 成人午夜精彩视频在线观看| 人体艺术视频欧美日本| 精品国产乱码久久久久久小说| 国产成人精品婷婷| 亚洲色图综合在线观看| 国产老妇伦熟女老妇高清| 欧美 日韩 精品 国产| 日韩电影二区| videosex国产| 日本与韩国留学比较| 精品国产一区二区久久| 成人影院久久| 亚洲婷婷狠狠爱综合网| 少妇人妻精品综合一区二区| 亚洲av在线观看美女高潮| 欧美激情国产日韩精品一区| 男女免费视频国产| 日本vs欧美在线观看视频| 日本av手机在线免费观看| 2021少妇久久久久久久久久久| 欧美成人午夜精品| 久久狼人影院| 精品国产一区二区三区久久久樱花| 国产一区二区在线观看日韩| 在线观看国产h片| 欧美 亚洲 国产 日韩一| 韩国精品一区二区三区 | 国产男人的电影天堂91| 精品久久久精品久久久| 久久精品国产a三级三级三级| 日韩熟女老妇一区二区性免费视频| 精品国产乱码久久久久久小说| 国产精品女同一区二区软件| 日韩 亚洲 欧美在线| 亚洲三级黄色毛片| 一级毛片黄色毛片免费观看视频| 亚洲熟女精品中文字幕| 青青草视频在线视频观看| 久久久精品区二区三区| 精品人妻在线不人妻| 中文天堂在线官网| 18禁裸乳无遮挡动漫免费视频| 啦啦啦视频在线资源免费观看| 亚洲高清免费不卡视频| 老熟女久久久| 亚洲,欧美,日韩| 午夜福利网站1000一区二区三区| 亚洲第一区二区三区不卡| 日韩制服骚丝袜av| 爱豆传媒免费全集在线观看| av一本久久久久| 麻豆乱淫一区二区| 国产毛片在线视频| 国产黄频视频在线观看| 国产日韩一区二区三区精品不卡| 欧美xxxx性猛交bbbb| 另类亚洲欧美激情| 欧美日韩成人在线一区二区| 97在线视频观看| 久久久久精品性色| 少妇精品久久久久久久| 你懂的网址亚洲精品在线观看| 草草在线视频免费看| 一二三四中文在线观看免费高清| 久久久久久久国产电影| 国产免费现黄频在线看| 国产极品粉嫩免费观看在线| 黄色怎么调成土黄色| 超碰97精品在线观看| 日韩,欧美,国产一区二区三区| 波多野结衣一区麻豆| 日本爱情动作片www.在线观看| 在线免费观看不下载黄p国产| 国产av一区二区精品久久| 99精国产麻豆久久婷婷| 日韩制服骚丝袜av| 99九九在线精品视频| 亚洲国产毛片av蜜桃av| 国产女主播在线喷水免费视频网站| 高清av免费在线| 一本大道久久a久久精品| 欧美精品高潮呻吟av久久| 午夜福利在线观看免费完整高清在| 国产一区二区在线观看av| av在线播放精品| av有码第一页| 免费人妻精品一区二区三区视频| 国产视频首页在线观看| 97人妻天天添夜夜摸| 少妇的逼水好多| 亚洲精品久久久久久婷婷小说| 免费av中文字幕在线| 韩国av在线不卡| 99国产精品免费福利视频| 只有这里有精品99| 91在线精品国自产拍蜜月| 内地一区二区视频在线| 精品国产露脸久久av麻豆| 精品国产一区二区久久| 一区二区三区乱码不卡18| 日本91视频免费播放| 在线看a的网站| 夫妻午夜视频| 一边摸一边做爽爽视频免费| 免费在线观看完整版高清| 亚洲色图 男人天堂 中文字幕 | 日韩 亚洲 欧美在线| kizo精华| 国产国语露脸激情在线看| 最近手机中文字幕大全| 午夜影院在线不卡| 22中文网久久字幕| 日本av手机在线免费观看| 丝瓜视频免费看黄片| 免费女性裸体啪啪无遮挡网站| 免费久久久久久久精品成人欧美视频 | 亚洲国产av影院在线观看| 国产免费一区二区三区四区乱码| 在线观看人妻少妇| 交换朋友夫妻互换小说| 两个人看的免费小视频| 久久久久久久大尺度免费视频| 全区人妻精品视频| 少妇的逼好多水| 日韩 亚洲 欧美在线| 欧美人与善性xxx| 精品一区二区免费观看| a级毛片黄视频| 欧美激情国产日韩精品一区| 久久精品久久久久久噜噜老黄| 一本久久精品| 考比视频在线观看| 最新中文字幕久久久久| 咕卡用的链子| 久久久欧美国产精品| 波多野结衣一区麻豆| 美国免费a级毛片| 午夜视频国产福利| 亚洲色图 男人天堂 中文字幕 | 国产精品一二三区在线看| 亚洲中文av在线| 国产精品国产三级国产av玫瑰| 熟女人妻精品中文字幕| 不卡视频在线观看欧美| 国产成人免费观看mmmm| 亚洲欧美中文字幕日韩二区| 18禁国产床啪视频网站| 国产精品秋霞免费鲁丝片| 我要看黄色一级片免费的| 久久这里有精品视频免费| 日本与韩国留学比较| 99久久综合免费| 国产成人精品一,二区| 男女啪啪激烈高潮av片| 看免费av毛片| 涩涩av久久男人的天堂| 欧美精品一区二区大全| 看免费av毛片| 国产色婷婷99| 超色免费av| 成人无遮挡网站| 久久久久久久亚洲中文字幕| 午夜日本视频在线| 国产成人a∨麻豆精品| 国产免费又黄又爽又色| 黄片无遮挡物在线观看| 亚洲综合精品二区| 日韩成人伦理影院| 亚洲国产av新网站| 久久久久精品性色| 亚洲国产日韩一区二区| 草草在线视频免费看| 亚洲精品成人av观看孕妇| 国产精品.久久久| 一级片'在线观看视频| 麻豆乱淫一区二区| 亚洲成人av在线免费| 久久久精品94久久精品| 韩国高清视频一区二区三区| 国产免费现黄频在线看| √禁漫天堂资源中文www| 久久人人97超碰香蕉20202| 永久免费av网站大全| 天美传媒精品一区二区| 丝袜喷水一区| 超色免费av| 日韩免费高清中文字幕av| 国产午夜精品一二区理论片| 国产麻豆69| 久久精品夜色国产| 美女脱内裤让男人舔精品视频| 午夜久久久在线观看| 午夜精品国产一区二区电影| 最近中文字幕2019免费版| 日韩av免费高清视频| 99热全是精品| 90打野战视频偷拍视频| 国产激情久久老熟女| 天堂俺去俺来也www色官网| 夫妻性生交免费视频一级片| 亚洲国产欧美在线一区| 亚洲国产欧美日韩在线播放| 99视频精品全部免费 在线| 国产亚洲欧美精品永久| 欧美3d第一页| 久久精品人人爽人人爽视色| 亚洲内射少妇av| 亚洲av综合色区一区| 国产高清不卡午夜福利| 亚洲内射少妇av| 99热这里只有是精品在线观看| 国产精品.久久久| 亚洲综合色网址| 亚洲精品第二区| 免费大片黄手机在线观看| 七月丁香在线播放| 亚洲人与动物交配视频| 日韩在线高清观看一区二区三区| 哪个播放器可以免费观看大片| 丝袜美足系列| 国产色婷婷99| 街头女战士在线观看网站| 咕卡用的链子| 另类亚洲欧美激情| 成人国产麻豆网| 国产成人aa在线观看| 免费看光身美女| 少妇的丰满在线观看| 亚洲欧美一区二区三区国产| 免费人成在线观看视频色| 欧美97在线视频| 最近中文字幕2019免费版| 日本午夜av视频| 久久久久人妻精品一区果冻| 热99国产精品久久久久久7| 免费黄网站久久成人精品| 岛国毛片在线播放| 91精品国产国语对白视频| 丰满迷人的少妇在线观看| 午夜福利视频在线观看免费| 99香蕉大伊视频| 99精国产麻豆久久婷婷| 黄片播放在线免费| 秋霞伦理黄片| 日韩中文字幕视频在线看片| 女的被弄到高潮叫床怎么办| 视频在线观看一区二区三区| av电影中文网址| 久久久国产精品麻豆| 自拍欧美九色日韩亚洲蝌蚪91| 久久精品久久久久久久性| 少妇高潮的动态图| 日韩制服骚丝袜av| 精品少妇内射三级| 狠狠精品人妻久久久久久综合| 精品99又大又爽又粗少妇毛片| 国产一级毛片在线| 久久精品国产鲁丝片午夜精品| 纯流量卡能插随身wifi吗| 美女脱内裤让男人舔精品视频| 精品第一国产精品| 少妇被粗大的猛进出69影院 | 久久热在线av| h视频一区二区三区| 亚洲精品色激情综合| 欧美成人午夜免费资源| 国产精品一国产av| 午夜激情av网站| 永久网站在线| 99久久精品国产国产毛片| 国产av码专区亚洲av| 嫩草影院入口| 大香蕉久久网| 久久青草综合色| kizo精华| 精品一区在线观看国产| 精品国产一区二区三区四区第35| 桃花免费在线播放| 纯流量卡能插随身wifi吗| 欧美人与性动交α欧美精品济南到 | 免费观看在线日韩| 午夜91福利影院| 90打野战视频偷拍视频| 亚洲 欧美一区二区三区| 免费看不卡的av| 一边亲一边摸免费视频| 欧美精品高潮呻吟av久久| 亚洲成国产人片在线观看| 乱人伦中国视频| 久久精品久久久久久噜噜老黄| 香蕉精品网在线| 69精品国产乱码久久久| 免费看av在线观看网站| 中文欧美无线码| 丰满饥渴人妻一区二区三| 欧美亚洲 丝袜 人妻 在线| 97精品久久久久久久久久精品| 天天操日日干夜夜撸| 嫩草影院入口| 男女无遮挡免费网站观看| 亚洲伊人久久精品综合| 纯流量卡能插随身wifi吗| 婷婷色综合大香蕉| 国产无遮挡羞羞视频在线观看| 国产女主播在线喷水免费视频网站| 看免费av毛片| 欧美老熟妇乱子伦牲交| 高清av免费在线| 18禁动态无遮挡网站| 2022亚洲国产成人精品| 少妇熟女欧美另类| 2022亚洲国产成人精品| 国产日韩欧美视频二区| 一级,二级,三级黄色视频| 欧美xxxx性猛交bbbb| 亚洲av中文av极速乱| 国产乱来视频区| 国产精品无大码| 丰满饥渴人妻一区二区三| 只有这里有精品99| 欧美3d第一页| av免费观看日本| 欧美亚洲日本最大视频资源| 哪个播放器可以免费观看大片| 啦啦啦啦在线视频资源| 日本午夜av视频| 久久久久久久国产电影| 精品国产一区二区久久| av在线app专区| 久久久久久伊人网av| 免费黄频网站在线观看国产| 黑人高潮一二区| 欧美97在线视频| 国产激情久久老熟女| 久久这里有精品视频免费| 黑人猛操日本美女一级片| 免费看不卡的av| 黑人欧美特级aaaaaa片| 亚洲成av片中文字幕在线观看 | 中文字幕人妻熟女乱码| 纯流量卡能插随身wifi吗| 免费播放大片免费观看视频在线观看| 久久人人97超碰香蕉20202| 国产一区二区在线观看av| 涩涩av久久男人的天堂| 大码成人一级视频| 久久婷婷青草| 狂野欧美激情性xxxx在线观看| 亚洲性久久影院| 看免费av毛片| 久久女婷五月综合色啪小说| 亚洲综合精品二区| 午夜福利网站1000一区二区三区| 最近2019中文字幕mv第一页| 91aial.com中文字幕在线观看| 晚上一个人看的免费电影| 国产精品秋霞免费鲁丝片| tube8黄色片| 少妇的逼好多水| av视频免费观看在线观看| 人人澡人人妻人| 深夜精品福利| 热re99久久国产66热| 国产有黄有色有爽视频| 春色校园在线视频观看| 日本黄大片高清| 一二三四在线观看免费中文在 | 久久久国产一区二区| 少妇的逼水好多| 日本av手机在线免费观看| 亚洲经典国产精华液单| 大片免费播放器 马上看| 久久人妻熟女aⅴ| 国产精品99久久99久久久不卡 | 亚洲天堂av无毛| 九九爱精品视频在线观看| 亚洲精品乱久久久久久| 久热这里只有精品99| 免费大片18禁| 色94色欧美一区二区| 国产在视频线精品| 午夜福利在线观看免费完整高清在|