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

    Experimental and numerical investigation of threedimensional vortex structures of a pitching airfoil at a transitional Reynolds number

    2019-12-19 02:07:38ZhenyoLILihoFENGHmidRezKARBASIANJinjunWANGKyungChunKIM
    CHINESE JOURNAL OF AERONAUTICS 2019年10期

    Zhenyo LI, Liho FENG, Hmid Rez KARBASIAN, Jinjun WANG,Kyung Chun KIM,*

    a Fluid Mechanics Key Laboratory of Education Ministry, Beihang University, Beijing 100083, China

    b School of Mechanical Engineering, Pusan National University, Busan 46241, Republic of Korea

    KEYWORDS

    DDES;

    Dynamic stall;

    Pitching airfoil;

    PIV;

    Transitional Reynolds number

    Abstract This research examines the vortex behaviors and aerodynamic forces in dynamic stall phenomena at a transitional Reynolds number (Re=90000) using experimental and numerical approaches. Periodic sinusoidal pitching motion at two different reduced frequencies is used to achieve the dynamic stall of a NACA 0012 airfoil. Several leading edge vortices form and detach in the dynamic stall stage. The flow then quickly transitions to a full separation zone in the stall stage when the angle of attack starts to decrease. There is discrepancy between the phaseaveraged and instantaneous flow field in that the small flow structures increased with angle of attack, which is a characteristic of the flow field at the transitional Reynolds number. The interaction between the streamwise vortices in the three-dimensional numerical results and the leading edge vortex are the main contribution to the turbulent flow.In addition,the leading edge vortex that supplies vortex lift is more stable at higher reduced frequency,which decreases the lift fluctuation in the dynamic stall stage. The leading edge vortex at higher reduced frequency is strong enough to stabilize the flow, even when the airfoil is in the down-stroke phase.

    1. Introduction

    The blades of a horizontal-axis wind turbine need to change pitch angle all the time to work efficiently in the case of unsteady incoming flow. However, this motion causes strong unsteady aerodynamic load on the blades. Dynamic stall happens during the motion of blades which is the key reason for the unsteady aerodynamic load.1One of the main characteristics of dynamic stall is that the maximum lift of an airfoil undergoing unsteady motion is far beyond the maximum lift that the static airfoil can reach. Other characteristics of dynamic stall include the formation and shedding of vortices of different size.2

    The aerodynamic characteristics of dynamic stall delay the static stall at high angles of attack, which has attracted many researchers because it can play an important role in wind turbines, aerospace, and energy harvesting. McCroskey3, Carr4,and Ekaterinaris and Platzer5comprehensively reviewed dynamic stall using experimental and numerical techniques.The dynamic stall of a pitching airfoil can be divided into four stages: no stall, stall onset, light stall, and deep stall6. These four stages can also be divided into three stages according to the development of the aerodynamics loads:pre-stall,dynamic stall, and stall.

    The pre-stall stage is a combination of no-stall and stall onset since the aerodynamic loads follow quasi-static aerodynamic theory in these two stages. The development of a Leading Edge Vortex (LEV) and large peak force are the main characteristics of the dynamic stall stage. In the stall stage,the aerodynamic loads drop quickly, and the flow transitions to a separation zone, which is similar to static stall. The separation behavior of flow is also divides dynamic stall into light and deep dynamic stall.

    In light dynamic stall,the LEV strengthens itself in the upstroke phase and terminates in the beginning of the downstroke phase.7In deep dynamic stall, the LEV transitions to severe flow separation in the early up-stroke phase.8LEV is thought to be the key flow structure in dynamic stall. Thus,understanding the development features of LEV can help explain the aerodynamic characteristics of dynamic stall.

    There are many parameters that influence the flow in dynamic stall. The Reynolds number and the reduced frequency were two important non-dimensional parameters to describe the pitching motion of an airfoil.3The reduced frequency is defined as k=πfc/U∞and is introduced to indicate the ratio of the convective time scale and the forced oscillation time scale.Here the f is the motion frequency of airfoil,the c is the chord length of the airfoil and the U∞is the free-stream velocity. The Reynolds number is defined as Re=ρU∞c/μ and is a common parameter in fluid dynamics. Here the ρ is the fluid density and the μ is the free-stream viscosity. The combination of these two parameters leads to complicated rules for vortex evolution and force dynamics.

    The Reynolds number is often used as a reference to distinguish different kinds of flow. Aircraft, helicopter rotors, and turbo-machinery often operate at high Reynolds number beyond the order of 106.Small to medium wind turbines work at medium Reynolds numbers on the order of 105. Low Reynolds numbers(Re ~O(104))are the operating range of Micro Air Vehicles (MAVs). With the different phenomena of dynamic stall at different Reynolds numbers,the LEV characteristics also exhibit great differences. Numerous researchers have studied the characteristics of a pitching airfoil at different Reynolds numbers. Some recent representative studies are summarized in Table 19-27.

    The Reynolds number used in the research in Table 1 are from 103to 106.These studies indicate that the vortex behavior strongly depends on the Reynolds number. At low Reynolds number,the LEV structure is more concentrated,and the maximum vortex strength is higher.At high Reynolds number,the LEV structure transitions to the full separation zone rapidly.More small-scale flow structures can also be observed during the motion process because it is much easier for the vortex structures to lose stability and separate into small threedimensional vortex structures.

    Table 1 Non-dimensional parameters of pitching airfoils in previous studies.

    The flow field shows more complicated flow structures at transitional Reynolds numbers, which are around 1.0×105.Studying a flow field that mixes the characteristics of these two regimes at transitional Reynolds numbers is important to understand the development of the flow field and aerodynamic forces at different Reynolds numbers. In addition,investigation of the evolution of flow structures, particularly three-dimensional ones, helps in understanding the development of aerodynamic forces during pitching motion since it is closely related to the vortex structure. However, threedimensional approaches are rarely used, including 3D Particle Image Velocimetry (3D-PIV) and 3D Computational Fluid Dynamics(3D-CFD),although there are some 3D CFD works at a low Reynolds number of 2×104.15,16Although there are many 2D-PIV studies at low Reynolds numbers, there are few comprehensive PIV studies at higher Reynolds numbers due to the difficulty of measurement.23-25

    Thus, this study mainly focuses on the evolutions of threedimensional flow structures and their influences on the aerodynamic forces during pitching motion at a transitional Reynolds number. The Reynolds number is fixed at 9×104, and two typical reduced frequencies are chosen as parameters. Both 2D time-resolved PIV and 3D CFD approaches are applied to reveal the physics of dynamic stall for a pitching NACA0012 airfoil.

    2. Experimental and numerical methods

    2.1. Experimental method

    An experiment was conducted in a recirculation water tunnel at Beijing University of Aeronautics and Astronautics. The size of the test section was 1000 mm×1200 mm in the spanwise and vertical directions, and the free-stream turbulence intensity was less than 0.8%. A schematic of the experimental setup is shown in Fig. 1. The experimental model was a NACA0012 airfoil with a chord length of c=300 mm and spanwise length of 800 mm. The distance between the top of the airfoil and the water surface was 2 mm, so the effective aspect ratio was extended to 5.33.14The airfoil was manufactured using rigid plastic to reduce the weight and black spray paint to reduce the surface reflection under laser illumination.The model was mounted vertically in the test section and maintained perpendicular to the flow direction, as shown in Fig. 1(a).The distance between the pivot point and the leading edge was one-quarter of the chord length, which follows most of previous studies about dynamic stall. It has been found that moving the pivot point forward is beneficial for the formation of the leading-edge vortex.28,29Thus, this study will not concern the effect of the pivot point.

    Airfoil pitching motion was produced by a servo motor(YASKAWA SGM7J) and a Panasonic AFPXHC60T PLC controller. A schematic of the pitching airfoil is shown in Fig. 1(b). The airfoil oscillated around the pivot point in the x-y plane. Since the rotation velocity was much smaller than the free-stream velocity,the effective angle of attack was equal to the angle between the free-stream velocity and chord-line.Given a sinusoidal oscillation as an input, the effective angle of attack α for pure pitching motion can be represented as follows:

    The mean effective angle of attack α0and the pitch oscillation amplitude αmwere fixed at 10° and 15°, respectively. The free-stream velocity was fixed at U∞=300 mm/s,and the corresponding Reynolds number was 9×104. The motion frequency was f=0.0318, which corresponds to k=0.1. Thus,the deflection is negligible at the half span location where the velocity fields were measured.

    The velocity field was measured by a time-resolved 2D PIV technique.A dual-head ND:YAG laser with 500 mJ/pulse was used to generate a light sheet with a thickness of 1.5 mm.Hollow glass beads with a median diameter of 20 μm and density of 1.05 g/cm3were used as tracer particles. To ensure the spatial resolution in a large field of view,two double-exposure 12-bit cameras (IMPERX B2520) with a resolution of 2058×2456 pixels and 50-mm Nikon lenses were installed in the stream direction. The field of view of one camera was approximately 178 mm×148 mm in the stream direction and vertical direction, respectively. The sampling frequency was 5 Hz, which corresponds to 155 snapshots in one cycle.

    Fig. 1 Schematic diagram of experimental setup and pitching airfoil.

    The PIV images were analyzed using the Multi-pass Iterative Lucas-Kanade (MILK) algorithm.30The convergent velocities were obtained by a pyramid of 4 and iterations of 7. The size of the interrogation windows was set to 32×32 pixels with 75% overlap. The window’s spatial resolution was 0.78% of the chord length. Erroneous velocity vectors were detected and replaced by peak detection, a normalized median test, and local interpolation.31To increase the measurement accuracy in the measurement of near-wall region,background subtraction techniques were used in the raw PIV pictures to get rid of the negative influence of laser sheet heterogeneity in PIV calculation. Meanwhile, the PIV meshes which contained the airfoil model boundary was eliminated.Phase average process was also done to over 5500 recorded velocity fields corresponding to 38 cycles to suppress the random noise in the flow field. The phase of motion cycle was identified depend on the output signal of servo motor and the angle of airfoil which was got by image identification techniques. The zero phase corresponded to the start position of the airfoil.

    2.2. Numerical method

    For 3D flow simulations,the flow field is solved using Delayed Detached Eddy Simulations (DDES) and the Finite Volume Method(FVM).The flow inside of boundary layer is modeled using the Unsteady Reynolds-Averaged Navier-Stokes(URANS) approach with the Shear Stress Transport (SST)k-ω turbulence model. However, outside of the boundary layer, the flow is resolved by the Large Eddy Simulation(LES) approach. The CFD flow solver package ANSYS Fluent version 16.1 (ANSYS Fluent, Academic Research, V16.1 2015)was used for the present study.In the case of spatial discretization, the bounced central differencing and second-order upwind formulations were selected for the momentum and turbulence, respectively. Additionally, a second-order implicit scheme was selected for the temporal discretization, and the time step was set as Δt=15×10-5T for all cases to make sure that the Courant-Friedrichs-Lewy (CFL) number is less than 2. The ‘‘Pressure Implicit with Splitting of Operator”(PISO) algorithm was used for the pressure-velocity coupling.

    Fig. 2 Computational domain and detailed mesh topology.

    Fig. 2 shows the 3D computational domain, which was designed with C-mesh topology. The mesh structure is magnified in different views.The blade is surrounded by a boundary mesh that is normal to the surface of the airfoil. The mesh is dense and fine enough near the leading and trailing edges to capture the strong shear layer and mixing layer in these regions appropriately. The profile section of the straight blade is a NACA0012 airfoil with a chord length of 0.1 m and aspect ratio of 2.The airfoil is enveloped by a boundary mesh normal to the wall with y+lower than 1, which meets the accuracy requirement of DDES solver. All outer boundaries were kept at about 12c away from the blade to make the blockage effect negligible.20The total number of nodes for the final design was 2.9×106, and 500 nodes were placed along the foil and clustered near the leading and trailing edges.

    To design a 3D computational domain for DDES simulations,a 2D domain in an x-y system was designed and run with the same boundary conditions. The boundary layer thickness was then computed for when the flow is still attached to airfoil surface. In the DDES approach, the flow is modeled inside of the boundary layer using URANS,while outside of the boundary layer, the solver switches to the LES approach and solve the flow in the wake region. Therefore, the boundary layer thickness allows us to design a mesh with appropriate cell dimensions to make sure that the flow is solved by the URANS solver inside the boundary layer. Next, the 2D computational domain is extruded uniformly in the z-direction with a spatial step of Δz=0.01c.

    For the boundary conditions, constant free stream velocity was applied upstream of the domain, and the upper, lower,and downstream boundaries were set to atmospheric pressure.The blade surface has no-slip conditions, and the side boundaries were set to periodic boundary conditions. A reduced frequency of k=0.1, Reynolds number of 1.35×105, and effective angle of attack of α=10°+15°sin(2πft) were chosen to verify CFD solver. Reduced frequencies of 0.1 and 0.2, a Reynolds number of 9×104, and the same effective angle of attack were used later in the study.

    2.3. Validation of CFD method

    Different computational domains were considered to validate the computational method.Four computational domains with different characteristics are provided in Table 2.Nzis the number of nodes in the spanwise direction. The number of nodes on the airfoil’s surface and the total number of nodes in the whole domain are also provided. In the case of H1, H2, and H3, the number of nodes in the spanwise direction was Nz=160,which was increased to Nz=200 for H4.In the case of the coarse grid(H1),the number of nodes on the surface of the airfoil was Nwall=200, while for the rest of the cases,Nwall=250. Eventually, H3 was chosen as an appropriate computational domain based on the average error computed for each computational domain with reference to H4. As shown in Table 2, the average error for H1 is very high(0.176), while the values for H2 and H3 are 0.0302 and 0.0196, respectively.

    The numerical results of lift and drag coefficients were compared with experimental results.21The lift and drag were normalized by chord length,fluid density and free-stream velocity.The experimental parameters are described above.All aerodynamic loads were computed based on integration over the entire surface of the airfoil. Fig. 3 which shows the development of lift and drag coefficients with pitching angle of attack is provided to investigate the accuracy of the numerical solver.As shown in Fig. 3(a), the lift coefficient increases in the prestall region with increasing angle of attack, and the numerical results are in good agreement with the experimental results of Lee and Gerontakos.21In the down-stroke phase, there are greater discrepancies between the experimental and numerical results due to the complexity of the flow in the post-stall region.

    Table 2 Comparison of different computational domains with specific characteristics.

    Fig. 3 Comparison of lift and drag coefficients for present numerical solver and experimental results (Lee and Gerontakos21) at Re=1.35×105, k=0.1, and α=10°+15°sin(2πft).

    In Fig. 3(b), the drag coefficient from the simulation is in good agreement with the experimental results. As the angle of attack increases, the numerical results follow the trend of the experimental results.The discrepancies between the numerical and experimental results might be due to the measurement methods since post processing of the pressure tap results rather than direct measurement was used in experiment. Additionally, the complex mechanism of the flow separation and the development of vortices make the numerical results more unstable in the post-stall region,which is still being considered by other researchers as well.32-34Overall, the numerical method is reliable, although experimental results are still needed as a reference.

    3. Results and discussion

    3.1. 2D PIV results

    Fig. 4 shows the phase-averaged spanwise vorticity (ωzc/U∞)field of the PIV results. The experimental measurement results have been ensemble-averaged over 38 cycles.In the early stage of the up-stroke phase,the flow remains attached to the airfoil.When the angle of attack increases to 15.5°,a vorticity concentration zone is observed at the leading edge of the airfoil,which means that the LEV forms and continually gains vorticity from the leading edge shear layer.With the development of the LEV, a weak vortex with an opposite sign forms between the shear layer and the LEV when the angle of attack increases to 21.67°. When the angle of attack increases to 22.73°, which is nearly the maximum angle, the LEV detaches from the airfoil surface, and a new LEV is formed.

    When the airfoil is in the down-stroke phase, the flow quickly transitions to a full separation zone, which is compressed with increases of the angle of attack. The detachment of the LEV is different from the phenomenon observed in experimental results at a lower Reynolds number from Gharali and Johnson18and Mackowski and Williamson14that the LEV detaches from the airfoil until the LEV reaches the trailing edge of the airfoil. In addition, the LEV in the phaseaveraged flow field is more stable than the experimental results at a higher Reynolds number from Taylor and Amitay23and Zanotti et al.25

    At high Reynolds number, the flow field is highly random.The phase-averaging process erases some of the flow features.Fig. 5 shows the instantaneous spanwise vorticity field at α=18.57° and 22.73° and the corresponding phase-averaged results to explore the instantaneous flow features. At low angles of attack up to α=18.57°, the phase-averaged fields are similar to the instantaneous fields. However, the discrepancy increases with the angle of attack. At high angles of attack around 22.73°, the main vortex structure of the phaseaveraged field is the same as the instantaneous field.However,small-scale structures form and shed during the pitching motion. These structures are not phase-locked since they are filtered out in the phase-averaging process.

    The reason for these differences is probably that there are three-dimensional structures produced by instability at high angles of attack. They interact with the LEV and generate the small-scale structures in the two-dimensional PIV results.Thus,there are different scales of the vortices in the flow field,which is the typical flow feature at transitional and high Reynolds numbers. The complex mutual process of vortices with different scales enlarge the noises of the flow information captured by the PIV,which further increases the noisy level of the calculated PIV data. Thus, the instantaneous flow fields seem to be disordered. This phenomenon has also been observed at high Reynolds numbers by other researchers,such as Taylor and Amitay23, Mulleners and Raffel24, and Zanotti et al.25

    Fig. 4 Phase-averaged spanwise vorticity fields at k=0.1 for PIV results.

    Fig. 5 Comparison between the phase-averaged and instantaneous vorticity fields at the same angle of attack.

    Fig. 6 Normalized turbulence statistics at different angles of attack.

    The experimental study of a pitching airfoil indicates that at the chosen transitional Reynolds number, the small-scale flow structures that are filtered by the phase-averaging process contribute substantially to the turbulence at a high angle of attack.To better study these phenomena, a 3D CFD method which performs well in capturing instantaneous small-scale flow structures is used to overcome the shortcomings of filtering such structures by phase-averaged process. Fig. 7 compares the instantaneous flow field between the PIV and CFD results.The CFD results are from the frame at the middle span. At α=19°, the LEV structures show good agreement. When α increases to 20.9°, the LEVs between PIV results and CFD results still show high similarity. However, the differences increase with the angle of attack. The main differences are the vortex strength and the secondary vortex with an opposite sign, as in the case of α=20.9°. When the angle of attack increases to 22.73°,the LEV of the CFD results is more intense than in the PIV results. The signal-to-noise ratio of the CFD results is also lower. The secondary vortex is not clear in the PIV results because the secondary vortex may lose stability and split into small-scale structures that are covered by the noise. Nevertheless, the formation and detachment of the LEV have good agreement between the CFD and experimental results.This proves that the 3D CFD can be used to study the details of the 3D flow field.

    3.2. 3D CFD results

    Fig.8 shows the spanwise vorticity fields in the mid-span of the airfoil at k=0.1 from the CFD results.In Fig.8(a),the angle of attack is α=19.8°,and the LEV is growing in strength and size.Von Karman streets also shed into the wake near the trailing edge because of trailing edge separation. In Fig. 8(b), the angle of attack is increased to 23.7°. Two LEVs have already detached from the leading-edge shear layer, and the first LEV induces a secondary vortex with the opposite sign beneath the LEV. Additionally, the Trailing Edge Vortex(TEV) begins to develop. One vortex that rotates clockwise can be observed above the TEV. These vortex structures are still stable and move downstream together.

    In Fig.8(c),the angle of attack is increased to 24.68°,which is near the maximum angle of attack. One obvious feature of the flow field at this angle of attack is that the flow structures start to lose their stability and deform. The secondary vortex has split into small-scale structures along the airfoil surface.In Fig.8(d),the airfoil just enters the down-stroke phase,and the angle of attack is 24.61°. The original LEV breaks down into complex small structures. Another new TEV is formed and moves the separation zone upward.

    Fig. 7 Comparison between CFD results and PIV results.

    Since the flow structures are highly unstable at high angles of attack, it is worthwhile to study the three-dimensional flow field, as shown in Fig. 9.In the three-dimensional case, the Qcriterion isosurface is used as an indicator of turbulent flow structures.Q-criterion is based on the calculation of the second invariant of the velocity gradient tensor.It has been proved as an efficient method to separate the vortex from a shear layer.35Additionally,the isosurface is colored according to the magnitude of the streamwise vorticity to show the strength of the flow structure.In the early stage of the up-stroke phase shown in Fig. 9(a),the LEV extends strictly along the span. With the increase of the angle of attack, numerous LEVs detach from the leading-edge shear layer and move downstream continually, as shown in Fig. 9(b). A cluster of streamwise vortices forms and arrange along the span direction in the rear part of the airfoil.At the same time,the spanwise vortices start losing their stability and have weak three-dimensional deformation.

    In Fig.9(c),the angle of attack increases to 24.68°,and the streamwise vortices move downstream and expand. At the same time, some streamwise vortices can also be observed between the first and second LEVs. When the angle of attack reaches the maximum and then decreases to 24.6°, the airfoil enters the deep stall process, and the flow field is completely turbulent.Those streamwise vortices are also observed in other researches to pitching airfoil like Buchner et al.9and Garmann and Visbal36by experimental and numerical simulation methods. However, the Reynolds number they used was on the order of 104which is quite lower than this research. In addition, the motion function they used was pitching at constant velocity not the harmonic function in this research. It states that the streamwise vortices are not special flow structures in this experimental parameters. The reason to its formation may due to the stability of the LEV.

    To better study the development of the streamwise vortices,the vorticity field in four y-z frames are shown in Fig.10.These frames are fixed at the present coordinates and do not move with the airfoil.Fig.11 shows the corresponding vorticity field.In Fig. 11(a), the angle of attack is 22.7°, and the streamwise vortices spread over all frames, especially the frame located at 1.25c.These vortices are arranged at the same height in pairs with a mushroom shape in one frame.The width of the vortex pair in one frame is nearly the same.In Fig.11(b),the angle of attack is 24.8°, and the airfoil starts to enter the down-stroke phase.The streamwise vortices keep developing,and the mushroom shape cannot be maintained. The vortices at the same height merge with each other and form a vorticity belt, as shown in the frames located at 0.5c and 1c. However, the streamwise vortices elongate in the vertical direction in the frame located at 0.75c.

    Fig. 8 Two-dimensional spanwise vorticity field for the CFD results at k=0.1.

    In Fig.11(c)the angle of attack decreases to 24.3°,the flow transitions to turbulence, and distinct vortices cannot be observed. In the development of the streamwise vortices, their characteristics include a spanwise-invariant width, attachment to the LEV, and merging. These substructures are assumed to be formed because of the centrifugal instability of the LEV.In Fig. 11(d), the angle of attack is decreases further, and the wake behind the airfoil shrinks as the vortex strength gradually loses strength. With dynamic stall, the core vortex is strong enough to remain as a coherent structure, while with decreasing of angle of attack,the vortices lose their energy and appear as large eddies, which do not have any order and diminish rapidly.

    3.3. Comparison of different k

    Fig. 10 Schematic of four chosen y-z frames.

    Fig. 11 Streamwise vorticity in four y-z frames for the CFD results.

    Fig. 12 shows the vortex development for a pitching airfoil at k=0.2 from the CFD results.In Fig.12(a),the airfoil is in the up-stroke phase, and the angle of attack increases to α=24.37°. In this instant, the LEV has already developed and separates from the airfoil surface. Additionally, the TEV also begins to grow. In Fig. 12(b), the airfoil is in the downstroke phase, and the angle of attack decreases to α=18°.In this instant, the LEV has already separated and sheds into the wake. Additionally, the TEV forms and is ready to leave the airfoil surface. The secondary LEV is still attached to the airfoil, even though the airfoil is in the down-stroke phase.This suggests that the response time of the flow field is slower at high k. With a further decrease of the angle of attack to α=11.9°, the flow mostly separates from the airfoil surface,but no full separation zone is observed on the suction side of the airfoil.

    Fig. 12 Two-dimensional streamwise vorticity field for the CFD results at k=0.2.

    Fig.13 shows the three-dimensional vortex development at k=0.2. When the angle of attack nearly reaches the maximum, as shown in Fig. 13(a), the LEV has detached from the shear layer. The LEVs and wake streets show good two-dimensional characteristics. In Fig. 13(b), a cluster of streamwise vortices can be observed between the primary and secondary LEVs.The leading-edge shear layer transforms into a streamwise vortex,which cannot be observed in the case of k=0.1.In Fig.13(c),streamwise vortices develop and multiply.At the same time,the TEV formats with streamwise vortices which is formed at the trailing edge of the airfoil and forces the streamwise vortex to separate from the airfoil surface. When the airfoil enters deep stall, the flow field is again full of small structures.

    Fig. 13 Q-criterion isosurface of the streamwise vorticity for the CFD results at k=0.2.

    There is a big discrepancy in the vortex development between the flow fields at the two reduced frequencies. The LEV detaches from the shear layer later and remains stable for a longer time at higher reduced frequency. Furthermore,three LEVs at most can be observed at k=0.1, while there are only two LEVs at k=0.2.The streamwise vortices remain longer in the case of k=0.2. This shows that the airfoil reaches stall later at higher k,which means the rotation velocity is higher. Besides that, it should be mentioned that the streamwise vortices format after the detachment and three dimensional deformation of the LEV at two different reduced frequencies.It can be assumed that the instability which occurs the deformation of the LEV is the reason to the formation of those streamwise vortices.

    Fig. 14 Lift and drag coefficients of the pitching airfoil at k=0.1 and 0.2 obtained by CFD.

    Fig. 14 shows the lift and drag coefficients of the pitching airfoil at different reduced frequencies. In Fig. 14(a), the lift coefficient increases linearly with increasing angle of attack in the pre-stall region at both k=0.1 and k=0.2. The lift coefficient at different k increases linearly in the pre-stall region, fluctuates in the dynamic stall region, and rapidly decreases in the stall region. In the pre-stall region, the lift coefficient at k=0.2 is a little bit higher than that at k=0.1. The maximum lift coefficient at different k is the same,while the corresponding angle of attack is different.This means that the airfoil reaches the stall stage later at k=0.2.The lift fluctuation at k=0.1 is also stronger than in the other case.

    The angle of attack where the lift coefficient starts to fluctuate is the same as the time when the LEV detaches from the airfoil.Since the LEV at k=0.1 detaches earlier and more LEVs are formed during the pitching motion, the vortex lift that influences the lift coefficient is much stronger. At higher k, the airfoil also reaches stall during the down-stroke stage.This means that the LEV is strong enough to keep the flow relatively attached and stable at higher k. Mackowssi and Williamson14found that the lift agrees well with Theodorsen linear theory37which states the amplitude of angle of attack determines the maximum lift at low k and low Reynolds number.The same maximum lift in both cases prove that this phenomenon also exists in transitional Reynolds number.

    Fig. 14(b) shows the drag coefficient hysteresis. In general,the drag coefficient for both reduced frequencies increases with the angle of attack.The wake behind the airfoil and separation of flow over the suction side of the airfoil are the reasons why the drag increases with the angle of attack. At k=0.1, the drag coefficient surpasses that at k=0.2 at high angle of attack. This occurs because the dynamic stall occurs right before the maximum angle of attack,which leads to high drag force on the airfoil. But in the case of k=0.2, the flow is not fully separated, which is why the drag coefficient is relatively low in the up-stroke phase. Once the airfoil starts the downstroke phase,the dynamic stall happens at k=0.2,and at this moment, the drag coefficient remains high,as opposed to that at k=0.1.

    4. Conclusions

    This research focused on the vortex behavior and aerodynamic force of a NACA0012 airfoil during the dynamic stall stage under two reduced frequencies at a transitional Reynolds number of Re=90000. To better study the flow field and aerodynamic force, the two-dimensional time-resolved PIV technique and three-dimensional CFD were used. The angle of attack varied from -5° to 25° and followed a sinusoidal function.

    The flow field of the pitching airfoil at k=0.1 and a transitional Reynolds number was first studied experimentally.The flow is attached to the airfoil in the early up-stroke phase.The LEV then takes shape, which indicates that the airfoil reaches dynamic stall. After the detachment of the LEV, the flow quickly transitions to a full separation zone. When comparing the phase-averaged results and the instantaneous results,there were small flow structures in the flow field with increasing the angle of attack. The results of the turbulent statistics also proves this.To better understand the flow field of the pitching airfoil, three-dimensional DDES was used to study the flow field and aerodynamic force at the same Reynolds number and two reduced frequencies of k=0.1 and k=0.2.

    Similar flow development was found in the CFD results.Streamwise vortices were also found in the three-dimensional results.The streamwise vortices were arranged along the spanwise direction.The interaction between the streamwise vortices and the LEV promotes the transition of the flow field to turbulence. When comparing the results between the two reduced frequencies, the flow was more stable and the airfoil reaches the stall stage later at higher reduced frequency.The vortex lift produced by the LEV contributes substantially to the lift in the dynamic stall stage. The formation and detachment of the LEV also cause fluctuation of the lift coefficient in the dynamic stall stage.

    Acknowledgments

    This work was supported by the National Natural Science Foundation of China (Nos. GZ 1280, 11722215 and 11721202). This work was also supported by the National Research Foundation of Korea (NRF) grant with funding from the Korean government (MSIT) (No. 2011-0030013,No. 2018R1A2B2007117).

    黑人高潮一二区| 中国国产av一级| 亚洲精品成人久久久久久| 欧美又色又爽又黄视频| 嫩草影视91久久| 日日啪夜夜撸| 亚洲精品成人久久久久久| 国产精品免费一区二区三区在线| 国产男人的电影天堂91| 成人特级av手机在线观看| 亚洲av第一区精品v没综合| 一区二区三区免费毛片| 成人性生交大片免费视频hd| 欧美成人精品欧美一级黄| 欧美一级a爱片免费观看看| 特级一级黄色大片| av在线观看视频网站免费| 99在线人妻在线中文字幕| 久久久久国内视频| 一卡2卡三卡四卡精品乱码亚洲| 久久久久久久久久成人| 亚洲成人中文字幕在线播放| 欧美zozozo另类| 成年女人毛片免费观看观看9| 亚洲成av人片在线播放无| 97碰自拍视频| 深爱激情五月婷婷| 国产v大片淫在线免费观看| 国产成人一区二区在线| 波多野结衣高清无吗| 有码 亚洲区| 精品欧美国产一区二区三| 亚洲精品色激情综合| av福利片在线观看| 日本三级黄在线观看| 亚洲图色成人| 欧美高清性xxxxhd video| 直男gayav资源| 九九爱精品视频在线观看| 简卡轻食公司| 国产成人aa在线观看| 丝袜美腿在线中文| 一本久久中文字幕| 国产人妻一区二区三区在| 一级av片app| 最近最新中文字幕大全电影3| 日韩成人伦理影院| 国产精品1区2区在线观看.| 亚洲成人久久性| 成人三级黄色视频| 香蕉av资源在线| 舔av片在线| 18禁在线无遮挡免费观看视频 | 亚洲av二区三区四区| 18禁黄网站禁片免费观看直播| 我要搜黄色片| 欧洲精品卡2卡3卡4卡5卡区| 午夜福利在线观看免费完整高清在 | 日韩人妻高清精品专区| 成年女人永久免费观看视频| 亚洲一区高清亚洲精品| 99久久九九国产精品国产免费| 欧美bdsm另类| 亚洲精品国产av成人精品 | 国产av不卡久久| 99久国产av精品国产电影| 亚洲av成人精品一区久久| 久久精品91蜜桃| 日韩三级伦理在线观看| 日韩欧美精品免费久久| 午夜视频国产福利| 18禁黄网站禁片免费观看直播| 两个人视频免费观看高清| 极品教师在线视频| 全区人妻精品视频| 成年女人永久免费观看视频| 又爽又黄无遮挡网站| 国产探花极品一区二区| 日日撸夜夜添| 人妻少妇偷人精品九色| 欧美激情国产日韩精品一区| 久久婷婷人人爽人人干人人爱| 国产精品久久视频播放| 99精品在免费线老司机午夜| 日本免费一区二区三区高清不卡| 日本免费a在线| 老司机午夜福利在线观看视频| 伊人久久精品亚洲午夜| 久久中文看片网| 最近手机中文字幕大全| 国产在线男女| 精品人妻视频免费看| av在线播放精品| 日本在线视频免费播放| 亚洲电影在线观看av| 国产av不卡久久| 少妇人妻精品综合一区二区 | 亚洲av二区三区四区| 国产私拍福利视频在线观看| 乱系列少妇在线播放| 成人精品一区二区免费| 午夜a级毛片| 亚洲婷婷狠狠爱综合网| 日日摸夜夜添夜夜添小说| 日韩精品青青久久久久久| 偷拍熟女少妇极品色| 蜜桃亚洲精品一区二区三区| 在线观看免费视频日本深夜| 九九热线精品视视频播放| 最近的中文字幕免费完整| 直男gayav资源| 亚洲欧美日韩东京热| av中文乱码字幕在线| 国产精品日韩av在线免费观看| 国产精品美女特级片免费视频播放器| 丝袜喷水一区| 成年版毛片免费区| 日韩欧美精品v在线| 欧美一区二区国产精品久久精品| 久久韩国三级中文字幕| 欧美3d第一页| 亚洲性久久影院| 成年女人看的毛片在线观看| 免费看a级黄色片| 精品午夜福利在线看| 日韩av在线大香蕉| 久久精品国产亚洲网站| 天堂网av新在线| 九九在线视频观看精品| av国产免费在线观看| 免费大片18禁| 99在线视频只有这里精品首页| 日本三级黄在线观看| 特级一级黄色大片| 久久欧美精品欧美久久欧美| 欧美精品国产亚洲| 春色校园在线视频观看| 九九爱精品视频在线观看| 大香蕉久久网| 国内揄拍国产精品人妻在线| 好男人在线观看高清免费视频| 亚洲av中文字字幕乱码综合| 3wmmmm亚洲av在线观看| 国产片特级美女逼逼视频| 亚洲人成网站在线观看播放| 日韩制服骚丝袜av| 久久精品国产亚洲网站| 亚洲经典国产精华液单| 淫秽高清视频在线观看| 亚洲第一电影网av| 国产av在哪里看| av视频在线观看入口| av黄色大香蕉| 此物有八面人人有两片| 一级毛片我不卡| 国产精品野战在线观看| 亚洲欧美成人综合另类久久久 | 丝袜美腿在线中文| 精品久久久久久久人妻蜜臀av| 又粗又爽又猛毛片免费看| 日本a在线网址| 中文字幕人妻熟人妻熟丝袜美| 12—13女人毛片做爰片一| 日韩欧美免费精品| 国产精品永久免费网站| 国产三级中文精品| av卡一久久| 免费观看精品视频网站| 欧美日韩在线观看h| 亚洲精品国产成人久久av| 人妻少妇偷人精品九色| 免费黄网站久久成人精品| av专区在线播放| 欧美三级亚洲精品| 我的女老师完整版在线观看| 欧美xxxx黑人xx丫x性爽| 极品教师在线视频| 亚洲人成网站在线观看播放| 国产成人精品久久久久久| 精品一区二区免费观看| 卡戴珊不雅视频在线播放| 成人二区视频| 91狼人影院| 午夜福利在线在线| 国产一区二区三区在线臀色熟女| 老司机影院成人| 国产综合懂色| 看十八女毛片水多多多| 成人国产麻豆网| 亚洲aⅴ乱码一区二区在线播放| 国产精品伦人一区二区| 又黄又爽又免费观看的视频| 精品久久久久久久久亚洲| 毛片女人毛片| 久久久久国产精品人妻aⅴ院| 亚洲av免费在线观看| 欧美丝袜亚洲另类| 精品一区二区三区人妻视频| 少妇猛男粗大的猛烈进出视频 | 色哟哟哟哟哟哟| 天美传媒精品一区二区| 麻豆一二三区av精品| 最新中文字幕久久久久| 一级黄片播放器| .国产精品久久| 变态另类丝袜制服| 欧美成人精品欧美一级黄| 国内精品美女久久久久久| 偷拍熟女少妇极品色| 观看免费一级毛片| 看非洲黑人一级黄片| 亚洲国产色片| 在线观看免费视频日本深夜| 岛国在线免费视频观看| 午夜日韩欧美国产| 男女做爰动态图高潮gif福利片| 狂野欧美白嫩少妇大欣赏| 国产黄片美女视频| 日韩欧美免费精品| 少妇的逼水好多| 色综合亚洲欧美另类图片| 在线免费观看不下载黄p国产| 国产高清不卡午夜福利| 岛国在线免费视频观看| 中文资源天堂在线| av免费在线看不卡| 中文字幕免费在线视频6| ponron亚洲| 国产成人a∨麻豆精品| 真人做人爱边吃奶动态| 婷婷六月久久综合丁香| 女人被狂操c到高潮| 国产成人精品久久久久久| 一区二区三区四区激情视频 | 国产精品综合久久久久久久免费| 卡戴珊不雅视频在线播放| 小蜜桃在线观看免费完整版高清| 99久久中文字幕三级久久日本| 国产亚洲欧美98| 国产一级毛片七仙女欲春2| 岛国在线免费视频观看| 在线播放无遮挡| 国产成年人精品一区二区| 欧美3d第一页| 午夜福利在线在线| 国产人妻一区二区三区在| 欧美最新免费一区二区三区| av视频在线观看入口| 麻豆一二三区av精品| 香蕉av资源在线| 干丝袜人妻中文字幕| 我的老师免费观看完整版| 欧美zozozo另类| 午夜日韩欧美国产| 国产一区二区在线av高清观看| 国产精品久久电影中文字幕| 深夜a级毛片| 简卡轻食公司| 亚洲一区二区三区色噜噜| av视频在线观看入口| 看黄色毛片网站| 久久久成人免费电影| 成年免费大片在线观看| 久久精品人妻少妇| 在现免费观看毛片| 中文亚洲av片在线观看爽| 午夜亚洲福利在线播放| 亚洲国产日韩欧美精品在线观看| 欧美一区二区国产精品久久精品| 啦啦啦韩国在线观看视频| 男女做爰动态图高潮gif福利片| 国产精品一二三区在线看| 午夜福利在线观看免费完整高清在 | 99久久无色码亚洲精品果冻| 久久鲁丝午夜福利片| 丰满乱子伦码专区| 国产伦精品一区二区三区四那| 国产精品一及| 又爽又黄a免费视频| 日韩欧美精品v在线| 熟女人妻精品中文字幕| 国产免费一级a男人的天堂| 夜夜爽天天搞| 黄色配什么色好看| 观看美女的网站| 亚洲人成网站在线观看播放| 国产高清视频在线观看网站| 精品久久久久久久久久免费视频| 日日撸夜夜添| 精品久久久久久久久久久久久| 18禁黄网站禁片免费观看直播| 非洲黑人性xxxx精品又粗又长| 亚洲中文字幕一区二区三区有码在线看| 丝袜美腿在线中文| 亚洲国产精品合色在线| 深爱激情五月婷婷| 久久久久免费精品人妻一区二区| 最近视频中文字幕2019在线8| 中文字幕久久专区| 长腿黑丝高跟| 免费观看精品视频网站| 长腿黑丝高跟| 精品国内亚洲2022精品成人| 国产男人的电影天堂91| 国产久久久一区二区三区| 国产真实伦视频高清在线观看| 99热6这里只有精品| 人妻夜夜爽99麻豆av| 久久久久九九精品影院| 国产熟女欧美一区二区| 18禁裸乳无遮挡免费网站照片| 黄色配什么色好看| videossex国产| 欧美激情久久久久久爽电影| 久久6这里有精品| 亚洲成av人片在线播放无| 亚洲内射少妇av| 精品一区二区三区av网在线观看| 亚洲综合色惰| 亚洲国产精品国产精品| 国产成人91sexporn| 女的被弄到高潮叫床怎么办| 国产一区二区在线观看日韩| 午夜精品国产一区二区电影 | 亚洲成人中文字幕在线播放| 欧美成人精品欧美一级黄| 欧美一区二区精品小视频在线| 亚洲av一区综合| 日韩大尺度精品在线看网址| 国产精品一及| 一级av片app| 色哟哟哟哟哟哟| 国产亚洲精品久久久com| 亚洲经典国产精华液单| 中国美女看黄片| 看十八女毛片水多多多| 丝袜美腿在线中文| 老熟妇乱子伦视频在线观看| 美女内射精品一级片tv| 亚洲,欧美,日韩| 在线免费观看的www视频| 午夜激情欧美在线| 不卡视频在线观看欧美| 亚洲自偷自拍三级| 国产中年淑女户外野战色| 大又大粗又爽又黄少妇毛片口| 亚洲乱码一区二区免费版| videossex国产| 少妇的逼水好多| 菩萨蛮人人尽说江南好唐韦庄 | 青春草视频在线免费观看| 九九久久精品国产亚洲av麻豆| 成年av动漫网址| 久久人人精品亚洲av| 99热这里只有是精品在线观看| 三级国产精品欧美在线观看| 人妻少妇偷人精品九色| 中文字幕av成人在线电影| 欧美色视频一区免费| 国产成人精品久久久久久| 午夜精品在线福利| 一级黄色大片毛片| 美女cb高潮喷水在线观看| 亚洲中文字幕日韩| 久久久色成人| 日韩一本色道免费dvd| 好男人在线观看高清免费视频| 欧美日韩一区二区视频在线观看视频在线 | 国产 一区精品| 乱系列少妇在线播放| 天天躁夜夜躁狠狠久久av| 最近2019中文字幕mv第一页| 神马国产精品三级电影在线观看| 最新在线观看一区二区三区| 日韩中字成人| 在线播放国产精品三级| 99热只有精品国产| 大又大粗又爽又黄少妇毛片口| 哪里可以看免费的av片| 久久久欧美国产精品| 最近视频中文字幕2019在线8| 日日摸夜夜添夜夜添av毛片| а√天堂www在线а√下载| 国产成人freesex在线 | 一区二区三区高清视频在线| 亚洲人与动物交配视频| 午夜影院日韩av| 日本黄色片子视频| 少妇丰满av| 精品久久久久久成人av| 免费av毛片视频| 国产高清三级在线| 国产一区二区在线av高清观看| 国产激情偷乱视频一区二区| 成人漫画全彩无遮挡| 亚洲无线观看免费| 欧美zozozo另类| 亚洲av二区三区四区| 最后的刺客免费高清国语| 欧美三级亚洲精品| 日韩 亚洲 欧美在线| 久久精品人妻少妇| 女人被狂操c到高潮| www日本黄色视频网| 欧洲精品卡2卡3卡4卡5卡区| 久久久久久久久久黄片| 看片在线看免费视频| 久久久精品欧美日韩精品| 国产亚洲精品av在线| 少妇裸体淫交视频免费看高清| 搞女人的毛片| 免费搜索国产男女视频| 欧美日韩一区二区视频在线观看视频在线 | 国产成人a∨麻豆精品| 在线观看免费视频日本深夜| 亚洲国产色片| 2021天堂中文幕一二区在线观| 在线看三级毛片| 日韩欧美精品v在线| 老司机福利观看| 欧美bdsm另类| 国内久久婷婷六月综合欲色啪| 色综合亚洲欧美另类图片| 老师上课跳d突然被开到最大视频| 丰满人妻一区二区三区视频av| 久久人妻av系列| 国产私拍福利视频在线观看| 免费观看精品视频网站| 国产一区二区亚洲精品在线观看| 尾随美女入室| 99国产精品一区二区蜜桃av| 亚洲精品一卡2卡三卡4卡5卡| 性色avwww在线观看| 免费观看在线日韩| 国产精品久久视频播放| 不卡视频在线观看欧美| 九九久久精品国产亚洲av麻豆| 国产精品av视频在线免费观看| 亚洲欧美日韩高清专用| 亚洲国产日韩欧美精品在线观看| 日本精品一区二区三区蜜桃| 中出人妻视频一区二区| 伦理电影大哥的女人| 国国产精品蜜臀av免费| 老师上课跳d突然被开到最大视频| 美女大奶头视频| 午夜爱爱视频在线播放| 国产熟女欧美一区二区| 日本爱情动作片www.在线观看 | 蜜桃久久精品国产亚洲av| 欧美日韩一区二区视频在线观看视频在线 | 乱人视频在线观看| 天堂影院成人在线观看| www日本黄色视频网| 久久精品夜色国产| 亚洲欧美精品自产自拍| 一区二区三区免费毛片| 伦精品一区二区三区| 亚州av有码| 99久国产av精品| 欧美日本亚洲视频在线播放| 美女cb高潮喷水在线观看| 国产精品乱码一区二三区的特点| 老司机影院成人| 天天一区二区日本电影三级| 亚洲精品久久国产高清桃花| 亚洲中文字幕一区二区三区有码在线看| 人妻少妇偷人精品九色| 三级国产精品欧美在线观看| 国产日本99.免费观看| 免费看a级黄色片| 特大巨黑吊av在线直播| 国产蜜桃级精品一区二区三区| 国内精品久久久久精免费| 91久久精品国产一区二区三区| 一个人看视频在线观看www免费| 国产成人a∨麻豆精品| 亚洲电影在线观看av| 在线免费观看不下载黄p国产| 3wmmmm亚洲av在线观看| 午夜精品一区二区三区免费看| 乱人视频在线观看| 变态另类成人亚洲欧美熟女| 午夜爱爱视频在线播放| 欧美日韩精品成人综合77777| 99久久成人亚洲精品观看| 亚洲精品色激情综合| 午夜免费男女啪啪视频观看 | 身体一侧抽搐| 国产美女午夜福利| 男人舔女人下体高潮全视频| 久久久久精品国产欧美久久久| 国产一区二区在线av高清观看| 国产精品不卡视频一区二区| 国产精品久久视频播放| 又粗又爽又猛毛片免费看| 成熟少妇高潮喷水视频| 天堂√8在线中文| 成人永久免费在线观看视频| 在线观看午夜福利视频| 欧美又色又爽又黄视频| 久久精品国产鲁丝片午夜精品| 成年女人永久免费观看视频| 99久久九九国产精品国产免费| 亚洲欧美日韩高清专用| 91av网一区二区| 国产在线精品亚洲第一网站| 免费看美女性在线毛片视频| 成人特级黄色片久久久久久久| 热99re8久久精品国产| 欧美日韩精品成人综合77777| 国产成人a区在线观看| 最近视频中文字幕2019在线8| 免费av不卡在线播放| 欧美xxxx性猛交bbbb| 最新中文字幕久久久久| 国产激情偷乱视频一区二区| 毛片一级片免费看久久久久| 大型黄色视频在线免费观看| 内地一区二区视频在线| 色综合站精品国产| 波野结衣二区三区在线| 国产精品美女特级片免费视频播放器| 日韩精品有码人妻一区| 亚洲丝袜综合中文字幕| 免费搜索国产男女视频| 91午夜精品亚洲一区二区三区| 九九爱精品视频在线观看| 色哟哟·www| 国产白丝娇喘喷水9色精品| 熟女人妻精品中文字幕| 国产精品av视频在线免费观看| 真人做人爱边吃奶动态| 久久天躁狠狠躁夜夜2o2o| 波多野结衣高清无吗| 尤物成人国产欧美一区二区三区| 乱人视频在线观看| 日本黄色视频三级网站网址| 国产亚洲欧美98| 亚洲av免费在线观看| 熟妇人妻久久中文字幕3abv| 久久国内精品自在自线图片| 国产av在哪里看| 久久人妻av系列| 综合色丁香网| 免费在线观看影片大全网站| 极品教师在线视频| 亚洲精品粉嫩美女一区| 搡老熟女国产l中国老女人| 亚洲av熟女| 亚洲av二区三区四区| 搡老岳熟女国产| 午夜a级毛片| 国内精品久久久久精免费| 99九九线精品视频在线观看视频| 国产精品一区二区三区四区免费观看 | 国产精品一区www在线观看| 国产一区二区三区av在线 | 日韩一区二区视频免费看| 我要搜黄色片| 精品不卡国产一区二区三区| 97热精品久久久久久| 一区二区三区高清视频在线| 国产中年淑女户外野战色| 深夜a级毛片| 精品午夜福利视频在线观看一区| 国产亚洲91精品色在线| 欧美不卡视频在线免费观看| 非洲黑人性xxxx精品又粗又长| 久久精品国产清高在天天线| 欧美绝顶高潮抽搐喷水| 午夜爱爱视频在线播放| 国产成人91sexporn| 国产精品亚洲美女久久久| 国产亚洲精品久久久com| 五月伊人婷婷丁香| 真人做人爱边吃奶动态| 18禁裸乳无遮挡免费网站照片| 亚洲欧美日韩无卡精品| 99热精品在线国产| 18禁在线无遮挡免费观看视频 | 在线播放国产精品三级| 久久综合国产亚洲精品| 12—13女人毛片做爰片一| 日本在线视频免费播放| 国产一区二区三区在线臀色熟女| av.在线天堂| 欧美另类亚洲清纯唯美| 国产黄a三级三级三级人| 国语自产精品视频在线第100页| 欧美色欧美亚洲另类二区| 国产精品一区二区三区四区久久| 免费无遮挡裸体视频| 国产精品一二三区在线看| 长腿黑丝高跟| 三级男女做爰猛烈吃奶摸视频| 干丝袜人妻中文字幕| 国产一区二区在线av高清观看| 欧美日韩综合久久久久久| 成年女人毛片免费观看观看9| 国产精品日韩av在线免费观看| 在线天堂最新版资源| 人妻久久中文字幕网| 中文字幕熟女人妻在线| 免费观看在线日韩| 亚洲人成网站在线播| 18禁在线播放成人免费| 91久久精品电影网| 中文在线观看免费www的网站| 日本免费a在线| 日韩一本色道免费dvd| 十八禁国产超污无遮挡网站| 69av精品久久久久久| 神马国产精品三级电影在线观看| 天堂网av新在线|