• <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| 高清黄色对白视频在线免费看 | 成人18禁高潮啪啪吃奶动态图 | 欧美一级a爱片免费观看看| 久久精品夜色国产| 午夜91福利影院| 嫩草影院新地址| 成人特级av手机在线观看| 六月丁香七月| 91精品国产国语对白视频| 国产精品久久久久成人av| av卡一久久| 久久97久久精品| 91在线精品国自产拍蜜月| 一区二区三区四区激情视频| 亚洲精品日本国产第一区| 亚洲精品国产av蜜桃| av在线老鸭窝| 久久久欧美国产精品| 五月开心婷婷网| 国产成人精品无人区| 欧美精品一区二区免费开放| 亚洲自偷自拍三级| 精品少妇久久久久久888优播| 一区二区三区乱码不卡18| 最近的中文字幕免费完整| 亚洲国产成人一精品久久久| 一级毛片aaaaaa免费看小| 亚洲性久久影院| 国产精品99久久久久久久久| 成人亚洲精品一区在线观看| 人妻夜夜爽99麻豆av| 亚洲av国产av综合av卡| av免费观看日本| 桃花免费在线播放| 偷拍熟女少妇极品色| 男男h啪啪无遮挡| 国产精品熟女久久久久浪| 水蜜桃什么品种好| av在线老鸭窝| 国产精品99久久久久久久久| av网站免费在线观看视频| 男女啪啪激烈高潮av片| 日韩大片免费观看网站| 免费大片18禁| 国产精品熟女久久久久浪| 久久精品久久久久久久性| 99久久人妻综合| 国产欧美日韩精品一区二区| 久久99一区二区三区| 69精品国产乱码久久久| 国产高清国产精品国产三级| 国产在视频线精品| 久久精品久久久久久久性| 日本91视频免费播放| 日日摸夜夜添夜夜添av毛片| 日韩欧美 国产精品| 如日韩欧美国产精品一区二区三区 | 精品一区二区免费观看| 国产国拍精品亚洲av在线观看| 一级毛片久久久久久久久女| av女优亚洲男人天堂| 免费看日本二区| 在线观看免费视频网站a站| 熟女av电影| 老司机影院毛片| freevideosex欧美| 欧美少妇被猛烈插入视频| 成人无遮挡网站| 亚洲欧美一区二区三区国产| 涩涩av久久男人的天堂| 女人久久www免费人成看片| 国产精品久久久久久久电影| 人人妻人人澡人人看| 伦理电影免费视频| 国产片特级美女逼逼视频| 国产一级毛片在线| 青春草视频在线免费观看| 少妇的逼水好多| 色94色欧美一区二区| 伊人亚洲综合成人网| 国产午夜精品久久久久久一区二区三区| 国产极品粉嫩免费观看在线 | 天天躁夜夜躁狠狠久久av| av福利片在线| 99九九线精品视频在线观看视频| 99九九在线精品视频 | 成人毛片60女人毛片免费| 黄片无遮挡物在线观看| 国产精品人妻久久久影院| 视频中文字幕在线观看| a 毛片基地| 在线免费观看不下载黄p国产| 狠狠精品人妻久久久久久综合| 国产熟女午夜一区二区三区 | 免费观看性生交大片5| 国产片特级美女逼逼视频| 女性被躁到高潮视频| 久久99热6这里只有精品| 国产色爽女视频免费观看| tube8黄色片| 欧美激情国产日韩精品一区| 亚洲精品视频女| 人人妻人人澡人人看| 久久人人爽人人爽人人片va| 蜜臀久久99精品久久宅男| 在线观看免费视频网站a站| 熟妇人妻不卡中文字幕| 欧美3d第一页| 欧美xxxx性猛交bbbb| 午夜免费鲁丝| 97超视频在线观看视频| 国产精品秋霞免费鲁丝片| 久久久久久久久久人人人人人人| 亚洲欧美成人精品一区二区| 精品亚洲成国产av| 九九爱精品视频在线观看| 亚洲av在线观看美女高潮| 纯流量卡能插随身wifi吗| 80岁老熟妇乱子伦牲交| 91精品一卡2卡3卡4卡| 欧美xxⅹ黑人| 国产欧美另类精品又又久久亚洲欧美| 热99国产精品久久久久久7| 熟妇人妻不卡中文字幕| 久久久久久久久久久免费av| 日韩精品有码人妻一区| 97超视频在线观看视频| 久久国产亚洲av麻豆专区| 三级国产精品片| 国产免费又黄又爽又色| 亚洲经典国产精华液单| 日韩伦理黄色片| 在线观看国产h片| 最后的刺客免费高清国语| 欧美3d第一页| 久久久精品免费免费高清| 三级经典国产精品| 亚洲综合色惰| 欧美日韩视频高清一区二区三区二| 一区二区三区乱码不卡18| 99热全是精品| 国产淫语在线视频| 国产成人精品无人区| 精品国产一区二区久久| 亚洲成人手机| 丝袜在线中文字幕| 又黄又爽又刺激的免费视频.| 99九九在线精品视频 | 亚洲av二区三区四区| 在线观看av片永久免费下载| 黑人高潮一二区| 亚洲高清免费不卡视频| 美女视频免费永久观看网站| 九草在线视频观看| 99热国产这里只有精品6| 日日撸夜夜添| 精品久久久精品久久久| 欧美成人精品欧美一级黄| 国产成人午夜福利电影在线观看| 国产成人精品无人区| 日韩三级伦理在线观看| 91精品国产国语对白视频| 免费看av在线观看网站| 99国产精品免费福利视频| 亚洲av国产av综合av卡| 成人18禁高潮啪啪吃奶动态图 | 国产精品一区二区在线不卡| 最近的中文字幕免费完整| 一级爰片在线观看| .国产精品久久| 少妇丰满av| 日韩伦理黄色片| 亚洲精品国产av成人精品| 人妻 亚洲 视频| 欧美日韩在线观看h| 91精品国产国语对白视频| 九九爱精品视频在线观看| 制服丝袜香蕉在线| 成人18禁高潮啪啪吃奶动态图 | 三上悠亚av全集在线观看 | 国产精品久久久久久精品电影小说| 18+在线观看网站| 69精品国产乱码久久久| 插阴视频在线观看视频| kizo精华| 日韩三级伦理在线观看| 少妇人妻精品综合一区二区| 久久久久久久久久久免费av| 精品久久久久久久久av| 简卡轻食公司| 日韩成人伦理影院| 亚洲精品日韩在线中文字幕| 丰满乱子伦码专区| 午夜激情福利司机影院| 免费看不卡的av| 国内少妇人妻偷人精品xxx网站| 视频区图区小说| 99九九线精品视频在线观看视频| 亚洲av综合色区一区| 如何舔出高潮| 久久精品国产亚洲av涩爱| 亚洲一区二区三区欧美精品| 夫妻午夜视频| 三级国产精品片| 97超碰精品成人国产| 蜜臀久久99精品久久宅男| 丰满迷人的少妇在线观看| 亚洲人与动物交配视频| 日本爱情动作片www.在线观看| 欧美精品国产亚洲| 成人国产av品久久久| 最近的中文字幕免费完整| 在线亚洲精品国产二区图片欧美 | 美女大奶头黄色视频| 日韩av不卡免费在线播放| 亚洲伊人久久精品综合| 久久久精品94久久精品| 狂野欧美激情性xxxx在线观看| 亚洲婷婷狠狠爱综合网| 亚洲三级黄色毛片| 五月开心婷婷网| 久久久久久久久大av| 久久精品久久精品一区二区三区| 中文在线观看免费www的网站| 99视频精品全部免费 在线| 国产欧美另类精品又又久久亚洲欧美| 水蜜桃什么品种好| 午夜福利,免费看| 街头女战士在线观看网站| 国产黄片视频在线免费观看| 免费看日本二区| 中文资源天堂在线| 欧美xxⅹ黑人| 亚洲欧洲精品一区二区精品久久久 | 一级黄片播放器| 91精品一卡2卡3卡4卡| 少妇人妻一区二区三区视频| 啦啦啦中文免费视频观看日本| 亚洲国产欧美日韩在线播放 | 亚洲精品国产av成人精品| 色5月婷婷丁香| 一级毛片 在线播放| 午夜激情福利司机影院| 日本黄大片高清| 天堂8中文在线网| 国产日韩欧美亚洲二区| 深夜a级毛片| 我要看黄色一级片免费的| 日韩视频在线欧美| 国产伦精品一区二区三区视频9| 午夜福利视频精品| 国产成人精品一,二区| 国产女主播在线喷水免费视频网站| 国产精品免费大片| 亚洲精品成人av观看孕妇| 亚洲无线观看免费| 欧美亚洲 丝袜 人妻 在线| 女人精品久久久久毛片| 色婷婷av一区二区三区视频| 亚洲欧美清纯卡通| 日本91视频免费播放| 欧美精品一区二区大全| 丰满迷人的少妇在线观看| 国产乱来视频区| 99热这里只有是精品50| av播播在线观看一区| 国产av国产精品国产| av黄色大香蕉| 美女大奶头黄色视频| 久久人人爽av亚洲精品天堂| 女的被弄到高潮叫床怎么办| 久久国产精品大桥未久av | 久久青草综合色| 涩涩av久久男人的天堂| 日韩人妻高清精品专区| 国产亚洲最大av| 五月玫瑰六月丁香| 精品一区二区三区视频在线| 亚洲精华国产精华液的使用体验| 全区人妻精品视频| 亚洲怡红院男人天堂| 亚洲综合精品二区| 国产亚洲精品久久久com| 亚洲国产欧美日韩在线播放 | 亚洲av国产av综合av卡| 久久6这里有精品| 午夜激情福利司机影院| 色网站视频免费| 午夜激情久久久久久久| 妹子高潮喷水视频| 国产av码专区亚洲av| 国产亚洲欧美精品永久| 男人狂女人下面高潮的视频| 精品国产露脸久久av麻豆| 国产男人的电影天堂91| 18禁在线无遮挡免费观看视频| 欧美成人午夜免费资源| 国产欧美亚洲国产| 啦啦啦在线观看免费高清www| 毛片一级片免费看久久久久| 精品一品国产午夜福利视频| 日韩制服骚丝袜av| 亚洲婷婷狠狠爱综合网| 一级爰片在线观看| 婷婷色麻豆天堂久久| 亚洲精品久久午夜乱码| 在线观看人妻少妇| av在线观看视频网站免费| 91久久精品国产一区二区成人| 久久国产精品男人的天堂亚洲 | 大片电影免费在线观看免费| 曰老女人黄片| 免费看不卡的av| 国产精品国产三级专区第一集| 卡戴珊不雅视频在线播放| 极品少妇高潮喷水抽搐| 精品国产一区二区三区久久久樱花| 国产精品三级大全| 一级毛片久久久久久久久女| 91成人精品电影| freevideosex欧美| 久久6这里有精品| 亚洲国产精品国产精品| 一级毛片aaaaaa免费看小| 成年美女黄网站色视频大全免费 | www.999成人在线观看| 亚洲一区二区三区欧美精品| 女人久久www免费人成看片| 国产精品影院久久| 人人妻,人人澡人人爽秒播| 日韩视频在线欧美| 黄色视频不卡| 亚洲精品一二三| 国产一区二区激情短视频 | 亚洲第一av免费看| 人人妻,人人澡人人爽秒播| 久久精品亚洲熟妇少妇任你| 国产精品久久久久久精品电影小说| 国产视频一区二区在线看| 婷婷成人精品国产| 超碰成人久久| 黄色 视频免费看| 丝袜在线中文字幕| 在线观看人妻少妇| 免费少妇av软件| 久久久精品免费免费高清| 久热这里只有精品99| 午夜福利免费观看在线| 欧美人与性动交α欧美精品济南到| 国产97色在线日韩免费| 美女扒开内裤让男人捅视频| 亚洲欧洲日产国产| 亚洲精品乱久久久久久| 狂野欧美激情性xxxx| 欧美另类一区| 免费在线观看日本一区| 夜夜骑夜夜射夜夜干| 国产av又大| 久久av网站| 巨乳人妻的诱惑在线观看| 国产一区二区三区综合在线观看| 久久精品人人爽人人爽视色| 欧美日韩亚洲综合一区二区三区_| 亚洲av美国av| 色94色欧美一区二区| 精品国产国语对白av| 色婷婷久久久亚洲欧美| 黄色视频在线播放观看不卡| 午夜两性在线视频| 国产免费av片在线观看野外av| 亚洲精品粉嫩美女一区| 国产在线免费精品| 国产一区二区三区综合在线观看| 一区二区三区激情视频| 九色亚洲精品在线播放| 人人澡人人妻人| 十八禁网站网址无遮挡| 美女高潮喷水抽搐中文字幕| 中文字幕人妻丝袜制服| 免费人妻精品一区二区三区视频| 成人国语在线视频| 一本色道久久久久久精品综合| 久久人人爽av亚洲精品天堂| 午夜激情av网站| a级毛片黄视频| 日本91视频免费播放| 亚洲欧美激情在线| 中文字幕人妻丝袜一区二区| 男女之事视频高清在线观看| 狠狠婷婷综合久久久久久88av| 国产免费福利视频在线观看| 国产淫语在线视频| 91麻豆精品激情在线观看国产 | 妹子高潮喷水视频| 天天躁夜夜躁狠狠躁躁| 男人添女人高潮全过程视频| 国产又爽黄色视频| 亚洲中文av在线| 精品人妻熟女毛片av久久网站| 国产一区二区三区在线臀色熟女 | 飞空精品影院首页| 亚洲专区国产一区二区| 国产成人影院久久av| 少妇的丰满在线观看| 欧美大码av| 国产老妇伦熟女老妇高清| 亚洲五月色婷婷综合| 日韩欧美国产一区二区入口| 国产成人欧美| 在线永久观看黄色视频| 久久中文字幕一级| 狠狠婷婷综合久久久久久88av| 热99re8久久精品国产| 成年女人毛片免费观看观看9 | 一区二区三区激情视频| 欧美日韩精品网址| 久久久欧美国产精品| 亚洲国产欧美一区二区综合| 黄色a级毛片大全视频| 国产精品免费视频内射| 亚洲第一av免费看| 亚洲美女黄色视频免费看| 精品高清国产在线一区| 日韩视频在线欧美| 满18在线观看网站| 制服诱惑二区| 高清av免费在线| 中文字幕精品免费在线观看视频| 亚洲精品一区蜜桃| 99re6热这里在线精品视频| 丝袜人妻中文字幕| 精品久久久精品久久久| 999久久久国产精品视频| 啪啪无遮挡十八禁网站| 国产成人一区二区三区免费视频网站| 久久热在线av| 精品亚洲乱码少妇综合久久| 久久人妻熟女aⅴ| 国产成人啪精品午夜网站| 国产亚洲av片在线观看秒播厂| 国产三级黄色录像| 国产97色在线日韩免费| 亚洲精品av麻豆狂野| 男人添女人高潮全过程视频| 黄色片一级片一级黄色片| 午夜日韩欧美国产| 欧美黑人精品巨大| 中国美女看黄片| 欧美精品av麻豆av| 热re99久久国产66热| 久久热在线av| 少妇 在线观看| 亚洲,欧美精品.| 亚洲欧美成人综合另类久久久| 纵有疾风起免费观看全集完整版| 亚洲伊人久久精品综合| 国产99久久九九免费精品| 啦啦啦免费观看视频1| 丝袜喷水一区| 久久人妻熟女aⅴ| 高清黄色对白视频在线免费看| 黑人欧美特级aaaaaa片| 在线亚洲精品国产二区图片欧美| 精品人妻在线不人妻| 精品国产超薄肉色丝袜足j| 精品人妻在线不人妻| 免费看十八禁软件| 熟女少妇亚洲综合色aaa.| 久久久久久人人人人人| 黄片小视频在线播放| 久久久久久久大尺度免费视频| videos熟女内射| 视频区欧美日本亚洲| 夫妻午夜视频| 黄色视频在线播放观看不卡| 国产淫语在线视频| 国产精品av久久久久免费| 大码成人一级视频| 50天的宝宝边吃奶边哭怎么回事| 亚洲三区欧美一区| 国产男女超爽视频在线观看| 欧美一级毛片孕妇| 在线亚洲精品国产二区图片欧美| 国产成人精品在线电影| 久久久久久久久免费视频了| 国产亚洲av片在线观看秒播厂| 中文欧美无线码| 亚洲中文日韩欧美视频| 丁香六月天网| 叶爱在线成人免费视频播放| 午夜福利视频精品| 亚洲欧美色中文字幕在线| 欧美成狂野欧美在线观看| 中文字幕av电影在线播放| 韩国精品一区二区三区| 99热国产这里只有精品6| 亚洲专区国产一区二区| 天天影视国产精品| 国产成人一区二区三区免费视频网站| a级毛片黄视频| 精品久久久久久久毛片微露脸 | 欧美精品人与动牲交sv欧美| 男男h啪啪无遮挡| 亚洲专区中文字幕在线| a在线观看视频网站| 热re99久久国产66热| 一二三四社区在线视频社区8| 夜夜夜夜夜久久久久| 蜜桃国产av成人99| 美女高潮到喷水免费观看| 久久热在线av| 在线观看人妻少妇| 日本撒尿小便嘘嘘汇集6| 久久午夜综合久久蜜桃| 午夜精品国产一区二区电影| 一个人免费在线观看的高清视频 | 欧美在线一区亚洲| 男人操女人黄网站| 日韩一区二区三区影片| 久久亚洲国产成人精品v| 99精国产麻豆久久婷婷| 制服诱惑二区| 王馨瑶露胸无遮挡在线观看| 黑丝袜美女国产一区| www.精华液| 伊人久久大香线蕉亚洲五| 制服诱惑二区| 黄色视频在线播放观看不卡| 久久中文字幕一级| 成人三级做爰电影| 午夜免费观看性视频| 亚洲欧洲精品一区二区精品久久久| 亚洲精品一卡2卡三卡4卡5卡 | 国产xxxxx性猛交| 日本一区二区免费在线视频| 老司机亚洲免费影院| 日韩三级视频一区二区三区| 欧美人与性动交α欧美精品济南到| 亚洲精品中文字幕在线视频| 日韩大码丰满熟妇| 国产精品99久久99久久久不卡| 亚洲视频免费观看视频| 久久久久国产一级毛片高清牌| 黄网站色视频无遮挡免费观看| 国产福利在线免费观看视频| 国产99久久九九免费精品| 黄片小视频在线播放| 亚洲国产欧美在线一区| 精品人妻一区二区三区麻豆| 女性生殖器流出的白浆| 一本色道久久久久久精品综合| 国产精品99久久99久久久不卡| 国产成人精品在线电影| xxxhd国产人妻xxx| 国产成人系列免费观看| 国产成人欧美在线观看 | 国产一区二区三区综合在线观看| 午夜精品久久久久久毛片777| 国产无遮挡羞羞视频在线观看| 亚洲中文av在线| 在线av久久热| av电影中文网址| 日本撒尿小便嘘嘘汇集6| 一级毛片电影观看| 激情视频va一区二区三区| 亚洲精品国产精品久久久不卡| 男女之事视频高清在线观看| 90打野战视频偷拍视频| 99久久国产精品久久久| 啦啦啦在线免费观看视频4| 亚洲熟女毛片儿| 在线观看人妻少妇| cao死你这个sao货| 制服人妻中文乱码| 欧美黄色淫秽网站| 免费日韩欧美在线观看| 日本黄色日本黄色录像| 国产一区二区在线观看av| 亚洲国产欧美日韩在线播放| 欧美精品一区二区免费开放| 丰满少妇做爰视频| 国产欧美日韩一区二区三 | 99国产精品免费福利视频| 精品久久久久久电影网| 老司机靠b影院| 最近中文字幕2019免费版| www.熟女人妻精品国产| 丝袜美足系列| 搡老熟女国产l中国老女人| 亚洲 国产 在线| 少妇粗大呻吟视频| 99re6热这里在线精品视频| 亚洲色图 男人天堂 中文字幕| 亚洲第一欧美日韩一区二区三区 | 男女国产视频网站| 18禁黄网站禁片午夜丰满| 一本久久精品| 国产免费av片在线观看野外av| 色94色欧美一区二区| 黑人猛操日本美女一级片| 日本一区二区免费在线视频| 日韩一区二区三区影片| 亚洲五月婷婷丁香| 一二三四在线观看免费中文在| www.av在线官网国产| 淫妇啪啪啪对白视频 | 18禁裸乳无遮挡动漫免费视频| 亚洲五月色婷婷综合| 欧美激情 高清一区二区三区| 国产男人的电影天堂91|