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

    Investigation of the self-induced magnetic field characteristics in a pulsed plasma thruster with flared electrodes

    2020-04-24 05:49:58QiLIU劉祺LeiYANG楊磊YupingHUANG黃玉平XuZHAO趙絮andZaipingZHENG鄭再平
    Plasma Science and Technology 2020年3期
    關(guān)鍵詞:楊磊黃玉

    Qi LIU(劉祺),Lei YANG(楊磊),Yuping HUANG(黃玉平),Xu ZHAO(趙絮)and Zaiping ZHENG (鄭再平)

    Beijing Research Institute of Precise Mechatronics and Controls,Beijing 100076,People’s Republic of China

    1 Author to whom any correspondence should be addressed.

    Abstract

    Keywords:pulsed plasma thruster,numerical investigation,magnetic flux density,flared electrodes

    1.Introduction

    Pulsed plasma thrusters (PPTs) are a kind of electromagnetic accelerator for spacecrafts.Due to their flexible power and robust structure,PPTs have been applied to several space missions to execute orbital and attitude adjustment.In 1964,the PPT in Zond-2 became the first electric propulsion system used in space[1].However,due to the lower performance of the PPT compared with other thrusters such as Hall thrusters and ion thrusters,PPTs were not recommended for space application.Actually,PPT has a typical thrust efficiency of 5%–15%,while other thrusters have an efficiency of 50%[2].In recent years,PPTs were optimized to higher efficiency and were selected as the main propulsion system for the lunar exploration mission BW1(Baden-Württemberg 1)[3].Much effort has been made to improve the performance of PPTs [3–8].It was found by experimental research that some approaches could optimize the performance,such as increasing the aspect ratio [3,9],which is the ratio of the gap between the electrodes and the width of the electrodes,and raising the discharge energy [9,10].So the configuration and the working process of this kind of PPT are worth investigating.Actually,such approaches may promote the thruster performance by influencing the magnetic field in the discharge channel,which could enhance the acceleration of the particles.

    Magnetic field significantly affects the behavior of plasma and plays an important role in the working process of a PPT.The magnetic field might also be responsible for the erosion phenomenon due to its confinement of the ions that bombard the electrodes.When it comes to numerical study,an accurate magnetic field model could be the foundation of the magnetohydrodynamics or particle-in-cell simulations[11,12].This work builds a model to describe the magnetic field in spatial and temporal view for a better understanding of the working process and thruster structural design.

    Several research works have been done in the selfinduced magnetic field in the past twenty years.Koizumi et al measured the magnetic flux density by a magnetic probe.However,in order to prevent the disturbance of the plasma flow,the magnetic probe was placed outside the sidewalls of the PPT,explicitly,and the probe was not in the discharge channel.So only qualitative results were presented [2].Scharlemann measured the magnetic field at three different axial positions in a Teflon-fed PPT and water-fed PPT with a magnetic probe.One of the positions was between the electrodes,the others were at the exit and in the plume.The distribution of the magnetic field between the electrodes was yet to be investigated [13].Lau et al measured the selfinduced magnetic field at several discrete points along the centerline of the discharge channel by inductive miniaturized probes.However,the magnetic field near the electrodes is difficult to acquire due to the fact that it would interfere with the plasma [14].They also presented that the magnetic flux density peaks at the area close to the Teflon surface,and the value is up to 0.82 T at 1300 V discharge voltage,with a total capacitance of 80 μF [5].

    It could be seen that several experimental investigations have been made on the magnetic field of PPTs,however,the insertion of inductive probes disturbs the flow of plasma in the discharge channel due to the fact that the size of the probes cannot be neglected.Considering the pulsed discharge characteristics of the PPT,the magnetic field might not be totally identical in different pulses.The magnetic flux density at different positions could not be measured synchronically only if multiple magnetic probes are applied,which may introduce more disturbance of the plasma.In addition,the probe could be easily destroyed by the particles in high energy in the discharge channel.Considering the facts mentioned above,studying the magnetic field by modeling and calculating is an efficient and economic way.Nawaz et al applied a refined magnetic field model to a PPT called SIMPLEX(Stuttgart impulsing magnetoplasmadynamic thruster for lunar exploration) with parallel electrodes.No plasma sheet was considered in their calculation [15].Nawaz et al then built a more complete model,including the slug model of this thruster,and compared the results with experimental data[16].As mentioned above,increasing the aspect ratio could promote the thruster performance,and the application of the flared electrodes could be a direct way to increase the aspect ratio.As for the flared structure,Krejci et al built an analytical model to assess the electrode configuration of a PPT[17].The electromagnetic impulse bit was calculated as a criterion of performance analysis.It was concluded that the electromagnetic impulse bit peaks when at a flared angle of electrodes is about 25 deg.

    Figure 1.The PPT with flared electrodes.(a) The schematic of the thruster.(b) The side view and computational domain.

    This work would build a magnetic field model and analyze the field quantitatively.Unlike previous research works,the structure of the electrodes would be discussed in terms of the magnetic field.Above all,this work is going to provide suggestions for thruster design and enhance the understanding of the working process of the PPT.

    2.Analytical model

    2.1.The model of the magnetic field

    In order to study the magnetic field in the PPT with flared electrodes,as shown in figure 1,a calculation model is established.In the figure,deis the distance between the parallel parts of the anode and the cathode,h is the gap between the electrodes,which increases as the position x moves downstream,φis the flared angle of the electrodes,dendingis the distance between the electrode endings,and lelectrodeis the total length of the electrodes.

    In this work,the model is built under the following assumptions:the plasma blob is assumed as the plasma sheet,which is perpendicular to the x axis at all times.In other words,the plasma sheet canting is not considered.The mass of the plasma sheet is assumed to increase as it moves downstream[18–20].It is also assumed that only one plasma sheet exists at every moment,and that the mass of the plasma concentrates on the plasma sheet,which is the same assumption appearing in the slug model.The discharge current is acquired from experiments and inputted into the simulation,while the displacement current is disregarded in the calculation of the magnetic field.In addition,the magnetic field induced by the current in the plasma sheet would be ignored in this simulation.Since the thickness of the plasma sheet would be much broader than that of the electrodes’,the induced magnetic flux density would be partly canceled out in the plasma sheet.On the other hand,the magnetic field induced by the current within the plasma sheet does not contribute to the acceleration process [17],so such ignorance may have little influence on the analysis of plasma acceleration via magnetic flux density distribution.However,the existence and the directional movement of the plasma sheet are still admitted.Moreover,the magnetic field in the thruster is simplified as a two-dimensional model in the x-y plane,and some characteristics in the third dimension are not investigated in this work.

    The magnetic field induced by the current could be expressed as

    The magnetic flux density in z direction Bzis the main issue considered in this work:

    where

    μ0is the vacuum permeability,I is the current,which induces the magnetic field,l is the direction vector of the current,eris the direction vector of the vector pointing from the current element to the calculation point,α is the included angle of l and er,r is the distance between the current element and the calculation point,and r0is the perpendicular distance between the current and the calculation point.

    Equation (2) is a concise description of magnetic field because it could be integrated analytically when the current flows in straight carriers:

    Figure 2.The detailed calculation elements in the PPT with parallel electrodes.

    The magnetic flux density could be calculated by the superposition of several parts where the current flows as shown in figure 2.Explicitly,the magnetic field is induced by the current in the two parts,I and II.In figure 2,the two parts of the magnetic field are going to be obtained by integrating the angles from α1to α2,and from β1to β2,respectively.Since only the magnetic flux density in the z-direction is considered in this work,the subscript z is omitted in the following equations.The total magnetic flux density of a traditional PPT with parallel electrodes could be expressed as

    where r1and r2are the perpendicular distance between the electrodes and the calculation point.

    The detailed expression of magnetic field calculated in PPTs with parallel electrodes is presented in equation(6)and in figure 2,however,the field of the flared electrode case is more sophisticated due to the fact that the structure of the flared angles in the electrodes should be considered.The magnetic field would be the superposition of the field induced by the current in the parallel parts and the flared parts of the electrodes,as shown in figure 3.The magnetic flux density induced by the current in the flared part(IV and V in figure 3)and the bending part (III in figure 3) of the electrodes,i.e.,Belectrodes_flaredand BAnode_bending,should be included as follows:

    where Belectrodes_parais the magnetic flux density induced by the current in the parallel part of the electrodes (I and II in figure 3).As shown in figure 3,the integration of equation(5)from α1to α2and from β1to β2are added to be Belectrodes_para:

    Figure 3.The detailed calculation elements in the PPT with flared electrodes.(a) The parallel part and (b) the flared part.

    where r1–r5are the perpendicular distances between the current in elements I–V and the calculation point.

    As mentioned in [18],for ablative PPTs,the mass of the plasma sheet would increase as it traverses,so a dynamic mass model is applied.A quasi-steady state of the mass flow is assumed and thus the mass of the plasma sheet could be derived.In the flow where the magnetic pressure is much higher than the plasma pressure and where the magnetic Reynolds number is high,a magnetosonic condition exists where the flow speed is equal to the local value of Alfvén sound speed[20].Based on the quasi-steady flow assumption and the law of conservation of momentum,the mass of the plasma sheet is estimated as

    where A is the area of the propellant surface,w is the width of the electrode and vcritis the Alfvén critical speed.The area A could be described as the product of the propellant height hp,which is related to de,and the propellant width wp.So themass of the plasma sheet is

    Table 1.Parameters of the validation case.

    The formula of mass calculation contains the information for electrode configuration,which is w and hp,and the measured current I,so it is believed that the model is compatible with the present thruster.

    During the operation of the PPT,the electromagnetic acceleration of the plasma sheet is produced by the Lorentz force,so the x-position of the plasma sheet xp(t) could be calculated by applying Newton’s second law and solving the kinematics equations as follows:

    where vpis the velocity of the plasma sheet.

    The discharge current I(t) is a crucial parameter for the description of the thruster performance,because it is related to the magnetic field and extensively the mass bit.In order to present a magnetic profile that is closed to its actual distribution,the current data is obtained experimentally for thrusters with different geometries.The Rogowski coil is placed around the cathode column of the capacitor.

    Considering the computer memory and the curve resolution,the time step of the magnetic field calculation is set as 10?9s,while the grid size is 2.5×10?4m.For the PPT with a flared angle of 10 deg and an electrode gap of 21 mm,the calculation area is 82.5 mm×33 mm and the total number of grids is 43560.The related electric and structural parameters of the thruster model are described in table 1.

    2.2.Model validation

    In the temporal axis,the calculated magnetic flux density at different positions on the centerline varying with time is compared with the experimental curve (figure 4(a)) in literature[14].The current data is acquired from the experiments of a PPT developed by the Beijing Research Institute of Precise Mechatronics and Controls.The flared angle of the thruster is 20 deg,deis 21 mm,and the capacitance is 80 μF,which are the same as the parameters in [14].As shown in figure 4(b),the maximum of the magnetic flux density at x=6 mm is

    Figure 4.The magnetic flux density at different axial positions on the centerline of the discharge channel.(a) Experimental results in [14](Reprinted from [14],copyright (2014),with permission from Elsevier).(b) Calculated results in this work.

    about 0.76 T,while at the exit of the discharge channel,x=70 mm and 80 mm,the peak values of magnetic flux density descend to 0.1–0.2 T,which coincide with the experimental data pretty well.In the middle area of the channel,the calculated magnetic flux density is higher than the experimental one.As a whole,the curves and the values of the magnetic flux density simulated by the model are in accordance with the experimental results.

    3.Results

    3.1.Magnetic profile

    The discharge current used in calculations was measured by a Rogowski coil experimentally.Thrusters equipped with electrodes in different geometric structures were tested and the discharge curves were acquired.The experimental current data of a thruster,of which the electrode flared angle is 10 deg and deis 21 mm,is applied to the magnetic field model.The magnetic profile is presented in figure 5.It should be mentioned that the magnetic field induced by the current in the bending part of the anode is considered in this calculation,although that mechanical part of structure is out of the computational domain and is not shown in the contour.Other parameters except for φ and deare the same as shown in table 1.It could be seen that the magnetic field flux density decreases from the upstream to the downstream in the discharge channel,ranging from 0.8 T to 0.1 T and below.

    Figure 5.Magnetic profile in the discharge channel.(φ=10 deg,de=21 mm,t=5 μs).

    Figure 6.Experimental discharge current and the magnetic fulx density Bz at different positions on the channel centerline.(φ=10 deg,de=21 mm).

    Since the magnetic field is induced by the discharge current,Bzshould be in close relation to the discharge current.According to equation (5),the current,which could be theoretically depicted as a function

    is a factor in the expression of Bz.In equation(15),parameters a,b and c are constants.As a consequence,the magnetic flux density curves are similar to the curve of a damped oscillation function.However,as shown in figure 6,the peaks of the magnetic flux density curves lag behind the current peak,and the lags are larger at the downstream than the ones at the upstream in the discharge channel,which coincides with the measured results in figure 4(a) [14].The time lags are significant,as shown in table 2.The time lag is due to the fact that there are time-related items,α1and α2,in equation (5).The movement of the plasma sheet is considered and the boundaries of the current in electrodes change at every moment.So at a fixed point in the discharge channel,the magnetic field changes with not only the value of the current in the electrodes but the position of the plasma sheet.The coincidence of the calculated and experimental results toward the time lag demonstrates that the plasma sheet literally moves in the discharge process.

    Figure 7.Bz with and without the magnetic flux density induced by the current in the bending part of the anode on the channel centerline.(t=5 μs,φ=10 deg,de=21 mm).

    Figure 8.Magnetic flux density in thrusters with different electrode gaps and flared angles (x=20 mm,in the midpoint between electrodes).

    Table 2.The time lag between the negative peaks of Bz and I.

    3.2.Magnetic field of different electrode structure

    3.2.1.Anode structure.The structure of the anode has been discussed in [21]from the point of view that the initial inductance could be reduced by shortening the distance between electrode endings,as can be seen in figure 1(b).As a consequence of this improvement,there is an extra part of the anode that is called the bending part in this paper (part III shown in figure 3(a)).From the view of the magnetic field,this anode structure is yet to be studied.The magnetic field is simulated with and without considering the bending part.As shown in figure 7,the magnetic flux density is higher when the magnetic field induced by the bending structure is considered,compared to the one without the bending part,especially at the upstream of the discharge channel ranging from x=0 m to x=0.02 m.It could be seen from the successive image of C+in the APPT firing [12]that the ionization process happens near the surface area of the propellant,then plasma is accelerated from the surface.Sothe increase of magnetic field at the beginning of the discharge channel could enhance the acceleration of the charged particles near the propellant.

    Table 3.The tested cases.

    3.2.2.Geometrical study.Four geometrical cases shown in table 3 are tested and calculated to investigate the influence of deand the flared angle on the magnetic field.

    As can be seen in figure 8,deexerts an evident influence on the magnetic field.For a thruster with a larger de,the greater exposure area of the propellant leads to a higher mass bit.However,the magnetic flux density on the centerline decreases with the increase of de.Consequently,the Lorentz force applied on the charged particles would diminish.

    The flared angle has a slight influence on the magnetic field,which only decreases a little when the flared angle increases from 10 deg to 20 deg.The ablation areas of the propellant are identical for the cases φ10 de21 and φ20 de21,so the mass bit is not affected by the flared angle.A positive effect of the flared angle on the thruster could be deduced that the ion bombardment is alleviated due to the fact that less ions could reach the electrodes with larger flared angles.

    Since the strong electromagnetic acceleration in the PPT,the electromagnetic impulse bit could represent the performance of the PPT,which is depicted as

    The electromagnetic impulse bit in cases φ10 de21 and φ20 de21 are computed as shown in table 4.Although the magnetic field in the case φ20 de21 is slightly lower than the counterpart in case φ10 de21,the electromagnetic impulse bit of the former case is higher than the other one by 0.11 mNs due to the fact that a higher length h is achieved.In a word,the structure of flared electrodes would increase h and promote IEM.

    Figure 9.The current curves of the PPT with φ=10 deg and de=21 mm.(a) Experimental current curve.(b) Simulated current curve.

    Table 4.The electromagnetic impulse bit in the two cases with different flared angles (Ud=1300 V).

    where U is the discharge voltage,R and L are the total resistance and the total inductance of the circuit and plasma,respectively.

    The simulated current curve is in good accordance with the experimental curve as shown in figure 9.

    Several cases of φ=10–40 deg are then simulated.The electromagnetic impulse bit is calculated according to equation (16).In figure 10,it could be seen that the flared angle has a positive effect on the electromagnetic impulse bit when the angle is in the range of 0 to 28 deg.The electromagnetic impulse bit reaches a maximum at φ=28 deg.However,for φ>28 deg,the impulse bit appears to decrease with the increase of the flared angle.Compared with the result of [17],in which the IEMincreases with φ and becomes maximum for φ beyond 25 deg,the

    In the calculation results of the PPT with flared angles of 10 deg and 20 deg,the tendency of IEMincreases when the flared angle is observed.In order to further analyze the effect of the flared angle,more cases of various flared angles are simulated.Due to the numerous cases with different electrode geometries,instead of using the experimental current data,the current is calculated as follows result of a maximum of performance at φ=28 deg in this work could provide a clear suggestion for the structural optimization of the electrodes.The result shows that IEMcould not monotonically increase with the flared angle.

    Figure 10.Change in electromagnetic impulse bit for PPTs with different flared angles (E0=68 J,de=21 mm).

    Actually the electromagnetic impulse bit is a comprehensive result affected by the magnetic field,discharge current and h.The electrode flared angle has a negative effect on the magnetic field,while it has a positive effect on h,which could be seen in figure 11.The absolute value of the slope of Bz,peakcurve increases with the flared angle exactly.In addition,the flared angle has little effect on the current.So it could be deduced that for the electrode flared angle smaller than 28 deg,the positive effect of the flared angle on h prevails over the magnetic field.Whereas,the negative effect of the flared angle on the magnetic field dominates when the angle is larger than 28 deg.Considering the working process of the thruster and the application to the spacecraft,the large flared angle may cause discontinuity and the separation of the plasma sheet and lead to a redundant configuration of the thruster.So the structural design of electrodes could be optimized referring to the simulation result presented above.

    Figure 11.Change in Bz,peak and h for PPTs with different flared angles(E0=68 J,de=21 mm,Bz,peak:the peak value of magnetic flux density at x=20 mm on the channel centerline).

    4.Conclusions

    The self-induced magnetic field in the PPT with flared electrodes is numerically investigated for the purpose of a better understanding of the working process and the structural design of the thruster.A two-dimensional model of the magnetic field in the discharge channel of the PPT is built and validated to be in accordance with the experimental results.The magnetic flux density in the discharge channel is presented and analyzed.The simulated magnetic field flux density decreases from the upstream to the downstream in the discharge channel,ranging from 0.8 T to 0.1 T and below.The peak of the simulated magnetic flux density lags behind the current peak,and the lags are larger at the downstream than the ones at the upstream in the discharge channel,which coincides with the experimental results in literature.By evaluating the expression of the magnetic field and the analytical discharge current,the existence of a moving plasma sheet in the discharge process could be deduced.The magnetic field induced by the current in the extra bending part of the anode is found to enhance the Lorentz force,which acts on the charged particles near the propellant.Finally,the geometric study indicates that the electromagnetic impulse bit could not monotonically increase with the flared angle,and that it reaches the maximum when the flared angle of the electrodes is 28 deg for the conditions in this work,which could provide significant suggestions for structural optimization.

    Acknowledgments

    This work is supported by National Natural Science Foundation of China (No.11602016).

    ORCID iDs

    猜你喜歡
    楊磊黃玉
    三里畈河濱采風(fēng)寄慨
    賞黃花(忘憂草)未開(kāi)花
    Empirical Likehood for Linear Models with Random Designs Under Strong Mixing Samples
    《溪山清遠(yuǎn)》
    一道屏障
    我家是個(gè)調(diào)色盤(pán)
    沉痛悼念黃玉民同志
    An integrated optimization design of a fishing ship hullform at different speeds *
    我的媽媽屬“炮仗”
    搬書(shū)
    久久久国产欧美日韩av| 人人妻人人澡人人爽人人夜夜| 亚洲精品国产av成人精品| 亚州av有码| 久久久久久久久大av| 日韩不卡一区二区三区视频在线| 满18在线观看网站| 国产熟女欧美一区二区| 日本黄色片子视频| 狂野欧美激情性xxxx在线观看| 超色免费av| 久久女婷五月综合色啪小说| 超碰97精品在线观看| 少妇人妻久久综合中文| 亚洲在久久综合| 国产精品久久久久久精品电影小说| 国产一区亚洲一区在线观看| 午夜激情久久久久久久| 边亲边吃奶的免费视频| 国产高清不卡午夜福利| 日韩视频在线欧美| 人妻少妇偷人精品九色| 一级毛片电影观看| 美女大奶头黄色视频| 丁香六月天网| 啦啦啦中文免费视频观看日本| 久久国产精品大桥未久av| 亚洲欧美日韩另类电影网站| 日本与韩国留学比较| 免费黄频网站在线观看国产| 亚洲精品av麻豆狂野| 成人毛片60女人毛片免费| 免费日韩欧美在线观看| 国产免费一级a男人的天堂| 久久久久久久国产电影| 91久久精品国产一区二区三区| 久久99一区二区三区| 丁香六月天网| 成人国产av品久久久| 亚洲av电影在线观看一区二区三区| 亚洲av福利一区| 精品国产一区二区久久| 国产成人精品久久久久久| 99精国产麻豆久久婷婷| 赤兔流量卡办理| 国产一级毛片在线| 搡老乐熟女国产| 国产av一区二区精品久久| 在线观看美女被高潮喷水网站| 久久亚洲国产成人精品v| 婷婷成人精品国产| 亚洲美女黄色视频免费看| 久久久亚洲精品成人影院| 99精国产麻豆久久婷婷| 国产白丝娇喘喷水9色精品| 久久久久久伊人网av| av免费观看日本| 欧美日韩精品成人综合77777| 欧美日韩精品成人综合77777| 看非洲黑人一级黄片| 亚洲色图 男人天堂 中文字幕 | 国产精品久久久久成人av| 9色porny在线观看| 中国三级夫妇交换| 最近中文字幕2019免费版| 国产精品久久久久成人av| 中国国产av一级| 免费不卡的大黄色大毛片视频在线观看| 午夜av观看不卡| 亚洲精品久久久久久婷婷小说| 天堂中文最新版在线下载| 不卡视频在线观看欧美| kizo精华| 午夜免费观看性视频| 国产深夜福利视频在线观看| 国产黄色视频一区二区在线观看| 亚洲国产最新在线播放| 国产av国产精品国产| 国产成人av激情在线播放 | 91成人精品电影| 日本黄色日本黄色录像| 国产永久视频网站| 草草在线视频免费看| 日韩强制内射视频| 国产 精品1| 国产乱人偷精品视频| 午夜福利,免费看| 我的老师免费观看完整版| 国产一区二区在线观看日韩| 久久综合国产亚洲精品| 一级黄片播放器| 天天影视国产精品| 欧美日本中文国产一区发布| av播播在线观看一区| 一边亲一边摸免费视频| 最近手机中文字幕大全| 日日啪夜夜爽| 国产成人精品在线电影| 久久精品国产自在天天线| 欧美少妇被猛烈插入视频| 一区二区三区免费毛片| 内地一区二区视频在线| 岛国毛片在线播放| 少妇被粗大猛烈的视频| 久久人人爽av亚洲精品天堂| 丰满饥渴人妻一区二区三| 中文字幕久久专区| 肉色欧美久久久久久久蜜桃| 3wmmmm亚洲av在线观看| 亚洲av电影在线观看一区二区三区| a级毛片黄视频| 亚洲av在线观看美女高潮| 国产毛片在线视频| 欧美日韩成人在线一区二区| 国产精品免费大片| 国产亚洲精品久久久com| 国产色婷婷99| 国产熟女欧美一区二区| 国产男人的电影天堂91| 久久国内精品自在自线图片| 蜜桃在线观看..| 国产成人精品无人区| 日本vs欧美在线观看视频| 最近中文字幕高清免费大全6| 老司机亚洲免费影院| 久久久精品免费免费高清| 国产有黄有色有爽视频| 最近最新中文字幕免费大全7| 国产综合精华液| 精品一区二区三区视频在线| 精品人妻熟女毛片av久久网站| 插逼视频在线观看| 色吧在线观看| 插阴视频在线观看视频| 久久精品国产鲁丝片午夜精品| 在线播放无遮挡| av网站免费在线观看视频| 日韩制服骚丝袜av| 亚洲精品av麻豆狂野| 伊人久久国产一区二区| 国产精品一区二区在线观看99| 欧美人与性动交α欧美精品济南到 | 国模一区二区三区四区视频| 亚洲少妇的诱惑av| 女人久久www免费人成看片| 如日韩欧美国产精品一区二区三区 | 九草在线视频观看| av播播在线观看一区| 亚洲成人手机| av在线老鸭窝| 人成视频在线观看免费观看| 少妇的逼好多水| 美女cb高潮喷水在线观看| 丝瓜视频免费看黄片| 不卡视频在线观看欧美| 中文欧美无线码| 免费av中文字幕在线| 人妻少妇偷人精品九色| 成人综合一区亚洲| 最近中文字幕2019免费版| 国产极品天堂在线| 国产高清国产精品国产三级| 简卡轻食公司| 亚洲人成网站在线观看播放| 男女啪啪激烈高潮av片| 18禁动态无遮挡网站| 成人免费观看视频高清| 看非洲黑人一级黄片| av不卡在线播放| 高清毛片免费看| 国产在线视频一区二区| 国产精品秋霞免费鲁丝片| 99久国产av精品国产电影| 国产欧美日韩一区二区三区在线 | 亚洲av.av天堂| 精品一区二区三区视频在线| 免费人成在线观看视频色| 中文天堂在线官网| 能在线免费看毛片的网站| 国产深夜福利视频在线观看| 午夜福利网站1000一区二区三区| √禁漫天堂资源中文www| 国内精品宾馆在线| kizo精华| 精品国产一区二区久久| 日韩免费高清中文字幕av| 亚洲精品一二三| 国产精品国产三级专区第一集| 国产高清国产精品国产三级| 国产视频内射| 成年美女黄网站色视频大全免费 | 日韩视频在线欧美| 日韩av不卡免费在线播放| 九色成人免费人妻av| 亚洲欧美中文字幕日韩二区| 国产精品人妻久久久久久| 精品国产一区二区久久| 国语对白做爰xxxⅹ性视频网站| 丁香六月天网| 日本-黄色视频高清免费观看| av免费在线看不卡| 中文字幕亚洲精品专区| 日韩,欧美,国产一区二区三区| 亚洲色图 男人天堂 中文字幕 | 赤兔流量卡办理| 免费观看无遮挡的男女| 少妇猛男粗大的猛烈进出视频| 亚洲色图综合在线观看| 青春草国产在线视频| 日本欧美视频一区| 少妇高潮的动态图| 国产永久视频网站| 成人国产av品久久久| 又黄又爽又刺激的免费视频.| 日本91视频免费播放| 久久毛片免费看一区二区三区| 在线精品无人区一区二区三| 大陆偷拍与自拍| 国产精品偷伦视频观看了| 日韩av免费高清视频| 亚洲色图综合在线观看| 国产黄色视频一区二区在线观看| 久久久久久人妻| 日韩成人伦理影院| 免费高清在线观看日韩| 永久免费av网站大全| 日本wwww免费看| 多毛熟女@视频| 丝袜美足系列| 国产精品一区www在线观看| 91精品国产九色| 黄色欧美视频在线观看| 国产精品99久久久久久久久| 国产乱来视频区| 国精品久久久久久国模美| 男男h啪啪无遮挡| 高清不卡的av网站| 特大巨黑吊av在线直播| 亚洲人成网站在线播| 精品人妻在线不人妻| 成年人免费黄色播放视频| 夜夜看夜夜爽夜夜摸| 欧美精品一区二区大全| 欧美三级亚洲精品| 国产熟女午夜一区二区三区 | 欧美人与性动交α欧美精品济南到 | 免费少妇av软件| 亚洲欧美成人精品一区二区| 最近中文字幕2019免费版| 日韩人妻高清精品专区| 午夜av观看不卡| 18禁裸乳无遮挡动漫免费视频| 久久综合国产亚洲精品| 婷婷色综合www| 国产片特级美女逼逼视频| 午夜激情av网站| 亚洲精品一区蜜桃| 一区二区日韩欧美中文字幕 | 少妇人妻久久综合中文| 自拍欧美九色日韩亚洲蝌蚪91| 丰满乱子伦码专区| 热re99久久国产66热| 美女国产高潮福利片在线看| 亚洲精品成人av观看孕妇| 老熟女久久久| 午夜免费观看性视频| 成年人免费黄色播放视频| 成人毛片a级毛片在线播放| 亚洲美女视频黄频| 啦啦啦在线观看免费高清www| 丝袜美足系列| 久久人人爽人人片av| 亚洲美女搞黄在线观看| 中文字幕久久专区| 久久av网站| 亚洲精品aⅴ在线观看| 97超视频在线观看视频| 国产精品嫩草影院av在线观看| tube8黄色片| 亚洲人成77777在线视频| 国产女主播在线喷水免费视频网站| 国产 一区精品| 一本—道久久a久久精品蜜桃钙片| 大片电影免费在线观看免费| 久久午夜福利片| 久久久久国产精品人妻一区二区| 熟妇人妻不卡中文字幕| 秋霞伦理黄片| 国产永久视频网站| 九九久久精品国产亚洲av麻豆| 美女xxoo啪啪120秒动态图| 久久久久精品性色| 欧美xxxx性猛交bbbb| 亚洲av综合色区一区| 少妇被粗大的猛进出69影院 | 久久99热这里只频精品6学生| 黄色视频在线播放观看不卡| 日本欧美视频一区| 天堂中文最新版在线下载| 美女大奶头黄色视频| 亚洲精品日韩av片在线观看| 一本—道久久a久久精品蜜桃钙片| 久久青草综合色| 免费黄色在线免费观看| 乱码一卡2卡4卡精品| 99精国产麻豆久久婷婷| 国产一区二区三区av在线| 亚洲国产精品一区三区| 水蜜桃什么品种好| 免费黄网站久久成人精品| 大香蕉久久网| 久久久午夜欧美精品| 免费久久久久久久精品成人欧美视频 | 少妇人妻久久综合中文| 少妇猛男粗大的猛烈进出视频| 99热6这里只有精品| 精品一区二区三卡| 日本与韩国留学比较| 亚洲av.av天堂| 久久ye,这里只有精品| 男女边吃奶边做爰视频| 亚洲国产色片| 国产乱人偷精品视频| 成年人免费黄色播放视频| 国产极品粉嫩免费观看在线 | 国产精品久久久久成人av| 一个人免费看片子| 精品人妻熟女毛片av久久网站| 91久久精品电影网| 日韩精品有码人妻一区| 天美传媒精品一区二区| 国产精品熟女久久久久浪| 18禁观看日本| 国产伦理片在线播放av一区| 欧美日韩一区二区视频在线观看视频在线| 视频在线观看一区二区三区| 国产精品国产三级国产av玫瑰| 精品卡一卡二卡四卡免费| 天堂俺去俺来也www色官网| 又黄又爽又刺激的免费视频.| 欧美最新免费一区二区三区| 国产黄色免费在线视频| 精品久久久久久电影网| 精品一区在线观看国产| 午夜老司机福利剧场| 欧美3d第一页| 久久狼人影院| 久久热精品热| 日韩中文字幕视频在线看片| 人妻少妇偷人精品九色| 国产日韩欧美在线精品| 伦精品一区二区三区| 熟女电影av网| 久久久久久伊人网av| 一区二区三区乱码不卡18| 久久久久久久国产电影| av视频免费观看在线观看| 久久99一区二区三区| 我的女老师完整版在线观看| 亚洲高清免费不卡视频| 飞空精品影院首页| 国产色爽女视频免费观看| 成年美女黄网站色视频大全免费 | 亚洲成色77777| 免费播放大片免费观看视频在线观看| 国产视频首页在线观看| 一级毛片我不卡| 人妻人人澡人人爽人人| 午夜福利影视在线免费观看| 亚洲美女搞黄在线观看| 伦理电影大哥的女人| 亚洲伊人久久精品综合| 午夜91福利影院| 999精品在线视频| 一本一本久久a久久精品综合妖精 国产伦在线观看视频一区 | 亚洲av综合色区一区| 久久av网站| 欧美精品一区二区大全| 一级二级三级毛片免费看| 久久久久久久大尺度免费视频| 性高湖久久久久久久久免费观看| 久久精品国产自在天天线| 精品久久久精品久久久| 伊人久久国产一区二区| 看十八女毛片水多多多| 久久久午夜欧美精品| 9色porny在线观看| 日本欧美国产在线视频| 久久精品夜色国产| 久久亚洲国产成人精品v| 制服诱惑二区| 久久99热这里只频精品6学生| 欧美精品高潮呻吟av久久| 午夜福利网站1000一区二区三区| 国产 一区精品| 精品少妇久久久久久888优播| 中文字幕制服av| 老司机影院毛片| 日韩不卡一区二区三区视频在线| 亚洲精品久久午夜乱码| 性色av一级| 国产伦理片在线播放av一区| 在线观看免费日韩欧美大片 | 美女cb高潮喷水在线观看| 精品亚洲成a人片在线观看| 免费大片黄手机在线观看| 午夜日本视频在线| 99热网站在线观看| 亚洲情色 制服丝袜| 亚洲欧美日韩另类电影网站| 色5月婷婷丁香| 国产精品秋霞免费鲁丝片| 一本久久精品| 香蕉精品网在线| 人人妻人人澡人人爽人人夜夜| 亚洲美女视频黄频| 亚洲四区av| 亚洲精华国产精华液的使用体验| 国产精品国产三级专区第一集| 欧美性感艳星| 人妻系列 视频| 久久毛片免费看一区二区三区| 国产成人免费观看mmmm| 夜夜骑夜夜射夜夜干| 日韩精品免费视频一区二区三区 | 欧美人与善性xxx| 国产国语露脸激情在线看| 美女xxoo啪啪120秒动态图| 中文乱码字字幕精品一区二区三区| 久久精品国产鲁丝片午夜精品| 99精国产麻豆久久婷婷| 日韩大片免费观看网站| 能在线免费看毛片的网站| 亚洲国产成人一精品久久久| 大香蕉久久网| 菩萨蛮人人尽说江南好唐韦庄| 欧美日韩国产mv在线观看视频| 国产深夜福利视频在线观看| 狠狠婷婷综合久久久久久88av| 国产精品免费大片| 国产成人aa在线观看| 麻豆成人av视频| 制服诱惑二区| 日韩成人av中文字幕在线观看| 高清毛片免费看| 精品国产一区二区三区久久久樱花| 七月丁香在线播放| 久久国产亚洲av麻豆专区| 欧美亚洲日本最大视频资源| 一级毛片黄色毛片免费观看视频| 简卡轻食公司| 亚洲国产毛片av蜜桃av| 国产极品粉嫩免费观看在线 | 国产淫语在线视频| 国产一区二区在线观看av| 妹子高潮喷水视频| 久久久亚洲精品成人影院| 老司机亚洲免费影院| 国产免费一级a男人的天堂| 久久久久网色| 边亲边吃奶的免费视频| 丝袜脚勾引网站| 22中文网久久字幕| 一区二区三区乱码不卡18| 欧美日本中文国产一区发布| 久久综合国产亚洲精品| 日韩欧美一区视频在线观看| 777米奇影视久久| 一本一本综合久久| 一本一本久久a久久精品综合妖精 国产伦在线观看视频一区 | 亚洲av欧美aⅴ国产| 国产日韩一区二区三区精品不卡 | 一区二区三区乱码不卡18| 亚洲精品一二三| 丰满迷人的少妇在线观看| 国产色爽女视频免费观看| 国产精品成人在线| 中文字幕精品免费在线观看视频 | 久久久久久久大尺度免费视频| 精品国产一区二区三区久久久樱花| 亚洲欧洲日产国产| 伦理电影免费视频| 亚洲中文av在线| 黄色毛片三级朝国网站| 亚洲av成人精品一区久久| 18在线观看网站| 亚洲精品久久午夜乱码| 美女视频免费永久观看网站| 国产精品久久久久久久电影| 97在线视频观看| 丝瓜视频免费看黄片| 交换朋友夫妻互换小说| 日本爱情动作片www.在线观看| tube8黄色片| 91aial.com中文字幕在线观看| 一区二区三区四区激情视频| 欧美精品高潮呻吟av久久| 3wmmmm亚洲av在线观看| a级毛片黄视频| 日日摸夜夜添夜夜添av毛片| 一本色道久久久久久精品综合| 97精品久久久久久久久久精品| 久久精品人人爽人人爽视色| 色哟哟·www| a级毛色黄片| 国产av码专区亚洲av| 久久久久人妻精品一区果冻| 亚洲图色成人| 日本av免费视频播放| 在现免费观看毛片| 国产精品不卡视频一区二区| 在线观看三级黄色| 久久久久久人妻| 国产精品国产三级专区第一集| 欧美日韩视频高清一区二区三区二| tube8黄色片| 欧美日韩在线观看h| 亚洲精品国产av成人精品| 91精品三级在线观看| 一边亲一边摸免费视频| 久久精品国产亚洲网站| 最近中文字幕2019免费版| 丰满乱子伦码专区| 中文字幕免费在线视频6| 交换朋友夫妻互换小说| 亚洲色图综合在线观看| 精品人妻在线不人妻| 国产一区二区在线观看日韩| 亚洲精品国产av成人精品| 男女无遮挡免费网站观看| 国产精品人妻久久久影院| 女的被弄到高潮叫床怎么办| 男女啪啪激烈高潮av片| 热99国产精品久久久久久7| www.av在线官网国产| 丝袜美足系列| 日韩,欧美,国产一区二区三区| 久久久精品94久久精品| 男女无遮挡免费网站观看| 七月丁香在线播放| 国产 精品1| 国产精品.久久久| 一边亲一边摸免费视频| 亚洲精品456在线播放app| 国产熟女欧美一区二区| 国国产精品蜜臀av免费| 在线亚洲精品国产二区图片欧美 | 天美传媒精品一区二区| 国产探花极品一区二区| 国产成人av激情在线播放 | 色网站视频免费| 青春草亚洲视频在线观看| 亚洲性久久影院| 久久97久久精品| 欧美日韩一区二区视频在线观看视频在线| 男人爽女人下面视频在线观看| 日韩人妻高清精品专区| 亚洲精品中文字幕在线视频| 黄片无遮挡物在线观看| 国产色婷婷99| 日本猛色少妇xxxxx猛交久久| 80岁老熟妇乱子伦牲交| 黑人猛操日本美女一级片| 欧美最新免费一区二区三区| 一本大道久久a久久精品| 精品人妻偷拍中文字幕| 精品99又大又爽又粗少妇毛片| 久久久久久伊人网av| 天天影视国产精品| 亚洲国产精品专区欧美| 麻豆成人av视频| 卡戴珊不雅视频在线播放| 99视频精品全部免费 在线| 久久99热6这里只有精品| 国产极品粉嫩免费观看在线 | 久久毛片免费看一区二区三区| 国产视频首页在线观看| 99国产综合亚洲精品| 简卡轻食公司| 亚洲精品亚洲一区二区| 国产无遮挡羞羞视频在线观看| 国产极品粉嫩免费观看在线 | 午夜影院在线不卡| xxx大片免费视频| 午夜视频国产福利| 国产男女内射视频| 高清视频免费观看一区二区| 国产免费视频播放在线视频| 亚洲在久久综合| 久久人人爽av亚洲精品天堂| 国产一区二区三区综合在线观看 | 一级二级三级毛片免费看| 亚洲国产精品999| 亚洲精品av麻豆狂野| 国产高清不卡午夜福利| 日本欧美国产在线视频| 日韩av不卡免费在线播放| 如何舔出高潮| 日本vs欧美在线观看视频| 青春草视频在线免费观看| 久久午夜综合久久蜜桃| 国产在线免费精品| 国产在视频线精品| 久久免费观看电影| 韩国高清视频一区二区三区| 五月伊人婷婷丁香| 老熟女久久久| 国产精品人妻久久久影院| 亚洲图色成人| 日韩大片免费观看网站| 午夜影院在线不卡| 五月伊人婷婷丁香| 欧美丝袜亚洲另类| 国产精品人妻久久久影院| 色94色欧美一区二区|