• <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.

    韩国av一区二区三区四区| 国产男靠女视频免费网站| 乱码一卡2卡4卡精品| 国产激情偷乱视频一区二区| 中国美白少妇内射xxxbb| 国产精品98久久久久久宅男小说| 身体一侧抽搐| 最好的美女福利视频网| 免费不卡的大黄色大毛片视频在线观看 | 听说在线观看完整版免费高清| a在线观看视频网站| 欧美3d第一页| 中文字幕久久专区| 欧美绝顶高潮抽搐喷水| 国产一区二区三区视频了| 国产真实伦视频高清在线观看 | 国产精品久久久久久久久免| 国产aⅴ精品一区二区三区波| 天堂网av新在线| 欧美一级a爱片免费观看看| 男插女下体视频免费在线播放| 欧美成人一区二区免费高清观看| 亚洲欧美精品综合久久99| 亚洲无线观看免费| 97人妻精品一区二区三区麻豆| 亚洲av日韩精品久久久久久密| 国产精品嫩草影院av在线观看 | 美女大奶头视频| 18禁在线播放成人免费| 国产精品伦人一区二区| 欧美日韩瑟瑟在线播放| 黄色女人牲交| 最近最新免费中文字幕在线| 观看免费一级毛片| 亚洲国产高清在线一区二区三| 国产免费av片在线观看野外av| av视频在线观看入口| 成人二区视频| 亚洲精品一卡2卡三卡4卡5卡| 男人舔女人下体高潮全视频| 成人无遮挡网站| 国产精品永久免费网站| 日韩中文字幕欧美一区二区| 免费电影在线观看免费观看| 久久久久九九精品影院| 免费观看精品视频网站| 欧美一区二区精品小视频在线| 精品午夜福利在线看| 禁无遮挡网站| 日韩一区二区视频免费看| 国产免费一级a男人的天堂| 国产精品人妻久久久影院| 麻豆国产97在线/欧美| 色综合亚洲欧美另类图片| 欧洲精品卡2卡3卡4卡5卡区| 国产精华一区二区三区| 国产中年淑女户外野战色| 欧美色欧美亚洲另类二区| 老司机福利观看| 色5月婷婷丁香| ponron亚洲| 深夜精品福利| 三级毛片av免费| 小蜜桃在线观看免费完整版高清| 亚洲欧美日韩无卡精品| 欧美黑人巨大hd| 亚洲美女搞黄在线观看 | .国产精品久久| 国产蜜桃级精品一区二区三区| 一个人观看的视频www高清免费观看| 精品日产1卡2卡| 国产精品国产高清国产av| 日本 av在线| 夜夜爽天天搞| 99热这里只有是精品50| 三级国产精品欧美在线观看| 色综合站精品国产| 免费看av在线观看网站| 不卡视频在线观看欧美| 热99在线观看视频| 性插视频无遮挡在线免费观看| 亚洲美女搞黄在线观看 | 日日撸夜夜添| 少妇猛男粗大的猛烈进出视频 | 午夜精品久久久久久毛片777| 国产色婷婷99| 亚洲真实伦在线观看| av在线天堂中文字幕| 成年免费大片在线观看| 国语自产精品视频在线第100页| 热99在线观看视频| 国产成人影院久久av| 久久国产精品人妻蜜桃| 大型黄色视频在线免费观看| 欧美一区二区亚洲| 日本在线视频免费播放| 我的老师免费观看完整版| 有码 亚洲区| 日韩 亚洲 欧美在线| 国产亚洲精品综合一区在线观看| 亚洲美女黄片视频| 两个人视频免费观看高清| 午夜福利成人在线免费观看| 全区人妻精品视频| 天堂√8在线中文| 国产成人aa在线观看| 午夜免费男女啪啪视频观看 | 少妇裸体淫交视频免费看高清| 国产单亲对白刺激| 欧美日韩综合久久久久久 | 亚洲美女搞黄在线观看 | 精品午夜福利在线看| 国产精品一区www在线观看 | 午夜福利视频1000在线观看| 国产av麻豆久久久久久久| 欧美又色又爽又黄视频| 两个人视频免费观看高清| 别揉我奶头 嗯啊视频| 99视频精品全部免费 在线| 欧美精品国产亚洲| 男人狂女人下面高潮的视频| av天堂中文字幕网| 国国产精品蜜臀av免费| 精品午夜福利视频在线观看一区| 国产极品精品免费视频能看的| 床上黄色一级片| 免费在线观看成人毛片| 老司机深夜福利视频在线观看| 嫩草影视91久久| 3wmmmm亚洲av在线观看| 男女下面进入的视频免费午夜| 18禁在线播放成人免费| 国产亚洲欧美98| 精品久久久久久久久久免费视频| 人妻制服诱惑在线中文字幕| 老熟妇仑乱视频hdxx| 我的女老师完整版在线观看| 日本a在线网址| 给我免费播放毛片高清在线观看| 亚洲人成伊人成综合网2020| 男插女下体视频免费在线播放| 51国产日韩欧美| 内地一区二区视频在线| 国产精品不卡视频一区二区| av在线蜜桃| 欧美又色又爽又黄视频| 男人舔女人下体高潮全视频| 啪啪无遮挡十八禁网站| 波多野结衣巨乳人妻| 深夜a级毛片| 亚洲五月天丁香| 天堂网av新在线| 麻豆久久精品国产亚洲av| 日韩强制内射视频| 国产精品久久电影中文字幕| 999久久久精品免费观看国产| 国产成人福利小说| 亚洲人成网站在线播| 亚洲中文字幕日韩| 午夜亚洲福利在线播放| 国产在视频线在精品| 99久久久亚洲精品蜜臀av| 乱人视频在线观看| 又黄又爽又刺激的免费视频.| 欧美+亚洲+日韩+国产| 91在线精品国自产拍蜜月| 又爽又黄a免费视频| 国产精品综合久久久久久久免费| 丝袜美腿在线中文| 看黄色毛片网站| 精品久久久久久,| 色精品久久人妻99蜜桃| 亚洲无线观看免费| 亚州av有码| 国产日本99.免费观看| 天天一区二区日本电影三级| 97超视频在线观看视频| 免费看日本二区| av福利片在线观看| 我要看日韩黄色一级片| 淫妇啪啪啪对白视频| 午夜影院日韩av| videossex国产| 成人欧美大片| 久久婷婷人人爽人人干人人爱| 1024手机看黄色片| 国内精品久久久久精免费| netflix在线观看网站| 黄色配什么色好看| 国产高潮美女av| 久久久久久久久大av| 深爱激情五月婷婷| 性欧美人与动物交配| 小蜜桃在线观看免费完整版高清| 精品人妻偷拍中文字幕| www.色视频.com| 亚洲aⅴ乱码一区二区在线播放| 99热这里只有是精品在线观看| 国产伦精品一区二区三区四那| 成年女人看的毛片在线观看| 香蕉av资源在线| 亚洲五月天丁香| 亚洲久久久久久中文字幕| 亚洲最大成人中文| 日本与韩国留学比较| 18禁在线播放成人免费| x7x7x7水蜜桃| 看黄色毛片网站| 22中文网久久字幕| 国产精品三级大全| 精品人妻熟女av久视频| 赤兔流量卡办理| 91久久精品电影网| 欧美黑人巨大hd| 小蜜桃在线观看免费完整版高清| 亚洲av免费在线观看| 亚洲成a人片在线一区二区| 麻豆成人av在线观看| 噜噜噜噜噜久久久久久91| 黄色丝袜av网址大全| 嫩草影视91久久| 欧美成人一区二区免费高清观看| 日日摸夜夜添夜夜添av毛片 | 日本三级黄在线观看| 嫩草影视91久久| 最近最新免费中文字幕在线| 男女那种视频在线观看| 联通29元200g的流量卡| 国产高清有码在线观看视频| 亚洲黑人精品在线| 成年免费大片在线观看| 亚洲av成人精品一区久久| 少妇熟女aⅴ在线视频| 麻豆成人av在线观看| 久久中文看片网| 久久久久国内视频| 精品久久久久久成人av| 在线观看舔阴道视频| 亚洲专区中文字幕在线| 色噜噜av男人的天堂激情| 观看美女的网站| 亚洲最大成人手机在线| 女的被弄到高潮叫床怎么办 | 日日摸夜夜添夜夜添小说| 成年版毛片免费区| av在线观看视频网站免费| 亚州av有码| 在线观看午夜福利视频| 最近在线观看免费完整版| 老司机午夜福利在线观看视频| 色尼玛亚洲综合影院| 99久久九九国产精品国产免费| av国产免费在线观看| 精品久久久久久久久av| 97热精品久久久久久| 欧美国产日韩亚洲一区| 亚洲,欧美,日韩| 国产一区二区三区av在线 | 中文字幕久久专区| 久久久久国产精品人妻aⅴ院| 国产高潮美女av| 婷婷六月久久综合丁香| 黄色视频,在线免费观看| 免费观看的影片在线观看| 国产中年淑女户外野战色| 久久国产精品人妻蜜桃| 亚洲性久久影院| 久久精品91蜜桃| 嫩草影院入口| 久久亚洲真实| 免费在线观看成人毛片| 国产色爽女视频免费观看| 18+在线观看网站| 国内精品美女久久久久久| 久久久久久久久大av| 黄色视频,在线免费观看| 日本一二三区视频观看| 久久久久性生活片| 成人美女网站在线观看视频| 看免费成人av毛片| 91av网一区二区| 欧美在线一区亚洲| 精品久久久久久成人av| 亚洲最大成人中文| 国内揄拍国产精品人妻在线| 69人妻影院| 99久久精品国产国产毛片| 国产 一区 欧美 日韩| .国产精品久久| 国产主播在线观看一区二区| 一个人观看的视频www高清免费观看| 日本撒尿小便嘘嘘汇集6| 小蜜桃在线观看免费完整版高清| 国产视频一区二区在线看| 亚洲自拍偷在线| 国产精品三级大全| 99在线视频只有这里精品首页| 欧美国产日韩亚洲一区| 看十八女毛片水多多多| 精品人妻视频免费看| 一级毛片久久久久久久久女| 国内精品久久久久久久电影| 中文资源天堂在线| 91午夜精品亚洲一区二区三区 | 免费在线观看成人毛片| 又黄又爽又刺激的免费视频.| 一个人看视频在线观看www免费| 91午夜精品亚洲一区二区三区 | 亚洲av电影不卡..在线观看| 免费人成在线观看视频色| 动漫黄色视频在线观看| 色综合亚洲欧美另类图片| 欧美成人性av电影在线观看| 免费看av在线观看网站| 国产亚洲精品av在线| 日本黄大片高清| avwww免费| 深爱激情五月婷婷| av女优亚洲男人天堂| 国产v大片淫在线免费观看| 国产精品永久免费网站| 亚洲成av人片在线播放无| 免费在线观看影片大全网站| 国产黄a三级三级三级人| 桃红色精品国产亚洲av| 国产av一区在线观看免费| 色综合色国产| 国产一区二区在线观看日韩| 91在线精品国自产拍蜜月| 直男gayav资源| 22中文网久久字幕| 亚洲av一区综合| 国产综合懂色| 男人舔奶头视频| 欧美一区二区亚洲| 国产中年淑女户外野战色| 中文字幕av在线有码专区| 97碰自拍视频| 国产精品日韩av在线免费观看| 好男人在线观看高清免费视频| 国产 一区 欧美 日韩| 国产精品1区2区在线观看.| 91av网一区二区| 国产精品精品国产色婷婷| 欧美xxxx性猛交bbbb| 免费人成在线观看视频色| 免费观看在线日韩| 久久精品国产清高在天天线| 欧美日韩乱码在线| 麻豆av噜噜一区二区三区| av在线亚洲专区| av中文乱码字幕在线| 成年女人永久免费观看视频| 男人的好看免费观看在线视频| 久久精品国产鲁丝片午夜精品 | 少妇人妻一区二区三区视频| 欧美性猛交黑人性爽| 简卡轻食公司| av专区在线播放| 亚洲第一区二区三区不卡| 99久久精品热视频| 性插视频无遮挡在线免费观看| 中出人妻视频一区二区| 国产午夜精品论理片| 村上凉子中文字幕在线| 美女xxoo啪啪120秒动态图| 夜夜夜夜夜久久久久| 亚洲精华国产精华液的使用体验 | 性色avwww在线观看| 国产在视频线在精品| 夜夜看夜夜爽夜夜摸| 美女黄网站色视频| 午夜精品久久久久久毛片777| 精品久久久久久久人妻蜜臀av| 又粗又爽又猛毛片免费看| 国产精品自产拍在线观看55亚洲| 女人被狂操c到高潮| 两性午夜刺激爽爽歪歪视频在线观看| 午夜a级毛片| 久久久久久久精品吃奶| 中文字幕免费在线视频6| 国产精品1区2区在线观看.| 男女做爰动态图高潮gif福利片| 精品欧美国产一区二区三| 可以在线观看毛片的网站| 久久欧美精品欧美久久欧美| 我要看日韩黄色一级片| 久久精品国产99精品国产亚洲性色| 女同久久另类99精品国产91| 国产精品久久久久久精品电影| 久久国内精品自在自线图片| 亚洲三级黄色毛片| 亚洲专区中文字幕在线| 国产高清视频在线观看网站| 色视频www国产| 亚洲色图av天堂| 亚洲av免费在线观看| 精品人妻熟女av久视频| av在线亚洲专区| 男女做爰动态图高潮gif福利片| 春色校园在线视频观看| 99久久精品一区二区三区| 午夜免费成人在线视频| 岛国在线免费视频观看| 国产精品久久久久久久久免| 国产黄片美女视频| 波多野结衣高清作品| 51国产日韩欧美| 伦精品一区二区三区| 久久婷婷人人爽人人干人人爱| 最后的刺客免费高清国语| 又黄又爽又刺激的免费视频.| 大型黄色视频在线免费观看| 男女视频在线观看网站免费| 最近最新中文字幕大全电影3| 国产女主播在线喷水免费视频网站 | 午夜亚洲福利在线播放| 欧美潮喷喷水| 亚洲国产精品sss在线观看| 嫩草影院入口| 免费在线观看日本一区| 波多野结衣高清无吗| 俺也久久电影网| 欧美日本视频| 美女免费视频网站| 国产亚洲av嫩草精品影院| a级一级毛片免费在线观看| 婷婷色综合大香蕉| 久久久久九九精品影院| 欧美高清成人免费视频www| 日本色播在线视频| 美女cb高潮喷水在线观看| 看黄色毛片网站| 国产主播在线观看一区二区| 精品人妻熟女av久视频| 国产探花极品一区二区| АⅤ资源中文在线天堂| 久久99热6这里只有精品| 又紧又爽又黄一区二区| 精品久久国产蜜桃| 国产v大片淫在线免费观看| 精品不卡国产一区二区三区| 国产精品自产拍在线观看55亚洲| 久久精品国产亚洲网站| 91在线观看av| 久久99热这里只有精品18| 亚洲成人久久爱视频| 国产精品99久久久久久久久| 一进一出好大好爽视频| 如何舔出高潮| 亚洲一区二区三区色噜噜| 免费看美女性在线毛片视频| 欧美一级a爱片免费观看看| 久久这里只有精品中国| 99在线视频只有这里精品首页| 一级黄色大片毛片| 偷拍熟女少妇极品色| 欧美一级a爱片免费观看看| 成人国产麻豆网| 日本黄色片子视频| 久久精品91蜜桃| 欧美日韩乱码在线| 人妻丰满熟妇av一区二区三区| 国产真实伦视频高清在线观看 | 又紧又爽又黄一区二区| 精品免费久久久久久久清纯| 大又大粗又爽又黄少妇毛片口| 欧美bdsm另类| 国产成人av教育| 麻豆成人av在线观看| 亚洲无线在线观看| 天堂影院成人在线观看| 日日啪夜夜撸| www.www免费av| 搡老熟女国产l中国老女人| 久久久成人免费电影| 欧美丝袜亚洲另类 | 日本爱情动作片www.在线观看 | 精品人妻一区二区三区麻豆 | 久久精品国产亚洲av涩爱 | 99热精品在线国产| 国产伦精品一区二区三区视频9| 国产真实乱freesex| 成年版毛片免费区| 国产精品女同一区二区软件 | 久久久久久伊人网av| 成人二区视频| 久久欧美精品欧美久久欧美| 午夜福利在线观看吧| 在线播放无遮挡| 国产高清视频在线观看网站| 国产v大片淫在线免费观看| 99久久久亚洲精品蜜臀av| 男人的好看免费观看在线视频| 欧美最新免费一区二区三区| 最近视频中文字幕2019在线8| 国内精品久久久久精免费| www.色视频.com| 男人狂女人下面高潮的视频| 又粗又爽又猛毛片免费看| 男人的好看免费观看在线视频| 国产免费男女视频| 99九九线精品视频在线观看视频| 真实男女啪啪啪动态图| www日本黄色视频网| 亚洲人成伊人成综合网2020| 日韩中文字幕欧美一区二区| 日本a在线网址| 中文字幕高清在线视频| 免费不卡的大黄色大毛片视频在线观看 | 看片在线看免费视频| 国内揄拍国产精品人妻在线| 如何舔出高潮| 2021天堂中文幕一二区在线观| av女优亚洲男人天堂| 国产高清不卡午夜福利| 色播亚洲综合网| 欧美激情国产日韩精品一区| 欧美一区二区国产精品久久精品| 欧美区成人在线视频| 国产精品女同一区二区软件 | 一a级毛片在线观看| 很黄的视频免费| 国产爱豆传媒在线观看| 给我免费播放毛片高清在线观看| .国产精品久久| 日本爱情动作片www.在线观看 | 99久久精品一区二区三区| 88av欧美| 国产精品一区二区免费欧美| 日日撸夜夜添| 在线播放国产精品三级| 精品国产三级普通话版| 欧美另类亚洲清纯唯美| av国产免费在线观看| 97超级碰碰碰精品色视频在线观看| 国产午夜福利久久久久久| 亚洲五月天丁香| 欧美+日韩+精品| 我要看日韩黄色一级片| 成人欧美大片| 欧美精品国产亚洲| 中文字幕人妻熟人妻熟丝袜美| 五月伊人婷婷丁香| 久久99热6这里只有精品| 日韩人妻高清精品专区| 午夜免费激情av| 国产视频一区二区在线看| 69av精品久久久久久| 国产精品99久久久久久久久| 九九热线精品视视频播放| 中文字幕高清在线视频| 亚洲av日韩精品久久久久久密| 日韩在线高清观看一区二区三区 | 亚洲精品久久国产高清桃花| h日本视频在线播放| 亚洲欧美日韩卡通动漫| 伦理电影大哥的女人| 男女之事视频高清在线观看| 久久久久久久久久成人| 18禁在线播放成人免费| 国产午夜精品论理片| 成人特级av手机在线观看| 免费在线观看日本一区| 国产精品日韩av在线免费观看| 啦啦啦观看免费观看视频高清| 亚洲一级一片aⅴ在线观看| 国产av一区在线观看免费| 男人舔女人下体高潮全视频| 成熟少妇高潮喷水视频| 久久国产乱子免费精品| 亚洲精品456在线播放app | 亚洲人与动物交配视频| 日韩人妻高清精品专区| 在线免费十八禁| 日本-黄色视频高清免费观看| 男人的好看免费观看在线视频| 亚洲欧美日韩东京热| 日本a在线网址| 成人国产综合亚洲| 成人午夜高清在线视频| 成人高潮视频无遮挡免费网站| 一夜夜www| av福利片在线观看| 久久久国产成人免费| 成人国产一区最新在线观看| 久久精品国产亚洲网站| 日本与韩国留学比较| 女人被狂操c到高潮| 亚洲成人中文字幕在线播放| 色av中文字幕| 一a级毛片在线观看| 最近最新中文字幕大全电影3| 亚洲av美国av| a级一级毛片免费在线观看| or卡值多少钱| 色5月婷婷丁香| 亚洲欧美日韩卡通动漫| 在线观看免费视频日本深夜| 国产亚洲精品久久久com| 成人av一区二区三区在线看| 乱系列少妇在线播放| 国内少妇人妻偷人精品xxx网站| 亚洲成人久久爱视频| 波多野结衣高清无吗| 国产v大片淫在线免费观看| 草草在线视频免费看| 最近最新中文字幕大全电影3| 精品午夜福利在线看| 美女被艹到高潮喷水动态| 18禁在线播放成人免费| 黄色一级大片看看| 最近最新免费中文字幕在线| 久久香蕉精品热| 亚洲国产色片|