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

    日韩高清综合在线| 久久6这里有精品| 国产男人的电影天堂91| eeuss影院久久| www.色视频.com| 婷婷六月久久综合丁香| 国产成人福利小说| 一级二级三级毛片免费看| 日本猛色少妇xxxxx猛交久久| 久久久成人免费电影| 欧美日韩在线观看h| 日韩欧美精品v在线| or卡值多少钱| 国产成人aa在线观看| 美女cb高潮喷水在线观看| 午夜精品在线福利| 欧美xxxx性猛交bbbb| 三级毛片av免费| 夜夜爽夜夜爽视频| 国产精品一区二区性色av| 亚洲美女视频黄频| 一个人观看的视频www高清免费观看| 久久久国产成人精品二区| 精品国内亚洲2022精品成人| 身体一侧抽搐| av免费在线看不卡| 丰满少妇做爰视频| 午夜精品一区二区三区免费看| 亚洲成人av在线免费| 26uuu在线亚洲综合色| 免费大片18禁| 桃色一区二区三区在线观看| 日韩成人伦理影院| av免费观看日本| 一级毛片久久久久久久久女| 久久精品国产99精品国产亚洲性色| 国产精品人妻久久久久久| 男女啪啪激烈高潮av片| 激情 狠狠 欧美| 国产成人aa在线观看| 可以在线观看毛片的网站| 国产极品精品免费视频能看的| 最近最新中文字幕免费大全7| 一区二区三区免费毛片| 中文在线观看免费www的网站| 欧美激情久久久久久爽电影| 精品久久久久久久久亚洲| 亚洲av一区综合| 青春草国产在线视频| 看黄色毛片网站| 日韩av不卡免费在线播放| 少妇的逼好多水| 1024手机看黄色片| 人人妻人人看人人澡| 欧美一区二区精品小视频在线| 黑人高潮一二区| 国产在视频线精品| 国产免费又黄又爽又色| 日韩欧美 国产精品| 久久精品影院6| 国产老妇伦熟女老妇高清| 啦啦啦观看免费观看视频高清| 国产成人精品久久久久久| 爱豆传媒免费全集在线观看| 国产精品一区www在线观看| 国产探花在线观看一区二区| 国产午夜精品久久久久久一区二区三区| 亚洲18禁久久av| 能在线免费看毛片的网站| 国产视频内射| 丝袜喷水一区| 波多野结衣巨乳人妻| av专区在线播放| 国产探花极品一区二区| 热99在线观看视频| 国产亚洲精品久久久com| 国产色爽女视频免费观看| 天堂中文最新版在线下载 | 特大巨黑吊av在线直播| av天堂中文字幕网| 欧美精品国产亚洲| 久久久久网色| 中文字幕av在线有码专区| 亚洲激情五月婷婷啪啪| 97在线视频观看| 99热6这里只有精品| 日韩欧美在线乱码| 黄片无遮挡物在线观看| 男女视频在线观看网站免费| 久久久精品大字幕| 亚洲av中文av极速乱| 91久久精品国产一区二区三区| 热99re8久久精品国产| 欧美高清成人免费视频www| 自拍偷自拍亚洲精品老妇| 欧美性猛交黑人性爽| 一级黄片播放器| 国产成人一区二区在线| 久久99精品国语久久久| 亚洲欧美日韩无卡精品| АⅤ资源中文在线天堂| 亚洲国产精品成人综合色| 高清毛片免费看| 三级男女做爰猛烈吃奶摸视频| 亚洲av男天堂| 少妇人妻精品综合一区二区| 国产伦一二天堂av在线观看| 国产成人免费观看mmmm| 最后的刺客免费高清国语| 99久久成人亚洲精品观看| 免费av不卡在线播放| 97超视频在线观看视频| 最近的中文字幕免费完整| 99久久精品一区二区三区| 午夜福利在线观看免费完整高清在| av又黄又爽大尺度在线免费看 | 白带黄色成豆腐渣| 国产av不卡久久| 少妇熟女欧美另类| 人人妻人人澡人人爽人人夜夜 | 久久久久久大精品| 国产精品伦人一区二区| 亚洲中文字幕日韩| 18禁在线播放成人免费| 精品久久久久久久久久久久久| 日本免费在线观看一区| 亚洲最大成人av| 亚洲欧美精品综合久久99| h日本视频在线播放| 日韩大片免费观看网站 | 插逼视频在线观看| 五月伊人婷婷丁香| 高清av免费在线| 免费电影在线观看免费观看| 综合色丁香网| 亚洲欧美精品综合久久99| av天堂中文字幕网| 国国产精品蜜臀av免费| 国产一级毛片在线| 亚洲电影在线观看av| 床上黄色一级片| 秋霞在线观看毛片| 三级经典国产精品| 亚洲熟妇中文字幕五十中出| 干丝袜人妻中文字幕| 赤兔流量卡办理| 亚洲国产高清在线一区二区三| 国产精品三级大全| 国产精品伦人一区二区| 午夜福利网站1000一区二区三区| 中文字幕熟女人妻在线| 国产在视频线在精品| av在线蜜桃| or卡值多少钱| 久久久a久久爽久久v久久| 91午夜精品亚洲一区二区三区| 综合色丁香网| 精品国内亚洲2022精品成人| 亚洲欧洲国产日韩| 日韩亚洲欧美综合| 日韩,欧美,国产一区二区三区 | 国产伦精品一区二区三区四那| 国产精品无大码| 午夜福利在线观看吧| 亚洲三级黄色毛片| 国产亚洲5aaaaa淫片| 性插视频无遮挡在线免费观看| 老师上课跳d突然被开到最大视频| 国产精品一区二区性色av| 亚洲av电影在线观看一区二区三区 | 26uuu在线亚洲综合色| 亚洲四区av| 欧美变态另类bdsm刘玥| 99久久精品热视频| 欧美丝袜亚洲另类| 成年女人看的毛片在线观看| 91久久精品国产一区二区三区| 午夜精品国产一区二区电影 | 久久人妻av系列| 国产免费一级a男人的天堂| 久久久精品欧美日韩精品| 日本三级黄在线观看| 身体一侧抽搐| 免费搜索国产男女视频| 日韩在线高清观看一区二区三区| 国产私拍福利视频在线观看| 成人鲁丝片一二三区免费| 插逼视频在线观看| 舔av片在线| 日本三级黄在线观看| 午夜福利高清视频| 亚洲精品国产av成人精品| 国产在线男女| 久久精品国产鲁丝片午夜精品| 欧美日本视频| 国产成人免费观看mmmm| 一级黄片播放器| av福利片在线观看| 色哟哟·www| 欧美极品一区二区三区四区| 日本五十路高清| 亚洲精品影视一区二区三区av| 成人综合一区亚洲| 国内揄拍国产精品人妻在线| 国产真实乱freesex| 九草在线视频观看| 青春草亚洲视频在线观看| 直男gayav资源| 国产乱人视频| 亚洲四区av| 一个人看视频在线观看www免费| 久久久久久久久久黄片| 亚洲av电影不卡..在线观看| 亚洲精品自拍成人| 99久久精品国产国产毛片| 亚洲无线观看免费| 国产一区亚洲一区在线观看| 免费看日本二区| 国产免费一级a男人的天堂| 国产高清有码在线观看视频| 午夜免费男女啪啪视频观看| 亚洲精品色激情综合| 婷婷色av中文字幕| 神马国产精品三级电影在线观看| 国模一区二区三区四区视频| 免费看光身美女| 国产一区亚洲一区在线观看| 亚洲精品乱码久久久v下载方式| 91久久精品电影网| 国语对白做爰xxxⅹ性视频网站| 18+在线观看网站| 精品一区二区免费观看| 久久精品国产99精品国产亚洲性色| 寂寞人妻少妇视频99o| 国产高清有码在线观看视频| 尾随美女入室| 97超视频在线观看视频| 91精品国产九色| 欧美一区二区精品小视频在线| av卡一久久| 国产又色又爽无遮挡免| 久久久久久久久大av| 美女高潮的动态| 精品一区二区三区人妻视频| 色播亚洲综合网| 午夜视频国产福利| 国产精品人妻久久久影院| 晚上一个人看的免费电影| 国产女主播在线喷水免费视频网站 | 欧美成人一区二区免费高清观看| 一个人看视频在线观看www免费| 国产成人精品婷婷| 精品99又大又爽又粗少妇毛片| 我的老师免费观看完整版| 岛国在线免费视频观看| 日韩成人伦理影院| 男女国产视频网站| 国产欧美日韩精品一区二区| 久久精品久久久久久噜噜老黄 | 亚洲精品成人久久久久久| 国产亚洲av嫩草精品影院| 亚洲在久久综合| 国产精品日韩av在线免费观看| 久久久久久久久久久丰满| 国产av在哪里看| 国产亚洲av嫩草精品影院| 亚洲最大成人中文| 国产精品爽爽va在线观看网站| 国产av一区在线观看免费| 免费看a级黄色片| 中文乱码字字幕精品一区二区三区 | 亚洲国产高清在线一区二区三| 日日干狠狠操夜夜爽| 99久久精品一区二区三区| 99久久中文字幕三级久久日本| 欧美3d第一页| 国产三级在线视频| 久久精品夜色国产| 中文字幕av成人在线电影| 久久久久九九精品影院| 精品国产一区二区三区久久久樱花 | 嘟嘟电影网在线观看| 男的添女的下面高潮视频| 能在线免费看毛片的网站| 看免费成人av毛片| 少妇裸体淫交视频免费看高清| 免费观看精品视频网站| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 一区二区三区四区激情视频| 亚洲无线观看免费| 丰满乱子伦码专区| 麻豆一二三区av精品| 九九久久精品国产亚洲av麻豆| 国产真实乱freesex| 亚洲国产日韩欧美精品在线观看| 日韩,欧美,国产一区二区三区 | 日本av手机在线免费观看| 国产视频内射| 精品国产一区二区三区久久久樱花 | 亚洲国产欧美人成| 特大巨黑吊av在线直播| 免费观看在线日韩| 变态另类丝袜制服| 99热这里只有是精品在线观看| 老司机影院毛片| 国产高清国产精品国产三级 | 男人狂女人下面高潮的视频| 欧美+日韩+精品| 搡女人真爽免费视频火全软件| 国产亚洲av嫩草精品影院| 亚洲人成网站高清观看| 亚洲激情五月婷婷啪啪| 最近中文字幕高清免费大全6| 欧美高清性xxxxhd video| 你懂的网址亚洲精品在线观看 | 久久久欧美国产精品| 亚洲精品成人久久久久久| 日韩欧美精品免费久久| 美女被艹到高潮喷水动态| 亚洲自拍偷在线| 一级毛片我不卡| 久久久久久久久中文| 欧美bdsm另类| 久久久国产成人免费| av在线亚洲专区| 欧美一区二区国产精品久久精品| 亚洲人成网站在线观看播放| 久久久久久久午夜电影| 久久久精品大字幕| 高清av免费在线| 国产成人精品一,二区| 国产成人91sexporn| 精品久久久久久成人av| 亚洲美女视频黄频| 尤物成人国产欧美一区二区三区| 亚洲av.av天堂| 99久久人妻综合| 国内揄拍国产精品人妻在线| 色尼玛亚洲综合影院| 亚洲av成人精品一二三区| 亚洲av电影不卡..在线观看| 热99在线观看视频| 国产毛片a区久久久久| 最近中文字幕高清免费大全6| 久久久久久久久中文| 国国产精品蜜臀av免费| 18+在线观看网站| 成人特级av手机在线观看| 天堂av国产一区二区熟女人妻| 亚洲国产欧美在线一区| 欧美97在线视频| 亚洲国产欧洲综合997久久,| 亚洲丝袜综合中文字幕| 长腿黑丝高跟| 精品久久久噜噜| 七月丁香在线播放| 丝袜喷水一区| 久久精品影院6| 啦啦啦观看免费观看视频高清| 精品午夜福利在线看| 观看免费一级毛片| 久久人人爽人人爽人人片va| 国产激情偷乱视频一区二区| 九九热线精品视视频播放| 波多野结衣巨乳人妻| 国产高清视频在线观看网站| 午夜精品一区二区三区免费看| 亚洲自偷自拍三级| 日韩一区二区视频免费看| 色哟哟·www| 国产激情偷乱视频一区二区| 乱人视频在线观看| 日本熟妇午夜| 亚洲五月天丁香| 国产高清不卡午夜福利| 大香蕉97超碰在线| 永久免费av网站大全| 精品欧美国产一区二区三| 成人毛片60女人毛片免费| 在线免费十八禁| 国产女主播在线喷水免费视频网站 | 精品无人区乱码1区二区| 午夜日本视频在线| 99热这里只有是精品在线观看| 自拍偷自拍亚洲精品老妇| 干丝袜人妻中文字幕| 97人妻精品一区二区三区麻豆| 99久国产av精品| 日韩精品有码人妻一区| 高清在线视频一区二区三区 | 九九爱精品视频在线观看| 日本wwww免费看| 高清毛片免费看| 99热网站在线观看| 黄色欧美视频在线观看| 亚洲欧洲日产国产| 欧美日韩国产亚洲二区| 久久精品久久精品一区二区三区| 夫妻性生交免费视频一级片| 国产精品熟女久久久久浪| 我的老师免费观看完整版| 熟女人妻精品中文字幕| 日本五十路高清| 国产精品伦人一区二区| 九九久久精品国产亚洲av麻豆| 天堂影院成人在线观看| 亚洲精品色激情综合| 免费人成在线观看视频色| 成人特级av手机在线观看| 色5月婷婷丁香| 一个人免费在线观看电影| 天天一区二区日本电影三级| 成人无遮挡网站| 亚洲欧美精品自产自拍| 熟女人妻精品中文字幕| 欧美日韩在线观看h| 精品99又大又爽又粗少妇毛片| 菩萨蛮人人尽说江南好唐韦庄 | 3wmmmm亚洲av在线观看| 中文字幕人妻熟人妻熟丝袜美| 美女黄网站色视频| 天堂av国产一区二区熟女人妻| 青春草亚洲视频在线观看| 两性午夜刺激爽爽歪歪视频在线观看| 免费搜索国产男女视频| 午夜福利在线观看免费完整高清在| av女优亚洲男人天堂| 国产爱豆传媒在线观看| 久久久国产成人精品二区| 久久精品久久久久久噜噜老黄 | 国产精品爽爽va在线观看网站| 高清视频免费观看一区二区 | 干丝袜人妻中文字幕| 26uuu在线亚洲综合色| 国产精品久久久久久久久免| 国产成人免费观看mmmm| 99久久精品一区二区三区| 美女xxoo啪啪120秒动态图| 欧美+日韩+精品| 少妇猛男粗大的猛烈进出视频 | 亚洲经典国产精华液单| 美女cb高潮喷水在线观看| 中文天堂在线官网| 国产亚洲一区二区精品| 亚洲av成人精品一二三区| 久久久久久久亚洲中文字幕| 一级毛片电影观看 | 男女边吃奶边做爰视频| 亚洲性久久影院| 一区二区三区四区激情视频| 日本黄色视频三级网站网址| 秋霞在线观看毛片| 少妇丰满av| 天堂影院成人在线观看| 久久精品综合一区二区三区| 少妇熟女aⅴ在线视频| 国产精品久久视频播放| 欧美一区二区精品小视频在线| 欧美日韩综合久久久久久| 一区二区三区乱码不卡18| 亚洲色图av天堂| av在线播放精品| 一个人看视频在线观看www免费| 青春草亚洲视频在线观看| 日韩欧美在线乱码| a级毛片免费高清观看在线播放| 国产免费又黄又爽又色| 国产精品一区二区三区四区久久| 久久久久久久国产电影| 蜜臀久久99精品久久宅男| 久久久久网色| 国产免费福利视频在线观看| 欧美xxxx黑人xx丫x性爽| 国产白丝娇喘喷水9色精品| 午夜亚洲福利在线播放| 免费av观看视频| 丰满乱子伦码专区| 最近最新中文字幕大全电影3| 婷婷六月久久综合丁香| 亚洲欧美精品综合久久99| 免费观看a级毛片全部| 国产亚洲一区二区精品| 波多野结衣高清无吗| 国产成人aa在线观看| 国产精品综合久久久久久久免费| 亚洲最大成人手机在线| 国产免费福利视频在线观看| 日韩亚洲欧美综合| 国产免费一级a男人的天堂| 亚洲精品国产av成人精品| 老司机影院毛片| videossex国产| 夫妻性生交免费视频一级片| 少妇猛男粗大的猛烈进出视频 | 国产黄色视频一区二区在线观看 | 1024手机看黄色片| 大话2 男鬼变身卡| 国产精华一区二区三区| 人体艺术视频欧美日本| 午夜免费激情av| 亚洲国产高清在线一区二区三| 禁无遮挡网站| 国产免费一级a男人的天堂| 久99久视频精品免费| 国产精品一区二区性色av| 日韩亚洲欧美综合| 国产精品一区二区性色av| 国产精品久久久久久精品电影小说 | 成人毛片a级毛片在线播放| 久久久国产成人免费| 国产三级中文精品| 国产成人aa在线观看| 在线观看av片永久免费下载| 日韩 亚洲 欧美在线| 国产极品天堂在线| 秋霞在线观看毛片| av播播在线观看一区| 亚洲四区av| 亚洲电影在线观看av| 在现免费观看毛片| 日本免费一区二区三区高清不卡| 少妇熟女欧美另类| 91av网一区二区| 精品99又大又爽又粗少妇毛片| 色吧在线观看| 欧美一级a爱片免费观看看| 男女边吃奶边做爰视频| 亚洲国产日韩欧美精品在线观看| 91在线精品国自产拍蜜月| 国产探花极品一区二区| 日韩av在线大香蕉| 夫妻性生交免费视频一级片| 性色avwww在线观看| 亚洲国产高清在线一区二区三| 久久久久九九精品影院| 午夜福利在线观看免费完整高清在| 老司机影院成人| 纵有疾风起免费观看全集完整版 | 一区二区三区四区激情视频| 女人久久www免费人成看片 | 国产老妇女一区| 国产精华一区二区三区| 两个人的视频大全免费| 欧美色视频一区免费| 久久热精品热| 久久精品国产亚洲网站| 国产一区二区三区av在线| 国产91av在线免费观看| 亚洲av.av天堂| 夜夜爽夜夜爽视频| 久久久久久久午夜电影| 久久久久免费精品人妻一区二区| 国产黄片视频在线免费观看| 国产精品国产高清国产av| 亚洲国产最新在线播放| 久久亚洲国产成人精品v| 国产精品爽爽va在线观看网站| 一个人看视频在线观看www免费| 国产精品久久久久久av不卡| 全区人妻精品视频| 男插女下体视频免费在线播放| 人体艺术视频欧美日本| 1024手机看黄色片| 日本免费一区二区三区高清不卡| 一级黄片播放器| 亚洲精品成人久久久久久| 国产av一区在线观看免费| 亚洲18禁久久av| 久久久欧美国产精品| 日本免费a在线| 亚洲av二区三区四区| 午夜福利视频1000在线观看| 国产免费一级a男人的天堂| 少妇熟女欧美另类| 日本午夜av视频| 在线播放国产精品三级| 久久人人爽人人爽人人片va| 日韩欧美在线乱码| 一区二区三区四区激情视频| 欧美区成人在线视频| 国产淫片久久久久久久久| 国内少妇人妻偷人精品xxx网站| 久久精品久久久久久久性| 亚洲自拍偷在线| 国产av一区在线观看免费| 天堂√8在线中文| 欧美日韩国产亚洲二区| 91精品伊人久久大香线蕉| 寂寞人妻少妇视频99o| 免费在线观看成人毛片| av免费观看日本| 国产大屁股一区二区在线视频| 国产av不卡久久| 国产伦理片在线播放av一区| 国产老妇伦熟女老妇高清| 婷婷六月久久综合丁香| 国产又色又爽无遮挡免| 亚洲四区av| 亚洲欧美成人综合另类久久久 | 男女视频在线观看网站免费| 国产av在哪里看| 大又大粗又爽又黄少妇毛片口| 亚洲五月天丁香| 国产一区二区三区av在线| 久久这里有精品视频免费| 国语对白做爰xxxⅹ性视频网站| 亚洲成人中文字幕在线播放| 赤兔流量卡办理| 亚洲四区av| 日本wwww免费看| 搞女人的毛片| 亚洲成人久久爱视频| 亚洲国产精品专区欧美| 免费看光身美女| 韩国av在线不卡|