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

    Pressure distribution guided supercritical wing optimization

    2018-09-27 07:08:06RunzeLIKaiwenDENGYufeiZHANGHaixinCHEN
    CHINESE JOURNAL OF AERONAUTICS 2018年9期

    Runze LI,Kaiwen DENG,Yufei ZHANG,Haixin CHEN

    School of Aerospace Engineering,Tsinghua University,Beijing 100084,China

    KEYWORDS Aerodynamic optimization;Man-in-loop;Pressure distribution;RBF-assisted;Supercritical wing

    Abstract Pressure distribution is important information for engineers during an aerodynamic design process.Pressure Distribution Oriented(PDO)optimization design has been proposed to introduce pressure distribution manipulation into traditional performance dominated optimization.In previous PDO approaches,constraints or manual manipulation have been used to obtain a desirable pressure distribution.In the present paper,a new Pressure Distribution Guided(PDG)method is developed to enable better pressure distribution manipulation while maintaining optimization efficiency.Based on the RBF-Assisted Differential Evolution(RADE)algorithm,a surrogate model is built for target pressure distribution features.By introducing individuals suggested by suboptimization on the surrogate model into the population,the direction of optimal searching can be guided.Pressure distribution expectation and aerodynamic performance improvement can be achieved at the same time.The improvements of the PDG method are illustrated by comparing its design results and efficiency on airfoil optimization test cases with those obtained using other methods.Then the PDG method is applied on a dual-aisle airplane’s inner-board wing design.A total drag reduction of 8 drag counts is achieved.

    1.Introduction

    Optimization design is more and more widely used to gain practical industrial designs in recent years.Multiple types of optimization algorithms have been well developed for more efficient and capable industrial optimal searching.However,although there have been many efforts using optimization algorithms to improve aerodynamic performances,practical design via optimization is still limited,and ‘‘cut-&-try” is still heavily relied on in the aircraft industry.1On one hand,even though multi-objective and multi-constraint optimization has become popular,2,3the complexity of industrial considerations makes it difficult to define all engineering-needed objectives and constraints for optimization algorithms.4,5On the other hand,engineers find their experiences,considerations,and judgments difficult to be introduced into an automatic optimization design process.1,4,6–8Therefore,in order to gain an engineering-acceptable design,optimization design not only needs a robust and flexible algorithm to endure a large amount of objectives and constraints,but also needs ways to transfer those in-brain requirements into what the optimization algorithm can accept or process.

    Experienced supercritical wing designers usually do not seek a wing with the highest lift/drag ratio,but a design that best compromises the performances of different disciplines and different flight conditions.What’s more,they emphasize much on design robustness,such as the drag divergence Mach number,the buffet onset lift coefficient,etc.Since both the performance and robustness are essentially the outcomes of flow structures as well as their evolutions,the flow physics is relatively clear.Engineers tend to judge a design through flow patterns and details.Currently,pressure distribution is the most cared flow structure in supercritical wing design.Many rules or criteria on pressure distribution have been proposed.For instance,the shock should be properly located to get a good robustness,and the aft loading of an airfoil should not be too large or else the nose-down pitching moment could be unacceptable.9–11

    By realizing these explicit rules of hint experiences on pressure distribution,a designer could improve a wing’s design point performances while achieving preferable off-design and multi-disciplinary properties.This is basically what they do in ‘‘cut-&-try”.From the 1980s,with the development of inverse design methods12–14researchers tried to gain a design by realizing a desired pressure distribution.However,due to the difficulties of ‘‘designing” a physically-existing pressure distribution,even for a simple geometry at a specific flow condition,11,15the application of these methods in the industry is still limited.

    Zhang et al.16studied the drag,moment,and especially robustness of three categories of typical pressure distributions for supercritical airfoils,i.e.,shock-free,double shock,and weak shock.The weak shock pressure distribution was evaluated as the best.Furthermore,other suggestions on pressure distribution have been proposed,such as that the shock location should range from 45%to 55%chord length for a single-aisle civil aircraft,the aft-loading should compromise the lift generation and a mild pitching moment,etc.To realize these suggestions,methods were developed to induce the pressure distribution to approach a desired pattern during optimization design.Some of them have been proven effective in industrial design.Optimal searching is no longer driven only by performances,but also by pressure distributions.Such methods are generalized in the present paper as Pressure Distribution Oriented(PDO)optimizations.

    The PDO method has also been applied to a dual-aisle airplane wing design.17The cruise performance and robustness were improved while the proposed weak shock pressure distribution was also achieved.The shock location was pushed slightly downstream that of a single-aisle civil aircraft to fit a higher Reynolds number and different cruise lift coefficients.Both the single-and dual-aisle design studies have shown that the location of shock wave is a critical factor for a supercritical wing’s balance of performance and robustness.

    Since the PDO method is characterized by the ability of manipulating pressure distributions,it can be consequently used to study the performance of a specified type of pressure distribution.By using the PDO method,Zhang et al.18achieved supercritical natural laminar airfoils with different pressure distributions characterized by assigned favorable pressure gradients and shock locations.Their gains on laminar friction reduction and penalties on the wave drag and robustness were then systematically compared.

    There are two types of PDO optimization developed in previous studies according to their methods of manipulating pressure distributions.The first one can be called Pressure Distribution Constrained(PDC)method.16–18Constraints are set to rule out or punish designs with an unsatisfying pressure distributions shape.It essentially posts restrictions on the optimal search direction.The optimization efficiency and global optimal searching capability are inevitably deteriorated.

    The other method is a manual or ‘‘man-in-loop”1,16one.Engineers can guide an optimization’s pressure distribution trend by manipulating the population.They need to introduce external individuals,which have the expected pressure distribution characteristics,into the population,or they need to eliminate unsatisfactory individuals from the population.This method demands a large amount of human labors and experiences.

    In this paper,a new approach of using the pressure distribution expectation to guide optimization is developed.Instead of manual population control,the so-called Pressure Distribution Guided(PDG)method uses a surrogate model to search potential individuals which best satisfy the expectation on pressure distribution and also have excellent performances.These individuals are introduced into the population of an evolutionary optimization algorithm.In this way,the optimization process is automatically guided to approach the expected pressure distribution feature,while the diversity of the population during the optimization is well preserved.To more rationally define the ‘‘expectation” on the pressure distribution,several physical or empirical relations are also studied.

    Unlike the PDC method,the pressure distribution is pursued by ‘‘guidance from good individuals” instead of ‘‘punishing bad individuals”.

    As a new branch of PDO optimization,in the present paper,the idea of the PDG method is firstly introduced.It is then tested and compared with previous methods by airfoil cases.Proven by results that it has better optimal searching efficiency and pressure distribution manipulation capabilities,the PDG method is applied to the supercritical wing design of a dual-aisle civil airplane.

    2.Optimization and modeling methods

    2.1.RBF assisted differential evolution algorithm

    An evolutionary optimization algorithm can use the assistance of surrogate models to improve efficiency.19In the present study,an RBF(radial basis function)Assisted Differential Evolution(RADE)algorithm is used as the primary optimization algorithm.20The basic flow chart of the RADE algorithm is outlined in Fig.1(a),in which k is the index of the current generation and P is short for population.The dash box contains the main optimization process of Differential Evolution(DE)optimization,and the upper part outside the box shows the RBF surrogate model’s behavior.By utilizing computed individuals’information,the RBF surrogate model can obtain approximation of CFD results.Then an optimal search is conducted on the RBF response surface to find individuals that are most likely to produce excellent objectives.Those individuals are added into the current candidate population to participate in the main DE optimization process.Since the search direction of an evolutionary algorithm is mainly based on the candidate population and the sorting strategy,this will form guidance to the search direction.In RADE,the surrogate model is not a substitution to an accurate CFD analysis.For individuals from both the main DE optimization process and the RBF surrogate model,a CFD objective evaluation is fully conducted.

    Fig.1 Flow charts of the RADE algorithm and PDC method and the PDG method.

    In an original RADE algorithm,the optimal search on the RBF surrogate model shares the same objectives and constraints as those of the main DE process.The purpose of the surrogate model is only to accelerate the main optimization process.The whole method can be used as a black-box multi-objective multi-constraints optimizer.The PDC and man-in-loop pressure distribution manipulation can both be directly realized on the RADE.The optimization can be directed to achieve the expected pressure distribution through manually adding/eliminating individuals by people,or through automatically eliminating individuals by constraints on the pressure distribution shape.

    There are drawbacks in these methods.For manual population manipulation,the basic idea is to manually judge all calculated individuals’pressure distributions and add or delete individuals to reconstruct the current population for evolutionary operators.On the contrary,the constraints of PDC optimization limit the optimization direction by eliminating non-feasible individuals,and then the population will meet the constraints eventually.Both methods alter the candidate population according to an engineer’s judgment,for example,their consideration on pressure distribution.However,eliminated individuals might have good objective function values,while artificially-introduced individuals usually do not have good enough performances to survive in the sorting procedure.Therefore,the optimization efficiency often suffers badly.If the introduced individuals could have competitive performances,the efficiency should be able to be maintained.

    Furthermore,proposed pressure distribution features may compromise the performance and robustness,and hence the method should only provide attempts on these features instead of forcing the entire population achieving them.

    The present PDG method is based on a modification to the RADE algorithm.The surrogate model contains the mapping between pressure distribution characteristic parameters and the design variable space.Objectives and constraints can be set different from those of the main DE process.A multiobjective optimal search on the RBF surrogate model is constructed to find potential individuals that have both good performances and expected pressure distribution features.These potential individuals are introduced to the DE’s candidate population to guide the optimization process,instead of forcing the entire population like the PDC method does.The constraints on the pressure distribution in the DE can therefore be relaxed or eliminated in the PDG process.The realization of the PDG method is shown in Fig.1(b).

    2.2.Modeling and numerical method

    The Class Shape Transformation(CST)method is widely used in wing design due to its capability of constructing a smooth airfoil with a few design variables.18In the present paper,a 6th-order Bernstein polynomial is used as the shape function for both the upper and lower surfaces of an airfoil.The number of design variables for an airfoil is 14.The ranges of the variables are the same as those in Ref.18.The wing surface is constructed by interpolation with several span-wise distributed airfoils.

    The aerodynamic performances of airfoils and wings are evaluated by an in-house solver,21,22which is a cell-centered finite volume compressible Reynolds Average Navier-Stokes(RANS)solver with multiple types of numerical schemes and turbulence models.The solver has been validated by various test cases and successfully applied in many industrial design activities.16–18,21–23In this paper,it employs the MUSCL scheme for the reconstruction,Roe’s scheme for the spatial discretization,the lower-upper symmetric Gauss-Seidel method for the time advancing,and the Shear Stress Transport(SST)model for the turbulence modeling.The accuracy of the solver has been validated by our previous study.23

    3.Method validation

    3.1.Pressure distribution consideration of a supercritical wing

    For a supercritical wing,the pressure distribution contains various information of aerodynamic characteristics.Since the off design performances are essentially the outcome of the flow structure evolution among different flight conditions,it is not necessary to evaluate the drag divergence Mach number or buffet onset point through CFD calculations of a series of flight conditions.Pressure distribution manipulation and performance optimization on several typical conditions can suffice.1,17

    On the other hand,many important characteristics are directly associated with certain flow patterns,such as the wave drag and the shockwave strength.Some experience of pressure distribution has been gathered to gain a fine compromise among lift,drag,moment,shockwave stability,trailing edge separation characteristic,and geometry considerations.16However,features are different for different flight conditions and different aircraft planforms.A dual-aisle aircraft usually has a higher flight Mach number and a higher Reynolds number than those of a single-aisle aircraft,and the shockwave is downstream with a higher strength.Therefore,the wave drag is a major party in drag reduction,and the off-design characteristic is more non-linear.In this paper,several pressure distribution considerations16of a wing of a dual-aisle aircraft at a flight Mach number of 0.85 are taken in the following optimizations.The pressure distribution considerations of an airfoil are transformed according to the cosine rule of a swept wing.23These considerations can be expressed as follows:

    (1)The pressure coefficient of the suction peak of a wing cannot be lower than-1.0;otherwise,the flow accelerates too much around the leading edge,which will make the aerodynamic performance like that of ‘‘peaky”airfoils.24

    (2)The slope of the pressure plateau on the upper surface(from the suction peak to the shockwave front)shall be between-0.2 and 0.5 to ensure an appropriate length of the suction platform for providing enough lift.

    (3)The pressure plateau shall be as smooth as possible to avoid undesirable robustness issues.

    (4)The slope of the pressure recovery zone near the trailing edge of a wing is kept lower than 3.0 to prevent flow separation.

    (5)The aft loading shall not be too large,or else the nosedown pitching moment will be too large.

    3.2.Relation between shockwave strength and drag coefficient

    The wave drag is directly related to the shockwave strength.Although the shockwave strength may have different definitions in previous studies,it can be roughly described by the pressure rise between the two sides of the shockwave ΔCp,as shown in Fig.2.According to the Oswatitsch theorem,25the wave drag of an airfoil can be described as in Eq.(1)under the assumption of isentropic flow before the shockwave,where Ma1is the Mach number at the wave front.

    The pressure coefficient can be defined aswhere γ is the specific heat ratio,Ma∞is the freestream Mach number,and p/p∞is the ratio of static pressure normalized by the freestream value.Then ΔCpcan be described as in Eq.(2)using the isentropic relation and the Rankine-Hugoniot equation.

    Fig.2 Definitions of Err and ΔCp.

    Then the wave drag can be expressed as in Eq.(3),where S is a coefficient associated with the airfoil shape and the flight condition.

    A large number of supercritical airfoils through optimization are gathered to verify the relation.The airfoils are evaluated by RANS.The free stream Mach number is 0.742,and the Reynolds number is 10 million.The lift coefficients of the airfoils are kept the same(CL=0.787).The airfoils all have similar flow patterns,which are a single shockwave and a smooth suction plateau,and all airfoils have normalized chord length,i.e.X ∈ [0,1].The smoothness function Err is defined as Cpfluctuations on the suction plateau,which is the area of the blue shadow region in Fig.2.The pressure distribution in Fig.2 has a smoothness function of Err=0.012.The drag coefficients of the sampled airfoils which have small fluctuations(Err<0.012)are shown in Fig.3.A curve fitting relation between the drag coefficient CDand the shockwave strength ΔCpis shown as the red curve in Fig.3,where the coefficient S equals to 0.05,and CD0is the smallest drag coefficient for non-shockwave airfoils.

    Fig.3 Relation between total drag and shockwave strength.

    The sampled airfoils fit the curve with a less than 0.1 derivation of ΔCp.Therefore,the shockwave strength is correlated with the total drag and can be used as an effective auxiliary parameter for drag reduction.It is also a parameter to distinguish between different pressure distribution types.As shown in Fig.3,when ΔCpis less than 0.3,the drag coefficient penalty caused by the shock wave is smaller than 0.0005,and the drag reduction caused by shockwave weakening is negligible when ΔCpis less than 0.2.Moreover,it demonstrates that airfoils with the same total drag coefficient can have different shockwave strengths,which reveals that while keeping the total drag unchanged,the shock strength can be adjusted to get a better robustness.Therefore,if the shockwave strength is used as an objective or sub-objective in a drag reduction optimization,it could have the potential to increase the overall optimization efficiency or robustness.

    3.3.Single-objective PDG optimization

    Pressure distribution features can be used as objectives in the main optimization process,or in the sub-optimization process.In this section,single-objective optimizations of a supercritical airfoil with different optimization settings are compared to demonstrate the acceleration effect of sub-optimization.The free stream Mach number is 0.742,the Reynolds number is 10 million,and the lift coefficient is kept 0.787 during the optimization.The design variables are 14 CST parameters.

    Five settings are studied in this case,as shown in Table 1,where RLEis the airfoil leading edge radius.Opt1 and Opt2 are baseline PDC optimizations using DE and RADE algorithms,respectively.Opt3 and Opt4 directly use the shockwave strength as an overall objective of RADE.Opt1 has no suboptimization process.Opt2,Opt3,and Opt4 have the same objectives and constraints in both the main optimization process and the sub-optimization;therefore,Opt3 and Opt4 are still called PDC optimizations.Opt5 uses the modified RADE with different objectives and constraints in the main optimization process and the sub-optimization,which is a demonstration of the present PDG method.According to the relation between the shockwave strength and the drag coefficient in Section 3.2,ΔCpis used as a sub-objective to guide the optimization generating airfoil designs for an inboard wing,of which ΔCpis expected to be smaller than 0.2.Pressure distribution considerations described in Section 3.1 are also applied as constraints in the main optimization.

    The population sizes of all the cases are 32,and 32 generations are carried out to compare the efficiency and effectiveness.All the cases have the same initial population,and the averaged drag coefficients of the final generation are used for comparison,as shown in Table 2.The convergence histories are shown in Fig.4.Table 2 shows that Opt2 has a slightly better final result than that of Opt1,and a result that has a similar performance to the final result of Opt1 can be achieved by Opt2 with only around 600 individuals,as shown in Fig.4,where ID is the index of each individual.Therefore,the RBF-assisted method could improve the efficiency of the optimization.

    Table 1 Optimization settings of 5 cases.

    Table 2 Final results of 5 cases.

    Fig.4 Convergence histories of the five cases.

    Fig.5 shows two typical pressure distributions of the results which have an unsatisfactory performance.The strong shockwave airfoil(Fig.5(a))whose CDequals to 0.01071 is the initial design.Although reducing the shockwave strength can reduce the wave drag,shock-free and weak-shock airfoils are not the only possibilities for reducing the shockwave.Due to a lack of proper constraints,Opt3 wastes part of its ability to search in a wrong direction,i.e.,airfoils with massive upper surface separation or double-shockwave airfoils(Fig.5(b)).Moreover,excessive constraints can have negative effects on the main process,as shown in Opt4,and the overall efficiency is low.Opt5 uses the shockwave strength as a sub-objective to provide an alternative searching direction for the main process.Because the computation cost for searching on the RBF is negligible,the efficiency of the main optimization process is not depressed.Table 2 shows that Opt5 has the best final result in all of the five optimizations.Fig.4 shows that a result similar to the final result of Opt2 can be achieved by Opt5 with around 500 individuals.Therefore,the optimization is accelerated by the sub-optimization,and the desired pressure distribution feature is achieved.

    Fig.6 shows the final result of Opt5,of which CDequals to 0.00958.The pressure distribution type is shock-free.The pressure distribution satisfies all of the constraints in the optimization,but the robustness is not good according to Ref.16.This case is only a demonstration of the PDG optimization.Multi-point optimization is carried out in the next section to improve the robustness.

    Fig.5 Typical pressure distributions with an unsatisfied performance.

    In summary,the PDG optimization is achieved by a combination of the main optimization process and sub-optimization.As shown by the results of Opt3 and Opt4,the pressure distribution features used as main process objectives might lead to undesired results,while adding more constraints in the main optimization process leads to a low efficiency.The PDG optimization,i.e.,Opt5,uses sub-objectives and sub-constraints to utilize the pressure distribution expectation,and the efficiency is not compromised.Moreover,the sub-optimization only serves as a source of potential individuals to the main process.Therefore,the sub-optimization can handleflexible expectations of engineers without interfering the entire optimization efficiency.

    3.4.Multi-point pressure distribution guided optimization

    Different methods of pressure distribution oriented optimization along with a baseline optimization are compared in this section.Cases are used to demonstrate the PDG method’s effectiveness of pressure distribution manipulation and optimization.A multi-point airfoil optimization is applied based on the RADE algorithm.The airfoil is expected to have a shockwave location at 50%chord length at the cruise condition.This case is to demonstrate the pressure distribution control ability of the present method.In order to make the comparison simple and clear,only the cruise drag(CD0)is used as the objective of the main optimization process.Drag divergence and buffet onset control are maintained through constraints.Flight conditions are shown in Table 3.The objectives and constraints of 4 different settings can be seen in Table 4.

    Fig.6 Pressure distribution of Opt5’s final result.

    Table 3 Flight conditions of three design points.

    Opt1 is the baseline optimization.Opt2 is a demonstration of manual population manipulation.Since the final effect can be highly influenced by different engineers’involvement,a fixed strategy is programed into the optimization sorting strategy to mimic the manual manipulation.The basic idea of manual manipulation is to select individuals for the DE operators based on personal judgement of performance and pressure distribution considerations.In this case,individuals with a low drag and a single shockwave in all 3 flight conditions as well as a small value of| XSW0-0.5|are more likely to be selected.The strategy for Opt2 is shown in Fig.7.Opt3 applies constraints of the shockwave position| XSW0-0.5|<0.05.Opt1,Opt2,and Opt3 share the same objectives and constraints in the main optimization process and the sub-optimization.Opt4 uses the PDG method,which means more objectives and constraints are applied in the sub-optimization to guide the search direction of the optimization.According to the relation between the shockwave strength and the drag coefficient in Section 3.2,the shockwave strength at the cruise condition is constrained by ΔCp0< 0.3 in the sub-optimization.

    Fig.8 shows the convergence histories of the optimizations.Fig.8(a)is the objective of the cruise condition.Fig.8(b)is the pressure distribution consideration,i.e.,the distance between the shockwave and the half chord length.Fig.8(c)and(d)show the drag increments for the drag divergence condition and the buffet onset point.The dash lines in Fig.8(b)–(d)show the critical constraint values.Table 5 shows some statistical results of the methods.Fig.9(a)shows the individuals within the requested shock wave region,i.e.,45–55%,and only the ones meeting all the constraints are included.A typical optimized result is shown in Fig.9(b).It is a weakshockwave airfoil with a cruise drag of 0.00988.

    As shown by the results in Fig.8(b),compared to Opt1,all the pressure distribution oriented methods are able to get shockwaves in the requested region.The results of Opt1 are almost equally distributed in different shock locations.By contrast,the results of Opt2 to Opt4 are getting near 50%chord length when the optimization proceeds.Table 5 shows that Opt2 and Opt3 have fewer results satisfying Opt1 constraints than those of Opt1,whereas Opt2 and Opt3 have more results satisfying Opt3 constraints.Therefore,although manual population manipulation and constraints in Opt3 and Opt4 have a better efficiency generating results with shockwaves in the requested region,the ability of generating results satisfying basic constraints,i.e.,Opt1 constraints,is compromised.By contrast,Opt4 has more results satisfying Opt1 and Opt3 constraints,and these results have shockwaves near the required 50%chord length location.

    According to the results of Opt2 and Opt3,we can see that the constraint in the main optimization process is more objective than manual manipulation,but also more in flexible.When the critical value for| XSW0-0.5|constraint is small,the constraint may eliminate most individuals in the early stage of the optimization,and a convergence is hard to achieve.However,if a loose constraint is applied,the final results cannot satisfy the expectation of the shockwave location.

    Opt4 applies more objectives and constraints in the sub-optimization to guide the searching direction.Thesub-optimization provides promising individuals for the main process.Since the sub-constraints are stronger and the shockwave strength along with| XSW0-0.5|are also used as subobjectives,the individuals provided by the sub-optimization have a higher possibility to reduce the drag and meet the constraints synchronously.Moreover,because the drag divergence and buffet characteristics are associated with shockwave induced boundary layer separation,the sub-constraints about the shockwave strength in these flight conditions can help improving the robustness,and then accelerating the main optimization process.Therefore,the PDG method could increase the optimization efficiency and generate more desirable results.

    Table 4 Four different optimization settings.

    Fig.7 Flow chart with an additional sorting strategy for Opt2.

    Table 5 Statistics of the 4 optimization settings.

    Fig.9 Final results of the 4 optimization settings.

    4.Method application

    Previous studies16,17have indicated that a weak-shockwave pressure distribution has a better balance between the cruise performance and the robustness than those of doubleshockwave and shock-free pressure distributions on a singleaisle airplane.However,since a dual-aisle airplane has a different planform and a different flight condition,the details of preferable features are unclear,especially for the inboard wing.The inboard wing has a much longer chord length than that of the outboard wing,and hence a small modification of the inboard wing can have a significant impact on the span-wise load distribution and a strong influence on the wing performance.Meanwhile,the inboard wing’s pressure distribution can directly influence the outboard wing due to the cross flow,and the interference of the wing,fuselage,and nacelle/pylon makes the non-linear phenomenon even more severe than on a single-aisle airplane.As from the discussion of the CD- ΔCprelation in Section 3.2,once ΔCpfor an airfoil is lower than 0.3,the drag increment caused by a shockwave is insignificant.Therefore,a weak shockwave of ΔCp=0.3 can be a good balance between cruise drag and robustness.

    For a typical supercritical airfoil,once the shockwave strength ΔCpis fixed,the section lift coefficient is mainly dependent on the suction plateau and the shockwave location.A higher suction peak causes a strong pressure recovery on the suction plateau,and may induce unsatisfying robustness.A further downstream shockwave location tends to have a slightly stronger pressure recovery behind the shockwave,which may increase the risk of shockwave-induced boundary layer separation,but it also provides more lift.It is especially important for the inboard wing,because an increment of the lift coefficient of inboard wing sections can lead to a large reduction of the outboard wing load,which is good for improving the buffet characteristics.Therefore,when the shockwave strength ΔCpis kept lower than 0.3,a further downstream shockwave location and a lower suction peak on the inboard wing may result in a more excellent lift/drag ratio and robustness.

    A dual-aisle airplane is used as the test case for the inboard wing optimization based on the PDG method.The cruise Mach number is 0.85,the Reynolds number is 40 million,and the lift coefficient is 0.48.The baseline wing was designed in our previous study,23which mainly focused on the outboard wing optimization.The baseline wing has an excellent robustness and a satisfying weak-shockwave pressure distribution on the outboard part.However,there is a strong shock on the inboard wing.In this paper,it is used to study different inboard wing pressure distributions and the influence on the overall performances.The inboard wing geometry is interpolated by 3 airfoils generated by the CST method.The locations are shown as Sections 1–3 in Fig.10.The design variables are the 42 CST variables of the 3 airfoils,and the computational mesh is the same as that in our previous study.23The black solid and dash lines in Fig.10 show the shockwave location of the baseline wing and the expected inboard wing shockwave location,respectively.

    Fig.10 Wing section locations and shockwave locations.

    The PDG method is used to manipulate the pressure distribution while reducing the cruise drag.Optimization settings are listed in Table 6.The cruise drag is set as the main objective and sub-objective.The distance between the actual shockwave location and the proposed position as well as the shockwave strength are also adopted as sub-objectives.The pressure distribution considerations mentioned in Section 3.1 are used as constraints,as well as the constraints about the geometry and the shockwave strength.The optimization population in each generation has a size of 40 individuals,and total 25 generations are carried out to gain final results.The convergence history is shown in Fig.11(a).Results with and without proposed features are used to compare the performances.Four typical results in Fig.11(b)are selected to compare the performances and pressure distributions.The pressure distributions are shown in Fig.12.The characteristics of the pressure distributions are shown in Table 7.The average values in the table are based on the values of Sections 1–4 in Fig.10.

    Design4 in Table 7 is the final result of the optimization.It has the minimum drag coefficient of all feasible results and a short distance to the proposed shock location.Although it has a weak double shock on the 17.5%half-span section,it does not violate the constraints of the main optimization pro-cess.Design3 is also a typical result for drag reduction.However,it cannot satisfy the constraint of the shockwave location,and its suction peak is high.Design2 shows that when the double shock in Design4 is eliminated,the drag reduction quantity is limited.If a lower suction peak is applied,i.e.,in Design1,the cruise drag is compromised.

    Table 6 Settings of dual-aisle airplane inboard wing optimization.

    Fig.11 Optimization history and results.

    Fig.12 Wing section pressure distributions of typical results.

    Table 7 Cruise performances of the baseline and typical results.

    Fig.13 shows the shockwave locations and span-wise loads of the wings.Basically,a further downstream shockwave has a larger inboard wing load(Design4);consequently,the maximum lift coefficient on the outboard wing is reduced,which usually means a weak shockwave strength(see Fig.12,17.5%,32.7%,and 45.0%sections)and a potentially better buffet characteristic.

    Fig.14(a)shows the Cmvs CLcurves of the 4 results.Fig.14(b)shows the CDvs Ma curve.Inflection in a Cmvs CLcurve can be used to define the buffet onset26.The CDvs Ma curve is used to determine the drag divergence Mach number.It can be seen that designs with a larger inboard wing load have a larger buffet onset lift coefficient.Although Design4 has a similar cruise drag and a similar drag divergence performance to those of Design3,it has a better buffet characteristic.Through the present optimization with the PDG method,the proposed shockwave location and the pressure distribution expectations of a dual-aisle airplane are achieved.The result shows that the optimized design has a low drag and good buffet and drag divergence performances.

    Fig.14 Buffet and drag divergence characteristics.

    5.Conclusions

    The pressure distribution of a supercritical wing contains comprehensive information about its performance and robustness.Pressure distribution oriented optimization can help achieve more engineering-applicable designs.This paper has developed a pressure distribution guided method of PDO optimization to introduce the considerations of pressure distribution into the optimization process.The work of this paper can be summarized as follows:

    (1)The PDG method is developed based on an RBF-assisted differential evolution algorithm.The suboptimization on the RBF response surface can have objectives and constraints different from those of the main optimization process.By generating candidate individuals more likely to have both satisfactory pressure distribution features and performances,the suboptimization enables the PDG method to better manipulate the pressure distribution while maintaining the aerodynamic optimization efficiency.

    (2)The relation between the shockwave strength and the drag coefficient is theoretically and statistically discussed in this paper.The result shows that the shockwave influence on the drag coefficient is insignificant when the shockwave strength is lower than 0.3.

    (3)Single-and multi-point optimizations of supercritical airfoils are carried out to study the efficiency of the PDG method.The results show that setting the shockwave strength as the objective for the sub-optimization can accelerate the drag optimization process,and the PDG method is able to more effectively achieve the pressure distribution expectation.

    (4)The PDG method is applied on the inboard wing design of a dual-aisle airplane.A proposed shockwave location of the inboard wing is achieved through different pressure distribution manipulation methods.The performance and robustness of the final result and some other typical results are compared.The results show that the PDG method can gain designs with the proposed pressure distribution while maintaining the aerodynamic optimization efficiency.The results also show that a weak-shockwave pressure distribution is a satisfactory balance between cruise drag and robustness.

    Acknowledgements

    This work was co-supported by the National Key Basic Research Program of China(No.2014CB744806)and Tsinghua University Initiative Scientific Research Program(No.2015Z22003).

    精品国产乱码久久久久久男人| 日本av免费视频播放| 国产麻豆69| 午夜激情久久久久久久| 在线观看免费高清a一片| 777久久人妻少妇嫩草av网站| 乱人伦中国视频| 精品乱码久久久久久99久播| 高清毛片免费观看视频网站 | 欧美+亚洲+日韩+国产| 国产精品久久久av美女十八| 性色av乱码一区二区三区2| 免费观看a级毛片全部| 久久久精品94久久精品| 手机成人av网站| 亚洲伊人色综图| 天堂动漫精品| 久久精品成人免费网站| 正在播放国产对白刺激| 成年人午夜在线观看视频| 精品亚洲乱码少妇综合久久| 国产一区二区三区在线臀色熟女 | 男女下面插进去视频免费观看| 亚洲国产av新网站| 免费在线观看日本一区| 欧美日韩亚洲高清精品| 中文字幕人妻丝袜一区二区| av福利片在线| 久久 成人 亚洲| 丰满饥渴人妻一区二区三| 久久精品熟女亚洲av麻豆精品| 俄罗斯特黄特色一大片| 欧美黄色淫秽网站| 一本综合久久免费| 国产一区有黄有色的免费视频| 亚洲av日韩精品久久久久久密| 午夜精品久久久久久毛片777| 久久久久久久久免费视频了| 天堂俺去俺来也www色官网| 人妻 亚洲 视频| 国内毛片毛片毛片毛片毛片| 精品国产乱码久久久久久小说| 久久久精品94久久精品| 国产精品99久久99久久久不卡| av天堂在线播放| 大香蕉久久网| 日韩有码中文字幕| 少妇猛男粗大的猛烈进出视频| 亚洲精品在线美女| 91麻豆精品激情在线观看国产 | 老司机影院毛片| 欧美日韩一级在线毛片| 亚洲欧美日韩高清在线视频 | 日本wwww免费看| 中文字幕av电影在线播放| 建设人人有责人人尽责人人享有的| 男女无遮挡免费网站观看| 精品免费久久久久久久清纯 | 一本大道久久a久久精品| 自线自在国产av| 19禁男女啪啪无遮挡网站| 久久午夜亚洲精品久久| 久久人妻av系列| 亚洲精品成人av观看孕妇| av视频免费观看在线观看| 精品免费久久久久久久清纯 | 老司机深夜福利视频在线观看| 午夜免费成人在线视频| 亚洲av电影在线进入| 一区二区av电影网| 午夜福利视频精品| 女人被躁到高潮嗷嗷叫费观| 老汉色av国产亚洲站长工具| 成人黄色视频免费在线看| 丰满迷人的少妇在线观看| 亚洲欧美日韩另类电影网站| 亚洲成人国产一区在线观看| 免费人妻精品一区二区三区视频| 久久久精品94久久精品| 欧美成人午夜精品| 精品乱码久久久久久99久播| 青草久久国产| 久久青草综合色| 三上悠亚av全集在线观看| 国产在线视频一区二区| 国产熟女午夜一区二区三区| 狠狠精品人妻久久久久久综合| 十八禁网站网址无遮挡| 亚洲熟女精品中文字幕| 亚洲欧美日韩高清在线视频 | 在线av久久热| 高清欧美精品videossex| 欧美精品高潮呻吟av久久| 亚洲九九香蕉| 亚洲伊人久久精品综合| 女人爽到高潮嗷嗷叫在线视频| 新久久久久国产一级毛片| 国产精品av久久久久免费| 亚洲精品乱久久久久久| 色婷婷久久久亚洲欧美| 一区福利在线观看| 电影成人av| 正在播放国产对白刺激| 天堂俺去俺来也www色官网| 十八禁网站网址无遮挡| 欧美人与性动交α欧美精品济南到| 免费一级毛片在线播放高清视频 | 高清av免费在线| 国产精品98久久久久久宅男小说| 久久精品aⅴ一区二区三区四区| 久久天躁狠狠躁夜夜2o2o| 国产人伦9x9x在线观看| 另类精品久久| 91精品三级在线观看| 精品国产乱码久久久久久小说| 男男h啪啪无遮挡| 国产日韩欧美亚洲二区| 精品国产乱码久久久久久男人| 人妻久久中文字幕网| 欧美日韩精品网址| 久久中文看片网| 欧美人与性动交α欧美精品济南到| 成人精品一区二区免费| 新久久久久国产一级毛片| 亚洲av日韩在线播放| 国产精品 欧美亚洲| 亚洲欧美日韩另类电影网站| 国产在线一区二区三区精| 欧美日韩中文字幕国产精品一区二区三区 | 黑丝袜美女国产一区| 亚洲第一欧美日韩一区二区三区 | 99国产综合亚洲精品| 国产单亲对白刺激| 欧美日韩av久久| 亚洲男人天堂网一区| 免费在线观看黄色视频的| 色播在线永久视频| 天天躁狠狠躁夜夜躁狠狠躁| 电影成人av| 精品一区二区三区av网在线观看 | 久久久久国内视频| 久久国产精品人妻蜜桃| 三上悠亚av全集在线观看| 香蕉丝袜av| 熟女少妇亚洲综合色aaa.| 亚洲欧美激情在线| 黄色视频不卡| av片东京热男人的天堂| 国产老妇伦熟女老妇高清| av视频免费观看在线观看| 免费观看a级毛片全部| 在线观看舔阴道视频| 国产真人三级小视频在线观看| 中亚洲国语对白在线视频| 精品一区二区三区av网在线观看 | 国产欧美日韩一区二区精品| 亚洲av片天天在线观看| 国产又色又爽无遮挡免费看| 久久性视频一级片| 91精品国产国语对白视频| 国产成人系列免费观看| av超薄肉色丝袜交足视频| 久久人人爽av亚洲精品天堂| 黑人巨大精品欧美一区二区蜜桃| 久久精品国产亚洲av高清一级| 国产男女内射视频| 搡老熟女国产l中国老女人| 一二三四在线观看免费中文在| 日本av免费视频播放| 久久国产精品人妻蜜桃| 国产免费视频播放在线视频| 超碰97精品在线观看| h视频一区二区三区| 十八禁网站免费在线| 性色av乱码一区二区三区2| 如日韩欧美国产精品一区二区三区| 免费av中文字幕在线| 不卡av一区二区三区| tube8黄色片| 丝袜美腿诱惑在线| 免费观看a级毛片全部| 777米奇影视久久| 一本综合久久免费| 99久久人妻综合| 美女高潮到喷水免费观看| 日本av免费视频播放| 久久午夜综合久久蜜桃| 久久国产精品影院| 午夜福利乱码中文字幕| 肉色欧美久久久久久久蜜桃| 午夜精品久久久久久毛片777| 亚洲情色 制服丝袜| 一级a爱视频在线免费观看| 国产精品久久久av美女十八| 18禁美女被吸乳视频| 一级毛片电影观看| 丝袜在线中文字幕| 狠狠狠狠99中文字幕| 一级毛片电影观看| 老熟妇仑乱视频hdxx| 男女之事视频高清在线观看| 国产高清国产精品国产三级| av天堂在线播放| 亚洲精品自拍成人| 午夜福利视频精品| 精品人妻在线不人妻| 99久久国产精品久久久| 757午夜福利合集在线观看| 国产欧美日韩精品亚洲av| 777米奇影视久久| 亚洲国产欧美日韩在线播放| 国产成人精品久久二区二区91| 国产成人av激情在线播放| 在线观看人妻少妇| 两性夫妻黄色片| 亚洲自偷自拍图片 自拍| 99re6热这里在线精品视频| 国产日韩一区二区三区精品不卡| 天堂俺去俺来也www色官网| 日韩熟女老妇一区二区性免费视频| 国产区一区二久久| 亚洲五月色婷婷综合| 国产成人精品在线电影| 成人18禁在线播放| 99久久人妻综合| 视频在线观看一区二区三区| 亚洲国产欧美网| 久久精品国产亚洲av香蕉五月 | 精品午夜福利视频在线观看一区 | 夜夜爽天天搞| 在线播放国产精品三级| 亚洲精品在线观看二区| 亚洲精品自拍成人| 在线观看舔阴道视频| 欧美激情 高清一区二区三区| 别揉我奶头~嗯~啊~动态视频| 久久久久视频综合| 纵有疾风起免费观看全集完整版| 一本色道久久久久久精品综合| 啦啦啦视频在线资源免费观看| 亚洲成人免费av在线播放| 国产精品av久久久久免费| 久久久久久亚洲精品国产蜜桃av| www日本在线高清视频| 丝袜喷水一区| 亚洲av日韩在线播放| 久久ye,这里只有精品| 12—13女人毛片做爰片一| 成人免费观看视频高清| 亚洲精品乱久久久久久| 亚洲中文日韩欧美视频| 午夜福利视频在线观看免费| 久久99一区二区三区| 国产成人免费观看mmmm| 热re99久久国产66热| 91成人精品电影| 久久九九热精品免费| 香蕉国产在线看| av不卡在线播放| 国产在线视频一区二区| 久久久精品94久久精品| 国产亚洲av高清不卡| 日韩视频在线欧美| 老司机影院毛片| 日韩一卡2卡3卡4卡2021年| 国产99久久九九免费精品| 99久久人妻综合| 最新的欧美精品一区二区| 亚洲午夜理论影院| 亚洲专区中文字幕在线| 女人久久www免费人成看片| 午夜精品久久久久久毛片777| 窝窝影院91人妻| 日韩制服丝袜自拍偷拍| 18禁美女被吸乳视频| 亚洲综合色网址| 激情在线观看视频在线高清 | 女同久久另类99精品国产91| 国产麻豆69| 亚洲精华国产精华精| 丰满迷人的少妇在线观看| 天天影视国产精品| 在线观看免费午夜福利视频| 激情在线观看视频在线高清 | 久9热在线精品视频| 下体分泌物呈黄色| 极品少妇高潮喷水抽搐| 久久精品91无色码中文字幕| 国产精品久久久久成人av| 女性生殖器流出的白浆| 黑人欧美特级aaaaaa片| 两个人免费观看高清视频| 可以免费在线观看a视频的电影网站| 欧美 日韩 精品 国产| 日韩大码丰满熟妇| 国产精品久久久久久精品电影小说| 99精国产麻豆久久婷婷| 中文字幕av电影在线播放| 亚洲第一青青草原| 欧美国产精品一级二级三级| 悠悠久久av| 这个男人来自地球电影免费观看| 两性午夜刺激爽爽歪歪视频在线观看 | 在线观看免费视频日本深夜| 欧美激情高清一区二区三区| 日日夜夜操网爽| 亚洲第一青青草原| 欧美国产精品一级二级三级| 国产真人三级小视频在线观看| av片东京热男人的天堂| 好男人电影高清在线观看| 亚洲色图 男人天堂 中文字幕| 一二三四在线观看免费中文在| 男女免费视频国产| 在线观看舔阴道视频| 国产精品.久久久| 国产精品一区二区免费欧美| 在线永久观看黄色视频| 国产男女内射视频| 免费少妇av软件| 极品人妻少妇av视频| 三上悠亚av全集在线观看| 欧美成狂野欧美在线观看| 黑人巨大精品欧美一区二区mp4| 亚洲国产精品一区二区三区在线| 美国免费a级毛片| 国产在线免费精品| 久久久久久久久免费视频了| 国产av一区二区精品久久| 成人三级做爰电影| 亚洲欧美色中文字幕在线| 亚洲成a人片在线一区二区| 成人国产一区最新在线观看| 最黄视频免费看| 婷婷成人精品国产| 欧美精品一区二区免费开放| 成人特级黄色片久久久久久久 | 91九色精品人成在线观看| 亚洲av日韩精品久久久久久密| 久久人人爽av亚洲精品天堂| a级片在线免费高清观看视频| 免费久久久久久久精品成人欧美视频| 国产精品久久久av美女十八| 一本—道久久a久久精品蜜桃钙片| 99久久人妻综合| 热99re8久久精品国产| 女人被躁到高潮嗷嗷叫费观| 免费观看a级毛片全部| 在线观看www视频免费| 亚洲人成电影观看| 亚洲精品成人av观看孕妇| 久久天躁狠狠躁夜夜2o2o| 777久久人妻少妇嫩草av网站| 男女下面插进去视频免费观看| 免费日韩欧美在线观看| 51午夜福利影视在线观看| 我要看黄色一级片免费的| 99国产综合亚洲精品| 在线观看免费高清a一片| 国产精品九九99| 午夜免费鲁丝| 18禁国产床啪视频网站| 多毛熟女@视频| 久久精品aⅴ一区二区三区四区| 日韩中文字幕视频在线看片| 黑人巨大精品欧美一区二区蜜桃| 国产欧美日韩综合在线一区二区| 国产日韩欧美视频二区| 桃红色精品国产亚洲av| 国产日韩欧美视频二区| 桃红色精品国产亚洲av| 久久久精品区二区三区| 成年人黄色毛片网站| 久久久久久免费高清国产稀缺| 成人手机av| 久久久久久免费高清国产稀缺| 欧美精品一区二区免费开放| 999精品在线视频| 男人操女人黄网站| 国产精品久久久久久精品古装| 男人操女人黄网站| 久久久久久久大尺度免费视频| 午夜精品久久久久久毛片777| 飞空精品影院首页| av福利片在线| 男人舔女人的私密视频| 99久久国产精品久久久| av天堂在线播放| 亚洲免费av在线视频| 久久久国产欧美日韩av| 精品欧美一区二区三区在线| 久久精品熟女亚洲av麻豆精品| 法律面前人人平等表现在哪些方面| 人妻 亚洲 视频| 中文字幕色久视频| 日韩视频在线欧美| 多毛熟女@视频| 国产欧美日韩一区二区三| www.熟女人妻精品国产| kizo精华| 午夜免费鲁丝| 国产亚洲欧美在线一区二区| 久久精品国产综合久久久| 久久人妻福利社区极品人妻图片| 欧美 亚洲 国产 日韩一| 在线看a的网站| 国产精品久久久久久精品电影小说| 老鸭窝网址在线观看| 少妇 在线观看| 黄片大片在线免费观看| 12—13女人毛片做爰片一| 99热国产这里只有精品6| 韩国精品一区二区三区| 国产精品亚洲av一区麻豆| 亚洲第一av免费看| 国产精品 欧美亚洲| 宅男免费午夜| 精品少妇久久久久久888优播| 美女高潮到喷水免费观看| netflix在线观看网站| 中文字幕人妻丝袜制服| 激情在线观看视频在线高清 | 伦理电影免费视频| 日韩欧美三级三区| 窝窝影院91人妻| 国产男女超爽视频在线观看| 国产精品久久久久成人av| 精品福利永久在线观看| 午夜两性在线视频| 99精品久久久久人妻精品| 国产精品影院久久| 精品久久蜜臀av无| av视频免费观看在线观看| 黑人欧美特级aaaaaa片| 777久久人妻少妇嫩草av网站| 亚洲精华国产精华精| 一区二区三区激情视频| 精品国产乱码久久久久久小说| 十八禁高潮呻吟视频| 91麻豆av在线| 欧美老熟妇乱子伦牲交| 成人精品一区二区免费| 亚洲国产精品一区二区三区在线| 国产日韩欧美视频二区| 视频区图区小说| 波多野结衣av一区二区av| 国产成人啪精品午夜网站| 久久av网站| 在线观看免费高清a一片| 免费观看人在逋| 777久久人妻少妇嫩草av网站| 国产精品偷伦视频观看了| 纯流量卡能插随身wifi吗| 久久久久久亚洲精品国产蜜桃av| 人妻 亚洲 视频| 黄色a级毛片大全视频| 90打野战视频偷拍视频| 美女国产高潮福利片在线看| 91国产中文字幕| 91字幕亚洲| 亚洲免费av在线视频| h视频一区二区三区| 国产精品一区二区免费欧美| 美女国产高潮福利片在线看| 中文欧美无线码| 日韩免费av在线播放| 一本色道久久久久久精品综合| 久久精品国产亚洲av香蕉五月 | 女人被躁到高潮嗷嗷叫费观| 国产成人欧美在线观看 | 亚洲人成电影观看| 嫩草影视91久久| 日本撒尿小便嘘嘘汇集6| 精品亚洲成国产av| 亚洲精品国产色婷婷电影| 国产在线视频一区二区| 免费在线观看视频国产中文字幕亚洲| 日韩大码丰满熟妇| 日韩欧美国产一区二区入口| 免费久久久久久久精品成人欧美视频| 久久久久精品人妻al黑| 国产精品久久久久成人av| 一边摸一边抽搐一进一小说 | 黄片小视频在线播放| 黄色a级毛片大全视频| 91精品三级在线观看| 咕卡用的链子| 亚洲av国产av综合av卡| 日韩欧美国产一区二区入口| √禁漫天堂资源中文www| 99在线人妻在线中文字幕 | 久久精品国产a三级三级三级| 久久久国产一区二区| 免费在线观看完整版高清| 99久久99久久久精品蜜桃| 国产成人精品无人区| 亚洲精品一卡2卡三卡4卡5卡| 汤姆久久久久久久影院中文字幕| 免费在线观看影片大全网站| 首页视频小说图片口味搜索| 久久天堂一区二区三区四区| 午夜福利视频精品| 51午夜福利影视在线观看| av免费在线观看网站| 国产熟女午夜一区二区三区| av线在线观看网站| 成人三级做爰电影| 精品人妻在线不人妻| 精品一区二区三区av网在线观看 | 精品欧美一区二区三区在线| 国产欧美日韩一区二区三| 脱女人内裤的视频| 国内毛片毛片毛片毛片毛片| 精品少妇一区二区三区视频日本电影| 一进一出抽搐动态| 免费少妇av软件| 搡老乐熟女国产| 啦啦啦 在线观看视频| 亚洲自偷自拍图片 自拍| 久久天躁狠狠躁夜夜2o2o| 老汉色∧v一级毛片| 精品一区二区三区四区五区乱码| 国产不卡一卡二| 欧美久久黑人一区二区| 国产男女超爽视频在线观看| 国产精品1区2区在线观看. | 国产精品久久久久久人妻精品电影 | 考比视频在线观看| 啦啦啦中文免费视频观看日本| 中文字幕另类日韩欧美亚洲嫩草| 国产精品免费视频内射| 99精品久久久久人妻精品| 无遮挡黄片免费观看| 91大片在线观看| 亚洲中文av在线| 岛国在线观看网站| 欧美精品一区二区免费开放| 亚洲少妇的诱惑av| 亚洲avbb在线观看| 久久久久精品人妻al黑| 一本大道久久a久久精品| 超色免费av| 国产深夜福利视频在线观看| 国产精品欧美亚洲77777| 99在线人妻在线中文字幕 | 黄网站色视频无遮挡免费观看| 视频区图区小说| 在线观看免费视频网站a站| 精品少妇久久久久久888优播| 精品国产一区二区三区久久久樱花| 欧美精品亚洲一区二区| 美女主播在线视频| 国产又爽黄色视频| 国产99久久九九免费精品| 欧美精品一区二区大全| 韩国精品一区二区三区| 欧美大码av| 中文字幕另类日韩欧美亚洲嫩草| 亚洲精品中文字幕在线视频| 大片电影免费在线观看免费| 视频区欧美日本亚洲| 一级a爱视频在线免费观看| 日本wwww免费看| 亚洲成人国产一区在线观看| 高清欧美精品videossex| 一区在线观看完整版| 99re在线观看精品视频| 久久国产亚洲av麻豆专区| 中国美女看黄片| 亚洲成人免费av在线播放| 高潮久久久久久久久久久不卡| 在线看a的网站| 亚洲 国产 在线| 一级毛片电影观看| 亚洲色图综合在线观看| 精品福利永久在线观看| 亚洲伊人色综图| 久久人人爽av亚洲精品天堂| 国产成人精品无人区| 韩国精品一区二区三区| 丰满迷人的少妇在线观看| 狠狠狠狠99中文字幕| 亚洲全国av大片| av一本久久久久| 国产亚洲精品久久久久5区| 老司机靠b影院| 夫妻午夜视频| 久久人妻福利社区极品人妻图片| 亚洲国产欧美日韩在线播放| 黑人猛操日本美女一级片| 欧美日韩亚洲综合一区二区三区_| 欧美亚洲日本最大视频资源| 免费久久久久久久精品成人欧美视频| 精品久久久久久久毛片微露脸| 国产亚洲一区二区精品| 高清av免费在线| 国产精品自产拍在线观看55亚洲 | 亚洲精品粉嫩美女一区| 麻豆成人av在线观看| 国产单亲对白刺激| 亚洲中文av在线| 亚洲熟女毛片儿| 久久久久视频综合| 国产精品.久久久| 一区二区三区激情视频| 国产亚洲欧美精品永久| 天天添夜夜摸| 丝袜在线中文字幕| 在线播放国产精品三级| 一夜夜www| 99久久人妻综合| xxxhd国产人妻xxx| 免费看十八禁软件| 天天影视国产精品| 亚洲国产欧美网|