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

    Decomposition of mean skin friction in incident shock wave/turbulent boundary layer interaction flows at Mach 2.25

    2023-10-25 12:12:24JunyiDUANFulinTONGXinliangLIHongweiLIU
    CHINESE JOURNAL OF AERONAUTICS 2023年9期

    Junyi DUAN, Fulin TONG, Xinliang LI, Hongwei LIU

    a LHD, Institute of Mechanics, Chinese Academy of Sciences, Beijing 100190, China

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

    c State Key Laboratory of Aerodynamics, China Aerodynamics Research and Development Center, Mianyang 621000, China

    KEYWORDS

    Abstract The evolution characteristics of the mean skin friction beneath the supersonic turbulent boundary layer that interacts with incident shock waves at Mach 2.25 are analyzed using Direct Numerical Simulation(DNS).The separated and attached boundary layers in the interaction region that respectively correspond to 33.2° and 28° incident shock angles are considered.The mean skin friction recovery rate for the separated boundary layer is much gentler and distinctly less than that for the attached case where the skin friction completes its recovery within one boundary layer thickness.The novel mean skin friction decomposition method for compressible flows proposed by the recent research is applied in the interaction region to investigate the internal evolution characteristics quantitatively.The results reveal that the three decomposition components are distinctly unequal between the two cases.The contributions of the turbulent motions at different scales to the associated term are focused on using empirical mode decomposition technology.It indicates that the outer large-scale structures dominate separation and reattachment regions,while contributions from inner small-scale structures are limited.In contrast, contributions from the outer largescale structures are dramatically reduced in the attached case,which results in the outer large-scale and inner small-scale motions being of equal importance.

    1.Introduction

    Owing to its ubiquitous presence and profound significance in the aeronautical field, Shockwave/Turbulent Boundary Layer Interactions (STBLIs) have been widely studied over the past half century.The encountered phenomena include the maximum surface fluctuating pressure, severe wall heat flux peak,and strong unsteadiness,which cause significant adverse effects on aircraft and component performance.1,2There has been some remarkable progress with these fundamental phenomena, especially over the last 20 years, through the fastdeveloping high-resolution experimental measurement facilities and High-Performance Computing (HPC) technologies.However, many challenges urgently need to be solved,3such as shock wave low-frequency oscillations, turbulence amplification, and quick and accurate prediction methods.

    Wall Shear Stress (WSS), together with Wall Heat Flux(WHF) and wall pressure, which are wall-quantities closely related to the aerodynamic design and thermal protection of the aircraft, is of obvious importance in studying wallbounded turbulent flows.4,5A considerable number of empirical or semi-theoretical formulations between mean WSS and Reynolds number in smooth-wall-bounded turbulent flows have been built up.6–9Besides, various studies have been conducted to improve understanding of the statistics and structure of fluctuating WSS.10–12Extremely high WSS events in the incompressible Turbulent Boundary Layer(TBL)were studied by Pan and Kwon,10and the relationships between the extreme positive WSS events and the strong sweeping events originating from the outer layer positive large-scale motions were detected.Statistical characteristics of WSS and WHF fluctuations in supersonic TBL were investigated by Tong et al.,11and the similar mechanism for WSS in compressible flows was found through conditional average analysis.

    The WSS has also been widely investigated in STBLI flows.Various kinds of methods, for instance, the buried wire gages,13laser interferometric instruments,14and global interferometry skin friction techniques,15were used to make measurements in the two- or three-dimensional interaction regions.The WSS was calculated theoretically for the normal shock wave/flat plate TBL interactions by asymptotic methods.16The statistical characteristics of WSS fluctuations,including the probability density function, weighted powerspectrum density,and space–time correlation, are significantly changed through the interaction region.17,18

    It is well known that the enhanced skin friction and heat transfer in the Reattached Boundary Layer (RBL) can be several times the values of the incoming undisturbed TBL,19which deteriorate significantly for hypersonic flows.20,21The excessive skin friction causes reduced aerodynamic efficiency and combat radius and increased fuel investments.To the authors’ knowledge, there is minimal systematic research about the skin friction evolution characteristics of the TBL under shock interactions.Recent research has proven that skin friction is related to turbulent motion within the boundary layer through mathematical derivations.22,23As noted by Tong et al.,24,25the turbulent vortex structures are significantly strengthened due to the shock interactions.Turbulent motion within the RBL is more complex and characterized by regenerated small-scale structures in the inner layer, and large-scale structures in the outer layer originating from the upstream interaction region.The relationship between the mean skin friction and complex vortical structures within the RBL needs further study.

    The primary objective of this study is to further improve the physical understanding of the skin friction relaxation process downstream of the interaction region in STBLI flows via Direct Numerical Simulation(DNS).To this end, the mean skin friction decomposition in the compressible boundary layer proposed by Li et al.26is combined with bidimensional Empirical Mode Decomposition (EMD) to analyze the DNS dataset of the 33.2°/28° incident shock that interacts with a turbulent boundary layer at Mach 2.25.This method can provide a physical interpretation for different contributing terms while quantifying the contributions associated with specific spanwise length scales.

    The remainder of this paper is organized as follows.Section 2 fully describes the simple flow configuration selected for the DNS and the adopted computational setup.Section 3 discusses the evolution characteristics of the mean skin friction in STBLI flows through analyses of the ensemble-averaged field, mean skin friction decomposition, and empirical mode decomposition.Conclusions are drawn in Section 4.

    2.Numerical strategy

    2.1.DNS setup overview

    The three-dimensional compressible conservative Navier-Stokes equations are nondimensionalized using the freestream velocity u∞, density ρ∞, temperature T∞, dynamic viscosity μ∞, and unit length L = 1 in-1(1 in = 0.0254 m), and are directly numerically solved using the high-order finite difference solver Opencfd-SC developed by Li et al.,27which has been widely applied in compressible turbulent flows and STBLI flows.28–32The subscript‘‘∞”represents the properties of the freestream flow.The perfect gas equation and Sutherland’s equation are employed to close the governing equations.The specific heat ratio and molecular Prandtl number are taken as γ=1.4 and Pr=0.7 respectively.For the spatial discretization of the governing equations, the convection terms are solved using a bandwidth-optimized WENO scheme combined with the limiters,33which has been generally used in STBLI problems.17,34The diffusion terms are solved using an eighth-order central scheme.In addition, a third-order Runge-Kutta method is used for time integration.

    A two-dimensional sketch of the computational domain is shown in Fig.1.The original point of the Cartesian coordinate system is at the nominal shock incident point xo, and the spatial coordinates are normalized by the turbulent boundary layer thickness δref=0.078 in-1at the reference station,where its non-dimensional coordinate is x= -4.0.Unless otherwise stated, the spatial coordinates (x, y, z) are normalized by δref.The laminar boundary layer, which has identical freestream parameters and wall temperature conditions, is prescribed at the inlet boundary.A fully developed TBL is generated through the transition of the incoming laminar boundary layer under the disturbance of wall blowing and suction.The region for these disturbances ranges from x= -47.5 to-41.1 downstream of the inflow laminar boundary layer, and the setup of wall-normal disturbance velocity vbsis the same as suggested by Pirozzoli et al.35and Fang et al.36Additionally,an isothermal non-slip boundary condition with a fixed wall temperature of Tw= 1.9T∞is imposed.Hereafter, the subscript ‘‘w”denotes the properties at the wall surface.The non-reflective boundary condition is adopted for the computational domain’s outlet and top far-field boundaries to protect it from staining by reflected disturbances.The single-point R-H relations are imposed at the top boundary around xs= - 9.8 or - 12.1 to generate an incident shock with shock angles of α=33.2°or 28°.The incident shock generated by this method is easier to control than the wedge angle used by Zhong et al.37Periodic boundary conditions are used in the spanwise direction.

    Fig.1 Sketch of computational domain in x-y plane.The grid points are plotted with only every tenth and fourth point in x and y directions, respectively.

    As shown in Fig.1,the extent of the computational domain is Lx×Ly×Lz=68.43δref×6.31δref×2.25δref.A structured Cartesian mesh composed of 3700 × 300 × 250 grid points in the x,y,and z directions discretizes the computational domain for both cases.The mesh in the streamwise direction is progressively refined in the transition region ranging from x = - 51.3 to - 12.8 (total of 600 grid points) and is uniformly distributed in the well-resolved region ranging from x = - 12.8 to 12.8 (total of 3000 grid points).In these two regions, the wall-normal mesh is clustered toward the wall using a hyperbolic tangent mapping with 193 grid points within the TBL and 280 grid points within the computational domain (y ≤3.85).In the top and outlet buffer regions, the mesh is rapidly coarsened in the wall-normal (20 grid points)and streamwise(100 grid points)directions respectively to further perfect the dominating computational region from being fouled by reflected disturbances.In addition, 250 grid points are equally distributed for the spanwise computational domain 0 ≤z ≤2.25.As evaluated at the reference station, the mesh resolution in the streamwise direction, wall-normal directions at wall surface and TBL edge, and the spanwise direction are Δx+=6?23,Δy+w=0?65,=11?66,and Δz+=6?54,respectively.Note that the superscript‘‘+”represents the variable in the inner scaling, i.e., it is normalized with the viscous scales computed at the reference station.Research by Poggie et al.38on the effects of mesh resolution on compressible turbulent flows indicates that the present simulation uses the standard mesh to resolve the finest structures in turbulent flows.

    The parameters of the fully developed TBL are similar to those considered by Pirozzoli et al.35,39and Fang et al.36These are a freestream Mach number of Ma∞= 2.25, freestream unit Reynolds number of Re = 635000 in-1, and freestream flow temperature of T∞= 169.44 K.At the reference station,the Reynolds number based on the momentum boundary layer thickness is Reθ= ρ∞u∞θ/μ∞= 3533, and the friction Reynolds number is Reτ= ρwuτδref/μw= 728.

    2.2.DNS validation

    The van Driest transformed mean streamwise velocityin the inner scaling is defined as

    and the Reynolds-averaged streamwise velocity u- in the outer scaling is plotted in Fig.2.Unless otherwise stated, the superscript bar ‘‘–” and tilde ‘‘~” represent the Reynolds average(i.e., time-spanwise average) and Favre average (i.e., densityweighted averagerespectively, with a prime and double prime for the corresponding fluctuations.As shown in Fig.2, the velocity profiles overlap well with the DNS data from Schlatter and O¨ rlu¨9(incompressible TBL at Reθ=2000)and Fang et al.36(similar compressible TBL at Reθ=3700),as well as the experimental data of Bookey et al.40(compressible TBL at Reθ= 2400).Specifically, the profile ofpresents a linear distribution in the viscous sublayer y+< 7 and obeys the log-law with a standard slope of 1/κ =2.44 (κ is the Karman constant) with a cutoff of 5.2 in the log layer 30 < y+< 100.

    where TKE is balanced by the terms of convection C,production of the mean velocity gradient P, turbulent transport T,pressure dilatation Π,viscous diffusion D,mass flux contribution associated with density fluctuations M, and viscous dissipation ε.Detailed formulae for these terms are found in the work of Pirozzoli et al.35and Tong et al.41As shown in Fig.4(a),the transport process for TKE is dominated primarily by the production P, viscous dissipation ε, turbulent transport T, and viscous diffusion D, while the other terms are negligible.In the region of y+< 30, the greatest generation of TKE is transported into the near-wall region via turbulent transport and viscous diffusion, and eventually dissipated as heat in the near-wall region.In the log layer for y+> 30,TKE production is balanced mainly by viscous dissipation.Furthermore, the production-to-dissipation ratio P/ε is illustrated in Fig.4(b), which shows that the ratio attains a maximum value of approximately 2.0 at y+≈10.These results are consistent with those of Schlatter and O¨ rlu¨9and Sun et al.42In particular,the ratio remains constant for 30

    Fig.2 Profiles of van Driest transformed mean streamwise velocity for inner scaling(a)and Reynolds-averaged streamwise velocity -u for outer scaling (b) at reference station.

    Fig.4 Profiles of turbulence kinetic energy budget terms(a)and production-to-dissipation ratio(b)at reference station.The symbols in(a)are the DNS results for Ma=2.25 zero-pressure-gradient turbulent boundary layer with the same flow conditions as Pirozzoli et al.35

    Fig.( 5 D)istribution of mean wall pressure,, with the incident shock angle of α = 33.2°.The dotted line represents the inviscid shock solution.

    3.Results and discussion

    3.1.Instantaneous and mean flow fields

    Contour maps of the instantaneous streamwise velocity u for α = 33.2° and 28° are plotted in Figs.6 and 7, respectively.An incident shock is generated through the single-point R-H relations imposed at the top boundary and interacts with the incoming TBL downstream.As shown in Fig.6(a),a transparent separation region can be observed around the incident shock point.According to Souverein et al.,44a strong interaction occurs for α=33.2°.The incoming TBL is blocked under the adverse-pressure-gradient effects caused by the incident shock, and a large area of low-momentum fluids accumulates in the downstream region of the reflected shock.Multiple small and separated elements, which are highlighted by solid white lines, can be found within the separation bubble.Notably,the reattached turbulent structures downstream of the interaction region are lifted by the low-momentum fluids within a long horizontal distance (at least 5δref).

    In the buffer layer,the canonical high-and low-speed alternating streamwise-elongated streak structures with spanwise widths of approximately Δz+= 80–100 are seen in Fig.6(b), which is consistent with the experimental observations of compressible TBL at Ma=2 by Ganapathisubramani et al.45Despite these streak structures disappear within the interaction region,and they begin to regenerate until 2δrefdownstream of xo.Furthermore, the regenerated streak structures, such as at x=3.0,have a much larger spanwise scale than the upstream TBL.These are associated with Go¨rtler vortices20,46–48and Kelvin-Helmholtz instability43,49in the detached shear layer,as previously observed by Pasquariello et al.50

    Fig.6 Contour maps of instantaneous streamwise velocity u for an incident shock angle α=33.2°in x-y plane(a),x-z plane(b),and y-z plane (c).The white solid lines are the isolines of u = 0.

    Fig.8 Instantaneous turbulence coherent structures around reference station (a) and downstream of interaction region (b) for an incident shock angle of α = 33.2°.

    Fig.9 Instantaneous turbulence coherent structures around reference station (a) and downstream of interaction region (b) for an incident shock angle of α = 28°.

    Fig.10 Distribution of mean wall pressure(a)and mean skin friction coefficient Cf (b).The dashed lines in Fig.10(a) are the inviscid shock solution.

    As reported in Fig.7,the interaction intensity for α=28°is weaker than that for α = 33.2°.The large area for lowmomentum fluids that is obviously observed in Fig.6(a) is absent, and the incoming boundary layer remains nearly attached when passing through the interaction region.Moreover,the flow field structure for α=33.2°,which is completely separated in the interaction region,is distinct from the interactions with α = 28°, where the fluids have minimal separation.In addition, a footprint of the incident shock brands on the incoming streamwise-elongated structures is shown in Fig.7(b), where the lower-speed streaks between × = - 1.0 and 0 are relatively weakened.Downstream of the interaction region, these streaky structures are essentially preserved and the spanwise length is slightly changed, which is utterly different from the 33.2° case.

    The instantaneous vortex structures beneath the fully developed TBL and RBL for both cases are visualized in Figs.8 and 9, respectively, using the Q criterion proposed by Jeong and Hussain.51Coherent structures appear as streamwise vortices upstream of the interaction region, which is known as the leg of horseshoe-like vortices and agrees with the numerical observations of Zhang et al.,52Wu et al.,53and Fang et al.,36as well as the experimental observations by Zhuang et al.54For α=33.2°,there are more vortex structures captured at higher altitudes around x = 3, with larger spatial scales than at the upstream TBL, which is consistent with observations of the streak structures in the near-wall region.However, differences in the vortex structures between the upstream and downstream zone are indistinctive for α = 28°, except that the vortices downstream of the interaction region are slightly more abundant and at approximately the same height as the upstream ones.

    Fig.10 illustrates the distribution of the mean wall pressureand mean )skin friction coefficient, defined as Cf=τw/ (1 /2),where τwis the streamwise components of WSS.Compared with the inviscid shock solution, we find that the distribution ofin Fig.10(a) is significantly impacted by the viscous effect,which spreads the wall pressure jump over a much larger area with a smoother gradient.Furthermore, the skin friction reported in Fig.10(b) drops as the wall pressure increases after entering the interaction region,which indicates that the boundary layer decelerates under the effects of the adverse pressure gradient.Specifically, for α = 28°, the skin friction quickly decreases to 38% of its incoming level before the incident shock point and then recovers quickly to its nominal value after the incident shock point.No negative skin friction is observed in the interaction region,indicating that the boundary layer is statistically fully attached across the interaction.However, for α = 33.2°, the skin friction drops sharply but becomes negative before rising to a minimal positive value.After this brief rise, the skin friction reduces to a negative value before increasing to a positive one.The distribution of the mean skin friction predicts two distinct separation bubbles,which is consistent with the observations of Fang et al.36

    Clearly, the streamwise distributions of the skin friction in the downstream region for both cases are distinct.To intuitively evaluate the streamwise evolution process of the skin friction, a theoretical estimate proposed by White55for the undisturbed TBL is plotted in Fig.10(b) with the dashed line,

    where x0is the ‘‘virtual” origin of the turbulent part of the boundary layer, and

    Pirozzoli et al.35showed good agreement between the above formula and its DNS results, thus providing a reference skin friction value for the downstream undisturbed TBL.For α = 28°, the dropped skin friction recovers sharply to a value of approximately 2.3×10-3,and then decreases gradually farther downstream.In contrast, for α = 33.2°, the entire recovery process for the skin friction is gentle.The skin friction experiences a consistent increase downstream of the interaction region and does not rise to the same level as the upstream TBL until x=7.4,which then presents a continuously increasing trend.In summary, the recovery process for the skin friction in the separated boundary layer is relatively distinct from that in the attached boundary layer.In addition,the skin friction of the downstream disturbed TBL for both cases are more significant than that of the undisturbed TBL calculated from Eq.(3), although under the blockage effects caused by the incident shock.

    3.2.Mean skin friction decomposition

    Three contributions in Eq.(5) are: Cf,V, which represents the direct viscous dissipation effect; Cf,T, which represents the power spent for turbulent kinetic energy production; Cf,G, which represents the spatial growth of the flow, including the flow convection, streamwise heterogeneity, and pressure gradient.

    The decomposition results for the fully developed TBL at the reference station of x = - 4 are given in Table 1, compared with those of compressible zero-pressure-gradient flat plate TBL from Fan et al.56at Ma = 2 and Reτ= 580, and Tong et al.11at Ma = 2.25 and Reτ= 769.All three components positively contribute to skin friction generation,which is closely consistent with the previous studies.The TKEproduction Cf,Tand direct dissipation Cf,Vare the predominant components, which yield up to 46.6% and 40.8% of Cf,respectively.Contribution coming from the spatial growth Cf,Gcan be neglected, only reaching 12.6% of Cf.Figs.11 and 12 show the decomposed components of the skin friction at different streamwise locations for both cases, respectively.The relative error in this research is calculated as 100 ×(Cf,V+ Cf,T+ Cf,G- Cf) /Cf%.The errors at all locations for both cases are negligibly small and generally confined to±3%,reflecting the accuracy of the decomposition method in STBLI flows.

    The mean skin friction decomposition results in the interaction region for α = 33.2° are plotted in Fig.11.Cf,Vand Cf,Tare positive contributions,while Cf,Gnegatively contributes to the skin friction generation in the interaction region, which closely resembles the decomposition results in adversepressure-gradient TBL as developed on flat-plates and airfoils.57In general,the absolute values of Cf,Tand Cf,Gincrease significantly in the separation region and become far larger than Cf, whereas a rapid decrease is observed in the reattachment region.The trend for Cf,Vexhibits an opposite behavior.

    Compared with the upstream TBL, a completely different scenario emerges in the separation region characterized by small negative skin frictions, where large positive TKEproduction and large negative spatial growth dominate.Specifically, the positive direct dissipation Cf,Vexperiences a sharp decrease,and its contribution is negligible, whereas the spatial growth Cf,Gbecomes negative and takes the leading role in the skin friction generation.These counteract the significantlyamplified positive TKE-production.For instance, at x =- 0.75, the absolute values of Cf,Tand Cf,Gare much greater than Cfby one order of magnitude.

    Table 1 Decomposed skin friction components for fully developed turbulent boundary layer.

    Fig.11 Decomposed skin friction components in interaction region for α = 33.2°.

    Fig.12 Decomposed skin friction components in interaction region for α = 28°.

    In the reattachment region, the predominance of Cf,Gis gradually overtaken by the increased positive Cf,Vand decreased Cf,T, which causes a consistently increased positive Cfthat appears from x = 1.69 to 11.81.Two observations are made regarding this behavior.First,despite the magnitudes of both Cf,Tand Cf,Gexperiencing a consistent decrease, the attenuation of Cf,Gis greater than that of Cf,T, whereby the positive Cf,Tplays a dominant role in the skin friction generation.As seen in Fig.11,at x =1.69,Cf,Tand Cf,Ghave comparable magnitudes only in the initial part of the reattachment region.However, at x = 11.81 and downstream of the reattachment region, Cf,Gis dramatically reduced with only 7.9% remaining, while Cf,Tdecreases to 46.4% of its value at x=1.69.Second,the positive contribution of Cf,Vbecomes apparent.At x = 11.81, the magnitude of Cf,Vis approximately 2.5 times larger than that for Cf,G, which contributes to 27.6% of Cf.

    Fig.12 shows similar trends for α = 28° at the direct viscous dissipation Cf,V, TKE-production Cf,T, and spatial growth Cf,G.Specifically, before x = - 0.33 and around the nominal shock incident point, Cf,Tincreases sharply and then consistently decreases in the downstream region, whereas Cf,Vand Cf,Ghold opposite behavior.However, the contributions of the three terms for the Cfgeneration change completely.As mentioned above,Cf,Vis negligible in the separation region and gradually takes a predominant role in the farther downstream reattachment region.For α = 28°, Cf,Vis comparable with Cf,Tat x = 0.3.It is the initial part of the recovering boundary layer and contributes to 33.4% of Cf, much greater than that in Fig.11.

    Moreover, Cf,Tand Cf,Gbecome relatively less predominant.For instance, the maximum value of Cf,Tin Fig.12 appears at x = - 0.33, which is only 1.6 times that of the upstream TBL and much less than 4.8 times that for α = 33.2°.Although Cf,Gis also gradually increased downstream of the interaction region, it has a small positive contribution to Cfgeneration for α = 28°.In general, the dramatically increased Cf,Vand relatively high Cf,Tlead to the rapid recovery of Cf, despite Cf,Ghaving a negative contribution.

    3.3.EMD analysis of TKE-production contributions

    Empirical mode decomposition is an adaptive mode decomposition method first proposed by Huang et al.58EMD directly extracts intrinsic mode functions based on the signal’s instantaneous features, unlike Fourier analyses that require a set of basis functions in advance.Thus, EMD is a data-driven and posteriori method for data analysis.Readers can refer to the work of Huang et al.58and Cheng et al.59for more details about EMD procedures.It is inherently suitable for nonstationary and non-linear processes like those in turbulent flows.The advantages in processing multiscale signals have led to EMD being applied extensively in fluid mechanics.10,11,59,60Recently, attached eddies with specific length scales in turbulent channel flows were identified by Cheng et al.59using bidimensional EMD.The contributions of these structures to skin friction generation were quantitatively evaluated.

    Here, the scale decomposition of the turbulent motion is achieved by employing bidimensional EMD technology to investigate the contribution of the turbulent fluctuations at different length scales in the interaction region to Cfgeneration in STBLI flows.The velocity fluctuations u” and v” in the y-z plane are decomposed into four EMD modes, which is different from that of Cheng et al.59who decomposed velocity fluctuations in the homogeneous directions(x-z)of channel flows.Furthermore,different from Tong et al.,24to have a more comprehensive understanding of the skin friction generation characteristics in the STBLI flows, the situation in the reattachment region is analyzed, and the separation region is included.In addition, an extra case without flow separation is conducted to evaluate the effect of the interaction strength on the results.characterizes the Large Structure Motions (LSM) and Very Large Structure Motions (VLSM) in the wall-bounded turbulent flows.Overall, the decomposed results in our study are consistent with the results of Cheng et al.59

    Fig.13 Spanwise premultiplied energy spectra of streamwise velocity fluctuations (a) and wall-normal velocity fluctuati-ons (b) at reference station.The black lines represent the original DNS case.The green, purple, orange and pink lines represent Modes 1 to 4,respectively.

    Fig.14 Spanwise premultiplied energy spectra at reference station:(a)-(d)spectra of decomposed streamw-ise velocity fluctuations in Modes 1 to 4, respectively; (e)-(h) spectra of decomposed wall-normal ve-locity fluctuations in Modes 1 to 4, respectively.

    To further assess the contributions of turbulent fluctuations to different length scales on the mean skin friction generation,the corresponding Reynolds Shear Stress (RSS) is calculated using u′′dand v′′d.As analyzed above, velocity fluctuations are separated into four EMD modes.Full velocity fluctuations have the relations:

    where the number subscript represents velocity fluctuations in the corresponding EMD mode.Hence, the Reynolds shear stress ~-u′′v′′can be decomposed into 16 components as four primary parts and 12 cross parts representing interactions among different scales.These 16 parts of the RSS, the reconstructed RSS, and the original RSS are plotted in Fig.15.The profiles of the original and reconstructed RSS collapse well with each other, which indicates the correctness of our EMD analysis.In addition,we find that the four primary parts dominate at different heights, whereas the 12 cross parts are almost negligible.

    Then, the Cf,Tterm related to the TKE-production can be recalculated based on the decomposed RSS components,which is associated with different-scale motions as

    where the subscripts i = 1, 2, ???, 16 have the following relationship with the mode number m and n: i = 4 m + n - 4.The 16 parts of the contribution computed by 100×Cf,Ti/Cf,T%are shown in Fig.16.The primary parts Cf,T1,Cf,T6,and Cf,T16provide the top three contributions to the Cf,Tterm,whereas the primary part Cf,T11is less significant.Moreover,the cross parts represent interactions between two adjacent modes, i.e., Cf,T2, Cf,T5, Cf,T7, Cf,T10, Cf,T12, and Cf,T15, which can provide limited contributions.However, the cross parts represent interactions between two nonadjacent modes, i.e.,Cf,T3, Cf,T4, Cf,T8, Cf,T9, Cf,T13, and Cf,T14, which can be neglected.

    Fig.16 Contributions of Reynolds shear stress shown in Fig.15 to Cf,T.

    In general,the majority of Cf,Thave contributions from the parts associated with the inner structures, including the primary parts Cf,T1and Cf,T6, as well as the cross parts Cf,T2and Cf,T5, which even includes Cf,T7and Cf,T10.These give a total contribution of 66.2%.Thus, we believe Cf,Tis dominated by near-wall small-scale turbulent motion upstream of the fully developed TBL.Additionally, secondary contributions are from Cf,T11, Cf,T12, Cf,T15, and Cf,T16at 32.6%.It indicates that large-scale motion at the outer layer and interactions between them with motion within the boundary layer are significant in high Reynolds number flows.Similarly,Li et al.26also found a secondary peak at the premultiplied integrands function of Cf,Twhen the Reynolds number is up to Reτ=800 when considering the impact of the Reynolds number on the mean skin friction in compressible channel flows.

    The full and decomposed spectra at the separation region of x = - 0.75 for α = 33.2° are respectively reported in Figs.17 and 18.The length scales are normalized by the viscous length scales calculated at the reference station in the following energy spectrum analysis.Two noteworthy phenomena are observed.First, as shown in Fig.17(a), the near-wall primary peak of the energy spectra for u", which is observed at the upstream TBL, is absent at the separation region.It indicates that the near-wall streak structures are destroyed in the separation region, consistent with the instantaneous streamwise velocity distribution reported in Fig.6(b).Second, most energy is concentrated at higher locations with a peak at y+≈220 (y ≈0.3, in the outer scaling) for a larger spanwise length scale.The outer peak of the energy spectra reflects the Fig.18 Spanwise premultiplied energy spectra at x = - 0.75 for α = 33.2°: (a)-(d) spectra of decomposed streamwise velocity fluctuationsin Modes 1 to 4, respectively; (e)-(h) spectra of decomposed wall-no-rmal velocity fluctuationsin Modes 1 to 4,respectively.characteristics of the intense velocity fluctuations within the detached shear layer.

    Fig.17 Spanwise premultiplied energy spectra of streamwise velocity fluctuations (a) and wall-normal velocity fluctuations (b) at x= -0.75 for α=33.2°.The black lines represent the original DNS case,and the green,purple,orange,and pink lines represent Modes 1 to 4, respectively.

    As shown in Fig.19, the RSS amplitude at the separation region for α=33.2°is amplified approximately two times that of the upstream TBL.The turbulence amplification phenomenon in the interaction region is one of the primary characteristics of STBLI flow, which is generally observed via numerical simulations and experimental research.36,53,61,62Research on the characteristics of RBL by Tong et al.24,25indicates that turbulent motion is associated with the TKEproduction term Cf,Tand plays an important role in Cfgeneration in STBLI flows.Section 3.2 shows that the magnified turbulence contributes greatly to the generation of Cfthrough the mean skin friction decomposition analysis, regardless of whether flow separation occurs.In Fig.19, the original RSS is dominated by the four primary parts, while the 12 cross parts are relatively insignificant.The peak of the four primary parts is approximately at the same height, consistent with the observation of energy spectra.Notably, the fourth primary part occupies most of the original RSS at a greater height far from the wall surface.

    Fig.19 Profile of Reynolds shear stress as calculated from decomposed velocity fluctuations at x = -0.75 for α=33.2°,(a)for primary terms and (b) for cross terms.

    The 16 parts that contribute to the Cf,Tterm as calculated from the decomposed RSS are plotted in Fig.20.The parts associated with small-scale structures are dominated at the upstream TBL, involving the primary parts Cf,T1and Cf,T6,as well as the cross parts Cf,T2and Cf,T5,but are not prominent in the separation region with a total contribution of only about 19.7%.In contrast,the parts associated with large-scale structures dominate.Even the single part Cf,T16makes up 38%,and the total contributions of the parts associated with large-scale motion,including the primary parts Cf,T11and Cf,T16and cross parts Cf,T12and Cf,T15,reach 65.4%.Therefore,we believe that the intense fluctuations within the detached shear layer,generally considered with large spanwise length scales,63dominate the Cf,Tterm contributions in the separation region.The contributions of small-scale motions within the separation bubble are relatively limited.

    We show the spanwise premultiplied energy spectra of the original and EMD decomposed velocity fluctuations at the reattachment region of x = 8.43 for α = 33.2° in Figs.21 and 22, respectively.The spectra for u” and v” are shown as black lines in Fig.21 with both primary peaks located at a higher position of y+≈275.The secondary peaks appear in the near-wall region, which resembles the primary peaks at the upstream TBL and suggests they are regenerated wallbounded vortical structures.However, the difference between the downstream RBL and upstream TBL is that the energy is concentrated primarily at the outer layer with a relatively large spanwise length scale instead of at the near-wall region with a small scale.

    Fig.20 Contributions of Reynolds shear stress shown in Fig.19 to Cf,T.

    The spectra of the velocity fluctuations in the EMD modes are given in Fig.22.For the streamwise velocity fluctuations,the characteristics (including the peak height and most energy-carried wavelength) of the first two modes at the reattachment region are similar to those at the TBL.We consider that the first two modes still represent the near-wall streak structures shown in Fig.6(b),but these structures are regenerated and closer to the wall with smaller length scales.Modes 3 and 4 are more complex.As analyzed previously, Mode 3 at the upstream TBL represents the self-similar attached vortices at the log layer,and we believe that this may be available at the reattachment region.However,the higher energy-carried spectra span crosses a large region in the wall-normal direction and extends up to y+≈260.A similar situation is found in mode 4.Thus, we believe that the streamwise velocity fluctuations captured by Modes 3 and 4 are not only regenerated attached vortices but also include large-scale motion as convected from the upstream interaction region.Such energy-spanned phenomena are stronger for the wall-normal fluctuations shown in Figs.22(e)-(h).Additionally, the most energy-carried wavelengths for different modes are separated as shown in Fig.21.

    The decomposition of the RSS and their contributions to Cf,Tterm at x = 8.43 are shown in Figs.23 and 24, respectively.The RSS reaches a peak value at a relatively high location of y+= 300 (y = 0.4 in the outer scaling) at the reattachment region.It agrees well with the primary peak of the energy spectra in the outer layer shown in Fig.21, which further supports the appearance of large-scale vortices in the outer layer.Moreover, the fourth primary part occupies most of the original distribution, which results in the single part Cf,T16contributing to 53.1%of Cf,T.However,the total contributions of the primary parts Cf,T1and Cf,T6,and cross parts Cf,T2and Cf,T5,which represent the near-wall small-scale structures,only occupy 18.7%.In summary,consistent with studies of the RBL conducted by Tong et al.,24motion with large-scale length scales in the outer plays the primary role in contributing to the Cf,Tterm at the reattachment region,while the contributions of the near-wall small-scale structures are limited.It contrasts with the upstream TBL, where contributions of the motion with small scales in the near-wall region dominate.

    Fig.21 Spanwise premultiplied energy spectra of streamwise velocity fluctuations (a), and wall-normal velocity fluctuations (b) at x =8.43 for α=33.2°.The black lines represent the original DNS case,and the green,purple,orange,and pink lines represent Modes 1 to 4, respectively.

    Fig.22 Spanwise premultiplied energy spectra at x =8.43 for α=33.2°:(a)-(d)spectra of decomposed streamwise velocity fluctuationsin Modes 1 to 4, respectively; (e)-(h) spectra of decomposed wall-normal velocity fluctuations in Modes 1 to 4, respectively.

    Finally, we consider α = 28° at the same station as analyzed in the reattachment region to compare the effects of flow separation on the skin friction generation.The spanwise premultiplied energy spectra for the full and decomposed velocity fluctuations beneath the boundary layer downstream of the incident shock are shown in Figs.25 and 26,respectively.Note that there are 19 contour levels of the spectra shown in Figs.14,18, 22, and 26 ranging from 10% to 100% of their maximum values.The RSS profile, including the original and decomposed components, and its contribution to the Cf,Tterm are plotted in Figs.27 and 28, respectively.

    As shown in Figs.25–28, the characteristics at this station are similar to those at the upstream TBL.However, there are two points worth noting.First, the outlines of the plotted spectra are similar to those at the upstream TBL, as shown in Fig.25.For streamwise velocity fluctuations,the near-wall primary peak is nearly unaffected, while the secondary peak is stronger and located at a much higher position with a larger length scale.It results in the corresponding changes in the fourth mode.It is plausible that the outer large-scale structures tend to be influenced by the incident shock compared with the inner small-scale turbulence.Second, the maximum value of the RSS is slightly less than that at the upstream TBL, but its trend remains flat across a large region between y+≈30–160.Cf,T16gives the greatest contribution to the Cf,Tterm at up to 26.4% at the downstream boundary layer, which is different from the upstream TBL where Cf,T1provides the greatest contribution at up to 23.2%.

    Fig.23 Profile of Reynolds shear stress as calculated from decomposed velocity fluctuationsat x=8.43 for α=33.2°,(a)for primary terms and (b) for cross terms.

    Fig.24 Contributions of Reynolds shear stress shown in Fig.23 to Cf,T.

    Contributions from near-wall small-scale structures,as represented by the sum of parts Cf,T1, Cf,T2, Cf,T5, and Cf,T6,reduce from 56.1%at the upstream TBL to 44%at the downstream boundary layer.Contributions from the large-scale structures in the outer layer, denoted by the sum of parts Cf,T11, Cf,T12, Cf,T15, and Cf,T16, increase from 32.5% at the upstream TBL to 44.2%.Cf,Tis no longer contributed dominantly by near-wall structures but jointly by the inner and outer structures, which is different from the situation in the reattachment region as analyzed in the previous section.However, we recall that the Cf,Tterm increases by 13.3%at x = 8.43, as analyzed in Section 3.The absolute contribution of the inner small-scale structures to Cfremains nearly constant, but the absolute contribution from the outer largescale structures is greatly increased.In general, the contributions of the inner and outer structures are comparable to α = 28°, which differs from α = 33.2° at which the outer motions dominate Cfgeneration.

    The decomposed results at the same station downstream of the interaction region differ for both cases.In the separated case,large-scale structures are generated within the interaction region through various instabilities, including Kelvin-Helmholtz and centrifugal instabilities, which contribute greatly to the Cf,Tterm generation.Such large-scale structures can rarely be observed in the attached case,as shown in Fig.9(b), and the contributions of the related term are dramatically reduced.In addition, both the Cf,Tterm and contributions of the outer large-scale motions to the Cf,Tterm at the downstream boundary layer are greater than those at the upstream TBL,which is attributed to the effects of turbulence amplification.Thus,the decomposition results between the two cases at the same downstream station are distinct.

    Fig.25 Spanwise premultiplied energy spectra of streamwise velocity fluctuations(a),wall-normal velocity fluctuations(b)at x =8.43 for α = 28°.The black lines represent the original DNS case, and the green, purple, orange, and pink lines represent Modes 1 to 4,respectively.

    Fig.26 Spanwise premultiplied energy spectra at x =8.43 for α=28°:(a)-(d)spectra of decomposed streamwise velocity fluctuationsin Modes 1 to 4, respectively; (e)-(h) spectra of decomposed wall-normal velocity fluctuationsin Modes 1 to 4, respectively.

    Fig.27 Profile of Reynolds shear stress calculated from decomposed velocity fluctuationsat x = 8.43 for α = 28°, (a) for primary terms and (b) for cross terms.

    Fig.28 Contributions of Reynolds shear stress shown in Fig.27 to Cf,T.

    4.Conclusions

    The mean skin friction evolution characteristics of supersonic shock wave/turbulent boundary layer interactions are investigated using data from the DNS of incident shock waves with two different incident angles of 33.2°and 28°that interact with a spatially developing boundary layer at Ma∞= 2.25 and Reτ= 728.The concluding remarks are given below.

    (1) The ensemble-averaged results reveal that the distributions of Cfthroughout the interaction region between the two shock angle cases are distinct.For α = 33.2°,Cfgradually recovers at the reattachment region and requires at least 7.4δrefto retain the same level as the upstream TBL.For α=28°,the TBL remains attached when crossing the incident shock, and Cfquickly recovers within 1δref.

    (2) The mean skin friction Cfin STBLI flows is decomposed into a viscous dissipation term Cf,V, TKE-production term Cf,T, and flow spatial growth term Cf,Gusing the decomposition method proposed by recent research.The decomposition results indicate that the contributions of Cf,Vare limited in the case of flow separation,but it still plays an important role when the boundary layer remains attached.The contributions of Cf,Tare strengthened in both cases.The maximal magnification factor for α = 33.2° is 4.75 in the separation region and is about 2–3 in the reattachment region.The maximal magnification factor for α = 28° is approximately 1.6 near the shock incident point and 1.1–1.4 downstream of the interaction region.To counter the dramatically enhanced Cf,Tterm, the absolute value of Cf,Gis increased in the case of flow separation, which can negatively contribute to Cfgeneration.

    (3) Bidimensional EMD technology is employed to separate turbulent motion within the interaction region into four modes, characterized by specific spanwise length scales,to investigate further the contributions of different scale motions to the Cf,Tterm in STBLI flows.For the fully developed TBL, the contributions of inner small-scale motion are dominant, while those of large-scale motion are much weaker.The energy-carried structures with a large spanwise length are observed in the detached shear layer and at the outer layer of the RBL, and their contributions occupy the vast majority of the Cf,Tterm generation.Those from the inner small-scale structures can be neglected.However, such large-scale structures are absent in the attached case, dramatically reducing their contributions.Even so, the outer and inner motions are equally important, possibly because of the turbulence amplification phenomena.

    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.

    91麻豆av在线| 婷婷色麻豆天堂久久| 欧美日韩综合久久久久久| 亚洲男人天堂网一区| 国产精品久久久人人做人人爽| 国产真人三级小视频在线观看| 国产亚洲av片在线观看秒播厂| 老司机亚洲免费影院| 欧美日韩精品网址| 男人爽女人下面视频在线观看| 老司机深夜福利视频在线观看 | 男的添女的下面高潮视频| 久久人妻熟女aⅴ| 99国产精品一区二区蜜桃av | 黄色片一级片一级黄色片| 99热国产这里只有精品6| 少妇人妻 视频| 99精品久久久久人妻精品| 免费观看a级毛片全部| 亚洲人成电影免费在线| 一区二区av电影网| 精品久久蜜臀av无| 欧美黄色淫秽网站| 国产精品人妻久久久影院| 高清av免费在线| 一边亲一边摸免费视频| 亚洲国产中文字幕在线视频| 国产片内射在线| av在线播放精品| 99久久综合免费| 国产成人a∨麻豆精品| 欧美国产精品va在线观看不卡| 日韩中文字幕欧美一区二区 | 大香蕉久久成人网| 亚洲美女黄色视频免费看| 99国产精品一区二区蜜桃av | 中文欧美无线码| 一区在线观看完整版| 麻豆国产av国片精品| 亚洲av成人精品一二三区| 国产成人一区二区三区免费视频网站 | 免费看不卡的av| 黄片小视频在线播放| 天天躁狠狠躁夜夜躁狠狠躁| 成人午夜精彩视频在线观看| 777米奇影视久久| 男女边吃奶边做爰视频| 精品久久久久久电影网| 日韩制服骚丝袜av| 亚洲伊人久久精品综合| 两人在一起打扑克的视频| 婷婷丁香在线五月| 在线观看国产h片| 啦啦啦在线免费观看视频4| 欧美人与性动交α欧美软件| 午夜两性在线视频| 亚洲专区中文字幕在线| av在线app专区| 久久亚洲精品不卡| kizo精华| 岛国毛片在线播放| 欧美日韩综合久久久久久| 男的添女的下面高潮视频| 别揉我奶头~嗯~啊~动态视频 | 日本wwww免费看| 色综合欧美亚洲国产小说| 黄片播放在线免费| 天天躁日日躁夜夜躁夜夜| 狂野欧美激情性bbbbbb| 成年av动漫网址| 欧美97在线视频| 国产欧美日韩精品亚洲av| 国产伦人伦偷精品视频| 一级片免费观看大全| 青草久久国产| 日韩大片免费观看网站| √禁漫天堂资源中文www| 国产成人免费无遮挡视频| 欧美人与性动交α欧美软件| 韩国高清视频一区二区三区| 日本91视频免费播放| 欧美激情 高清一区二区三区| 啦啦啦 在线观看视频| 视频区欧美日本亚洲| 51午夜福利影视在线观看| 亚洲激情五月婷婷啪啪| 久久九九热精品免费| 99久久精品国产亚洲精品| 国产精品欧美亚洲77777| 久久青草综合色| 一边摸一边抽搐一进一出视频| 国产亚洲av高清不卡| 国产日韩欧美亚洲二区| 亚洲,欧美精品.| 99精国产麻豆久久婷婷| 男的添女的下面高潮视频| 又大又黄又爽视频免费| 亚洲欧洲日产国产| 精品一区二区三区av网在线观看 | 久久鲁丝午夜福利片| 久久久国产精品麻豆| 精品福利永久在线观看| 国产国语露脸激情在线看| 国产人伦9x9x在线观看| 日韩免费高清中文字幕av| 日韩中文字幕欧美一区二区 | 最新在线观看一区二区三区 | 秋霞在线观看毛片| 亚洲图色成人| 男女无遮挡免费网站观看| 久久久久精品人妻al黑| 日韩av在线免费看完整版不卡| 人人妻人人澡人人看| av一本久久久久| 欧美乱码精品一区二区三区| 欧美 日韩 精品 国产| 欧美黑人精品巨大| 成人亚洲欧美一区二区av| 免费人妻精品一区二区三区视频| e午夜精品久久久久久久| 黑人欧美特级aaaaaa片| 高清欧美精品videossex| 老司机靠b影院| 久久亚洲精品不卡| 久久国产精品大桥未久av| 超色免费av| 极品人妻少妇av视频| 久久ye,这里只有精品| 2021少妇久久久久久久久久久| 妹子高潮喷水视频| 老熟女久久久| 久久久欧美国产精品| 成人亚洲欧美一区二区av| 欧美日韩成人在线一区二区| av片东京热男人的天堂| 日日爽夜夜爽网站| 精品少妇一区二区三区视频日本电影| 亚洲欧美日韩另类电影网站| 嫩草影视91久久| 午夜福利一区二区在线看| 欧美日韩综合久久久久久| 亚洲一区二区三区欧美精品| 精品国产一区二区久久| 国产日韩一区二区三区精品不卡| 美女国产高潮福利片在线看| 好男人电影高清在线观看| 午夜福利,免费看| 久久99一区二区三区| 亚洲熟女精品中文字幕| 一级毛片 在线播放| 男女高潮啪啪啪动态图| 亚洲精品av麻豆狂野| a级毛片在线看网站| 一级片免费观看大全| 一本色道久久久久久精品综合| 亚洲色图综合在线观看| 欧美人与性动交α欧美软件| 老汉色∧v一级毛片| h视频一区二区三区| 国产有黄有色有爽视频| 波多野结衣一区麻豆| 亚洲欧洲日产国产| 久久精品亚洲熟妇少妇任你| 国产爽快片一区二区三区| 亚洲,欧美精品.| 天天躁夜夜躁狠狠久久av| 成人国产av品久久久| 老司机午夜十八禁免费视频| 国产精品一区二区在线不卡| 咕卡用的链子| 飞空精品影院首页| 久久亚洲国产成人精品v| 大型av网站在线播放| 国产成人av教育| 黄片播放在线免费| 777久久人妻少妇嫩草av网站| videos熟女内射| 国产一区二区三区综合在线观看| 一级,二级,三级黄色视频| 中文字幕色久视频| 丰满少妇做爰视频| 丝袜在线中文字幕| 国产不卡av网站在线观看| 亚洲一码二码三码区别大吗| 午夜老司机福利片| 色视频在线一区二区三区| 老司机在亚洲福利影院| 精品人妻一区二区三区麻豆| 欧美人与善性xxx| 亚洲一区二区三区欧美精品| 一边亲一边摸免费视频| 老司机深夜福利视频在线观看 | 国产亚洲av片在线观看秒播厂| 老汉色av国产亚洲站长工具| 国产一区二区三区av在线| 亚洲精品国产色婷婷电影| 无限看片的www在线观看| 女人高潮潮喷娇喘18禁视频| 国产又色又爽无遮挡免| 香蕉丝袜av| 久久精品久久久久久噜噜老黄| 一级毛片女人18水好多 | 亚洲欧美一区二区三区黑人| 久久这里只有精品19| 亚洲精品av麻豆狂野| 热re99久久精品国产66热6| a级片在线免费高清观看视频| 一级片免费观看大全| 日韩av免费高清视频| 亚洲自偷自拍图片 自拍| 极品少妇高潮喷水抽搐| 老司机影院毛片| 99精品久久久久人妻精品| 国精品久久久久久国模美| 国产一级毛片在线| 高清视频免费观看一区二区| 亚洲精品久久久久久婷婷小说| 亚洲国产中文字幕在线视频| 热99国产精品久久久久久7| 97在线人人人人妻| 亚洲av日韩精品久久久久久密 | 99国产精品99久久久久| 精品国产超薄肉色丝袜足j| 亚洲第一青青草原| 欧美亚洲 丝袜 人妻 在线| 久久影院123| 一级黄色大片毛片| 国产精品国产av在线观看| 国产精品 国内视频| svipshipincom国产片| 亚洲午夜精品一区,二区,三区| 新久久久久国产一级毛片| 午夜福利免费观看在线| 婷婷色麻豆天堂久久| 波多野结衣一区麻豆| 亚洲av电影在线观看一区二区三区| 久久热在线av| 一级黄片播放器| 男女午夜视频在线观看| 亚洲九九香蕉| 99热网站在线观看| 亚洲av片天天在线观看| 十八禁高潮呻吟视频| 黑丝袜美女国产一区| 七月丁香在线播放| 一区二区三区乱码不卡18| 国产欧美日韩一区二区三区在线| 亚洲,欧美,日韩| 久久99精品国语久久久| 丰满迷人的少妇在线观看| 亚洲av综合色区一区| 18禁国产床啪视频网站| 久久精品成人免费网站| 国产精品一区二区在线不卡| 午夜日韩欧美国产| 另类精品久久| 精品少妇黑人巨大在线播放| 后天国语完整版免费观看| 国产精品一区二区免费欧美 | 久久国产精品影院| 男女无遮挡免费网站观看| 精品福利永久在线观看| 夫妻午夜视频| 麻豆乱淫一区二区| 熟女av电影| 十八禁人妻一区二区| 制服人妻中文乱码| 激情视频va一区二区三区| 午夜福利,免费看| 在线观看免费午夜福利视频| 51午夜福利影视在线观看| 电影成人av| 国产人伦9x9x在线观看| 黄网站色视频无遮挡免费观看| 免费看不卡的av| 香蕉国产在线看| av国产精品久久久久影院| 悠悠久久av| 国产精品av久久久久免费| 久久久久久亚洲精品国产蜜桃av| 日韩免费高清中文字幕av| 亚洲第一av免费看| 最近最新中文字幕大全免费视频 | 无遮挡黄片免费观看| 中文字幕人妻丝袜制服| 一本色道久久久久久精品综合| 亚洲国产最新在线播放| 国产亚洲av片在线观看秒播厂| 国产成人精品久久久久久| 亚洲精品国产一区二区精华液| 久久国产精品影院| 国产av一区二区精品久久| 亚洲精品久久久久久婷婷小说| 人人澡人人妻人| 国产精品国产三级国产专区5o| kizo精华| 国产精品成人在线| 黄色a级毛片大全视频| 少妇粗大呻吟视频| 最近中文字幕2019免费版| 久久久久国产一级毛片高清牌| 国产精品99久久99久久久不卡| 宅男免费午夜| 免费观看a级毛片全部| netflix在线观看网站| 丁香六月天网| 亚洲九九香蕉| 国产激情久久老熟女| 秋霞在线观看毛片| 黄色一级大片看看| 黄色片一级片一级黄色片| 好男人视频免费观看在线| 精品亚洲成国产av| 亚洲欧美一区二区三区黑人| 亚洲精品国产一区二区精华液| 一级黄片播放器| 国产欧美亚洲国产| 亚洲av男天堂| 下体分泌物呈黄色| 国产伦理片在线播放av一区| 后天国语完整版免费观看| 久久青草综合色| 黑人巨大精品欧美一区二区蜜桃| 天天添夜夜摸| 亚洲国产中文字幕在线视频| 久久精品久久久久久久性| 啦啦啦 在线观看视频| 日本wwww免费看| 国产老妇伦熟女老妇高清| 熟女av电影| 中文字幕人妻熟女乱码| 亚洲国产最新在线播放| 亚洲第一av免费看| 1024视频免费在线观看| 人人澡人人妻人| 新久久久久国产一级毛片| 亚洲av欧美aⅴ国产| 天堂8中文在线网| 久久精品国产a三级三级三级| 亚洲七黄色美女视频| 日韩av免费高清视频| 建设人人有责人人尽责人人享有的| 国产一区二区三区综合在线观看| 国产一区亚洲一区在线观看| 777久久人妻少妇嫩草av网站| 亚洲成av片中文字幕在线观看| 少妇人妻久久综合中文| 热re99久久国产66热| 国产真人三级小视频在线观看| 香蕉国产在线看| 美女福利国产在线| 一区福利在线观看| 成人黄色视频免费在线看| 中文字幕精品免费在线观看视频| 日韩电影二区| 两性夫妻黄色片| 日日摸夜夜添夜夜爱| 男人爽女人下面视频在线观看| 国产无遮挡羞羞视频在线观看| 美女午夜性视频免费| 精品免费久久久久久久清纯 | 久热这里只有精品99| 亚洲精品一二三| 你懂的网址亚洲精品在线观看| 精品国产一区二区久久| 中文字幕av电影在线播放| 欧美成人午夜精品| 少妇被粗大的猛进出69影院| 国产精品偷伦视频观看了| netflix在线观看网站| 日本五十路高清| 在线看a的网站| 色视频在线一区二区三区| 国产一区有黄有色的免费视频| 丝袜人妻中文字幕| 午夜91福利影院| 成年动漫av网址| 别揉我奶头~嗯~啊~动态视频 | 亚洲国产最新在线播放| 国产精品香港三级国产av潘金莲 | 亚洲欧美一区二区三区国产| 高清视频免费观看一区二区| 国产精品亚洲av一区麻豆| 麻豆乱淫一区二区| 丝袜美腿诱惑在线| 欧美日韩视频高清一区二区三区二| 悠悠久久av| 久久中文字幕一级| 国产片内射在线| 亚洲国产精品国产精品| 久久精品成人免费网站| 高清不卡的av网站| av欧美777| 国产极品粉嫩免费观看在线| 欧美精品亚洲一区二区| 欧美激情极品国产一区二区三区| √禁漫天堂资源中文www| 久久人人97超碰香蕉20202| 国产福利在线免费观看视频| 亚洲色图综合在线观看| 国产精品香港三级国产av潘金莲 | 久久精品久久久久久久性| 肉色欧美久久久久久久蜜桃| 我的亚洲天堂| 欧美av亚洲av综合av国产av| 成人手机av| 久久精品久久精品一区二区三区| 国产淫语在线视频| 菩萨蛮人人尽说江南好唐韦庄| 国产精品久久久久久精品古装| 最近最新中文字幕大全免费视频 | 看免费成人av毛片| 一级毛片 在线播放| 欧美日韩黄片免| 国产精品一区二区在线观看99| 成人影院久久| 天天躁狠狠躁夜夜躁狠狠躁| 午夜av观看不卡| www日本在线高清视频| 国产在线视频一区二区| 免费久久久久久久精品成人欧美视频| 纯流量卡能插随身wifi吗| 国产成人精品无人区| netflix在线观看网站| 深夜精品福利| 亚洲欧美成人综合另类久久久| 亚洲一区中文字幕在线| 国产精品一区二区在线观看99| 国产视频首页在线观看| 新久久久久国产一级毛片| 久9热在线精品视频| 侵犯人妻中文字幕一二三四区| 国产在线视频一区二区| 在线av久久热| 晚上一个人看的免费电影| 日韩av在线免费看完整版不卡| 亚洲av片天天在线观看| 国产精品久久久人人做人人爽| 首页视频小说图片口味搜索 | 亚洲精品久久久久久婷婷小说| 精品高清国产在线一区| a级片在线免费高清观看视频| 成人国产av品久久久| 欧美日韩亚洲综合一区二区三区_| 精品第一国产精品| 日韩精品免费视频一区二区三区| 少妇人妻久久综合中文| 女人被躁到高潮嗷嗷叫费观| 老鸭窝网址在线观看| 国产精品久久久久成人av| 女人久久www免费人成看片| 伦理电影免费视频| 欧美97在线视频| 十分钟在线观看高清视频www| 久久久久国产精品人妻一区二区| 亚洲免费av在线视频| 99香蕉大伊视频| 国产亚洲av高清不卡| 国产福利在线免费观看视频| 丰满少妇做爰视频| 啦啦啦 在线观看视频| 成人18禁高潮啪啪吃奶动态图| 国产在线观看jvid| 日本一区二区免费在线视频| 亚洲精品在线美女| 1024视频免费在线观看| 欧美国产精品va在线观看不卡| 久久女婷五月综合色啪小说| 久久久久久久精品精品| 国产一卡二卡三卡精品| 精品欧美一区二区三区在线| 久久精品国产a三级三级三级| 男人操女人黄网站| 欧美成人午夜精品| av网站免费在线观看视频| 欧美人与性动交α欧美软件| videosex国产| 久久久久久久精品精品| 欧美人与善性xxx| 久久午夜综合久久蜜桃| 亚洲精品久久成人aⅴ小说| 日韩av免费高清视频| 天天躁狠狠躁夜夜躁狠狠躁| 操美女的视频在线观看| 午夜免费成人在线视频| 99热国产这里只有精品6| 亚洲中文字幕日韩| 日韩av在线免费看完整版不卡| 老司机影院成人| 天天添夜夜摸| 日日夜夜操网爽| 亚洲精品美女久久av网站| 久久国产精品人妻蜜桃| 赤兔流量卡办理| 国产精品免费大片| 少妇精品久久久久久久| 精品国产一区二区久久| 亚洲国产成人一精品久久久| 国产一区二区三区av在线| a级毛片黄视频| 国产精品av久久久久免费| 男男h啪啪无遮挡| 久久精品人人爽人人爽视色| 午夜免费鲁丝| 热99久久久久精品小说推荐| 首页视频小说图片口味搜索 | 后天国语完整版免费观看| 热re99久久精品国产66热6| 18在线观看网站| 老司机在亚洲福利影院| 欧美精品一区二区大全| 免费日韩欧美在线观看| 男女免费视频国产| 老司机影院成人| 亚洲欧美一区二区三区黑人| 我要看黄色一级片免费的| 夫妻性生交免费视频一级片| 美女福利国产在线| 亚洲中文av在线| 高潮久久久久久久久久久不卡| 亚洲精品一二三| 久久午夜综合久久蜜桃| 啦啦啦中文免费视频观看日本| 久久久久视频综合| 精品亚洲乱码少妇综合久久| 中文字幕人妻熟女乱码| 婷婷色综合www| 日韩,欧美,国产一区二区三区| 国产高清videossex| 午夜福利一区二区在线看| 桃花免费在线播放| 欧美成人精品欧美一级黄| 国产精品成人在线| 啦啦啦啦在线视频资源| 丝袜人妻中文字幕| 建设人人有责人人尽责人人享有的| 免费高清在线观看视频在线观看| 亚洲五月色婷婷综合| 国产女主播在线喷水免费视频网站| 日本欧美国产在线视频| 久久女婷五月综合色啪小说| 久久 成人 亚洲| 日韩,欧美,国产一区二区三区| 亚洲欧美日韩高清在线视频 | 一边摸一边抽搐一进一出视频| 免费在线观看完整版高清| 国产精品免费大片| 一区二区日韩欧美中文字幕| 午夜激情久久久久久久| 久久久国产欧美日韩av| 操美女的视频在线观看| 精品亚洲成国产av| 日本vs欧美在线观看视频| 欧美+亚洲+日韩+国产| 国产国语露脸激情在线看| 欧美 日韩 精品 国产| 后天国语完整版免费观看| 亚洲精品国产区一区二| 亚洲精品自拍成人| 亚洲成人免费电影在线观看 | netflix在线观看网站| 国产三级黄色录像| 午夜91福利影院| 精品久久久久久久毛片微露脸 | 国产成人一区二区三区免费视频网站 | 中文字幕色久视频| 天天躁狠狠躁夜夜躁狠狠躁| 亚洲伊人久久精品综合| 精品亚洲乱码少妇综合久久| 日韩伦理黄色片| 悠悠久久av| 免费在线观看黄色视频的| 亚洲专区中文字幕在线| 久久精品久久久久久噜噜老黄| 亚洲 国产 在线| 国产成人一区二区三区免费视频网站 | 国产黄色视频一区二区在线观看| 欧美在线一区亚洲| 国产精品99久久99久久久不卡| 亚洲精品日韩在线中文字幕| 久久久久国产精品人妻一区二区| 久久天躁狠狠躁夜夜2o2o | 国产成人精品无人区| 老司机在亚洲福利影院| 天天躁夜夜躁狠狠久久av| 亚洲欧美精品综合一区二区三区| 亚洲人成电影观看| 亚洲七黄色美女视频| 91精品国产国语对白视频| 老司机影院毛片| 一本久久精品| 免费女性裸体啪啪无遮挡网站| 在线 av 中文字幕| 电影成人av| 黄色 视频免费看| 日本vs欧美在线观看视频| 国产日韩欧美亚洲二区| 自拍欧美九色日韩亚洲蝌蚪91| 精品一区二区三卡| 国产精品欧美亚洲77777| 国产不卡av网站在线观看| 国产人伦9x9x在线观看| 亚洲av美国av| 日日爽夜夜爽网站| 黑丝袜美女国产一区| 99国产精品一区二区三区| 亚洲黑人精品在线| 国产精品一区二区精品视频观看| 最新的欧美精品一区二区| 国产无遮挡羞羞视频在线观看| 日本一区二区免费在线视频| 丁香六月天网| 美女扒开内裤让男人捅视频| 久久天堂一区二区三区四区| 久久久久久久精品精品|