• <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ū)
    18禁黄网站禁片免费观看直播| 欧美绝顶高潮抽搐喷水| 亚洲精品456在线播放app | 精品人妻偷拍中文字幕| 国产老妇女一区| 三级男女做爰猛烈吃奶摸视频| 久久精品亚洲精品国产色婷小说| 怎么达到女性高潮| 男女那种视频在线观看| 国产精品 欧美亚洲| 人妻久久中文字幕网| 欧美国产日韩亚洲一区| 国产色婷婷99| 免费一级毛片在线播放高清视频| 99国产精品一区二区三区| 亚洲国产精品合色在线| 两性午夜刺激爽爽歪歪视频在线观看| 国产免费男女视频| 欧美激情久久久久久爽电影| 亚洲av五月六月丁香网| 日本成人三级电影网站| 黄色成人免费大全| 一进一出抽搐gif免费好疼| 神马国产精品三级电影在线观看| 婷婷精品国产亚洲av在线| 国产老妇女一区| 国产精品亚洲一级av第二区| 欧美日韩黄片免| 成人永久免费在线观看视频| 婷婷亚洲欧美| 级片在线观看| 啦啦啦观看免费观看视频高清| x7x7x7水蜜桃| 国产成人欧美在线观看| 中国美女看黄片| 18美女黄网站色大片免费观看| 琪琪午夜伦伦电影理论片6080| 久久久久久久午夜电影| 国产精品一区二区三区四区久久| 亚洲精品国产精品久久久不卡| 偷拍熟女少妇极品色| 岛国在线免费视频观看| 亚洲精品久久国产高清桃花| 久久久久九九精品影院| 久久九九热精品免费| 欧美日韩一级在线毛片| 99精品久久久久人妻精品| 嫩草影院入口| 欧美激情在线99| 成年人黄色毛片网站| 好男人在线观看高清免费视频| 国产精品亚洲av一区麻豆| 午夜激情欧美在线| 国产亚洲欧美98| 精品日产1卡2卡| 国产高清videossex| 99热精品在线国产| 99久国产av精品| 久久久久精品国产欧美久久久| 国产色爽女视频免费观看| 中文字幕av在线有码专区| 亚洲一区二区三区不卡视频| 熟妇人妻久久中文字幕3abv| 国产精品久久久久久久电影 | 18禁在线播放成人免费| 一个人观看的视频www高清免费观看| 国产成年人精品一区二区| 亚洲18禁久久av| 国产91精品成人一区二区三区| 在线免费观看的www视频| av欧美777| 国产av不卡久久| 亚洲激情在线av| АⅤ资源中文在线天堂| 搡女人真爽免费视频火全软件 | 国产精品久久久人人做人人爽| 中文字幕人妻丝袜一区二区| 久久6这里有精品| 欧美日韩一级在线毛片| 精品久久久久久久久久久久久| 欧美高清成人免费视频www| 亚洲成人免费电影在线观看| 欧美日本亚洲视频在线播放| 欧美黄色片欧美黄色片| 亚洲aⅴ乱码一区二区在线播放| 无人区码免费观看不卡| 国产在视频线在精品| 成人永久免费在线观看视频| 精品一区二区三区视频在线观看免费| 午夜福利欧美成人| 亚洲国产欧美人成| 性色avwww在线观看| 色视频www国产| 免费在线观看日本一区| 变态另类成人亚洲欧美熟女| 色在线成人网| 丰满乱子伦码专区| 国产欧美日韩一区二区精品| 成人亚洲精品av一区二区| 看片在线看免费视频| 国产亚洲精品久久久久久毛片| 91av网一区二区| 欧美性猛交黑人性爽| 精品欧美国产一区二区三| 亚洲男人的天堂狠狠| 中文资源天堂在线| 亚洲成人久久性| 国产蜜桃级精品一区二区三区| 性色av乱码一区二区三区2| 毛片女人毛片| 欧美极品一区二区三区四区| 好看av亚洲va欧美ⅴa在| 国内精品美女久久久久久| 亚洲精品成人久久久久久| 国产黄a三级三级三级人| 最近最新免费中文字幕在线| 国产精品久久久久久久久免 | 尤物成人国产欧美一区二区三区| 一级毛片高清免费大全| 欧美区成人在线视频| 亚洲国产精品999在线| a级毛片a级免费在线| 亚洲欧美日韩高清专用| 精品久久久久久成人av| 国产成+人综合+亚洲专区| 中文字幕av成人在线电影| 色吧在线观看| av视频在线观看入口| 久久精品国产99精品国产亚洲性色| 可以在线观看毛片的网站| 女同久久另类99精品国产91| 日韩 欧美 亚洲 中文字幕| 听说在线观看完整版免费高清| 亚洲最大成人手机在线| svipshipincom国产片| 少妇熟女aⅴ在线视频| 午夜免费成人在线视频| 老汉色av国产亚洲站长工具| 中亚洲国语对白在线视频| 性欧美人与动物交配| 成人精品一区二区免费| 国产精品综合久久久久久久免费| 一进一出好大好爽视频| 亚洲欧美日韩卡通动漫| 性欧美人与动物交配| 91在线观看av| 国产高清视频在线播放一区| 久久久久久久久久黄片| 天天躁日日操中文字幕| 少妇的丰满在线观看| 欧美日本视频| 亚洲专区中文字幕在线| 日本精品一区二区三区蜜桃| 免费看十八禁软件| 校园春色视频在线观看| 欧美精品啪啪一区二区三区| 亚洲精品在线观看二区| 长腿黑丝高跟| 亚洲最大成人手机在线| 中文在线观看免费www的网站| 婷婷亚洲欧美| 国产真人三级小视频在线观看| 久久久精品欧美日韩精品| 美女 人体艺术 gogo| 国产激情欧美一区二区| 1000部很黄的大片| 国产高清视频在线观看网站| 桃红色精品国产亚洲av| 制服人妻中文乱码| 国产精品99久久99久久久不卡| 五月伊人婷婷丁香| 麻豆国产97在线/欧美| 一级作爱视频免费观看| 国内精品久久久久精免费| 成人特级黄色片久久久久久久| 99久久精品热视频| 老鸭窝网址在线观看| 欧美极品一区二区三区四区| 国产毛片a区久久久久| 日本与韩国留学比较| 偷拍熟女少妇极品色| 首页视频小说图片口味搜索| 欧美日韩国产亚洲二区| 免费看a级黄色片| 99国产综合亚洲精品| 亚洲欧美激情综合另类| 亚洲中文日韩欧美视频| 又粗又爽又猛毛片免费看| 免费高清视频大片| 一本一本综合久久| 欧美一级a爱片免费观看看| 动漫黄色视频在线观看| 亚洲av美国av| 成年人黄色毛片网站| 国产色婷婷99| 日本免费一区二区三区高清不卡| 久久香蕉国产精品| 亚洲一区二区三区色噜噜| 18禁国产床啪视频网站| 禁无遮挡网站| 免费人成在线观看视频色| 1000部很黄的大片| 老汉色av国产亚洲站长工具| 国产淫片久久久久久久久 | 亚洲国产精品久久男人天堂| 亚洲精品一卡2卡三卡4卡5卡| 国产免费av片在线观看野外av| 母亲3免费完整高清在线观看| 蜜桃久久精品国产亚洲av| 国语自产精品视频在线第100页| 国产欧美日韩一区二区三| 久久久久国内视频| 国产色爽女视频免费观看| 色哟哟哟哟哟哟| 天堂√8在线中文| 成人av在线播放网站| 看片在线看免费视频| 人人妻人人澡欧美一区二区| 亚洲av一区综合| 又爽又黄无遮挡网站| 熟妇人妻久久中文字幕3abv| 久久精品国产自在天天线| 亚洲精华国产精华精| eeuss影院久久| 少妇裸体淫交视频免费看高清| 狂野欧美激情性xxxx| 国产激情偷乱视频一区二区| 欧美一区二区精品小视频在线| 国产一区二区亚洲精品在线观看| 国产成人av教育| 久久精品国产综合久久久| 村上凉子中文字幕在线| 91在线精品国自产拍蜜月 | 他把我摸到了高潮在线观看| www.www免费av| 久久国产精品人妻蜜桃| 久久久久久人人人人人| 午夜福利免费观看在线| 97超级碰碰碰精品色视频在线观看| 亚洲国产日韩欧美精品在线观看 | 亚洲一区高清亚洲精品| 在线十欧美十亚洲十日本专区| 在线观看一区二区三区| 天天躁日日操中文字幕| 日本一本二区三区精品| 两个人的视频大全免费| 90打野战视频偷拍视频| 久久久国产成人精品二区| 婷婷精品国产亚洲av| 欧美午夜高清在线| 国产欧美日韩一区二区精品| 午夜日韩欧美国产| 欧美在线黄色| 精品久久久久久久末码| 久久久久久久亚洲中文字幕 | 国产高清激情床上av| 亚洲五月婷婷丁香| 亚洲欧美日韩高清在线视频| 亚洲专区中文字幕在线| 亚洲av一区综合| 日韩欧美在线二视频| 熟女电影av网| 亚洲内射少妇av| 欧美色欧美亚洲另类二区| 欧美激情久久久久久爽电影| 丰满人妻一区二区三区视频av | 国产精品99久久久久久久久| avwww免费| 午夜免费激情av| av专区在线播放| 无遮挡黄片免费观看| 最近最新免费中文字幕在线| 国内精品久久久久精免费| 亚洲精品在线观看二区| 欧美在线一区亚洲| 淫妇啪啪啪对白视频| 久久精品人妻少妇| 国产黄色小视频在线观看| 国产在视频线在精品| 国内少妇人妻偷人精品xxx网站| 精品久久久久久久毛片微露脸| 狂野欧美激情性xxxx| 国产老妇女一区| 国产精品日韩av在线免费观看| 国产v大片淫在线免费观看| 女同久久另类99精品国产91| 国产一区二区亚洲精品在线观看| 成年女人毛片免费观看观看9| av天堂中文字幕网| 中文字幕人妻丝袜一区二区| 两人在一起打扑克的视频| 波野结衣二区三区在线 | 欧美色欧美亚洲另类二区| 国产精品99久久99久久久不卡| 国产欧美日韩精品亚洲av| 欧美丝袜亚洲另类 | 婷婷丁香在线五月| 成人特级黄色片久久久久久久| 舔av片在线| 日本三级黄在线观看| 亚洲一区二区三区色噜噜| 亚洲av一区综合| 中文在线观看免费www的网站| 亚洲av成人不卡在线观看播放网| 狠狠狠狠99中文字幕| 国产成人a区在线观看| 日本与韩国留学比较| 久久天躁狠狠躁夜夜2o2o| 国产精品野战在线观看| avwww免费| 熟女人妻精品中文字幕| 757午夜福利合集在线观看| 99精品欧美一区二区三区四区| 免费在线观看亚洲国产| 亚洲专区国产一区二区| 久久香蕉国产精品| 熟女少妇亚洲综合色aaa.| 日本成人三级电影网站| 国产精品自产拍在线观看55亚洲| 一卡2卡三卡四卡精品乱码亚洲| 国内精品美女久久久久久| 欧美成人性av电影在线观看| 最近最新中文字幕大全电影3| 神马国产精品三级电影在线观看| 亚洲精华国产精华精| 波多野结衣巨乳人妻| 久久99热这里只有精品18| 亚洲成人久久爱视频| 精品福利观看| 国产亚洲欧美在线一区二区| 内地一区二区视频在线| 麻豆久久精品国产亚洲av| 丰满乱子伦码专区| 老司机午夜福利在线观看视频| 精品国产超薄肉色丝袜足j| 亚洲av一区综合| 久久久久久久亚洲中文字幕 | 久久久久久久久久黄片| 国产97色在线日韩免费| 国产老妇女一区| 听说在线观看完整版免费高清| 黄色丝袜av网址大全| 此物有八面人人有两片| 18美女黄网站色大片免费观看| 亚洲国产欧美网| 国产精品久久久久久精品电影| 国产成人av教育| 国产精品亚洲一级av第二区| 欧美一区二区亚洲| 国产色爽女视频免费观看| 色在线成人网| 好男人在线观看高清免费视频| 又黄又爽又免费观看的视频| 免费高清视频大片| 精品99又大又爽又粗少妇毛片 | 丁香六月欧美| 99久久精品热视频| 99久久精品国产亚洲精品| 在线观看免费视频日本深夜| 超碰av人人做人人爽久久 | 亚洲精品在线观看二区| 俺也久久电影网| 搞女人的毛片| 波多野结衣高清作品| 精品久久久久久成人av| 国产精品久久久久久亚洲av鲁大| www.www免费av| 两性午夜刺激爽爽歪歪视频在线观看| 99在线人妻在线中文字幕| 日本三级黄在线观看| 日本与韩国留学比较| 男人和女人高潮做爰伦理| 色老头精品视频在线观看| 国产一区二区激情短视频| 久久国产精品人妻蜜桃| 18禁黄网站禁片免费观看直播| 国产真实乱freesex| 在线观看66精品国产| av在线天堂中文字幕| 高潮久久久久久久久久久不卡| 美女免费视频网站| 狂野欧美白嫩少妇大欣赏| 波野结衣二区三区在线 | 黄色丝袜av网址大全| 精品久久久久久,| bbb黄色大片| 日本黄色视频三级网站网址| 最近最新中文字幕大全免费视频| 丁香六月欧美| 宅男免费午夜| 老司机福利观看| 99国产极品粉嫩在线观看| 国产亚洲av嫩草精品影院| 亚洲国产精品合色在线| 噜噜噜噜噜久久久久久91| 欧美一级a爱片免费观看看| 国产精品av视频在线免费观看| aaaaa片日本免费| 夜夜看夜夜爽夜夜摸| 亚洲av第一区精品v没综合| 可以在线观看毛片的网站| 精品国产美女av久久久久小说| 成年女人毛片免费观看观看9| 男人舔奶头视频| 欧美色视频一区免费| 国产精品一区二区三区四区免费观看 | 国产免费av片在线观看野外av| 国产精品亚洲av一区麻豆| 又粗又爽又猛毛片免费看| 男人和女人高潮做爰伦理| 亚洲专区国产一区二区| 亚洲七黄色美女视频| 18+在线观看网站| 变态另类丝袜制服| 国产av一区在线观看免费| 3wmmmm亚洲av在线观看| 亚洲成a人片在线一区二区| 久久久久久九九精品二区国产| 久久久久九九精品影院| 成人高潮视频无遮挡免费网站| 嫩草影院精品99| 99精品在免费线老司机午夜| 精品乱码久久久久久99久播| 精品久久久久久成人av| 久久国产乱子伦精品免费另类| bbb黄色大片| 悠悠久久av| 成年人黄色毛片网站| 国产精品99久久99久久久不卡| 性色av乱码一区二区三区2| 欧美高清成人免费视频www| 床上黄色一级片| 99热这里只有是精品50| 一进一出好大好爽视频| 人妻久久中文字幕网| 一个人看视频在线观看www免费 | 国产亚洲精品一区二区www| 欧美黄色片欧美黄色片| 色尼玛亚洲综合影院| www日本黄色视频网| 色视频www国产| 婷婷精品国产亚洲av在线| 国产真实乱freesex| 国产一区二区三区视频了| 男人的好看免费观看在线视频| 亚洲av二区三区四区| 日韩av在线大香蕉| 特级一级黄色大片| 乱人视频在线观看| 国产激情偷乱视频一区二区| 亚洲国产欧洲综合997久久,| 真实男女啪啪啪动态图| 香蕉av资源在线| 国产av一区在线观看免费| 国产高清激情床上av| 精品99又大又爽又粗少妇毛片 | 亚洲欧美日韩东京热| 国产精品精品国产色婷婷| a级一级毛片免费在线观看| 亚洲18禁久久av| 两个人看的免费小视频| 精品国产美女av久久久久小说| 国产久久久一区二区三区| 亚洲片人在线观看| 99热这里只有精品一区| 亚洲片人在线观看| 亚洲不卡免费看| 男女做爰动态图高潮gif福利片| 亚洲男人的天堂狠狠| 少妇丰满av| 熟女少妇亚洲综合色aaa.| 国产一区在线观看成人免费| 美女高潮的动态| 亚洲乱码一区二区免费版| 亚洲欧美激情综合另类| 午夜福利在线在线| 在线看三级毛片| 免费看a级黄色片| 成年免费大片在线观看| av中文乱码字幕在线| 一级黄片播放器| 亚洲精品在线观看二区| 亚洲欧美激情综合另类| 免费看a级黄色片| 欧美3d第一页| 国产三级中文精品| 精品一区二区三区视频在线 | 国产欧美日韩一区二区精品| 亚洲成人久久性| 成人性生交大片免费视频hd| 色哟哟哟哟哟哟| 成人av在线播放网站| 亚洲专区中文字幕在线| 18禁黄网站禁片午夜丰满| 日韩精品青青久久久久久| 亚洲人与动物交配视频| 久久欧美精品欧美久久欧美| 免费在线观看亚洲国产| 国产av在哪里看| 最新在线观看一区二区三区| 在线视频色国产色| 精品日产1卡2卡| 国产一区二区在线观看日韩 | 1000部很黄的大片| 99精品久久久久人妻精品| 日韩人妻高清精品专区| 老司机午夜十八禁免费视频| 午夜免费成人在线视频| 性欧美人与动物交配| 欧美在线黄色| 国产精品 欧美亚洲| 久久精品国产清高在天天线| 级片在线观看| 亚洲av美国av| 免费人成在线观看视频色| 九九热线精品视视频播放| 91在线观看av| 亚洲av第一区精品v没综合| 亚洲色图av天堂| 午夜老司机福利剧场| 亚洲 欧美 日韩 在线 免费| 麻豆成人午夜福利视频| 此物有八面人人有两片| 国产亚洲欧美98| 欧美在线黄色| www.熟女人妻精品国产| 色综合亚洲欧美另类图片| 最新美女视频免费是黄的| 精品乱码久久久久久99久播| 中文字幕久久专区| 90打野战视频偷拍视频| 国产亚洲精品一区二区www| 日日干狠狠操夜夜爽| 波野结衣二区三区在线 | 免费av观看视频| 亚洲av电影在线进入| 国产免费男女视频| 久久亚洲真实| 国产亚洲精品久久久久久毛片| 欧美色视频一区免费| 国产三级中文精品| 欧美中文日本在线观看视频| 精华霜和精华液先用哪个| 久久伊人香网站| 在线免费观看的www视频| 国产av麻豆久久久久久久| av在线天堂中文字幕| 又黄又爽又免费观看的视频| 最后的刺客免费高清国语| 亚洲电影在线观看av| 亚洲人与动物交配视频| 欧美av亚洲av综合av国产av| 国产久久久一区二区三区| 日韩欧美三级三区| 麻豆一二三区av精品| 欧美性感艳星| 岛国在线观看网站| 俺也久久电影网| 精品人妻偷拍中文字幕| 桃色一区二区三区在线观看| 亚洲美女黄片视频| 亚洲avbb在线观看| 美女黄网站色视频| 天堂av国产一区二区熟女人妻| av在线天堂中文字幕| 亚洲av五月六月丁香网| 久久性视频一级片| 久久精品人妻少妇| 91麻豆av在线| 国产精品一及| 免费人成在线观看视频色| 18禁在线播放成人免费| 国语自产精品视频在线第100页| 毛片女人毛片| 婷婷精品国产亚洲av在线| 亚洲人成网站在线播| 天堂动漫精品| 香蕉久久夜色| 午夜福利在线观看吧| 久久久久久人人人人人| 首页视频小说图片口味搜索| 国产蜜桃级精品一区二区三区| 脱女人内裤的视频| 俄罗斯特黄特色一大片| 黄色成人免费大全| 亚洲欧美日韩高清专用| 午夜影院日韩av| 尤物成人国产欧美一区二区三区| 床上黄色一级片| netflix在线观看网站| 老鸭窝网址在线观看| 亚洲中文字幕一区二区三区有码在线看| 久久久久性生活片| 婷婷六月久久综合丁香| 美女被艹到高潮喷水动态| 亚洲av第一区精品v没综合| 成人特级黄色片久久久久久久| x7x7x7水蜜桃| 黄色女人牲交| 真人做人爱边吃奶动态| 男女之事视频高清在线观看| 一区二区三区激情视频| 亚洲人成伊人成综合网2020| 国产主播在线观看一区二区| 99热精品在线国产| 亚洲欧美一区二区三区黑人| 午夜福利在线观看免费完整高清在 | 欧美黑人欧美精品刺激| 在线观看舔阴道视频| 久久精品91无色码中文字幕| ponron亚洲| 精品午夜福利视频在线观看一区| 亚洲18禁久久av| 伊人久久精品亚洲午夜| 成年人黄色毛片网站|