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

    Variance-based Sensitivity Analyses of Piezoelectric Models

    2015-12-13 06:10:14LahmerIlgandLerch

    T.Lahmer,J.Ilgand R.Lerch

    Variance-based Sensitivity Analyses of Piezoelectric Models

    T.Lahmer1,J.Ilg2and R.Lerch2

    In the recent years many publications appeared putting emphasis on the simulation-based identification of piezoelectric material parameters from electrical or mechanical measurements and combinations of them.By experience,one is aware of the importance of a single input parameter.However,it is not yet fully understood and in particular quantified to which extend missing knowledge in the single parameters(parameter uncertainty)influences the quality of the model’s prognosis.In this paper,we adapt and apply variance-based sensitivity measures to models describing the piezoelectric effect in the linear case and derive global information about the single input parameter’s sensitivities.

    Piezoelectricity,sensitivity analysis,uncertainty analysis,material parameter determination,Finite Element Method.

    1 Introduction

    Piezoelectric ceramics are nowadays widely used in many sensor and actuator applications,like microphones,sonar technologies or injection-valves.The increasing number of the application of numerical simulations in design,development and control of piezoelectric devices requires reliable simulation outcomes,see,e.g.,Zhang,Xu,and Wang(2014);Dziatkiewicz and Fedelinski(2007)for different formulations or Olyaie,Razfar,Wang,and Kansa(2011);Olyaie,Razfar,and Kansa(2011)for topics based on the design optimization of piezoelectric structures.

    Simulation responses can be assumed to be reliable,when the underlying mathematical model itself can be validated and the numerical treatment verified,but also when all input parameters are known precisely.This,however,is often not the case,as parameter identification methods are costly,time consuming and underlay other challenges and restrictions.Further,inherent nonlinear behavior of the material parameters needs to be taken into account,e.g.temperature or ex-citation frequency dependency.These non-linear parameter relationships are not always fully understood and therefore often neglected in the analysis.Many publications can be found which deal with the characterization of piezoelectric materials,e.g.Du,Wang,and Uchino(2004);IEEEStandard(1985);Kaltenbacher,Lahmer,and Mohr(2006);Kwok,Chan,and Choy(1997);Sherrit,Wiederick,and Mukherjee(1997,1992);Jonsson,Andersson,and Lindahl(2013),where different strategies are applied,e.g.nonlinear optimization,genetic programming or regularizing methods,see examples in Poteralskia,Szczepanika,Dziatkiewicza,W.,and Burczy′nskiab(2011);Lahmer(2009);Rupitsch,Sutor,Ilg,and Lerch(2010).Also contributions differ in the types of measurements evaluated,e.g.impedance,charges or displacement measurements and combinations of them,see,e.g.,IEEEStandard(1985);Lahmer(2009);Rupitsch and Lerch(2009);Rupitsch,Sutor,Ilg,and Lerch(2010).All these approaches might have theirpros and consbut all have in common that the uncertainty in the parameter estimates cannot be reduced completely to zero.The sources of these uncertainties are twofold:Firstly,there is the aleatoric uncertainty caused by the randomness in the material composition or by Gaussian noise in the measurements used to estimate the parameter.Particularly,in the case of sparse measurements,the statistical uncertainty rises.Secondly,pure epistemic uncertainty is present due to systematic errors in either the model setup or the con figuration of the experiment.In this case,systematic errors are directly mapped to the parameter estimates.Regularization methods may reduce this effect a little,but a certain amount of bias remains.As a consequence,one has to consider the parameter estimates as random variables each following a specific probability distribution which is characterized by its statistical properties like mean value and standard deviation(parameter uncertainty).In particular the latter might be high due to low sensitivities of the parameters,see,e.g.,Lahmer,Kaltenbacher,Kaltenbacher,Leder,and Lerch(2008).An uncertainty quantification of MEMS(micro-electromechanical systems)derived in a unified framework is reported in Alwan and Aluru(2011).Generally,in order to reduce the parameter uncertainty,further experiments need to be conducted,probably combined with a preceding design of experiments,see Rus,Palma,and Pérez-Aparicio(2011);Lahmer,Kaltenbacher,and Schulz(2008).

    In order to answer the question if this additional experimental and calibration efforts are justified,a sensitivity analysis should be conducted.Doing so,it can be quantified to which extend uncertainties in the input parameters in fluence the uncertainty in the model’s response or by the de finition of Saltelli,Ratto,Andres,Campolongo,Cariboni,Gatelli,M.,and Tarantola(2008)“how uncertainty in the output of a model can be apportioned to different sources of uncertainty in the model input”.Sensitivity analyses for piezoelectric models based on partial derivatives are reported in Lahmer(2009);Lahmer,Kaltenbacher,Kaltenbacher,Leder,and Lerch(2008),a visual type of sensitivity analysis is reported in Rupitsch,Wolf,Sutor,and Lerch(2012).

    The aim of this article is to quantify the sources of uncertainty using global sensitivity measures for different typical piezoelectric problems,e.g.to clarify how a coefficient of variation of 10%of any input parameter in fluences the frequency of the first eigenmode.

    From this information,the following questions might be answered:

    ·Which parameters are mainly responsible for the system’s behavior?

    ·Which parameter has a completely non-in fluential variance and can therefore be considered deterministically?

    ·Is my simulation reliable or does the model answer depend strongly on weak assumptions on the input parameters?

    ·To which type of measurements are my parameters sensitive?(Design of Experiments)

    ·Prioritization of further research:Which input factors require further investigations and which not?

    Finally,results of a carefully conducted sensitivity analysis supports the user’s understanding of the model and the process of decision-making based on numerical simulation results,as e.g.done in Most(2011).In particular,this point should be discussed intensively as it(de-)motivates further time and money consuming research activities.

    The sensitivity measures applied in this research are variance-based methods,i.e.the variance is considered as the measure of uncertainty.Further,the measures are of a global character and can handle both non-linearities and interactions between pairs of input parameters.In the case that models are prohibitively expensive for global measures,sensitivity indices based on their correlations or computed by response surface methods are shown to be an alternative.

    The paper is organized as follows:In Section 2 the model,i.e.balance equations and piezoelectric constitutive equations are introduced and the Finite Element discretisation is briefly discussed.In Section 3 variance-based sensitivity measures are introduced and efficient computing strategies are presented.Section 4 contains results of this analysis and Section 5 closes this article with a discussion of the results obtained.

    2 Piezoelectric Effect And Equations

    The material law describing the piezoelectric effect under small signal loading is given by

    It relates the mechanical stress tensor[σ]and the electrical displacementD,respectively,to the mechanical strain tensor[S]and the electric fieldE.Due to the symmetry of the mechanical tensors[σ]and[S],we may rewrite them in Voigt notation as follows

    In the electrostatic case,the relation between electric fieldEand the electric potentialφis given by

    Mechanical strainsSare related to the mechanical displacementsuby

    whereBis the first order differential operator relating strains and displacements.The material tensors[cE],[εS],and[e],appearing in(1),(2)are the elasticity coefficients,the dielectric constants,and the piezoelectric coupling coefficients,respectively.According to the crystal structure and polarization of the piezoelectric material,these matrices show a certain symmetry and sparsity pattern(cf.IEEEStandard(1985)).For the 6mm crystal class we have,e.g.

    The mechanical behavior of piezoelectric materials is described by Newton’s law

    where ρ denotes the mass density.Since piezoelectric materials are insulating,i.e.,do not contain free volume charges

    Plugging the constitutive laws(1),(2)in the balance equations(5)and(6),one obtains a fully determined set of partial differential equations which are solved for the primary unknowns u and φ inside a piezoelectric domain ?

    The following boundary conditions are according to the typical experimental setting of vanishing normal stresses at the boundaries,and two electrodes being applied at opposite positions Γgand Γe,one of them loaded by a prescribed electric potential φe

    The application of the Fourier Transformation and Finite Element discretization to equations in(7)gives the following set of coupled algebraic equations in frequency domain,see e.g.Lerch(1990);Kaltenbacher(2004)

    Herein Kuu,and Muudenote the mechanical stiffness,and mass matrix,respectively,Kφφand Kuφare the dielectric stiffness-and the piezoelectric coupling matrix.The variable ω denotes the angular frequency.The Finite Element solutions are the nodal vector of displacementsand the nodal vector of the electric potential

    In order to take energy-dissipation of the material into account,complex valued material parameters are considered,which is equivalent to a velocity proportional damping in time domain.For each tensor,a constant scaling factor between the real and imaginary part is introduced assuming that an equal damping of the different entries in the material tensors is a valid approach,i.e.

    In(10),jdenotes the imaginary unit,i.e.j2=?1 in the complex number system.The eigenvalues problem in(9)is solved by the aid of the ARPACK package,see Sorensen,Lehoucq,Yang,and Maschhoff(1996),yielding the angular eigenfrequenciesω1,ω2,....The eigenfrequencies correspond to all possible mechanical eigenmodes,even the ones that cannot be excited electrically.

    3 Global Sensitivity Analysis

    Sensitivity analyses can be conducted in different ways.The most simple approach and mathematically motivated by the definition of a sensitivity in calculus is to compute partial derivatives.This strategy is justified for a first understanding of the problem,however,one needs to be aware that it is just alocalmeasure as it computes a parameter’s sensitivity when the values of all other parameters are fixed,i.e.assumed to be known exactly.The parameter’s sensitivity however might change drastically when the other parameters in the model contain different values.This can be exemplified regarding the sensitivity of the mechanical deflection with respect to an infinitesimal change of the mechanical damping parameter.If we are away of resonance,no significant influence will be seen.However,a change of stiffness properties may shift the resonance’s positions and the damping parameter might become visibly influential as it is now mainly responsible for the magnitude of the deflection in resonance.

    Due to these interrelations,we considerglobalmeasures,where all parameters are varied while determining the sensitivity of a parameter of interest.Those methods where firstly introduced in 1973 by Cukier et.al.Cukier,Fortuin,Shuler,Petschek,and Schaibly(59)and were extended by Homma,Imam,Sobol and Saltelli,Homma and Saltelli(1996);Sobol(1990);Saltelli,Tarantola,and Chan(1999)and applied in various models in economy,engineering and science,see e.g.Sin,Gernaey,Neumann,Loosdrecht,and Gujer(2009);Keitel,Karaki,Lahmer,Nikulla,and Zabel(2011);Kami′nski and Lauke(2013).According to our knowledge,global variance-based sensitivity measures have not been applied to any qualitative assessment of models for smart materials,like piezoelectric ones.Let us assume,that in the followingFis the forward operator of our model

    which maps thenmodel input parametersp1,...,pnto the model outputy=F(p1,...,pn).WithXandY,we denote the finite-dimensional parameter and measurement spaces.WithE(y)andV(y),we denote the mean and the variance of the model’s response,respectively,and byE(pi)andV(pi)fori=1,...,nthe mean and variances of the parameters.

    3.1 Variance-based Sensitivity Measures

    In this context,all parameters are assumed to be uncertain,i.e.they are allowed to vary according to an assumed probability distribution.Now the question arises,what happens to the variance of the model’s response,if we were able to know one parameter precisely,i.e.

    The above is the variance ofyunder the condition thatpiis known to have exactly the value.As the true value ofpicannot be assumed to be not known precisely,the mean of the conditional variance is further regarded.This measures the expected amount by which the uncertainty inyis reduced if we were able to determinepiprecisely on average.By this,the dependence of the measure onp?iis removed.Now

    is on average the expected amount of variance which is removed if we were able to knowpiprecisely.According to the decomposition of variance(law of total variance)

    either a large value ofVpi(Ep~i(y|pi))or a small valueEpi(Vp~i(y|pi))implies thatpiis an important(sensitive)parameter.Normalization gives finally the definition of a first order sensitivity index

    The index Siis therefore a measure of the exclusive influence of the input parameter pi.If the sum of all Siis close to one,the model is additive with respect to its variances and no remarkable interactions between the parameters seem to exist.

    However,for the piezoelectric equations it is not fully understood to which extent choices of single parameters may influence each other so that higher order sensitivity measures should be evaluated in addition.The term

    gives now the expected amount of variance which remains if we were able to know all parameters except piprecisely,i.e.piis allowed to vary over its uncertainty range while all other parameters are kept fixed.After normalization and applying(15)the total effects sensitivity measure is de fined as

    see Sobol(1990);Homma and Saltelli(1996).The total effect index STiincludes both,the in fluence of the first order effects and all higher order(interaction)effects.Consequently,

    If both sums are almost equal,then all input parameters seem to be independent.The quantities 1?∑Sior∑STi?1 are measures of the degree of interactions,which is in particular a valuable hint concerning the unique identifiability of the parameters in the associated inverse problem of parameter identification.

    3.2 Computational Aspects

    A brute force Monte-Carlo-based computation of the variance of the mean or vice versa would actually require N×N evaluations of the piezoelectric problem under investigation,where N denotes the number of samples in a Monte Carlo type simulation.Therefore,for the computation of the senstivity indices we follow the strategy poposed in Saltelli,Ratto,Andres,Campolongo,Cariboni,Gatelli,M.,and Tarantola(2008)going back to earlier approaches by Sobol(1990)and Homma and Saltelli(1996).Our implementation is based on a Latin-Hypercube Sampling strategy with which two matrices A and B of dimensions N×n are filled with random entries.The samples are generated according the assumed probability distribution of the single model parameters.The number of samples N is required to be chosen sufficiently high in order to remove any statistical uncertainty from the results.

    Further matricesCi,i=1,...,nare set up then,which are equal toAexcept theith column,which is taken fromB.Now the model output,i.e.the Finite Element solution of the piezoelectric problem is computed for all rows inA,BandCiyielding 2+nvectors of lengthN.The entries are given by solving

    Where as the first order sensitivity indices might rather be used in the context of research prioritization(which parameters need a more careful investigation/identification),the total order sensitivity indices provide information to which extent models can be simplified as it may be shown that some parameters are completely non-in fluential.

    The random input parameters are in all cases assumed to be normally distributed and are generated using a Latin Hypercube sampling strategy,see e.g.Iman and Conover(1980).

    The Finite Element simulations have been conducted with the noncommercial toolCFS++Kaltenbacher(2010)with the calculation timetFfor one model outputy=F(p1,...,pn).Pre-and post-processing is performed with MATLAB?(The MathWorks,Inc.).The effort to compute(19)and(20)is of orderO(N(n+2)).By using theParallel Computing Toolboxof MATLAB?and an appropriate computer with minimumn+2 CPU’s the computation time reduces toN·tF.

    3.3 Alternative Strategies to Assess Sensitivities

    Variance-based sensitivity measures provide global information,however,under sometimes prohibitive high computational costs.We therefore present two alternatives that approximate the information obtained from the variance-based sensitivity analysis.The first one measures simply correlations between input and output signals;the second alternative discusses the use of variance-based methods on surrogate models.

    3.3.1 Correlation-Based Sensitivity Analysis

    To anticipate,the numerical results show a rather linear relation between the system responses and uncertainties in the input parameters.Therefore,working with weighted linear correlation coefficients may give good approximations of the sensitivity indices.However,as the weighted correlation coefficients are based on Pearson’s definition of correlation,they only measure primary effects and no interactions,thus they are of the same quality as the first order effects.Correlation coefficients are de fined as

    withμpi,μybeing the mean values of the model’s input and output quantities andσpi,σyare the standard deviations ofpiandy.A weighted sensitivity index based on linear correlation can now be de fined as

    Besides their cheaper computation correlation coefficients provide the additional advantage that the sign of correlation(positive/negative)can be retrieved.

    3.3.2 Sensitivity Analysis on Surrogate Models

    As high numbers of samples are generally required to approximate the conditional variances sufficiently precise,uncertainty analysis based on response surfaces is proposed by several authors,see e.g.for the analysis of other engineering problems Marrel,Iooss,Veiga,and Ribatet(2012);Karaki(2012).Here,the true model is approximated by a response surface,also called a meta-model,which allows the main functional dependencies of the model output from model input to be represented.For rather smooth model responses,multivariate polynomial regression is generally a good choice.Any polynomial regression for models withninput parameters requires the definition of an appropriate basis,e.g.with terms of linear and quadratic order

    Higher order terms could be included in this basis,however,then the problem of over fitting could become severe.According to the law of parsimony,often less complex approaches are preferable.The total number of entries in the basis is denoted bym.

    The regression model is given now by

    for which the regression parametersneed to be estimated.This is generally done by minimizing all residuals,i.e.the difference between the regression model and the true model’s response for every sample pointj=1,...,Nse.One method to do so is the ordinary least-squares minimization(sum of squared errors),i.e.

    Minimization of this terms leads to a set of normal equations from which the coefficients can be retrieved as follows

    where all regression coefficients are stored in the vector

    The matrix X is the so-called Vandermonde matrix whose columns are composed by thembasis entries(see e.g.the basis in(23))which vary for every sample over theNserows,for details see,e.g.,Forrester,Sobester,and Keane(2008),Chapter 2.As long as all regression variables are linearly independent this is a well posed problem.In some contexts,a regularized version of the least squares solution could reduce the risk of non-invertability of the matrix(X┬X).Applying Tikhonov regularization one adds a constraint in the form of a penalty term to the minimization problem(24)so that the regression coefficientswill not exceed a given value.Equivalently,one may solve an unconstrained minimization problem of the sum of least-squares approach with a penalty of the formαreg‖β‖2added to the diagonal entries of(X┬X),whereαreg>0 is a regularizing constant.In a Bayesian setting,this strategy is equal to de fining a zero-mean normally distributed prior on the parameter vector.

    The so called coefficient of determinationR2

    measures how well the regression model fits to the true model since it relates the explained variance(variance of the model’s predictions)

    Table 1:Assumed mean values for the parameters of a Pz27 ceramicin N/m2;in As/Vm;exyin As/m2;density ρ in kg/m3).

    Table 1:Assumed mean values for the parameters of a Pz27 ceramicin N/m2;in As/Vm;exyin As/m2;density ρ in kg/m3).

    Hereyiare the model outputs,is the approximated output andis the mean value of the model output evaluated.Values ofR2close to one(at least>0.8)are required for an acceptable approximation of the true model by the regression model.Sensitivity measures which are computed on such response surfaces need to be reduced by the factorR2as some behavior of the true model might remain unexplained,i.e.

    Certainly,there are many more issues to discuss when applying regression models(reduced accuracy,over- fitting,confidence,local/global regressions,...),which however,are beyond the scope of this paper and we refer to Armstrong(2012);Cumming(2012)for further reading on these issues.

    4 Results

    The uncertainty analysis is conducted for two different piezoelectric problems,namely the computation of the electrical impedance of a piezoelectric disc and the computation of displacements of the tip of a piezoelectric unimorph actuator.

    4.1 Electrical impedance-eigenvalue problem

    Figure 1:Electrical impedance curve of a piecoelectric disc with diameter D=25mm and height h=2mm computed with the parameters as given in Table 4.1.The radial(left box)and thickness(right)resonance of the disc are visualized and can be found at~80kHz and~1MHz,respectively.

    The fastest way to compute resonant frequencies is to solve the eigenvalue problem.By doing so,one obtains a list of eigenfrequencies with the corresponding eigenvectors.This list contains all mechanical eigenfrequencies of the given model,even the ones that cannot be electrically excited with the potential Φ de fined by the electrodes and the applied voltageU(see Fig.1).The major problem for this approach is the identification of the eigenfreqency that belongs to the sought-after resonance.An analysis of the eigenfrequencies and the corresponding eigenvectors shows,that every electrically excitable resonance leads to very low values of the impedance in the associated eigenvector.Since the radial resonance is the first one that can be excited electrically,it can be easily found by means of this impedance value.The frequency range around the thickness resonance also contains higher order radial modes that do not allow for a unique identifiability of the thickness resonance by means of an eigenfrequency analysis.

    The second column in Table 4.1 lists the first and total order sensitivity indices with respect to the location of the radial resonance for selected parameters.The selection of the parameter is motivated by a previous parameter study Rupitsch,Sutor,Ilg,and Lerch(2010).As material a Pz27 ceramic from FERROPERM was chosen with assumed material mean-values according to Table 4.1.A check of statistical convergence revealed that after computing 100000 samples a good approximation of the sensitivity indices could be achieved.From this calculations(see Tab.4.1),we can draw directly following conclusions

    ·For the location of the radial resonator’s resonance,only the entries in thestiffness tensor show a significant sensitivity,especially

    Table 2:Sensitivity indices for the resonance frequencies of a radial resonator.

    ·The decomposition of the system’s variance is nearly linear.98%of the uncertainty can be explained by primary effects and roughly 4%are due to interaction effects between certain variables.

    ·The interaction is mainly affecting the parameter

    4.2 Electrical impedance-harmonic analysis

    In the following,the sensitivity analysis was applied to a harmonic analysis of a piezoelectric disc where both the fundamental radial and the thickness resonance are investigated.Simulations have been conducted for 400 linearly distributed frequencies around the resonances.In detail,the frequency ranges 70–90 kHz and 850–1050 kHz are considered and the resonances are given by searching for the minimum impedance within these ranges.Table 4.1 shows that the harmonic analysis results in the same sensitivities as w.r.t.the solution of the eigenvalue problem(compare columns two and three).The calculated sensitivity indices for the thickness resonance allow for the following conclusion:

    ·For the location of the thickness resonance onlyshows a significant sensitivity.

    Table 3:Coefficients of Determination Using different Bases for the Regression Model.

    ·The decomposition of the system’s variance is also approximately linear.98%of the uncertainty can be explained by primary effects and roughly 13%are due to interaction effects between certain variables.

    ·The interaction is mainly affecting the parameters

    4.3 Piezoelectric unimorph-static analysis-UA based on response surfaces

    The strategy of applying sensitivity measures on regression models is now applied to a piezoelectric unimorph modeled in 3D as it is assumed to be?too costly"for an analsyis directly on the model.The unimorph consists of two plates,a piezoelectric ceramic,Pz27,and copper sheet,which are glued together.The piezoelectric layer is polarized in the third dimension(z),the electrodes are always inx?yplane,see Figure 2.For the copper,a Young’s modulus of 9.6e+10Paand a Poisson ratio of 0.3 are assumed.

    Different bases for setting up a polynomial regression model have been tested leading to coefficients of determination as given in Table 4.3,which allows for a very good approximation of the true model behavior.For setting up the regression models,500 samples have been generated with a coefficient of variation of 5%.The sensitivity measures computed on the regression model with linear,quadratic and mixed quadratic terms are reported in the Tables 4.3 and 4.3.For the generation of the random input parameters,a coefficient of variation of 2%was assumed and 200000 samples have been evaluated on the regression model.The values in Table 4.3 belong to a model where theZdirection corresponds with the 3 direction of the material.In Table 4.3 a reorientation of the material is assumed so that theZdirection of the model is now the 1 direction of the material.As for setting up the regression models,already a series of model responses needs to be computed,a sensitivity analysis based on correlation can be conducted additionally,where the results are reported in Tables 4.3 and 4.3 in the third last column.

    Figure 2:Displacement field and electric potential of a piezoelectric unimorph polarized in z-direction.The unimorph is mechanically clamped at the left-hand side(y=0).

    The findings have shown that sensitivity analysis on regression models give quite accurate approximations of the real sensitivities,however,with a significantly reduced computational time.The analysis of the unimorph with a reoriented material serves mainly academic considerations as it shows sensitivities of the parameters for less common models where experiences about the main input parameters are less available.As we see,indices based on correlation give a fair approximation to the variance-based indices and may also be used for the purpose of factor prioritization.

    Remark:By comparing the changes in the results while evaluating new samples,one can de fine a stopping criterion.With this criterion an upper limit of samples can be determined.See Figure 3 for the convergence of the sum of all sensitivity measures for the unimorph example assuming uniformly distributed input parameters.As the sensitivity indices are approximated by stochastic sampling schemes,the results are only to be regarded as estimates of the true values.The accuracy of these estimates depends mainly on the number of simulation runs,the assumed scatter in the input parameters and the number of parameters in the model.A too small number of simulations or poorly chosen samples may give inaccurate results.

    Table 4:Sensitivity Indices of a 3 D Model unimorph for the displacements in Z(3)direction under static analysis.Z-Polarization on a response surface with mixed quadratic terms.

    Table 5:Sensitivity Indices of a 3 D Model unimorph for the displacements in Z(1)direction under static analysis.Z-Polarization on a response surface with mixed quadratic terms.

    Figure 3:Convergence of first and total order sensitivity indices,logarithmic scale.

    These are visible,e.g.when sensitivity indices close to zero have a negative sign or when someSivalues are larger than theSTivalues.If these effects are rather visible,it is recommended to repeat the sampling with a higher number of samples.Results are only reliably with respect to the second position after the decimal point.Repeated analyses and averaging would yield more robust results.

    5 Conclusion and Outlook

    From the numerical results one can draw the conclusion that for the applications discussed here,the variance in the responses of piezoelectric models might be separated in an almost additive manner according to the variances of the single input parameters.The relatively small differences in the values of first-order and totalorder sensitivity measures allow such conclusions.However,the result might be different when looking at further objectives,e.g.the quality factor or the emitted waves of a piezoelectric ultrasound generator.Variance-based methods are proven to be a valuable tool in analysing and studying piezoelectric models of moderate complexity.For larger models,the high computational costs may limit the applicability of the proposed methods.In those cases,indices computed with regression models or correlation-based measures have been shown to be an acceptable approximation.In order to obtain more precise statistical information about the model input parameter’s,a parameter estimation in the framework of Bayesian model updating should be conducted.The result of those inversions is a so calleda posterioriprobability from which all statistical moments can be retrieved.Those results will mainly be dependent on the number of measurements and the measurements’noise.

    Acknowledgement:The second and third author gratefully acknowledges the support of the German Research Foundation(DFG)in context of the Collaborative Research Centre/Transregio 39 PT-PIESA,subproject C6.

    The authors additionally want to thank the reviewers and Dr.Holger Keitel,Dr.Stephan Rupitsch and Peter Olney for their careful reading of the manuscript.

    Alwan,A.;Aluru,N.(2011):Uncertainty quantification of MEMS using a datadependent adaptive stochastic collocation method.Comput.Methods Appl.Mech.Engrg,vol.200,pp.3169–3182.

    Armstrong,J.S.(2012):Illusions in regression analysis.International Journal of Forecasting(forthcoming).

    Cukier,R.I.;Fortuin,C.M.;Shuler,K.;Petschek,A.;Schaibly,J.(59):Study of the sensitivity of coupled reaction systems to uncertainties in rate coefficients.Theory.J.Chem.Phys,vol.1973.

    Cumming,G.(2012):Understanding The New Statistics:Effect Sizes,Con fidence Intervals,and Meta-Analysis.Routledge,New York.

    Du,X.-H.;Wang,Q.-M.;Uchino,K.(2004):An accurate method for the determination of complex coefficients of single crystal piezoelectric resonators I:Theory.IEEEUFFC,vol.51,pp.227–235.

    Dziatkiewicz,G.;Fedelinski,P.(2007):Dynamic analysis of piezoelectric structures by the dual reciprocity boundary element method.Computer Modeling in Engineering&Sciences,vol.17,no.1,pp.35–46.

    Forrester,A.;Sobester,A.;Keane,A.(2008):Engineering Design via Surrogate Modelling.John Wiley and Sons Ltd.

    Homma,T.;Saltelli,A.(1996):Important measures in global sensitivity analysis of nonlinear models.Reliability Engineering And System Safety,vol.52,pp.1–17.

    IEEEStandard(1985):IEEE Standard on Piezoelectricity,1985.IEEE/ANSI.

    Iman,R.;Conover,W.(1980): Small sample sensitivity analysis techniques for computer models with an application to risk assessment.Communications in Statistics,vol.17,no.9,pp.1749–1842.

    Jonsson,U.;Andersson,B.;Lindahl,O.(2013):A FEM-based method using harmonic overtones to determine the effective elastic,dielectric,and piezoelectric parameters offreely vibrating thick piezoelectric disks.IEEE Transactions on Ultrasoncis,Ferroelectrics and Frequency Control,vol.60,no.1,pp.243–235.

    Kaltenbacher,B.;Lahmer,T.;Mohr,M.(2006):PDE based determination of piezoelectric material tensors.European Journal of Applied Mathematics.

    Kaltenbacher,M.(2004):Numerical Simulation of Mechatronic Sensors and Actuators.Springer Berlin-Heidelberg-New York.ISBN:3-540-20458-X.

    Kaltenbacher,M.(2010):Advanced simulation tool for the design of sensors and actuators.Proc.Eurosensors XXIV,vol.5,pp.597–600.

    Kami′nski,M.;Lauke,B.(2013):Parameter sensitivity and probabilistic analysis of the elastic homogenized properties for rubber filled polymers.Computer Modeling in Engineering&Sciences,vol.93,no.6,pp.411–440.

    Karaki,G.(2012):Assessment of coupled models of bridges considering timedependent vehicular loading.PhD thesis,Bauhaus University Weimar,Germany,DFG-Research Training Group 1462 ?Modelquality”,2012.

    Keitel,H.;Karaki,G.;Lahmer,T.;Nikulla,S.;Zabel,V.(2011):Evaluation of coupled partial models in structural engineering using graph theory and sensitivity analyses.Engineering Structures,vol.33,pp.3726–3736.

    Kwok,K.W.;Chan,H.L.W.;Choy,C.L.(1997):Evaluation of the material parameters of piezoelectric materials by various methods.IEEEUFFC,vol.44,no.4,pp.733–742.

    Lahmer,T.(2009): Modified landweber iterations in a multilevel algorithm applied to inverse problems in piezoelectricity.Journal of Inverse and Ill-Posed Problems,vol.17,no.1.

    Lahmer,T.;Kaltenbacher,B.;Schulz,V.(2008): Optimal experiment design for piezoelectric material tensor identification.Inverse Problems in Science and Engineering,vol.6,pp.369–387.

    Lahmer,T.;Kaltenbacher,M.;Kaltenbacher,B.;Leder,E.;Lerch,R.(2008):FEM-based determination of real and complex elastic,dielectric and piezoelectric moduli in piezoceramic materials.IEEE Transactions on Ultrasonics,Ferroelectrics,and Frequency Control,vol.55,no.2,pp.465–475.

    Lerch,R.(1990): Simulation of piezoelectric devices by two-and threedimensional finite elements.IIEEE Transactions on UFFC,vol.37,no.3,pp.233–247.ISSN 0885-3010.

    Marrel,A.;Iooss,B.;Veiga,S.;Ribatet,M.(2012):Global sensitivity analysis of stochastic computer models with joint metamodels.Statistics and Computing,vol.22,no.3,pp.833–847.

    Most,T.(2011):Assessment of structural simulation models by estimating uncertainties due to model selection and model simplification.Computers&Structures,vol.89,no.17-18,pp.1664–1672.

    Olyaie,M.S.;Razfar,M.R.;Kansa,E.J.(2011): Reliability based topology optimization of a linear piezoelectric micromotor using the cell-based smoothed finite element method.Computer Modeling in Engineering&Sciences,vol.75,no.1,pp.43–88.

    Olyaie,M.S.;Razfar,M.R.;Wang,S.;Kansa,E.J.(2011): Dynamic analysis of piezoelectric structures by the dual reciprocity boundary element method.Computer Modeling in Engineering&Sciences,vol.82,no.1,pp.55–82.

    Poteralskia,A.;Szczepanika,M.;Dziatkiewicza,G.;W.,K.;Burczyn′skiab,T.(2011):Immune identification of piezoelectric material constants using BEM.Inverse Problems in Science and Engineering,vol.19,no.1,pp.103–116.

    Rupitsch,S.;Sutor,A.;Ilg,J.;Lerch,R.(2010): Identification procedure for real and imaginary material parameters of piezoceramic materials.Ultrasonics Symposium(IUS),2010 IEEE,pp.1214–1217.

    Rupitsch,S.;Wolf,F.;Sutor,A.;Lerch,R.(2012): Reliable modeling of piezoceramic materials utilized in sensors and actuators.Acta Mechanica,pp.1–13.10.1007/s00707-012-0639-7.

    Rupitsch,S.J.;Lerch,R.(2009): Inverse method to estimate material parameters for piezoceramic disc actuators.Applied Physics A:Materials Science and Processin,vol.97,no.4,pp.735–740.

    Rus,G.;Palma,R.;Pérez-Aparicio,J.L.(2011): Experimental design of dynamic model-based damage identification in piezoelectric ceramics.Mechanical Systems and Signal Processing,vol.26,pp.268–293.

    Saltelli,A.;Ratto,M.;Andres,T.;Campolongo,F.;Cariboni,F.;Gatelli,D.;M.,S.;Tarantola,S.(2008):Global Sensitivity Analysis.J.Wiley&Sons.Ltd.

    Saltelli,A.;Tarantola,S.;Chan,K.(1999): Quantitative model-independent method for global sensitivity analysis of model output.Technometrics,vol.41,no.1,pp.39–56.

    Sherrit,S.;Wiederick,H.D.;Mukherjee,B.K.(1992): Non-iterative evaluation of the real and imaginary material monstants of piezoelectric resonators.Ferroelectrics,vol.134,pp.111–119.

    Sherrit,S.;Wiederick,H.D.;Mukherjee,B.K.(1997): A complete characterization of the piezoelectric,dielectric,and elastic properties of Motorola PZT 3203 HD,including losses and dispersion.SPIE Proceedings,vol.3037-22,pp.158–169.

    Sin,G.;Gernaey,K.;Neumann,M.;Loosdrecht,M.;Gujer,W.(2009):Uncertainty analysis in WWTP model applications:A critical discussion using an example from design.Water Research,vol.43,no.11,pp.2894–2906.

    Sobol,I.M.(1990): Sensitivity estimates for nonlinear mathematical models.Matematicheskoe Modelirovanie,vol.2,pp.112–118.

    Sorensen,D.;Lehoucq,R.;Yang,C.;Maschhoff,K.(1996):ARPACK Homepage-http://www.caam.rice.edu/software/arpack/,1996.

    Zhang,X.;Xu,X.;Wang,Y.(2014): Wave propagation in piezoelectric rods with rectangular cross sections.Computer Modeling in Engineering&Sciences,vol.100,no.1,pp.1–17.

    1Bauhaus University Weimar,Germany

    2Friedrich-Alexander University,Erlangen-Nuremberg,Germany

    欧美日韩av久久| 97在线人人人人妻| 亚洲精品aⅴ在线观看| 伊人久久大香线蕉亚洲五| 91在线精品国自产拍蜜月| 国产黄色免费在线视频| 久久精品国产亚洲av涩爱| 欧美精品亚洲一区二区| 久久国产亚洲av麻豆专区| 最新中文字幕久久久久| av免费观看日本| 日本色播在线视频| 日韩一区二区视频免费看| 久久久久久久久免费视频了| 制服人妻中文乱码| 如日韩欧美国产精品一区二区三区| 两个人免费观看高清视频| 久久久国产欧美日韩av| 少妇精品久久久久久久| 亚洲五月色婷婷综合| 黑人猛操日本美女一级片| 又大又黄又爽视频免费| 寂寞人妻少妇视频99o| h视频一区二区三区| 一本大道久久a久久精品| 欧美另类一区| 香蕉国产在线看| 伊人久久大香线蕉亚洲五| 精品第一国产精品| 国产 一区精品| 国产成人精品福利久久| 久久狼人影院| 亚洲伊人色综图| 乱人伦中国视频| 亚洲国产欧美日韩在线播放| 精品一区在线观看国产| 免费女性裸体啪啪无遮挡网站| 一边摸一边做爽爽视频免费| 婷婷色综合www| 免费黄色在线免费观看| 国产 精品1| 建设人人有责人人尽责人人享有的| 美女视频免费永久观看网站| 久久久a久久爽久久v久久| 亚洲av综合色区一区| 久久精品久久精品一区二区三区| 黑人欧美特级aaaaaa片| 久久青草综合色| 搡老乐熟女国产| 欧美少妇被猛烈插入视频| 国产男女内射视频| av.在线天堂| 国产白丝娇喘喷水9色精品| 午夜福利一区二区在线看| 黄色 视频免费看| 成人国语在线视频| 欧美亚洲 丝袜 人妻 在线| 国产精品无大码| 在线观看免费高清a一片| 成年女人毛片免费观看观看9 | 91成人精品电影| 91在线精品国自产拍蜜月| 免费av中文字幕在线| 七月丁香在线播放| 亚洲伊人久久精品综合| 日本vs欧美在线观看视频| 欧美bdsm另类| 亚洲美女视频黄频| 777米奇影视久久| 国产激情久久老熟女| 国产片内射在线| 国产成人免费观看mmmm| 日韩中字成人| 日韩中字成人| 超色免费av| 高清黄色对白视频在线免费看| 男女免费视频国产| 男女免费视频国产| 国产毛片在线视频| 欧美+日韩+精品| 亚洲精品国产av蜜桃| 麻豆精品久久久久久蜜桃| 在现免费观看毛片| 一级片'在线观看视频| 熟女少妇亚洲综合色aaa.| 国产97色在线日韩免费| 日本vs欧美在线观看视频| av网站免费在线观看视频| 亚洲欧美成人精品一区二区| 韩国av在线不卡| www日本在线高清视频| 啦啦啦在线免费观看视频4| 久久人人爽av亚洲精品天堂| 男女国产视频网站| 一级毛片 在线播放| 中文字幕人妻熟女乱码| 欧美日韩成人在线一区二区| 边亲边吃奶的免费视频| 在线天堂中文资源库| 深夜精品福利| 少妇的逼水好多| 美女午夜性视频免费| 精品亚洲乱码少妇综合久久| 久久久久久久国产电影| 久久久久久久大尺度免费视频| 男人爽女人下面视频在线观看| 国产成人a∨麻豆精品| av卡一久久| 国精品久久久久久国模美| 欧美精品av麻豆av| 成人国产麻豆网| 久久久亚洲精品成人影院| 亚洲第一区二区三区不卡| 国产精品av久久久久免费| 考比视频在线观看| 欧美成人午夜精品| 老司机影院成人| 欧美国产精品va在线观看不卡| 国产精品久久久久成人av| 七月丁香在线播放| 国产女主播在线喷水免费视频网站| 少妇的逼水好多| 不卡av一区二区三区| 日韩欧美精品免费久久| 日本爱情动作片www.在线观看| 亚洲色图综合在线观看| 久久久久视频综合| 欧美亚洲 丝袜 人妻 在线| 亚洲第一区二区三区不卡| 午夜福利一区二区在线看| 国产淫语在线视频| 午夜福利乱码中文字幕| 午夜福利乱码中文字幕| 少妇猛男粗大的猛烈进出视频| 一级毛片电影观看| 欧美日韩亚洲国产一区二区在线观看 | 天天躁夜夜躁狠狠久久av| 只有这里有精品99| 婷婷色综合大香蕉| 蜜桃在线观看..| 999精品在线视频| www日本在线高清视频| 成年人午夜在线观看视频| 国产精品久久久久久av不卡| 久久 成人 亚洲| 毛片一级片免费看久久久久| 久久久久网色| 国产乱人偷精品视频| 日韩视频在线欧美| 久久精品aⅴ一区二区三区四区 | 国产精品久久久av美女十八| 欧美 日韩 精品 国产| 啦啦啦中文免费视频观看日本| 亚洲精品国产色婷婷电影| 日韩av免费高清视频| 中文字幕人妻熟女乱码| 性少妇av在线| 午夜久久久在线观看| 日韩av免费高清视频| 高清av免费在线| 久久久久国产一级毛片高清牌| 久久影院123| 久久精品久久精品一区二区三区| 亚洲国产色片| 黄色 视频免费看| 久久久国产欧美日韩av| 伦精品一区二区三区| 国产av精品麻豆| 婷婷色麻豆天堂久久| 99香蕉大伊视频| 欧美精品av麻豆av| 在线观看免费高清a一片| 久久午夜福利片| 国产成人免费无遮挡视频| 亚洲精品国产色婷婷电影| 一边亲一边摸免费视频| 国产成人91sexporn| 国产亚洲欧美精品永久| 色哟哟·www| 一级毛片电影观看| 午夜老司机福利剧场| 亚洲欧美一区二区三区国产| 999久久久国产精品视频| 一区二区日韩欧美中文字幕| 女性生殖器流出的白浆| 国产黄色免费在线视频| av免费在线看不卡| 久久久久久久精品精品| 曰老女人黄片| 国产日韩欧美视频二区| 日韩制服丝袜自拍偷拍| 日日啪夜夜爽| 亚洲三区欧美一区| 99九九在线精品视频| 校园人妻丝袜中文字幕| 久久久精品94久久精品| 国产片内射在线| 熟女少妇亚洲综合色aaa.| 日本欧美国产在线视频| 国产成人欧美| 亚洲精品第二区| 日韩 亚洲 欧美在线| 超碰97精品在线观看| 欧美日韩亚洲高清精品| 一级a爱视频在线免费观看| 日本欧美国产在线视频| 三级国产精品片| 男人添女人高潮全过程视频| 亚洲男人天堂网一区| 九草在线视频观看| 国产男女超爽视频在线观看| 老司机影院毛片| 亚洲av欧美aⅴ国产| 最黄视频免费看| xxx大片免费视频| 国产视频首页在线观看| av电影中文网址| 亚洲国产精品成人久久小说| 精品亚洲成a人片在线观看| 精品一品国产午夜福利视频| 国产精品一区二区在线观看99| 亚洲精品自拍成人| 麻豆av在线久日| 精品国产一区二区久久| 色吧在线观看| av国产精品久久久久影院| 侵犯人妻中文字幕一二三四区| 一级黄片播放器| 在线观看免费视频网站a站| videosex国产| 久久精品亚洲av国产电影网| 两性夫妻黄色片| 最新中文字幕久久久久| 欧美另类一区| 久久国产亚洲av麻豆专区| av一本久久久久| 毛片一级片免费看久久久久| 精品国产超薄肉色丝袜足j| 国产色婷婷99| 亚洲经典国产精华液单| 亚洲一区二区三区欧美精品| 午夜久久久在线观看| 国产免费视频播放在线视频| 蜜桃国产av成人99| 丝袜喷水一区| 国产日韩欧美亚洲二区| a级片在线免费高清观看视频| 老鸭窝网址在线观看| 欧美日韩亚洲国产一区二区在线观看 | 18禁观看日本| 亚洲视频免费观看视频| 你懂的网址亚洲精品在线观看| 视频在线观看一区二区三区| 日日撸夜夜添| av国产精品久久久久影院| 精品久久久精品久久久| 建设人人有责人人尽责人人享有的| 男人舔女人的私密视频| 国产成人精品一,二区| 亚洲国产精品999| 亚洲欧美一区二区三区黑人 | 午夜91福利影院| 美女国产高潮福利片在线看| 国产爽快片一区二区三区| 精品视频人人做人人爽| kizo精华| 人体艺术视频欧美日本| 两性夫妻黄色片| 久久久久国产一级毛片高清牌| 搡老乐熟女国产| 免费日韩欧美在线观看| 国产精品嫩草影院av在线观看| 国产精品国产三级国产专区5o| 美女视频免费永久观看网站| 亚洲欧美色中文字幕在线| 国产探花极品一区二区| 男的添女的下面高潮视频| 久久久久久伊人网av| 国产1区2区3区精品| 精品久久久精品久久久| 婷婷色av中文字幕| 国产xxxxx性猛交| 精品酒店卫生间| 亚洲av中文av极速乱| 亚洲av成人精品一二三区| 欧美国产精品一级二级三级| 日本欧美国产在线视频| 欧美 日韩 精品 国产| 欧美97在线视频| 久久久久国产网址| 天天躁夜夜躁狠狠久久av| 极品人妻少妇av视频| 一级毛片电影观看| av天堂久久9| 午夜激情久久久久久久| 两个人看的免费小视频| 最近手机中文字幕大全| 少妇人妻精品综合一区二区| 三级国产精品片| 国产免费又黄又爽又色| 亚洲精华国产精华液的使用体验| 伦理电影大哥的女人| 天天操日日干夜夜撸| 欧美日韩视频精品一区| 国产成人91sexporn| 亚洲综合色网址| 日日摸夜夜添夜夜爱| 亚洲综合精品二区| 一级毛片我不卡| 少妇熟女欧美另类| 熟妇人妻不卡中文字幕| 在线观看一区二区三区激情| 久久这里有精品视频免费| 国产一区二区激情短视频 | 日韩视频在线欧美| 在线精品无人区一区二区三| 免费女性裸体啪啪无遮挡网站| 国产97色在线日韩免费| 久久久久网色| 91成人精品电影| 欧美激情极品国产一区二区三区| 午夜日本视频在线| 婷婷色综合大香蕉| 久久综合国产亚洲精品| 国产精品一二三区在线看| 免费看不卡的av| 99精国产麻豆久久婷婷| 制服诱惑二区| 国产片内射在线| 亚洲欧美日韩另类电影网站| 日韩一区二区视频免费看| 国产av精品麻豆| 七月丁香在线播放| 大片免费播放器 马上看| 伊人久久国产一区二区| 欧美+日韩+精品| 青草久久国产| 在线观看国产h片| 久久久久久人人人人人| 人体艺术视频欧美日本| 七月丁香在线播放| 精品国产国语对白av| 成年av动漫网址| 人人妻人人爽人人添夜夜欢视频| 波多野结衣一区麻豆| 亚洲国产欧美网| 久久久a久久爽久久v久久| 中文天堂在线官网| av天堂久久9| www.自偷自拍.com| 在线看a的网站| 黄片播放在线免费| 女人高潮潮喷娇喘18禁视频| 精品福利永久在线观看| 男女国产视频网站| 久久鲁丝午夜福利片| 成人毛片60女人毛片免费| 婷婷色综合大香蕉| 丰满乱子伦码专区| a级毛片黄视频| 国产成人午夜福利电影在线观看| 欧美精品一区二区免费开放| 欧美人与性动交α欧美软件| 国产淫语在线视频| 亚洲欧美成人综合另类久久久| 天天躁日日躁夜夜躁夜夜| 欧美日韩亚洲高清精品| 国产白丝娇喘喷水9色精品| 哪个播放器可以免费观看大片| 人妻人人澡人人爽人人| 老女人水多毛片| 制服丝袜香蕉在线| 久久ye,这里只有精品| 精品国产超薄肉色丝袜足j| 久久国产精品大桥未久av| 国产精品三级大全| 国产一级毛片在线| 亚洲一码二码三码区别大吗| 欧美精品av麻豆av| 免费黄网站久久成人精品| 成人亚洲精品一区在线观看| 欧美少妇被猛烈插入视频| 亚洲美女搞黄在线观看| 久久久亚洲精品成人影院| 青青草视频在线视频观看| 欧美成人午夜免费资源| 欧美精品人与动牲交sv欧美| 啦啦啦中文免费视频观看日本| 亚洲少妇的诱惑av| 亚洲激情五月婷婷啪啪| 丰满少妇做爰视频| 熟女av电影| 久久久精品国产亚洲av高清涩受| 丝袜美足系列| 老熟女久久久| 国产成人精品无人区| 久久久久国产一级毛片高清牌| 青春草国产在线视频| 国产精品二区激情视频| freevideosex欧美| 黑人欧美特级aaaaaa片| 晚上一个人看的免费电影| 久久这里有精品视频免费| 只有这里有精品99| 99九九在线精品视频| 日韩在线高清观看一区二区三区| 看非洲黑人一级黄片| 最近手机中文字幕大全| 免费看不卡的av| 国产一区二区在线观看av| 亚洲四区av| 亚洲在久久综合| 亚洲av男天堂| 最近2019中文字幕mv第一页| 亚洲男人天堂网一区| 美女国产视频在线观看| 在现免费观看毛片| 1024香蕉在线观看| 汤姆久久久久久久影院中文字幕| 亚洲精品av麻豆狂野| 在线观看人妻少妇| 午夜福利网站1000一区二区三区| 国产精品 欧美亚洲| 麻豆精品久久久久久蜜桃| 一本一本久久a久久精品综合妖精 国产伦在线观看视频一区 | av免费观看日本| 欧美人与性动交α欧美软件| 国产成人一区二区在线| 少妇人妻 视频| 女性被躁到高潮视频| av.在线天堂| 久久精品国产综合久久久| 美女高潮到喷水免费观看| 免费观看无遮挡的男女| 午夜免费男女啪啪视频观看| 日韩av免费高清视频| 在现免费观看毛片| 大香蕉久久网| 中文精品一卡2卡3卡4更新| 亚洲国产成人一精品久久久| 在线天堂最新版资源| 成人18禁高潮啪啪吃奶动态图| 极品少妇高潮喷水抽搐| 国产成人欧美| 久久精品国产亚洲av高清一级| 久久久国产一区二区| a 毛片基地| 亚洲中文av在线| 久久久久国产一级毛片高清牌| 欧美日韩精品网址| 黄色配什么色好看| 亚洲av男天堂| 午夜免费鲁丝| 午夜福利影视在线免费观看| 国产不卡av网站在线观看| 国产精品一区二区在线观看99| 伊人亚洲综合成人网| 国产精品女同一区二区软件| av女优亚洲男人天堂| 亚洲美女搞黄在线观看| 欧美xxⅹ黑人| 亚洲色图 男人天堂 中文字幕| 男女高潮啪啪啪动态图| 精品国产超薄肉色丝袜足j| 日本欧美视频一区| 人妻一区二区av| 免费女性裸体啪啪无遮挡网站| 一级毛片电影观看| a级毛片黄视频| 精品国产国语对白av| 午夜福利一区二区在线看| 18禁国产床啪视频网站| 国产成人精品久久二区二区91 | 中文字幕人妻丝袜制服| 一二三四中文在线观看免费高清| 亚洲伊人色综图| 国产在线视频一区二区| 爱豆传媒免费全集在线观看| 国产片特级美女逼逼视频| 亚洲国产日韩一区二区| 国产极品天堂在线| 欧美中文综合在线视频| 咕卡用的链子| 国产欧美亚洲国产| 成人毛片a级毛片在线播放| 18禁裸乳无遮挡动漫免费视频| 国产野战对白在线观看| 你懂的网址亚洲精品在线观看| av天堂久久9| 国产精品99久久99久久久不卡 | 激情五月婷婷亚洲| 99久久精品国产国产毛片| 日韩在线高清观看一区二区三区| 大话2 男鬼变身卡| 久久亚洲国产成人精品v| 日韩人妻精品一区2区三区| 丰满迷人的少妇在线观看| 一区二区三区激情视频| 男女啪啪激烈高潮av片| 色播在线永久视频| 黄色配什么色好看| 日韩av免费高清视频| 午夜影院在线不卡| 国产 一区精品| 另类精品久久| 99精国产麻豆久久婷婷| 亚洲在久久综合| 色婷婷av一区二区三区视频| 免费在线观看视频国产中文字幕亚洲 | 成年人免费黄色播放视频| 黄片播放在线免费| 亚洲伊人色综图| 欧美精品一区二区大全| 国语对白做爰xxxⅹ性视频网站| 免费在线观看完整版高清| 国产一区二区 视频在线| 丝袜喷水一区| 久久这里只有精品19| 水蜜桃什么品种好| 纵有疾风起免费观看全集完整版| 观看美女的网站| 国产毛片在线视频| 日日撸夜夜添| 国产精品99久久99久久久不卡 | 午夜福利乱码中文字幕| 波多野结衣一区麻豆| 国产一级毛片在线| 免费观看av网站的网址| 欧美少妇被猛烈插入视频| 国产精品久久久av美女十八| 丝袜人妻中文字幕| 国产精品熟女久久久久浪| 午夜福利一区二区在线看| av免费在线看不卡| 夜夜骑夜夜射夜夜干| 伦理电影大哥的女人| 国产精品亚洲av一区麻豆 | 国产熟女欧美一区二区| 精品一区二区免费观看| 99久久综合免费| 亚洲国产毛片av蜜桃av| 伊人亚洲综合成人网| 国产片内射在线| 国产免费又黄又爽又色| 欧美日韩成人在线一区二区| 免费av中文字幕在线| 亚洲av综合色区一区| 欧美日韩综合久久久久久| 精品第一国产精品| 晚上一个人看的免费电影| 视频在线观看一区二区三区| 综合色丁香网| 精品亚洲成国产av| 日韩伦理黄色片| 在现免费观看毛片| 色网站视频免费| 精品午夜福利在线看| 一级黄片播放器| 亚洲精品自拍成人| 亚洲成人一二三区av| 在线精品无人区一区二区三| 人人妻人人添人人爽欧美一区卜| 黄片小视频在线播放| 亚洲美女视频黄频| 丝袜在线中文字幕| 国产男女超爽视频在线观看| 国产黄频视频在线观看| 亚洲精品成人av观看孕妇| 国产成人一区二区在线| 最黄视频免费看| 精品国产超薄肉色丝袜足j| 高清视频免费观看一区二区| 男女国产视频网站| 两个人免费观看高清视频| 国产精品人妻久久久影院| 久久这里有精品视频免费| 欧美日韩一级在线毛片| 天天躁夜夜躁狠狠躁躁| 精品一区二区免费观看| 黑人巨大精品欧美一区二区蜜桃| 欧美精品一区二区大全| 欧美人与性动交α欧美软件| 亚洲,欧美精品.| 丰满乱子伦码专区| 久久精品熟女亚洲av麻豆精品| 日韩av不卡免费在线播放| 色婷婷久久久亚洲欧美| 在线精品无人区一区二区三| 国产视频首页在线观看| 久久精品国产a三级三级三级| 国产国语露脸激情在线看| 国产精品久久久久久精品电影小说| 狠狠婷婷综合久久久久久88av| 伦理电影免费视频| 999久久久国产精品视频| 亚洲美女黄色视频免费看| 国产精品久久久久成人av| 香蕉丝袜av| 天美传媒精品一区二区| 伦精品一区二区三区| 少妇人妻久久综合中文| 欧美精品亚洲一区二区| 午夜福利网站1000一区二区三区| 亚洲国产精品一区三区| 少妇被粗大的猛进出69影院| 日韩不卡一区二区三区视频在线| 久久久久久人人人人人| 老鸭窝网址在线观看| 女性被躁到高潮视频| 国产精品二区激情视频| 亚洲欧美成人综合另类久久久| 999久久久国产精品视频| 欧美精品一区二区大全| 午夜激情av网站| 香蕉丝袜av|