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

    Sparse grid-based polynomial chaos expansion for aerodynamics of an airfoil with uncertainties

    2018-05-17 10:06:20XiaojingWUWeiweiZHANGShufangSONGZhengyinYE
    CHINESE JOURNAL OF AERONAUTICS 2018年5期

    Xiaojing WU,Weiwei ZHANG,Shufang SONG,Zhengyin YE

    School of Aeronautics,Northwestern Polytechnical University,Xi’an 710072,China

    1.Introduction

    A vast amount of uncertainties exist in the practical aircraft design and application,which can cause fluctuations of aircraft performance.Therefore,it is important to take these uncertainties into account at the beginning of aircraft design.1,2Recently,many researches have concerned the topics.3–8The Uncertainty quantification(UQ)and uncertainty sensitivity analysis of aerodynamics are concerned in the paper.

    Computational Fluid Dynamics(CFD)technology has been widely used to solve problems in fluid mechanics with the development of computer technology.Traditional CFD simulation is deterministic.However,a variety of uncertainties are inevitably introduced into CFD simulation due to the increasing complexity of the problems.This leads to the mismatch between the results of CFD simulation and the actual results.The UQ in CFD simulation has gained extensive attention in Ref.9.Sources and classifications of uncertainty in CFD were described in Ref.10.Several UQ strategies have been used in CFD,including Monte Carlo Simulation(MCS)method,moment method and Polynomial Chaos(PC)in Ref.11.MCS is a statistical method,which needs many samples to accurately quantify uncertainty.Moment methods are suitable to solve the problem of small parameter uncertainty or linear model.Recently,PC which is based on the spectral representation of the uncertainty has been adopted in UQ for fluid problems in Refs.12,13.PC methods can be divided into intrusive and non-intrusive ones according to the coupling ways with CFD solvers.In general,an intrusive approach is adopted to obtain unknown polynomial coefficients by projecting resulting equations into basis functions for different modes,and it requires the modification of CFD codes,which may be difficult and time-consuming for complex problems such as Navier–Stokes equations.To overcome the shortcomings of intrusive polynomial chaos,Non-Intrusive Polynomial Chaos(NIPC)has been developed.The CFD is regarded as a black box model without changing the CFD program codes in non-intrusive methods.There are two different sampling approaches to build NIPC:random sampling and deterministic sampling.Random sampling methods use MCS to evaluate the unknown coefficients,but their convergence rate is low.Deterministic sampling methods use the quadrature to evaluate the unknown coefficients.The quadrature-based methods are more efficient than random sampling methods for low-dimensional problems.However,they are inefficient for relatively high-dimensional problems because of the exponential rising of quadrature points with the increasing dimensions.The UQ based on PC and its applications in fluid mechanics were comprehensively reviewed in Refs.14,15.

    Recently,the NIPC is sufficiently used for stochastic aerodynamic analysis with operational uncertainties.A subsonic aerodynamic analysis was conducted on a NACA0012 airfoil with an uncertain free stream velocity using a commercial flow solver in Ref.16.They proved that an uncertain free stream velocity led to the highest variation in pressure on the upper surface near the leading edge.Transonic stochastic response of two-dimensional airfoil to parameter uncertainty(Mach number Ma and angle of attack α)is focused using generalized Polynomial Chaos(gPC)in Ref.8.Two kinds of non-linearities are critical to transonic aerodynamics in their study:the shock characteristics and boundary-layer separation.A stochastic investigation of flows about NACA0012 airfoil was conducted at transonic speeds in Ref.17.A point-collocation NIPC method was used to quantify uncertainty of aerodynamic characteristics with uncertain variables Ma and α in the transonic-wing case in Ref.18.A stochastic fluid analysis on 3D wind blades considering the wind speed as an uncertain parameter was conducted in Ref.19.In their studies,when the flow separation appears,the separation vortex region corresponds to the maximum variation area which extends to the trailing edge and even to the whole suction side.However,the stochastic aerodynamic analysis considering geometric uncertainty was rarely involved.The geometric uncertainty on an aerodynamic surface resulting from manufacturing errors has significant effect on the aerodynamic performance.It is impractical to remove the impact of these geometric variations by improving the manufacturing tolerance due to the high cost of the precise surface manufacturing technique.In other words,the geometric uncertainty due to manufacturing errors is unavoidable.Therefore,it is important and necessary to conduct a stochastic aerodynamic analysis considering geometric uncertainty.Nevertheless,the description of the geometric uncertainty is difficult in a computing environment and many random variables are needed to represent aerodynamic surface fluctuations.Therefore,the quadrature-based NIPC method is inefficient in solving this high-dimensional stochastic problem.

    To improve the computational efficiency of traditional quadrature-based NIPC method for high-dimensional problems,the sparse grid numerical integration is introduced to solve the coefficients of PC.The sparse grid technique extends one-dimensional formulae to higher dimensions by tensor product and then selects sampling points according to Smolyak theory in Ref.20.It has been widely used in numerical integration and interpolation21,22as well as data mining.23Recently,a new sparse grid-based method has been developed for UQ in Ref.24.From their research,it can be known that when the dimension of the uncertain variables is larger than 5,the computational cost required by the sparse grid method is much smaller than that required by tensor product method.A sparse grid interpolant was developed to solve the high-dimensional stochastic partial different equations in Ref.25.Sparse grid collocation schemes were applied to stochastic natural convection problems in Ref.26.Sparse grids-based stochastic approximations with applications to subsonic steady flow about a NACA0015 airfoil in the presence of geometrical and operational uncertainties with both simplified aerodynamics model and Reynolds-Averaged Navier–Stokes(RANS)simulation was presented in Ref.27.

    In this paper,a Sparse Grid-based Polynomial Chaos(SGPC)method is constructed to UQ and sensitivity analysis for transonic aerodynamics considering airfoil geometric and operational uncertainties in detail.The paper is structured as follows:Section 2 introduces the mathematical formulations of the sparse grid technique;in Section 3,the SGPC is built;in Section 4,a stochastic aerodynamic analysis considering geometric and operational uncertainties is conducted in detail;Section 5 outlines several useful conclusions.

    2.Sparse grid numerical integration

    Sparse grid technique selects sampling points under Smolyak theory,which uses a weighted linear combination of special tensor products to reduce the grid size.The sparse grid has been successfully used for numerical integration.Locations and weights of the univariate quadrature points with a range of accuracy are determined in each dimension,and then the univariate quadrature point sets are extended to form a multi-dimensional grid using the sparse grid theory.The introduction of Sparse Grid Numerical Integration(SGNI)is as follows:

    where?denotes the tensor product rule.

    Quadrature points of multivariate integrals are all possible combinations of one-dimensional quadrature nodes(Nfull=nm1nm2···nmn).Thus,it is time-consuming for a problem with relatively high dimensions by the tensor product algorithm.

    Since the full tensor product is inefficient for high dimensions,Smolyak theory is adopted because it can reduce the grid size with a weight linear combination of special tensor products.For the n-dimensional sampling points with p-level accuracy generated by the sparse grid technique,the tensor product rule based on the sparse grid is shown as

    where|m|=m1+m2+ ···+mn.The bounded sum makes sure that tensor products exclude from full grids points that contribute less to the improvement of the required integration accuracy.The integration of multivariate nonlinear functionin terms of variable ξ by sparse grid can be computed by

    where Psis the sum of all possible combinations of the multiple indices that satisfy q=p+n,n is the number of random dimensions,and ξ is the selected collocation points by sparse grid method.The weight alcorresponding to the lth collocation point is computed by

    Table 1 shows the change of sparse grid collocation points(Nsg)from 2 to 10 dimensions with the 2-order accuracy Gaussian quadrature(Ntp).In a low-dimensional problem,with the same level p,the sparse grid approach needs more collocation points than the full grid approach.When n increases(e.g.n≥5),the number of collocation points of the sparse grid method is much smaller than that of the full grid.It can be observed that with the same polynomial exactness p the number of collocation points produced by the full grid scheme is(p+1)nwhich increases exponentially with dimensionality,whereas the sparse grid technique remarkably reduces the number of collocation points.

    3.Sparse grid-based polynomial chaos

    PC has been widely used in UQ,but the method is inefficient for high-dimensional problems.Therefore,the SGPC method is constructed to alleviate the computational burden for relatively high-dimensional problems.The method is introduced in detail in this section.

    PC is a stochastic method based on the spectral representation of uncertainty.According to the spectral representation,the random function can be decomposed into deterministic and stochastic components.For example,a random variable(X)can be represented by

    where αj(x)is the deterministic component and Ψj(ξ)is the random basis function corresponding to the jth mode.From Eq.(6),the random variable X is the function of the deterministic independent variable vector x and the n-dimensional standard random variable vector ξ =(ξ1,ξ2,···,ξn).PC expansion given by Eq.(6)should contain an in finite number of terms.In a practical computational context,PC is truncated in both order p and dimension n.The number of truncated terms is finite as follows:

    The random basis function Ψjis chosen according to the type of the input random variable.For example,if the input uncertainty obeys Gauss distribution,the basis function is the multidimensional Hermite polynomial.If ξ is chosen to be uniform with the random variable,the basis function must be the Legendre polynomial.A complete description of gPC scheme is introduced in Ref.13.

    The purpose of the PC method is to obtain unknown polynomial coefficients.Eq.(6)can be transformed to Eq.(8)by the inner product:

    Table 1 Number of collocation points with different dimensions by 2-order accuracy of Gaussian quadrature(tensor product and sparse grid).

    〈·〉represents the inner product which can be expressed by

    Because of orthogonality,Eq.(8)can be transformed to

    And then it can be derived:

    There are two ways to solve the coefficients:intrusive and non-intrusive methods.The intrusive method computes unknown polynomial coefficients by projecting the resulting equations into basis functions for different modes.It requires the modification of the deterministic codes,which may be difficult,expensive,and time-consuming for complex computational problems.The non-intrusive method treats the CFD as a black box without changing the program code when propagating uncertainty.Steps of the non-intrusive method are as follows:

    Step 1.Adopting relevant sampling methods to produce sample

    Step 2.For each sample ξj,evaluate the uncertainty parameter λj.

    Step 3.Use the selected N samples to determine the expansion coefficients by Galerkin projection.

    Step 4.Given the computed coefficient αk,a polynomial approximation model can be built∑

    Sampling approaches can be divided into random and deterministic sampling methods in Step 1.Random sampling methods use MCS to compute projection integrals,where the convergence rate is low.Deterministic sampling methods use the numerical quadrature to evaluate unknown coefficients.Using n-dimensional Gauss quadrature with q points in each dimension,the unknown coefficients can be computed by

    One-dimensional integral is expanded into a highdimensional integral by tensor product in Eq.(12),the calculation times of which are qn(requiring(p+1)npoints for pth order chaos).For low-dimensional problems,the efficiency of the deterministic sampling method has been greatly improved compared with that of the random sampling method.However,the computational cost grows exponentially with the increasing dimensions.Thus,the NIPC method is inefficient in high-dimensional problems.Such shortcomings motivate us to improve NIPC method.If we use much less integral points than the conventional tensor product to solve the multi-dimensional integral Eq.(8),the computational cost of UQ can be reduced.It is well known that the SGNI can improve the computational efficiency to solve the multidimensional integral Eq.(8).From the description in Section 2,the sparse grid method is more efficient than the tensor product for high-dimensional integration(e.g.n≥5).Hence,we use the SGNI to replace the tensor product to solve Eq.(8)as follows:

    4.Stochastic aerodynamic analysis

    A stochastic aerodynamic analysis is conducted considering airfoil geometric and operational uncertainties.The operational uncertainty generally contains two parameters:Ma and α.However,the relatively high-dimensional random parameters are needed to represent aerodynamic surface fluctuations.In this section,SGPC and MCS methods are applied to the stochastic transonic aerodynamics analysis around a RAE2822 airfoil.Firstly,the deterministic flow computations are briefly introduced.

    4.1.Deterministic flow computations

    The deterministic steady-state flow solutions are computed by the RANS equations or method combined with the Spalart–Allmaras turbulence model.The cell-centered finite volume method is used for spatial discretization and AUSM-up scheme is used to evaluate the numerical flux.The implicit dual-time stepping method is used for temporal discretization.The symmetric Gauss–Seidel iterative time-marching scheme is applied in the pseudo-time step,and the second order accurate full implicit scheme is used to solve the equation in the physical time step.The developed flow solver is also used for unsteady flows.28,29

    The computational mesh surrounding a RAE2822 airfoil is based on a structured C-grid.The 2D mesh is composed of two blocks,the sizes of which are 300×150(C block surrounding the airfoil)and 200×150.The chord of the airfoil is c=1 m and the far- field boundary is placed at a distance d=20c from the airfoil.Fig.1 shows the computing grids of RAE2822 airfoil and the reliability verification of the CFD program,in which Cpis the pressure coefficient.The geometry of airfoil keeps changing in the process of UQ,and Radial Basis Function(RBF)is used to mesh deformation.30

    4.2.Geometric uncertainty

    In aerospace engineering,despite advances of manufacturing engineering techniques,airfoilsvery often exhibitsome deviation from their intended shape due to the noisy manufacturing process.Moreover,it is often impractical to remove the impact of these geometric variations by improving the manufacturing tolerance due to the high cost of precise surface manufacturing techniques.In other words,the geometric uncertainty resulting from manufacturing errors is unavoidable.How to describe this geometric uncertainty in computing environment should be concerned.

    Main geometric variation modes are obtained by Principal Component Analysis(PCA)with a large amount of geometric statistical measurement data of an airfoil in Refs.31,32.The description of the airfoil geometric uncertainty by PCA is shown as

    where gnis the nominal geometry;ˉg is the average geometric variation;viis the geometric mode shape;nsis the number of mode shapes used to represent the variation in geometry.The geometric mode can be computed by PCA based on measurement samples.σiis the ith singular value of the measurement snapshot matrix,which represents the geometric variability attributable to the ith mode.ziis a random parameter which obeys the standard normal distribution,and thus the product σiziis the stochastic contribution of the ith mode.

    It is difficult to describe this geometric variation in the computing environment.A Gaussian random process simulation is used to obtain the geometric data in Ref.33.Then PCA is used to obtain main geometric modes.In addition,a variety of parametric methods have been employed to describe the geometric variation in aerodynamic design so far,such as PARSEC-11 geometry parameterization,34Class-Shape function Transformation(CST)method,35Free-Form Deformation(FFD)method,36and Hicks-Henne bump functions.37,38By changing the parameter of these parametric methods,the geometric variation is realized.Generally,these parametric methods need a lot of parameters to represent the airfoil shape.To reduce the dimensions of the variables,the PCA technology combined with airfoil parameterization is used to describe the geometric variation39,40used in the paper.Firstly,a parametric method is used to generate a set of sample data,and then the PCA is used to obtain the main deformation modes based on the generated sample data.In this paper,we use a parameterized representation of the airfoil ARE2822 by CST method with 12 parameters.The data of the measurement points on the airfoil surface are obtained by random perturbation of CST parameters.In this way,the PCA based sample data are conducted.Fig.2 shows the eigenvalue with the number of modes, which represents the geometric variation attributable to each mode.The smaller the eigenvalue of one mode is,the smaller the proportion of the mode to geometric variation is.The first 6 modes obtained by PCA are showed in Fig.3.It can be observed that these modes are global geometric deformation modes and present some typical geometric deformation.specifically,Mode 1 and Mode 2 are the scale modes in the thickness direction;Mode 3 and Mode 4 are translation modes of the maximum thickness in the axial direction;Mode 5 and Mode 6 are the extrusion modes of the upper surface.

    4.3.Stochastic transonic aerodynamics analysis

    In the paper,the uncertainties are specified by means of PDFs.The input uncertainties,which PDFs should obey,accord to the corresponding statistical information,which are found in the methods of statistical inference.In the current study,the uncertain geometric variables are zi(i=1,2,...,6)in Eq.(14)which obey the standard normal distribution according to Ref.26.The operational uncertainties(Ma and α)obey truncated Gaussian distribution.The Mach number has a 0.73 mean value and a±0.01 variability,and the angle of attack has a 2.5°mean value and a± 0.3°variability.The steady flow state is selected in a transonic region(Ma=0.73,α=2.5°,Re=3.0× 106),and the range of geometric uncertainty is showed in Fig.4.

    4.3.1.Verification of SGPC

    Table 2 Results of uncertainty quantification of aerodynamic coefficients.

    Table 3 Comparison of CFD calls among selected methods.

    A convergence study of SGPC has been performed by the stochastic transonic aerodynamic analysis.The polynomial order p is increased to enhance the accuracy of SGPC.Moreover,MCS is also used in this analysis,which aims to validate the accuracy and efficiency of SGPC.The number of CFD calls of MCS is 5000.

    Figs.5 and 6 show the standard deviations of the pressure coefficient Cpdistribution and skin-friction coefficient Cfdistribution along the chord and they also re flect the convergence of SGPC method.It can be seen that the standard deviation distributions of Cpand Cfcoincide well with the results of MCS when SGPC reaches the 3-order accuracy.Table 2 also gives the convergence of SGPC for lift coefficient CLand drag coefficient CD.It can be observed that when p reaches 2-order accuracy,the Standard Deviation(StD)of aerodynamic coefficients is convergent and coincides well with the results of MCS.The convergence study of SGPC verifies the accuracy and correctness of SGPC to the stochastic aerodynamic analysis.

    The computational cost of SGPC should also be concerned for the stochastic transonic aerodynamic analysis.From Section 3,it can be learned that the computational cost of SGPC is related with parameters n and p.Table 3 shows computational costs of SGPC,NIPC and MCS for the 8-dimensional problem.It can be observed that SGPC is more efficient than the other two.It must be emphasized here that NIPC method is not used to the stochastic aerodynamic analysis because of the prohibitive computational cost.The CFD calls of NIPC is obtained according to the formula(p+1)n.

    4.3.2.Uncertainty quantification of transonic aerodynamics

    Now,the in fluence of input uncertainties on transonic flows is studied.Fig.7 shows the mean and fluctuations of Cpon the airfoil surface.It is well known that the shock wave exists in the transonic flow.The flow characteristics before and after the shock wave change dramatically with the input uncertainties.Therefore,it can be observed that a strong spurious effect appears near the shock wave foot.Fig.8 shows the mean and fluctuation of Cfon the airfoil surface.The skin-friction behavior displays discrepancies from the pressure coefficient distribution.The variation of Cfis similar to that of Cpbefore and in the shock disturbance region.However,the difference occurs in the downstream region of shock wave;the standard deviation has larger magnitudes below the shock disturbance region,and then it gradually decreases.This indicates that the flow characteristics not only near the shock wave but also in the boundary-layer separation are sensitive to geometric and operational uncertainties.

    Moreover,Figs.7 and 8 show the corresponding local PDF pro files at the selected five chord locations.It is well known that once the mean position of the shock wave is constrained by the separated shear layer,the probability density function of the solution may exhibit a bifurcation corresponding to a jump in the solution.Therefore,the stochastic response shows a bimodal nature at x/c=0.56.

    4.3.3.Global sensitivity analysis of transonic aerodynamics

    Global Sensitivity Analysis(GSA)should be applied to study how uncertainty of the model output can be apportioned to different sources of uncertainty of the model input.A sensitivity analysis is performed by using the Sobol’s analysis,which analyzes the individual and coupled effects between the random parameters(a deep introduction of the Sobol’s analysis can be referred to in Ref.41).When a PC model is built,the Sobol’s analysis based on the PC can be directly calculated from the expansion coefficients.41In the paper,the aerodynamic load distribution and aerodynamic coefficients are included for Sobol’s analysis.

    Figs.9 and 10 show the results of Sobol’s analysis of Cpand Cf.The partial standard deviation and standard deviation caused by the coupling effect for each uncertainty variable are given.It can be observed that Ma is important to the flow characteristics both in the vicinity of the shock wave and in the boundary-layer separation.Now,we concern the contribution of the geometric deformation modes.It can be observed that deformation modes of the upper surface are much important than those of the lower surface to aerodynamic load distribution.Specifically,the scale mode(Mode 1)and extrusion mode(Mode 5)are important to the flow characteristic of shock wave and boundary-layer separation,but the translation mode(Mode 3)is only important to the flow characteristic of shock wave.Moreover,it can also be observed that the coupling effect among geometric deformation modes is not negligible near the shock wave.

    Figs.11 and 12 show the results of Sobol’s analysis of aerodynamic coefficients.It can be observed that α is the most important to lift coefficients.To drag characteristic,Ma is the most important;geometric modes of the upper surface and angle of attack also have obvious influences.

    4.4.Stochastic aerodynamic investigation considering geometric uncertainty in different flow regimes

    A further stochastic aerodynamic investigation is conducted in different flow regimes.Three flow states are selected:Case A,subsonic flow with a small angle of attack(Ma=0.5,α=2.5°);Case B,subsonic flow with a large angle of attack(Ma=0.5,α =20°);and Case C,transonic flow with a small angle of attack(Ma=0.73, α =2.5°).SGPC with 3-order accuracy is used in the stochastic aerodynamic analysis.At first,the three kinds of flow states have been employed in the steady flow computation,respectively.The streamline of the flow field is shown in Fig.13.It can be observed that the flow is smooth and there is no separation along the airfoil in Case A.There is a separation vortex above the upper surface of the airfoil in Case B.In Case C,the separation appears on the aft part of the airfoil because of the boundary layer separation.

    Fig.14 shows error bars of pressure and friction coefficient distribution of RAE2822 airfoil in the three states.In case A,variations of Cpand Cfare the same,and the maximum fluctuations mainly lie in the upper surface near the leading edge.In Case B,we can also observe that the maximum fluctuations mainly display on the upper surface near the leading edge.Fluctuations in other regions are even smaller than those in Case A.This indicates that flow characteristics in the separation region are intensive to this geometric uncertainty.In Case C,variation regularities of Cpand Cfare different.For Cp,the variation mainly appears in the shock wave disturbance region,and for Cf,besides this region,it also exists in the boundary layer separation region.

    Table 4 shows uncertainty analysis results of aerodynamic coefficients.The Cov of CDis relatively small for Case A and Case B while it is relatively large for Case C.The Cov of CLis small for all the three cases.This means that subsonic aerodynamic coefficients are not sensitive to the geometric uncertainty.The drag is sensitive to the geometric uncertainty,but the lift is not in the transonic region.

    Table 4 Uncertainty quantification of aerodynamic coefficients in three different flow states.

    5.Conclusions

    In the paper,the sparse grid-based polynomial chaos expansion was used for UQ and GSA for transonic aerodynamics with geometric and operational uncertainties.The method is much more efficient than NIPC and MCS in dealing with relatively high-dimensional stochastic aerodynamic problems.The accuracy and efficiency of the SGPC are verified by the stochastic transonic aerodynamic analysis.It is proved that the method is more efficient than NIPC and MCS.It is observed that the flow characteristics in the vicinity of the shock wave and the boundary-layer separation region are sensitive to uncertainties because of the nonlinear flow characteristics of shock wave and boundary-layer separation.By uncertainty sensitivity analysis(Sobol’s analysis),it can be learnt that Ma is important to the flow characteristics both in the vicinity of the shock wave and in the boundary-layer separation;the geometric deformation modes of the upper surface are much important than those of the lower surface to aerodynamic load distribution.Specifically,the scale mode(Mode 1)and extrusion mode(Mode 5)are important to the flow characteristic of shock wave and boundary-layer separation,but the translation mode(Mode 3)is only important to the flow characteristic of shock wave.The coupling effect is not negligible in the vicinity of the shock wave foot.

    Moreover,a stochastic aerodynamic investigation considering geometric uncertainty is conducted by the SGPC in different flow states.It can be concluded that fluctuations of subsonic aerodynamic characteristics are mainly concentrated in the airfoil head.When the flow is in the transonic condition,the affected region shifts from the upper surface near the leading edge to the shock wave and boundary-layer separation region because of the nonlinear characteristics of the shock wave and boundary-layer separation behind the shock wave.

    Acknowledgements

    This work was supported by the National Natural Science Foundation of China(No.11572252),the ‘111” Project of China(No.B17037),and the National Science Fund for Excellent Young Scholars(No.11622220).

    References

    1.Zang TA,Hemsch MJ,Hilburger MW,Kenny SP,Luckring JM,Maghami P,et al.Needs and opportunities for uncertainty-based multidisciplinary design methods for aerospace vehicles.Washington,D.C.:NASA;2002.Report No.:NASA/TM-2002-211462.

    2.Yao W,Chen XQ,Luo WC,Tooren MV,Guo J.Review of uncertainty-based multidisciplinary design optimization methods for aerospace vehicles.Prog Aerosp sci 2011;47(6):450–79.

    3.Song SF,Lu ZZ,Zhang WW,Ye ZY.Reliability and sensitivity analysis of transonic flutter using improved line sampling technique.Chin J Aeronaut 2009;22(5):513–9.

    4.Tang J,Wu Z,Yang C.Epistemic uncertainty quantification in lf utter analysis using evidence theory.Chin J Aeronaut 2015;28(1):164–71.

    5.Zhu H,Tian H,Cai GB,Bao WM.Uncertainty analysis and design optimization of hybrid rocket motor powered vehicle for suborbital flight.Chin J Aeronaut 2015;8(3):676–86.

    6.Dai Y,Yang C.Methods and advances in the study of aeroelasticity with uncertainties.Chin J Aeronaut 2014;27(3):461–74.

    7.Pagnacco E,Souza DCE,Sampaio R.Subspace inverse power method and polynomial chaos representation for the modal frequency responses of random mechanical systems.Comput Mech 2016;58:129–49.

    8.Simon F,Guillen P,Sagautn P,Lucor D.A gPC-based approach to uncertain transonic aerodynamics.Comput Method Appl Mech Eng 2010;199:1091–9.

    9.Roache PJ.Quanti fication of uncertainty in computational fluid dynamics.Annu Rev Fluid Mech 1997;29:123–60.

    10.Pelletier D,Turgeon E,Lacasse D.Adaptivity,sensitivity,and uncertainty:Toward standards of good practice in computational fluid dynamics.AIAA J 2003;41(10):1925–32.

    11.Walters RW,Huyse L.Uncertainty analysis for fluid mechanics with applications.Washington,D.C.:NASA;2002.Report No.:NACA/CR-2002-211449.

    12.Mathelin L,Hussaini MY,Zang TZ.Stochastic approaches to uncertainty quantification in CFD simulations.Numer Algorithms 2005;38(1):209–36.

    13.Xiu DB,Karniadakis GE.Modeling uncertainty in flow simulationsvia generalized polynomialchaos.J ComputPhys 2003;187:137–67.

    14.Knio OM,Maitre OPL.Uncertainty propagation in CFD using polynomial chaos decomposition.Fluid Dyn Res 2006;38:616–40.

    15.Najm HN.UQ and polynomial chaos techniques in computational lf uid dynamics.Annu Rev Fluid Mech 2009;41:35–52.

    16.Loeven GJA,Witteveen JAS,Bijl H.Probabilistic collocation:An efficient non-intrusive approach for arbitrarily distributed parametric uncertainties.45th AIAA aerospace sciences meeting and exhibit;2007 Jan 8-11;Reno,Nevada.Reston:AIAA;2007.

    17.Chasseing JC,Lucor D.Stochastic investigation of flows about airfoils at transonic speeds.AIAA J 2010;48(5):918–49.

    18.Hosder S,Walters RW,Balch M.Point-collocation nonintrusive polynomial chaos method for stochastic computational fluid dynamics.AIAA J 2012;48(12):2721–30.

    19.Liu ZY,Wang XD,Shun K.Stochastic performance evaluation of horizontal axis wind turbine blades using non-deterministic CFD simulations.Energy 2014;73:126–36.

    20.Bungartz HJ,GriebelM.Sparse grids.Acta Numerica 2004;13:147–269.

    21.Thomas G,Michael G.Numerical integration using sparse grids.Numer Algorithms 1998;18:209–32.

    22.Novak E,Ritter K.High dimensional integration of smooth functions over cubes.J Numer Math 1996;75:79–97.

    23.Garchke J,Griebel M,Thess M.Data mining with sparse grids.Computing 2001;67:225–53.

    24.Xiong FF,Greene S,Chen W.A new sparse grid based method for uncertainty propagation.Struct Multidisc Optim 2010;41:335–49.

    25.Xiu D,Hesthaven J.Higher-order collocation method for differential equations with random inputs.SIAM J Sci Comput 2005;27(3):1118–39.

    26.Baskar G,Nicholas Z.Sparse grid collocation schemes for stochastic natural convection problems.J Comput phys 2007;225:652–85.

    27.Resmini A,Peter J,Lucor D.Sparse grids-based stochastic approximations with applications to aerodynamics sensitivity analysis.Int J Numer Meth Eng 2016;1(106):32–57.

    28.Zhang WW,Li XT,Ye ZY,Jiang YW.Mechanism of frequency lock-in in vortex-induced vibrations at low Reynolds numbers.J Fluid Mech 2015;783:72–102.

    29.Zhang WW,Gao CQ,Liu YL,Ye ZY,Jiang YW.The interaction between flutter and buffet in transonic flow.Nonlinear Dynam 2015;82(4):1851–65.

    30.Jakobsson S,Amoignon O.Mesh deformation using radial basis functions for gradient-based aerodynamic shape optimization.Comput Fluids 2007;36:1119–36.

    31.Garzon V,Darmofal D.Impact of geometric variability on axial compressor performance.J Turbomac 2003;125(4):692–703.

    32.Thanh TB,Willcox K.Parametric reduced-order models for probabilistic analysis of unsteady aerodynamic applications.AIAA J 2008;46(10):2520–9.

    33.Chen H,Wang QQ,Hu R,Paul C.Conditional sampling and experiment design for quantifying manufacturing error of transonic airfoil.49th AIAA aerospace sciences meeting including the new horizons forum and aerospace exposition;2011 Jan 4–7;Orlando,USA.Reston:AIAA;2011.

    34.Dodson M,Parks GT.Robust aerodynamic design optimization using polynomial chaos.J Aircraft 2009;46(2):635–45.

    35.Kulfan B,Bussoletti J. ‘Fundamental” parameteric geometry representations for aircraft component shapes.11th AIAA/ISSMO multidisciplinary analysis and optimization conference;2006 Sep 6–8;Portsmouth,USA;Reston:AIAA;2006.p.1–45.

    36.Padulo M,Maginot J,Guenov M.Airfoil design under uncertainty with robust geometric parameterization.50th AIAA/ASME/ASCE/AHS/ASC structures,structural dynamics,and materials conference;2009 May 4–7;Palm Springs,USA.Reston:AIAA;2009.

    37.Duan YH,Cai JS,Li YZ.Gappy proper orthogonal decomposition-based two-step optimization for airfoil design.AIAA J 2012;50(4):968–71.

    38.Padulo M,Campobasso MS,Guenov MD.Novel uncertainty propagation method for robust aerodynamic design.AIAA J 2012;49(3):530–43.

    39.Wu XJ,Zhang WW,Song SF.Uncertainty quantification and sensitivity analysis of transonic aerodynamics with geometric uncertainty.Int J Aerospace Eng 2017;8107190:1–16.

    40.Wei PF,Lu ZZ,Song JW.Variable importance analysis:A comprehensive review.Reliab Eng Syst Saf 2015;42:399–432.

    41.Bruno S.Global sensitivity analysis using polynomial chaos expansions.Reliab Eng Syst Saf 2008;93:964–79.

    久久精品国产99精品国产亚洲性色| 亚洲在线自拍视频| 在线观看舔阴道视频| 在线观看www视频免费| 搡老岳熟女国产| 国产一区二区激情短视频| 一区二区三区精品91| 亚洲成国产人片在线观看| 窝窝影院91人妻| 日本成人三级电影网站| 亚洲一区二区三区色噜噜| 亚洲免费av在线视频| 1024视频免费在线观看| 十八禁人妻一区二区| 18禁黄网站禁片午夜丰满| x7x7x7水蜜桃| 久久草成人影院| 午夜福利成人在线免费观看| 正在播放国产对白刺激| 我的亚洲天堂| 欧美乱色亚洲激情| 婷婷亚洲欧美| 国产精品亚洲一级av第二区| 欧美三级亚洲精品| 一二三四在线观看免费中文在| 亚洲成人免费电影在线观看| ponron亚洲| 久久久久久免费高清国产稀缺| 一区二区三区激情视频| 久久久久久九九精品二区国产 | 久久精品aⅴ一区二区三区四区| 国产一区二区三区在线臀色熟女| 美女 人体艺术 gogo| 国产片内射在线| 精品久久久久久成人av| 在线看三级毛片| 成年女人毛片免费观看观看9| 久久这里只有精品19| 午夜福利在线观看吧| 一边摸一边抽搐一进一小说| 在线观看免费日韩欧美大片| 亚洲五月婷婷丁香| 成年人黄色毛片网站| 18美女黄网站色大片免费观看| 久久久久九九精品影院| 50天的宝宝边吃奶边哭怎么回事| 18禁观看日本| 日韩欧美免费精品| 国产激情久久老熟女| 女性被躁到高潮视频| 欧美av亚洲av综合av国产av| 久久亚洲真实| 人人澡人人妻人| 亚洲三区欧美一区| 国产精品香港三级国产av潘金莲| 18禁国产床啪视频网站| 搡老岳熟女国产| 国产又色又爽无遮挡免费看| 亚洲自拍偷在线| 97碰自拍视频| 一进一出抽搐gif免费好疼| 最好的美女福利视频网| 宅男免费午夜| 国产av不卡久久| 欧美黄色片欧美黄色片| 丰满的人妻完整版| 女性生殖器流出的白浆| 欧美zozozo另类| 成人手机av| aaaaa片日本免费| 成人18禁在线播放| 亚洲av中文字字幕乱码综合 | 91成年电影在线观看| 国产三级在线视频| 亚洲片人在线观看| 国产一区二区激情短视频| 亚洲黑人精品在线| 99热6这里只有精品| 欧美另类亚洲清纯唯美| 日日爽夜夜爽网站| 波多野结衣高清无吗| 在线天堂中文资源库| 久久久国产成人免费| 国产激情欧美一区二区| 此物有八面人人有两片| 国产高清有码在线观看视频 | 久久人妻av系列| 在线观看www视频免费| 一二三四社区在线视频社区8| 麻豆国产av国片精品| 精品国内亚洲2022精品成人| 大型黄色视频在线免费观看| 51午夜福利影视在线观看| 十分钟在线观看高清视频www| 成人永久免费在线观看视频| 日韩欧美一区二区三区在线观看| 日韩高清综合在线| 国产激情偷乱视频一区二区| 美女午夜性视频免费| 亚洲国产精品成人综合色| 精品久久久久久成人av| 熟女电影av网| 久久久久久大精品| √禁漫天堂资源中文www| 变态另类成人亚洲欧美熟女| 亚洲精品国产精品久久久不卡| 无人区码免费观看不卡| 亚洲午夜精品一区,二区,三区| 国产主播在线观看一区二区| 欧美日韩一级在线毛片| 精品久久久久久久人妻蜜臀av| 午夜免费成人在线视频| 亚洲精华国产精华精| 男女那种视频在线观看| 最好的美女福利视频网| 国产精品久久视频播放| 精品久久久久久久人妻蜜臀av| 日本五十路高清| 性色av乱码一区二区三区2| 欧美日韩亚洲国产一区二区在线观看| 欧美成人午夜精品| 日韩视频一区二区在线观看| 亚洲第一电影网av| 国产一级毛片七仙女欲春2 | 麻豆成人av在线观看| 亚洲男人天堂网一区| 久9热在线精品视频| 精品久久久久久久人妻蜜臀av| 国产亚洲精品第一综合不卡| 99riav亚洲国产免费| 午夜免费鲁丝| 女人高潮潮喷娇喘18禁视频| 高清毛片免费观看视频网站| 热re99久久国产66热| 国产精品久久久av美女十八| 日韩高清综合在线| 女人爽到高潮嗷嗷叫在线视频| 亚洲第一av免费看| 1024手机看黄色片| 99国产精品一区二区三区| 国产精品久久久久久精品电影 | 91成人精品电影| 亚洲av日韩精品久久久久久密| 亚洲精品久久国产高清桃花| 99精品久久久久人妻精品| 美女扒开内裤让男人捅视频| 国产激情欧美一区二区| 中文亚洲av片在线观看爽| 国产高清视频在线播放一区| 亚洲精品久久成人aⅴ小说| 久久久国产精品麻豆| 视频在线观看一区二区三区| 十分钟在线观看高清视频www| 中文字幕av电影在线播放| 成年版毛片免费区| 99久久久亚洲精品蜜臀av| 18禁国产床啪视频网站| 成人午夜高清在线视频 | 99精品欧美一区二区三区四区| 国产麻豆成人av免费视频| 无限看片的www在线观看| 少妇被粗大的猛进出69影院| 亚洲欧美激情综合另类| 国产精品一区二区精品视频观看| 视频区欧美日本亚洲| 国产视频内射| 亚洲午夜理论影院| 国产三级在线视频| 狠狠狠狠99中文字幕| 特大巨黑吊av在线直播 | 大型av网站在线播放| 欧美成人免费av一区二区三区| 757午夜福利合集在线观看| 亚洲国产精品sss在线观看| av视频在线观看入口| 成人18禁在线播放| 变态另类成人亚洲欧美熟女| 88av欧美| 俺也久久电影网| 精品第一国产精品| 欧美日韩黄片免| 久久99热这里只有精品18| 久久精品国产亚洲av高清一级| 热99re8久久精品国产| 老司机福利观看| 国产免费男女视频| 国产精品影院久久| 亚洲精品美女久久久久99蜜臀| 色精品久久人妻99蜜桃| 熟女少妇亚洲综合色aaa.| 国产熟女xx| 亚洲专区中文字幕在线| 亚洲国产欧洲综合997久久, | 变态另类成人亚洲欧美熟女| 色在线成人网| a级毛片a级免费在线| 99精品久久久久人妻精品| 高潮久久久久久久久久久不卡| 免费电影在线观看免费观看| 日韩欧美国产在线观看| 国产人伦9x9x在线观看| 精品国产亚洲在线| 日本精品一区二区三区蜜桃| 欧美zozozo另类| 精品久久蜜臀av无| 黑人操中国人逼视频| 亚洲国产毛片av蜜桃av| 精品免费久久久久久久清纯| 91麻豆av在线| 美国免费a级毛片| 免费搜索国产男女视频| 久久久久久大精品| 亚洲一区高清亚洲精品| av电影中文网址| 最好的美女福利视频网| 淫秽高清视频在线观看| 99在线视频只有这里精品首页| 淫妇啪啪啪对白视频| 一级黄色大片毛片| 亚洲熟妇中文字幕五十中出| 欧美人与性动交α欧美精品济南到| 窝窝影院91人妻| 久久精品人妻少妇| av视频在线观看入口| 人妻久久中文字幕网| 搡老熟女国产l中国老女人| 手机成人av网站| 亚洲av成人不卡在线观看播放网| 满18在线观看网站| 亚洲av成人一区二区三| 久久99热这里只有精品18| 很黄的视频免费| 亚洲五月色婷婷综合| 黄色视频,在线免费观看| 欧美大码av| 亚洲真实伦在线观看| 午夜福利视频1000在线观看| av免费在线观看网站| 亚洲性夜色夜夜综合| 制服人妻中文乱码| 1024视频免费在线观看| 这个男人来自地球电影免费观看| 亚洲片人在线观看| 香蕉丝袜av| 99精品久久久久人妻精品| 免费观看精品视频网站| 在线观看www视频免费| 午夜福利在线在线| 国产免费av片在线观看野外av| 十八禁网站免费在线| 亚洲欧美一区二区三区黑人| 99久久综合精品五月天人人| 久久久久免费精品人妻一区二区 | 草草在线视频免费看| 一本久久中文字幕| 91老司机精品| 欧美 亚洲 国产 日韩一| 女人高潮潮喷娇喘18禁视频| 少妇裸体淫交视频免费看高清 | 97人妻精品一区二区三区麻豆 | 国产精品久久久久久人妻精品电影| 欧美一级毛片孕妇| 亚洲自偷自拍图片 自拍| 亚洲自偷自拍图片 自拍| 色综合婷婷激情| 人人澡人人妻人| 丰满人妻熟妇乱又伦精品不卡| 在线av久久热| 欧美日韩福利视频一区二区| 欧美成人一区二区免费高清观看 | 老司机午夜福利在线观看视频| av福利片在线| 午夜福利视频1000在线观看| 最近最新中文字幕大全免费视频| 国产精品电影一区二区三区| 婷婷丁香在线五月| 免费在线观看亚洲国产| 国产av又大| 这个男人来自地球电影免费观看| 18禁美女被吸乳视频| 久久久国产精品麻豆| 99精品欧美一区二区三区四区| 久热这里只有精品99| 丁香欧美五月| www.自偷自拍.com| 丝袜在线中文字幕| 在线观看免费视频日本深夜| 欧美日韩瑟瑟在线播放| 亚洲av片天天在线观看| 国内精品久久久久精免费| 中文在线观看免费www的网站 | 嫁个100分男人电影在线观看| 99国产精品一区二区蜜桃av| 丝袜美腿诱惑在线| 亚洲无线在线观看| 免费高清视频大片| 精品乱码久久久久久99久播| 草草在线视频免费看| 日韩 欧美 亚洲 中文字幕| 精品人妻1区二区| 亚洲全国av大片| 身体一侧抽搐| 妹子高潮喷水视频| 一级黄色大片毛片| 国产精品一区二区精品视频观看| 少妇熟女aⅴ在线视频| 久久精品亚洲精品国产色婷小说| 在线国产一区二区在线| 成人一区二区视频在线观看| 淫妇啪啪啪对白视频| 男女那种视频在线观看| 久久热在线av| 最近最新中文字幕大全电影3 | 成人精品一区二区免费| 欧美最黄视频在线播放免费| 日韩大尺度精品在线看网址| 高潮久久久久久久久久久不卡| 男女那种视频在线观看| 黑人欧美特级aaaaaa片| 黄片播放在线免费| 两个人看的免费小视频| 九色国产91popny在线| 最新美女视频免费是黄的| 黑人欧美特级aaaaaa片| 人妻丰满熟妇av一区二区三区| 国产精品一区二区免费欧美| 精品久久久久久久久久免费视频| 黑人欧美特级aaaaaa片| 夜夜躁狠狠躁天天躁| 精品国内亚洲2022精品成人| 脱女人内裤的视频| 久久青草综合色| 窝窝影院91人妻| 波多野结衣巨乳人妻| 亚洲国产日韩欧美精品在线观看 | 国产精品美女特级片免费视频播放器 | 精品国产乱码久久久久久男人| 淫秽高清视频在线观看| 俄罗斯特黄特色一大片| 黄频高清免费视频| 男女视频在线观看网站免费 | 日本精品一区二区三区蜜桃| 亚洲精品美女久久av网站| 亚洲中文av在线| 露出奶头的视频| 国产伦人伦偷精品视频| 亚洲aⅴ乱码一区二区在线播放 | 久久久久久久久久黄片| 十八禁网站免费在线| videosex国产| av在线天堂中文字幕| 欧美zozozo另类| 久久久久九九精品影院| 欧美日韩一级在线毛片| 国产亚洲欧美98| 午夜福利成人在线免费观看| 久久午夜综合久久蜜桃| 一级毛片女人18水好多| 老司机深夜福利视频在线观看| 免费电影在线观看免费观看| 别揉我奶头~嗯~啊~动态视频| 我的亚洲天堂| 国产一区二区三区视频了| 免费在线观看完整版高清| 国内精品久久久久久久电影| 国产1区2区3区精品| 久久精品人妻少妇| 91麻豆精品激情在线观看国产| av有码第一页| 亚洲一码二码三码区别大吗| 亚洲成av人片免费观看| 成年免费大片在线观看| 日本在线视频免费播放| 成人欧美大片| 亚洲最大成人中文| 亚洲精品在线美女| 2021天堂中文幕一二区在线观 | 欧美日韩一级在线毛片| 精品午夜福利视频在线观看一区| 级片在线观看| 久久久久久久久中文| 国产精品免费视频内射| 亚洲男人天堂网一区| 欧美人与性动交α欧美精品济南到| 精品国内亚洲2022精品成人| 黄色 视频免费看| 在线天堂中文资源库| 欧美日韩福利视频一区二区| 免费一级毛片在线播放高清视频| 日本免费一区二区三区高清不卡| 久久中文字幕一级| 97超级碰碰碰精品色视频在线观看| 级片在线观看| 中文字幕精品免费在线观看视频| 一区二区三区激情视频| 亚洲自拍偷在线| 亚洲av成人av| 亚洲国产精品成人综合色| 欧美日韩福利视频一区二区| 亚洲男人的天堂狠狠| 丁香六月欧美| 变态另类成人亚洲欧美熟女| av天堂在线播放| 1024香蕉在线观看| 中文字幕另类日韩欧美亚洲嫩草| 黄片大片在线免费观看| 精品久久蜜臀av无| 日日爽夜夜爽网站| 亚洲专区国产一区二区| 神马国产精品三级电影在线观看 | 欧美亚洲日本最大视频资源| 丁香欧美五月| 久久久久久国产a免费观看| 最近最新免费中文字幕在线| 热99re8久久精品国产| 久久午夜综合久久蜜桃| 熟妇人妻久久中文字幕3abv| 欧美亚洲日本最大视频资源| 国产又色又爽无遮挡免费看| 成人免费观看视频高清| 国产成人影院久久av| 成人国语在线视频| 亚洲精品av麻豆狂野| 久热爱精品视频在线9| 91成人精品电影| 日韩三级视频一区二区三区| 中国美女看黄片| 老司机午夜福利在线观看视频| 波多野结衣巨乳人妻| 欧美三级亚洲精品| 久久亚洲精品不卡| 俄罗斯特黄特色一大片| 99在线视频只有这里精品首页| 欧洲精品卡2卡3卡4卡5卡区| 女人爽到高潮嗷嗷叫在线视频| 丝袜在线中文字幕| 91麻豆精品激情在线观看国产| 亚洲精品中文字幕一二三四区| 黄色片一级片一级黄色片| 人人澡人人妻人| 一级a爱视频在线免费观看| 久久国产乱子伦精品免费另类| 国产精品九九99| 国产色视频综合| 黄色丝袜av网址大全| 久久久久国产精品人妻aⅴ院| 日韩精品青青久久久久久| 一二三四社区在线视频社区8| 久久久久国内视频| 亚洲国产中文字幕在线视频| 又黄又粗又硬又大视频| 婷婷精品国产亚洲av在线| av有码第一页| 麻豆国产av国片精品| 两个人看的免费小视频| 在线观看66精品国产| 最近最新免费中文字幕在线| 久久精品国产99精品国产亚洲性色| 精品高清国产在线一区| ponron亚洲| 亚洲精品美女久久av网站| 91国产中文字幕| 中文字幕最新亚洲高清| 亚洲五月色婷婷综合| 国产蜜桃级精品一区二区三区| 亚洲av熟女| 成人三级做爰电影| www日本在线高清视频| 亚洲国产精品sss在线观看| 欧美日韩亚洲综合一区二区三区_| 国产成人一区二区三区免费视频网站| 成人一区二区视频在线观看| 老司机靠b影院| av欧美777| 国产成人精品无人区| 欧美成人免费av一区二区三区| 免费在线观看视频国产中文字幕亚洲| 国产三级在线视频| 免费看美女性在线毛片视频| 首页视频小说图片口味搜索| 91九色精品人成在线观看| 国产精品久久久久久亚洲av鲁大| 琪琪午夜伦伦电影理论片6080| 中文字幕人妻熟女乱码| 脱女人内裤的视频| 日韩免费av在线播放| videosex国产| 欧美色欧美亚洲另类二区| 91在线观看av| 色在线成人网| 日日干狠狠操夜夜爽| 欧美中文综合在线视频| 一级a爱视频在线免费观看| 香蕉av资源在线| 91国产中文字幕| 国产99久久九九免费精品| 每晚都被弄得嗷嗷叫到高潮| 国产午夜精品久久久久久| 中文字幕精品亚洲无线码一区 | 老司机在亚洲福利影院| 欧美黑人精品巨大| 黄片大片在线免费观看| 99精品在免费线老司机午夜| 97碰自拍视频| 黄频高清免费视频| 欧洲精品卡2卡3卡4卡5卡区| 亚洲熟女毛片儿| 国产黄a三级三级三级人| 国产国语露脸激情在线看| 一本综合久久免费| 久久亚洲精品不卡| 国产亚洲av嫩草精品影院| 日日爽夜夜爽网站| 亚洲美女黄片视频| 老鸭窝网址在线观看| 国产伦在线观看视频一区| 欧美日韩中文字幕国产精品一区二区三区| 久久久久精品国产欧美久久久| 亚洲成人精品中文字幕电影| 久久香蕉激情| 精品不卡国产一区二区三区| 男女之事视频高清在线观看| 国产精品野战在线观看| 91麻豆精品激情在线观看国产| 成人国语在线视频| 亚洲五月色婷婷综合| 琪琪午夜伦伦电影理论片6080| 国产v大片淫在线免费观看| 亚洲熟女毛片儿| 亚洲av电影在线进入| 亚洲avbb在线观看| 在线观看www视频免费| 天天添夜夜摸| 成人18禁在线播放| 日韩欧美 国产精品| 18禁美女被吸乳视频| 日本一本二区三区精品| 99国产极品粉嫩在线观看| 日韩av在线大香蕉| 日日夜夜操网爽| 精品国产一区二区三区四区第35| 久久九九热精品免费| 91九色精品人成在线观看| 日韩欧美一区二区三区在线观看| 免费在线观看黄色视频的| 日韩中文字幕欧美一区二区| 色播亚洲综合网| 日本熟妇午夜| 久久中文字幕一级| 看黄色毛片网站| 久久香蕉国产精品| 人妻丰满熟妇av一区二区三区| 欧美午夜高清在线| 国产亚洲欧美精品永久| 无遮挡黄片免费观看| 很黄的视频免费| 好男人电影高清在线观看| 在线观看日韩欧美| 18禁观看日本| 在线观看免费午夜福利视频| 白带黄色成豆腐渣| 嫩草影院精品99| 两个人看的免费小视频| 亚洲国产毛片av蜜桃av| 国产欧美日韩一区二区精品| 欧美激情久久久久久爽电影| 99在线人妻在线中文字幕| www.精华液| 人人妻,人人澡人人爽秒播| 午夜免费激情av| 色婷婷久久久亚洲欧美| 亚洲中文av在线| 国产高清videossex| 精品久久久久久,| 亚洲一卡2卡3卡4卡5卡精品中文| 曰老女人黄片| 久久久精品欧美日韩精品| 两个人看的免费小视频| www日本黄色视频网| 国产麻豆成人av免费视频| 中文亚洲av片在线观看爽| 狠狠狠狠99中文字幕| 国产成人精品久久二区二区免费| 1024香蕉在线观看| 特大巨黑吊av在线直播 | 老鸭窝网址在线观看| 一本综合久久免费| 国产精品九九99| 久久国产精品人妻蜜桃| 啦啦啦免费观看视频1| 99热6这里只有精品| 999久久久精品免费观看国产| 欧美zozozo另类| 搡老熟女国产l中国老女人| 亚洲熟妇熟女久久| 最近最新中文字幕大全电影3 | 99精品欧美一区二区三区四区| 99热这里只有精品一区 | 午夜免费成人在线视频| 欧美色视频一区免费| 精品不卡国产一区二区三区| 国产成人精品无人区| 亚洲一区二区三区不卡视频| 国产亚洲欧美在线一区二区| 男女那种视频在线观看| 黄色女人牲交| 亚洲精品国产一区二区精华液| 色在线成人网| 亚洲欧美激情综合另类| 18禁观看日本| 99精品久久久久人妻精品| 久久久久久免费高清国产稀缺| 久久精品91无色码中文字幕| 国产精品国产高清国产av| 51午夜福利影视在线观看|