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

Transitional wave configurations between Type III and Type IV oblique-shock/bow-shock interactions

2023-04-22 02:04:30JunPENGShuaiLIFanYANGMingyueLINGuilaiHANZongminHU
CHINESE JOURNAL OF AERONAUTICS 2023年3期

Jun PENG, Shuai LI, Fan YANG, Mingyue LIN,*, Guilai HAN,Zongmin HU

a State Key Laboratory of High-temperature Gas Dynamics, Institute of Mechanics, Chinese Academy of Sciences, Beijing 100190, China

b The System Design Institute of Mechanical-Electrical Engineering, Beijing 100854, China

c School of Engineering Science, University of Chinese Academy of Sciences, Beijing 100049, China

KEYWORDS Shock interactions;Transitional configuration;Aerodynamic heating;Shear layer;Mach interaction

Abstract The interactions of oblique/bow shock waves are the key flow phenomena restricting the design and aerothermodynamic performance of high-speed vehicles.Type III and Type IV Shock/Shock Interactions(SSIs)have been extensively investigated,as such interactions can induce abnormal aerodynamic heating problems in hypersonic flows of vehicles.The transition process between these two distinct types of shock/shock interactions remains unclear.In the present study,a subclass of shock/shock interaction configuration is revealed and defined as Type IIIa.Type IIIa interaction can induce much more severe aerodynamic heating than a Type IV interaction which was ever reported to be the most serious in literature.The intense aerodynamic heating observed in this configuration highlights a new design point for the thermal protection system of hypersonic vehicles.A secondary Mach interaction between shock waves in the supersonic flow path of a Type III configuration is demonstrated to be the primary mechanism for such a subclass of shock/shock interaction configuration.

1.Introduction

The interactions of oblique/bow shock waves,as key flow phenomena around hypersonic vehicles, have been investigated in the past few decades.The overall interaction configuration changes diversely according to the strengths and the geometric parameters of the intersecting shock waves.As a consequence,the aerodynamic loads over the vehicle surface in the vicinity of the interaction region may vary significantly when the interaction pattern changes.1Edney2–4conducted tunnel tests as well as theoretical analyses, and defined six types of oblique/bow Shock-Shock Interaction (SSI), i.e., Type I to Type VI.His research revealed abnormally high aerothermal loads induced by the Type IV SSI.2,5,6The aforementioned studies confirmed that the most severe pressure load and heat transfer rate occur at the point where the jet impinges at the surface of the blunt-body in a Type IV shock-shock interaction.Based on the experimental data for such an interaction configuration,the enhancement factor of the peak heat transfer rate had been approximately correlated with the peak pressure amplification,jet width, and the blunt radius.2–5Further experimental investigations conducted by Grasso et al.7,8indicated that peak heat flux arises when the jet strikes the model surface perpendicularly.

Both experimental studies5–9and numerical simulations10–16revealed that Type IV interactions feature intrinsic unsteadiness.The transition between Type IV and Type IVa interaction patterns was visualized in a run of tunnel tests.7The underlying mechanism was not clear and was supposed to be either the hysteresis phenomenon associated with the transition between regular and Mach reflections of shock waves17or the transient phenomenon occurring with the development of the test flow.Detailed numerical simulations conducted by Windisch et al.13–16indicated that the unsteadiness in wall pressure or heat flux is mainly correlated with the jet unsteady mechanism.Non-monotonic dependence of peak heat flux on the obliquebow shock intersecting point was revealed in the tests.8,18When Type IV interaction turns into Type III accompanied by the disappearance of expansion fan in the jet structure, a second heat flux maximum is formed.Unsteady phenomena in the overall Type V SSI configuration were revealed by computational investigations on hypersonic double-wedge flows.19–21Here, the transition between regular reflection and Mach reflection of shock waves and the shock wave-shear layer interaction are the key issues to give rise to the unsteadiness aforementioned.When the transition between different SSI patterns occurs, high-temperature effects may lead to remarkable variation in the transition conditions.22–28

Since the starting studies of Edney2–4on the oblique/bow shock-shock interactions which lead to the definition of six types of SSI, extensive investigations either by experiment or computation have been reported in the literature.However,to the best of the authors’knowledge, the transition criteria among different SSI types have not been obtained to date due to the complexity of the flow field.A preliminary study using a machine learning method from Peng et al presented the SSI transition criteria for Type III-IV, and Type IV-IVa,respectively.29However, Type IV and Type IVa interactions are hard to distinguish due to the inherent unsteady features.In addition,there are significant differences in the overall wave configurations of Type III and Type IV interactions as shown in Fig.12.In Fig.1,IS repesents incident shock wave,BS repesents bow shock wave, TS repesents transmitted shock wave,SL repesents shear layer, NS repesents normal shock wave.It is clear that Type III and Type IV configurations cannot turn from one into the other like an on-off switch and there should be a transition process between them.The nonmonotonic phenomenon revealed in the heat flux measurements conducted by Boldyrev et al.8implies that the SSI configuration within the transition path should be significantly different from either Type III or Type IV interactions.

In this paper,the transitional flow structures between Type III and Type IV oblique shock/bow shock interactions are studied by numerical simulation.The main objective of the present work is to reveal the distinctive flow features during the transition process as well as the resulting surface quantities such as pressure and heat flux.

2.Computation setting

2.1.Model definition

The present work focuses on shock/shock interactions during the transition process between Type III and Type IV SSI configurations.According to the parameter-correlation study conducted by Peng et al.,29the control parameters for such a shock-dynamics dominant problem includes three variables,i.e., the freestream Mach number, Ma∞, the incident shock angle, β, and the geometrical parameter Irwhich determines the spatial relation between the incident shock wave (IS) and the bow shock wave (BS).All the control parameters are schematically depicted in Fig.2.Two of the abovementioned parameters, i.e., Ma∞and β can uniquely determine the flow solutions in regions (2), (3) and (4) based on the classically two- and three-shock theories proposed by von Neumann.30φ represents the angle of circumferential position on the circular surface, and R is the radius of the cylinder.Parts of the remaining flow regions can be theoretically obtained when Iris given simultaneously.However,there are no theoretical solutions to the coordinates of triple points A and B by far to the best of the authors’knowledge.Thus,such SSI problems cannot be solved theoretically, but can be worked out by numerical simulation.10–16,29

Fig.1 Type III and Type IV shock/shock interaction configurations according to Edney’s definition.2

Fig.2 Control parameters for flows of Type III and Type IV shock-shock interactions.

Fig.3 Computational domain, boundary specifications and overall SSI configuration shown by Mach number contour (Grid:1600 × 800, Ma∞= 8, β = 21°, Ir = 0.17).

Two-dimensional laminar compressible Navier-Stokes equations are supposed to be the governing equations for the flows studied hereafter.For the numerical algorithms, the second-order TVD (Total Variation Diminishing) scheme based on a non-linear Riemann solver named HLLC(Harten-Lax-van Leer Contact)31is applied to solve the convective terms.This computational code based on the finitevolume method has been validated in the previous investigations29for the simulations of hypersonic flows especially Type IV oblique shock/bow shock interactions.It is verified that the distributions of heat flux calculated by the code used in this paper are in good agreement with the experimental data,including the locations and magnitudes of peaks of heat flux.

2.2.Convergence study

Fig.4 Grid independence evaluation for simulations of oblique shock/bow shock interactions (Ma∞= 8, β = 21°, Ir = 0.17).

The computational domain used in the following simulations is given in Fig.3.The domain is discretized by structured grid while the grid is refined both in the shock-shock interaction region and along the cylinder surface.The SSI configuration for the following conditions (Ma∞= 8, β = 21°, Ir= 0.17)is depicted in the figure by the contour of flow Mach number.For all the following computations, the domain should be adjusted carefully to maintain that the sonic lines downstream the bow shock waves should be completely accommodated inside.Boundary specifications are also shown in Fig.3.The boundary conditions are set to user-defined inflow and outflow at the left and right of the computational domain respectively.Nonslip and isothermal conditions are used for the solid wall at a fixed temperature of 300 K.For all the computations in the present study, the static temperature and pressure of the freestream airflow are given as T∞=200 K, P∞=428 Pa,respectively.The freestream flow conditions are selected according to the test flow conditions of a large-scale detonation-driven shock tunnel located at the State Key Laboratory of High-temperature Gas Dynamics of Institute of Mechanics of Chinese Academy of Sciences in Beijing.32–34The unit Reynolds number of the flow is 1.27 × 106/m.When the model is small and the boundary layer is not fully developed, the flow is generally laminar.Thermal equilibrium gas is utilized in this study, thermodynamic and transport properties of which are calculated by polynomials of temperature.35The polynomials reflect the vibration excitation of air molecules at high temperatures.Of all the cases studied in this paper, the maximum temperature of the flow field is about 2800 K.The degree of the dissociation of molecular is small,if any.Therefore,turbulence and other dissipative phenomena are not considered.

Fig.6 Peak heat flux and peak pressure during the transition of Type III-IV shock/shock interaction (Ma∞= 8, β = 21°).

Fig.5 SSI transition from Type III to Type IV with rising geometrical parameter Ir (Ma∞= 8, β = 21°).

Fig.7 Sketch and zoom-in flow structure of Type IIIa shock/shock interaction (corresponding to Fig.5(c) where Ma∞= 8, β = 21°,Ir = 0.16).

3.Results

3.1.Type IIIa shock/shock interaction

Different SSI types can be obtained by changing Irwhile maintaining Ma and β constant.As aforementioned,the flow states in each region around triple point A(see Fig.1 or Fig.2)can be determined once Ma and β are given.The flow states include the thermodynamic and kinematic parameters in regions (2),(3) and (4).However, triple point A must move as Irchanges to feel the disturbance which is generated by impingement of the shear layer, SL1, over the cylinder surface and propagates upstream in the subsonic flow of region(4).Consequently,the SSI configuration can vary from Type I to Type VI with the increase of Ir.Fig.5 shows the wave structures of different SSI types when Irincreases from 0.12 to 0.5 while maintaining Ma∞= 8, and β = 21°.

Fig.8 Shock polar system for Type IIIa interaction (corresponding to Fig.5(c)) (Ma∞= 8, β = 21°, Ir = 0.16).

Fig.9 Typical profiles of heat flux and pressure) along the cylinder surface for Type IIIa (Ir = 0.2) and Type IV (Ir = 0.35)interactions (Ma∞= 8, β = 21°).

According to the sketch in Fig.1,it is easy to find a typical Type III configuration of SSI as depicted in Fig.5(a) for the simulation when Ir=0.12.Fig.5(b)follows the main features of Type III interaction configuration while a little different from Fig.5(a) appears in the vicinity of point B (see labels in Fig.1) which becomes a triple point of a Mach reflection.In Fig.5(a), B is a point of intersection of two oblique shock waves of the same family.On the other hand,Type IV interactions for Ir= 0.3, 0.4 and 0.5 are given in Fig.5(g), Fig.5(h)and Fig.5(i), respectively.As Type IV SSI is reported to be unsteady mainly due to jet impingement,13–16the present simulations for the aforementioned three conditions appear unsteady too.Each frame here shows a transient SSI structure for each condition of Ir.

Fig.10 Peak heat flux and pressure for Type IIIa interaction with oscillation (Ir = 0.25, Ma∞= 8, β = 21°).

The SSI structures depicted in the remaining frames, i.e.,Fig.5(c)-Fig.5(e), are different from either Type III or Type IV in the vicinity of point B (see labels in Fig.1) when the supersonic flow in region (3) arrives at the cylinder surface.This kind of SSI structure occurs within a relatively small parameter domain of Ir= 0.16 - 0.2 for the given conditions β=21°and Ma∞=8.Here,this special situation is classified as Type IIIa shock/shock interaction to make difference from Edney’s terminology of Type III.Obviously, the overall structure of Type IIIa interaction looks similar to that of a general Type III interaction while appears subtle differences in the vicinity of interaction region near the cylinder surface from the latter.Therefore, Type IIIa configuration is a subclass of Type III interaction patterns.It was reported in the literature2–8that Type IV interference results in the most critical dynamic and thermal loads associated with perpendicular jet impingement on the surface.From the present study as depicted in Fig.6, the studied SSI structure of Type IIIa can lead to more severe aerodynamic heating than Type IV.In Fig.6,the peak heat flux and pressure are nondimensionalized by the heat flux, Qst(3.5 × 105W/(m2?s)) and pressure, Pst(3.58×104Pa),at the stagnation point of the cylinder surface without shock interference, respectively.

Fig.12 Sketch of double-InMR structure of unsteady Type IIIa SSI at t = 2800 μs (Ir = 0.25, Ma∞= 8, β = 21°).

Fig.11 Transient flow structures of Type IIIa interaction with oscillation (Ir = 0.25, Ma∞= 8, β = 21°).

Fig.13 Peak heat flux and its location along the cylinder surface(labels a to j correspond to each frame in Fig.11).

Fig.14 Peak heat flux and pressure for Type IV interaction with oscillation (Ir = 0.4, Ma∞= 8, β = 21°).

The distinctive feature of Type IIIa SSI is schematically depicted in Fig.7 with a zoom-in illustration (see Fig.7(b))for the local wave pattern in the vicinity of triple point B.The theoretical solution of this special kind of wave system can be obtained by shock polar analysis as shown in Fig.8 for the case of Ir= 0.16.When the supersonic flow between the shear layer SL1 and the transmitted shock TS, both of which originate from triple-point A approaches the cylinder surface, an oblique shock wave, i.e., CE (see Fig.7(a)) is formed to match the boundary conditions along the surface.The interaction of the oblique shock CE and the reflected shock BD originating from triple-point B results in a bi-Mach reflection structure and forms new triple-points, D and E.Two shear layers, SL3 and SL4, come into being and form a convergent-divergent flow tube.Such a bi-Mach reflection structure is defined as double-DiMR by Li et al.36and can maintain steady in steady supersonic flows.As a consequence,the flow structures of Type IIIa interaction as depicted in Fig.5(c)-Fig.5(e)are always steady unlike Type IV interaction which is substantially unsteady.

Shear layers SL1 and SL3 assemble a supersonic jet developing along the cylinder surface.The reflected shock wave of the triple-wave structure at point E is reflected repeatedly between the two shear layers.As schematically depicted in Fig.8,the flow in region(7)is close to the sonic flow condition;a Mach disc appears at point F(see Fig.7(b))in the supersonic jet.The interaction between the aforementioned Mach disc with the shear layer SL1 leads to an extremely high heating load shortly downstream point F.As the parameter Irincreases with the translational moving of the incident shock,the aforementioned triple-wave structure will shrink and the supersonic jet flow surrounded by SL1 and SL3 will get thinner.When Ir= 0.2, the peak heat flux becomes around twenty-six times the corresponding value at the stagnation point,as can be seen in Fig.6,i.e.,Qpeak(Ir=0.2)/Qst≈26.Such an amplification factor of aerodynamic heating is much higher than the maximum amplification factor induced by the Type IV interaction when Ir= 0.35.The distribution profiles of heat flux and pressure along the cylinder surface for Type III, Type IIIa and Type IV interactions are combined in Fig.9 for a direct comparison.

3.2.Type IIIa shock/shock interaction with oscillations

Fig.15 Transient flow structures of Type IV interaction with oscillation (Ir = 0.4, Ma∞= 8, β = 21°).

The Type IIIa SSI structures as depicted in Fig.5(c)-Fig.5(e)are steady with the aid of a convergent-divergent flow tube downstream the Mach stem DE (see Fig.7(b)).With the increase of the geometric factor Ir, the tripe-wave structure at point E will shrink and disappear completely at a critical condition,e.g.,Ir=0.2 for Ma∞=8,β=21°.At such a critical point, the shock wave DE intersects with the shear layer SL1 perpendicularly and triggers unsteadiness.As a consequence, the overall structure of such a special Type IIIa interaction undergoes approximately periodic oscillation which can be reflected by the history of the peak heat flux or the peak pressure along the cylinder surface as given in Fig.10.The frequency of the oscillation is around 11 kHz.The transient flow structures within one cycle (corresponding to the time range enclosed by the dash-dotted box in Fig.10) are depicted in Fig.11.

The unsteadiness of Type IIIa structure as shown in Figs.10 and 11 are primarily induced by the interaction of the shock wave DE and the shear layer SL1.Fig.12 schematically shows the transient wave configuration corresponding to the transient structure of Type IIIa interaction at t=2800 μs depicted by the first frame of Fig.11.The aforementioned shock-shear layer interaction results in a λ-shaped structure or a Mach reflection structure at triple-point E.The two shear layers,i.e.,SL3 and SL4,assemble a divergent flow tube locally,indicating that the transient double-Mach reflection structure connected by DE is unsteady in nature.Such a shock structure is denoted as double-InMR by Li et al.36which is not physical in steady flows.Shocklets or perturbation waves propagate upstream during the unsteady evolution of the double-Mach reflection structure and influence the triple-points A and B.However, the influences are so gentle that the trajectories of points A and B during an oscillation period are confined in tiny domains as demonstrated in Fig.12 by the line segment and triangle, respectively.

The peak heat flux and the corresponding location (represented by parameter φ) at each time transient associated with each frame in Fig.11 are depicted in Fig, 13.The peak varies from around ten to twenty-nine times of the stagnation heat flux, while its location varies within 6 degrees of φ.The peak heat flux reaches a maximum at time transient of 2840 μs shortly after the birth of a new double-I nMR structure.It should be noted that the peak heat flux given in Fig.6 for the case of Ir= 0.25 is the temporal average as it varies with time.Oscillations of the wave structure leads to a decrease in the averaged heat flux at each point of the cylinder surface.

The fact that a Type IV interaction is substantially unsteady is widely reported in the literature.5–16Unlike a Type III or IIIa interaction in which the shear layer SL2 does not reach the cylinder, both primary shear layers of a Type IV interaction, i.e., SL1 and SL2 (see the definition in Fig.1(b)),strike the cylinder surface and induce a supersonic jet impingement.Such an impingement can result in severe heating and pressure loads to the surface accompanied by flow oscillations as depicted in Fig.14 for the case of Ir=0.4 for instance.The frequency of the oscillation is around 10 kHz for this case,within which the main flow structures are shown in Fig.15.As the supersonic jet flow in region (3) is approximately perpendicular to the cylinder surface, a normal shock wave, as denoted by NS in Fig.1(b), appears and interacts with the shear layers.The normal shock wave varies in length and location slightly and becomes the shortest and nearest to the surface at a time transient of about 2865 μs (see Fig.15) when the peak heat flux reaches maximum.

As discussed hereinbefore, the unsteady behavior of Type IIIa interaction is different from that of Type IV interaction.The primary difference is associated with the local wave structures at the spot where the supersonic jet impinges the cylinder surface.The historic variation of peak heat flux and surface pressure as respectively shown in Figs.10 and 14 reflect the aforementioned difference.

Fig.16 Profiles of surface heat flux for different types of shockshock interaction configuration.

3.3.Discussion

The Type IIIa SSI studied here is a subclass wave configuration of Type III interaction that generally results in the most severe aerodynamic heating in the interference region along the cylinder surface.Such a wave configuration occurs during the transition domain between Type III and IV SSI,and combines the primary flow features of the latter two interaction patterns.The peak heat flux induced by the Type IIIa SSI within a special range of parameter Irmay become much larger than that in a Type IV SSI with the same Ma∞and β as can be seen in Fig.16.The most serious aerodynamic heating conditions may be misestimated if one applies interpolation among predicted solutions of Type III and IV interaction.However,a Type IIIa SSI occurs within a very small parameter domain of Iras shown in Fig.16.Therefore, the flow structures of Type IIIa interaction can hardly be visualized by hypersonic tunnel tests and have not been paid careful attention to in open literature, according to the best of the authors’knowledge.

As can be seen in Figs.7 and 8, the Mach type interaction between oblique shock waves CE and BD gives birth to a Type IIIa SSI configuration.Here,shock CE is induced by the interaction of the supersonic flow in region (3) while BD is the reflected shock wave of the Mach reflection wave configuration associated with Triple-point B.When shear layers SL3 and SL4 assemble a convergent-divergent flow tube as depicted in Fig.7, the resulting Type IIIa SSI configuration maintains steady, and is not vice versa as illustrated by Figs.10 to 13.Nevertheless, a regular interaction instead of a Mach interaction between shock waves CE and BD,is theoretically accepted as shown in Fig.17 where a smaller incident shock angle is given,i.e.,β=15°.Solution states(8,9)in the shock polar diagram represent the regular interaction aforementioned.Such a Type IIIa SSI configuration always maintains steady since the flow is always supersonic downstream the regular interaction.Fig.16 (c) indicates that the peak heat flux in a Type IIIa SSI with regular shock/shock interaction in the supersonic flow path between shear layers can exceed that of a Type IV configuration.

Borovoy et al.revealed a non-monotonic distribution of the peak heat flux with the peak position by experimental measurements using specially arranged thermal couples.8,18The first maximum peak heat flux is associated with the supersonic jet impingement in a Type IV shock/shock interaction.A great number of studies have been concentrated on such kind of interaction and thrown light on the aerodynamic mechanisms.Borovoy et al.stated that the second maximum is primarily due to the disappearance of an expansion fan in the supersonic jet when Type IV flow turns into Type III flow.According to the present study, the primary mechanism of the second maxima of peak heat flux is related to subclass of Type III shock/shock interaction, i.e., Type IIIa, occurring during the transition between Type III and Type IV.In addition, the main inducement of a Type IIIa configuration is the Mach interaction or regular interaction between two shock waves of different families in the supersonic flow path between shear layers as depicted in Fig.5(d)-Fig.5(f)or Fig.17(b),respectively.However, due to the insufficient resolution of flow visualization in experiments or of numerical simulation.8,18around twenty years back, such a subtle flow structure could not be well captured.

4.Conclusions

In the present study, the flow patterns during the transition between Type III and Type IV oblique shock/bow shock interactions over a cylinder model are numerically simulated.A geometric ratio Iris introduced to represent the spatial relation between the incident shock wave and the cylinder surface.Iris one of the control parameters of the studied problems along with the incident shock angle,β,and the freestream flow Mach number,Ma∞.By carefully adjusting parameter Irwhile maintaining the rest control parameters fixed, different types of shock/shock interaction configurations are obtained in the simulations.SSI configuration can vary from Type I to Type VI when Irchanges,which has been widely reported.However,a transition SSI configuration between Type III and Type IV which has not been paid careful attention to is found in this work,when Iris located inside a very small range of parameter domain.

Fig.17 Type IIIa SSI configuration with regular shock-shock interaction in supersonic jet flow (Ir = 0.125, Ma∞= 8, β = 15°).

This subclass of Type III shock/shock interaction configuration is revealed and defined as Type IIIa SSI.Such an interaction configuration occurs during the transition process between Type III and Type IV SSI configurations and appears different flow features in a subtle region from both of the latter.A Mach interaction or regular interaction between two shock waves of different families in the supersonic flow path between shear layers of the primary triple-points is found to be the main mechanism for the birth of a Type IIIa SSI configuration.Such a local shock/shock interaction can be either steady or unsteady depending on whether a convergentdivergent flow tube is formed or not downstream the local interaction.The aforementioned local shock/shock interaction results in a secondary shock wave/shear layer interaction along the cylinder surface which may induce an extremely high heat flux.The most severe aerodynamic heating condition occurs when a Type IIIa configuration is steady.In fact,the transient peak of heat flux in the unstable Type IIIa interaction(such as the case given in Fig.11) can be higher than a Type IV interaction.However, the location of the peak heat flux varies along the cylinder surface which may significantly decrease the time-averaged peak of heat flux.

This is a primary study using numerical techniques for high-temperature gas flows without considering turbulence in the simulations.In the future, turbulent simulation or experimental visualization will be conducted for further studies on such special type of shock/shock interaction.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

This study was co-supported by the National Key Research and Development Plan of China (No.2019YFA0405204),and the National Natural Science Foundation of China(Nos.12172365, 12072353 and 12132017).

主站蜘蛛池模板: 欧美www在线观看| 亚洲AV无码精品无码久久蜜桃| 2021天堂在线亚洲精品专区| 亚洲中文精品人人永久免费| 国内99精品激情视频精品| 99成人在线观看| 凹凸精品免费精品视频| 亚洲色图欧美视频| 91小视频版在线观看www| 国产成人精品2021欧美日韩| 欧美在线精品一区二区三区| 亚洲人成人伊人成综合网无码| 精品国产电影久久九九| 日本不卡视频在线| 日韩高清无码免费| 人妻丝袜无码视频| 老司机精品一区在线视频| 成人亚洲国产| 国产在线自揄拍揄视频网站| 色偷偷综合网| 欧美性猛交一区二区三区| 中文字幕调教一区二区视频| 国产激爽爽爽大片在线观看| 亚洲av无码片一区二区三区| 一本大道无码高清| 热99精品视频| 亚洲无码91视频| 久久综合激情网| 97se亚洲| 国产精品视频观看裸模 | 再看日本中文字幕在线观看| 影音先锋丝袜制服| 2020精品极品国产色在线观看 | 中文无码精品A∨在线观看不卡 | 国产激情影院| 中文字幕在线一区二区在线| 伊人丁香五月天久久综合 | 久久伊人久久亚洲综合| 黄色免费在线网址| 精品国产一区91在线| 欧美色图第一页| 国产成人无码Av在线播放无广告| 免费人欧美成又黄又爽的视频| 亚洲av日韩av制服丝袜| 亚洲成人www| 99久久精品无码专区免费| 国产91色在线| 亚洲第一天堂无码专区| 日本一本在线视频| 久视频免费精品6| 亚洲国产AV无码综合原创| 亚洲欧美日本国产综合在线| 国产乱子伦无码精品小说| 91色爱欧美精品www| 日韩中文无码av超清| 久久96热在精品国产高清| 日韩欧美网址| 国产综合网站| 欧美成一级| 98超碰在线观看| 99热最新网址| 国产精品漂亮美女在线观看| 999福利激情视频| 亚洲一区二区在线无码| 精品视频在线观看你懂的一区 | 精品伊人久久久香线蕉 | 国产亚洲精品yxsp| 美女被躁出白浆视频播放| 在线中文字幕日韩| 国产簧片免费在线播放| 伊人久久久久久久| 成人中文字幕在线| 亚洲国产成人久久精品软件| 亚洲精品第五页| 日本久久网站| 五月天福利视频| 亚洲熟女偷拍| 欧美一级视频免费| 色噜噜在线观看| 亚洲精品桃花岛av在线| 在线永久免费观看的毛片| 国产69精品久久久久孕妇大杂乱|