999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

A simulation study of a windowless gas-stripping room in an E//B neutral particle analyzer

2021-08-05 08:22:02YuanLuoWeiPingLinIDPeiPeiRenGuoFengQuJingJunZhuXingQuanLiuIDXiaoBingLuoZhuAnIDRoyWadaLinGeZangYuFanQuZhongBingShi
Nuclear Science and Techniques 2021年7期

Yuan Luo Wei-Ping Lin ID Pei-Pei Ren Guo-Feng Qu Jing-Jun Zhu Xing-Quan Liu ID Xiao-Bing Luo Zhu An ID Roy Wada Lin-Ge Zang Yu-Fan Qu Zhong-Bing Shi

Abstract The neutral particle analyzer(NPA)is one of the crucial diagnostic devices in a Tokamak facility.The stripping unit is one of the main parts of the NPA.A windowless gas-stripping room with two differential pipes has been constructed in a parallel electric and magnetic fields(E//B)NPA.The pressure distributions in the stripping chamber are simulated by ANSYS Fluent together with MolFlow+.Based on the pressure distributions obtained from the simulation,the stripping efficiency of the E//B NPA is studied using GEANT4.Hadron reaction physics is modified to track the charge state of each particle in a cross-section-based method in GEANT4. The transmission rates (R) and stripping efficiencies f+1 are examined for particle energies ranging from 20 to 200 keV with the input pressure (P0),ranging from 20 to 400 Pa.According to the combined global efficiency,R×f+1,P0=240 Pa is obtained as the optimum pressure for the maximum global efficiency in the incident energy range investigated.

Keywords Neutral particle analyzer · Windowless gasstripping chamber · Stripping efficiency · Ansys Fluent ·MolFlow+ · GEANT4

1 Introduction

Tokamak is a toroidal device used in nuclear fusion research for magnetic confinement plasmas. It provides a place to test integrated technologies,materials,and physics regimes necessary for future commercial production of fusion-based electricity [1]. The neutral particle analyzer(NPA)is a crucial diagnostic device in the Tokamak facilities.It is used to determine the bulk ion temperature,isotopic ratio,and fast ion distribution of the plasma by measuring the charge exchange neutral particles escaping from the plasma.Different types of NPAs have been built in Tokamak facilities worldwide [2–13], such as the parallel direction of electric and magnetic fields (E//B) NPA on the Tokamak Fusion Test Reactor (TFTR) [2], compact neutral particle analyzer (CNPA) on the Wendelstein 7-AS stellarator [3],low- and high-energy NPAs (LENPA and HENPA) on the International Thermonuclear Experimental Reactor(ITER)[10], solid-state NPA (ssNPA) on the Experimental Advanced Superconducting Tokamak(EAST)[11],and the CP-NPA on the HuanLiuqi-2A(HL-2A)[12].

The stripping unit of the NPA plays an important role in the analysis of neutral particles, except in flux measurement NPAs, such as ssNPA [11]. It provides a place to reionize the charge-exchange neutral particles. According to the state of the stripping material,the stripping unit can be divided into two types: stripping foil and gas chamber.When a stripping foil is used in the NPA for low-energy neutrals, an additional accelerating or focusing voltage is required for the secondary ions [3, 9, 10]. A carbon foil with a thickness of 100 A?is commonly used as the stripping foil. In contrast, a gas chamber requires a differential pumping system when stripping gas is used. Typically, an integrated target thickness of the order of 1016atoms/cm2for H2gas is used in the Joint European Torus (JET) NPA[13], and 1015atoms/cm2for He gas is used in the E//B NPA on TFTR [2].

The energetic particles,also known as fast,superthermal,hot,and high-energy particles,are expected to play a critical role in plasma heating,current drive,momentum transport,energy transfer,and plasma stability[14,15].Many experimental and theoretical studies have been conducted in this field[16–20],and other related fields[21–26].To study the frontier physics of energetic particles and to measure the fuel ratio, a new E//B NPA was designed in the present experimental devices[27].This E//B NPA is a tandem-type NPA—similar to the CNPA built at the Ioffe Physicotechnical Institute, Russia [3]. It provides mass resolution (H and D resolution)for particles in the energy range of 20–200 keV.The magnetic field is generated with a permanent magnet for a smaller size of the NPA and easy maintenance.The upper limit energy of the E//B NPA was determined from the negative ion source neutral beam heating on the Huanliuqi-2M (HL-2M) device. The lower limit was set to 20 keV because our study focused on the fast ions instead of the background ions.For more details,refer to our previous work in Ref.[27].

In this study,a gas-stripping chamber of a new E//B NPA was designed and studied. A windowless gas-stripping chamber was adopted to avoid the replacement of the stripping foils and ensure their easy maintenance during actual operation. The performance of the gas-stripping chamber was investigated using ANSYS Fluent [28, 29] and Mol-Flow+[30],together with GEANT4[31,32].This article is organized as follows:The design and pressure calculation of the gas-stripping chamber is presented in Sect. 2.The results of the GEANT4 simulation and discussions are presented in Sect. 3.A brief summary is provided in Sect. 4.

2 Design and pressure calculation of the gasstripping chamber

The stripping unit is one of the main parts of the NPA.The electrons of the escaped neutral particles are stripped in the stripping unit. A windowless gas-stripping chamber was adopted in the design to avoid replacement of the stripping foil and for easy maintenance. To obtain a specific high pressure inside the stripping room and high vacuum in the outside vacuum chamber at the same time,two differential pipes with small flow conductance were used for the stripping room. Figure 1 shows the schematic layout of the stripping chamber.The stripping room(1)with two differential pipes(2)36 mm in length and 4 mm in diameter was placed inside a vacuum chamber(3).Two holes with diameters of 6 mm(7 and 8)were made on the entrance and exit flanges to limit the beam size and maximize the vacuum isolation from the upstream pipe and downstream chamber.H2gas was used as the stripping gas to avoid polluting the Tokamak fuel.It was filled in the stripping room from the top flange (5) of the vacuum chamber through a bellow (4). A machinery-bearing molecular pump with a pumping speed of 340 L/s was used at the bottom of the vacuum chamber together with a gate valve(11).

The pressure distribution inside the stripping chamber is one of the main concerns in our design. A pressure of dozens of Pa is required in the stripping room to obtain sufficient stripping efficiency for high-energy hydrogen(H)and deuteron(D)atoms.In this pressure region,the gas flow state in the stripping room remains in a viscousmolecular flow state[33].The pressure in the gas inlet and bellow was higher than that in the stripping room. On the contrary, a pressure, two or three orders of magnitude lower, was estimated in the vacuum chamber. The mean free path of the gas molecule inside the vacuum chamber is larger than the size of the chamber, and the motion of the gas molecule can be treated as collision-less. The Monte Carlo code of MolFlow+ is often used to calculate the pressure distribution of collision-less gas in a high-vacuum system.However,it is not accurate for all gas regions in the gas-stripping chamber. The computational fluid dynamics(CFD)software,which can evaluate the nonlinear effect of viscous fluid,has a better performance in the high-pressure region. However, the calculations of the CFD software in the high-vacuum region show an unnatural bump at the corners. Therefore, the gas pressure distribution of the stripping chamber was calculated by combining the CFD software of ANSYS Fluent[28,29](for the viscous region in the bellow, stripping room, and differential pipes) and the Monte Carlo software of MolFlow+ [30] (for the lowpressure collision-less region in the vacuum chamber).

Three-dimensional CFD calculations were performed using the ANSYS Fluent software. The fluid region was established according to the structure of the stripping chamber,as shown in Fig. 1.A laminar viscous model was adopted in the calculations. A pressure-type gas inlet was present at the top flange (5) of the vacuum chamber, and the pressure at the gas inlet P0ranged from 20 to 400 Pa with a step of 20 Pa set in the simulation. Three pressuretype gas outlets were present at the entrance and exit holes(7 and 8), and at the bottom of the gate valve (11). A pressure of 10-3Pa was assumed for all three outlets.Owing to the large flow conductance of the two differential pipes (2), a small change in the outlet pressure does not affect the pressure distribution in the stripping room.Moreover,the pressure distribution in the vacuum chamber was changed based on the results from MolFlow+.Therefore,a pressure of 10-3Pa at the outlets was used for all the ANSYS Fluent calculations.Stainless steel was used as the wall material. A room temperature of 300 K was used for the calculations. A typical gas flow rate of 9.97 Pa L/s was obtained at the bottom of the bellow for the input pressure P0= 100 Pa.

In the low-pressure region in the vacuum chamber(3)in Fig. 1 and inside the differential pipe (2) at |z|>42 mm,where z=0 was set at the center of the stripping room,the pressure distribution was simulated by MolFlow+ at the same temperature of 300 K.For the MolFlow+simulation,the outgassing rate was adopted from the ANSYS Fluent calculations at z=± 42 mm inside the differential pipe, 4 mm from the pipe exit. For pumping, it was assumed that the gas molecules are absorbed when they hit the surfaces of the entrance and exit holes(7 and 8);that is,the sticking factor was set to 1 on the surfaces. This resulted in a pumping speed of 12.4 L/s through the entrance and exit holes.A pumping speed of 340 L/s was set at the bottom of the gate valve (11).

The simulated two-dimensional (2D) pressure distribution at P0=100 Pa in the y-z plane at x=0 is shown in Fig. 2a. The detailed pressure distribution around the stripping room is shown in Fig. 2b at a magnified scale.The pressure distribution along the beam line is shown in Fig. 2c plotted on a logarithmic scale. Using the two differential pipes in the design, the desired high pressure was achieved inside the stripping room,and a linear decrease in pressure was observed inside the two differential pipes and the entrance and exit holes. A sharp change in pressure at the entrance of the two differential pipes was observed in the ANSYS Fluent calculations.However,it was not found when MolFlow+was used to simulate the entire gas region in the stripping chamber.

To evaluate the pressure distribution change as the pressure at the gas inlet changes, the pressures at four typical positions were examined for all the P0values investigated. Figure 3 shows the pressure inside the stripping room (P1), at the entrance of the differential pipes(P2), at the vacuum chamber (P3), and at the outside surfaces of the entrance and exit holes (P4), as shown in Fig. 2c as a function of pressure at the gas inlet P0.Linear changes in P1, P2, P3, and P4on P0were obtained,although slight fluctuations were observed. By using the two differential pipes in the design, the pressure in the vacuum chamber was observed to be more than 500 times lower than that in the stripping room.Through the entrance and exit holes, pressure, approximately one order of magnitude lower, was obtained for the upstream pipe and downstream chamber.

Fig. 3 (Color online) Pressure P1 (solid circles), P2 (solid squares),P3 (solid triangles), and P4 (solid triangles) as a function of P0. The lines were obtained using linear fits

The integrated target thickness (nT) is an important quantity in gas-stripping chambers.It is commonly used to evaluate the efficiency of the stripping chamber. Because the pressures at the four typical positions were used to construct the pressure distribution in GEANT4 in the next section,a comparison of the exact nTand the nTcalculated from the pressure distribution used in GEANT4 is necessary. Figure 4 shows nTas a function of P0for the results of ANSYS Fluent and MolFlow+(solid circles)and that of GEANT4 (open circles). A good agreement was found between them.

3 Results of GEANT4 simulation and discussion

Utilizing the obtained P1, P2, P3, and P4, the stripping efficiencies of the stripping room were further studied using GEANT4 [31, 32]. GEANT4 is a toolkit for simulating the passage of particles through matter. It has been applied in various studies, such as our previous studies on the average neutron detection efficiency for the DEtecteur MOdulaire de Neutrons detectors (DEMON) [34] and the module test of the Collision Centrality Detector Array(CCDA) [35], as well as many other studies of electron backscattering [36], neutron time-of-flight spectrometer system at HL-2M[37],and performance of a large-size CsI detector [38].

Fig. 4 (Color online) nT as a function of P0. Solid circles are the exact nT values calculated from the pressure distributions obtained from the ANSYS Fluent and MolFlow+simulations.Open circles are the nT values calculated from the pressure distributions constructed in the GEANT4 simulation

In the GEANT4 simulation, the physics list includes electromagnetic physics [39] and hadronic physics [40],in which ion transportation and electromagnetic, nuclear elastic, and inelastic processes are activated—although some of the processes may not be used in the actual simulation. This is because the integrated target thickness is approximately 1016atoms/cm2, and the scattering probability of the incident particles and the target atoms is small enough such that multi-scattering is negligible. Therefore,the G4ScreenedNuclearRecoil class [41, 42] is included in the standard electromagnetic physics for incident energies ranging from 10 eV to 100 MeV.

The charge states of the H and D atoms are the key variables in this study.However,the original GEANT4 could not properly handle the charge state.To simulate the charge state variation,the hadron reaction physics was modified to trace the charge state of each particle in a cross-section-based method.By introducing a global charge state variable in the hadron reaction physics,the charge states of H and D were recorded when the charge exchange reaction occurred in the gas-stripping chamber. During the last century, many charge-exchange cross-section measurements have been performed for H on H2gas [43–54]. Figure 5 shows the electron loss cross-sections of H0(σ0,1), electron capture cross-sections of H+(σ1,0), electron capture cross-sections of H0(σ0,-1),and electron loss cross-sections of H-(σ-1,0)on H2gas as a function of the incident energy(E)in(a),(b),(c),and(d),respectively.Solid circles,solid squares,solid up triangles, solid down triangles, open circles, open squares,and open triangles represent the data from Gealy [45, 46],Stier[51],Barnett[43],Sanders[48],Smith[50],McClure[47],and Van Zyl[52],respectively.Solid curves represent the ORNL-recommended cross-section[53,54].Because the experimental σ-1,0does not cover the high-energy region at E greater than 30 keV,the ORNL recommended cross-section of the H atom on H2gas is used in the simulation.For a given velocity,the charge exchange cross-sections of D on Cs[55]or Rb[56]vapor are the same as those of H.Incident energies per nucleon(E/A)of H and D are the same for the investigated energy range when they have the same velocity.Therefore,the charge exchange cross-sections of H are also used for D at the same E/A in the simulations.

Fig. 5 (Color online) Charge exchange cross-section of σ0,1 (a), σ1,0(b),σ0,-1 (c),and σ-1,0 (d)in the order of of 10-16 cm2 per molecule as a function of the incident energy (E) of the H atom or ion. Solid circles, solid squares, solid up triangles, solid down triangles, open circles, open squares, and open triangles represent the data from Gealy [45, 46], Stier [51], Barnett [43], Sanders [48], Smith [50],McClure [47], and Van Zyl [52], respectively. The solid curves represent the ORNL-recommended cross-section [53, 54], which is also used in the simulation

The simulations were performed for H and D on H2gas with incident energies ranging from 20 to 200 keV in steps of 20 keV and with P0ranging from 20 to 400 Pa in steps of 20 Pa.H and D atoms are generated at the entrance hole(7)of the vacuum chamber,corresponding to the z position of -120 mm, and distributed uniformly on the entrance hole surface with the same diameter of 6 mm. The momentum direction was assumed to be parallel to the zaxis. One million events were generated for each run. The energy losses(ΔE/E)of H and D at 20 keV as a function of P0are shown in Fig. 6. A slightly lower energy loss is observed for D because the mass of D is twice that of H,which causes less energy loss during the collisions. It can be seen from Fig. 6 that H and D increase linearly as P0increases. The maximum energy loss of 20 keV for H and D was less than 4% for all P0investigated. Compared to the energy resolution of this NPA, the energy losses of H and D after passing through the stripping chamber are small and can be neglected.

Because the diameters of the two differential pipes are smaller than these of the entrance and exit holes, some of the incident particles will be stopped by the geometry of the stripping chamber. Moreover, owing to the Coulomb scattering of the incident particles and the target atoms,some of the incident particles are scattered away from their original directions. Therefore, the transmission rate (R) of the incident particles is important, especially for low-energy particles that suffer more Coulomb scattering when they pass through the stripping chamber. In this study, the transmission rate is defined as the ratio between the number of particles reached at the exit hole (8) after passing through the stripping chamber with H2gas (P0>0 Pa)and that without the gas (P0= 0 Pa, vacuum). Thus, the particle loss caused by the geometry of the stripping chamber was canceled out. Figure 7 shows the transmission rate as a function of the incident energy for P0=20 Pa(circles), 100 Pa (squares), and 400 Pa (triangles) in (a);and as a function of P0for the incident energy E =20 keV(circles), 40 keV (squares), and 100 keV (triangles) in (b).The solid and open symbols represent H and D,respectively.A slight increasing trend was observed for the transmission rate as the incident energy increased, but it showed an opposite trend as P0increased. The scattering loss was small (less than 3%) for all the investigated incident energies and input pressures.

The stripping efficiency in the stripping chamber was evaluated using the charge fraction variable (f) for the incident particles after the stripping area. The evolution of the charge fraction inside the stripping area is the most important aspect of this study. Figure 8 shows the charge state fraction as a function of the z position in the stripping chamber for H and D atoms at 20, 100, and 200 keV. The pressures used in the simulation are those from the results of P0=100 Pa,as shown in Sect. 2.The solid,dashed,and long-dashed curves correspond to the fractions of the charge states 0, +1, and -1, respectively. Owing to the small cross-section of σ0,-1,the fraction of charge state-1 is less than 2%for H and D at an incident energy of 20 keV and becomes negligible for higher incident energies. The fraction of charge states 0 and +1 shows a sharp change starting from around z -40 mm, which is because the stripping room is located at-46 mm <z<46 mm.Owing to the energy dependence of the charge exchange crosssection,particles with lower incident energies have smaller saturation thicknesses for the charge fractions. The higher stripping efficiency for particles with larger energies is mainly caused by the sharp decrease in σ1,0as the incident energy increases.

Fig.8 (Color online)Charge state fraction as a function of z position in the stripping chamber for H at 20 keV in(a),D at 20 keV in(b),H at 100 keV in(c),D at 100 keV in(d),H at 200 keV in(e),and D at 200 keV in(f).The pressures used in the simulation were those from the results of P0 =100 Pa.The solid,dashed,and long-dashed curves correspond to the fractions of charge states 0, +1, and -1,respectively

Figure 9 shows the fraction of charge state +1 (f+1) at z=120 mm as a function of E/A for P0= 20 Pa (circles),100 Pa (squares) and 400 Pa (triangles) in (a), and as a function of P0for the incident energy E=20 keV(circles),100 keV(squares)and 200 keV(triangles)in(b).Solid and open symbols represent H and D,respectively.One can see from Fig. 9a that f+1increases at lower incident energies for different P0up to E/A around 100 keV. After reaching the maximum value,f+1decreases for P0=20 Pa,and stays flat for P0= 100 Pa, but keeps slowly increasing for P0=400 Pa as E/A increases.This indicates that the thickness of the stripping gas is not enough for higher energy particles at P0=20 Pa.No noticeable difference between H and D is observed, indicating that the results can also be applied to neutral tritium particles when E/A is used. As shown in Fig. 9b, the fractions of charge state +1 quickly reach a maximum value and keep the maximum as P0increases at lower E.For larger E,the fractions increase faster at lower P0and reach the maximum at a pressure of around P0=240 Pa.

Fig.9 (Color online)a f+1 as a function of E/A.Circles,squares,and triangles represent the values at P0 = 20, 100, and 400 Pa,respectively. b f+1 as a function of P0.Circles,squares, and triangles represent the values at E = 20, 100, and 200 keV, respectively. The solid and open symbols represent H and D, respectively

To verify the GEANT4 results, f+1is also calculated using the gas-integrated target thickness (nT) as σ01and σ10are the stripping (electron loss) and charge exchange (electron capture) cross-sections of H(D) and H+(D+), respectively. The small amount of particle loss due to the electron capture of H(D), σ0-1, and σ-10, is neglected in Eq. (1).The results are shown in Fig. 10.The symbols are the same f+1values from the GEANT4 simulation in Fig. 9b but plotted as a function of nT, and the solid and dashed lines are those from Eq. (1) for H and D,respectively. Good agreement was found between the calculations and GEANT4 simulations. These good agreements originate from the following facts: In the GEANT4 simulation, the evaluation of f+1is obtained from a Monte Carlo sampling of the charge state along the particle track according to the cross-sections of σ01, σ10, σ0-1, and σ-10.As mentioned earlier, in Eq. (1), only a part of the crosssections (σ01and σ10) is used, neglecting σ0-1and σ-10,because the latter values are orders of magnitude smaller.

To determine the optimum condition, the global efficiency, a combination of transmission rate and fraction of charge state +1, R×f+1, was further studied. Figure 11 shows the global efficiency of R×f+1as a function of P0for the incident energies E = 20, 100, and 200 keV in (a),(b), and (c), and as a function of nTin (d), (e), and (f),respectively. The global efficiency of H and D decreases gradually as the pressure P0increases for E = 20 keV. On the other hands,the global efficiency shows a trend similar to that of f+1as the pressure P0increases for E = 100 and 200 keV.For pressures P0>240 Pa,the global efficiency becomes flat for all conditions of H and D at the incident energy E ≥100 keV. Considering the low temperature of the plasma in HL-2A/M, the number of high-energy particles is orders of magnitude lower than that of the lowenergy particles. Therefore, P0= 240 Pa was obtained as the optimum pressure for the maximum global efficiency in the incident energy range investigated. At this P0, the pressure in the vacuum chamber is less than 0.1 Pa, which is within the operating pressure range of the molecular pump.The simulation results would provide a useful guide for practical applications.

Fig. 10 (Color online) The fraction of H+ (solid symbols) and D+(open symbols)as a function of nT for E=20 keV(circles),100 keV(squares), and 200 keV (triangles). The solid and dashed lines represent the results of Eq. (1). More details are provided in the text

Fig.11 Global efficiency of R×f+1 as a function of the pressure P0 for incident energies E=20 keV(a),100 keV(b),and 200 keV(c–f)is the same as those in a–c—but as a function of nT. The solid and open circles represent H and D, respectively.

4 Summary

The neutral particle analyzer (NPA) is a crucial diagnostic device in Tokamak facilities. The stripping unit is one of the main parts of NPA. A windowless gas tripping room with two differential pipes was adopted to maintain a certain pressure in the parallel electric and magnetic fields(E//B) NPA. The gas pressure distribution of the striping chamber was calculated by combining the computational fluid dynamics (CFD) software ANSYS Fluent, and the Monte Carlo software MolFlow+ for the low-pressure collision-less region in the vacuum chamber. The pressure distribution along the beam direction was obtained for different input pressures. A certain high pressure was achieved inside the stripping room,and a linear decreasing pressure was obtained inside the differential pipes and the entrance and exit holes. A pressure more than two orders smaller was obtained in the vacuum chamber than that inside the stripping room.

Based on the pressure distributions calculated by ANSYS Fluent and MolFlow+, the stripping efficiency of the stripping chamber for H and D atoms at incident energies ranging from 20 to 200 keV was studied using GEANT4. The energy losses of H and D after passing through the stripping chamber are small and can be neglected for all the investigated incident energies and input pressures.The scattering loss of H and D atoms on H2gas was studied using the transmission rate (R) of the incident atoms. A slight increasing trend was observed for R as the incident energy increased,but showed an opposite trend as the input pressure (P0) increased. The scattering loss is small(less than 3%)for all the investigated incident energies and input pressures.

A charge state variable was introduced to track the charge state of the particles in the GEANT4 simulation.Adopting the ORNL-recommended charge exchange crosssections with modified hadron reaction physics, the charge state of each particle is traced in the simulation. The behaviors of the charge fractions along the beam direction(z-axis)in the H2gas were investigated for E=20,100,and 200 keV H and D atoms. The stripping efficiency was obtained as the fraction of charge state +1 at the exit hole of the vacuum chamber (z=120 mm). After reaching the maximum value,f+1decreases for P0=20 Pa,and remains flat for P0=100 Pa,but continues to increase slowly for P0= 400 Pa as the incident energy per nucleon increases. f+1quickly reaches a maximum value and maintains the maximum as P0increases at lower E. For larger E, the fractions increase faster at lower P0and reach a maximum at an input pressure of approximately P0= 240 Pa.

According to the combined global efficiency, R×f+1,P0= 240 Pa was found to be the optimum pressure for the maximum global efficiency in the incident energy range investigated.The simulation results would provide a useful guide for practical applications.

Author ContributionsAll authors contributed to the study conception and design. Material preparation, data collection, and analysis were performed by Yuan Luo and Wei-Ping Lin.The first draft of the manuscript was written by Yuan Luo and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.

主站蜘蛛池模板: 久久特级毛片| 手机永久AV在线播放| 永久免费精品视频| 亚洲色图在线观看| 91久久偷偷做嫩草影院精品| 亚洲一区二区精品无码久久久| 午夜不卡视频| 欧美在线一二区| 欧美成人a∨视频免费观看| 好紧好深好大乳无码中文字幕| h网站在线播放| 成年网址网站在线观看| 再看日本中文字幕在线观看| 国产精品中文免费福利| 日韩欧美中文字幕在线韩免费| 国产高清不卡| 国产精品网址你懂的| 久热这里只有精品6| 毛片一区二区在线看| 一本无码在线观看| 久久一本日韩精品中文字幕屁孩| 久久大香伊蕉在人线观看热2| 久操线在视频在线观看| 国产日韩精品一区在线不卡| 免费无码AV片在线观看国产| 亚洲国模精品一区| 激情综合五月网| 狠狠色狠狠色综合久久第一次| 伊人久久婷婷| 国内毛片视频| 国产色伊人| 亚洲中字无码AV电影在线观看| 国产福利拍拍拍| 福利国产微拍广场一区视频在线| 91精品小视频| 丰满人妻久久中文字幕| 色天天综合久久久久综合片| 日韩在线成年视频人网站观看| 亚洲乱码在线视频| 日韩成人午夜| 91人人妻人人做人人爽男同| 欧美综合中文字幕久久| 国产亚洲欧美日韩在线一区| 国产国产人免费视频成18| 在线观看91精品国产剧情免费| h网站在线播放| 欧美亚洲激情| 成年人久久黄色网站| 亚洲网综合| 天天色天天综合网| 亚洲色无码专线精品观看| 国产成人一区免费观看| 国产成人乱无码视频| 国产精品亚洲va在线观看| 国产精品无码影视久久久久久久| 久久亚洲国产一区二区| 中文字幕va| 婷婷开心中文字幕| 麻豆国产精品| 91精品国产自产在线老师啪l| www.91在线播放| 毛片久久久| 久久精品国产一区二区小说| 久久综合色88| 亚洲黄网在线| 日韩欧美中文在线| 自慰网址在线观看| 日本欧美一二三区色视频| 国产91九色在线播放| 欧美一级99在线观看国产| 在线观看精品自拍视频| 国产日韩AV高潮在线| 国产高清国内精品福利| 婷婷六月综合网| 亚洲不卡无码av中文字幕| 一级毛片在线播放免费观看| 青青草国产精品久久久久| 婷婷色丁香综合激情| 综合色在线| 精品欧美一区二区三区久久久| 国产精品香蕉| 真实国产乱子伦高清|