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

    Aerodynamic performance of hovering micro revolving wings in ground and ceiling effects at low Reynolds number

    2023-02-09 08:59:06JinjingHAOYanlaiZHANGChaoZHOUSongtaoCHUJianghaoWU
    CHINESE JOURNAL OF AERONAUTICS 2023年1期

    Jinjing HAO, Yanlai ZHANG, Chao ZHOU, Songtao CHU, Jianghao WU

    School of Transportation Science and Engineering, Beihang University, Beijing 100083, China

    KEYWORDS Ceiling effect;Flapping wing;Ground effect;Micro air vehicle;Revolving wing

    Abstract Numerous investigations have been conducted to understand the wall effects on rotors.The purpose of this study is to further investigate the aerodynamic performance of revolving wings,especially when it is very close to the ground and ceiling (i.e., less than half the wingspan) at low Reynolds numbers. Hence, the ground and ceiling effect for hovering micro revolving wings at low Reynolds numbers are investigated by improving the theoretical models.The theoretical model for the ground effect is established based on the wall-jet assumption,and that for the ceiling effect is improved by considering the uneven spanwise distribution of induced velocity. These two models are validated by comparing the results of experiments and CFD simulations with the Lattice-Boltzmann Method(LBM).Both ground and ceiling effects are found helpful to enhance the thrust,especially with small wing-wall distances, by making a difference to the induced velocity and the pressure distribution. By comparing the thrust generation and aerodynamic efficiency between the ground and ceiling effects, the former is found more helpful to the thrust augmentation, and the latter is more beneficial for the aerodynamic efficiency promotion.

    1. Introduction

    In the past two decades, Micro Air Vehicles (MAVs) have attracted great research interests due to their potential in military,security,and civil applications,including but not limited to intelligence gathering,fixed-point patrol,precise attack,and disaster relief.MAV prototypes in rotary-wing(e.g.,the Black Hornet1equipped by British scouts in 2016, DJI PHANTOM series,2and Parrot Bebop3) and flapping-wing (e.g., the RoboBee,4,5DelFly Nimble,6and KUBettle7,8) configurations especially stand out because of their hovering capability.These small and light-weighted flying robots are desirable to perform missions in confined environments, confronted with energy savings and reliable operation challenges.

    One of the prominent problems in the application of MAVs is the short mission time.9,10The shrinking size inevitably leads to increased viscous loss and reduced efficiency conversion at low Reynolds numbers.10,11Flying creatures such as birds and insects take advantage of the ground effect to save energy.Researchers exploit this bioinspired method to extend the endurance of MAVs by perching on or hovering near surfaces.12As a result, studies on the influence of interaction between wings (including rotary and flapping wings) and surrounding walls on the aerodynamic performance of MAVs emerge in large numbers.

    It is commonly accepted that the effect of lateral walls is negligible,13,14and lots of attention are paid to the ground and ceiling effects. Jardin et al.13extracted aerodynamic loads and velocity fields by experiments to check the influence of ground and ceiling effects on the aerodynamic performance of a hovering micro-rotor. In their experiments,the disk loading and pitching angle of the rotors were kept unchanged. In this condition, the rotor power loading will rapidly increase when the wing-wall distance is within one rotor diameter,indicating an increase in thrust at a given power.Furthermore,the increase in power loading in the ground effect is more obvious than that in the ceiling effect. They associated the ground effect with the change in the axial induced velocity at the rotor disk, and the ceiling effect with the distortion of the flow toward the rotor. Prothin et al.15carried out a numerical and experimental analysis on the flow field around a revolving wing in the ground effect.They validated that a decreasing distance between the rotor and the ground makes the flow recirculate upward and impinge the rotor from underneath, which changes the force distribution on the blade over the wingspan.Kocer et al.16tested the flight of micro-rotors on close proximities of ceilings and collected the real-time flight data,including rotor speed, rotor thrust, and battery current. The results showed a decrease in the current drawn by the battery with the ceiling effect,suggesting that the thrust efficiency increases while the rotor approaches the ceiling.

    Apart from revolving wings, flapping wings in the wall effects also attracted much attention. Lu et al.17-19proposed that the ground effect on flapping wings consists of ‘‘force enhancement-force reduction-force recovery”, unlike rotary wings. They claimed that the force enhancement is because the ground blocks the downward convection of the trailingedge vortex shed from the previous stroke;the reduction might be due to the induced jet flow pointing toward the leeward side of the wing; the recovery is associated with the leading-edge vortex reattachment.However,in the ceiling effect,the aerodynamic force of flapping wings is found to increase monotonically with a closer position to the ceiling.14,20Meng20interpreted this phenomenon as the increasing relative velocity and the pitching angle, together with an interaction between the leading-edge vortex and its mirror vortex actuated by the ceiling.

    Experimental evidence shows that ground and ceiling effects are not avoidable in the operation of MAVs, such as takeoff and landing, and must be treated with caution.21,22For example, MAVs are very likely to hit the ceiling due to the suction of the additional thrust, resulting in a crash without proper control.16,23Nakata et al.24developed a surface detector for autonomous flying vehicles. This collisionavoidance mechanism is independent of vision, mediated by perceiving modulations of induced flow caused by wall effects.Because control strategies accommodating to realistic situations call for accurate quantitative data about the change in aerodynamic force and power, it is ideal to have a theoretical model to assist design.

    Most existing analytical models toward wall effects are based on the situation where a helicopter rotor is in the ground proximity,dating back to the early 20th century.In the 1950s,Zbrozek25implemented experiments on a hovering helicopter and determined that the thrust of a helicopter rotor operating at constant power increases as it approaches the ground.Betz26then built a theoretical model in which the rotor is replaced by a sink and the ground by an image sink with equal strength situated below the revolving wing. However, Cheeseman and Bennett27later argued that the flow pattern beneath the rotor resembles a source rather than a sink.Following this new assumption, they deduced an empirical evaluation of the ground effect for hovering revolving wings, i.e., the renowned Cheeseman-Bennett (C-B) model. Although the C-B model agrees well with flight test results, it fails when Dg/R ≤0.25(Dgrepresents the distance between the wing and the ground,and R represents the radius of the rotor) because the thrust is predicted to increase exponentially to infinity within this range.Apart from the C-B model, there are also some theoretical models, such as Hayden’s model and Johnson’s model(Table 127-29). However, these models hugely overvalued the thrust when wings are extremely close to the ground (Dg/R ≤0.25). As for the ceiling effect, helicopters do not have to bother with this issue,and related models are scarce.Conyers et al.30used the C-B model to predict the thrust of rotors in the ceiling effect. It turns out that the predicted results agree with the experimental results in the case of Dc/R >1 (Dcrepresents the distance between the wing and the ceiling) but diverge obviously when Dc/R ≤1. In doing so, the limitation of the C-B model mentioned above cannot be avoided, and the distinct physical mechanisms behind the ground and ceiling effects are ignored. Afterward, Hsiao and Chirarattananon9further proposed an original model for especially dealing with the ceiling effect based on the blade element momentum theory, which can provide good descriptions of experimental observations even when Dc/R ≤0.25.But this model neglected the uneven distribution of the induced velocity along the wingspan and confused the differential thrust for a blade element with the integral thrust over the entire wing.

    In summary, rotary-wing or flapping-wing MAVs could benefit from reduced energy consumption and extended flight endurance with a closer position to the ground and the ceiling,such as perching on surfaces or hovering in the close vicinity of walls. To remain stable in these situations, the flight control system must adjust the altitude and attitude of MAVs. As a result, it is necessary to acquire a priori information about how MAVs perform with the presence of walls. Because the existing models are inadequate to predict change in aerody-namic performance with small wing-wall distances, this paper first modifies those aforementioned theoretical aerodynamic models by: (A) introducing the wall-jet assumption into the ground effect model to deal with small ground clearances;(B) improving Hsiao’s ceiling effect model by considering the uneven induced velocity distribution along the wingspan, and using the correct form of integral thrust over the entire wing.Then, modified models are compared with experimental measurements and CFD calculations. Finally, comparisons between the ground and ceiling effects in terms of thrust,aerodynamic power, and aerodynamic efficiency are conducted to analyze how to maximize the advantages of these two effects.

    Table 1 Theoretical aerodynamic models for a single rotarywing in ground effect.

    2. Model and methods

    2.1. Model

    This study makes use of a simple rectangular wing planform.The wing section is a flat plate, the thickness of which is 1%of the chord length of the wing.The constant chord length(c = 0.08 m) and changing wingspan (R) lead to changing aspect ratio(AR =R/c). AR = 1,3,5, and 7 are considered,basically covering the aspect ratio range of wings in MAVs.

    The global coordinate system O-XYZ and the local coordinate frame o-xyz are introduced to describe the 3D motion and forces of a wing, as shown in Fig. 1 (a) and 1(b). These two coordinate systems share the same origin located on the 1/2 chord position of the wing root. The OY-axis, that is, the axis of revolution, is perpendicular to the horizontal plane O-XZ.The coordinate system o-xyz is fixed on the wing. The oxaxis is spanwise and points from the root to the tip. The oyaxis is chordwise and points from the trailing edge to the leading edge.The oz-axis is determined according to the right-hand rule.

    The thrust (T), aerodynamic power consumption (P), and aerodynamic efficiency (η)are calculated to quantify the aerodynamic performance of the revolving wing. As plotted in Fig.1,the thrust is defined as the projection of resultant aerodynamic force along the OY-axis,and it can be nondimensionalized as CT=2 T/(ρUr2efS);the horizontal force(H)is defined as the projection of resultant aerodynamic force in the O-XZ plane,and it can be nondimensionalized as CH=2H/(ρUr2efS).Aerodynamic power is consumed to overcome the horizontal force, defined as P = Qa·˙φ, where Qais the aerodynamic torque around the OY-axis, and the power consumption can be nondimensionalized as CP= 2P/(ρUr3efS). The efficiency of thrust generation is described as η = CT3/2/CP.35

    Wing-wall distance (Di(i = g, c)) ranging from 0 to 2R is considered, where the ground and the ceiling are specified by the subscript ‘‘g” and ‘‘c”, respectively. In particular, CTand CPof cases with Di= c, 0.5c, and 0.25c are computed using the CFD method to quantitively validate the theoretical model when the revolving wing is extremely close to the wall.

    2.2. Ground effect modeling

    This section establishes a theoretical model to explain the ground effect on the aerodynamics of a revolving wing in hovering. When the wing is far from the ground, that is, Out of the Ground Effect (OGE), the upstream and downstream flows can be regarded as vertical. When the wing is In the Ground Effect(IGE),the wingtip vortex(TV)first sheds from the trailing edge, then shrinks along the radius direction, and extends again when approaching the ground. In this process,the wake of the wing transits from the vertical direction to the radial direction. Meanwhile, the wing TV gets stretched due to dissipation in the presence of the ground. As a result,Fig.1 Schematics of wing motion and aerodynamic forces of a revolving wing in ground effect and in ceiling effect,and an airfoil section of revolving wing at a location r from root.the radial velocity in the wake first increases and then decreases sharply as the wing-ground distance (Dg) decreases,behaving like a wall-jet.36To describe the feature of the wake structure and the radial velocity, the wall-jet assumption is introduced into the theoretical model to ensure applicability in small ground clearances.

    Fig. 2 Schematic of airflow around a revolving wing in ground effect.

    As shown in Fig. 2, it is reasonable to build the hypothesis that the air vertically passes through the wing and generates the axial induced velocity vi. Assuming that viis radially uniform, which yields pressure difference between the upper surface (p-) and the lower surface (p+) of the wing, we have.

    where pr(r) is the pressure of the air near the ground, p0is the atmospheric pressure, and v∞is the velocity of the downwash far from the wing.

    Combining Eq. (1) and Eq. (4), we have.

    Eqs. (3)-(5) associate Δp with vr(r). For simplicity, higherorder terms caused by a6, a5,and a0are omitted due to a negligible contribution to velocity vr(r).The dynamic pressure can be expressed as.

    where CLαand CDαrepresent the slope of the lift and drag curve, respectively; Nbrepresents the number of revolving wings; θ represents the pitch angle (Fig. 1 (c)). CLαand CDαare determined by the quasi-steady model for flapping wings39to address for the cases with high angles of attack in this paper:

    T and P can be obtained by substituting Eq. (17) into Eq.(14) and Eq. (15), respectively.

    2.3. Ceiling effect modeling

    In this section,a theoretical model proposed by Hsiao et al.9to explain the ceiling effect on the aerodynamics of a revolving wing in hovering is improved by considering an uneven distribution of vialong the wingspan. When the wing is ICE, the ceiling mainly obstructs the incoming air instead of the downwash. Therefore, the wall-jet hypothesis in the ground effect model does not work. It is assumed that the incoming airflow horizontally enters the region between the ceiling and the wing,and vertically passes through the wing, yielding axial vi. The horizontal velocity vr(r) is regarded as a function of the radial position r, independent of axial distance from the ceiling.

    Let’s consider a blade element with a width dr at the radial position r, and a tube that connects the horizontal incoming flow beneath the ceiling and the downwash far from the wing(Fig. 3). Along the tube, the flow rate must remain constant:

    After passing through the wing, the airflow continues to develop and finally becomes uniform. According to the Bernoulli equation, the pressure of the air near the wing can be described as.

    Fig. 3 Schematic of airflow around a revolving wing in ceiling effect.

    In particular, when the wing is far from the ceiling, that is,out of the ceiling effect(OCE),Dc→∞,γ(r )→1,Eq.(28)can be further expressed as.

    which is exactly the same as the expression for the induced velocity of hovering revolving wings without wall effects.40

    The above-mentioned analysis of the theoretical model for the ceiling effect is mainly based on the hypothesis and method proposed by Hsiao et al.9Compared with their model, an uneven distribution of vialong the wingspan (Eq. (28)) is considered in this paper, aiming at providing a more accurate description of airflow passing through the wing.

    Substitute Eq.(29)into Eq.(10)and Eq.(11),and T and P over the entire wing are the integral of them along the wingspan, respectively.

    2.4. CFD method

    The validation of theoretical models concerning small wingwall distances,and the comparison of the aerodynamic performance between IGE and ICE wings are conducted based on the CFD results. In this study, the CFD simulation is performed using a Lattice-Boltzmann Method (LBM), whose capability to solve flows for revolving wings at low Re has already been demonstrated.41,42The LBM solves the discrete Boltzmann equation, a statistical equation for the kinetics of gas molecules, instead of directly solving the Navier-Stokes equations. Compared to traditional CFD methods solving N-S equations, in the LBM, fluid particles perform consecutive propagation and collision processes over a discrete lattice Cartesian grid. A Bhatnagar Gross and Krook (BGK) model with the single-relaxation-time approximation is introduced into the collision step for relaxation to equilibrium, and in the BGK model,the probability fi(x,t)with velocity eiat location × and time t is described as.

    where ωiis the weight factor,and csis the speed of sound.The fluid density ρfand velocity u are computed according to the D3Q19 model and defined as.

    The relaxation time τ is related to the viscosity by τ = 3ν + 0.5, where ν is the kinematic viscosity measured in lattice units. To ensure a non-negative kinematic viscosity, τ should be larger than 0.5. Thus, stability problems arise as the relaxation time approaches this limiting value at high Re.43To overcome this problem, the Recursive Regularized BGK (RR-BGK) collision model44is introduced to ensure a special regularization procedure. To treat the boundary condition on the wings, the immersed boundary method45is applied. The body forces are applied on lattice points near the boundary to enforce the no-slip condition on the boundary. The simulation in this study is conducted using LES with Wall-Adapting Local Eddy(WALE) model, which is assumed to have good properties both near to and far from the wall and for both laminar and turbulent flows.

    Convergence checks are performed using the AR=5 wing far from the ground or ceiling as a typical case (c = 0.08 m,AR=5, ˙φ=1.5π rad/s,α=45°,Re≈6000).Results of three grids (Grids 1-3 in Table 2) with different lattice spacings but the same computational domain size (30c) are shown in Table 2. The time-averaged aerodynamic forces from Grid 2 and Grid 3 agree with each other approximately after five revolving cycles,indicating that the computed results are independent of the lattice spacing.Fig.4 shows that there are some obvious differences between Grid 1 and Grid 3 in the vorticity and iso-vorticity surface, but the discrepancy between Grid 2 and Grid 3 is slight. Because the computational time for Grid 3 is much longer than that for Grid 2,Grid 2 is ideal.Similarly,the influence of the computational domain size (20c, 30c, 40c)is also checked separately. The averaged thrust coefficients of Grid 2 with 20c, 30c, and 40c domain size are 1.4440,1.4217, and 1.4189, respectively. The difference in thrust between 30c and 40c domain size is around 0.2%,but the full domain elements of the latter are four times more than the former,leading to a longer calculation time.Specifically,the computational domain of 30c, the finest lattice spacing of 0.01c around the wing, and the coarsest lattice spacing of 0.25c are chosen.The time step is indirectly determined,keeping the stability parameter U(Δt/Δx)(where U is the velocity magnitude,Δt is the time step,and Δx is the finest lattice spacing)between 0.1 and 0.3.

    Table 2 Lattice spacing and calculation time of three different grids (computational domain size of 30c).

    Fig.4 Vorticity fields around the wing with different grids at φ=360°:(upper) contours of vorticity at R2 position (the magnitude of vorticity at outer contour is 500 and contour interval is 20);(down)iso-surface of vorticity(the magnitude of non-dimensional vorticity is 100).

    A quantitative comparison with the existing results in Ref.34 is made.The contrast of the averaged force coefficients is shown in Fig. 5, where CFD results agree with the experimental measurements in the literature. Therefore, it is reasonable to believe that the CFD method is suitable for this study.

    3. Results and discussion

    3.1. Validation of theoretical model

    Fig.5 Contrast of averaged CT and CH between CFD results in this study and experimental measurements in Ref.34.

    Firstly, the accuracy of theoretical models in Section 2 is validated by comparing their estimations with experimental data.For the ground effect, experimental data are taken from Light’s experiment,46in which the thrust of a four-bladed rotor with constant chord and untwisted blades is measured(Fig.6).Parameters in Light’s experiment are shown as follows:R = 1.105 m, c = 0.18 m, ˙φ ≈172.3 rad/s, θ = 13°-17°,Dg= 0.25R-2R, and corresponding Re≈7 × 106. In their results, the ratio of IGE thrust (TIGE) to OGE thrust (TOGE),i.e.,TIGE/TOGE,is used to describe the enhancement or degradation of the thrust with the presence of the ground.13Accordingly,the curve-fitted relationship between TIGE/TOGEand Dg/R is presented, and the analytical results based on the C-B model,27Hayden’s model,29and Johnson’s model30are also given for comparison. As shown in Fig. 6, all the five models show a consistent tendency with the experimental data, and the thrust gets augmented with the decrease of Dg. Quantitatively, predictions agree well with measurements when Dg/R ≥1, except Johnson’s low loading model. However, when Dg/R <1, only the IGE model in this paper can match the experimental results, while others all overestimate the thrust augmentation caused by the ground effect.

    Similarly, for the ceiling effect, measurements taken from Conyers’ experiment28are compared with model estimation in terms of the ratio of ICE thrust(TICE)to OCE thrust(TOCE)with respect to Dc/R (Fig. 6). The Conyers’ experiment used the Gemfan 9 × 4.7 propellers with R = 0.11 m and ˙φ-≈ 1000 rad/s, where Re ≈ 1 × 105. In the range of 0.25 <Dc/R <2,the ICE model in this paper agrees well with the experimental results. Overall,the theoretical model in Section 2 effectively solves the problem of inaccurate prediction of thrust augmentation in the ground and ceiling effects with relatively small clearances,particularly when Di/R ≤0.5,and also performs well with Diranging from R to 2R.

    Fig. 6 Contrast between experimental measurements and model estimations.

    Table 3 Comparison of thrust between theoretical model estimations and CFD calculations.

    Table 4 Comparison of aerodynamic power between theoretical model estimations and CFD calculations.

    Besides, theoretical estimations are also compared with CFD results to further validate the model in this paper with relatively small wing-wall distances and low Re (Re≈1000-9000). TIGE/TOGEand TICE/TOCE, as well as PIGE/POGEand PICE/POCE, of the wing with different AR (AR = 1,3,5, and 7)at varying Di(Di=0.25c,0.5c,and c)from theoretical estimations and CFD calculations are listed in Table 3 and Table 4. The model in this paper gives predictions close to CFD results,either IGE or ICE cases.For the thrust,the maximum errors between the estimation and calculation are - 6.45 % IGE (AR = 1, Dg= 0.5c) and - 10.52 % ICE(AR=7,Dc=0.5c),respectively;for the aerodynamic power,the numbers are - 6.02 % IGE (AR = 7, Dg= 0.25c)and-8.72%ICE(AR=7,Dc=0.5c),respectively.It seems that the effect of the aspect ratio cannot be ignored when it comes to the applicability of our theoretical model. In the range of AR = 3-5, the model performs well in the whole range of 0 <Di<2R, and the maximum error is no more than 4 %; when AR ≤1, predictions are relatively accurate when 0.5R <Di<2R, but gradually deviate from calculations with smaller Di, especially for ICE cases; when AR ≥7, the error gets even larger, and the model tends to underestimate the increase in thrust and aerodynamic power.As for the effect of Re, our improved model runs well compared with experiments with Re exceeding 105and calculations with Re around 103.Overall,the improved theoretical model is applicable within the range of Re≈103-106,covering the operation of most rotary-wing MAVs.

    Fig.7 Contrast of time-averaged flow field and iso-surface of Q criterion around the wing at φ=1800°(AR =1, Di=0.5c)(arrows represent velocity vector, and color represents the velocity magnitude, the radial velocity and distance from the wall are nondimensionalized by Uref and R, respectively).

    The flow field around the wing is depicted to validate the flow assumptions,which are the basis of the theoretical model.Fig.7 shows the velocity field and iso-surface of vorticity based on Q-criterion(Q=20)around the AR=1 wing in three situations: far from the ground and ceiling (OG(C)E), in the ground effect, and in the ceiling effect. For the OG(C)E wing,LEV migrates to the wingtip and then sheds into the wake together with a strong TV. The wake gradually develops into a spiral path below the wing surface. The spiral wake is obvious within the range of 2R apart from the wing along the axis.

    As shown in Fig. 7(b), in contrast, the ground forced the wake integrated with LEV and TV to expand within a small axial distance,behaving like a wall jet.This jet-like flow structure is consistent with our assumption in the theoretical model,assuming that the radial velocity in the wake first increases and then decreases sharply in the process of approaching the ground (Fig. 7 (d)). Meanwhile, the loop-shaped wake rebounds upward from the ground. Due to the reboundedflow, on one hand, the spiral wake is of larger diameter; on the other hand, the vortices accumulate underneath the wing root and resemble a fountain, resulting in additional lifting force acting on the wing.47,48Simultaneously,the flow affected by the ‘‘fountain effect” is too complex to simplify as a walljet, which may be one of the reasons for the discrepancy between the model estimation and CFD calculation.

    Fig. 8 Contrast of CFD-calculated and model-predicted of wings with different AR at different Di.

    As shown in Fig.7(c),For the ICE wing,the incoming flow is approximately horizontal above the wing due to the obstruction of the ceiling, consistent with the flow assumption in the theoretical model.Influenced by the horizontal flow,the spanwise and tangential flow carries the LEV and TV away from the wing surface with higher velocity, collecting near the leading edge. The spiral wake in the stable spiral structure of the wake cannot be kept; instead, it quickly dissipates and breaks into many pieces within the range of R apart from the wing along the axis.

    3.2. Effects of Di on aerodynamic performance

    This part discusses the influence of Dion the aerodynamic performance of IGE and ICE wings based on CFD results. The contrast of CFD-calculated and model-predicted TIGE/TOGEand TICE/TOCEof wings with different AR at varying Diare presented in Fig. 8 (a) and 8(b). With a specific AR, either TIGEor TICEkeeps increasing as Didecreases; the smaller Diis, the higher the increasing rate is. Within the range of Diand AR studied, TIGEor TICEcan increase by 20 % at most.This indicates that the wing would benefit from thrust enhancement in the close vicinity of the ground or ceiling.The trend in the aerodynamic power is quite similar to that in the thrust, as shown in Fig. 8 (c) and 8(d).

    With a fixed Di,the tendency of thrust enhancement weakens as AR increases for the ground effect. Particularly, when Di≤0.2R, thrust experiences a sharp augmentation in the AR = 1 wing, while only minor differences occur among AR ≥3 wings. However, the thrust enhancement in AR = 1 and AR = 7 cases is noticeable for the ceiling effect, but approximately remains a relatively low level in AR = 3 and AR = 5 cases.

    According to the blade element theory(Eq.(14)),if ˙φ and α of the wing are kept unchanged,the thrust is highly dependent on the reduced velocity (vi). Taking the ceiling effect as an example,Fig.9(a)shows the distribution of vialong the wingspan of the AR=3 wing at different Dc.CFD-calculated and model-predicted viagree well with each other in the middle wingspan (approximately from r/R = 0.3 to r/R = 0.7), but show distinct discrepancy near the wing root and tip.In particular, the model predicts that vifor the OCE wing increases monotonously from the root to the tip,while CFD shows that it drops near the wingtip. Besides, in situations where Dc≤0.15, CFD gives viin opposite sign with a large change rate near the root area,which is not seen in model predictions.Such a discrepancy of vithat occurred near wingtip and root might be caused by the wing root vortex and tip vortex.However,both CFD and model results reveal that the magnitude of videcreases with the decrease of Dc.In the range of Dc/R ≤0.4,even a tiny fluctuation in Dccould lead to a noticeable change in vi. Besides, the spanwise distribution of vialso varies. With relatively large Dc, there is an obvious increase in vifrom r/R = 0.2 to r/R = 0.8. As Dcdecreases to a value less than 0.4R,the viprofile along the wingspan gets flattered within this region. The peak position of vigradually moves towards the wing root as Dcdecreases, resulting in a more severe change in vinear the distal end than the proximal end of the wing.Quantitatively, the peak viis -0.28Urefat r/R = 0.8 for the OCE wing, compared with - 0.12Urefat r/R = 0.4. Apart from vi, differences are also seen in the pressure coefficient Cpdistribution on the wing (Fig. 9 (b)). The enlarged area of negative pressure on the upper surface is the main reason for the thrust augmentation with decreasing Dc.Overall,the smaller the Dcis,the smaller the viis,and the more benefit the ICE thrust gains.A similar phenomenon is also true for the ground effect,reported by Prothin et al.15The closer the wing is placed to the ground,the more uniform the induced velocity distribution becomes, and the more significant the influence of it on the aerodynamic pressure distribution is.

    3.3. Comparison of aerodynamic performancebetween IGE and ICE wings

    Both ground and ceiling effects are found beneficial to thrust(T) augmentation for wings. Not only the thrust but also the power consumption (P) and aerodynamic efficiency (η) are of great importance in practical application. Table 5 shows the comparison of CFD results based on criteria including T, P,and η between ground and ceiling effects of different AR wings. Both T and P increase with the decrease of Di. And since T increases faster than P,49η gets enhanced as the wing approaches the ground or the ceiling. Table 5 also shows that the IGE wing gains more thrust augmentation than the ICE wing except for the AR=7 case.TIGE/TOGEis approximately 2%-8%larger than TICE/TOCEwhen AR <7.In general,the ceiling effect is more aerodynamically effective,with ηICE/ηOCEexceeding ηIGE/ηOGEby 2%-7%on average(see Fig.8(e)and 8(f)).

    Fig. 9 Distribution of induced velocity along wingspan, and static pressure on upper and lower surfaces (AR = 3) at different Dc.

    The AR = 5 wing is taken as an example. As shown in Fig. 10 (a) and 10(b), the sectional CTand CPincrease across the entire wingspan as Didecreases. Comparing the two wall effects, we can see that the IGE wing expresses a lower rate of growth in sectional CTand CPover the proximal end(0 < r/R < 0.5), but a higher rate over the distal end(0.5 <r/R <1) than the ICE wing. As for η, both ground and ceiling effects show enhancement in magnitude but withdifferent distribution regimes. The IGE wing has a higher η over the middle part of the wingspan, while the ICE wing has a higher η over root and tip ends, as shown in Fig. 10 (c).

    Table 5 Contrast of CFD-calculated T, P, and η of wings with different AR at different Di (averages over four revolving cycles).

    Fig. 10 Contrast of CFD-calculated distribution along wingspan of AR = 5 revolving wing (averages over four revolving cycles).

    Differences in T and P are reflected more intuitively by the distribution of the static pressure. Fig. 10 (d) depicts the average static pressure coefficient(Cp=2p/(ρUr2ef))over the entire upper and lower surfaces of the AR=5 wing.Quantitatively,the average Cpon the upper and lower surfaces of the OG(C)E wing are - 1.4694 and 0.3848, respectively; for the IGE wing,numbers are - 1.7135 and 0.7846; for the ICE wing, they are - 1.7483 and 0.4879. The ground effect creates a larger area of positive pressure with a higher magnitude on the lower surface. In comparison, the ceiling effect creates a larger area of negative pressure on the upper surface, but the peak values of the negative pressure in the two effects are at the same level.Overall, the integration of the pressure difference for the IGE wing(about 0.4928)is larger than that for the ICE wing(about 0.4687),resulting in greater aerodynamic forces.However,the corresponding power consumption inevitably increases because it is proportional to the drag in the revolving motion.Accordingly,MAVs with revolving wings can be applied in the close vicinity to walls for improved aerodynamic performance.Priority can be given to the ground effect for carrying a high payload and the ceiling effect to maintain long flight endurance.

    4. Conclusions

    (1) The theoretical models estimating the thrust and aerodynamic power of the revolving wing in two wall effects are improved.The model for the ground effect is established based on the wall-jet consumption, and that for the ceiling effect is improved by considering the uneven spanwise distribution of induced velocity.These two models are both validated by comparing the experimental measurements with CFD simulations in LBM methods. It is suggested that the model developed in this paper performs well when the wing-wall distance is within the range of 0 <Di≤2R. Although the applicability of our theoretical model is aspect-ratio-dependent, its predictions agree well with the existing experimental measurements. The error between theoretical estimations and CFD calculations is no more than 10 %, and keeps as low as 3 % on average in the range of AR = 3-5. The wall-jet assumption and the blade element momentum theory used to establish the model are validated to be reasonable by analyzing the flow field structures.

    (2) Based on the theoretical models and CFD simulation,the influence of wall clearance is addressed. Both ground and ceiling effects are found helpful for thrust enhancement,which increases by at least 20%compared to hovering cases without wall effects. The small wing-wall distance can highlight the aerodynamic benefit because of the reduced induced velocity and enlarged negative pressure area over the upper surface of the revolving wing.

    (3)The aerodynamic performance in the ground and ceiling effects is compared in terms of thrust generation and aerodynamic efficiency.It is found that operating MAVs with revolving wings in the close vicinity of the ground is more helpful to enhance thrust for carrying a high payload, while the ceiling effect is more helpful to improve aerodynamic efficiency for long flight endurance.

    Declaration of Competing Interest

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

    Acknowledgements

    This research was supported by the National Natural Science Foundation of China (No. 11902017) and the China Postdoctoral Science Foundation (Nos. 2020T130043 and 2019M650418). The authors would like to thank Dr. Long CHEN of Northeastern University for the helpful discussion regarding this project and the valuable comments on the manuscript.

    亚洲精品日韩av片在线观看| 成年美女黄网站色视频大全免费 | 亚洲在久久综合| 国产日韩欧美亚洲二区| 国产 精品1| 黄色视频在线播放观看不卡| 搡老乐熟女国产| 卡戴珊不雅视频在线播放| 国产高清有码在线观看视频| tube8黄色片| 高清毛片免费看| 国产久久久一区二区三区| 久热这里只有精品99| 精品酒店卫生间| 国产精品国产三级专区第一集| 丰满人妻一区二区三区视频av| 日韩强制内射视频| 天堂8中文在线网| 久久国产亚洲av麻豆专区| 五月玫瑰六月丁香| 黄色一级大片看看| 99久久人妻综合| 久久久久久人妻| 国产精品久久久久久精品古装| 国产精品久久久久久精品古装| 美女高潮的动态| 少妇 在线观看| 下体分泌物呈黄色| av女优亚洲男人天堂| 黄色视频在线播放观看不卡| 热re99久久精品国产66热6| 亚洲精品中文字幕在线视频 | 亚洲av.av天堂| 亚洲伊人久久精品综合| 亚洲国产av新网站| 少妇被粗大猛烈的视频| 亚洲成人中文字幕在线播放| 日本午夜av视频| 26uuu在线亚洲综合色| 少妇人妻久久综合中文| 一本—道久久a久久精品蜜桃钙片| 十分钟在线观看高清视频www | 国产熟女欧美一区二区| 国产成人精品福利久久| 少妇高潮的动态图| 欧美日韩综合久久久久久| 在线天堂最新版资源| 国产成人精品久久久久久| 日本-黄色视频高清免费观看| 天天躁夜夜躁狠狠久久av| 亚洲真实伦在线观看| 自拍偷自拍亚洲精品老妇| 国产午夜精品久久久久久一区二区三区| 欧美亚洲 丝袜 人妻 在线| 日韩人妻高清精品专区| 亚洲美女黄色视频免费看| 99热这里只有是精品50| 蜜桃久久精品国产亚洲av| 亚洲av在线观看美女高潮| 国产 一区 欧美 日韩| 少妇精品久久久久久久| 久久国产精品男人的天堂亚洲 | 久久av网站| 欧美+日韩+精品| 伦理电影大哥的女人| 国产毛片在线视频| 中文资源天堂在线| 久久毛片免费看一区二区三区| h视频一区二区三区| 亚洲人成网站在线观看播放| 亚洲国产精品999| 麻豆成人av视频| 国产在线免费精品| 国产在线免费精品| 人人妻人人澡人人爽人人夜夜| 免费av不卡在线播放| 国产熟女欧美一区二区| 天堂中文最新版在线下载| 国产高潮美女av| 国产精品一区二区性色av| 国产爱豆传媒在线观看| 亚洲av国产av综合av卡| av播播在线观看一区| av播播在线观看一区| 亚洲欧洲国产日韩| 97超视频在线观看视频| 亚洲国产欧美人成| 超碰av人人做人人爽久久| 大片免费播放器 马上看| 久久韩国三级中文字幕| 99久久精品国产国产毛片| 日韩 亚洲 欧美在线| 国产69精品久久久久777片| 亚洲国产成人一精品久久久| av在线app专区| 午夜免费男女啪啪视频观看| 久久久久久久久久成人| 久久国产精品大桥未久av | 你懂的网址亚洲精品在线观看| 纵有疾风起免费观看全集完整版| 男人舔奶头视频| 亚洲精品视频女| 精品国产三级普通话版| 国产伦精品一区二区三区视频9| 一本—道久久a久久精品蜜桃钙片| 亚洲国产欧美人成| 久久99热这里只频精品6学生| 在线观看免费高清a一片| 色婷婷久久久亚洲欧美| 国产精品嫩草影院av在线观看| 我的女老师完整版在线观看| 亚洲精品色激情综合| 中文精品一卡2卡3卡4更新| 一边亲一边摸免费视频| 校园人妻丝袜中文字幕| 亚洲欧美清纯卡通| 热99国产精品久久久久久7| 丝袜脚勾引网站| 国产亚洲一区二区精品| 18禁动态无遮挡网站| 精品亚洲成国产av| 中文资源天堂在线| 老司机影院毛片| 我的老师免费观看完整版| 成人黄色视频免费在线看| 中国三级夫妇交换| 久久ye,这里只有精品| 一二三四中文在线观看免费高清| 在现免费观看毛片| 一级毛片久久久久久久久女| 久久综合国产亚洲精品| 欧美成人a在线观看| 国产精品久久久久久av不卡| 亚洲精品久久午夜乱码| 深爱激情五月婷婷| 一级毛片久久久久久久久女| 丰满少妇做爰视频| 欧美xxxx性猛交bbbb| 精品99又大又爽又粗少妇毛片| 精品视频人人做人人爽| 黄色配什么色好看| 交换朋友夫妻互换小说| 国产一区二区在线观看日韩| 国产精品人妻久久久影院| 欧美少妇被猛烈插入视频| 伦理电影免费视频| 综合色丁香网| 欧美成人一区二区免费高清观看| 欧美xxxx性猛交bbbb| 精品99又大又爽又粗少妇毛片| 一级爰片在线观看| 国产成人午夜福利电影在线观看| 精品酒店卫生间| 汤姆久久久久久久影院中文字幕| 国产亚洲91精品色在线| 亚洲一级一片aⅴ在线观看| 亚洲一级一片aⅴ在线观看| 丰满迷人的少妇在线观看| av女优亚洲男人天堂| 亚洲精品乱码久久久v下载方式| 国产精品人妻久久久久久| 国产男女超爽视频在线观看| 国产黄色免费在线视频| 草草在线视频免费看| 一级a做视频免费观看| 99久久综合免费| 亚洲成人一二三区av| 日韩电影二区| 国产毛片在线视频| 国产av码专区亚洲av| 男女啪啪激烈高潮av片| 欧美日韩在线观看h| 久久久久视频综合| 久久ye,这里只有精品| 精品久久久精品久久久| 蜜桃亚洲精品一区二区三区| 国产淫语在线视频| h日本视频在线播放| 欧美亚洲 丝袜 人妻 在线| 久久久久久久久大av| 国产精品一区二区在线不卡| 久久午夜福利片| a级毛色黄片| 日韩电影二区| 欧美最新免费一区二区三区| 国产亚洲5aaaaa淫片| 一级毛片 在线播放| 久久综合国产亚洲精品| 成人漫画全彩无遮挡| 日本一二三区视频观看| 精品99又大又爽又粗少妇毛片| 国产精品不卡视频一区二区| 高清黄色对白视频在线免费看 | 日本av免费视频播放| 各种免费的搞黄视频| 欧美一区二区亚洲| 熟女人妻精品中文字幕| 老师上课跳d突然被开到最大视频| 国产欧美亚洲国产| 欧美亚洲 丝袜 人妻 在线| 激情五月婷婷亚洲| 亚洲国产最新在线播放| 国产亚洲一区二区精品| 男女无遮挡免费网站观看| 国产成人freesex在线| 国产女主播在线喷水免费视频网站| 免费播放大片免费观看视频在线观看| a级毛片免费高清观看在线播放| 午夜老司机福利剧场| 中文字幕制服av| 午夜免费男女啪啪视频观看| 寂寞人妻少妇视频99o| 少妇人妻久久综合中文| 91在线精品国自产拍蜜月| 97精品久久久久久久久久精品| 99热这里只有是精品50| 国产白丝娇喘喷水9色精品| 蜜臀久久99精品久久宅男| 精品午夜福利在线看| 王馨瑶露胸无遮挡在线观看| 一级av片app| 久久精品熟女亚洲av麻豆精品| 啦啦啦在线观看免费高清www| 成人18禁高潮啪啪吃奶动态图 | 午夜视频国产福利| 精品国产一区二区三区久久久樱花 | 久久精品久久久久久久性| 国产伦精品一区二区三区四那| 插逼视频在线观看| 草草在线视频免费看| 久久青草综合色| 少妇人妻精品综合一区二区| 大又大粗又爽又黄少妇毛片口| 精品国产露脸久久av麻豆| 视频区图区小说| 欧美日韩视频高清一区二区三区二| 色婷婷久久久亚洲欧美| 久久精品国产鲁丝片午夜精品| 成人亚洲欧美一区二区av| 日韩中字成人| 在线观看美女被高潮喷水网站| 欧美xxxx性猛交bbbb| 成人综合一区亚洲| 亚洲国产高清在线一区二区三| 色婷婷av一区二区三区视频| 久久精品久久精品一区二区三区| 天堂俺去俺来也www色官网| 久久这里有精品视频免费| 最近中文字幕2019免费版| 日韩精品有码人妻一区| 久久国产亚洲av麻豆专区| 久热久热在线精品观看| 日韩国内少妇激情av| 国产真实伦视频高清在线观看| 久久久国产一区二区| 国产精品嫩草影院av在线观看| 美女中出高潮动态图| 精品熟女少妇av免费看| 国产精品国产三级国产av玫瑰| 午夜老司机福利剧场| 日日啪夜夜爽| 欧美+日韩+精品| 亚洲欧美精品专区久久| 多毛熟女@视频| 国产成人aa在线观看| 人妻 亚洲 视频| 22中文网久久字幕| 亚洲国产欧美人成| 在现免费观看毛片| 免费少妇av软件| 成人毛片a级毛片在线播放| 欧美日韩在线观看h| 在线观看美女被高潮喷水网站| 三级国产精品片| 国产在线男女| 这个男人来自地球电影免费观看 | 国产精品女同一区二区软件| 秋霞在线观看毛片| 国产成人一区二区在线| 日本av免费视频播放| 日韩成人av中文字幕在线观看| 99热全是精品| 久久久精品免费免费高清| 欧美精品一区二区大全| 国产精品人妻久久久影院| 亚洲va在线va天堂va国产| 免费黄网站久久成人精品| 国产成人a区在线观看| 不卡视频在线观看欧美| 欧美bdsm另类| 丰满乱子伦码专区| 啦啦啦啦在线视频资源| 十八禁网站网址无遮挡 | 日本av免费视频播放| 久久97久久精品| 婷婷色av中文字幕| 国产精品99久久久久久久久| 国产精品福利在线免费观看| 国产精品精品国产色婷婷| 妹子高潮喷水视频| 尾随美女入室| 日本色播在线视频| 国产精品精品国产色婷婷| 国产精品久久久久久精品古装| 日韩成人伦理影院| 2022亚洲国产成人精品| 国产精品精品国产色婷婷| 色视频www国产| 一个人看视频在线观看www免费| 日韩中字成人| 一二三四中文在线观看免费高清| 三级国产精品片| 久久久久国产网址| 成人国产麻豆网| 亚洲精品国产av成人精品| 国产 一区 欧美 日韩| 国产免费福利视频在线观看| 亚洲内射少妇av| 韩国av在线不卡| 久久久久久久久久成人| 五月伊人婷婷丁香| 嫩草影院入口| av在线老鸭窝| 男的添女的下面高潮视频| 激情 狠狠 欧美| 中文字幕亚洲精品专区| 高清午夜精品一区二区三区| 日韩av免费高清视频| 女性被躁到高潮视频| 国产毛片在线视频| 日韩伦理黄色片| 好男人视频免费观看在线| 尤物成人国产欧美一区二区三区| 高清黄色对白视频在线免费看 | 亚洲精品国产成人久久av| 国产精品久久久久久久久免| 国产av码专区亚洲av| 舔av片在线| 亚洲人成网站在线播| 永久网站在线| 免费大片18禁| 直男gayav资源| 国产成人精品久久久久久| 成人漫画全彩无遮挡| 久久久久精品性色| 久久毛片免费看一区二区三区| 视频区图区小说| 成年人午夜在线观看视频| 成人高潮视频无遮挡免费网站| 免费不卡的大黄色大毛片视频在线观看| 久久久久国产精品人妻一区二区| 在线精品无人区一区二区三 | 亚洲人成网站在线播| 高清午夜精品一区二区三区| 一级a做视频免费观看| 美女高潮的动态| 国产精品.久久久| 亚洲av.av天堂| 国产乱人视频| 少妇裸体淫交视频免费看高清| 高清毛片免费看| 边亲边吃奶的免费视频| 免费大片18禁| 久久久久久久亚洲中文字幕| 汤姆久久久久久久影院中文字幕| 午夜福利视频精品| 久久婷婷青草| 国产一区二区在线观看日韩| 亚洲精品国产av成人精品| 国产精品人妻久久久久久| 色婷婷av一区二区三区视频| 简卡轻食公司| 尾随美女入室| 亚洲av成人精品一二三区| 99热这里只有是精品在线观看| 精品酒店卫生间| 午夜激情久久久久久久| 国产欧美日韩一区二区三区在线 | 97精品久久久久久久久久精品| 成人亚洲欧美一区二区av| 丰满乱子伦码专区| 高清在线视频一区二区三区| 日韩人妻高清精品专区| 毛片一级片免费看久久久久| 久久国产精品男人的天堂亚洲 | 在线观看免费日韩欧美大片 | 身体一侧抽搐| 观看美女的网站| 在线观看一区二区三区激情| 亚洲精品视频女| 精品久久久精品久久久| 国产午夜精品一二区理论片| 菩萨蛮人人尽说江南好唐韦庄| 狂野欧美白嫩少妇大欣赏| 成人免费观看视频高清| 国产精品免费大片| 晚上一个人看的免费电影| 大话2 男鬼变身卡| 啦啦啦视频在线资源免费观看| 日韩电影二区| 蜜臀久久99精品久久宅男| 毛片一级片免费看久久久久| a级一级毛片免费在线观看| 最新中文字幕久久久久| 亚洲av免费高清在线观看| 十八禁网站网址无遮挡 | 日韩av免费高清视频| 免费黄网站久久成人精品| 久久热精品热| 精品久久久久久久末码| 赤兔流量卡办理| 美女国产视频在线观看| 一级黄片播放器| 九九在线视频观看精品| 97在线视频观看| 九九爱精品视频在线观看| 我的老师免费观看完整版| av在线播放精品| 小蜜桃在线观看免费完整版高清| 26uuu在线亚洲综合色| 日韩 亚洲 欧美在线| 老司机影院毛片| 国产精品女同一区二区软件| 国产成人精品一,二区| 成人国产麻豆网| 韩国av在线不卡| 亚洲电影在线观看av| 嫩草影院新地址| 国产成人a区在线观看| 在线观看一区二区三区| 少妇人妻精品综合一区二区| 99视频精品全部免费 在线| 久久人人爽av亚洲精品天堂 | 高清日韩中文字幕在线| a级毛色黄片| av网站免费在线观看视频| 亚洲综合精品二区| 伊人久久精品亚洲午夜| 国产精品免费大片| 精品亚洲成国产av| 成人黄色视频免费在线看| 男人和女人高潮做爰伦理| 久久久精品免费免费高清| 国产精品一区二区性色av| 亚洲怡红院男人天堂| 91久久精品国产一区二区三区| 国产免费一区二区三区四区乱码| 97超碰精品成人国产| 国产无遮挡羞羞视频在线观看| 亚洲图色成人| 国产精品福利在线免费观看| 国产老妇伦熟女老妇高清| 插逼视频在线观看| 久久99精品国语久久久| 日韩中文字幕视频在线看片 | 国产午夜精品一二区理论片| 黑人猛操日本美女一级片| 老师上课跳d突然被开到最大视频| 欧美精品一区二区免费开放| 麻豆精品久久久久久蜜桃| 夜夜看夜夜爽夜夜摸| 在现免费观看毛片| 欧美精品人与动牲交sv欧美| 中文字幕精品免费在线观看视频 | 亚洲国产欧美在线一区| 丝瓜视频免费看黄片| 日本wwww免费看| 欧美3d第一页| 国产精品偷伦视频观看了| 久久久色成人| 中文天堂在线官网| 国产精品一区www在线观看| 欧美日韩综合久久久久久| 午夜免费鲁丝| 高清视频免费观看一区二区| 亚洲av电影在线观看一区二区三区| 成人一区二区视频在线观看| 一本—道久久a久久精品蜜桃钙片| 国产精品无大码| 嫩草影院入口| 女人久久www免费人成看片| 久久韩国三级中文字幕| 国产精品一区二区在线观看99| 王馨瑶露胸无遮挡在线观看| 国产有黄有色有爽视频| 久久国内精品自在自线图片| 亚洲欧美精品自产自拍| 一本久久精品| 老司机影院毛片| 欧美成人一区二区免费高清观看| av天堂中文字幕网| 夜夜骑夜夜射夜夜干| 精品亚洲成国产av| 亚洲电影在线观看av| 最近手机中文字幕大全| 亚洲中文av在线| 久久国产亚洲av麻豆专区| 女人久久www免费人成看片| 99久久综合免费| 插逼视频在线观看| 国产老妇伦熟女老妇高清| 亚洲一区二区三区欧美精品| 少妇熟女欧美另类| av网站免费在线观看视频| 日本一二三区视频观看| 99精国产麻豆久久婷婷| 99久久精品国产国产毛片| 国产爱豆传媒在线观看| 国产视频内射| 汤姆久久久久久久影院中文字幕| 日韩一区二区三区影片| 久久精品国产自在天天线| 国产欧美日韩一区二区三区在线 | 欧美高清成人免费视频www| 国产高清国产精品国产三级 | 女人十人毛片免费观看3o分钟| 久久鲁丝午夜福利片| 久久6这里有精品| 黄色一级大片看看| 一边亲一边摸免费视频| 国产精品麻豆人妻色哟哟久久| 国产精品三级大全| 国产精品不卡视频一区二区| 深爱激情五月婷婷| 国产欧美亚洲国产| 18禁在线无遮挡免费观看视频| 亚洲第一区二区三区不卡| 丰满少妇做爰视频| 久久av网站| 欧美一区二区亚洲| 九草在线视频观看| 美女内射精品一级片tv| 全区人妻精品视频| 人妻少妇偷人精品九色| 亚洲精品,欧美精品| 97精品久久久久久久久久精品| 亚洲天堂av无毛| 99久久人妻综合| 少妇人妻久久综合中文| 黄色欧美视频在线观看| 久久ye,这里只有精品| 亚洲最大成人中文| 成人无遮挡网站| 精品熟女少妇av免费看| 狂野欧美白嫩少妇大欣赏| 精品人妻视频免费看| 亚洲国产精品国产精品| 久久人人爽人人爽人人片va| 欧美日本视频| 国精品久久久久久国模美| 夫妻性生交免费视频一级片| 女性被躁到高潮视频| 男人舔奶头视频| 国产成人aa在线观看| 日本与韩国留学比较| 久久久久久久大尺度免费视频| 亚洲欧洲国产日韩| 老师上课跳d突然被开到最大视频| 国产欧美另类精品又又久久亚洲欧美| av专区在线播放| 超碰av人人做人人爽久久| 久久精品国产亚洲av涩爱| 一区二区三区四区激情视频| 九草在线视频观看| 99久久综合免费| 夜夜看夜夜爽夜夜摸| 人妻少妇偷人精品九色| 久久国产乱子免费精品| 精品亚洲成国产av| 国产视频首页在线观看| 成人免费观看视频高清| 日本一二三区视频观看| 免费观看在线日韩| 国产精品99久久99久久久不卡 | 2021少妇久久久久久久久久久| 亚洲av国产av综合av卡| 国产精品伦人一区二区| 搡女人真爽免费视频火全软件| 最近最新中文字幕大全电影3| 欧美少妇被猛烈插入视频| 亚洲国产日韩一区二区| 高清欧美精品videossex| 男女边摸边吃奶| 亚洲人成网站在线播| 大码成人一级视频| 尾随美女入室| 国产免费一区二区三区四区乱码| 亚洲第一区二区三区不卡| 黄色视频在线播放观看不卡| 永久免费av网站大全| 国产69精品久久久久777片| 日韩制服骚丝袜av| 精品视频人人做人人爽| 久久久久久久精品精品| 日本免费在线观看一区| 少妇人妻精品综合一区二区| 在线观看一区二区三区| 美女中出高潮动态图| 久久综合国产亚洲精品| 成人亚洲欧美一区二区av| 观看美女的网站| 免费av不卡在线播放| 色婷婷久久久亚洲欧美| 精品人妻视频免费看| 国产色爽女视频免费观看| 国产国拍精品亚洲av在线观看| 日韩国内少妇激情av| 国语对白做爰xxxⅹ性视频网站| 久久精品人妻少妇| 内地一区二区视频在线| 麻豆国产97在线/欧美| 欧美激情极品国产一区二区三区 | 国产毛片在线视频| 丰满少妇做爰视频| 免费高清在线观看视频在线观看|