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

    Frequency domain analysis of pre-stressed elastomeric vibration isolators

    2023-07-31 13:30:16SomnthMrimuthuShnkrKrishnpilli
    Defence Technology 2023年7期

    S.Somnth ,R.Mrimuthu ,Shnkr Krishnpilli

    a Indian Space Research Organisation, Anthariksh Bhavan, New BEL Road, Bengaluru, 560231, India

    b STR/VSSC, Indian Space Research Organisation, Thiruvananthapuram, 695022, India

    c Department of Mechanical Engineering, IIT Madras, Chennai, 600036, India

    Keywords:Elastomer Isolator Axi-symmetric Frequency Random Response Strain amplitude Complex modulus Transmissibility

    ABSTRACT Two types of elastomeric vibration isolators used for equipment vibration isolation in aerospace vehicles are considered for the present study.These isolators are constructed using elastomers mounted in steel encasings.These isolators are initially deformed statically and dynamic loads are applied on the deformed configuration.To capture the static deformation,equivalent static load corresponding to its load rating and specified displacements are created.Static deformation is computed using Finite Element methods with four node axi-symmetric element which include the geometric non-linear effect for steel and with standard Yeoh hyper-elastic material model for elastomers(Muhammed and Zu,2012)[1].Yeoh material constants are derived from uni-axial tension test data of the elastomer specimen.These isolators are subjected to harmonic and random excitations in the pre-deformed state.For numerical analysis,elastomeric constants at dynamic conditions are obtained as complex function of frequency using Dynamic Mechanical Analyzer (DMA) for a range of frequencies.The standard material model of Yeoh is modified incorporating frequency dependant material characteristics and damping in the range of frequencies of interest.A multiplicative non-separable variables law is derived for Yeoh material model to include the effect of static pre-stress,based on the methodology given in literature(Nashif et al.,1985;Beda et al.,2014)[2,3].The modifications of Yeoh model suitable for frequency domain analysis is the novelty in the present study.In the analysis,while dynamic loads are applied,the configuration is updated considering initial static loading.The frequency response of the isolators is computed using material properties evaluated at progressive dynamic strains until a match in natural frequency is observed.Appropriate damping corrections are then incorporated to match the test observed transmissibility.Then updated material properties are used to compute the random response which showed good agreement with results of experiments,validating the approach taken for the development of this model.

    1.Introduction

    Rubber based vibration isolators play an important role in vibration control,and have been widely used in machine components for their low stiffness and suitable damping properties at relatively low costs.They have been employed in numerous applications including aerospace,civil and railway structures,automotive parts[4,5] and seismic load mitigation [6].Rubber-based vibration isolation systems being passive in nature,are less complex,cheaper and easier to manufacture than their active counterparts.Rubbers are essentially long chain molecules or polymers,and exhibit elastic behavior.Hence,the term ‘Elastomer’ is used synonymously with rubber [7].Elastomers generally exhibit large deformations when subjected to large quasi-static loads,which exceed linear elastic theory,and also stress softening [8].In addition,elastomers also exhibit loss of stiffness under the effect of cyclic loading,due to which the stress-strain curve of an elastomer varies after each cycle of loading.This is called the Mullins effect[9].Den Hartog[10]was one of the first researchers to suggest the use of isolators,absorbers and Dampers for vibration isolation.Finite element method plays an increasingly important role in the design of such isolation systems.Due to their behavior,Non-linear finite element method is required to perform a computational analysis of components containing elastomers [1,11-17].Over the years,many hyper-elastic material models have been formulated to model the large deformation characteristics of elastomers.These models define the strain energy density for an elastomeric material in terms of the strain invariants and material constants [18-21].In the present study,the Yeoh hyper-elastic material model is used to represent the rubber material.Finite element formulation for hyper-elastic material is represented through Total Lagrangian (TL) or Updated Lagrangian (UL) frame of reference.In TL formulation,the deformation is always with respect to the original undeformed configurations.The stress and strain measures are computed with reference to undeformed frame,which can be transformed to the present deformed configuration using continuum mechanics principles.UL formulation is derived with respect to previously converged deformed configuration.Mathematically both the formulations are equivalent and used depending upon convenience.Presently TL formulation is considered to study the non-linear static deformation of the isolator.

    The dynamic properties of elastomers are seen to depend highly on the frequency of excitation,temperature and strain amplitude.Many experimental studies have been conducted to study this dependence and have been well documented,for example,by De Wilde et al.[22] and Jurado et al.[23].Elastomers exhibit a combination of viscous fluid behavior and elastic solid behavior,termed as visco-elasticity.Visco-elastic materials show a time-dependent behavior where the applied load does not cause an instantaneous deformation,but there is a time lag between the application of load and the resulting deformation.This behavior is modelled mathematically using different approaches,the simplest among which,is to represent the constitutive relation as a combination of several spring and dashpot elements connected in series and parallel,each with different stiffness and viscosity respectively.This form of representation is explained in detail by Betz [24].

    Another approach is to relate the time-dependent stresses and strains using differential operators.To model material with frequency dependence a large number of higher-order time derivatives and parameters are required,making this model timeconsuming.This gave rise to the use of fractional derivative model for the frequency dependent properties of elastomers,which uses fewer parameters,reducing the complexity [25,26].Enelund et al.[27] studied the time domain response of viscoelastic structures governed by constitutive relations involving fractional calculus operators,and compared the results obtained with analytical solution.Fredette et al.[28] propose a novel spectral element approach to determine the dynamic stiffness of elastomeric isolators with frequency dependent damping characteristics.

    To simulate the linear viscoelastic behavior of elastomers,the complex modulus model has been widely used,as described in detail in Refs.[29,30].The frequency-dependent complex modulus of an elastomeric sample is obtained by conducting Dynamic Mechanical Analyzer (DMA) testing.Gil-Negrete et al.[31],presented the variation in shear modulus of several natural rubbers of different hardness with frequency and strain amplitude.It is generally observed that the storage modulus of elastomeric materials increases with increasing excitation frequency,and decreases with increasing dynamic strain amplitude[2].

    Additionally,the complex modulus of elastomeric materials is also dependent on the static preload.Nashif et al.[2] documented that with increasing preload,the storage modulus increases,while the loss factor decreases at any frequency of excitation.It thus becomes imperative to model the elastomeric properties as a function of both static and dynamic strains.The properties of elastomers under combined static and dynamic loads can be obtained experimentally,but requires a large number of tests due to the number of parameters involved.Kim et al.[32,33] proposed a viscoelastic constitutive equation for pre-stressed rubbers,derived through linearization of Simo's nonlinear viscoelastic constitutive model and reference configuration transformation [34].They further compared the results for different rubber specimens from their model with extensive experimental data and reported good prediction of elastomeric characteristics.Ahn &Kim [35] further provided a simple semi-empirical method to estimate the dynamic stiffness of pre-deformed elastomers using their static deformation characteristics.Beda et al.[3] later proposed a multiplicative non separable variables law for characterizing highly preloaded viscoelastic materials subjected to small-amplitude vibrations.This involved splitting the material property dependency on static stretch and time/frequency as though they are independent of each other,allowing the modelling of elastomeric materials with independent static and dynamic tests.They further validated their model with the experimental results obtained by Kim&Yuon[32],and reported a good match.This model representing the combined static and dynamic effects on an elastomer was proposed by Nashif et al.[2] and is considered in the current study,though it has not been presented for Yeoh material model in literature.Moreover,application of this method for obtaining accurate dynamic response of pre-stressed vibration isolators using Finite Element(FE)analysis and comparison with experimental results is the specific outcome of the present work.

    In the design of vibration isolators for aerospace vehicles,inertial acceleration is generated by the thrusting of the vehicle and random vibration excitation are considered applied over the preloaded systems.The present study is on vibration isolators used in aerospace applications subjected to dynamic loads in a preloaded state,which is hitherto not addressed in the referred literature.This paper presents an improved approach to generate a mathematical model for the elastomeric materials and a corresponding FE formulation to satisfactorily predict the random vibration response of typical isolators.

    Axisymmetric Finite Element Formulation is developed using Total Lagrangian approach to predict the static deformation of the isolator.The deformed state under static load is considered for further dynamic studies.The elastomeric constants are obtained from uniaxial tension tests and shear mode tests using Dynamic Mechanical Analyzer (DMA).The modified material properties of the elastomers to account for the pre-deformation and frequency dependence are used to compute the harmonic response and random response of two types of vibration isolators with different elastomer formulations.Quadratic quadrilateral element is developed and used effectively for the analysis of rocket motor solid propellant grains[36] whose linear version with pressure variable condensation is used for the present study.

    2.Material characterization

    The vibration isolators considered in the present study consist of parts made from steel and two types of elastomers,viz.VPM2-225A and VPM2-450B,where VPM2 is the nomenclature of isolators used in VSSC which is based on Polymeric(Silicone rubber)and Molded construction,whereas 225A and 450B represent the load rating applicable.These materials are selected depending on where it is being used: in general elastomers with good heat,weather and wear resistance with good damping characteristics are preferred.For the present case it is chosen with good damping characteristic to reduce the tramissibility to have good vibration control in the load transfer end.These molded vibration isolators are used for mounting of equipment in a launch vehicle.The properties of steel are taken from literature and elastomers are statically characterized to study the non-linear static deformation.To numerically simulate these elastomers under dynamic loads,the elastomeric material samples are tested in a DMA for different dynamic strain amplitudes and the properties are obtained as complex frequency dependent functions over a range of frequencies.The properties obtained from DMA are modified to account for the preload effect by including the static stretch and non-linear material constants.The procedure used to obtain these properties is provided in this section.

    2.1.Static material characterization

    The behaviour of hyper-elastic materials is derived from the strain energy density function,U,which is a function of the three invariants,I1,I2andI3of the Green deformation tensor.In the current study first,the suitability of Mooney-Rivlin and Yeoh material models are initially considered for comparison.The constants for these models are obtained by performing a uniaxial tension test.

    A uniaxial tension test was conducted on dumbbell shaped specimens prepared according to ASTM D412 standards.The dimensions of the specimen used for the test are shown in Fig.1.It is generally accepted that tensile test data is sufficient to compute Yeoh constants and this can represent deformation in other modes as well.

    Fig.1.Uniaxial tension test standard specimen dimensions.

    The test was conducted at ambient conditions for both the elastomeric material samples of VPM2-225A and VPM2-450B,and the measured stress vs.strain plots are presented.The hyper-elastic curve fitting capability in ABAQUS software is used to compute the hyper-elastic material constants for Mooney-Rivlin and Yeoh material models by providing the nominal stress vs.nominal strain data obtained from the tension test.By iterating through the constants using a least-squares fit method,the relative error between the predicted and experimental values is reduced by ABAQUS[37].

    The material constants calculated along with the corresponding goodness of fit values are presented in Table 1.The Fit obtained with both the material models are shown in Fig.2 and Fig.3 for materials VPM2-225A and VPM2-450B respectively.Brief description of the procedure to obtain Mooney-Rivlin and Yeoh material constants provided in the table below.These are from a linear fit between Cauchy's stress vs principal stretch and quadratic fit between Cauchy's stress vs(I1-3).Cauchy's stress can be obtained from strain energy density function for uniaxial tension which is provided in Section 3 for Yeoh material model.For Mooney-Rivlin model,strain energy density function is given as

    Table 1 Material constants obtained from ABAQUS.

    Fig.2.Moony Rivlin &Yeoh models fit with Tension Test data for VPM2-225A.

    Fig.3.Moony Rivlin &Yeoh models fit with Tension Test data for VPM2-450B.

    From Fig.2 and Fig.3 it is clearly evident that Yeoh material model fit much better compared to Mooney-Rivlin model.Hence Yeoh material constants will be used for numerical study.

    It can be seen from Fig.2 and Fig.3 that the Yeoh material model matches well with the observed non-linear behavior very closely.It can also be noted that the Mooney-Rivlin model fails to predict the“upturn” or stiffening at higher strains.Yeoh material model is further used for model improvement and response studies as the baseline.

    2.2.Dynamic material characterization

    The dynamic characterization for the elastomers under study is done using the TA Instruments Q800 Dynamic Mechanical Analyzer(DMA).The force resolution of the DMA is 0.0001 N with a maximum force of 18 N,while the strain resolution is 1 nm.The measurement range of modulus for the Q800 DMA is 10E3 to 30E12 Pa with a precision of ±1% [38].

    The Q800 DMA can be used to conduct dynamic property evaluation in different modes such as single/dual cantilever,thin film tension,compression etc.For the current study,the dynamic characterization is done in the shear mode by using the TA instruments shear sandwich clamp.This mode of test is preferred,because the elastomer in the isolators under study experience predominately shear deformation due to their geometric configuration.The test is conducted at room temperature for different strain amplitudes(Υd)with the frequency ranging from 1 to 60 Hz.The frequency dependent storage modulus (shear) and loss factor are measured at different dynamic strain amplitudes and the results are presented in Fig.4 and Fig.5.

    Fig.4.Frequency dependent properties of VPM2-225A.

    Fig.5.Frequency dependent properties of VPM2-450B.

    It can be seen that at a constant dynamic strain amplitude,the storage modulus increases with excitation frequency,while the loss factor remains almost constant.At constant frequency however,the storage modulus reduces with increasing strain amplitude while the loss factor increases.For VPM2-225A,the testing is done till Υdof 6.6%,and for Υdabove this value,the material sample slips off the clamps.

    2.3.Geometry and material properties

    The 3-D geometric configuration of type-A and type-B vibration isolators used for present analysis,are shown in Fig.6 and Fig.7 respectively.The isolator consists mainly of four parts,the pin,body,elastomer and cover.The pin,body and cover are made of AISI 304 steel.

    Fig.6.Configuration of Type-A (VPM2-225A) isolator.

    Fig.7.Configuration of Type-B (VPM2-450B) isolator.

    Fig.8.4-Node Quad element.

    The material properties considered for the analysis are given in Table 2.

    Table 2 Material properties used for non-linear static analysis.

    3.Material model formulation for incorporating staticstretch dependence

    In this section,a novel multiplicative non-separable variables law is derived,accounting for the static deformation and its effect on dynamic material characteristics,for Yeoh material model.The method used is similar to that given in Refs.[2,3].The assumption involved is that the effects of static stretch and frequency of excitation on the dynamic material properties are independent of each other.For elastomeric materials subjected to a dynamic load superimposed on a non-linear static deformation,the stress can be factored into a function of frequency (ω) and a function of stretch ratio (λ) as

    Strain energy density function as given by the standard Yeoh material model,

    where,I1is the first invariant of the green deformation tensor,

    where,λ1,λ2,λ3are the principal stretches in 3 directions.For uniaxial tension,

    From Eq.(4),Eq.(3) becomes

    The stress is given as

    Simplifying,

    Taking the right-hand side of Eq.(7) asF(λ)in Eq.(1),

    Eq.(9) can be used to represent the combined static and dynamic properties given after simplification,

    Using the complex modulus representation of the elastomeric material properties,

    Separating the real and imaginary parts of the above equation,real part is provided in Eq.(10) and imaginary part is given as

    To determine η1,η2,η3and H,it is to be noted that from Ref.[27],η2and η3can be assumed to be zero.As λ→1,

    From Eq.(10) to Eq.(16),

    Thus,Eq.(10) and Eq.(15) become

    whereY(ω)and η(ω)are the frequency dependent storage modulus and loss factor obtained from DMA testing.Using Eq.(18) and Eq.(19),the material properties are modified for at each frequency value while computing the response.Eq.(18)is modified with shear modulus and bulk modulus using elasticity relations with the material constants obtained from DMA tests for numerical computation.It is evident that with increase in static stretch (increase in load) the storage modulus increases and the damping reduces.Especially at higher static preloads,the modulus increases significantly.In the present study,significant effects on elastomeric constants are observed due to pre-load so as to change the dynamic characteristics of the isolator.

    4.Finite element (FE) formulation for response studies

    The axi-symmetric four node quadrilateral element for geometric non-linear analysis using Total Lagrangian approach is provided along with Newton-Raphson solution technique.The statically deformed configuration is updated to study frequency domain response.For this purpose,axi-symmetric four node element is developed using mixed u/p formulation with condensation of pressure degrees of freedom.Two by two Guassian quadrature is used for numerical computation of load vector and finite element matrices.Throughout this section,formulation is provided for one element and load vector and other assembly and solution procedure follows the usual finite element methods.

    4.1.Non-linear axisymmetric FE formulation for hookean material

    The formulation for 3D finite elements is given by Zienkiewicz[39],and is suitably modified for four-node axisymmetric finite elements,as presented in this section.The residual force for axisymmetric four-node quadrilateral element is given as

    whereArepresents area of the element,is the incremental strain displacement matrix and can be obtained from Green-Lagrange strain vector Ev,for axi-symmetric body.Svis the second P-K stress vector and Feis the externally applied load vector which is deformation independent.For the present study,only traction is applied,the details of the load vector for the traction are given as

    where l represents the line of the edge where pressurepris applied andm,nare the direction cosines of the edge and N is the shape function matrix.Relation between second P-K stress vector and Green-Lagrange strain vector is given as

    where D is the stress-strain relation.Axisymmetric quadrilateral four-node element is shown in Fig.8.

    The non-linear Eq.(20) is linearized through Taylor's series expansion considering only the linear term as provided in Zienkiewicz [39].Equilibrium solution which satisfies Eq.(20) is obtained in an iterative process using the linearized equation starting with zero displacement.To obtain the iterative solution,tangent stiffness matrix is computed as given below having relevance for Yeoh material model provided in the next section.

    where KT(q),G,MSare tangent,geometric stiffness,stress matrices respectively and dψ,dq are incremental residual force and displacement vectors respectively.

    4.2.Axisymmetric FE formulation for yeoh material model

    The relations for second P-K stress vector and constitutive matrix for axisymmetric mixed p-type element is given by Muhammed et al.[1],and is modified for Yeoh material model as presented in this section.The residual force vector and tangent stiffness matrix provided can be directly used for Yeoh material model with proper definition of second P-K stress vector and the constitutive relation.For that purpose,strain energy density function for Yeoh material model is used.The first and the second partial derivative of strain energy density function with respect to Green Lagrange strain provides second P-K stress vector and the constitutive matrix respectively.Strain energy density function for Yeoh material model is given as

    where,C10,C20,C30are Yeoh material constants,D1,D2,D3are penalty parameter to impose the incompressibility constraints,

    For axi-symmetric element,

    Second P-K stress tensor is defined as

    Green-Lagrange strain tensor is defined as

    where I is a 3 × 3 Identity matrix,F is the deformation gradient given as

    Second P-K stress vector is given as

    where,

    The constitutive matrix for Yeoh material is given as

    Using the above formulation non-linear static analysis is carried out.The details of loading and displacement boundary conditions are provided in the numerical study.The statically deformed state is used to carry out further the frequency domain response analysis.

    4.3.Axisymmetric FE formulation for frequency domain response

    Static mixed u-p based formulation for quadratic elements to address incompressibility effects is given by Marimuthu et al.[36],whose linear element version is used in the present study with pressure condensation.In this section,four-node quadrilateral finite element formulation is presented for incompressible material for static load.The same formulation is used for compressible material as well,and necessary modifications are made to study the frequency domain response.The total potential for an axisymmetric solid by imposing incompressibility constraint is given as

    where ε,σ,εvare the linear strain,stress vectors and volumetric strain respectively,p is the hydrostatic pressure at the element centroid,k is the bulk modulus,b is the body force,t is the traction,u is the displacement vector at any point within the element,and q is the nodal displacement vector.Minimization of total potentialw.r.to q,p after substituting strain vector as a sum of deviatoric and volumetric parts and stress vector as a sum of deviatoric and hydrostatic part results in the following with pressure condensation given as

    where Gdis a diagonal matrix given as{2G,2G,2G,G}which relates deviatoric stress and strain vectors,Bdis the deviatoric strain displacement matrix and Bvis a vector which relates voulmentic strain and displacement vector.

    4.3.1.FE formulation for frequency response analysis

    The static finite element equation for an incompressible material can be modified and written for dynamics by using the d’Alembert's principle [39].In dynamic analysis,the displacement vector at any point within the element u and nodal displacement vector q be complex and henceforth will be represented with a × superscript.When displacements of an elastic body vary with time,two sets of additional forces are called into effect.The first is the inertia force,which for an acceleration characterized by ü*can be replaced by its static equivalent,-ρü*where ρ is the mass density.The second force is that due to resistance opposing motion.A linear,viscous type resistance can be characterized by an equivalent static force of magnitude -c˙u*,wherecis damping property per unit volume.The equivalent static problem,at every instant of time,can be discretized in a manner similar to previous section,by replacing the distributed body force b by its equivalent

    In the present study,damping is expressed in the material property itself,hence viscous damping is neglected.Assuming

    Static finite element formulation provided in Eq.(43) get modified for frequency response as

    The above matrix is a complex function of frequency and static stretch dependent is due to the material properties.The diagonal matrix Gdin equation and bulk modulus in which was provided in static formulation get modified as={2G*(ω,λ),2 G*(ω,λ),2 G*(ω,λ),G*(ω,λ)} andk*(ω,λ).

    whereG(ω,λ)andk(ω,λ)are storage shear and bulk modulus respectively and η(ω,λ)is the loss factor and static stretch λ is a geometric mean of the eigen values of matrix CM.Mass matrix M is given as

    4.3.2.FE formulation for random response analysis

    A random process is said to be stationary when its mean does not change over time.This essentially means that a stationary process is one whose probability distribution does not change with time.In addition,if the mean obtained for a sample of an ensemble,is equal to the mean across the ensemble at any arbitrary time ti,the process is called ergodic[40].For the numerical analysis presented in this study,all random processes are assumed to be stationary ergodic.The finite element formulation for obtaining the random response is an extension of the formulation presented in the previous section for obtaining the frequency response.The RMS value of the reactive response in the Z direction to a random excitation with a spectral densitySF(fn) can be obtained as

    wherefnis the excitation frequency in Hertz,Nis the number of equal intervals of excitation frequency,andSRZ(fn)is the spectral density of the response given as

    where,|Hn| is the absolute value of the net reactive frequency response at frequency inZdirection atfn.

    4.4.Implementation algorithm

    The procedure implemented for obtaining the frequency response of the isolators is described in Fig.9 is as follows:

    Fig.9.Implementation algorithm for frequency response analysis.

    The non-linear static analysis of the isolators shown in Fig.6 and Fig.7 is carried out under static loading conditions at the rated load and the deformed configuration thus obtained is considered for the dynamic analyses.For the present study,the dynamic material properties of the elastomers have been obtained at different dynamic strain amplitude(Υd)values as a function of frequency.These properties are modified to include the effect of static stretch as described in section 3.

    The frequency response analysis is first carried out using the material properties at lowest Υdand the resonant frequency match between the experiment and numerical result is checked.If a match is not observed,the material properties at a higher Υdare used for obtaining the response.This process is continued until a close match is observed in the resonant frequency obtained from experiment and simulations.Further,a suitable material damping correction is provided to match the numerical transmissibility to the experiments.

    5.Experimental setup for vibration response studies

    The experimental study of isolators for harmonic response and random response is conducted using the electrodynamic shaker consists of an electromagnetic circuit with a stationary (field) coil and a moving (drive/armature) coil which is a part of the “Head”that vibrates.The force rating of the shaker is 3.5 ton with 50 mm double amplitude capability having frequency range is up to 2000 Hz shown Fig.10.The isolator used for the present study is fixed on top of the moving head.To simulate the load rating an equivalent to mass of 2.25 kg and 4.5 kg for Type-A and Type-B isolators,a dummy mass is made in hollow cylindrical form with bottom closed whose inverted position is placed on top of the isolator shown in Fig.10.Due to the placement of mass on the pin an equivalent load is transferred to the pin and futher a precompression is created due to specified displacement on the cover.This pre-stressed state is obtained using non-linear analysis presented in numerical study section.This pre-compressed state is used for the dynamic experimentation by fixing an accelerometer on top of dummy mass as shown in Fig.10 which is close to the flexural centre to avoid unwanted noise pick up.The experimental observations are compared in a graphical form with the present FE formulation are detailed in subsection 6.3.

    Fig.10.Dummy mass placed on top of the isolator for test with accelerometer attached.

    The data acquisition system used is the Spider-81 supplied by Crystal Electronics.The Integrated Electronics Piezoelectric (IEPE)accelerometer used is of Endevco make with a sensitivity of 10 mV/g.Data analysis is performed using the Engineering Data Management(EDM)software,compatible with Spider-81.The above test is conducted at Saraswathi Dynamics at the Vibration Test Facility in Liquid Propulsion Systems Centre(LPSC),ISRO,Trivandrum.

    6.Numerical study

    The results from static and frequency domain analyses are presented in this section.The appropriate material properties are taken from Table 2.The numerical results are presented in contour and graphical form.

    6.1.Finite element mesh

    The axisymmetric finite element mesh with the applied boundary conditions for type-A and type-B isolators are shown in Fig.11 and Fig.12 respectively.The nodes with red circles are fully fixed,whereas the nodes with blue circles are constrained in the radial direction.The arrows represent the location and direction of the applied load.The finite element mesh connectivity and nodal coordinates are generated using ‘PreWin’,a pre and post processing software which is developed in ISRO for the Finite Element Program“FEAST”.The data generated is then exported to a MATLAB code for further computations.The elements colored green and white in Fig.11 and Fig.12 correspond to the elastomer and the steel parts respectively.The FE details of the models used are provided in Table 3.

    Table 3 Finite Element details.

    Fig.11.Finite Element mesh of Type -A isolator.

    Fig.12.Finite element mesh of Type -B isolator.

    In Figs.11 and 12,the interface between the pin and elastomer are assumed as bonded and modelled using common nodes and this assumption is valid due to high friction coefficient for elastomers.

    6.2.Non-linear static analysis

    For the isolators under study,an initial pre-compression isprovided when the isolator body and cover are riveted together.A mass equal to the rated mass of the isolator is also placed on top during its operation and testing.For both the isolators,the precompression provided is 1.3 mm.For type-A isolator which has a mass rating of 2.25 kg,the static load applied is 22.5 N as edge load on the pin.For type-B isolator with a mass rating of 4.5 kg,the static load applied is 45 N.The results of the static analysis carried out using the code are presented in Table 4.The response is given in terms of displacement of the pin.The resulting deformation contour is given in Fig.13 and Fig.14 for type-A and type-B isolators respectively.

    Table 4 Pin displacements obtained from non-linear static analysis.

    Fig.13.Deformation contour of Type -A isolator.

    Fig.14.Deformation contour of Type -B isolator.

    6.3.Frequency domain analysis

    The statically deformed configuration is considered for the numerical and experimental studies with the addition of point mass equivalent to dummy mass at the top most node in axis of the pin.Five and three experiments respectively are conducted for frequency and random response which are bought out in tabular form provided in Table 5 and Table 6 for Type-A and Type-B isolator.The difference in numerical values from one experiment to another is attributed to variation in material properties.

    Table 5 Frequency response data obtained from 5 trials with Type -A,B isolators.

    Table 6 Response RMS obtained from random excitation of isolators.

    6.3.1.Experimental frequency response

    Five experiments are conducted on different isolators whoseresonance frequency and amplitude are given in Table 5.

    6.3.2.Comparison of frequency response results

    For the numerical comparison one of the experiments is chosen which is trial-1 from Table 5.Numerically observed frequency response is a net reactive response computed from eliminated system matrix,with multiplication of displacement vector computed for frequency ranging from 10 to 60 Hz in the interval of one Hz.Fig.15 shows the comparison between experimental observation and numerically computed values as the properties obtained for different Υdvalues as shown in Fig.4.

    Fig.15.Comparison of frequency response with material properties at varying dynamic strain amplitudes.

    It can be seen from Fig.15 that resonance frequency and amplitude are high compared to the experimental value for lower strain amplitude due to the dynamic elastomeric properties such as storage modulus and loss factor which are progressively lower and higher respectively as the frequency increases.This phenomenon is well presented in Fig.4 and explanation are provided thereafter.From the numerical comparisons shown in Fig.15 has a decent match in natural frequency is achieved for 6.6%Υd,while the peak transmissibility is only 3.2.The resonance frequency is indicated against each response curve.Here,a damping correction factor is introduced into the mathematical model to correct and match the numerically obtained transmissibility with that of the experiment.Since the random response computations are dependent on the response at each frequency component,applying this correction factor is essential for accurate response predictions.This will be evident with the results presented in the subsequent section.Transmissibility shown in Fig.15 is a magnification factor for the input provided on the isolator.Elastomers are designed such a way that this factor should be the least to ensure minimum response for sinusoidal excitation.For random excitation lower resonance frequency with lower transmissibility helps in overall responsereduction for acoustic excitations,stage separation etc.,to protect mission critical electronic packages mounted in launch vehicles.

    An important aspect to note here regarding hysteretically damped systems is that their natural frequency is not affected by the damping present in the system but rather only by the stiffness.A change in damping will only affect the transmissibility of the system as shown in Fig.16.

    Fig.16.Experimental vs.numerical frequency response with updated damping.

    The response shown in Fig.16 is obtained by scaling down the damping to 50% of the value obtained from DMA testing.The variation in natural frequency between FE analysis and experiment is 4.3% while the variation in peak transmissibility after scaling down the damping is 1.9%.The same scaling factor for damping is applied for random response study of Type-A isolator.The experimental frequency response of the isolator shows a small peak at around 860 Hz.Resonance at such a high frequency suggests that it might be due to the highly stiff components of the isolator such as the steel base,cover and pin.A sharp dip is seen in the numerical frequency response as shown in Fig.17.The cut-off frequency is under-predicted numerically at 51.5 Hz against that observed experimentally at 56 Hz.

    Fig.17.Frequency Response of Type -A Isolator: Numerical vs.Experiment.

    Fig.18 shows the real and imaginary parts of the response and it can be seen that the real part of the response changes sign and goes to zero at resonance and the response is affected only by the imaginary part i.e.,the damping associated response.

    Fig.18.Real and imaginary parts of numerical frequency response.

    In the present mathematical model,the damping properties of steel are not considered,hence both the real and imaginary parts of response go to zero at the second resonance point,as shown in Fig.19.

    Fig.19.Real and imaginary parts of numerical frequency response at 840 Hz.

    This leads to the response at this frequency being close to zero.In spite of this,as can be seen from Fig.17,that the response at such high frequencies is very low and this mismatch in transmissibility is considered acceptable here.This behaviour of the numerical result is further noticed in the random response studies as well.

    The material properties for the compound VPM2-450B are obtained up to 5% Υd.It is noticed that there is good match in the natural frequency with properties at 1%Υdand the higher dynamic strain data is not used.The numerical frequency response is shown in Fig.20.

    Fig.20.Numerical vs.experimental frequency response of type -B isolator.

    The correction factor for damping applied in this case is 1.25.The variation in natural frequency is0.92%and that in peak transmissibility is4.5%.The numerical response is compared with that from the experiment and is presented in Fig.21.The cut-off frequency obtained numerically is 44.5 Hz,while the experimentally observed value being 47 Hz.It can be inferred that the cut-off frequency is under-predicted with the current mathematical model.

    Fig.21.Frequency response of type -B isolator: numerical vs.experiment.

    6.3.3.Experimental random response

    Experiments were done with three different isolators of Type-A and Type-B with excitation having same power spectral density and an RMS value of 13.5 g.The response RMS is listed in Table 6.6.3.4.Comparison of random response results

    The power spectral density of the excitation,with an RMS value of 13.5 g,and the corresponding response are shown in Fig.22 and Fig.23 in comparison with the numerical results,where,‘Control Input’is the input PSD.The FE model used for obtaining the random response is the same as described in frequency response.

    Fig.22.Random Response of Type -A Isolator: numerical vs.experiment.The response RMS obtained numerically is 0.8538 g with an error of 2.7% with respect to the experimental mean,which is 0.8629 g.The comparison of numerical and experimental random response for type-B isolator is shown in Fig.23.

    Fig.23.Random Response of Type -A Isolator: numerical vs.experiment.

    The response RMS obtained numerically is0.673with an error of 14.5%with respect to the experimental mean,which is0.7872 g.

    7.Conclusions

    An improved material model for elastomers and a finite element formulation for frequency domain analysis of the vibration isolators is presented.The dynamic analyses are carried out using a mathematical model based on Yeoh hyper-elastic material model modified to handle preload and frequency dependence.The material data obtained from uniaxial tension test and DMA tests are used for the mathematical model generation.The frequency response of both type-A and type-B isolators are computed numerically and with material properties obtained at all Υdvalues.It is observed that for type-A isolator,the properties obtained at 6.6% Υdgive a close match in natural frequency of the isolator obtained experimentally.For type-B isolator,a match in natural frequency is obtained with material properties at 1% ΥdHowever,in both cases the error in the computed peak transmissibility is high,and hence a suitable damping correction factor is applied to the loss factor measured with the DMA,to match the numerically computed transmissibility with the experiments.This is done in order to compute the random response of the isolators with good accuracy.With the corrected loss factor,the random response of the isolators is computed numerically and compared with experiments.A good match in the response RMS between the numerical and experimental values is observed.Frequency response and random response matches well between 10 and 60 Hz after a suitable damping correction.The reason for this match can be attributed to the dynamic elastomeric constants that are obtained only up 60 Hz,but for higher frequencies the same material constants are being used and the deviation in the numerical simulations are well observed in the response graph.The system matrix is close to singularity at higher frequencies,which is due to assumption of zero damping of steel parts.

    The approach presented in this paper can be used for predicting the vibration response of equipment mounted using elastomeric vibration isolators in aerospace vehicles,where random excitations are generally encountered on statically loaded systems.Similar such approach can be suitably adopted for other FE model depends on as is where is basis.

    Declaration of competing interest

    The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

    Acknowledgements

    The authors like to gratefully acknowledge the contributions of Shri S Sai Sandeep Kumar,Dr RS Rajeev and Shri S Chandramouli of VSSC,ISRO.The authors would like also to thank Dr S Unnikrishnan Nair,Director,VSSC,ISRO for his unwavering support in carrying out this work.

    99国产极品粉嫩在线观看| 成年人黄色毛片网站| 午夜久久久在线观看| 丝袜美足系列| 日韩免费av在线播放| 80岁老熟妇乱子伦牲交| 99国产精品一区二区蜜桃av | 精品久久久久久电影网| 午夜日韩欧美国产| 国产欧美日韩一区二区三区在线| 欧美日本中文国产一区发布| 最黄视频免费看| 亚洲av美国av| 亚洲一区二区三区欧美精品| 高清在线国产一区| 日本黄色视频三级网站网址 | 久久天躁狠狠躁夜夜2o2o| av又黄又爽大尺度在线免费看| 中文字幕av电影在线播放| 在线观看www视频免费| 日韩欧美一区二区三区在线观看 | 国产精品 国内视频| 精品欧美一区二区三区在线| 一区二区三区乱码不卡18| 91字幕亚洲| 亚洲精品美女久久久久99蜜臀| 狠狠狠狠99中文字幕| 日韩大片免费观看网站| 国产在视频线精品| 成人国产一区最新在线观看| 亚洲av日韩精品久久久久久密| 亚洲熟女精品中文字幕| 久久久久久久久免费视频了| 国产成人欧美| 欧美日韩中文字幕国产精品一区二区三区 | 国产欧美日韩一区二区精品| 一级毛片精品| 亚洲 国产 在线| 在线天堂中文资源库| 成人手机av| av天堂久久9| 在线观看66精品国产| 精品一区二区三卡| 最近最新免费中文字幕在线| 成人18禁在线播放| 久久国产精品影院| 高清视频免费观看一区二区| 日本一区二区免费在线视频| 一区二区三区乱码不卡18| 亚洲少妇的诱惑av| 狠狠婷婷综合久久久久久88av| 国产区一区二久久| 亚洲午夜精品一区,二区,三区| 巨乳人妻的诱惑在线观看| 国产亚洲精品久久久久5区| 国产97色在线日韩免费| 免费日韩欧美在线观看| 久久精品国产亚洲av香蕉五月 | 欧美国产精品一级二级三级| 精品少妇一区二区三区视频日本电影| 美女高潮喷水抽搐中文字幕| 一进一出抽搐动态| 丰满少妇做爰视频| 日韩视频在线欧美| 69av精品久久久久久 | 亚洲精品美女久久久久99蜜臀| 后天国语完整版免费观看| 国产精品一区二区免费欧美| 久久久久久人人人人人| 黄色视频,在线免费观看| 菩萨蛮人人尽说江南好唐韦庄| 日日夜夜操网爽| 中文字幕制服av| 国产在线视频一区二区| 交换朋友夫妻互换小说| 亚洲国产中文字幕在线视频| 国产av国产精品国产| 国产精品久久久久久精品电影小说| 亚洲第一av免费看| 99riav亚洲国产免费| 黑人欧美特级aaaaaa片| 亚洲 国产 在线| 亚洲人成77777在线视频| 极品人妻少妇av视频| 亚洲国产欧美网| 男女午夜视频在线观看| 欧美精品人与动牲交sv欧美| 91大片在线观看| 成人特级黄色片久久久久久久 | 多毛熟女@视频| 少妇 在线观看| 国产一卡二卡三卡精品| 热99re8久久精品国产| 亚洲国产欧美网| 久久久国产欧美日韩av| 男女边摸边吃奶| 亚洲精品国产色婷婷电影| 黄网站色视频无遮挡免费观看| 人妻一区二区av| 免费观看a级毛片全部| 一级毛片精品| 日韩视频在线欧美| 男女边摸边吃奶| 两个人看的免费小视频| 黑人巨大精品欧美一区二区蜜桃| 国产有黄有色有爽视频| www.熟女人妻精品国产| 超碰97精品在线观看| 捣出白浆h1v1| 性少妇av在线| av免费在线观看网站| 黄片小视频在线播放| 欧美成狂野欧美在线观看| 黄色毛片三级朝国网站| 美女午夜性视频免费| 啪啪无遮挡十八禁网站| 亚洲精品av麻豆狂野| 午夜免费成人在线视频| 久久天躁狠狠躁夜夜2o2o| 久久久久久久久免费视频了| 波多野结衣av一区二区av| 99国产精品一区二区蜜桃av | 天堂俺去俺来也www色官网| 老汉色av国产亚洲站长工具| 三级毛片av免费| 日韩免费高清中文字幕av| 亚洲国产看品久久| 一区在线观看完整版| 美女福利国产在线| 男女午夜视频在线观看| 又大又爽又粗| 久久精品国产综合久久久| 国产在线观看jvid| 1024视频免费在线观看| 色综合婷婷激情| tube8黄色片| 一级,二级,三级黄色视频| 欧美久久黑人一区二区| 久久久久久久大尺度免费视频| 黄色视频,在线免费观看| 国产亚洲欧美精品永久| 亚洲国产欧美日韩在线播放| 男女免费视频国产| xxxhd国产人妻xxx| 亚洲专区国产一区二区| 美女高潮喷水抽搐中文字幕| 99久久国产精品久久久| 热99国产精品久久久久久7| 午夜福利一区二区在线看| 久久午夜亚洲精品久久| 免费日韩欧美在线观看| 日韩三级视频一区二区三区| 性少妇av在线| 成人特级黄色片久久久久久久 | 淫妇啪啪啪对白视频| 亚洲国产成人一精品久久久| 国产精品久久久久久人妻精品电影 | 亚洲成人免费av在线播放| 高清黄色对白视频在线免费看| 国产精品一区二区在线不卡| 国产亚洲精品第一综合不卡| 成人手机av| 欧美日韩国产mv在线观看视频| 日韩视频在线欧美| bbb黄色大片| 久久精品亚洲熟妇少妇任你| 欧美老熟妇乱子伦牲交| 欧美黄色淫秽网站| 欧美日韩黄片免| 国产av精品麻豆| 夫妻午夜视频| 亚洲全国av大片| 国产成人精品久久二区二区91| 国产精品一区二区免费欧美| 男女午夜视频在线观看| 成人黄色视频免费在线看| 精品国产国语对白av| 精品人妻在线不人妻| 波多野结衣一区麻豆| 欧美乱码精品一区二区三区| 91精品三级在线观看| 桃花免费在线播放| 亚洲五月色婷婷综合| 视频区欧美日本亚洲| 色综合欧美亚洲国产小说| 久久久久久免费高清国产稀缺| 人人妻人人澡人人爽人人夜夜| 午夜久久久在线观看| 久久久久久久大尺度免费视频| 亚洲欧美一区二区三区黑人| 真人做人爱边吃奶动态| a级片在线免费高清观看视频| 99久久国产精品久久久| 成人特级黄色片久久久久久久 | 久久99热这里只频精品6学生| 天天躁狠狠躁夜夜躁狠狠躁| 欧美日韩福利视频一区二区| 日韩一区二区三区影片| 18禁美女被吸乳视频| 天堂8中文在线网| 精品亚洲成a人片在线观看| 亚洲精品自拍成人| 我要看黄色一级片免费的| a级片在线免费高清观看视频| 亚洲人成伊人成综合网2020| 久久久久久免费高清国产稀缺| 中文字幕另类日韩欧美亚洲嫩草| 午夜免费鲁丝| av欧美777| 操出白浆在线播放| 亚洲欧美激情在线| 在线看a的网站| 日本av手机在线免费观看| 欧美精品高潮呻吟av久久| 嫁个100分男人电影在线观看| 亚洲国产欧美网| 国产不卡av网站在线观看| 亚洲伊人久久精品综合| 飞空精品影院首页| 欧美精品人与动牲交sv欧美| 无限看片的www在线观看| 美女高潮喷水抽搐中文字幕| 国产精品欧美亚洲77777| 18禁观看日本| 久久久久国内视频| 欧美性长视频在线观看| 老司机影院毛片| 最近最新中文字幕大全电影3 | 亚洲,欧美精品.| 国产区一区二久久| 久久久久久久国产电影| 桃红色精品国产亚洲av| 下体分泌物呈黄色| 久久人妻福利社区极品人妻图片| 看免费av毛片| 在线观看免费视频网站a站| 国产欧美日韩综合在线一区二区| 精品一区二区三区视频在线观看免费 | 757午夜福利合集在线观看| 在线天堂中文资源库| 亚洲欧洲日产国产| 国产野战对白在线观看| 精品国产乱码久久久久久小说| 黄色视频在线播放观看不卡| 国产亚洲精品久久久久5区| 一级毛片女人18水好多| 国产不卡av网站在线观看| 18禁裸乳无遮挡动漫免费视频| 国产精品久久久久久精品古装| 高清黄色对白视频在线免费看| 啦啦啦视频在线资源免费观看| 国产免费av片在线观看野外av| 亚洲成a人片在线一区二区| 动漫黄色视频在线观看| 色精品久久人妻99蜜桃| 成年人午夜在线观看视频| 性高湖久久久久久久久免费观看| 青青草视频在线视频观看| 免费少妇av软件| 在线观看免费午夜福利视频| 国产精品久久久久久精品古装| cao死你这个sao货| 又黄又粗又硬又大视频| 精品久久久久久久毛片微露脸| 国产日韩一区二区三区精品不卡| 精品人妻在线不人妻| 国产成人免费观看mmmm| 亚洲av成人不卡在线观看播放网| 69精品国产乱码久久久| 天堂俺去俺来也www色官网| 在线看a的网站| 午夜日韩欧美国产| 19禁男女啪啪无遮挡网站| 在线播放国产精品三级| 国产精品香港三级国产av潘金莲| 18禁裸乳无遮挡动漫免费视频| 中文字幕制服av| 这个男人来自地球电影免费观看| 亚洲色图av天堂| 国内毛片毛片毛片毛片毛片| 亚洲中文日韩欧美视频| 女人久久www免费人成看片| 一级毛片电影观看| 啦啦啦免费观看视频1| 三级毛片av免费| 一本大道久久a久久精品| 日韩精品免费视频一区二区三区| 精品欧美一区二区三区在线| 热re99久久精品国产66热6| 亚洲精品美女久久av网站| 黑人猛操日本美女一级片| 日韩成人在线观看一区二区三区| 三级毛片av免费| 亚洲精品中文字幕在线视频| 两性夫妻黄色片| 国产精品一区二区精品视频观看| 99国产精品一区二区蜜桃av | 国产在线观看jvid| 成人国产av品久久久| 久久久国产成人免费| a级毛片黄视频| 久久影院123| videosex国产| 在线av久久热| 婷婷成人精品国产| 欧美在线一区亚洲| 亚洲成人手机| 精品一品国产午夜福利视频| 午夜精品久久久久久毛片777| 国产成+人综合+亚洲专区| 亚洲男人天堂网一区| 国产国语露脸激情在线看| 久久久久网色| 欧美日本中文国产一区发布| 国产精品欧美亚洲77777| 国产野战对白在线观看| 两性夫妻黄色片| 一边摸一边抽搐一进一出视频| 亚洲精品粉嫩美女一区| 久久久精品94久久精品| 亚洲欧美精品综合一区二区三区| 国产成人精品久久二区二区91| 亚洲专区字幕在线| 中文字幕另类日韩欧美亚洲嫩草| 成人永久免费在线观看视频 | 国产精品电影一区二区三区 | 香蕉国产在线看| 国产精品亚洲一级av第二区| 超色免费av| 侵犯人妻中文字幕一二三四区| 狠狠狠狠99中文字幕| 成人永久免费在线观看视频 | 少妇 在线观看| 欧美精品一区二区大全| 国产精品久久久久久精品电影小说| 免费观看人在逋| 咕卡用的链子| 叶爱在线成人免费视频播放| 欧美人与性动交α欧美软件| 一夜夜www| 欧美乱码精品一区二区三区| 日韩欧美免费精品| 热99国产精品久久久久久7| 日韩欧美免费精品| 国产精品熟女久久久久浪| 91av网站免费观看| 免费日韩欧美在线观看| 久久久久久久精品吃奶| 亚洲色图av天堂| 天天操日日干夜夜撸| 欧美性长视频在线观看| 欧美在线一区亚洲| 亚洲精品一卡2卡三卡4卡5卡| 久久亚洲真实| 老熟妇乱子伦视频在线观看| 色精品久久人妻99蜜桃| 亚洲熟女毛片儿| 视频在线观看一区二区三区| 日韩中文字幕欧美一区二区| 欧美日韩视频精品一区| 黄色 视频免费看| 成人黄色视频免费在线看| 蜜桃国产av成人99| 亚洲人成电影免费在线| 久久精品国产a三级三级三级| 免费在线观看黄色视频的| 欧美日韩亚洲高清精品| 久久久精品国产亚洲av高清涩受| 91成人精品电影| 日本撒尿小便嘘嘘汇集6| www.自偷自拍.com| 日本vs欧美在线观看视频| 日韩欧美国产一区二区入口| 在线av久久热| 精品视频人人做人人爽| 国产片内射在线| 久久久欧美国产精品| 国产精品久久电影中文字幕 | 性高湖久久久久久久久免费观看| 精品国产国语对白av| 成在线人永久免费视频| 超碰成人久久| 男女下面插进去视频免费观看| 国产精品久久久久久精品电影小说| 人成视频在线观看免费观看| 一级a爱视频在线免费观看| 一边摸一边抽搐一进一小说 | 热99久久久久精品小说推荐| 人人澡人人妻人| 亚洲男人天堂网一区| 高潮久久久久久久久久久不卡| 精品人妻熟女毛片av久久网站| 精品一区二区三区av网在线观看 | 男女高潮啪啪啪动态图| 午夜福利,免费看| 我的亚洲天堂| 欧美日韩中文字幕国产精品一区二区三区 | 在线亚洲精品国产二区图片欧美| 老司机午夜十八禁免费视频| 亚洲成人免费av在线播放| 日韩中文字幕视频在线看片| 人妻 亚洲 视频| 中文字幕最新亚洲高清| 99在线人妻在线中文字幕 | 91字幕亚洲| 老司机午夜十八禁免费视频| 国产又爽黄色视频| 伦理电影免费视频| 99国产精品免费福利视频| 免费在线观看日本一区| 老司机影院毛片| 一本大道久久a久久精品| 亚洲七黄色美女视频| 色老头精品视频在线观看| 最近最新免费中文字幕在线| 不卡一级毛片| 热99re8久久精品国产| 国产日韩欧美亚洲二区| 国产男女内射视频| 欧美黄色片欧美黄色片| 久久 成人 亚洲| 亚洲人成77777在线视频| 高清欧美精品videossex| 大型av网站在线播放| 少妇 在线观看| 两性夫妻黄色片| 老熟妇仑乱视频hdxx| 性少妇av在线| 精品一品国产午夜福利视频| 久久人人爽av亚洲精品天堂| 90打野战视频偷拍视频| 老司机亚洲免费影院| 啪啪无遮挡十八禁网站| 999久久久精品免费观看国产| 久久香蕉激情| 午夜福利视频在线观看免费| 午夜成年电影在线免费观看| 99精品欧美一区二区三区四区| 啦啦啦在线免费观看视频4| 国产成人精品在线电影| 成年人午夜在线观看视频| 又大又爽又粗| 天天躁夜夜躁狠狠躁躁| 亚洲成国产人片在线观看| 久久精品熟女亚洲av麻豆精品| 欧美精品一区二区免费开放| 老司机靠b影院| 一区二区av电影网| 天天影视国产精品| 中文字幕最新亚洲高清| 女人久久www免费人成看片| 日韩免费高清中文字幕av| 亚洲av欧美aⅴ国产| 一区二区三区激情视频| 无人区码免费观看不卡 | 亚洲国产中文字幕在线视频| 深夜精品福利| 免费在线观看日本一区| 91av网站免费观看| 在线观看免费视频日本深夜| 18禁黄网站禁片午夜丰满| 波多野结衣av一区二区av| 夫妻午夜视频| 成年人免费黄色播放视频| 亚洲av成人一区二区三| 成人国产av品久久久| 亚洲综合色网址| 在线永久观看黄色视频| 日日摸夜夜添夜夜添小说| 99国产综合亚洲精品| 国产精品久久久久久精品电影小说| 大香蕉久久成人网| av天堂在线播放| 亚洲综合色网址| 欧美黄色片欧美黄色片| 在线观看舔阴道视频| 国产福利在线免费观看视频| 亚洲精品国产区一区二| 国产成人免费观看mmmm| 久久毛片免费看一区二区三区| 国精品久久久久久国模美| 免费观看av网站的网址| 欧美精品亚洲一区二区| 成年人黄色毛片网站| 精品亚洲乱码少妇综合久久| 黑人猛操日本美女一级片| 欧美日本中文国产一区发布| 99国产精品一区二区三区| a在线观看视频网站| 成人国产一区最新在线观看| 午夜久久久在线观看| 精品人妻熟女毛片av久久网站| 香蕉久久夜色| 欧美亚洲 丝袜 人妻 在线| 91字幕亚洲| 大片电影免费在线观看免费| 老汉色av国产亚洲站长工具| 欧美性长视频在线观看| 香蕉国产在线看| videos熟女内射| 一级毛片电影观看| 久久精品人人爽人人爽视色| 国产精品99久久99久久久不卡| tocl精华| 免费少妇av软件| 91老司机精品| 欧美日韩亚洲综合一区二区三区_| 少妇 在线观看| 黄色片一级片一级黄色片| 午夜福利视频在线观看免费| 色尼玛亚洲综合影院| 国产精品.久久久| 老司机亚洲免费影院| 成年人黄色毛片网站| 免费黄频网站在线观看国产| 久久ye,这里只有精品| 久久久水蜜桃国产精品网| 欧美人与性动交α欧美软件| 成年人黄色毛片网站| 日韩成人在线观看一区二区三区| 国产在线精品亚洲第一网站| 国产亚洲精品久久久久5区| 日本黄色视频三级网站网址 | 两性午夜刺激爽爽歪歪视频在线观看 | 国产精品免费一区二区三区在线 | 麻豆成人av在线观看| 热99re8久久精品国产| 精品人妻在线不人妻| 真人做人爱边吃奶动态| 日韩人妻精品一区2区三区| 18禁观看日本| 久久久精品94久久精品| 亚洲人成77777在线视频| 久久国产精品影院| 亚洲伊人久久精品综合| 国产一卡二卡三卡精品| 一本久久精品| 日本av免费视频播放| 亚洲欧美色中文字幕在线| 黄片播放在线免费| 欧美激情高清一区二区三区| 国产一区二区三区视频了| 欧美一级毛片孕妇| 在线天堂中文资源库| 一边摸一边抽搐一进一出视频| av福利片在线| 久久久精品免费免费高清| 国产免费视频播放在线视频| 一级,二级,三级黄色视频| 99久久人妻综合| 搡老乐熟女国产| 极品教师在线免费播放| 亚洲精品在线美女| 精品午夜福利视频在线观看一区 | 最近最新中文字幕大全电影3 | 久久精品亚洲av国产电影网| 国产一区有黄有色的免费视频| 亚洲欧美日韩高清在线视频 | 久久婷婷成人综合色麻豆| 午夜激情久久久久久久| 日韩熟女老妇一区二区性免费视频| 国产精品免费视频内射| av片东京热男人的天堂| 深夜精品福利| 色老头精品视频在线观看| bbb黄色大片| 久久毛片免费看一区二区三区| 久久精品国产综合久久久| 国产黄频视频在线观看| 啦啦啦免费观看视频1| 日日摸夜夜添夜夜添小说| 麻豆av在线久日| 丁香六月天网| 久久99热这里只频精品6学生| 啪啪无遮挡十八禁网站| 午夜免费成人在线视频| 黄色a级毛片大全视频| 久久久久久久精品吃奶| av线在线观看网站| 国产精品一区二区在线不卡| 51午夜福利影视在线观看| 黄色 视频免费看| 国产精品欧美亚洲77777| 成在线人永久免费视频| 不卡一级毛片| 国产男女内射视频| 老司机影院毛片| 热re99久久精品国产66热6| 国产亚洲欧美精品永久| 午夜福利在线免费观看网站| 在线观看www视频免费| 高清av免费在线| 国产一区二区三区在线臀色熟女 | 老司机亚洲免费影院| 女警被强在线播放| 久久天堂一区二区三区四区| 国产精品熟女久久久久浪| 人人妻人人澡人人看| 欧美老熟妇乱子伦牲交| 色婷婷久久久亚洲欧美| 久久久精品国产亚洲av高清涩受| 亚洲av美国av| 手机成人av网站| 免费不卡黄色视频| 人人妻人人爽人人添夜夜欢视频| 岛国在线观看网站| 国产成人av教育| 深夜精品福利| 国产精品一区二区精品视频观看| 国产区一区二久久| 欧美激情极品国产一区二区三区| 欧美成人午夜精品| 成人三级做爰电影| 国产精品自产拍在线观看55亚洲 | 啦啦啦 在线观看视频|