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

    A methodology for simulating 2D shock-induced dynamic stall at flight test-based fluctuating freestream

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

    Khider AL-JABURI, Daniel FESZTY, Fred NITZSCHE

    Department of Mechanical and Aerospace Engineering, Carleton University, Ottawa K1S 5B6, Canada

    KEYWORDS

    Mach fluctuation;

    Pitching airfoil;

    Shock-induced dynamic stall;

    Transonic flow;

    Unsteady freestream velocity

    Abstract A comprehensive methodology for simulating 2D dynamic stall at fluctuating freestream is proposed in this paper. 2D CFD simulation of a SC1095 airfoil exposed to a fluctuating freestream of Mach number 0.537±0.205 and Reynolds number 6.1×106 (based on the mean Mach number)and undergoing a 10°±10° pitch oscillation with a frequency of 4.25 Hz was conducted.These conditions were selected to be representative of the flow experienced by a helicopter rotor airfoil section in a real-life fast forward flight. Both constant freestream dynamic stall as well as fluctuating freestream dynamic stall simulations were conducted and compared.The methodology was carefully validated with experimental data for both transonic flow and dynamic stall under fluctuating freestream.Overall,the results suggest that the fluctuating freestream alters the dynamic stall mechanism documented for constant freestream in a major way,emphasizing that inclusion of this effect in the prediction of dynamic stall related rotor loads is imperative for rotor performance analysis and blades design.

    1. Introduction

    Dynamic stall is considered as one of the most complicated flow phenomena in the field of rotary-wing aerodynamics and in particular for helicopters. For a helicopter in fast forward flight, dynamic stall appears on the retreating rotor blades due to the high frequency periodic variation of the effective angle of attack of a blade section.1It is because dynamic stall is associated with the increase of lift, drag and negative pitching moment,as well as the occurrence of a characteristic hysteresis in these variables. The combination of all these effects leads to the appearance of excessive vibrations in the rotor hub as well as the control system of the helicopter,which not only deteriorates crew and passenger comfort but also limits the forward flight speed of the helicopter.1-3However, dynamic stall on a helicopter rotor blade occurs at the simultaneous variation of the relative freestream as well as of the effective angle of attack. At moderate and fast advance ratios, the amplitude of the relative freestream fluctuation can be as high as Mach number 0.3, the amplitude of the angle of attack as high as 10° and all these occur at the rotational frequency,i.e.at about 3-5 Hz.Unfortunately,reproduction of a fluctuating freestream variation of Mach number 0.3 amplitude and at a frequency of 3-5 Hz is extremely challenging in ground-based experiments. Therefore, the vast majority of dynamic stall studies have historically focused on considering only the variation of the pitch angle, while neglecting the fluctuation in the freestream.4-8

    There have been only a handful of researchers who conducted experimental or CFD studies of dynamic stall by considering fluctuating freestream as well. Gosselin and Feszty9has conducted a CFD study on dynamic stall with fluctuating freestream, which was dedicated to exploring transonic effects on dynamic stall in fast forward flight, but with validating the simulation only at the lower extreme of the freestream fluctuation as well as with a narrow scope of the influence of numerical parameters,such as the size of the computational domain.Table 1 provides a summarized literature survey.9-26

    Kerho20Gharali et al.22as well as Glaz et al.27,28have conducted CFD simulations of dynamic stall with fluctuating freestream,but either without validating the results to experiments and/or without exploring the transonic effects that could appear and affect the flow hysteresis.20,27,28In the case of Gharali et al.22the conditions were representative of wind turbines and not helicopters. Fernie29and Hird24,25et al. conducted experimental studies, but with freestream fluctuations corresponding to very modest advance ratios. Recently, Zhao et al.23conducted a CFD study with Mach number ranged from 0.2 -0.6 in an effort to study the effect of such a freestream velocity fluctuation on the aerodynamic characteristics of an airfoil subjected to oscillation. However, although the Mach number range suggested the occurrence of transonic flow on the advancing blade section, its effect was absent in the hysteresis loops and flow visualizations.

    It is clear from the above, that there is limited information regarding the methodology of simulating dynamic stall at fluctuating freestream as well as there is an absence of experimental results at the Mach number amplitudes corresponding to the real-life fluctuations of the relative freestream acting on a blade section. Therefore, the objective of this work is: (A) to provide a guideline on how to simulate 2D dynamic stall at real helicopter forward flight conditions,i.e.at fluctuating freestream yielding both transonic flow as well as dynamic stall within the same cycle, as well as (B) to determine whether the occurrence of transonic flow on the advancing side of the rotor affects the dynamic stall phenomenon on the retreating side of the rotor disc.

    2. Test case

    To provide realistic flight conditions as much as possible, the present study considered the conditions of UH-60A Blackhawk helicopter in fast forward flight.Table 2 lists the simulation parameters and conditions. The operating parameters of UH-60A as well as the flight conditions were selected with the help of references.30-33The simulation conditions were chosen to be similar to that in Coleman and Bousman34(which is also based on UH-60A flight tests), i.e. an advance ratio of 0.33 at radius r/R=0.865 of the main rotor radius was selected.

    Table 1 Summary of prior work on fluctuating freestream dynamic stall.

    Regarding the angle of attack variation, Kerho20with the help of Bousman35predicted the blade pitching history of UH-60A using the comprehensive rotorcraft analysis code CAMRAD II at advance ratio of 0.33 and radius r/R=0.865. According to Kerho,20a sinusoidal function of 10°±10° can be viewed as an acceptable approximation for the effective angle of attack history and was therefore used in the current study too.According to Table 2,the Mach number fluctuation shall be set to 0.537±0.205. Eqs. (1) and (2)provide the formulation of the angle of attack α oscillation and Mach number Ma fluctuation,while Fig.1 illustrates their resulting combination, ψ is the azimuth angle.

    Table 2 Simulation parameters and conditions.

    Fig. 1 Freestream Mach number fluctuation versus angle of attack oscillation as a function of the rotor azimuth angle.

    where Ω is the rotational frequency.

    Note that there is a phase difference of 180° between the periodic fluctuation of the freestream Mach number and the angle of attack.

    Fig. 2 UH-60A main rotor blade and SC1095 airfoil.36.

    3. Airfoil geometry

    The current study mimics the flight conditions, which an UH-60A main rotor blade section at r/R=0.865 would experience in forward flight at the flight Mach number of 0.205.The UH-60A main rotor blade consist of two airfoils, SC1095 and SC1094 R8. At r/R=0.865, the airfoil is the SC1095, hence,this was selected for this study. Fig. 2 illustrates the UH-60A main rotor blade and SC1095 airfoil.36

    4. Numerical method

    For the size of the computational domain, mesh generation,initialization, verification and validation, the techniques of Al-Jaburi and Feszty37were used in the current study. The novelty is the addition of the fluctuating freestream by setting the freestream at the inlet to Mach number 0.537±0.205. It was anticipated that this case will yield transonic flow at the extreme of the Mach number and thus it offers the possibility to examine the effects of transonic flow on dynamic stall.

    Introducing a fluctuating freestream velocity simultaneously with the airfoils’ oscillation can be achieved either by producing some form of velocity vector changes with time9or by a significantly more complicated procedure.20,23Table 3 summarize the differences between the current method and that used in the literature so far, c is the airfoil chord length.

    The methodology proposed in this work does not involve mesh deformation nor does not need to consider extra measures for the velocity vectors changes, i.e. the current method is simpler and is anticipated to be faster.

    Table 3 CFD literatures for fluctuating freestream dynamic stall in forward flight.

    The ANSYS FLUENT CFD code38was used for all simulations, i.e. for both the constant and fluctuating freestream cases. This discretizes the unsteady Reynolds-Averaged Navier-Stokes equations in the computational domain via a Finite Volume Method (FVM), which has 2nd order accuracy in space.The time derivatives are discretized via Implicit Time Integration method,with 2nd order accuracy in time.Two turbulence models were used in the study;Spalart-Allmaras(SA)and Shear Stress Transport (SST) k-ω model. Although the SST k-ω model is generally more accurate compared than the SA model when simulating constant freestream dynamic stall,39it was found by Richter et al.40Klein et al.41and Al-Jaburi et al.37that it generates two extra peaks in the aerodynamic characteristics (CL, CD, Cm) around the peak angle of attack. These appear to occur due to the SST k-ω model producing non-physical vortices. Therefore, the SA turbulence model was chosen for all constant freestream as well as fluctuating freestream simulations in this study, since it is expected to capture the flow physics more credibly as well as to reduce the computational cost. The one-equation SA turbulence model proved to be effective and robust for a variety of 2D airfoil flows with significant flow separation.21,42,43

    4.1. Computational domain

    The domain and mesh developed by Al-Jaburi et al.37for constant freestream dynamic stall were adopted for the current,fluctuating freestream study. The computational domain design and size justifications were discussed in detail in reference,37here only the key parameters of the domain are discussed. The domain is of rectangular shape, with 1000c by 1000c size and the airfoil placed in the middle. A sliding mesh technique is applied,38and thus a circular domain of 200c diameter is placed inside the rectangular domain, which allowed it to rotate and thus to introduce the angle of attack variation. Fig. 3 provides a sketch to the domain used by Al-Jaburi and Feszty.37

    4.2. Computational domain reconfiguration

    Fig. 3 Computational domain sketch used by Al-Jaburi and Feszty,37 the ‘‘black” boundaries represent pressure far field and the ‘‘gray” boundaries are mesh interface.

    For the present study,the computational domain was reconfigured to enable a fluctuating freestream simulation. Instead of using the pressure far field boundary condition, an inlet boundary condition was introduced at the left side of the domain, pressure outlet for the right side, while the upper and lower sides were set to be free-slip surfaces as shown in Fig.4.This modification would provide a one directional flow for the fluctuating velocity, i.e. no change in the direction of the vectors is needed.

    Fig. 4 Computational domain sketch after reconfiguration, the‘‘black”lines represent inlet and outlet boundaries with inlet to the left, the ‘‘dashed” lines were set to free-slip boundaries and the‘‘gray” boundaries are mesh interface.

    An in-depth analysis to the physics of the fluctuating freestream at the inlet suggests that using the domain37(inlet at 500c or 265 m from the airfoil) is too large because the wave length of the velocity fluctuation is merely a 44.45 m, i.e. for the velocity fluctuation to reach the airfoil,6 cycles are needed,and this would imply too high computational cost. Therefore,based on the mean velocity, the inlet was proposed to be located at only 25% of the wave length (≈18.9c in this case).

    As a result, the sizes of the domain from Ref.37remained intact but the inlet location was reduced to 20c, a distance,which proved to be the minimum sufficient inlet distance for a typical 2D CFD simulation.44

    4.3. Computational time delay

    Because changes at the inlet boundary take time to propagate to the airfoil,the resulting velocity profile at the airfoil location was in a phase delay relative to the inlet. Thompson45,46suggested a solution for time-dependent boundary conditions(the inlet in this case), which is capable of addressing the time-difference behind this phase-delay accurately and it is given as:

    where td(t) is the time delay in seconds at each instantaneous velocity U t( ) and a is the speed of sound. Once the time delay is calculated,the maximum time will be selected(usually at the minimum velocity for a uniform velocity profile) and then it will be implemented in the equations of motion. It was found that adding the time delay to Eq. (2) would increase the complexity of the simulation. Therefore, the best solution was to modify Eq. (1) to accommodate the time shift, leading to:

    Fig.5 provides a detailed presentation for the velocity profiles and the final angle of attack and Mach number at the airfoil location. As seen in the figure, the intended formulation illustrated previously in Fig. 1 was successfully achieved.

    Fig. 5 Mach number at inlet and airfoil location and the final Ma-α formulation at airfoil location after applying time delay.

    4.4. Mesh generation

    The computational mesh was of a hybrid structuredunstructured type, with a total of 252,321 cells. 650 cells were distributed along the surface of the airfoil.A structured rectangular mesh was used in the‘‘inflation layer”around the airfoil,so that the boundary layer can be effectively captured. 40 layers of structured mesh were employed inside the inflation layer,with the first spacing around the wall kept well below y+=1.Therefore, no wall function was used in the turbulence model.An unstructured mesh of 2D triangles was used everywhere outside of the inflation layer. The unstructured mesh was refined in the area above and downstream of the airfoil as shown in Fig. 6, where the transonic shockwave and dynamic stall vortices are expected to appear.

    4.5. Initialization

    The initialization of the transient simulations in this study was accomplished following the technique described in Ref.37, i.e.that first an initial steady-state simulation was generated with the airfoil set to the conditions corresponding to 0° azimuth angle in Fig.5,i.e.the minimum angle of attack and maximum Mach number. The steady-state simulation was run until the residuals of the continuity and energy equations converged to a level of 10-12. Then, the fluctuating freestream boundary condition was introduced in the fashion discussed above.Fig. 7 illustrates the time histories of the aerodynamic loads for 3 cycles. Note that with the help of the steady-state flow initialization, no initial transient is needed, periodicity can be observed directly from the 1st cycle. Thus, following the methodology of this paper, only one cycle is enough to simulate the fluctuating freestream dynamic stall.

    4.6. Verification

    To increase the confidence in the above methodology, it was important to conduct a new series of grid and other simulation parameters dependence studies to see if the modified domain and mesh is appropriate for a fluctuating freestream dynamic stall simulation; this is true because the inlet is now closer to the airfoil (20c) and the radius of the rotating domain was decreased from 100c in Ref.37to only 10c in the current work.Time-step independence was achieved by using 2000 time-steps per pitching cycle and 100 inner iterations with a Courant number 200 and a residual of 10-4between each consecutive time-steps. Fig. 8 shows the results of this grid convergence analysis.

    4.7. Validation

    The validation process was divided into two parts. In the first one, a transonic steady-flow was compared to experiment to validate the ability of the current mesh to capture transonic flow on the advancing helicopter blades featuring a shockwave on the airfoil surface. In the second part the methodology of simulating dynamic stall under fluctuating freestream conditions was compared to an experiment which, however, did not feature transonic flow.

    To the knowledge of the authors,no other paper in the past have presented such validation.

    4.7.1. Transonic validation

    In this case a SC1095 airfoil was subjected to a transonic steady flow of Mach number 0.8 and Reynolds number of 5.65×106. These conditions correspond to the wind tunnel data of Ref.47for angles of attack of 2.1°and 6.2°.The experiment in Ref.47was designed to test full scale rotorcraft airfoils at full scale Reynolds numbers and Mach numbers in the NASA Ames Eleven-Foot Transonic Wind Tunnel. Fig. 9 provides the pressure coefficient Cpat the designated angles of attack. As illustrated in the figure, an excellent agreement was reached using the proposed mesh in this study.Hence,this mesh was considered to be fine enough for simulating fluctuating freestream dynamic stall where shock-induced flow separation might occur.

    4.7.2. Fluctuation freestream validation

    The validation in this section was accomplished with the help of references.24-26,48The experiment involves a SSC-A09 airfoil undergoing a dynamic stall in a fluctuating freestream environment,with the reduced frequency of 0.05 for the simultaneous pitching oscillation of 8.5°±13° and Mach number fluctuation of 0.4±0.08. The airfoils oscillation is given by Eq. (5) while the Mach number fluctuation is approximated by Eq. (6). Note that Eq. (6) is adopted from Ref.24, where Δφ represents the phase difference.

    Fig. 6 Magnified images of mesh used in current study.

    Fig. 7 Time history of aerodynamic loads for fluctuating freestream dynamic stall simulation (SC1095 airfoil, frequency=4.25 Hz,Re=6.1×106, Ma=0.537±0.205 and α=10°±10°).

    Fig. 8 Results of grid dependent study using two different grids (SC1095 airfoil, frequency=4.25 Hz, Re=6.1×106, Ma=0.537±0.205 and α=10°±10°).

    Although this experiment has a very low Mach number amplitude (0.08) compared to that of the current study(0.205),as well as it focused on recording only lift and pitching moment(drag was omitted),it was decided to be used here for validation purposes because so far, this is the only experiment that involves dynamic stall representative of helicopter forward flight as well as at some moderate freestream fluctuation.Also, the frequency of pitch and Mach number fluctuations is representative of helicopters, up to 17 Hz.26However, the experiment still has some considerable limitations, which can be summarized as:

    (1) Pressure tabs limitation: the airfoil model in Ref.26was outfitted with 53 surface pressure taps that covers approximately 80c%of the airfoil from the leading edge,i.e. the rest of the 20c% of the airfoil was left without pressure taps due to geometrical limitation. Hence, and the pressure at the trailing edge vicinity was estimated via a ghost tap in post-processing, which might lead to a possible post-processing error. Moreover, due to the insufficient resolution of the pressure taps near the leading edge,26the accuracy of the results was significantly affected, especially when the leading-edge vortex was convected downstream on the airfoil upper surface with the increase of angle of attack.

    Fig. 9 Transonic flow validation (SC1095 airfoil,Re=6.65×106, Ma=0.8).

    (2) Test section aspect-ratio: the test section used in Ref.26had an aspect-ratio of 1. According to Ref.49, at such low aspect-ratio, the two-dimensionality of the flow was likely not perfect,which could lead to an artificially high velocity in the midspan region because of the end wall boundary layers constricting the flow. Also, the low aspect-ratio was responsible for delaying the onset of dynamic stall and promoting an early dynamic reattachment, weakening the leading-edge vortex. As a result, the pitching moment peak was smaller than expected (lower strength) and the effective angles of attack were changed due to the downwash caused by the end walls of the test section.

    (3) Pitch oscillation assembly vibration: when the experiment of Ref.26was designed, the intended maximum angle of attack was 20°. However, with the unpowered pitching amplitude, and the significant vibration in the oscillation assembly, this made the measured amplitude angle of attack to be increasing with frequency. Therefore,when the pitching amplitude increased,the angular acceleration and the effective frequency at the onset of stall or reattachment were slightly higher than the design values.

    Figs. 10-12 illustrates the results of the entire validation process. As one can see, simulation matches the experimental data quite well, especially when one considers the limitations of the experiment above and the expected overshoots in the peak values of the hysteresis loops,along with the discrepancy in the downstroke phase of moment, which are very common in the literature and is linked to the turbulence model used.

    5. Results and discussion

    CFD simulations of 2D dynamic stall for SC1095 airfoil were completed with the numerical parameters described above.Two simulations were conducted.The first one was a dynamic stall case with a constant freestream velocity corresponding to the mean velocity of the fluctuating freestream case. This will serve as the baseline case,to which fluctuating freestream cases will be compared to.The second case is dynamic stall with fluctuating freestream.The discussion of results below is therefore divided into two sections. However, before presenting the results,a note needs to be made about the scaling factors used to express the coefficients of aerodynamic forces and moments.

    Fig. 10 Simulation methodology validation validated with Ref.26 (Re=3×106, Δφ=13.3°, k=0.05, α=8.5°±13° and Ma=0.4±0.08).

    In the current study,for both the Constant Freestream Simulations(CFS)and Fluctuating Freestream Simulations(FFS),lift,drag and pitching moment were scaled by the mean velocity at each instant in time to express the aerodynamic coefficients. This is in contrast to the literature, where forces and moment of FFS cases are scaled by the instantaneous freestream velocity at each instant in time,for example,Gosselin,9Kerho20Hird24,25Zhao23and Gregory26et al. This type of scaling leads to FFS results to be higher than those under CFS, as illustrated for example in Fig. 13, which reproduces the experimental results of Ref.26. However, this gives a false impression about how the FFS and CFS loads relate to each other in reality, i.e. in dimensional terms. Fig. 14 illustrates that the actual loads,i.e.the dimensional forms of the aerodynamic forces and moment,behave in the opposite manner and that for these the FFS values are actually smaller than the CFS values. The authors of this paper believe that the nondimensionalized aerodynamic loads as well their dimensional equivalents should show the same trends. Therefore, it is proposed in this paper that the aerodynamic loads are scaled by the mean velocity and density, which yield the correct trends between the FFS and CFS results.

    Fig. 11 CFS dynamic stall validated with Ref.26 (SSC-A09 airfoil, Re=3×106, k=0.05, α=8.5°±13° and Ma=0.4).

    A detailed flow visualization of the same case of Ref.26at angle of attack 15.14° is also provided in Fig. 15 to support the fact above. One can clearly see the increased size of separation in FFS compared to that of the CFS, and hence, lift should be smaller in the FFS case compared to the CFS case,not the opposite.

    Fig. 12 FFS dynamic stall validated with Ref.26 (SSC-A09 airfoil, Re=3×106, Δφ=13.3°, k=0.05, α=8.5°±13° and Ma=0.4±0.08).

    This is further proof that the recommended scaling, based on the mean values appears to be correct.This was the method used in the current study.

    5.1. Compressible constant freestream dynamic stall: UH-60A case

    For comparing the CFS and FFS dynamic stall results,as it is intended in this study,one must choose a common mean freestream Mach number for both cases.For this study,this common freestream Mach number was selected by taking the rotational velocity at 86.5%radius of the blade,which according to Table 2 corresponds to Mach number 0.537.Therefore,this will be the freestream for the constant freestream dynamic stall simulation as well as the mean value used for nondimensionalizing the aerodynamic loads in the fluctuating freestream cases too. According to Table 2, the amplitude of the freestream fluctuation will be Mach number 0.205, yielding a fluctuating freestream of 0.537±0.205.

    Fig. 13 CFS (Ma=0.4) vs FFS (Ma=0.4±0.08) dynamic stall of Ref.26 experiment (lift and pitching moment coefficients,SSC-A09 airfoil, Re=3×106, k=0.05, α=8.5°±13°).

    Results for the CFS simulations are shown in Fig.16.Note that the main influence of exposing an airfoil to dynamic stall conditions (in terms of frequency and angle of attack fluctuation)at the constant freestream of Mach number 0.537 leads to the domination of compressibility effects. This is best manifested by the dramatic decrease of the stall angle of attack from 20°(as was seen at Mach number 0.4 in Fig.11)to only about 12°, i.e. lift stall occurs closer to the static stall value, rather than at the usual overshoot associated with a dynamic stall vortex. This is due to the occurrence of transonic flow and shock-induced boundary layer separation on the upper surface well below the ‘‘classical” dynamic stall angle of attack. This means that compressibility effects dominate the stall mechanism, a fact supported also in the literature, for example in references.49,50

    5.2. Fluctuating freestream dynamic stall: UH-60A case

    The FFS simulation conditions of the main rotor blade of UH-60A in a forward flight are sketched in Fig.17.The selected r/R=0.865 section is represented by the dashed-arrows in the velocity distribution plot. Because of the fluctuating relative freestream of 0.537±0.205, the advancing blade will be subject to transonic flow (Ma=0.742) while the retreating blade to dynamic stall(Ma=0.332). A detailed analysis of the flow via flow visualizations is provided in Fig. 18 and according to the selected frames (boxed letters) shown in Fig. 19.

    From Fig.19,there are four substantial differences between the ‘‘real life” dynamic stall at FFS and the ‘‘representative dynamic stall” at CFS:

    (1) Stall in the aerodynamic loads in the FFS case occur earlier than in the CFS case, i.e. at lower angles of attack.

    (2) The peak values at stall are smaller in the FFS case than in the CFS case.

    (3) Before stall, the aerodynamic loads are generally larger in the FFS case than in the CFS case.

    (4) Beyond stall,the aerodynamic loads are generally lower in the FFS case than in the CFS case.

    Fig.14 CFS(Ma=0.4)vs FFS(Ma=0.4±0.08)dynamic stall of Ref.26 experiment(lift and lift coefficient,SSC-A09,Re=3×106,k=0.05, α=8.5°±13°).

    Moreover, due to the compressibility effects in the FFS cases, there are two large peaks during the upstroke phase(Fig.19).The presence of these can be explained by examining the flow visualisation frames on Fig. 18. From these it can be seen that the first peak (Frame C) is generated by the shockformation as the angle of attack increases, while the second one likely because of the stationary shockwave at the leading-edge causing shock wave - boundary layer interaction and a consequent vortex shedding, which forms alternating leading and trailing-edges vortices (shock-induced dynamic stall). This process is visible from Fig. 19, Frames C, D, E,F and G.

    Fig. 15 CFS vs FFS dynamic stall of Ref.26 experiment (White dashed-line represents maximum vorticity, SSC-A09 airfoil,Re=3×106, k=0.05, α=15.14° and Ma=0.4±0.08).

    Fig. 17 Top view sketch of UH-60A main rotor disk showing Mach number distribution along the blade in a forward flight phase according to Table 2.

    In the FFS case,there is a large shock wave at the quarterchord position of the airfoil,starting from the beginning of the upstroke phase (for example, see Frame A, Fig. 18 at α=2°and ψ=36°).This diminishes in strength and moves upstream along the airfoil as the angle of attack increases, until the first leading-edge vortex is shed (Frame C, Fig. 18 at α=10° and ψ=90°). During this process, the flow upstream and downstream of the shock appears to be attached.

    The lift during upstroke phase was higher compared to the CFS case from angle of attack 0° to 7°. This is because of the higher speed on the upper surface of the airfoil during this range.However,lift starts to decrease during the next upstroke phase due to the decrease in speed on the airfoil upper surface.This is true because this will lead to higher pressure on the upper surface compared to that on the lower surface. This can be illustrated by Frame C on Fig.18,where the shockwave is moving upstream during the upstroke phase. The separated flow region starts to increase and becomes large enough to cause the lift coefficient for the FFS case to be less than that for the CFS case, which is logical since separated flow generates less lift than attached flow.This further supports the argument, that the aerodynamic loads shall be nondimensionalized by the mean velocity and not the instantaneous velocity to reflect the flow physics.

    Fig. 16 CFS dynamic stall (SC1095, Re=6.1×106, f=4.25 Hz, α=10°±10° and Ma=0.537).

    Fig. 18 Mach number contours and instantaneous vorticities (SC1095 airfoil, Re=6.1×106, f=4.25 Hz and α=10°±10°).

    Fig. 18 (continued)

    On the other hand, because of the shockwave was present from the beginning of the upstroke phase, between angles of attack 0°to 10°the drag was higher in the FFS case compared to that in the CFS case (Fig. 19). For the rest of the upstroke,the drag in the FFS case becomes lower than that in CFS case.This is likely due to the shockwave moving upstream. As a result, the wave-drag didn’t have enough time to develop into a value that could affect the total drag. Moreover, the instantaneous velocity starts to decrease after Frame C (Mach number decreases) and this leads to a decrease in the speed, and hence, the drag.

    Fig. 18 (continued)

    Regarding the FFS pitching moment,at angle of attack 12°(Frame E in Fig. 19), the pitching moment was remarkably lower in the FFS case than in the CFS case. This behaviour is due to the trailing-edge vortex, which has the effect of decreasing the pitching moment as well as the drag.

    It is noteworthy that in the FFS case, when the angle of attack increases, the leading-edge vortex convected aft to the trailing-edge is shrunk in size and appears to be closer to the airfoil.On the other hand,the trailing-edge vortex starts stronger at angle of attack 12° and becomes weaker as the freestream Mach number decreases during the upstroke phase(see Frames E to H, Fig. 18). Further in the upstroke phase of the FFS case, at angle of attack 20° (Frame H, Fig. 18),the wake of the airfoil features a weaker trailing-edge vortex(divided into two vortices) with a leading-edge vortex significantly lower in its strength compared to that in the CFS case.

    The downstroke phase begins with a leading-edge vortex shedding. As the vortex is shed downstream over the airfoil,a significant enhancement in lift and a reduction in the magnitude of the pitching moment can be observed due to the increasing freestream Mach number and the flow reattachment process. For example, from Fig. 19 Frame I, the pitching moment is higher due to these reasons,and as a result,the drag is significantly reduced compared to that in the CFS case. Also, in Fig. 19 Frame K, although the drag has increased, lift has greatly increased compared to that in CFS case. This is due to the presence of the shockwave that is clearly shown in Fig. 18, Frame K, which it was absence in the CFS case.

    Fig.19 Comparison of aerodynamic loads at FFS dynamic stall(Ma=0.537±0.205)and at CFS dynamic stall(Ma=0.537)with the selected azimuthal positions (SC1095, Re=6.1×106, f=4.25 Hz, and α=10°±10°).

    6. Conclusions

    A methodology for simulating real life 2D shock-induced dynamic stall under fluctuating freestream was described in this paper, which is among the first ones to provide such detailed methodology and validation provided in the literature.The technique is based on a previous study conducted by the same authors for dynamic stall under constant freestream condition, thus generating a ‘‘fluctuating freestream” methodology to accommodate the required unsteady conditions in the freestream.

    From the analysis of the current study,it was found that for both of CFS and FFS, lift, drag and pitching moment should all scaled by the mean velocity and density at each instant in time to drive the correspondent coefficients. The current approach considered more accurate because,the resulted coefficients are in consistent with the unscaled forces and moment.

    It was shown that dynamic stall under FFS significantly differs from the typical dynamic stall at CFS in many ways.Most comprehensive rotor codes are based on using dynamic stall data neglecting the fluctuating freestream and these represent an optimistic approach since these over-predict the loads when compared to that seen in dynamic stall using FFS. In general,FFS is characterized by shock-induced flow separation and as such, stall will occur much earlier than in constant freestream dynamic stall.Also,one can notice that the phenomena on the advancing blades do affect the phenomena on the retreating blades, emphasizing the need to consider fluctuating freestream when transonic flow is achieved on the advancing blades of a helicopter.

    The effect of compressibility is significant when simulating dynamic stall under CFS conditions as this will shift the location of the peak aerodynamic loads to occur before the airfoil static-stall angle of attack, which decreases the overall values of the aerodynamic loads.

    Future works shall involve the examination of these effects in 3D rotor blade simulations. It would be important to conduct an experimental study of the test case presented in this research to validate the conclusions further. However,setting-up an experimental model for the case studied in this work has yet to be carried out due to the practical challenges of reproducing such amplitudes of freestream fluctuation in a wind tunnel.

    大话2 男鬼变身卡| 麻豆成人av视频| 亚洲四区av| 中国国产av一级| 欧美一区二区亚洲| 国内精品宾馆在线| 熟妇人妻不卡中文字幕| 三级国产精品片| 日日啪夜夜爽| 久久久精品免费免费高清| 日本熟妇午夜| 国产黄片美女视频| www.av在线官网国产| 免费av不卡在线播放| 国产乱人视频| 一级毛片 在线播放| 国国产精品蜜臀av免费| 69av精品久久久久久| 欧美zozozo另类| 99久国产av精品国产电影| 韩国高清视频一区二区三区| 国产色爽女视频免费观看| 亚洲人成网站在线观看播放| 大陆偷拍与自拍| 亚洲精品乱码久久久久久按摩| 日本午夜av视频| 国产精品国产三级国产av玫瑰| 听说在线观看完整版免费高清| 欧美日韩国产mv在线观看视频 | 九九久久精品国产亚洲av麻豆| 大码成人一级视频| 一个人观看的视频www高清免费观看| 最近中文字幕高清免费大全6| 国语对白做爰xxxⅹ性视频网站| 少妇人妻精品综合一区二区| 丝袜喷水一区| 三级经典国产精品| 免费看a级黄色片| 久久久成人免费电影| 一级a做视频免费观看| 校园人妻丝袜中文字幕| 在线免费观看不下载黄p国产| 国产久久久一区二区三区| 日韩一区二区视频免费看| www.av在线官网国产| 纵有疾风起免费观看全集完整版| 精品国产三级普通话版| 免费大片18禁| 水蜜桃什么品种好| 国产高清三级在线| 久久久久久久久久成人| 久久久亚洲精品成人影院| 国产一区二区在线观看日韩| 国产高清三级在线| 三级男女做爰猛烈吃奶摸视频| 男的添女的下面高潮视频| 99热国产这里只有精品6| 九色成人免费人妻av| 久久6这里有精品| 欧美xxxx性猛交bbbb| 日本三级黄在线观看| 午夜视频国产福利| 免费观看av网站的网址| 国内揄拍国产精品人妻在线| 舔av片在线| 日韩成人av中文字幕在线观看| 日韩亚洲欧美综合| 欧美高清性xxxxhd video| 人妻 亚洲 视频| 国产黄片美女视频| 少妇被粗大猛烈的视频| 超碰av人人做人人爽久久| 亚洲最大成人手机在线| 涩涩av久久男人的天堂| 成人免费观看视频高清| 亚洲精品日韩av片在线观看| 激情五月婷婷亚洲| 久久久久精品性色| 国产色婷婷99| 免费av不卡在线播放| 2021天堂中文幕一二区在线观| 高清午夜精品一区二区三区| 午夜福利在线在线| 国产精品成人在线| 亚洲内射少妇av| 亚洲精品影视一区二区三区av| 免费在线观看成人毛片| 成人一区二区视频在线观看| av天堂中文字幕网| 五月伊人婷婷丁香| 精品少妇久久久久久888优播| 高清欧美精品videossex| 性色av一级| 日韩电影二区| 天堂中文最新版在线下载 | 国产 一区精品| 久久久国产一区二区| 卡戴珊不雅视频在线播放| 免费黄色在线免费观看| a级一级毛片免费在线观看| 国产老妇女一区| 中文欧美无线码| 一级毛片电影观看| 精品一区二区三区视频在线| 亚洲人成网站高清观看| 插逼视频在线观看| www.色视频.com| 午夜免费鲁丝| 成年免费大片在线观看| 亚洲,一卡二卡三卡| 成年免费大片在线观看| 免费av不卡在线播放| 中国三级夫妇交换| 亚洲欧美一区二区三区国产| 搞女人的毛片| 久久久久久伊人网av| 国产在视频线精品| 国产在线男女| 99热国产这里只有精品6| 美女内射精品一级片tv| 中文字幕久久专区| 久久国产乱子免费精品| 777米奇影视久久| 欧美日韩一区二区视频在线观看视频在线 | 国产免费福利视频在线观看| 九色成人免费人妻av| 少妇人妻精品综合一区二区| 色综合色国产| 亚洲自偷自拍三级| 天堂俺去俺来也www色官网| 亚洲,欧美,日韩| 又大又黄又爽视频免费| 亚洲欧美中文字幕日韩二区| 99re6热这里在线精品视频| 欧美日韩精品成人综合77777| 一个人观看的视频www高清免费观看| 久久精品国产鲁丝片午夜精品| 日本免费在线观看一区| 搡老乐熟女国产| 超碰97精品在线观看| 精品酒店卫生间| 高清日韩中文字幕在线| 少妇人妻精品综合一区二区| 国产黄片美女视频| 男女边吃奶边做爰视频| 亚洲人成网站在线播| 成人亚洲欧美一区二区av| 国产男人的电影天堂91| 国产精品伦人一区二区| 国产有黄有色有爽视频| 噜噜噜噜噜久久久久久91| 超碰97精品在线观看| 网址你懂的国产日韩在线| 18禁裸乳无遮挡动漫免费视频 | 精品久久久久久久久亚洲| 亚洲国产最新在线播放| 久久久欧美国产精品| 天天躁夜夜躁狠狠久久av| 狂野欧美激情性xxxx在线观看| 一级av片app| 成年版毛片免费区| 精品久久久久久久久av| 99热这里只有精品一区| 亚洲激情五月婷婷啪啪| 久久韩国三级中文字幕| 久久久a久久爽久久v久久| 日韩av在线免费看完整版不卡| 18禁在线播放成人免费| 亚洲天堂国产精品一区在线| 国产午夜福利久久久久久| 综合色av麻豆| 99久久人妻综合| 亚洲第一区二区三区不卡| 久久久精品94久久精品| 尤物成人国产欧美一区二区三区| 久久久亚洲精品成人影院| 校园人妻丝袜中文字幕| 国产色婷婷99| 欧美日韩在线观看h| 性插视频无遮挡在线免费观看| 亚洲久久久久久中文字幕| 永久网站在线| 少妇裸体淫交视频免费看高清| 亚洲av在线观看美女高潮| 麻豆久久精品国产亚洲av| 国产高清国产精品国产三级 | 国产成人精品福利久久| 亚洲精品乱久久久久久| 国产一级毛片在线| 在线看a的网站| 免费大片18禁| 大又大粗又爽又黄少妇毛片口| 一区二区三区精品91| 视频中文字幕在线观看| 国产成人精品一,二区| av网站免费在线观看视频| 午夜亚洲福利在线播放| 男女无遮挡免费网站观看| 蜜臀久久99精品久久宅男| 天天躁日日操中文字幕| 成人一区二区视频在线观看| 日韩制服骚丝袜av| 亚洲成人久久爱视频| 免费大片黄手机在线观看| 一本色道久久久久久精品综合| 久久久久国产精品人妻一区二区| 成人亚洲欧美一区二区av| 夫妻午夜视频| 国产大屁股一区二区在线视频| 亚洲av中文av极速乱| 六月丁香七月| 亚洲丝袜综合中文字幕| 日本色播在线视频| 简卡轻食公司| 美女脱内裤让男人舔精品视频| 精品一区在线观看国产| 乱码一卡2卡4卡精品| 亚洲精品乱码久久久久久按摩| 69人妻影院| 在线观看一区二区三区激情| 少妇被粗大猛烈的视频| 麻豆精品久久久久久蜜桃| 99视频精品全部免费 在线| 搡老乐熟女国产| 国产乱人偷精品视频| 久久久亚洲精品成人影院| 毛片一级片免费看久久久久| 麻豆成人av视频| 爱豆传媒免费全集在线观看| 亚洲国产成人一精品久久久| 纵有疾风起免费观看全集完整版| 精品久久久久久久末码| 久久亚洲国产成人精品v| 免费电影在线观看免费观看| 国产精品蜜桃在线观看| freevideosex欧美| 一级毛片我不卡| 国产伦在线观看视频一区| 色吧在线观看| 热re99久久精品国产66热6| 日本-黄色视频高清免费观看| 18禁在线播放成人免费| 青春草视频在线免费观看| 五月伊人婷婷丁香| 色婷婷久久久亚洲欧美| 国产欧美亚洲国产| 搡老乐熟女国产| 国产大屁股一区二区在线视频| av播播在线观看一区| 毛片女人毛片| 亚洲经典国产精华液单| 超碰av人人做人人爽久久| 丝袜喷水一区| 成人毛片a级毛片在线播放| av一本久久久久| 黄片无遮挡物在线观看| 女人久久www免费人成看片| 国产一区亚洲一区在线观看| 亚洲精品成人av观看孕妇| 五月开心婷婷网| 啦啦啦在线观看免费高清www| 麻豆乱淫一区二区| 三级经典国产精品| 一级a做视频免费观看| 建设人人有责人人尽责人人享有的 | 午夜激情久久久久久久| 涩涩av久久男人的天堂| 亚洲欧美成人综合另类久久久| 99热全是精品| 日本-黄色视频高清免费观看| 午夜免费鲁丝| 另类亚洲欧美激情| 一本色道久久久久久精品综合| 极品教师在线视频| 欧美亚洲 丝袜 人妻 在线| 亚洲四区av| 看十八女毛片水多多多| av在线播放精品| 午夜免费观看性视频| 菩萨蛮人人尽说江南好唐韦庄| 热re99久久精品国产66热6| 两个人的视频大全免费| 2021少妇久久久久久久久久久| 五月开心婷婷网| 听说在线观看完整版免费高清| 欧美一级a爱片免费观看看| 在线观看一区二区三区| 国产乱来视频区| 身体一侧抽搐| 婷婷色综合www| 五月开心婷婷网| 日韩强制内射视频| 亚洲在线观看片| 少妇丰满av| 99热这里只有是精品在线观看| 欧美xxⅹ黑人| 能在线免费看毛片的网站| 欧美激情国产日韩精品一区| 18禁在线播放成人免费| 午夜视频国产福利| 欧美国产精品一级二级三级 | 男人狂女人下面高潮的视频| 亚洲精品,欧美精品| 亚洲精品日本国产第一区| 一个人看的www免费观看视频| 青青草视频在线视频观看| 久久精品熟女亚洲av麻豆精品| av免费在线看不卡| 新久久久久国产一级毛片| 欧美bdsm另类| 国产v大片淫在线免费观看| 大话2 男鬼变身卡| 最近手机中文字幕大全| 免费av不卡在线播放| 丰满乱子伦码专区| 中文字幕久久专区| 久久人人爽人人爽人人片va| 国产成人精品一,二区| 国产成人91sexporn| 亚洲国产成人一精品久久久| 91aial.com中文字幕在线观看| 简卡轻食公司| 国产精品一区二区性色av| videossex国产| 女人久久www免费人成看片| 观看美女的网站| 日韩在线高清观看一区二区三区| 日韩强制内射视频| 美女被艹到高潮喷水动态| 亚洲一级一片aⅴ在线观看| 久热久热在线精品观看| 免费看a级黄色片| 久久久久国产精品人妻一区二区| 亚洲精品视频女| 在线 av 中文字幕| 永久免费av网站大全| 精品久久久精品久久久| 国产精品一区二区在线观看99| 嫩草影院新地址| 麻豆成人av视频| 我要看日韩黄色一级片| 看黄色毛片网站| 国产亚洲午夜精品一区二区久久 | 中国国产av一级| 精品久久久久久久人妻蜜臀av| 久久人人爽人人爽人人片va| 精品久久久噜噜| 亚洲精品乱久久久久久| 亚洲天堂国产精品一区在线| 九九久久精品国产亚洲av麻豆| 国产91av在线免费观看| 免费播放大片免费观看视频在线观看| 欧美日韩亚洲高清精品| 免费观看a级毛片全部| 成人鲁丝片一二三区免费| 亚洲真实伦在线观看| 久久久久国产精品人妻一区二区| 日日啪夜夜撸| 久久精品久久久久久久性| 国产精品国产三级专区第一集| 乱系列少妇在线播放| 麻豆久久精品国产亚洲av| 亚洲av免费高清在线观看| 欧美日韩一区二区视频在线观看视频在线 | 国产乱来视频区| 亚洲av国产av综合av卡| 麻豆精品久久久久久蜜桃| 国产乱人视频| 在线免费观看不下载黄p国产| 久久精品人妻少妇| 久久人人爽人人片av| 亚洲精品日韩在线中文字幕| 国产乱来视频区| 日产精品乱码卡一卡2卡三| 一级黄片播放器| 午夜福利高清视频| 色播亚洲综合网| 老司机影院成人| 日本av手机在线免费观看| 最近2019中文字幕mv第一页| 国产成人精品久久久久久| 国产免费视频播放在线视频| 国产精品人妻久久久久久| 国产精品一及| 超碰97精品在线观看| 纵有疾风起免费观看全集完整版| 最近中文字幕高清免费大全6| 人人妻人人看人人澡| 又粗又硬又长又爽又黄的视频| 国产大屁股一区二区在线视频| 国产免费一级a男人的天堂| 国产av国产精品国产| 欧美变态另类bdsm刘玥| 99久久中文字幕三级久久日本| 亚洲最大成人手机在线| 尾随美女入室| 亚洲自偷自拍三级| 99久久精品一区二区三区| 丰满乱子伦码专区| 亚洲,欧美,日韩| 免费观看性生交大片5| 欧美激情久久久久久爽电影| 成人无遮挡网站| 免费观看在线日韩| 水蜜桃什么品种好| 国产黄频视频在线观看| 香蕉精品网在线| 精华霜和精华液先用哪个| 国产爽快片一区二区三区| 精品久久久久久电影网| 国产免费又黄又爽又色| 久久ye,这里只有精品| 人人妻人人爽人人添夜夜欢视频 | 欧美日韩综合久久久久久| 亚洲综合精品二区| 制服丝袜香蕉在线| 狂野欧美白嫩少妇大欣赏| 汤姆久久久久久久影院中文字幕| 欧美成人精品欧美一级黄| 搡女人真爽免费视频火全软件| 涩涩av久久男人的天堂| 九色成人免费人妻av| 国产成人午夜福利电影在线观看| 最近中文字幕高清免费大全6| 大香蕉久久网| 一级毛片我不卡| 国内精品美女久久久久久| 亚洲欧美一区二区三区国产| 亚洲熟女精品中文字幕| 精品亚洲乱码少妇综合久久| 777米奇影视久久| 波多野结衣巨乳人妻| 免费大片18禁| 欧美最新免费一区二区三区| 99久国产av精品国产电影| 天天躁夜夜躁狠狠久久av| 韩国高清视频一区二区三区| 日韩中字成人| 国产色爽女视频免费观看| 色哟哟·www| 日韩在线高清观看一区二区三区| 免费观看无遮挡的男女| 免费观看在线日韩| 一级黄片播放器| 精华霜和精华液先用哪个| 久久精品国产鲁丝片午夜精品| av又黄又爽大尺度在线免费看| 80岁老熟妇乱子伦牲交| 91在线精品国自产拍蜜月| 国产探花极品一区二区| 久久人人爽av亚洲精品天堂 | 夜夜看夜夜爽夜夜摸| 一本久久精品| .国产精品久久| 卡戴珊不雅视频在线播放| 国产av码专区亚洲av| 国产精品久久久久久精品电影小说 | 三级国产精品片| 少妇人妻一区二区三区视频| 国产日韩欧美亚洲二区| 男女边吃奶边做爰视频| 亚洲国产日韩一区二区| 免费观看的影片在线观看| 91午夜精品亚洲一区二区三区| 亚洲精品456在线播放app| 久久久成人免费电影| 联通29元200g的流量卡| 视频中文字幕在线观看| 九九在线视频观看精品| 午夜激情福利司机影院| 少妇丰满av| 国产日韩欧美在线精品| 欧美性感艳星| 99精国产麻豆久久婷婷| 国产黄色视频一区二区在线观看| 国产 一区精品| 中文字幕人妻熟人妻熟丝袜美| 午夜激情久久久久久久| 在线精品无人区一区二区三 | 中文在线观看免费www的网站| 啦啦啦在线观看免费高清www| 91精品一卡2卡3卡4卡| 69av精品久久久久久| 欧美xxⅹ黑人| 国产69精品久久久久777片| 九九久久精品国产亚洲av麻豆| 精品视频人人做人人爽| 69人妻影院| 亚洲电影在线观看av| 久久人人爽人人片av| 国内揄拍国产精品人妻在线| 天堂网av新在线| 国产精品女同一区二区软件| 亚洲自偷自拍三级| 亚洲无线观看免费| 国语对白做爰xxxⅹ性视频网站| 在线观看av片永久免费下载| 亚洲精品国产av蜜桃| 亚洲国产欧美在线一区| 欧美激情国产日韩精品一区| 日本一二三区视频观看| 国产黄频视频在线观看| 最近最新中文字幕大全电影3| 国产视频首页在线观看| 看十八女毛片水多多多| 精品人妻熟女av久视频| 看非洲黑人一级黄片| 欧美一级a爱片免费观看看| 免费大片18禁| 精品熟女少妇av免费看| 亚洲精品中文字幕在线视频 | 欧美xxⅹ黑人| 成人二区视频| 亚洲精品色激情综合| 亚洲性久久影院| 日韩免费高清中文字幕av| 亚洲精品久久久久久婷婷小说| 别揉我奶头 嗯啊视频| 久久久久久国产a免费观看| 国产精品一区二区三区四区免费观看| 精品一区二区三区视频在线| 精品久久久久久久久亚洲| 精品一区二区三卡| 亚洲色图综合在线观看| 日韩欧美 国产精品| 91精品国产九色| 自拍欧美九色日韩亚洲蝌蚪91 | 国产在线男女| 嫩草影院入口| 婷婷色麻豆天堂久久| 国产在线一区二区三区精| 男女国产视频网站| 久久久久久久久久人人人人人人| 久久精品国产亚洲网站| 我的女老师完整版在线观看| 日日摸夜夜添夜夜添av毛片| 亚洲精品国产av成人精品| 亚洲精品,欧美精品| 最近最新中文字幕大全电影3| 欧美性感艳星| 欧美xxⅹ黑人| 国产精品蜜桃在线观看| 国产男女内射视频| 网址你懂的国产日韩在线| 狂野欧美白嫩少妇大欣赏| 日韩制服骚丝袜av| 欧美老熟妇乱子伦牲交| 男插女下体视频免费在线播放| 黄片无遮挡物在线观看| 成人美女网站在线观看视频| 国产精品精品国产色婷婷| 中文资源天堂在线| 亚洲国产精品成人综合色| 制服丝袜香蕉在线| 18禁动态无遮挡网站| 色哟哟·www| .国产精品久久| 国产乱来视频区| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 99热网站在线观看| 亚洲国产精品999| 亚洲四区av| 国产一区二区在线观看日韩| 国产伦精品一区二区三区四那| 嫩草影院入口| 国产精品偷伦视频观看了| 99热这里只有是精品50| 久久精品久久精品一区二区三区| 男的添女的下面高潮视频| av线在线观看网站| 久久久久久国产a免费观看| 国产高清不卡午夜福利| 久久精品久久久久久久性| 中文精品一卡2卡3卡4更新| www.av在线官网国产| 狂野欧美白嫩少妇大欣赏| 一级毛片 在线播放| 国产成人午夜福利电影在线观看| 久久久久九九精品影院| 国产午夜精品久久久久久一区二区三区| av线在线观看网站| 久久人人爽av亚洲精品天堂 | 久久人人爽人人爽人人片va| 白带黄色成豆腐渣| 久久久久性生活片| 五月玫瑰六月丁香| 白带黄色成豆腐渣| 三级国产精品片| av又黄又爽大尺度在线免费看| 国产精品一区二区三区四区免费观看| 狂野欧美激情性bbbbbb| 中国三级夫妇交换| 欧美激情久久久久久爽电影| 99热6这里只有精品| 日韩三级伦理在线观看| 精品酒店卫生间| 韩国av在线不卡| 大码成人一级视频| 国产白丝娇喘喷水9色精品| 青春草国产在线视频| 亚洲精品乱久久久久久| 大香蕉久久网| 麻豆精品久久久久久蜜桃| 精品午夜福利在线看| 一个人看的www免费观看视频| 性色av一级| 精品久久国产蜜桃| 插逼视频在线观看| 成人亚洲精品av一区二区| 嫩草影院精品99| 两个人的视频大全免费| 亚洲精品国产av蜜桃| 久久久久久国产a免费观看| 国产精品久久久久久久电影| 亚洲精品国产av蜜桃|