D I S C R E T E D Y N A M I C V I S C O E L A S T I C S Y S T E M S A N D V I B R A T I O N A N A L Y S I S O F A N E N G I N E S U P P O R T E D O N V I S C O E L A S T I C M O U N T S By Alexander Muravyov M.Eng. (Hons), Chelyabinsk Polytechnic Institute, Russia, 1982 MA.Sc., UBC, Department of Mechanical Engineering, 1994 A T H E S I S S U B M I T T E D I N P A R T I A L F U L F I L L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F D O C T O R O F P H I L O S O P H Y in T H E F A C U L T Y O F G R A D U A T E S T U D I E S D E P A R T M E N T O F M E C H A N I C A L E N G I N E E R I N G We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y O F B R I T I S H C O L U M B I A June 1997 © Alexander Muravyov, 1997 a In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of ^eyeing n L ca & • En&n een'lA The University of British Columbia Vancouver, Canada . Date Ju*e , IP. , 1997' DE-6 (2/88) Abstract An analysis of linear dynamic viscoelastic discrete systems is presented, and its application to the vibration of an engine supported on viscoelastic mounts is discussed. A new procedure is developed which allows exact (closed form) homogeneous solutions in the time domain to be derived for a dynamic system consisting of isotropic viscoelastic components, for which the relaxation kernels are represented as a sum of exponentials. The developed procedure (which is given a name "substitution method") for determination of closed form solutions is extended to the solution of boundary value problems. The application of the substitution method is also extended to the case of periodic loading. Based on this method, a numerical investigation of free and forced vibration responses of some viscoelastic systems is presented. Several approximation techniques are developed in this study which allow the parameters of the relaxation kernels (represented as a sum of exponentials) to be determined from exper-imental data. Also a numerical procedure for determination of complex moduli of isotropic viscoelastic materials is developed, in which certain experimental data related to the material specimen are required as input information. A hereditary (viscoelastic) stiffness matrix-operator is obtained by replacement of the elastic constants in the elastic stiffness matrix by the corresponding viscoelastic operators, or by complex moduli (for steady-state response problems). Comparison of experimental results (in terms of steady-state responses) with the numerical ones is presented. The sufficient conditions of diagonalization of discrete viscoelastic systems are formulated in this study. Analysis of the conditions of overdamping of a simple viscoelastic single-degree-of-freedom ii system is conducted, and new results concerned with this analysis are demonstrated. A particular case of a dynamic viscoelastic system (an internal combustion engine on elastomeric mounts) is given special consideration. A new dynamic model of an engine-mount system is developed where rotating and reciprocating parts lead to the mass matrix and velocity matrix (matrix-coefficient at the velocity vector) as periodic functions of time. The derivation of the equations of motion on the basis of Lagrange's equations is demonstrated. An analysis of parametric resonance phenomena for some examples of engine-mount systems is conducted. A method for steady-state response calculations for the case of time-dependent (periodic) matrices in the equation of motion is developed and some numerical results are presented. An optimization problem is posed and solved with the associated constraints and the objective function reflecting the optimum criteria of the performance of an engine-mount system. As a result, the optimum parameters of the mount material are determined. iii Table of Contents Abstract ii List of Tables vii List of Figures viii List of Symbols xiii Acknowledgement xviii 1 Introduction 1 1.1 Objectives of the research 1 1.2 Background and review of the literature 2 1.2.1 Constitutive laws of linear viscoelasticity 3 1.2.2 Viscoelastic discrete dynamic systems 13 2 Constitutive properties of viscoelastic media 16 2.1 The complex modulus and relaxation kernel concepts 16 2.1.1 Complex modulus in terms of relaxation kernel 16 2.1.2 Relaxation kernel in terms of complex modulus 17 2.1.3 Multiplication and division of operators and calculation of complex moduli 18 2.2 Approximation of experimental relaxation and creep curves 22 iv 2.2.1 Determination of the parameters in the constitutive law and calculation of complex moduli 28 2.2.2 Use of two-part kinematic loading 29 2.3 Approximation of complex modulus curves 31 3 Steady-state solutions for discrete viscoelastic systems 37 3.1 Homogeneous material systems 37 3.2 Procedure for determination of the complex moduli 43 3.3 Inhomogeneous material systems 52 3.3.1 Numerical results and comparison with experimental data 52 4 Closed form solutions for discrete viscoelastic systems 57 4.1 Single-degree-of-freedom system 57 4.1.1 Application of the Laplace transform 58 4.1.2 The substitution method 58 4.2 Multi-degree-of-freedom system. Poisson's ratio is constant 61 4.2.1 Laplace transform method 62 4.2.2 The substitution method 63 4.3 Periodic loading case. Application of the substitution method 67 4.4 Numerical results 68 4.5 The substitution method. Poisson's ratio as a viscoelastic operator 79 4.6 Application of modal analysis technique 82 4.6.1 Conditions of diagonalization of two matrices. Analysis of eigenvalues 82 4.6.2 Conditions of diagonalization of three matrices 85 4.6.3 Diagonalization of discrete viscoelastic systems 89 4.7 Investigation of the conditions of overdamping of a SDOF system 91 v 5 Dynamics of an engine-mount system 100 5.1 Derivation of equations of motion 100 5.2 Parametric resonance investigation 116 5.2.1 Parametric resonance analysis of an engine-mount system 120 5.3 Steady-state response for the case of periodic matrices M and D 123 5.4 Optimization problem for an engine-mount system 128 5.4.1 Numerical results for the optimization problem 133 6 Conclusions and future work 137 Bibliography 141 A Transfer functions, complex moduli 148 B Dimensions of the vibration rig 162 C Derivation of the transformation matrix R (for section 5.1) 164 D Input parameters of an engine-mount system 167 vi List of Tables 4.1 Nondimensional eigenvalues for the selected points 97 5.1 Eigenvalues of monodromy matrix (multipliers). Example 1, stable case . . . 121 5.2 Eigenvalues of monodromy matrix (multipliers). Example 2, unstable case . . 122 D.l Input parameters of the engine-mount system. Units: length - [m], mass -[kg], angle - [degree], Jtj, J?° - [kg * m2], K - [N/m] 169 vii List of Figures 1.1 Maxwell and Voigt models 3 1.2 Standard linear solid model 4 1.3 Kelvin chain and Biot models 4 2.1 1-term exponential approximation. Line with squares - experiment, line with triangles - approximation 26 2.2 2-term exponential approximation. Line with squares - experiment, line with triangles - approximation 27 2.3 Two-part kinematic loading 30 2.4 Complex shear modulus by an 1-term exponential approximation, Gi=real part, (j2=imaginary part, experiment, approximation 35 2.5 Complex shear modulus by a 2-term exponential approximation, Gi=real part, (?2=imaginary part, experiment, approximation . 35 3.1 Viscoelastic brick clamped at the hatched sides 40 3.2 Periodic loading 41 3.3 Steady-state responses: solid line - i/* as frequency-dependent, dashed line -v. = 0.4995 42 3.4 Steady-state responses: solid line - i/„ as frequency-dependent, dashed line -v. = 0.4995 42 3.5 Two types of base excitation test; tension type - 1) and shear type - 2) . . . 45 3.6 3-D finite element model of the specimen with two types of loading 46 3.7 Vibration rig with viscoelastic mounts 53 viii 3.8 Response function (point 5). Material - DPNR. — Numerical, Experimental 54 3.9 Response function (point 6). Material - DPNR. — Numerical, Experimental 55 3.10 Response function (point 5). Material - EAR C-1002. — Numerical, Experimental 55 3.11 Response function (point 6). Material - EAR C-1002. — Numerical, Experimental 56 4.1 Viscoelastic beam with fixed ends 69 4.2 Imaginary part of eigenvalues 70 4.3 Real part of eigenvalues 71 4.4 The 31st eigenvector (1st underdamped) 71 4.5 The 29th eigenvector (analog to the 31st one) 72 4.6 The 2nd underdamped eigenvector 72 4.7 Overdamped eigenvector (analog to the 2nd underdamped) 73 4.8 Free vibration response, case i) 74 4.9 Free vibration response, case ii) 74 4.10 Contribution of eigenvectors, case i) 75 4.11 Contribution of eigenvectors, case ii) 76 4.12 Initial and terminal shapes 76 4.13 Solution for conditions (4.20), T=0.02 s 77 4.14 Solution for conditions (4.20), T=0.016 s . . . 77 4.15 Forcing function 78 4.16 Forced response with zero initial conditions 78 4.17 The boundaries of the region where D < 0 (all roots are real) 95 4.18 Overdamped case region, D < 0 96 ix 4.19 Overdamped case regions (D < 0) with different viscous damping: | = 0;0.5;1 96 4.20 Selected points in the plane of nondimensional parameters 97 4.21 Imaginary part of the complex root in terms of f, | 98 4.22 Real part of the complex root in terms of | , | 99 4.23 Real root in terms of f, f 99 5.1 Schematic view of engine framework and crankshaft 101 5.2 Determination of velocity vector V B . 104 5.3 Determination of kinetic energy of the pistons 106 5.4 Assumed location of the mounts . Ill 5.5 Example of a SDOF system 118 5.6 Regions of stability and instability for the SDOF system (Fig. 5.5) 120 5.7 Steady-state responses. Crankshaft with = J^3°, elastic mounts 126 5.8 Steady-state responses: — M,D constant, — M,D time-dependent. Crankshaft with ^ J^°, elastic mounts 127 5.9 Steady-state responses: — M,D constant, — M,D time-dependent. Crankshaft with J220 7^ J330, viscoelastic mounts 127 5.10 Mount model 132 5.11 Profiles of E2{u>) for different unbalance radii p (m): 1 - p - 0.00014, 2 -p = 0.00016 , 3 - p = 0.00018 , 4 - p = 0.00020 134 5.12 Transmitted force factor T(u>) for different unbalance radii p (m): 1 - p — 0.00014 ,2 -p = 0.00016 , 3 - p = 0.00018, 4 - p = 0.00020 135 5.13 Profiles of E2(u>) for different stiffness ratios: solid line - 72/71 = 1, long-dash line - 72/71 = 0.5. Unbalance radius - p = 0.00018 m 136 5.14 Transmitted force factor T(u>) for different stiffness ratios: solid line - 72/71 = 1, dash-dotted line - 72/71 = 0.5, dashed line - rigid mounts. Unbalance radius - p = 0.00018 m 136 A.l Transfer function in the 1st test for PECH 148 A.2 Transfer function in the 2nd test for PECH 149 A.3 Transfer function in the 1st test for EAR C-1002 149 A.4 Transfer function in the 2nd test for EAR C-1002 150 A.5 Transfer function in the 1st test for DPNR 150 A.6 Transfer function in the 2nd test for DPNR 151 A.7 Transfer function in the 1st test for CR 151 A.8 Transfer function in the 2nd test for CR 152 A.9 Transfer function in the 1st test for NBR 152 A.10 Transfer function in the 2nd test for NBR 153 A.11 Complex Young's modulus of PECH 154 A.12 Complex Poisson's ratio of PECH 154 A.13 Complex shear modulus of PECH 155 A.14 Complex Young's modulus of EAR C-1002 155 A.15 Complex Poisson's ratio of EAR C-1002 156 A.16 Complex shear modulus of EAR C-1002 156 A.17 Complex Young's modulus of DPNR 157 A.18 Complex Poisson's ratio of DPNR 157 A.19 Complex shear modulus of DPNR 158 A.20 Complex Young's modulus of CR 158 A.21 Complex Poisson's ratio of CR 159 A.22 Complex shear modulus of CR 159 xi A.23 Complex Young's modulus of NBR 160 A.24 Complex Poisson's ratio of NBR 160 A. 25 Complex shear modulus of NBR 161 B. l Vibration rig with viscoelastic mounts, dimensions in [inch] 162 B.2 Vibration rig, View A - A 163 xii List of Symbols Roman C viscous damping matrix of system det(...) determinant of a matrix D(t) velocity matrix of engine-mount system exp(...) exponential function E instantaneous Young's modulus E*(<jj) complex Young's modulus E*(ijj) conjugate of complex Young's modulus Ei(u>) real part of complex Young's modulus (storage modulus) E2(<x>) imaginary part of complex Young's modulus (loss modulus) E viscoelastic Young's modulus operator F(t) forcing function acting on system G instantaneous shear modulus G+(LO) complex shear modulus Gi(u>) real part of complex shear modulus G2(w) imaginary part of complex shear modulus G viscoelastic shear modulus operator i = 1, n i=l,2,...,n I unit matrix 3A tensor of moments of inertia of engine framework JA matrix of components of tensor 3A Jf- ij-th component of matrix J A xiii 3B tensor of moments of inertia of crankshaft JB matrix of components of tensor 3B J$ ij-th component of matrix JB JBo matrix of components of tensor in rigidly connected (with the body B) coordinate system Jf° ij-th component of matrix JBo k*stor(w) real part of complex stiffness of viscoelastic spring Kioss(u) imaginary part of complex stiffness of viscoelastic spring K instantaneous stiffness matrix of system KQ matrix-coefficient associated with Young's modulus Ki matrix-coefficient associated with Lame's coefficient A K2 matrix-coefficient associated with Lame's coefficient G Kj ij-th component of matrix K K viscoelastic stiffness matrix operator of system K(t — T) creep kernel I length of crank l2 length of connecting rod M mass matrix of system Mij ij-th component of matrix M Pi ith eigenvalue (or characteristic root) QM vector of forces and moments transmitted through the mounts to the engine framework t time TA kinetic energy of engine framework TB kinetic energy of crankshaft TP kinetic energy of pistons x iv 7\ transfer function value in the 1st test T2 transfer function value in the 2nd test T(u>) transmitted force factor as function of frequency UIA x-displacement of the mass centre of engine framework U2A y-displacement of the mass centre of engine framework UZA z-displacement of the mass centre of engine framework VA velocity vector of the mass centre of engine framework Vjg velocity vector of the mass centre of crankshaft V G velocity vector of point G (the mass centre of piston) x(t) x-displacement (scalar) as function of time x Laplace transform of x-displacement X(t) vector of displacements as function of time X Laplace transform of vector of displacements Greek ai zth viscoelastic material constant in exponential relaxation kernel jpi angle of crankshaft rotation associated with ith crank (or piston) T(t — r) relaxation kernel ro(...) gamma function e uniaxial strain e Laplace transform of uniaxial strain eij strain tensor components A*(u>) complex Lame's coefficient A viscoelastic Lame's operator v instantaneous Poisson's ratio i/»(u>) complex Poisson's ratio v viscoelastic Poisson's ratio operator xv p unbalance radius of the crankshaft c r uniaxial stress a Laplace transform of uniaxial stress t T j j stress tensor components 4>\ angular displacement of engine framework about axis X <f>2 angular displacement of engine framework about axis Y </>3 angular displacement of engine framework about axis Z \& error of approximation of experimental relaxation, creep and complex modulus curves ui frequency <X>A vector of angular velocity of engine framework U>B vector of angular velocity of crankshaft Subscripts A associated with body A (engine framework) B associated with body B (crankshaft, or rotor) p associated with pistons * associated with complex moduli (or coefficients) Superscripts A associated with body A (engine framework) B associated with body B (crankshaft, or rotor) Bo associated with body B and rigidly connected coordinate system M associated with the engine mounts T transpose of a matrix Abbreviations for viscoelastic materials CR chloroprene DPNR natural rubber xvi EAR isodamp NBR acrylonitrile-butadiene rubber PECH polyepichlorohydrin Other abbreviations SDOF single-degree-of-freedom M D O F multi-degree-of-freedom Acknowledgement The author would like to offer sincere thanks to Professor Stanley G. Hutton for his insight and advice in the preparation of this thesis. The author would like also to acknowledge the financial support provided by Defense Research Establishment Atlantic. xviii Chapter 1 Introduction 1.1 Objectives of the research The following were the primary objectives of this study: i) Survey of the literature which is concerned with the treatment of viscoelastic dynamic systems. Review of different constitutive models of linear viscoelasticity, formulations of dy-namic problems for discrete viscoelastic systems, methods of solution, optimization problems concerned with the performance of engine-viscoelastic mount systems. ii) Analysis of the mechanical properties of viscoelastic materials and substantiation of the application of exponential kernels for modelling of their constitutive properties. Inclusion of two parameters in the constitutive stress-strain relation for a viscoelastic isotropic material, namely, Young's modulus and Poisson's ratio, which will be considered as operators of the hereditary type. iii) Experimental determination of complex Young's modulus and complex Poisson's ratio as functions of frequency for a series of viscoelastic materials. Design of an experimental model and creation of the appropriate software for this purpose. iv) Further development of the existing applied methods to treat discrete dynamic viscoelastic systems in the framework of the finite element formulation. Particularly, this includes the development of a method which allows the solutions to be obtained in the time domain and which is more effective than methods using numerical integration schemes. Creation of the appropriate software based on developed numerical methods in order to conduct investigation 1 Chapter 1. Introduction 2 (numerical experiments) of free and forced vibration responses of some viscoelastic systems and to investigate the effect of variation of material's properties. Comparison of the obtained numerical results with experimental ones. v) Formulation of a dynamic model of an engine-mount system where the contribution of rotating and reciprocating parts are taken into account with accuracy. Creation of the ap-propriate software in order to investigate (numerically) free and forced vibration responses for this model of an engine-mount system. vi) Optimization of the performance of engine-mount system in terms of reduction of the force transmissibility through the mounts, which are considered as viscoelastic. Formulation and solution of an optimization problem with the associated constraints and the objective function reflecting the optimum performance criteria for an engine-mount system. Creation of the appropriate software for this purpose. Determination of the optimum properties of the material used in vibroisolation mounts. 1.2 Background and review of the literature Viscoelastic damping provides an energy dissipation mechanism and the use of passive damp-ing mounts (isolators) is expedient, because it is inexpensive and a reliable way of vibroisola-tion, and noise reduction. The use of polymeric (elastomeric) materials for these purposes is quite extensive in the industry - [1], [2], [3], [4], [5]. Hereditary (viscoelastic) properties are present in many elastomeric materials, for example, chloroprene, polyepichlorohydrin, acrylonitrile-butadiene and many others. The definition of a hereditary medium will be given below. In this study, isotropic homogeneous viscoelastic materials and systems which utilize them are considered. It is known that polymers possess linear properties at quite significant strains Chapter 1. Introduction 3 [6], [7]. It will be assumed that the theory of linear viscoelasticity is adequate for the de-scription of constitutive laws of considered materials. It is also assumed that the material is supposed to be in an isothermal state, i.e., the change of temperature is small. Vibration mounts are used to minimize the levels of forces that are transmitted from operating machinery to the supporting structure. In seismic applications the mounts are used to isolate the structure from motions produced by the ground. The problem of choosing optimum mounts is not straightforward. Viscoelastic mounts have complex material properties which are dependant upon the frequency of excitation. The determination of optimum mount properties depends also on optimality criteria which needs to be specified. 1.2.1 Constitutive laws of linear viscoelasticity Elementary viscoelastic models are shown schematically in Figs, l . la .b, 1.2a,b, where K, Elt E2, E0, Eoo are stiffnesses of the springs and C, nE2, bE^ are viscous coefficients of the dashpots. c c K a) b) Figure 1.1: Maxwell and Voigt models Combining such models one can create more complicated models of viscoelastic materials particularly like those shown in Fig. 1.3 a), b) which were discussed in [8]. Chapter 1. Introduction 4 E i 7]E 2 -VW-E ^ b 00 oo •AAV E 0 a) b) Figure 1.2: Standard linear solid model a) Kelvin chain model b) B io t model Figure 1.3: Kelvin chain and Biot models The term "viscoelastic material" has quite a broad meaning, e.g., in the literature there is use of this term if the equation of motion includes a viscous damping term, i.e., the equation of motion is written in terms of current instant values of displacement, velocity, acceleration. In this study the equation of motion will include the integral term and a history of strain (or stress) is required for its formulation. Materials yielding such a constitutive relation (requiring history) are also called as viscoelastic, but the term "hereditary" will be used in this study as well to distinguish them. Models using the viscous linear damping mechanism can be considered as simplest models Chapter 1. Introduction 5 of viscoelastic constitutive laws. Note that these models do not contain hereditary properties. Such viscous damping models were investigated, e.g., in [9], [10], [11], [12], [13], [14], [15], [16]. The application of the component modes synthesis (CMS) along with the finite element procedure is demonstrated in [17], [18], [19], [20], [21], where the solution of the eigenvalue problems (free vibrations) is considered for systems with nonclassical (nonproportional) viscous damping. The case of arbitrary periodic excitation is considered in [19], [20] with the use of a modal representation of the solution. The necessary and sufficient conditions of classicall damping are considered in [13], [20], [22] [23], where an explicit expression of the damping matrix in terms of the mass and stiffness matrices is presented. In the theory of viscoelasticity, one of the models of the stress (tr)-strain (e) relation [24], [25], [26], [27], [28], [29] is a constitutive law of the form: der (Per de cPe aDv + «!-£ + ... + a n — = bae + bx- + ... + bn— (1.1) Note that a constitutive relation in such form is not convenient to use, thus it is expedient to represent it in other equivalent form. Taking the Laplace transform of equation (1.1) one obtains (aa + axp + ... + anp n)a- + i>n-i((ro,<ro\ —,(r0 /l~ 1\a1, ...,an,p) = (K + hp + ... + bnpn)e + £„-i(e 0, $ \ e j , n - 1 ) , h,bn,p) or in abbreviated form , A4p)cr + n^.1(...,p) = Bn{p)e-rCn-i(...,p) (1-2) where V'ri_1(...,p), £n_i(...,/>) are polynomials of order n — 1 , coefficients of which depend on the initial values of cr,...,a^ n~ 1\ e e(n_1) and coefficients at, bk and Laplace transform variable p. If the roots of polynomials An and Bn are real, distinct and negative the material is called an hereditary-elastic (or viscoelastic) medium [24]. Chapter 1. Introduction 6 From (1.2) one obtains It is a customary assumption that for t < 0 the material is assumed unstressed and undeformed, i.e., e(t) = 0, cr(t) = 0 for t < 0 and a history of deformation starts at t = 0 with zero conditions (l) O- 1 ) r> (!) ("-!) n In this case (zero conditions) the polynomials in (1.3) £n-i(...,p) = V'n-i(•••>/>) = 0. Here the conditions at £ = 0 will be generalized, i.e., it will be still assumed that for t < 0 the material has e(t) = 0, a(t) = 0, but at t = 0 the conditions are not absolutely zero ones. The fraction of polynomials Bn(p) and An(p) can be represented in the following form ^ \ = Et1_y2^—) (1.4) where a; > 0, i = l,n (recall that the roots of polynomial An were assumed negative). Using the inverse Laplace transformation and the convolution theorem one obtains from (1.3), (1-4): cr{t) = E(e(t) - [tl£ciexp[-ai{t-T)]e{T)dT) + 0(...,t) (1.5) where function J2ciexp[-ai{t-r)] = V{t-T) (1.6) is called the relaxation kernel. The function ®(...,t) originated from the second term in (1.3), and it depends on o-0 0 o " - 1 \ e0 e0n_1^ a n c ^ coefficients a*., It was verified in this study that if the initial conditions a0 tfo"-1^. e0,...,e0n_1^ are related according to * = E{e - f J2 Ciexp[-ai(t - r)}e{r)dr) (1.7) J o i=i Chapter 1. Introduction 7 then in equation (1-3) £ n - l ( - , p ) -•0n-l(".,p) = 0 and consequently the term Q(...,t) in (1.5) vanishes too. This fact will be illustrated on an example of standard linear solid model below. In [24] it is mentioned that the constitutive relation (1.1) is equivalent to (1.7) under the appropriate initial conditions, though they were not discussed. Here these appropriate initial conditions are specified as those which satisfy relation (1.7). Note that all necessary relations between derivatives of tr(t) and e(t) at t = 0 are obtained from (1.7) by differentiation and then substituting t = 0. One can rewrite relation (1.7) in the inverted form e = i ( < 7 + j f ^ - ^ M i - ) * - ) (i-s) where K{t-r) = '£diezp[-pi{t-T)} i=l is called the creep kernel. Coefficients di, fa can be expressed through coefficients Ci, ai of the relaxation kernel T(t — r) and vice versa. An example of the standard linear solid model of a hereditary-elastic medium is shown in Fig. 1.2a. This model yields the following relation between stress and strain: CT + ao~ = Ei(e + pe) where p, = 1/rj, a = Elr^ 2• Taking the Laplace transform of this equation one obtains , (p + a)cr - a0 = Ei(p + p)e - -Eien which can be rewritten as -a = El?±±e+ = EL(1 - SLZ»re+«LZ^ ( 1 . 9 ) p + a p + a p + a p + ct Chapter 1. Introduction 8 Applying the inverse Laplace transformation and the convolution theorem to (1.9) the follow-ing constitutive relation is obtained tr = EAe — \ (a — p,)exp[—a(t — T)]e(r)dr) Jo Note that the second term in (1.9) vanishes upon the assumption that the inital values are related according to (1.7), which in this case means <r0 = Eie0. The latest relation has a clear physical interpretation (see Fig. 1.2a). The dashpot cannot undergo instantaneous deformation, because of the infinite strain rate. Therefore the spring Ei takes all instantaneous deformation, i.e., a0 = Eie0-Constitutive relations (1.7), or (1.8) of a hereditary material have so-called difference type relaxation and creep kernels, i.e., T(i, r) = T(t - T) K(t, T) = K(t - T) which is an appropriate assumption for polymers, however it is not possible to say the same about some ageing materials, for example, concrete which has properties that vary in time. To justify that these kernels have this property consider the steady-state response case when both functions (stress and strain) are periodic in time with period T, (the fact that this is possible can be verified for any elastomeric material). To have an exact steady-state response at time t < oo, it is necessary to start the loading at time — oo (unless the initial conditions which correspond to the steady-state response at time t = 0 are prescribed). This is why instead of 0 as a lower limit of integration in the constitutive law the value — oo will be used. Thus we have at time t Ee{t) = a(t) + f K(t,T)<r(T)dT J — oo at time t + T rt+T Ee(t + T) = a(t + T) + / K(t + T, r)a(r)dT J — oo Chapter 1. Introduction 9 Recall that cr(t + T) = <r(t) e(t + T) = e{t) Then replacing r by r + T in the integral for time t + T one obtains f K(t,r)cr(T)dT= f K(t + T,T + T)(r(T)dT J—oo J—oo With an arbitrary T this relation is possible if, and only if, K(t,r) can be represented as a function K(t — r), i.e., as a function of one argument K(£), where f = t — r. This condition is also called the Volterra principle of closed cycle. Consider now the concept of a complex modulus. Substituting the harmonic strain: e = eQexp(iijji) i = in the constitutive stress-strain relation <r{t) = E(e(t) - f Y(t - r)e(r)dT) J — oo one obtains cr(i) = Ee0exp(iu>t)[l - f T(t - T)exp(iu}(r — t))dr] J—oo Replacing t — r by £ one obtains /•oo CT = Ee0exp(iut)[l — / r(£)exp(-iu>£)ci£] or in another form a = e0exp(iu>t)[Ei(u>) + iJ52(u;)] where E1(u) = E[l- [°° T{t)cos{u>t)dt] Jo Chapter 1. Introduction 10 and Jo The quantity Ei(u) + iE2(uj) is called the complex modulus, Ex is called the storage modulus, E2 = the loss modulus, E2/Ei = the loss factor (or loss tangent). Analogous formulae are mentioned in [30], [24]. Note that one can find the relaxation kernel r(£) from a given complex modulus E\(tjj) + iE2(w) (this formula will be given in the next chapter). The coefficients C j , ai in (1.6) can be found from an experimentally obtained Some experimental data and descriptions of procedures for determination of viscoelastic characteristics can be found in [7], [31], [32], [33], [34], [35], [36], [37], [38], [39]. In [40] viscoelastic materials such as (re-entrant foams) with possible negative Poisson's ratio were investigated. An example of expressions for the complex modulus is presented below. For a material with properties as in Fig. 1.2b (standard linear solid model), the complex modulus can be expressed as: E-I (<jj) + iE->(u}) — h i- — It has been found that this model gives a qualitative description of the behaviour of some elastomeric materials, but not a quantitative one [41]. In [42], [41] empirical 5-parameter for-mulae for Ei, E2 were proposed, which yield some quantitative agreement with experimental data. Viscoelastic models with a constitutive relation based on fractional calculus were consid-ered in [30], [24], [6], [3], [43], [44], [45]. The definition of the fractional integral of function y(t) as it follows from [46] is: D " » M = lWi7^ M* *>0 (L10) where r0(...) = gamma function. Chapter 1. Introduction 11 The definition of the fractional derivative of order vartheta (i? > 0) of y(t) is D*y{t) = Dn[D-"y{t)] = Dn \{n) Jo (t - r ) i - " n > 0 £ > 0 To where n is the smallest integer greater than and i9 = n — n. For example, for $ C (0,1) the fractional derivative of order •& will be I _d /* y(r) Y0(l - ti) dt Jo (t-ry As an example of a fractional model, the constitutive relation can be expressed as [43]: * » M = rasjfii*7F* t > 0 ( M 1 ) cr(£) + feD^cr(i) = E0e(t) + E^eit) (0,1) where 5 parameters b, Eo,Ei,"&,8 are involved. Another model based on fractional derivative and considered in [45] was <r{t) = ke{t) + cD*e(t) 0C(O,1) where 3 parameters k,c, i? are involved. The fractional derivative (1.11) can be expressed in equivalent form [46] as follows "M-T^tt^+fjhrf t>0 d.i2) This form has only the operator of integration, while form (1.11) involves integration and differentiation, which is less convenient to use. Note that according to definitions (1.11), (1.12) the fractional derivative is not defined for t = 0. Consider some fractional models of the constitutive law which use the concept of the fractional integral (1.10). For example, one can assume the following constitutive relation: a(t) = Ee(t) - e jf* -^-^e(r)dr * > 0 (1.13) Using the convolution theorem and finding the corresponding Laplace images [47] one obtains Chapter 1. Introduction 12 This model corresponds to the Davidson-Cole model according the terminology from [42]. The complex modulus for this model is F(u,\-F CT°{*] One can generalize the relaxation kernel in (1-13) further, introducing a summation on the quantity j3, i.e., *{t) = Ee(t) E Cj { t _ r)1_^(r)dT The complex modulus for this model will be which can be compared with the complex modulus of the Havrilak-Negami model [42]: E-^ = E-((WTW ( L 1 4 ) It has not yet been established what kind of relaxation kernel can yield such a complex modulus as (1.14). In [30], [24] so-called fractional exponential functions E#((3,t — r) were proposed. With their use the constitutive relation can be presented as follows: <r(t) = Ee(t) - j* cE4(3, t - T)e(r)dT d C (0,1) where Note that with i? = 0 the function E#((3,t — r) is a regular exponential function E0{P,t - T) = E^"T) and the first term (n = 0) in the above sum (1.15) constitutes the following kernel Chapter 1. Introduction 13 As a conclusion to the review of fractional models, one can say, that these models produce a good fit to the experimental data with a few parameters (3,4 or 5) in the constitutive relation. However the advantage of using of regular exponential kernels instead of fractional exponential (or fractional) models, as will be shown later, is that analytical (closed form) solutions in the time domain can be derived for discrete viscoelastic systems. 1.2.2 Viscoelastic discrete dynamic systems According to the correspondence (Volterra) principle a solution of the viscoelastic problem can be obtained from a solution of the corresponding elastic problem by replacement of elasticity constants (Young's modulus, Poisson's ratio) in the solution by their hereditary analogs (operators). However this works only in those cases [48], [24] when the elastic solution is found analytically in the space domain, and not applicable when the boundary conditions depend on time. The application of the finite element method to elastic systems leads to the formulation of dynamic problems in terms of mass and stiffness matrices. In the case of viscoelastic systems it is necessary to replace the material constants E, v, or A, G (in the stiffness matrix) by their viscoelastic analogs (operators) E, v, or A, G, or by complex moduli Et, i/*, A*, G* in steady-state response problems. It is a fact that very little data is available for complex Poisson's ratio i/» as function of frequency [2]. Some data related to the complex moduli E„ and G* of polymeric materials can be found in [2] and [7]. Finite element method applications for dynamic viscoelastic systems are usually described in the literature in the context of step-by-step numerical integration (in time) schemes [49], [50], [51], or based on numerical inversion algorithms of the Laplace transformed solution. Description of these numerical methods can be found, e.g., in [52], [53]. It may be noted that the use of numerical integration for the boundary value problem (when conditions are Chapter 1. Introduction 14 prescribed at different points in time) is much more complicated than for the initial value problem. In practice one of the widely used models of the constitutive viscoelastic law is when the relaxation kernel is represented by the sum of exponentials. The existence of homogeneous analytical solutions to the free vibration problem for this case has been shown [54], using the Laplace transform method. In papers [55], [56] a different approach is developed to determine the unknown parameters which are involved in the analytical solution. This approach also yields the formulation of an eigenvalue problem. In papers [4], [57] an introduction of "mini-oscillators" was shown incorporated into the finite element formulation, where a quite restrictive combination of certain exponential terms for the relaxation kernel was assumed in order to obtain symmetrical matrices in the Laplace transformed equation of motion, written in some state-space form (additional dissipative coordinates were introduced). Each "mini-oscillator" must consist of two exponential terms: T(t — T) = aiexp[—oti(t — r)] + a2exp[—a2(t — r)] where the coefficients in this expression must satisfy the following condition a,\Oti - — a2a2 and the relaxation kernel must consist of some number of such couples which leads to even number of terms. This may impose a significant restriction on the behavior of viscoelastic models. To illustrate this circumstance, it is sufficient to consider a material, the relaxation kernel of which is described by an odd number of exponential terms. Multilayered elastic-viscoelastic structures were considered in [58], [59], [60] along with optimization of damping properties (loss modulus) of the material. In [2] the damping layer treatment was considered along with application of tuned dampers. The application of vibration mounts (isolators) for multicylinder engines is quite a common Chapter 1. Introduction 15 practice [61], [62]. In this study some examples involving optimization of the mounts for such type of engines will be considered. In practice engines are attached to structures through multipoint isolation systems, and there are coupling effects in the engine and structure between the various mounting locations. Therefore a traditional optimization based on a single isolator model is not adequate [63]. In [64], [65] the authors described an optimization technique, which is aimed at the mini-mization of forces transmitted through engine mounts, with constraints imposed on the static engine displacements. The engine was modelled as a rigid body subjected to periodic load-ings. Steady-state responses were computed with the following evaluation of the transmitted forces. The Voigt element was used to simulate the properties of elastomeric material of the mounts, hence some inadequacy of this model should be expected (it is known that hereditary properties cannot be described by the Voigt element even qualitatively). The design variables were the stiffness and orientation of each individual engine mount. In [66] the objective of the optimization was the removal of natural frequencies of an engine-mount system from the undesirable frequency range. Voigt elements were also used to model the constitutive properties of the mount material. Some comparison of the numerical results with laboratory measurements was produced in terms of optimized frequencies. Variations of the optimization problem depend upon the choice of an objective function, imposition of different constraints on the design variables, or the level of vibrations. Different numerical techniques [67] can be employed to accomplish the optimization problem. In the present study the force transmissibility through the mounts will be considered as a primary factor. A dynamic model of engine-mount system will be shown where contributions of rotating and reciprocating parts are taken into account exactly. The elastomeric material of the mounts, namely, its complex modulus as function of frequency will be viewed as a design variable. Chapter 2 Constitutive properties of viscoelastic media In the theory of linear viscoelasticity, one of the commonly used hereditary models [30], [24], [27], [68] is a constitutive law of the form: where the scalar function T(t — r) is called the relaxation kernel and E is instantaneous modulus. An alternative form of representing (2.1), which was used by [6], [7], [29] is as follows: which can be reduced to (2.1) by integration by parts. R(t — r) is called the relaxation function. Note that ^p- = ET(t), where £ = t - r, E = R(0). In this study relation (2.1) will be used. 2.1 The complex modulus and relaxation kernel concepts 2.1.1 Complex modulus in terms of relaxation kernel Consider the constitutive stress-strain relation (2.1) for the particular case when the strain is harmonic: (2-1) (2.2) e = 16 Chapter 2. Constitutive properties of viscoelastic media 17 and by considering the steady-state response (i.e., assuming that deformation started at t = — oo) one obtains <r(t) = E{e0eiut - f T{t - r)e0eiurdr) J—oo = E(e 0 e i w t + e0e iut f* T(t - r)e- iu,^d(t - r) J—oo = eoe^Eil - / r(Oe-^ie) = E.(cu)e0e tut Jo where the complex quantity E. = E^u) + iE2(u) = E(l - f°° r(Oe-^ de) (2-3) J O will be called the complex modulus. Note that from (2.3) it follows that E;(u;) = E.(-u) (2.4) i.e., Ei(u>) — Ei(—u>) - an even function, E2{u) = — E2( — u) - an odd function. These properties will be used below. The quantity E\ is called the storage modulus, E2 — the loss modulus, | £ = the loss factor (or the loss tangent). 2.1.2 Relaxation kernel in terms of complex modulus Now one can show how to derive the relaxation kernel T(£) from a given (known) complex modulus. Using (2.3) and assuming T(£) = 0 for £ < 0 one obtains £r(t)e-**dt = l-±{E1+iE2) Taking the inverse Fourier transform of both sides one obtains Chapter 2. Constitutive properties of viscoelastic media 18 1 f°° Ei E2 Ei E2 = — / [COSLVLi —cosuti + —sinwl; + i(sinu>ti —sinwi; + —cosu>c;)}dw 2TT J-OO E E E E Taking care of the odd and even terms, i.e., using relation (2.4), one obtains: m = \ r°((i • it)cosujt+^sinu)^ (2-5) In the literature, e.g. in [24], there are alternative formulae which were derived by using Fourier cosine, or sine transforms: 2 Ei r(0 = - f°° (1 - ^coscoCdw (2.6) TT JO hi no = - r (2-7) 7T Jo ill Note that formula (2.5) follows from (2.6) and (2.7) as well. When an experimentally obtained complex modulus E+ = Ei+iE2 is substituted into (2.6) and (2.7), it can result in different values of T(£). This can occur, because the constitutive law (2.1) is just an approximate model. Thus it may be expedient to take the average of the values of T(£), which is done in formula (2.5). 2.1.3 Multiplication and division of operators and calculation of complex moduli The stress-strain relation for an isotropic elastic material has the form ( T j - j = \Sij9 + 2 ( j £ j j where e ~ V e k k X - ( l + u)(l-h (l + )(l-2u) 2(1 + 1/) For an isotropic viscoelastic material this relation can be rewritten as o~ij = \Sij9 + 2Geij Chapter 2. Constitutive properties of viscoelastic media 19 where viscoelastic Lame's operators are A = , ^ r G E (1 + *)(1 - 2£) 2(l + i/) In the development below the relaxation kernels of viscoelastic operators are assumed in the form of a sum of exponential functions. Note that there is no restriction on the number of terms in this sum. For example, for the operator E: Ex = E(x(t) - f T(t - r)x(r)dr) the kernel r(£) will be represented as r(fl = »=1 Introducing the auxiliary definition of operator M*(r): M*(r)x= f er^x{T)dr J—oo Assuming, for now, a kernel T(£) = ae"a£, one can rewrite the expression for the operator E: E = E(l - aM*{-ct)) Assume an analogous expression for the operator v u = u(l- bM*{-8)) Consider their product Ev = E(l - aM*(-a))u(l - bM*(-8)) According to the theorem of multiplication of such kind of operators [24], there exists an equality: M . W M > ) = * » - * » ( , ) ( 2 . 8 ) r — s Chapter 2. Constitutive properties of viscoelastic media 20 Thus Ev = Eu[l - bNT(-P) - aM*(-a) + ° b AM*(-a) - M*(-/?))] —a + p Now let x = e i u , t and calculate (Eu)x, it will be b a a b where (Eu)x = Evil - -r^-. ^ r - + -.—" , w , " . J e " " = P.x v y v /3 + %LO a + icu (a + tu) {(3 + *«)7 JP .E (1 ^ ^ *^ ^ (2 9) * (3 + iu> a + iu> (a + iu>) (/3 + iu) will be called the complex modulus of the operator P = Ev. From the other side, consider separately the complex moduli of operators E, and v. One has Ee™1 = Ell %-)eiwt OL + IU) thus the complex modulus of operator E is a and thus (3 + iu Comparing (2.10), (2.11) with (2.9) one concludes that P. = E.v, E. = E(l — ) (2.10) a + tu> veiut = 1/(1 - ———-,—) e'w* Therefore the following property holds. Property 1. The complex modulus of an operator, which itself is a product of two operators, is equal to the product of complex moduli of these two operators. Chapter 2. Constitutive properties of viscoelastic media 21 The extension of this property to the case of multiplication of several operators (each of them can contain several exponential terms) is straightforward and for the sake of brevity is not shown here. The case, when a = 0, deserves some attention. In this case formula (2.8) is not applicable (r = s). There is another one, mentioned in [24], which can be used, but for our purposes of proving Property 1, it is enough to consider the limit transition at a — 0 —> 0: aM*(-a)bM*{-a)e iu,t = ^Q[aM*(-a)bM*(-0)e iu' t] r r a h i 1 1 M iut ,. r a h ( -a + /3) 1 i u t = lim -( : —)le — lim [- —- . . , . r-r|e a _ ^ o l - a + /3va + to; 0 + tw" a-/3-+oL(-a + 0) (a + IUJ){0 + tw)J ab j„t — -e (a + iu>)(a + iw) Now by considering (2.9), (2.10), (2.11), it becomes obvious that Property 1 will be also valid for this particular case when a = 0. For future purposes it is also necessary to consider the division of operators. Let - _ E _ E(l-aM*(-a)) ~ ~Z ~ i/(l - bM*(-0)) Using the multiplication rule (2.8) one can derive [1 - &M*(r)] _ 1 = l + bM*(b + r) Using this result one obtains: - E E ab D = - = - (1 + bM*(-0 + b)- aM*{-a) — — (AT{-a) - M*{-0 + b))) V V —OL + 0 — 0 Now applying operator D to e Df* = (^i + b- a- —E* _)e*" v 0 — b-\-iu) a-\-iu> (a + iuj)(0 — b + iu>) Thus the complex modulus of D is E . b a ab . ^\ * ~ + 0 - b + iuj~a + iu> (a + iu)(0 - b + iuy ^' ' Chapter 2. Constitutive properties of viscoelastic media 22 From the other side, consider the ratio of complex moduli * = ^ f c ) ( 2 - 1 3 ) After some rearrangement of fractions in (2.12), (2.13) one can see that Therefore the following property holds. Property 2. The complex modulus of an operator, which itself is a fraction of two operators, is equal to the fraction of complex moduli of these two operators. The extension of this property to the general case of operator D = f ^ - f ' (2.14) BiB2...Bm (each of the operators can contain several exponential terms) is straightforward and for the sake of brevity is not shown here. Proposition. The complex modulus of operator D in (2.14) is expressed in terms of the complex moduli of participating operators as following: BuB2*...Bmt The proof follows from Properties 1 and 2. Remark. The multiplication rule (2.8) has been proven in [24] for so-called exponential fractional operators with the same indices. In this study a particular case (subclass) of such exponential fractional operators is utilized, namely, the kernels of operators are assumed to be a sum of exponential functions. 2.2 Approximation of experimental relaxation and creep curves The objective of this section is to show that the use of exponentials for the relaxation kernel in a constitutive relation allows a good modelling of viscoelastic properties to be obtained, and Chapter 2. Constitutive properties of viscoelastic media 23 also to show some methods of approximation developed in this study. These methods allow the parameters of the constitutive relation to be determined from experimental relaxation, or creep curves. Assume that an experimental stress relaxation curve is given on a set of points in time: tk with values <rfc, k = l , m . Assume an approximation as a sum of exponential functions: (TkKij^die-^x + B i=i where unknowhs are di, ji, B, and n is a chosen number of terms. One can set the approximation problem as follows: m n * = Y^Wk - die~litk - B]2 -»• Minimum (2.15) fe=i j=i There are many numerical algorithms of unconstrained minimization [67], for example one of them the conjugate-gradient method, which can be applied to the problem (2.15). However due to the fact that analytical expressions for the derivatives on the design variables are available Newton's method can be employed. At first the minimization problem is posed with respect to parameters di, B, while pa-rameters 7 ; are assumed known (their prescription is discussed below). In this case the minimization problem can be reduced to the solution of a system of linear equations with respect to unknowns dit B. Namely, differentiating (2.15) with respect to di and B and setting these derivatives to zero (necessary condition of an extremum point) one obtains: £ f c e - ( 7 i + 7 i ) t * £ f c e-(7i+72)t* ... Y,ke~[~ll+ln)tk Efee - 7 1 '* dx Efee-^-H' 1)** £ f e e - ( 7 i + 7 2 ) t * ... £ f c e-(n+-rr,)tk ^ e"7'4* di J2ke-^"+^tk Efe e _ ( 7" + 7 2 ) t* ... Efe e~ ( 7 n + 7 n ) t* E f e e-7"'* dn £ f e e " 7 l t i Eke'12'" .» Efee-7"4* n \ [ B (2.16) E ^ e " 7 1 Efe Cke"1' E f e c T f e e - 7 ' Efe f^e Chapter 2. Constitutive properties of viscoelastic media 24 If the determinant of the system (2.16) is not zero, then an optimum point exists and is unique. Obviously it will be the global minimum point (not a maximum). The solution of this linear system yields the optimum values of di, B for prescribed values of 7;. As an illustration consider original data for the relaxation Young's modulus curve (which is actually a stress relaxation curve divided by the value of the constant strain) of polyisobutylene (at 25° C) taken from [7]. This curve is shown in Fig. 2.1 by the line with the squares. With an 1-term approximation, a search of optimum values di and B was undertaken for the interval 71 C [0,0.le + 06] and the increment A71 = 0.25e + 04, i.e., with a given set of values 71 the optimum values for di and B were calculated from the system (2.16) and the best approximation point (minimum \P) was selected. These values were obtained as follows This approximate set of parameters was then used as an initial point for Newton's method which will be described below. Remark. An analogous technique can be employed when an experimental curve is a creep curve. Then the experimental values will be the values of strain on a set of points in time. In this case the values <rfc should be replaced by values of strain ek in (2.16). Another numerical method which is proposed to use for approximation purposes is New-ton's method. Presume that a system of nonlinear equations is given as following: B = 0.86775 di = 2.90729 71 = 2500 with * = 0.62405 (2.17) /(*) = 0 where X = [xi,XN]T, and f(X) is a vector-function f(X) = [f1(X),...,fN(X))T Chapter 2. Constitutive properties of viscoelastic media 25 or in abbreviated form where matrix Ri has the form -Xi+i = X{ — R{ 1f(Xi) ajx aix. dxi 8x2 dx„ ah ah ah dxi 8X2 dxn au au dxi 9x2 dxn (2.18) and derivatives are taken at point i. For our case vector X will consist of unknowns parameters [di,...,<i n , 7 i , . . . , 7 n , £ ] T . The necessary conditions of an extremum point of (2.15) give a system of N = 2n + 1 nonlinear equations as follows: f* = dd; = 0 * = 1 , n fs+n = ^ — = 0 s = l,ra * - m - r> In more detail: /. = E[e_7*tfc(Erfie_7jtk + * - = 0 «= 1, k=l j = l /.+» = E[(-^ e~7'ti)(E^e_7itfc + ^ - c T f e ) ] = 0 * = l,n fc=i j = i m n / 2 n + l = £ ( £ ^ e _ 7 j t f c + B - C T f c ) = 0 fe=i j = i The expressions for the components of matrix Ri are obtained by straightforward partial differentiation of the above formulae for f s . This algorithm was implemented in a program and some results are presented below. Chapter 2. Constitutive properties of viscoelastic media 26 ro 5 2 CD 0.000 0.005 0.010 0.015 0.020 Time, s Figure 2.1: 1-term exponential approximation. Line with squares - experiment, line with triangles - approximation For example, the set of parameters (2.17) was taken as an initial point. Then it required four Newton's iterations to yield the optimum point: B = 0.82902 dx = 2.82834 7 l = 1996.9 with * = 0.57584 The approximation graph (line with the triangles) for this set of parameters is presented in Fig. 2.1. Analogously this combination of two methods was undertaken for a 2-term approximation. With a 2-term approximation, a search of optimum values di, d2, B was undertaken for the interval of changing of 71 C [0,5000], 72 C [0,5000] and increments A71 = A72 = 100, i.e., with a given set of values 71, 72 the optimum values for d\, d2, B were calculated from the system (2.16) and the best approximation point (minimum \P) was selected. These values were obtained as follows B = 0.6405 di = 2.5230 7 l = 3300 d2 = 0.6192 72 = 100 with * = 0.1036 This set of parameters was then used as an initial point for Newton's method. It required 27 0.020 Figure 2.2: 2-term exponential approximation. Line with squares - experiment, line with triangles - approximation three iterations to yield the optimum point: B = 0.6462 dr = 2.4753 71 = 3449.9 d2 = 0.667 72 = 132.38 with * = 0.1017 The approximation graph (line with the triangles) for this set of parameters is presented in Fig. 2.2. Analogously one can conduct approximations with an arbitrary number of exponential terms. Note that using the 1st method (system (2.16)) one can determine the optimum point with a good accuracy, i.e., if small increments on the assigned values of 7; are taken. This would mean that the use of Newton's method is not really necessary. However if the number n of exponential terms is large, then it will require a large number of solutions of system (2.16). Denote the number of increment points on each parameter 7; as L, then the number of solutions of (2.16) will be Ln. For example for n = 5, L = 20 it would require 3.2e + 06 times to solve the system (2.16) (in this case 6x6) which would take a significant amount of computer time. Therefore the use of Newton's method is expedient for the case when number of exponential terms is chosen > 3. Chapter 2. Constitutive properties of viscoelastic media co 5 CL 2 =& 0.000 0.005 0.010 0.015 Time, s Chapter 2. Constitutive properties of viscoelastic media 28 2.2.1 Determination of the parameters in the constitutive law and calcula-tion of complex moduli The constitutive stress-strain relation is assumed in the following form: <r(t) = E(e - ake-ak{t-T)e{T)dr) (2.19) where the parameters E, ak, ak are to be determined, and re is a chosen number of terms. Having represented the relaxation (or creep) curve as a sum of exponential functions, one can without difficulty determine the parameters in the constitutive law and compute the corresponding complex modulus. Assume that the variation of stress upon the application of constant strain e0 is as follows 43 = * + e° i=i From the other side using the constitutive law (2.19) one obtains 43 = E(l - f £ aie-^-^dr) = E(l-J2- + it ~e~ait) e o Jo ttt fri fri ai Thus the following relations between the parameters in the expression for relaxation modulus and parameters in the constitutve law hold: ai = 7 i a i = diai/E E{1-Y,—) = B (2.20) i=i ai Therefore having determined di, 7 ^ B from experiment, the parameters of the constitutive law (2.19) are determined from (2.20). Then the corresponding complex modulus can be computed using (2.3), which produces the following ai ^[ai + iu Assume now that an experimental creep curve is given, which is already represented as e(t) (TO = J5 + E* e _ 7 <* (2-21) i = l Chapter 2. Constitutive properties of viscoelastic media 29 From the other side using the constitutive law in the form: e(t) = i(<7 + T £ cke-W-T)a{T)dT) (2.22) h J o fc=i parameters E, Bk, ck are to be determined. Substituting a(t) = <r0 into (2.22) one obtains e(t) Co JO . = 1 f = 1 # - = 1 /Ji Comparing it with (2.21) the following relations follow: 0i = l i ci = -Edi0i E-\\ + Yt^) = B (2.23) »=i Therefore having determined di, 7; , B from experiment, the parameters of the constitutive law (2.22) are determined from (2.23). Then the corresponding complex compliance can be computed analogously to (2.3) and the corresponding complex modulus is just the inverse of the compliance: E,(w) = ./.(a;)"1 2.2.2 Use of two-part kinematic loading Consider a kinematic loading (Fig. 2.3) in a uniaxial test. This graph has two portions: 1st one when the strain is a linear function of time, an 2nd one when it is constant. This type of loading can be realized experimentally. Note that relaxation curves consid-ered in section 2.2 imply that the strain is applied instantaneously and kept constant. An experiment to realize the instantaneous application of strain is, strictly speaking, impossi-ble. Therefore relations are developed in this section which can be used to determine the Chapter 2. Constitutive properties of viscoelastic media 30 0 S Time Figure 2.3: Two-part kinematic loading parameters of the constitutive relation by using the type of kinematic loadings shown in Fig. 2.3. Note that by the "strain" here is implied either tensile, or shear strain, depending on the type of experiment which is conducted with a specimen. The strain as a function of time is assumed to be t > s with the assumption that for t < 0 the specimen was in natural state (zero stress and strain histories). According to (2.19) for any arbitrary instant t > s one obtains <r(t) = E(e„ - f E ake-a^e(r)dr) J o fe=i = E(eo- rJ2ake-akit-T)-TdT- f ake~a^^e0dr) J o fe=i s J" k=i E[e0 - eo E - - E e-ak\-ekak - e0^e«»)} fe=i afe fe=i ak Chapter 2. Constitutive properties of viscoelastic media 31 where k = l , n Having the experimental stress-time dependence for t > s approximated by <r(t) = B + Yidke-'»*t k=i one can express the parameters of the constitutive relation (2.19) in terms of experimentally determined quantities B, dk, jk- Namely, they will be as follows dk = Ik (2.24) where E is viewed as the instantaneous Young's modulus for the case of a tensile experiment, or the instantaneous shear modulus for the case of a torsion experiment. The formulae (2.24) along with experimental input data can be used to evaluate the parameters of the constitutive relation ai, ait E. Note that with the limit transition s —>• 0 formulae (2.24) converts to formulae (2.20) which corresponds to instantaneous (step-type) kinematic loading. To calculate this limit it is necessary to use L'Hopital's rule. In the following section another approach will be presented which allows the constitutive parameters to be obtained through the approximation of complex modulus curves. 2.3 Approximation of complex modulus curves If the experimental data are given in terms of a complex modulus then an analogous approx-imation technique can be employed as it was described in the previous section. Chapter 2. Constitutive properties of viscoelastic media 32 The objective of this section is to show some methods of approximation developed. These methods allow the parameters of the constitutive relation to be determined from experimental complex modulus curves. Assume that experimental complex modulus is given on a set of points in the frequency domain: u)k with values E*k = Elk +iE2k, k = 1,N. This complex discrete function will be approximated here as following E.k « EU = E ( l - ± — £ - ) (2.25) where parameters aj, aj, E to be determined, and n is a chosen number of terms. Note that it follows from (2.3), that representation (2.25) is equivalent to representation of the relaxation kernel in the form r(t-r)=E«ie-°'M The function to be minimized is introduced as follows * = E(E:k-E.k)(W,k-E,k) = J 2 [ E ( l - ± - ^ ) - E . k ] [ E ( l - ± - ^ ^ )-E.k] This is a real-valued function, because the multiplication of conjugate quantities is involved. One will introduce new notations bj = Ea,j and in the first method of approximation, it will be assumed that quantities aj are known and quantities bj, E are to be determined. Setting the derivatives to zero "L = j;w-^ + —^±-K- + —!±-ftd) = o m = i , » obm £r[ a m - iiok am + iujk ~ aj - iuk am - ivk | | = Y:(Re[E - ± - \ - - E*]) = 0 (2.26) one obtains a linear system of n + 1 equations. Note that each of the quantities Jj^ , and | § yields a sum of two conjugate terms. Here in expressions (2.26) the real part of one of these Chapter 2. Constitutive properties of viscoelastic media 33 terms is set to zero. The solution of this linear system yields the optimum values of bj, E for prescribed values of ctj. As an illustration consider original data for the complex shear modulus G* = G\ + iG2 of styrene-butadiene rubber at 25° C taken from [7]. This curve is shown in Fig. 2.4 by the solid line. With an 1-term approximation, a search of optimum values 6 1 ( E was undertaken when the interval of changing of a i was [0,5000] and the increment A a i = 100, i.e., with a given set of values ct\ the optimum values for blt E were calculated from the system (2.26) and the best (minimum \P) approximation point was selected. These values were obtained as follows E = 0.803128 h = 39.87521 a i = 100 with = 7.0172e - 02 (2.27) Note that with one term approximation the search interval for a can be evaluated by noticing that for the complex modulus Em(u>) = E i + iE2 = E(l the following fact follows ^ = El (a 2 -w 2 ) du> {a2+u;2)2K ' i.e., ^ = 0 if a = u>. Thus if E2(w) attains its maximum at some point com, then having this maximum point from the experimental curve, one can assign the appropriate search interval for a near u>m to seek the best approximation. The approximation point (2.27) was then used as an initial point for Newton's method. The procedure of Newton's iterations was described in the previous section in detail. Here it is applied to the solution of a nonlinear system of equations: dO> N —E 1 n b 1 If- = EM— -^ + — h - E — ^ - + — = o ™ = i.» obm ^[ am - tuh am + iujk ~ ay - iu}k ocm - tujk » = £ ( & [ J ^ _ ^ + ± _ J L _ + 1 - ^ E ^ ] ) - o Oam (am-iu)k)2 (am + iuk)2 £ J aj - iuk [am - tuky Chapter 2. Constitutive properties of viscoelastic media 34 dvp N n b-°E fc=i j=i ao ~ lU)k Note that am (m = l ,n) are now unknown parameters. Newton's method provides a fast convergence when an initial point is chosen in a small enough vicinity of the solution point and may diverge if this condition is not satisfied. Thus it is expedient at first to apply the 1st method (linear system (2.26) solution) to obtain an initial point and then to apply Newton's method. With the initial point (2.27) five Newton's iterations were required to obtain the optimum point with values: E = 0.7869 ^ = 29.683 « i = 72.956 with * = 6.87901e - 02 The approximation graph is presented in Fig. 2.4, where the solid line is original data, and the dashed line is approximation results. The same combination of two methods was then used to obtain a 2-term approximation. For the 1st method, a search of optimum values bi,b2, E was undertaken when the interval of changing of cti and a2 was [0,1500] with the increments Acti = A a 2 = 30. This method gave: E = 0.93851 bx = 8.2443 c*i = 30 b2 = 317.8787 a2 = 1200 with * = 1.6619e - 02 This result was then used as an initial point for Newton's method, and six iterations were required to obtain the optimum point with values E = 0.91863 fei = 4.8562 a i = 17.724 b2 = 230.60 a2 = 830.18 with * = 1.4631e - 02 The approximation graph (dashed line) for this set of parameters is presented in Fig. 2.5. Chapter 2. Constitutive properties of viscoelastic media 35 100 150 200 Frequency, Hz 300 Figure 2.4: Complex shear modulus by an 1-term exponential approximation, Gi=real part, (?2 =imaginary part, experiment, approximation 100 150 200 Frequency, Hz Figure 2.5: Complex shear modulus by a 2-term exponential approximation, (7i=real part, (?2=iniaginary part, experiment, approximation Chapter 2. Constitutive properties of viscoelastic media 36 As a conclusion one can say that the use of exponentials for the relaxation kernel provided a good fit to experimental data for both curves (relaxation modulus and complex modulus) in the considered examples. The constitutive properties of many elastomeric materials can be satisfactorily described by relation (2.19) and the use of exponentials one can find in many applications, see, for example, [49], [24], [25], [69], [70]. In the American literature a sum of exponentials used in the relaxation kernel is sometimes referred as Prony's series. Chapter 3 Steady-state solutions for discrete viscoelastic systems 3.1 Homogeneous material systems The application of the finite element method to an elastic system yields the mass matrix of the system M and the stiffness matrix K. The derivation of the global system stiffness matrix starts with an element stiffness matrix, where the constitutive law (Hooke's law) is involved (3-1) where [a] = [an &22 033 °"i2 0"i3 0 2 3 ] T is the column of stress tensor components, and [e] = [ e n 622 e 3 3 e 1 2 e i 3 e 2 3 ] T is the column of strain tensor components. For an isotropic material matrix [E] has the form: [E] A 0 0 0 where Lame's coefficients A, G are A + 2G A A A + 2G A 0 0 0 A 0 A 0 A + 2G 0 0 0 0 0 0 0 0 0 0 A vE 2G 0 0 0 2G 0 0 0 2G E (1+ i/)(l-2i/) 2(l + i/) 37 Chapter 3. Steady-state solutions for discrete viscoelastic systems 38 The hereditary analog of this matrix is obtained by replacement of E, v by the corre-sponding viscoelastic operators, i.e., E —> E, v -¥ u, or by another replacement A —> A, G -> G. One can see that relation (3.1) can be presented in the two-part form: [a} = [\I1 + 2GI][e} (3.2) where h = 1 1 1 0 0 0 1 1 1 0 0 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 and / is the unit matrix. Then the global elastic system matrix, due to the previous form (3.2), will have the following two-part form [54]: K = \Ki + 2GK2 (3.3) Note that matrices Kx, K2 do not depend on material properties. The equation of motion for harmonic excitation of a viscoelastic system can be written as: MX + (XKi + 2GK2)X = Foe^* (3.4) where viscoelastic operators A and G for the case of steady-state response have the following form: A(.) = A[1(.)- f R(t-r)(.)dr] G(.) = C7[l(.) - f S(t-r)(.)dr] Chapter 3. Steady-state solutions for discrete viscoelastic systems 39 Recall that in the steady-state response consideration the deformation is assumed to begin at t = —oo. This is why as a lower limit of integration in the expressions above the value —oo is used. One can seek the steady-state response solution of (3.4) in the form X = X0eiu,t Substituting it into (3.4) one obtains ( _ o ; 2 M + x,Ki + 2G.K2)XQei"t = F0eiut (3.5) where A*, G* are complex Lame's moduli and can be computed using the rule (2.3) and the properties 1,2 of section 2.1.3, namely, _ vmE. E. (1+ ! / , ) ( ! -2 i / . ) 2(1 + 1/.) Therefore the steady-state solution (in real form) is X(t) = Re {[-Lo2M + A . i f i + 2G. / i : 2 ] - 1 Foe i ' j t ] Extension to the general case of periodic loading with a period T = ^ is accomplished by using the Fourier series and by application of the principle of linear superposition. Namely, for F(t) = Y, FkeiUkt Lok = kco, (3.6) k=0 the solution is X(t) = Re ] T [ - ^ M + M<"*)*i + 2Gt(uk)K2]-1Fkei^t lk=0 A program was written which calculates the steady-state response according to the for-mulae presented above. Matrices M, K\, K2 used in the equation of motion (3.5) can be built by a finite element code written for elastic dynamic systems. Here the program [71] will Chapter 3. Steady-state solutions for discrete viscoelastic systems 40 Figure 3.1: Viscoelastic brick clamped at the hatched sides be used, which was modified by the author of this study for this purpose. As an illustration some numerical results are given below. Numerical examples. Numerical results are presented for a viscoelastic brick loaded by periodic concentrated forces F(t) and clamped at the hatched sides (Fig. 3.1). The parameters of the brick's cross section were 0.1 x 0.1 m, the length = 1 m. 20-node solid finite elements [71] were used to mesh the brick. The periodic forcing function F(t) is shown in Fig. 3.2 for the duration of one period (T = 0.1 s), thus the 1st harmonic frequency u>i in (3.6) is 27r X 10 rad/s (or 10 Hz). The Fourier series of this loading function has a noticeable contribution of harmonics up to the order of 9 (frequency 90 Hz). As far as Young's complex modulus = Ei(oj) + iE2(w) is concerned the experimentally obtained values (Fig. A.11, Appendix A) for polyepichlorohydrin (PECH) were used. The Chapter 3. Steady-state solutions for discrete viscoelastic systems 41 0.00 0.02 0.04 0.06 0.08 0.10 Time, s Figure 3.2: Periodic loading procedure used to obtain these values from certain experimental input data will be discussed in the next section. To illustrate the effect of Poisson's ratio the computations of steady-state response were produced for two models. In the 1st model complex Poisson's ratio was frequency-dependent i/* = i/,(a;) and its values were taken from experimentally obtained data (Fig. A. 12, Appendix A). In the 2nd model Poisson's ratio was assumed as frequency-independent and equal v = 0.4995. Note that Young's complex modulus was the same for both models. The density prescribed to the material was equal 300 kg/m3. Calculated steady-state responses are shown in Fig. 3.3, where F-displacements of the force application node (node 1) and node 2 (Fig. 3.1) are presented as functions of time for the duration of one period T = 0.1 s. The solid line corresponds to the displacement for varying v* model, and dashed line to v* = const model. One can see some difference between the steady-state responses for these two models. Particularly the 2nd model predicts larger maximum displacement. The results when the period T of the same type of forcing finction was 0.05 s (the 1st Chapter 3. Steady-state solutions for discrete viscoelastic systems 42 -0.0010 0.00 0.02 0.04 0.06 Time, s 0.10 Figure 3.3: Steady-state responses: solid line - as frequency-dependent, dashed line i/, = 0.4995 0.0010 £ 0.0005 k c OJ E CD " r ro c Q . CO T3 -0.0005 0.020 0.030 Time, s Figure 3.4: Steady-state responses: solid line - i/» as frequency-dependent, dashed line -v. = 0.4995 Chapter 3. Steady-state solutions for discrete viscoelastic systems 43 harmonic frequency u>i is 2ir x 20 rad/s, or 20 Hz) are presented in Fig. 3.4. The Fourier series of this forcing function had noticeable contribution of harmonics up to frequency 180 Hz. Again one can see noticeable difference between the steady-state responses for these two models. For this case the 1st model predicts larger maximum displacement. It should be noted that analogous results were also obtained for such material as CR. Its material properties such as E+, v* are presented in Fig. A.20, A.21. One can see (Fig. A.21) that Poisson's ratio is almost constant, and as it was expected, the difference between steady-state responses for these two models was negligible. As a conclusion one can state the following. The use of viscoelastic models where Poisson's ratio is constant can be inappropriate in some cases. The experimental data, for example, for such material as PECH indicate that for the frequency range 10 to 250 Hz its value varies in the range between 0.4995 to 0.406, and differences in the calculated steady-state responses for the considered example (Fig. 3.3, 3.4) confirm this conclusion. 3.2 Procedure for determination of the complex moduli A procedure will be presented below which was used to experimentally obtain the complex Lame modulus A, and shear modulus G*. Note that given these two moduli all the other char-acteristic moduli such as complex Young's modulus, or complex Poisson's ratio, or complex bulk modulus can be computed. Some experimental procedures were described in [31], [32] for specimens that had the shape of a long bar which was subjected either to tension, or torsion vibration tests. The assumption of uniform tensile, or shear stress state for such a long rod specimens was made. In this study the specimen's shape is a 3-D cubic one, and no simplifying assumptions about stress (or strain) state will be made. As it will be seen from the procedure an arbitrary shape of specimen can be used. In this approach it is necessary to have an accurate finite element Chapter 3. Steady-state solutions for discrete viscoelastic systems 44 model of the specimen. For viscoelastic systems, Lame coefficients A and G should be replaced by their hereditary analogs (operators). In the case of the finite element method this implies the replacement of material constants A and G in the stiffness matrix by their viscoelastic analogs (operators) A, G, and for a particular case (harmonic deformation) this replacement leads to corresponding complex moduli A», G> which depend on frequency of excitation. To obtain steady-state responses in analytical form for periodic forcing functions, the forces acting on a viscoelastic system are represented in the form of the complex Fourier series. A steady-state response (also in the form of the complex Fourier series) can then be determined, provided two complex moduli for the material of the system are known. The determination of these complex moduli is considered in this section. Consider an example of harmonic excitation of the systems shown in Fig. 3.5. Base harmonic excitation tests are conducted and transfer functions fl(u) = y O l / z t U = y02/«02 are measured for both tests. At first, assume the 3-D element connecting the mass m and the base (Fig. 3.5a) is modelled as a one-dimensional viscoelastic spring of length L and cross-section area A. Considering a single-degree-of-freedom (SDOF) equation, one obtains: m'z + k*z = —mx where z — y — x, x — xoexp[iu)i\. Assume z = zoexp[iu>t] and substitute into the equation of motion. The term k*z equals to b(Ei(u) + iE2(w))z, where b = A/L. Thus mxoU>2 Z ° = - m w 2 + i ( £ i H + i £ 2 H ) and y = x + z = (x0 + Zo)exp[iu>t]. Chapter 3. Steady-state solutions for discrete viscoelastic systems 45 m BASE 1 1 m BASE 1 a) 1st test b) 2nd test reference system Figure 3.5: Two types of base excitation test; tension type - 1) and shear type - 2) The transfer function will be: T = y = XQ + ZQ = bjE^u) + iE2(u)) Now instead of a viscoelastic spring consider a finite element model which will be used to simulate the experimental results. This model includes as parameters the complex Lame moduli A , and G*. 20-node solid finite elements are used to mesh a given 3-D domain (Fig. 3.6). One side of the specimen is assumed fixed, and for the opposite side some displacements will be prescribed (side T). Two types of prescribed kinematic loadings (displacements on side V) will be considered according to two experimental tests: 1) ux = 1, uy = uz = 0 and 2) uy = 1, ux — uz = 0. As was shown before the stiffness matrix of a viscoelastic homogeneous system based on application of 3-D elements can be presented in the following form: K = XK, + 2GK2 Chapter 3. Steady-state solutions for discrete viscoelastic systems 46 a) 1st loading y s s / b) 2nd loading Figure 3.6: 3-D finite element model of the specimen with two types of loading For a static loading we have [\KX + 2GK2]U = F where U = vector of displacements, F = vector of forces. Neglecting the mass of the specimen in comparison with the attached rigid mass, one can write the following expression for a harmonic loading with frequency u>, namely, [A.(w)J»:i + 2G.{w)K2}U0el"t = F0e (3.8) where the force vector on the right-hand side represents the action from the attached rigid mass on the specimen. Note that force vector F0 will have only non-zero components which correspond to degrees of freedom on the side T. Also note that in the vector of displacements U0 the components which correspond to degrees of freedom on the side T are known (they are prescribed according to the two types of loading). Chapter 3. Steady-state solutions for discrete viscoelastic systems 47 Assume that the size of matrices K\ and K2 is N, and the number of degrees of freedom pertaining to the side T is 3ra, then denote the rest N — 3n as 3m degrees of freedom which correspond to internal degrees of freedom. It assumed here that the elements have 3 linear degrees of freedom per node. Considering the 1st loading, the vector of displacements can be subdivided as follows: - i T r Uo = I 1 0 0 ... 1 0 0 I ulx uly ulz ... u3rnx u3my u3mz = UQJ Analogously one can represent the force vector F0 = fix fly flz ••• / 3 n x f3ny f. 3 n z 0 0 0 ... 0 0 0 Renumerating the degrees of freedom of the system and eliminating the time-dependent coefficient on both sides, one can represent equation (3.8) in the following form: / 1 K r i \ For A» + 2G. — — — — = — V K n . 1 Kh 1 Uoi 0 where the subvector of displacements U0i (size 3m) is unknown, and subvector of forces For (size 3n) is also unknown. Thus one can solve the following system of linear equations -I + 2G*Krj For —\tK^TUor — 2G+KYT-UOT 0 + 2G*K2n Uoi . —\+K}TUor — 2GtKj^Uor and subvectors U0i and F 0r can be determined. Then one can find the total resulting force which is neccessary to apply to create the prescribed vector of displacements on T, namely, for the 1st loading it will be Fx = [l 0 0 ... 1 0 0 ] ^ , ! ^ ^ + A , t f ^ / + 2G,^ (3.9) Chapter 3. Steady-state solutions for discrete viscoelastic systems 48 Having the non-zero components of U0r equal to 1, from (3.9) one obtains the complex stiffness Sx: Sx = A,(Ax + A2) + 2G.{BL + B2) where Ai = RK{\TUQT A2 = RK^Uoi £ x = RK^RUOR B2 = RK^UQ! (3.10) and R = [1 0 0 . . . 1 0 0] Introducing additional notations C 1 = A 1 + A2 C2 = BL + B2 (3.11) the previous expression for the complex stiffness in the 1st type loading becomes: Sx — A , C i + 2G*C2 For the 2nd type of loading the displacement vector will look like U0 0 1 0 ... 0 1 0 I ^ l x ^ l y ^ l z ^ 3 m x ^ 3 m y ^ 3 m z and the auxiliary vector R used in (3.10) will be R = [0 1 0 ... 0 1 0] Applying an analogous procedure one can obtain the complex stiffness for the 2nd type of loading (Fig. 3.6b): Sy = 2{\.D1 + 2G.D2) where coefficient 2 stems from the fact that in the experiment two specimens are used to sup-port the mass m 2 . Note that these complex stiffnesses depend on coefficients C i , C2,Di,D2 Chapter 3. Steady-state solutions for discrete viscoelastic systems 49 (3.11) which in turn depend on the distribution of internal displacements U0i (3.10), which in turn depends on \*,G*. Thus the dependence will be nonlinear, namely, S»(A.,G.) = A.Ci(A„G.) + 2G.G2(A.,G.) 5 t f(A„ G.) = 2(\.D1(\., G.) + 2G.ZMA., <?.)) Now consider the base excitation test (Fig. 3.6a) with a certain frequency u> which will correspond to the 1st type of loading. Assuming that the mass at the top is mi one can write the expressions for transfer functions analogous to (3.7): j/oi A»Ci(A», G*) + 2G„,C2(A», G») T l ( a ; ) xoi -m^ 2 + A.Gi(A.,G.) + 2G,G2(A.,G.) and analogous for the 2nd type of loading (Fig. 3.6b) 2(A,£>1(A.,G.) + 2G.ZMA.,G.)) (3.12) T2(UJ) = — = x 0 2 -m2Lo* + 2(A.D1(A.,G.) + 2G.D2(A.,G.)) One can rewrite these previous nonlinear equations in the following form (3.13) d(A. ,a) (Ti - -1) G 2(A.,G,)(Zi- 1) " ' A. 2Dx(\.,G.)(T2 -1) 2D2(A., G.)(T2--1). 2G. T2m2u 2 H Introduce notations for vector Z — [A»,2G»]T and matrix Gi(A.,G.)(Ti - 1) G2(A.,G.)(Ti - 1) 2Ui (A„G. ) (T 2 - l ) 2£> 2 (A.,G0(T 2 -1) If Tx, T2 are known and A.i, G,i are assumed one may write Txmxu 2 Zi+1 = [H{Zi)}-1 i = 1,2,3,, (3.14) T2m2u}2 The iterations are repeated until convergence is reached \ \Zi+1 — ZA\ < e, (where e is a smal appropriate quantity), or another equivalent criteria can be used, e.g., \\Tle-T1\\<e and | |T 2 c -T 2 | |<e Chapter 3. Steady-state solutions for discrete viscoelastic systems 50 Z0 = [H(Zo)}-1 where Tlc,T2c = calculated values of the transfer function (3.12),(3.13), and TltT2 = exper-imental values of the transfer function in the 1st test, and 2nd test respectively. In other words a fixed point Z0 of the mapping (3.14) is determined Txmxu}2 T2m2u2 The convergence of this iterative process can be analyzed from a mathematical point of view using the contraction mapping principle (see for example [72]), this analysis is omitted here. The transfer functions Tx, T2 (available from experiment [73]) for such material as PECH (polyepichlorohydrin) are presented as functions of frequency in Appendix A, Fig. A . l for the test 1, and in Fig. A.2 for the 2nd test. For test 1 the value of the top mass was mi = 1.91 kg, for the 2nd test it was m2 = 0.54 kg. The dimensions of the specimen were 0.03 x 0.03 x 0.03 m. The temperature of all specimens in all tests was 20° C. Analogous transfer functions for such materials as EAR C - 1002 are presented in Fig. A.3.A.4, DPNR (natural rubber) (Fig. A.5.A.6), CR (chloroprene) (Fig. A.7.A.8) and NBR (acrylonitrile-butadiene rubber) in Fig. A.9,A.10. These transfer functions Tx{w), T2(w) were utilized in equation (3.14). The criteria of convergence for the iterative procedure (3.14) was l ^ i c - T i l , \T2c-T2\ ^ — < re and — < re \J-l\ \ 12\ where TXc,T2c are calculated values, and the relative error re was assumed 0.5 %. Having determined the complex moduli A*, G* one can determine E* and according to the following formulae: j2 _ 3A* + 2<j* v _ E* ^ A* + Gt 2G* Convergence of the procedure (3.14) was observed for all cases (no matter what the starting point was), though the number of iterations required to reach the prescribed level of Chapter 3. Steady-state solutions for discrete viscoelastic systems 51 re was different for different materials and frequencies. The fastest convergence was observed in the case of EAR C — 1002 material, which yielded values of Poisson's ratio between 0.35 - 0.2 (real part) in the considered frequency range. In Appendix A the values of complex Young's modulus, Poisson's ratio, shear modulus are presented which were obtained using the above procedure for five materials PECH, EAR C - 1002, DPNR, CR, NBR (Fig. A.11-A.25). For EAR C - 1002 the results are compared with results from [2] (Fig. A.16), where the real part of the complex shear modulus Gi (storage modulus), and the loss factor ^ are shown as functions of frequency. One can notice a good agreement between these results. The values from [2] correspond to the temperature 24° C, and the experimental results in this study are presented for the temperature 20° C. As a conclusion one can say that in the considered frequency range 10 — 250 Hz for such materials as CR,NBR complex Poisson's ratio can be considered as real and constant, because its values fluctuate between 0.4995 and 0.4991, and it has a negligible imaginary part. The same conclusion is not true for such materials as PECH, EAR, DPNR which demonstrate a dependence of v* on frequency in the range of 10 — 250 Hz. At low frequencies the measured transfer functions is close to 1 for all materials (Fig. A.1-A.10). Thus for frequencies below about 10 Hz additional experiments with different geometric parameters of the specimens are required to provide transfer function values which are not too close to 1 to cause any accuracy distortion. One can also see that the value of the transfer functions tends to zero (real and imaginary part) as the frequency increases (Fig. A.1-A.10). Therefore to use its values for higher frequencies does not seem appropriate due to a possible accuracy distortion of experimental data. The behaviour of E, for a greater frequency range > 250 Hz requires additional experiments with different geometric parameters of the specimens to provide transfer function values which are not too small to cause any accuracy distortion. Chapter 3. Steady-state solutions for discrete viscoelastic systems 52 3.3 Inhomogeneous material systems Consider now a system consisting of a number of subsystems with different material properties. The stiffness matrix of the system can be presented in this case as N K = E ( A i ^ i i + *GjK2i) where N is the number of homogeneous subsystems. The equation of motion (3.5) can be written as N ( - ^ M + E ( A * ^ i i + 2G.jK2j))X0ei"t = F0e Therefore the steady-state solution (in real form) is X(t) = Re N [-u>2M + + 2G,jK2j)]-1F0e iu>t Extension to the general case of periodic loading with a period T = ^ is accomplished by using the complex Fourier series fc=o Using the linear superposition principle one can find the solution in the following form X(t) = Re N fe=o i=i 3.3.1 Numerical results and comparison with experimental data A program was written to treat the case of inhomogeneous material viscoelastic systems and some numerical results are demonstrated below. A vibration rig (Fig. 3.7) was used as an example of an inhomogeneous viscoelastic system. The steel box (mass 34.2 kg) was mounted on polymeric pads which in turn were placed on two cantilevered steel beams. The dimensions of the vibration rig are presented in Appendix B. Chapter 3. Steady-state solutions for discrete viscoelastic systems 53 5 Y X Figure 3.7: Vibration rig with viscoelastic mounts Four pads (mounts) which are indicated in Fig. 3.7 by 1-1',2-2',3-3',4-4' have a cubic shape 3 x 3 x 3 (cm) and such polymeric materials as DPNR and EAR C — 1002 were prescribed for these mounts. Note that the complex moduli for these materials (which were prescribed to the mounts in the numerical model) are presented in Appendix A. An electromagnetic shaker provided a sinusoidal (of frequency a>) vertical force applied at point A. The following quantity (which will be further called the response function) was measured in the experiments [73]: where | | V | | is the amplitude of the velocity of point 5 or 6, and is the amplitude of the external harmonic force at point A (Fig. 3.7). The comparison of experimental results and numerical ones is produced in terms of this response function. The numerical results are represented by a series of calculated steady-state responses, each of them yields a point in the graph (Fig. 3.8). The experimental results were produced by using a frequency sweep (0.1 Hz/s) with the measurement of the response function for each (3.15) Chapter 3. Steady-state solutions for discrete viscoelastic systems 54 Figure 3.8: Response function (point 5). Material - DPNR. — Numerical, Experi-mental passing frequency. The experimental response function was measured on a set of discrete points in the frequency domain with a frequency increment 0.25 Hz. A frequency increment used in the numerical model was 0.1 Hz. For such material of the mounts as natural rubber (DPNR) the results are shown in Fig. 3.8 (for point 5) and Fig. 3.9 (for point 6), where the solid line represents the numerical results in terms of the response function (3.15), and the dashed line - experimental results. For such material of the mounts as EAR C — 1002 analogous results are presented in Fig. 3.10, 3.11 for points 5 and 6 respectively. One can see a good correspondence of the numerical results with the experimental ones for both materials. Some disagreement can be explained by the following facts. Firstly, the frequency sweep rate used in the experiments was 0.1 Hz/s, which is not slow enough for the assumption to be made that at each passing frequency a steady-state response takes place. Recall that numerical graphs in Fig. 3.8-3.11 are obtained by evaluation of the ratio of the steady-state response velocity amplitude to the force amplitude at each frequency. Secondly, Chapter 3. Steady-state solutions for discrete viscoelastic systems 55 0.012 Frequency, Hz Figure 3.9: Response function (point 6). Material - DPNR. mental Numerical, Experi-Frequency, Hz Figure 3.10: Response function (point 5). Material - E A R C-1002. — Numerical, Experimental Chapter 3. Steady-state solutions for discrete viscoelastic systems 56 Figure 3.11: Response function (point 6). Material - E A R C-1002. — Numerical, Experimental the values of complex moduli used in the numerical model may differ to some extent from the real values. Recall that these complex moduli (see Appendix A) were calculated from certain experimental input data according to the procedure presented in section 3.2. Thirdly, the frequency increment in the experiments was 0.25 Hz (in the numerical modelling it was 0.1 Hz), therefore some points in the frequency domain near peaks are not presented in the experimental graphs. Chapter 4 Closed form solutions for discrete viscoelastic systems Two methods of obtaining homogeneous solutions for the vibration of discrete viscoelastic systems are compared. The first [54] uses the Laplace transform; the second, a new method [56], [55] which is given a name "substitution method", does not. The advantage of the substitution method is indicated for the case of multi-degree-of-freedom (MDOF) systems. Analytical results computed by using the substitution method are demonstrated on an example of a viscoelastic beam. The free and forced vibration responses with various initial conditions are presented. 4.1 Single-degree-of-freedom system Consider a mass m connected to the base by a one-dimensional massless element which combines the properties of hereditary and viscosity in parallel. The equation of unforced motion is as follows mx (4.1) where the relaxation kernel is assumed to be of the form: n -<*i{t-r) i = l The initial conditions are: x(0) = x0 x(0) = x0 57 Chapter 4. Closed form solutions for discrete viscoelastic systems 58 The presence of the term cx arises from placing a viscous dashpot in parallel with a viscoelastic element. It should be noted that for some elastomeric materials the term cx can be introduced into the constitutive law of the material. Here the integral term plus the instantaneous part kx(t) will be referred as a viscoelastic element, and term cx as a viscous element. 4.1.1 Application of the Laplace transform Applying the Laplace transform to (4.1), one obtains: n ^_ m(p2x — x(0) — px(0)) + c(px — x(0)) + kx — fciV — = 0 where p is the Laplace transform variable. From the expression above the Laplace transform of the x-displacement will be m(x(0) +px(0)) + ca(0) _ R(p) X ~ p>TO + p c + * ( l - 5 - i ; j S L _ ) - T(p) Knowing the roots, in general complex, of denominator T(p) (denote them pi), and assuming they are simple (multiplicity of 1) one can obtain a solution in the form (see the theorem in [74]): Thus the free vibration solution in complex form is given by (4.2). If the initial conditions x0 and xo are assumed real, then the solution (4.2) will be real, because the roots pi are either complex conjugate, or real. In general one can assume that x0 and x0 are complex, then the real part of the complex solution (4.2) will be a solution for the real parts of x0, x0, and imaginary part of (4.2) will be a solution for the imaginary parts of x0, x0. 4.1.2 The substitution method The solution of (4.1) will be sought in the following form [56], [55]: n+2 x(t) = £ c,e*" (4.3) Chapter 4. Closed form solutions for discrete viscoelastic systems 59 where Cj, pj are complex. Substituting it in (4.1) one obtains n+2 n n n+2 E W + cpj + fe-fcE drH c; e P y t + * EE _ , a. j=l i=l "•" "» t=l j=i Po ^ "* Cj]e-ait = 0 In this equation there is a sum of exponential functions. To satisfy this equation one can set all the coefficients of these functions to zero, namely: mp) + cpi + fc(l-E ) = ° i=i Pj i a* and The initial conditions: n+2 E ai ±i Pj + a •CJ = 0 n+2 E CJ = X° i = l,n n+2 E CjPj = x0 (4.4) (4.5) (4.6) provide two more equations. Equation (4.4) can be called as the characteristic equation with respect to n + 2 complex (in general) roots pj. They can be determined, for example, by using Newton's method, or by reduction to an eigenvalue problem (see the Remark below). Equations (4.5), (4.6) yield a linear system of n + 2 equations with n + 2 unknown complex constants Cj. In matrix form one obtains: <*1 Pl+<*1 a i P 2 + a i ai Pn+2+Ol On On a„ Pl+°tn P2+an 1 1 1 Pi P2 Pn+2 0 Cl ... ... 0 Cn+2 x0 XQ (4.7) Note that ai ^ 0 for i = l,n, otherwise it does not make sense to introduce a term in the constitutive law which is a priori 0. It is interesting to note that the requirement that roots Pj are simple (multiplicity of 1) is necessary here to provide the nonsingularity of the matrix in Chapter 4. Closed form solutions for discrete viscoelastic systems 60 (4.7). It may be noted that if there is a repeated root then the corresponding columns in (4.7) will be identical and the determinant of the matrix vanishes. This simple roots requirement is not a restriction imposed on a system, it is just an indicator which can be used to determine if a given viscoelastic system (4.1) assumes a solution in the form (4.3), or not. If it assumes, then according to the theorem of uniqueness [53] this solution is unique. The satisfaction of this requirement is expected for most viscoelastic systems, although if this requirement is not satisfied, then the analytical homogeneous solution may differ from (4.3), and this requires an additional investigation. Remark. The reduction of the characteristic equation (4.4) to an eigenvalue problem is as follows. Equation (4.4) can be rewritten as pn+2bn+2 +Pn+1bn+1 + ... +p6i + b0 n?=1(p + ai) = 0 or Pn+2bn+2 + pn+1bn+1 + ... +pb1 + bo = 0 (4.8) where coefficient bn+2 = m, and the other coefficients can be easily evaluated (expressions for them are omitted here). Introducing the state-space vector Q = [q0 pqo p2qo ••• Pn+1<lo]T, the following eigenprob-lem will correspond to the characteristic equation (4.8): P V bi b2 63 1 0 0 0 1 0 0 0 1 bn+2 0 0 0 0 0 or in abbreviated form \ + 60 0 0 0 - 1 0 0 0 - 1 0 0 0 0 0 0 Q = 0 0 (pA + B)Q = 0 (4.9) Chapter 4. Closed form solutions for discrete viscoelastic systems 61 The eigenvalues of (4.9) will be the characteristic roots of (4.4). Example. For the purpose of illustration, consider the equation with initial conditions: a:(0) = 1 + iO, x(0) = 0 + iO. The complex solutions obtained by application of the Laplace transform (4.2) and by using the substitution method (4.9), (4.7), (4.3) are identical: x(t) = (0.6184 - i0.03739)ePlt + (0.6184 + i0.03739)ep2' - 0.2368eP3< where the roots of the characteristic equation are Pi = -0.1761 + i0.86071 p2 = -0.1761 - iO.86071 p3 = -0.64779 whence the free vibration response comprises a combination of a decaying oscillatory mode and one overdamped mode. 4.2 Multi-degree-of-freedom system. Poisson's ratio is constant. Application of the finite element method to an elastic system yields the mass matrix M and the stiffness matrix K. Represent the stiffness matrix as It is assumed here that Poisson's ratio operator is elastic (i.e. constant) v = v = const. The approach when Poisson's ratio operator is not constant will be presented in section 4.5. K = EK0 (4.10) Young's modulus operator E has form: (4.11) Chapter 4. Closed form solutions for discrete viscoelastic systems 62 where the relaxation kernel is assumed as a sum of exponentials: r ( * - r ) = t ^ " Q i M (4.12) t=i The equation of free motion can be written as MX + CX + EK0X = 0 (4.13) The term CX takes into account the presence of viscous damping arbitrary distributed over the system. 4.2.1 Laplace transform method Applying the Laplace transform to (4.13) and taking into account (4.11), (4.10), (4.12) one obtains: n „ ai P * M + P C + K ( I - J : — ^ ) X = M[X{0)+pX(0)) + CX{0) In abbreviated form, denoting the matrix coefficient of X as S: S{P)X = M{X(0) + PX{0)) + CX(0) (4.14) where the matrix S(p) can be written as Sip) = ^—r1 ;D(P) The elements of matrix D(p) are polynomials of degree n + 2. Inverting the matrix S(p) S-1(p) = II»=1(P + ai)D-1(p) and introducing the adjoint matrix A(p) = D~l(p)det(D(p)), equation (4.14) can be written as follows X = UUiP + ^ O d e f f p ) ) [ M ( x ( o ) + PX(0)) + CX(0)} Chapter 4. Closed form solutions for discrete viscoelastic systems 63 The roots of the characteristic equation det(D(p)) = 0 can be determined numerically in an analogous way to that shown in (4.9). The extension of (4.9) to the M D O F case will be shown in the following section. Denoting the characteristic roots as pk, the determinant of D(p) will be det(D(p)) = gIl»=1(P - Pk) where N = m(n + 2), m is the number of degrees of freedom, n = number of terms in (4.12), and g = a constant complex coefficient. Now one can express the solution of the free vibration problem in a form analogous to (4.2): = E n * = i ( f t + « i ) T T * ^ -AM(X(0)+PiX(0)) + CX(0)}e^ (4.15) This is the extension of the theorem from [74] (the scalar case was considered there) to the matrix case. It should be noted that a formula similar to (4.15), but for the case of an elastic undamped system is mentioned in [9]. Remark. Note that in the use of formula (4.15), it is not an easy task to derive the analytical expressions for D_1{p) and consequently for A(p). Also the computation of A(pi) by A(Pi) = D~l (pi)det(D(Pi)) requires the calculation of the limit A(Pi)=limi[D-1(p)det(D(p))] where limp_^.p. ||D - 1(/>)|| = oo and limp^.Pi det(D(p)) = 0. The substitution method described below, yields a better way to obtain the closed form solution. 4.2.2 The substitution method The solution of (4.13) is sought in the form [56], [55]: m(n+2) X(t)= £ CiXjf* (4.16) Chapter 4. Closed form solutions for discrete viscoelastic systems 64 where Xj = a complex vector (m x 1), Cj, pj are complex, m=number of degrees of freedom (size of matrices, M, C, K). Substituting (4.16) into (4.13) and taking into account (4.11), (4.10), (4.12) one obtains: n(n+2) E 3=1 ptM + PiC + Kil-Z—^;) i=l f3 ' "* , „ m(n+2) a; i=i j=i Pj + <*i c0Xj\e-ait = 0 Therefore the following equations need to be satisfied: and (with a nonsingular matrix K): m(n+2) Xj = 0 j = l,m(n + 2) (4.17) E ai -CjXj = 0 i = l,n i=\ Pj + <*i If matrix K is singular then instead of (4.18) one should consider m(n+2) (4.18) « E ai -CjXj = 0 i = l,n j7i Pi + ai Equation (4.17) represents a so-called characteristic equation, which can be reduced to the eigenvalue problem (see the Remark below) and the characteristic roots (eigenvalues) pj along with vectors Xj can be determined. The initial conditions of the problem are: a(n+2) E C3^3 = Xo i=i m(n+2) E C3P3X3 = X 0 (4.19) Note that initial conditions (4.19) can be replaced by other conditions in the following ways: i) E CjXj = X0 X; CjXj^=X(T) (4.20) m(n+2) E C 3 3 ~ X * 3=1 or ii) m(n+2) E C 3 X 3 ~ Xo ^(n+2) E CjPjXj<?>T = X(T) Chapter 4. Closed form solutions for discrete viscoelastic systems 65 or or or m(n+2) iii) Y cjPixi = i=i m(n+2) iv) coPoX3 = Xo v) m(n+2) E c i X i f T = X{T) 3=1 m(n+2) E C J A > ^ = X(T) 3=1 m(n+2) 3=1 m(n+2) E c ^ e ^ T = X(T) 3=1 where T represents some given instant of time (T > 0). Thus instead of the initial value problem, the terminal value problem or its variations can be posed, and this will only require a change in the two last matrix rows in the system (4.21) accordingly to conditions i) — v). Proceeding with conditions (4.19), relations (4.18), (4.19) constitute a system of linear m(n + 2) equations with respect to m(n + 2) unknowns Cj. In matrix form this system can be written as P I +<*! P2+C1 X\ X2 PiXx p2X2 P m ( n + 2 ) + « l %(n+2) P m ( n a + 2 ) + < * r > X m ( n + 2 ) (n+2) cm(n+2) J 0 0 Xo Xo (4.21) Pm(n+2)Xm(n+2) Remark. Introducing a common denominator for all terms, equation (4.17) can be rewritten as follows: [n?=1(p; + aAtfM + PjC + K) - o&Liwbi + = 0 i=l Collecting the matrix coefficients of p" + 2 ,P j + 1 , and denoting them as B n + 2 , B „ + i , B 0 Chapter 4. Closed form solutions for discrete viscoelastic systems 66 respectively, equation (4.17) is equivalent to: Br B2 B3 ••• Bn+2 Bo 0 0 .. . 0 I 0 0 ... 0 0 - / 0 .. . 0 p 0 / 0 ... 0 + 0 0 - / .. . 0 \ 0 0 • • • / 0 0 0 0 .. . - / • \ xj 0 Pjxj 0 PjXj 0 .) 0 or in abbreviated form where j th eigenvector will be (pA + D)Q = 0 Qj j = l,m(ra + 2) (4.22) (4.23) PJXJ P)xJ Pj+'XJ The matrix coefficient B n + 2 = M, and other coefficients can be easily evaluated. The eigenvalues of (4.23) will be the characteristic roots of (4.17). The basic part Xj of eigenvector Qj is used in (4.21), (4.16). Here it is assumed that coefficients ^ 0 for i = l,n (otherwise the determinant of the matrix in (4.21) would be 0). As far as the eigenvalues pj are concerned, they are allowed to be multiple (for the M D O F system case only). However they must have linearly independent corresponding eigenvectors, to provide nonsingularity of matrix in (4.21). It is easy to note that if there is a root, for example of multiplicity 2, which has only one eigenvector (the second one is linearly dependent), then the two corresponding columns in (4.21) will be linearly dependent and the determinant of the matrix vanishes. In other words, in order to Chapter 4. Closed form solutions for discrete viscoelastic systems 67 have a solution of (4.13) in the form (4.16), it is required that the eigenproblem (4.22) yields m(n + 2) linearly independent eigenvectors. It should be noted that the size of the eigenvalue problem (4.22) can be quite large, so it is necessary to have an effective eigensolver for matrices of the type in (4.22). 4.3 Periodic loading case. Application of the substitution method The general case of periodic loading with a period T — ^ is treated by using a complex Fourier series. The forcing function is represented as F(t) = J2Fkei"kt where Fk are vectors of size m. The equation of motion is written as follows: u>k = hu MX + CX + EK0X = E Fke k=o (4.24) with initial conditions X{0) = Xo X{0) = Xo The general solution is now sought as a sum of homogeneous and particular solutions, namely: (4.25) Substituting (4.25) in (4.24) and taking into account (4.11), (4.10), (4.12) one obtains: Xj(?i*+ m(n+2) X(t)= £ CiXif* + £ Zke^ 3=1 L E fc=o i(n+2) E 3=1 p'M + pjC + Kil-Y,-^) i = l + 3 1 . L + E k=0 -UJ, \M + iu3kC + K(l - £ -—;—) Chapter 4. Closed form solutions for discrete viscoelastic systems 68 n m(n+2) +*£[ E 0; fc=o From this the characteristic equation (4.17), and subsequently the eigenproblem (4.22) follows yielding values for pj and Xj. The quantities Zk are determined as follows: - l -LO\M + iwkC + K(1-Y ; ) k = 0,L The linear system of equations (4.21) with respect to unknowns Cj in this case will have the modified right-hand side: P l + « 1 1 P2+<*1 * -^-xL -^-x2 X l X 2 P l - X l P2^T 2 4.4 Numerical results P m ( n + 2 ) + « ! X m ( n + 2 ) P m ( n + 2 ) + « « m ( n + 2 ) -^m(n+2) Pm(n+2)-^m(n+2) L Cm(n+2) J XQ — Efe=0 ^fc A program was written which calculates analytical solutions according to the substitution method and was used for the examples below. Numerical results are presented for a viscoelastic beam of square cross section with fixed ends (Fig. 4.1). The parameters of the beam's cross section were 0.01 x 0.01 m, the length = 0.12 m, instantaneous Young's modulus E was 0.15e+08 Pa, v = 0.3, the density of the material was 0.141e+04 kg/mz. The beam was meshed by 6 general finite beam elements. Each node of a beam element had 6 degrees of freedom (3 linear and 3 rotational). Thus the size of the problem (number of degrees of freedom) was m = 30. Chapter 4. Closed form solutions for discrete viscoelastic systems 69 The relaxation kernel in (4.11) was taken as T(t -T) = a i e- a i(<-T) (4.26) where a\ = 150, c*i = 200 (n = 1). Thus the size of eigenvalue problem (4.22) and of linear system (4.21) was m(n + 2) = 90. The matrix coefficients in (4.22) for this case are computed as follows: B0 = K(a1-a1) BX = Cax + K B2 = Max + C BZ = M One can see that if ai is negligible [ax « 0), then the characteristic matrix equation (4.17) is reduced to the usual complex eigenvalue problem of viscously damped system. To investigate the influence of the hereditary part (term EKQX) in (4.13), it was assumed that the damping matrix C = 0. Figure 4.1: Viscoelastic beam with fixed ends A standard subroutine "DREIGN" [75] was used for the eigenproblem (4.22) and sub-routine "CDSOLN" for (4.21). Consider the eigenvalue problem (4.22) in more detail. For the system in Fig. 4.1 the chart of computed eigenvalues is presented in Fig. 4.2 (imaginary parts) and in Fig. 4.3 (real parts). Numeration of the eigenvalues was done on the basis of their absolute values. Chapter 4. Closed form solutions for discrete viscoelastic systems 70 Eigenvalue number Figure 4.2: Imaginary part of eigenvalues It is found that for the model (4.26) defined by only one exponential term there are 30 repeated (note this number is equal to the number of degrees of freedom) real eigenvalues which correspond to "overdamped" eigenvectors (also real). The other 30 pairs of com-plex conjugate eigenvalues which have nonzero imaginary part correspond to "underdamped" eigenvectors (complex conjugate). Although the terms "overdamped", "underdamped" are usually used in the context of viscously damped system, their use in the case of a viscoelastic (hereditary) system may be adopted as well. For the purpose of illustration, some of the eigenvectors are presented below in Figs. 4.4-4.7. The 1st of the underdamped eigenvectors (corresponding to eigenvalue No. 31) has a shape of the 1st vibration mode in plane X - Y (Fig. 4.4). For eigenvalue No. 29 (note that numeration of repeated real eigenvalues is arbitrary to some extent) the eigenvector has the same shape (Fig. 4.5) as the eigenvector corresponding to eigenvalue No. 31. A second vibration mode in plane X — Y is represented by the underdamped eigenvector No. 35 (Fig. 4.6). The analog to it from the spectrum of the overdamped eigenvectors (corresponding to eigenvalue No. 28) is in Fig. 4.7. Chapter 4. Closed form solutions for discrete viscoelastic systems 71 CO CD I 1 1 1 ' I 1 1 1 1 I 1 1 1 1 I 1 1 1 1 I 1 1 1 ' I 1 1 1 1 I 1 1 1 1 I -20 •c co 10 20 30 40 50 60 70 80 90 Eigenvalue number Figure 4.3: Real part of eigenvalues Figure 4.4: The 31st eigenvector (1st underdamped) Chapter 4. Closed form solutions for discrete viscoelastic systems 72 5.00E-5 Node Figure 4.5: The 29th eigenvector (analog to the 31st one) Node Figure 4.6: The 2nd underdamped eigenvector Chapter 4. Closed form solutions for discrete viscoelastic systems 73 1.00E-4 1 2 3 4 5 6 7 Node Figure 4.7: Overdamped eigenvector (analog to the 2nd underdamped) The numerical results in terms of displacements are presented in Figs. 4.8 and 4.9. The vertical displacements and velocities were imparted to the middle node 4 at t = 0. The following variations of initial conditions were considered for node 4: case i) y(0) = 0.005 m, y(0) = 0; case ii) y(0) = 0, y(0) = 10 m/s. All the other degrees of freedom had zero initial conditions. The vertical displacements of the node 4 (Fig. 4.1) were chosen to illustrate the response. The results of free vibration response are presented in Fig. 4.8 for the case i), Fig. 4.9 (case ii). The contribution of vectors CjXj to the solution, namely, their norms | | C J X , - | | are presented in Fig. 4.10, 4.11 for the initial conditions i) and ii) respectively. Note that the complex conjugate eigenvectors contribute equally, for example, pairs of eigenvectors No. 31 and 32, 43 and 44, 59 and 60, 71 and 72, 85 and 86 in Fig. 4.10. One can see that the contribution of overdamped eigenvectors No. 24, and 29 is quite noticeable for the case of initial conditions ii) (Fig. 4.11). It is clear that relative contribution of the underdamped eigenvectors (with respect to each other) will not change in time, because the real parts of their eigenvalues are Chapter 4. Closed form solutions for discrete viscoelastic systems 74 0.000 0.005 0.010 0.015 Time, s 0.020 0.025 Figure 4.8: Free vibration response, case i) Figure 4.9: Free vibration response, case ii) Chapter 4. Closed form solutions for discrete viscoelastic systems 75 i . . . . i . . . .mill. . . i .III. • i i i i .III. . . . .Ill 111,1 10 20 30 40 50 60 70 80 90 Eigenvector number Figure 4.10: Contribution of eigenvectors, case i) the same (Fig. 4.3). One can see that the contribution of the overdamped eigenvectors No. 24 and 29 is quite noticeable for the case of initial conditions ii) (Fig. 4.11). Note that their relative contribution with respect to the underdamped terms will increase in time, because real parts of their eigenvectors are greater (less negative) (Fig. 4.3). As it should be (one can see it from system (4.21)) this contribution depends on initial conditions. A boundary value problem (4.20) was then considered. The conditions shown in Fig. 4.12 were assumed, where for a chosen instant of time T = 0.02 s , all the displacements were prescribed to be zero. The solution (the vertical displacements of nodes 2,3,4) for this case is shown in Fig. 4.13. Then an instant T = 0.016 s was chosen, and the response is shown in Fig. 4.14. One can observe the difference in the initial velocities which are not prescribed, but computed. For the forced vibration a vertical periodic force F(t) (Fig. 4.15) is applied at the node 4 (Fig. 4.1). This loading has a period T — 0.1 s and the Fourier series (30 harmonics) 0.070 0.060 h _ 0.050 : <T ,j- 0.040 [-0.030 0.020 0.010 0.000 Chapter 4. Closed form solutions for discrete viscoelastic systems 76 0.08 h 10 20 30 40 50 60 70 80 Eigenvector number Figure 4.11: Contribution of eigenvectors, case ii) 5.00E-5 > -CO X < 0.00E0I 3 4 5 Node, axis X Figure 4.12: Initial and terminal shapes Chapter 4. Closed form solutions for discrete viscoelastic systems 77 Figure 4.13: Solution for conditions (4.20), T=0.02 s Figure 4.14: Solution for conditions (4.20), T=0.016 s Chapter 4. Closed form solutions for discrete viscoelastic systems 78 0.04 0.06 Time, s Figure 4.15: Forcing function 0.0020 CD i« o « CD-GO 0.00 0.02 0.04 0.06 Time, s 0.08 0.10 Figure 4.16: Forced response with zero initial conditions Chapter 4. Closed form solutions for discrete viscoelastic systems 79 was used to represent the forcing function F(t). In Fig. 4.15 this function is shown for the duration of one period t C [0,T]. The solution (4.25) with zero initial conditions (for all degrees of freedom) is presented in Fig. 4.16 for t C [0,T], where the vertical displacements of the nodes 2,3,4 (Fig. 4.1) were chosen to illustrate the response. 4.5 The substitution method. Poisson's ratio as a viscoelastic operator The finite element method application yields a mass matrix of the system M and a stiffness matrix K. The derivation of the global system stiffness matrix starts with an element stiffness matrix, where the constitutive law is involved. For an elastic material: [<r] = [E][e] where [IT] = [ ( T n C T 2 2 0 3 3 C T 1 2 <T 1 3 <T23]T, X + 2G X X X + 2G where [E} = X 0 0 0 A 0 0 0 A = uE r l r i3 eJ = L E H E 2 2 633 612 6l3 ^23. A 0 A 0 X + 2G 0 0 0 0 0 0 0 0 0 0 2G 2G 0 0 0 2G 0 0 0 2G E {l + u){l-2u) ~ l + v and E=Young's modulus, i/=Poisson's ratio. The hereditary analog of this matrix is obtained by replacement of E, u by corresponding operators, i.e., E —> E, v —> u, or by another replacement A —>• A, G —> G. Chapter 4. Closed form solutions for discrete viscoelastic systems 80 Viscoelastic Lame's operators A and G are assumed in the following form ft n Xy(t) = X(y(t) - / V a,iexp[-ai(t - r)]y(r)dT) J o i=i Gy(t) = G(y(t) - f E biexp[-7i(t - T)}y(r)dr) J o i=i Now the equation of free vibrations can be written in the following form MX + CX + (\Ki + 2GK2)X = 0 (4.27) (4.28) (4.29) The application of the substitution method in this case is analogous to one described in section 4.2.2. The solution of (4.29) is sought in the form: m(2n+2) X(t)= E CiXj*** (4.30) where Xj = a complex vector (m x 1), Cj, Pj are complex, m=number of degrees of freedom, n = number of exponential terms in (4.27) and (4.28). Denoting N = m(2n + 2) and substituting (4.30) into (4.29), using (4.27),(4.28) one obtains: N 3=1 p>M + PjC + \ K 1 ( l - £ : ^ ) + 2GK2(l-£ h i tr( Pj + a{ fri PJ + ii CjXj^-T-N + £ Z - ^ - ^ X ^ + Y\*GKt E -^-CjXj]e-^ = 0 r - Oi tri' 'jTiPj + ati-' - " j ^ f t + 7 i Therefore the following equations need to be satisfied: n n 1 P)M + PjC + A ^ ( l - E - f - ) + 2 G K & ~ E T-±-) and tri Pj + a{ jTi Pj + ai triPs+V J CjXj = 0 i = l , n (4.31) (4.32) Chapter 4. Closed form solutions for discrete viscoelastic systems 81 and N —CjXj = 0 i = l,n (4.33) Equation (4.31) represents a so-called characteristic equation, which can be reduced to an eigenvalue problem and the characteristic roots (eigenvalues) p j along with vectors Xj can be determined (see section 4.2.2). Relations (4.32), (4.33) and initial conditions (4.19) constitute a system of linear N = m(2n + 2) equations with respect to N unknowns Cj analogous to (4.21): ~^~X2 P l + a i P 2 + « * l P i + 7 1 * P2+71 P / V + 7 1 ' XN -%-XX - ^ - X 2 P l + a n P 2 + « n P l + 7 n X -J*-x2 P2+7n " -^xN -^S—XN PN+fn " Cl CN 0 X0 (4.34) Xi X2 ... XN PIXI piX2 ... PNXN Thus the homogeneous solution in the form (4.30) can be found for this case when Poisson's ratio is considered as operator. Note that the number of exponential terms in the solution (4.30) is taken for this case as N = m(2n + 2) in comparison with m(n + 2) in (4.16) when only Young's modulus was considered as a viscoelastic operator. Remark. Note that the extention of the substitution method to the case of in homoge-neous material systems is accomplished in an analogous way, i.e., by representing the system stiffness matrix-operator as L K = ^ 2(XjKij + 2GjK2j) i=i where L is the number of homogeneous subsystems. Then the solution is sought in the Chapter 4. Closed form solutions for discrete viscoelastic systems 82 following form m(2Lxn+2) X(t) = £ CjXje^ 3=1 where n = number of exponentials in relaxation kernels of operators A j , Gj, m = number of degrees of freedom (size of matrices M, C, Kij, K2j, j = l,L). 4.6 Application of modal analysis technique 4.6.1 Conditions of diagonalization of two matrices. Analysis of eigenvalues Consider an eigenvalue problem of the following form: (XA + B)x = 0 (4.35) where x=eigenvector, A=eigenvalue. A, B are assumed to be real symmetric matrices (a special case of Hermitian matrices) and one of them, say, A is positive definite. Then according to a theorem of linear algebra, see for example [76] (theorem 12.7), there exists a nonsingular matrix (in our case it can be real) P such that PTAP = I (4.36) and PTBP = A (4.37) where A is a real diagonal matrix. Introducing a change of variable in (4.35): x = Pq and premultipling (4.35) by PT one obtains (XI + A)q = 0 (4.38) Chapter 4. Closed form solutions for discrete viscoelastic systems 83 where relations (4.36), (4.37) were used. From this it follows that the eigenvalues of the problem (4.35) are the elements of —A and therefore all real. Also all eigenvectors x are real. It is obvious that the eigenvectors Q = [ c j i . . . e j n ] of problem (4.38) are linearly independent, hence det{Q] ^ 0. Actually one can put Q = I, which leads to the conclusion that the matrix of eigenvectors X of problem (4.35) is X = [Xl...xn] = PQ = P and these eigenvectors are linearly independent, because matrix P is nonsingular as was mentioned above. Consider now the case when neither of the matrices A, B is positive definite. Sometimes in practical applications, e.g. viscously damped dynamic systems, the matrices A, B can be (so-called state-space form): C M K 0 A = B = M 0 0 -M where M=mass matrix (positive definite), -ftT=stiffness matrix (positive semi definite), C=damping matrix (symmetric). One can see that matrices A and B are not positive definite, because according to a necessary condition [76] of positive definiteness all the diagonal elements must be positive. It is expedient to establish whether A is semidefinite, indefinite (without zero eigenvalues), or negative definite. Represent A in the following form: ' c M J C ' 0 M M 0 0 M J 0 From this representation one obtains: det[A] = {-l)mdet[M]det[M} ^ 0 because M is positive definite. The number m in the formula above is the size of matrices M and C. Therefore A cannot be semidefinite, because in this case at least one of its Chapter 4. Closed form solutions for discrete viscoelastic systems 84 eigenvalues would be 0, and the determinant of A (as the product of all eigenvalues) would be 0 as well. Thus one comes to the conclusion that matrix A is either indefinite (without zero eigenvalues), or negative definite. Consider the following eigenproblem: Ap = rjp or (A — nl)p = 0 where p — an eigenvector, and n is an eigenvalue. It is known that the eigenvalues of a general real symmetric matrix are real and eigenvectors (can be taken real) are linearly independent and orthonormal [76]. Upon the normalization of eigenvectors pfPi = 1, the following holds: pfApi = T)i As it was mentioned above at least one of the eigenvalues rji is negative (matrix A is either indefinite without zero eigenvalues, or negative definite). One can renormalize eigenvectors in P to get PTAP = / , namely, by introducing the new eigenvectors as follows pi —^ •^Pi- Thus it is clear that some of these new eigenvectors P=\pi...pn] corresponding to the negative eigenvalues 77, will be complex. Now introducing a change of variable x = Pq and premultiplying (4.35) by PT one obtains: (AJ + PTBP)q = 0 (4.39) where PTBP will be a complex symmetric matrix (that is not Hermitian). Therefore to say that the eigenvalues and eigenvectors of problem (4.39) (and respectively of problem (4.35)) for the case when neither of the matrices A, B is positive definite are real is not possible. In general they can be complex. The conditions of diagonalization of matrices A, B for this case were discussed, for example, in [77]. Chapter 4. Closed form solutions for discrete viscoelastic systems 85 4.6.2 Conditions of diagonalization of three matrices Consider at first the system of equations of free vibration of a discrete viscously damped system MX + CX + KX = 0 (4.40) where mass matrix M is symmetric positive definite, J^=stiffness matrix (symmetric positive semidefinite), C=damping matrix (symmetric). In the section above it was shown that the diagonalization of two real symmetric matrices (one of them is positive definite) by means of a nonsingular real matrix transformation is always possible. It is known that the mode shape (eigenvectors) matrix of the undamped system can make this transformation with respect to the M and K matrices. In this section a necessary and sufficient condition of diagonalization of C by this matrix of undamped mode shapes is considered. For a general M, K, C system such condition has been formulated, for example, in [78], [13] as M~lC commutes with M~XK, or in equivalent form CM~XK = KM~lC It should be noted that a necessary and sufficient condition for diagonalization of matrix C in terms of undamped system modes $ $ r C$ = [Of], or C = ($ r)" 1[a i]$- 1 (4.41) ( [cii] is a diagonal matrix) implies that N real coefficients a; is enough to describe all the variety of possible damping matrices for the given system defined by M, K. Starting consideration with systems which have distinct (no repeated) undamped natural frequencies the necessary and sufficient condition of proportionality imposed on matrix C will be formulated in terms of matrices M, K in the proposition below. Chapter 4. Closed form solutions for discrete viscoelastic systems 86 Proposition. For the given matrices M,K of order N in (4.40) the matrix C will be diagonalized simulta-neously with M , K if, and only if, it is represented as N C = Mj^gtiM^K)1-1 (4.42) i=i where = 1,2,...,AT are arbitrary coefficients. Proof. Necessity. There is such a matrix $ consisting of undamped mode shapes that <Srtf$ = [bi] $TC$ = [ai\ where [bi] = b, [ai] = a are diagonal matrices and the matrix $ is mass normalized (first equation), so the elements of b will be the natural frequencies squared. Rewrite the above expressions as C = ($ r)- 1a$" 1 (4.43) K = ($T)"16$-1 (4.44) and M = ($$ T ) _ 1 (4.45) or A T 1 = <MT (4.46) It is obvious that there is one to one mapping, namely, a —>•<— C, b —»<— K. Express matrix a as a = g1I + g2b+... + gNbN-1 (4.47) Chapter 4. Closed form solutions for discrete viscoelastic systems 87 Equation (4.47) can be rewritten as a system of linear equations for coefficients gk,k 1,2, ...,N N - l 1 W b\ ... b» 1 b2 bl ... of" 1 1 bN b 2N iN-l bN 9i 92 = Q-2 9N aN The determinant of this system (the Vandermonde determinant) is non-zero because there are no repeated values bi (natural frequencies squared), so the solution is guaranteed and this representation (4.47) holds. Rewrite (4.47) using (4.43),(4.44) N \ fe - l $ r c $ = J2 gk($ TK$) k Therefore any matrix C can be represented as C = ($ r)- 1[5l/ + g2$TK$ + g3$TK$$ TK$ + ...+ +gN$TK$...$ TK$}$- 1 and using (4.46) one can obtain N c^M^miM-'K)1-1 1=1 as required. Sufficiency. The damping matrix is given as in (4.42). Using (4.45),(4.46),(4.44) one can obtain c = ( $ $ r ) - 1 f > ( $ $ r ( $ r ) - 1 0 $ - 1 ) ' - 1 = i=i = ($ T)- 1(9i+92b+... + gNbN- 1)$- 1 Therefore Chapter 4. Closed form solutions for discrete viscoelastic systems 88 1-1 1=1 represents a diagonal matrix as a sum of diagonal matrices, which was required to prove. It may seen from (4.42) that if the first two terms are taken then C = giM + g2K which is Rayleigh damping. If more terms taken, e.g. three, then C = giM + g2K + g3KM~lK and so on. Remark. Note that for systems with repeated undamped natural frequencies expression (4.42) will represent a sufficient condition, but not the necessary one. This means that there might be matrices C, which are not represented as in (4.42), but are diagonalizable. In the case of repeated eigenvalues a formula for general classical damping matrix will be as follows N c = M ^ < / , ( M - 1 i i r ) ' - 1 + * -1=1 0 0 0 0 0 0 0 a; 0 0 cci+i 0 0 a i + s - l 0 0 0 0 0 (4.48) Here it is assumed that jth eigenvector is the j'th column of the matrix $, and the eigenvalues numerated from i to i + s are repeated (multiplicity = 5 + 1). In (4.48) s quantities a;, ai+1, Chapter 4. Closed form solutions for discrete viscoelastic systems 89 ...,a1+9-! are arbitrary real coefficients which are allocated along the principal diagonal in columns i,..., i + s — 1 respectively. Formula (4.48) is the necessary and sufficient condition of diagonalization for the case when repeated undamped natural frequencies are present. In this representation the matrix of mode shapes (eigenvectors) $ is involved, i.e., they need to be calculated. Therefore in this case there is actually no advantage of formula (4.48) in comparison with formula (4.41). 4.6.3 Diagonalization of discrete viscoelastic systems Now one can consider the conditions of transformation of the equation of motion of a discrete viscoelastic system to a uncoupled form. Assume that the system consists of one viscoelastic homogeneous material. If Poisson's ratio of the material is assumed constant, then the stiffness matrix-operator K can be expressed in the following form K = EK0 where Young's modulus operator E is taken out as a common factor for all components of matrix-operator K, and K0 is a symmetric matrix. Therefore the equation of free motion can be written in the following form: MX + CX + EK0X = 0 (4.49) with the initial conditions X{0) = Xo X{0) = Xo Matrix C (viscous damping) is assumed here to be diagonalizable together with matrices M and Ko. Therefore introducing the change of variable X = $Q Chapter 4. Closed form solutions for discrete viscoelastic systems 90 and premultiplying (4.49) by $ r , the decoupled system of equations is obtained: rmqi + Ciqi + Ekiqi = 0 i=l,N (4.50) where matrix $ is nonsingular, and (as it was shown in section 4.6.1) is equal to the matrix of undamped modes of the system MX + K0X = 0 These N scalar integro-differential equations (4.50) can now be solved by the Laplace transform method, or by the substitution method. The initial conditions in terms of variable Q = [<?i(0)...<7n(0)]r are following Q(0) = $-^(0) Q(0) = $-^(0) which can be used directly in the case of the Laplace transform method application (see (4.2) in section 4.1.1). In the case of the substitution method application the solution of each equation (4.50) is sought as 3=1 where Cj, pj are complex. The procedure of determination of Cj and pj is described in section 4.1.2. As a conclusion to this section, one can state, that the equations of motion for a discrete viscoelastic system (of a homogeneous material and with the proportional viscous damping) are always diagonalizable (i.e., can be decoupled) if Poisson's ratio is assumed constant and real. If Poisson's ratio is considered as a viscoelastic operator, or the system consists of several different materials (even with constant Poisson's ratios), then decoupling of equations, in general, is not possible. An additional investigation is required for this case. Chapter 4. Closed form solutions for discrete viscoelastic systems 91 4.7 Investigation of the conditions of overdamping of a SDOF system Assume that a M D O F viscoelastic system is diagonalizable, for example, consider a viscoelastic system of a homogeneous material with constant Posson's ratio and with the proportional viscous damping term (if present). Then it is expedient to consider a S D O F case [79] which will correspond to one of the modes associated with the system. The equation of motion is as follows: mx + cx + k{x - f T(t- T)x(r)dr) = 0 (4.51) Jo where the relaxation kernel will be assumed as for the standard linear solid model: T{t -T) = ae-"^ In equation (4.51) the quantity k is the instantaneous stiffness, e.g., for an axial element it is calculated as k — ^4, where ^^instantaneous Young's modulus, A—area of the cross section and L=length of the element. The eigenvalue problem which arises for this case is as follows / ' 6i b2 bs " b0 0 0 \ 0 p 1 0 0 + 0 -1 0 Q = 0 V 0 1 0 0 0 -1 ) 0 where bi = ca + k 62 = ma + c b3 = m b0 = fc(a — a) The corresponding characteristic equation is p3 + clP2 + c2p + c3 = 0 (4.52) where ci = a + 0 c2 = 0a + X2 c3 = A2(a - a) (4.53) Chapter 4. Closed form solutions for discrete viscoelastic systems 92 where A 2 = k/m, fi = c/m. Introduce some auxiliary quantities: Q = (3c2 - c2)/9 R = (9clC2 - 27c3 - 2c?)/54 D = Q3 + R2 (4.54) The roots pi,P2,P3 of (4.52) are determined according to known Cardan's formulas [80]. If D < 0 then all three roots are real and equal Pi = 2^Qcos{6lZ) - ci/3 P2 = 2^Qcos{\ir + 6lZ) - ci/3 (4.55) Pz = 2y[^Qcos{^K + 0/3) - CJ3 where 0 = arcosiR/sf-Q3) For the case D > 0 there are two complex conjugate roots and one real: Pi = -\(V + T)-Cl/3 + ^i(V-T) P2 = -\(V + T)- d/3 - ^ t ( V - T) (4.56) P l = V + T-Cl/Z where Note that the imaginary part of the complex roots is reduced to zero when D = 0. Therefore the cases D < 0, D = 0 can be combined to represent the overdamped case D < 0 (all three roots are real). Consider the conditions which are implied for parameters k, m, c, a, a of equation (4.51) to satisfy the condition D < 0. Substituting (4.53) in (4.54) one obtains Chapter 4. Closed form solutions for discrete viscoelastic systems 93 and Note that the quantities | , f, | are nondimensional. Thus the expression for D will be D = Q3 + R2 = X6(Q3 + R2) = A 6 D ( p ^ 2) (4.57) The sign of Z> is the same as of D, which depends on nondimensional parameters ^ , | , f. Denote them as a a 8 7 i = Also introduce notations 7 1 = A 7 2 = A 7 3 = A s = 7 i 7 3 + 1 q = 7 i + 73 T = 7 i - 72 The expression for /3 will be (intermediate derivations are omitted): f> ^ 3 - ^ 2 2 ^ , ^ 2 , ^ 3 D = —s s q sqr A—r H q r 27 108 y 6 y 4 27 y The boundary between the overdamped case (all roots are real) and underdamped case is expressed by the equation: D{s,q,r) = 0 (4.58) Consider a special case when the viscous damping is absent, i.e., | = 73 = 0, then s — 1, g = 7 ! = f and equation (4.58) becomes D = r2 + r(—q3 - -q) - —q2 + — = 0 K27q 3 y ; 27 y 27 From which one can determine *>=5[-F«'+W& ,-i« ) ,-4 (-F« , +£ ) 1 Recall that r = 71 - 72 = q - f , then <5>u=*4«,+i>^v^,-5«'+8'-1 ( 4 ' 5 9 ) Chapter 4. Closed form solutions for discrete viscoelastic systems 94 Thus the quantity | gets real values when S = ^q6 — |g 4 + q2 — 1 > 0. The case S = 0 corresponds to the minimum value of q = 71 = | . The value qmin is easily determined from the solution of cubic algebraic equation S(z) = ±fz*-1-z2 + z-1 = 0 which yields z i = z2 = z3 = 3 = q^in and 3Wn = v^ 3- This value is minimum one because S'(z) is which is semipositive, therefore S(q2) is a monotonically increasing function of argument q2. For the value gm;n = (f )min = \/3 one can calculate from (4.59) the corresponding value .a. , 1 2 1 \ ^ • ( ^ ) m m = 2 g m i n ( — g m i n + —) = g V O Thus for any 47 = | > i j m i n one can determine from (4.59) two values ( f )i,2 which yield D = 0 in (4.57). In this way the boundary of the region D < 0 can be determined. Two lines demonstrate this boundary (D = 0) in Fig. 4.17. The region between them corresponds to D < 0 (overdamped case). It should be noted that the condition a - < 1 (4.60) a must be taken into account, which implies that the stress in the material cannot relax to less than zero. To verify this, it is sufficient to substitute e(t) = e0 in the constitutive relation cr(t) = E(e(t) - f a exp[-a(t - T)]e(r)dT) Jo and to integrate to t = 00. If ^ > 1 then, e.g., for uniaxial extension, the stress will relax to a negative value which is not possible. Taking into account the constraint (4.60) results in the region shown in Fig. 4.18. Thus a viscoelastic system whose parameters j , | fall within Chapter 4. Closed form solutions for discrete viscoelastic systems 95 3.0 2.5 2.0 a IX 1.5 1.0 0.5 0.0 0.0 0.5 1.0 1.5 2.0 2.5 3.0 a IX Figure 4.17: The boundaries of the region where D < 0 (all roots are real). the region D < 0 will not be able to sustain any oscillatory motion. Further it may be noted that such overdamped behaviour is not possible when or in another form, when the relaxation time 1/a satisfies the following I T where T is the period of the undamped elastic system. For the case when the viscous damping term is present (c ^ 0 in (4.51)) the region of overdamping is broadened. For the values /3/A = 0.5, (3/\ = 1.0 these regions are presented in Fig. 4.19. For a series of prescribed points in the plane f, | (Fig. 4.20) the eigenvalues, namely, their ratios to quantity A = were computed and presented in Table 4.1 for the case when the viscous damping term is neglected (3=^ = 0. It follows from formulas (4.55), (4.54) that the fractions y (i = 1,2,3) will depend only Chapter 4. Closed form solutions for discrete viscoelastic systems 96 Figure 4.19: Overdamped case regions (D < 0) with different viscous damping: f = 0;0.5;l Chapter 4. Closed form solutions for discrete viscoelastic systems 97 a/X 1.0 1.5 2.0 , 2.5 3.0 OL IX Figure 4.20: Selected points in the plane of nondimensional parameters Point Coord-s ^, ^ Pi /A p 2 / A Ps/A 1 (1.5;1.0) -1.0 -0.25+i0.6614 -0.25-i0.6614 2 (1.5;1.5) 0.0 -0.75+i0.6614 -0.75-i0.6614 3 (v/3,1-2) -1.275 -0.2285+i0.6042 -0.2285-i0.6042 4 (V3,fV3) -V3 /3 -v^/3 -V3 /3 5 0.0 -V3/2+i0.5 -v /3/2-i0.5 6 (2.0;1.45) -1.588 -0.2058+i0.5513 -0.2058-i0.5513 7 (2.0;50/27) -1/3 -1/3 -4/3 8 (2.0;1.92) -0.0984 -0.6489 -1.253 9 (2.0;2.0) 0.0 -1.0 -1.0 10 (3.0;2.91134) -0.1835 -0.1835 -2.633 11 (3.0;2.96) -0.0463 -0.3288 -2.625 12 (3.0;3.0) 0.0 -0.382 -2.618 Table 4.1: Nondimensional eigenvalues for the selected points Chapter 4. Closed form solutions for discrete viscoelastic systems 98 Figure 4.21: Imaginary part of the complex root in terms of ^, | on nondimensional parameters j, j. The same holds for the case when complex roots are present (4.56). The quantities y will be called nondimensional eigenvalues. Table 4.1 can be expanded to present more points in the plane ~, | . Thus if it is necessary to prescribe certain eigenvalues pi, p2, P3 to the system (4.51), one can choose from Table 4.1 the quantities | , f to yield the most suitable nondimensional eigenvalues and then to rescale them by assigning the appropriate factor A. For the underdamped region (case D > 0), the imaginary part of one of the complex conjugate roots was calculated in terms of viscoelastic parameters | , | and it is illustrated in Fig. 4.21, where the lines of constant value are presented. The real part of the complex conjugate roots is analogously presented in Fig. 4.22. The third root (real root) is presented in Fig. 4.23. Chapter 4. Closed form solutions for discrete viscoelastic systems 99 Figure 4.23: Real root in terms of | C h a p t e r 5 D y n a m i c s of an eng ine -mount s y s t e m 5.1 D e r i v a t i o n of equat ions of m o t i o n The system which is under consideration consists of the following rigid bodies: 1) body A -the engine framework (it includes all the non-moving parts of the engine) with mass mA and tensor of moments of inertia 3 A taken at the mass centre (point A, Fig. 5.1); 2) body B - the crankshaft with the mass TUB and tensor of moments of inertia 3B taken at the mass centre (point B); 3) the pistons, each with the mass ra^, i = l,iV. Note that connecting rods are not included, because their mass is assumed to be distributed between the crankshaft and the pistons. It is a customary assumption [61] to distribute the mass of the connecting rod m„ into two concentrated masses mi, m 2 using the following relations miLi = m2L2 m\ + m 2 = mcr where L\ = distance from the axis of the crank joint to the e.g. of the connecting rod, and L2 = distance from the axis of the piston joint to the e.g. of the connecting rod. One can introduce the following generalized coordinates for the system: displacements uiAi^AtUsA of the body A mass centre in a ground based system of coordinates e i , e 2 , e 3 (Fig. 5.1) and angular rotations <t>i,tj>2,<f>3 about axes e i , e 2 , e 3 respectively. The angular ro-tations are presumed to be very small, so that they can play the role of generalized coordinates [81]. The angular speed of the crankshaft in rotation about axis b—b is presumed to be constant (u>) with respect to the ground based system of coordinates e i ,e 2 ,e 3 . 100 Chapter 5. Dynamics of an engine-mount system 101 Figure 5.1: Schematic view of engine framework and crankshaft One can begin with a derivation [82] of the kinetic energy of the system. For the body A it will be: TA = ^rnAvA + • 3A • u)A (5.1) where the second term is the kinetic energy of rotational motion about the mass centre (point A), and is calculated as a double scalar product of the vector of angular velocity and tensor of moments of inertia. The vector of angular velocity is expressed in terms of the generalized coordinates as follows UA — </>iei + 4)2^2 + <f>3e3 Tensor 3A will be also expressed in terms of projections onto the axes parallel to ei,e2,e3, and therefore will be time-dependent. However assuming that the angles 4>x,<f>2,(l>3 are very small, it will be assumed that 3A is constant. Expression (5.1) can be rewritten as follows TA = \™>AvA + ~[4>i,<j>2,JA [4>1,4>2,4>3]T Chapter 5. Dynamics of an engine-mount system 102 where JA is the matrix of components of tensor 3A. The kinetic energy of body B (the crankshaft) will be: TB = ^mBvB + \JJB -JB-UB (5.2) where the second term is the kinetic energy of rotational motion about the mass centre of body B. The vector of angular velocity is expressed in terms of generalized coordinates as follows U>B = wei + <f>2e2 + <f>3e3 The tensor JB is time-dependent when expressed in terms of ground basis vectors ei,e 2,e 3. Consider a set of basis vectors gi,g2,g3 rigidly connected with the body B, in which tensor of moments of inertia 3 B has constant components 3B = JS'&gj (5.3) where summation on repeated indices is assumed. Also introduce an auxiliary rotating set of basis vectors i i , i2» 13. which has the angular velocity vector wei. These basis vectors , i 2 ,13 can be expressed in terms of the ground basis vectors as follows ii 1 0 0 ei 12 = 0 C 0 S 7 siny e 2 h 0 —sinj cosj e 3 where 7 = u>t + <j>0 (angle of crankshaft rotation). Given the angular displacements <f>i,<j>2,<f>s at a given moment of time t, the basis vectors gi>g2,g3 rigidly connected with the body B can be expressed in terms of the ground basis vectors as follows gl 1 0 0 0 -<t>2 ei g2 = { 0 cosy siny + fasiny — (f>zCosy —cpxsiny (j)\COSy } e 2 . S3 . 0 —sinj cosy (j)zsiny + (f)2cosy —<j>iCos*y —(pisin-y e3 Chapter 5. Dynamics of an engine-mount system 103 or introducing notation T for the first matrix term and $ for the second, one obtains g i = (Tik + $ik)ek = Rikek i, k = 1,2,3 (5.4) where summation on repeated indices is assumed. Matrix R will be called the transformation matrix. Note that components of matrix T in (5.4) are functions of y only, whereas components $ i k are functions of angular coordinates cf>i,cf>2,<t>3 as well as of 7. Thus the basis gi,g2,g3 is related with the basis e i ,e 2 , e 3 through matrix R: 1 (j>Z ~<i>2 R — fasiny — cpzcos'-j cos-y — cpisinj siwy + cpicosy cpzsiwy + (f>2cosy —siwy — (j>\cos^ cos~f — cpisiwy The details of the derivation of the matrix R are presented in Appendix C. The contribution of components containing 4>i W 'H he neglected, because the angular displacements are assumed very small (|</>;| << 1, i = 1,2,3). Note that comparing the C-norms [72] of these functions, one obtains ||(/>iC067||c = ||^>i5m7||c = |0iI << ||co37||c = ||sin7|| c = 1 7 C [0,27r] Thus the final form of the transformation matrix R will be further assumed as 1 <f>3 -<f>2 R = (j>2siwy — cpzcosy cosj siwy 4>zsiwy + <f)2cosy —siwy cosj Substituting (5.4) into (5.3) one obtains Jfi = Jij''Rik^kRjm^m = RikJfj"Rjm^k^m (5-5) where summation on repeated indices is assumed. Using matrix form, expression (5.5) will correspond to JB = RTJBoR (5.6) Chapter 5. Dynamics of an engine-mount system 104 Note that matrix JB is a function of time, and strictly speaking, of angular coordinates <j>i,4>2,(f>3, because matrix R is a function of these variables. Now one can rewrite the expression for the kinetic energy of the body B in matrix form. Namely, TB = ^rnBv2B + JB{t) [w, <£2, <£ 3] r (5.7) Figure 5.2: Determination of velocity vector \B It is left to express the velocity v>3 in terms of the generalized coordinates. The following vector relation hold (Fig. 5.2): rs = *A + rAD + rDB (5.8) where *A = UlA^l + U2A*2 + U3AG3 Vectors YAD< *DB are represented in the coordinate system P i , P 2 5 P 3 (which is rigidly con-nected with body A, Fig. 5.1) as follows *AD = -SlPl + S 2 p 2 + 5 3 P 3 *DB = pCOS-yp2 + psiwypz Chapter 5. Dynamics of an engine-mount system 105 where point D is the projection of mass centre B on the axis of shaft rotation, and quantities si,s2,s3 are constant and assumed known. Differentiating (5.8), one obtains V B = v,4 + rAD + rDB (5.9) Differentiation of TAD yields TAD = uA x *AD = (^2*3 - ^s^ex + (cp3si - <j>is3)e2 + {-<t>2si + <£is 2)e 3 where vector TAD in the multiplication above was represented by the components Si,s2,s3 in the basis ei,e 2,e 3 instead of basis pi, p 2 , P 3 (assumption about smallness of angles <j>\,(j>2,<l>3 was used). Differentiation of TDB yields TDB = ">A X rDB — pusin-yp2 + ptocosjp3 Neglecting the first term and replacing basis vectors p 2,p3 by e 2,e 3 (assumption about smallness of angles </>i,</>2,</>3 is used again) one obtains TDB = —pusin'je2 + pu>cos~/e3 Thus from (5.9) the velocity of mass centre of body B will be vs = {u1A+if>2S3-<i>3S2)e1 + (u2A+<i>3S1-<p1s3-pu;siny)e2 + and the kinetic energy of linear motion of body B will be as follows TBI = ^rnB([u1A + <j)2s3 - cp3s2}2 + [u2A + (j>3Si ~(j>is3 - pusinf]2+ + [U3A - 4>2Sl + <j)\S2 + pUJCOSj]2) Consider now the determination of the kinetic energy of the pistons (Fig. 5.3). It will be presumed that pistons are moving in the direction of p 3 (Fig. 5.1), i.e., the engine cylinders Chapter 5. Dynamics of an engine-mount system 106 Figure 5.3: Determination of kinetic energy of the pistons are parallel to each other. The derivation will be shown for one cylinder and then generalized for the case of several cylinders. The velocity vector of point G (mass centre of the piston) can be calculated as V G = V J 4 + W A X rAG + vGrp3 (5.10) where vor is the velocity of linear motion of the piston with respect to the cylinder along the vertical direction. One can find that this velocity VQT will satisfy the following relation: VQT — Iwsin^ptanS + luicosyp (5-H) where 7 P = 7 + 7 p 0 , 7Po is the initial phase of the piston, and sinS — —cosfp h The radius-vector TAG can be represented as follows *AG = *AD' + RD'G = (si + / ) p i + S2P2 + (-83 + lsinjp + l2cosS)p3 (5.12) Chapter 5. Dynamics of an engine-mount system 107 where / is the offset of the cylinder line (point D') with respect to point D. Introduce the notation: si + f = a\ which will be used below. Upon the calculation of u>A X rAG in (5.10), vector rAG is represented by the same components as in (5.12), but in the basis e i , e 2 , e 3 instead of basis P i , P 2 , P 3 (assumption about smallness of angles <f>i,<f>2,<f>3 is used) and analogously the term vGrp3 is replaced by *>Gre 3. Substituting relations (5.12), (5.11) in (5.10) and taking into account the mentioned above replacements, one obtains Now consider the case of N cylinders, introduce the following auxiliary notations which will be used later: N Mp = 53 mpi hi = s3 + Isin^pi + l2cosSi di = Icos^piU) — l2sinSiSi v G = (xA + 4>2{s3 + Isinip + l2cos8) — s2<f>3)ei + + 0/A - <}>i{s3 + Isinjp + l2cosS) + ai</>3)e2+ +(zA + luisin^ptanS + lujcos'jp + s2(f>\ — ai<^2)e3 N A = ^rripihi N P = Y,m.pihl (5.13) where i=l i=l i=l mass of ith cylinder. The kinetic energy of all pistons will be: 1 N . . . . Tp = - 53 m-piK^A + fahi - s2<j>3)2 + (yA - ipihi + au4>3)2+ i=l +(zA + s2<f>i — au(f>2 + Iwsin-ypitanSi + lucosjpi)2] (5.14) where hi were defined in (5.13). Chapter 5. Dynamics of an engine-mount system 108 The total kinetic energy of the system will be T = TA + TB + TP where TA, TB, Tp are given by (5.1), (5.2), (5.14) respectively. Introduce a new notation for the generalized coordinates, namely, qi = u1A q2 = u2A q3 - u3A qi = 4>\ 95 — 02 ?6 = 03 Using Lagrange's equations [81]: one obtains the following equation in matrix form: M(t)q + D(t)q = Q(t) (5.16) where the inertia matrix M(t) is symmetric, time-dependent and has components Mij defined by: M n = mA + mB + Mp M12 = 0 M 1 3 = 0 M i 4 = 0 M i 5 = mBs3 + A Mie = - m B s 2 - Mps2 M22 — mA + mB + Mp M23 = 0 M 2 4 = - m B s 3 - A M25 = 0 M 2 6 = m B si + 5 M 3 3 = + m B + M p M 3 4 = mB*2 + Mps2 M35 = -mB&\ - S M 3 6 = 0 (5.17) M 4 4 = Jn + mB(sl + s22) + P + MpS22 M 4 5 = - mBs1s2 - Ss2 M 4 6 = Jf3 - mBs1s3 - W Chapter 5. Dynamics of an engine-mount system 109 M 5 5 = J22 + mB{s\ + s\) + Jg + P + U M56 = J23 ~ mBs2s3 + Jg - As2 M 6 6 = J* + mB{s\ + s22) + jg + Mps\ + U where JA = components of the tensor J ^ , and J^(t) = components of the tensor 3B(t). The expressions for Jg,Jg<Jg are as follows 1 jg(t) = jg°sin2~f + jg°cos2j + Jg°sin2'y The quantities Mp, A, P, S, W, U used in (5.17) were defined in (5.13). The matrix coefficient D(t) associated with the velocity vector will be: D(t) = 0 0 0 0 B 0 0 -B 0 0 0 0 0 0 0 -B 0 2 £ £ 1 T J V M i 0 0 0 0 0 - E i l i TUpidudi B 0 0 0 0 0 0 Jg + 2 E i l i rripihidi Jg - S 2 B E i l i rripiaudi jg ~ s2B J33 + G(t) (5.18) where the additional term G(t) arises from (j = 5,6) and will be presented below. The quantity B used in (5.18) was defined in (5.13). Below the assumptions which were made in the process of the derivation of term G(t) are presented. These assumptions allow the obtained matrices M, D to be independent of generalized coordinates. Namely, upon the calculation jt{^), {j = 5,6), the assumption that the angular dis-placements (j>\,(j)2,4>3 are small is used, so the second term (components $;&) in (5.4) are neglected in comparison with components T^, i.e., transformation matrix R is assumed to Chapter 5. Dynamics of an engine-mount system 110 be a function of 7 only. However upon the calculation (j = 5,6) this term (matrix $) is taken into account, i.e., using (5.7) -Q-; = 2 ^ ft * » J - ^ l w ft 96] J = 5,6 where (using (5.6)) ^ = f r t + H V f ,' = 5,6 where for matrix R (when it is not under differentiation) is assumed dependence only on 7. Therefore the calculation of the term yields: B dq5 2 ' where the terms containing 95175, q5q6, q6q6 have been neglected in comparison with the terms containing 950;, q6uj, u2. Analogously dTi - = liG^w2 + 2U,G<8K + 2wG<$qe) dq6 2 where G{$ = l « n 2 7 ( J*° - J*°) + J3B2°cos21 G{$ = \sin2j(jg° - J2B2°) - J 3 B 2 ° ^ 2 7 = - Jn + Jf2°8in2i + Jgfcos2'! + J^air^ Therefore the additional matrix G(t), which can be called a gyroscopic matrix, will be Chapter 5. Dynamics of an engine-mount system 111 Figure 5.4: Assumed location of the mounts expressed as follows G(t) = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 One can see that with an axisymmetric body rotating about the axis of symmetry, the terms — Gvz — 0 a r | d G$ = —G^ which yields a skew-symmetric addition to the velocity matrix D, but in the general case when the body is not axisymmetric this addition to the velocity matrix will not be skew-symmetric. Now consider the generalized forces Qj in (5.15), or vector Q(t) = [ Q i , Q 6 ] T in (5.16). One contribution to these generalized forces are the forces arising from interaction of body A with the mounts (Fig. 5.1). The location of the mounts is assumed as shown in Fig. 5.4, where point A is the mass centre of the engine framework, and k3,k2,ki are viscoelastic springs parallel to axes e i , e 2 , e 3 respectively. Chapter 5. Dynamics of an engine-mount system 112 Each of the mounts is modelled as a combination of 3 viscoelastic uniaxial elements with the following relation between the force F and elongation S: F(t) = k[8(t) - /"*r(< - T)S{r)dr] Jo where parameter k = E 7 , E = instantaneous Young's modulus, and 7=the stiffness factor depending on the mount's geometry and Poisson's ratio. The stiffness factor in each direction (for each spring) can be, for example, determined by application of the finite element method to a real mount 3-D model. The derivation of the stifness matrix operator will be conducted as if springs were purely elastic, and only in the final form the stiffness matrix will be converted to a matrix-operator by replacing Young's modulus E by the operator E. The relation between the generalized coordinates and generalized forces arising from forces and moments (denote them as vector QM) transmitted through the mounts to the engine framework can be written as follows QM = -Kq (5.19) where matrix K will be called the stiffness matrix. To derive this matrix the definition of generalized forces is used. According to this definition for a system with N points at which the external forces are applied, the generalized force will be [81]: Qi = ,tir- F> ( 5- 2°) *=i d $ where r*. is the radius-vector of the point k, and Fk is the external force vector acting at the point k. Given generalized coordinates q = [qu ...,qe], the following resulting forces and moments (in projections on axes e i , e 2 , e 3 ) are transmitted from the mounts to the engine framework: Fx = -4*391 + 4&3#g5 Chapter 5. Dynamics of an engine-mount system 113 Fy = -4k2q2 - 4k2HqA Fz = -4kxq3 MxA = -4k2Hq2 - 4(k2H2 + hA2)qA MyA = -4{kxA\ - k3H2)q5 + 4k3Hqi MzA = {-4k2A\ - 4k3Al)q6 where index A means that the moments are calculated with respect to the axes passing through the point A. Now one can calculate the generalized forces contributed from the mounts using (5.20), namely, for this case the components of the vector QM will be: = FX Q? = Fy Qf = Fz Qf = MxA Qf = MyA Q f = MzA From these expressions the stiffness matrix K mentioned in (5.19) can be obtained. It will be symmetric and with the following components: Kn = 4k3 K12 = K13 = K1A = K16 = 0 K15 = -4k3H K22 = 4k2 K23 = K25 = K26 = 0 K2i = 4k2H K33 = 4fcx K3i = K35 = K36 = 0 (5.21) K44 = 4(k2H2 + kUl) K45 = Ki6 = 0 K55 = 4{hAl + k3H2) K56 = 0 K66 = \k2A\ + 4k3A\ Note that the generalized force contribution from gravity does not depend upon the gener-alized coordinates. The same situation exists regarding to the generalized force contributions from the gas pressure in the cylinders, except that they will be functions of time. Derivation Chapter 5. Dynamics of an engine-mount system 114 of them is straightforward and is omitted here. Also note that modeling of mounts in the form of springs is not a unique approach which can be appropriate, in general the components of matrix K can be derived according to any finite element model of the mounts and their orientation, or can also be obtained analytically using other elements, e.g., beam elements instead of bar elements. Now having collected the terms resulting from (5.15) which depend on generalized coor-dinates, velocities, accelerations on the left-hand side of the equations and having transfered all the other terms to the right-hand side, the equation of motion of the engine-mount system in matrix form can be written as follows M{t)q +D{t)q +Kq = F{t) (5.22) where matrices M and D were given respectively by (5.17), (5.18). The matrix-operator K follows from (5.21) by replacing ki k\ i — 1,2,3. Note that matrices M, D are periodic functions of time: M(t) = M(t + T), D(t) = D(t + T),T = 2TT/UJ. The right-hand side of (5.22) will be as follows F(t) = m) F2(t) ... F6(t)]T where Ncyl m = £ f t ) Ncyl F2{t) = mBpu;2cosut + £ Ff°s(t) i=i Ncyl Ncyl Fa(t) = mBpu?sinut - (mA + mB + Mp)g - £ m^cw + £ Ff?s(t) (5.23) i=l j=l Ncyl Ncyl F^t) — —mBps3u)2coswt + mBps2w2sinwt — ^ mpiS2VGri + £ MxA[Ffa°(t)] i=l j=l Ncyl Ncyl FB(t) = -G{$u2 - JxB2(t)u - mBpsxu2sinujt + £ m^vcri + £ MyA[Ffa3(t)] ^ i=i j=i Chapter 5. Dynamics of an engine-mount system 115 I Ncyl Mt) = ^iV - J*(t)v + rnBpSlu,2cosu,t + £ MzA[F?aa(t)) 1 i=i where vori is given by (5.11), and = Ji2°sin-f + Ji3°cos~f = -J^cos^ + J^siwy Ff°', Fj"*, Fjzs = X,Y,Z projections of the resulting gas force in jth cylinder, and MxA[Ffas], MyA[Ff], MzA[Ff] = resulting moments of the gas forces in j th cylinder with respect to the axes X,Y,Z passing through the point A. The expressions for J ^ ( £ ) > ^ ( * ) are obtained by differentiation of ^(O = J i 2 c o s l ~ Jg°sin~i jf3(t) = J^siwy + J^3°cosj Remark 1. Note that forces and moments related with the gas pressure should be only taken into acount when there is a gradient of pressure in the combustion chamber. If the pressure is uniformly distributed inside the chamber at each instant of time, then the resulting force is zero at each instant of time, and all terms with superscript "gas" in the right-hand side of equations (5.23) vanish. Remark 2. Note that for some engines (slow, or medium speed marine engines) there is a propeller drive shaft which joins the crankshaft and the propeller. In the dynamic model considered here the interaction between the crankshaft and the propeller drive shaft is reduced to a torque moment (Mx) upon the assumption that a flexible coupling will minimize all other forces and moments of interaction. Therefore its contributions to the generalized forces are neglected. The moment Mx does not contribute into the generalized forces (5.20), because the rotation angle of the crankhaft is prescribed 7 = uit + </>0, i.e., 7 is not considered as a generalized coordinate. Chapter 5. Dynamics of an engine-mount system 116 5.2 Parametric resonance investigation Consider a general equation M(t)q + D(t)q + Kq = 0 (5.24) where matrices M, D of order n, T-periodic, and q is n x 1 vector. In general K can be T-periodic as well. To investigate parametric resonance phenomena (dynamic instability) [83], [84], [85] it is sufficient to consider the conditions that give rise to unbounded solutions q(t) of (5.24). It is assumed in this section that no viscoelastic properties are present in the system which means K — K (elastic case). This represents the case when parametric resonance arises more easily. A general theory of linear differential equations with periodic coefficients [83] allows a treatment of equations like (5.24). However the case when viscoelastic (hereditary) elements are present brings a new type of equation, namely, integro-differential where some coefficients are time periodic. The analysis of resonance conditions for this case can be a subject for future work. The equation (5.24) can be represented in equivalent state-space form A(t)X + B(t)X = 0 where A(t) = D(t) M{t) B(t) = K 0 X = q M(t) 0 0 -M(t) . q . or in abbreviated form X = P(t)X (5.25) where P = -A~XB, 2n x 2n matrix, and P(t + T) = P(t). Note also that matrix A is nonsingular, upon the presumption that the mass matrix M is nonsigular for all t C [0,T]. Chapter 5. Dynamics of an engine-mount system 117 Consider now X as a 2n x 2n matrix as well. Assuming the initial conditions X(0) = I, the matrix-function X(t) will be called a matrix of fundamental solutions (or matrizant), and the matrix X(T) is called the monodromy matrix. According to the Floquet-Lyapunov theorem [83] the matrix of fundamental solutions can be expressed as X(t) = S{t)etR where 5(0) = I, R = £ZnX(T), S(t + T) = S{t). Eigenvalues of R, rn are called charac-teristics exponents, eigenvalues of X(T), & are called multipliers of the system. Note the relation between them & = er>iT t = l,2n Therefore X(T) the monodromy matrix (namely its eigenvalues) yields all information required to analyse the stability of trivial solutions of (5.25), and consequently of (5.24), namely [83]: i) all solutions bounded on [0,oo], if multipliers are inside, or on the unit circle (the latest case with simple elementary divisors), ii) asymptotic stability, if multipliers are inside of unit circle, iii) instability of solution, if at least one multiplier is either outside the unit circle, or on the unit circle with multiple elementary divisor. There are several methods to obtain regions of stability and instability for systems de-scribed by differential equations with periodic coefficients. For example, for the Mathieu equation such methods as straightforward expansion (in power series of parameter), the method of strained parameters, Whittakers's method, and the method of multiple scales are described in the literature [86], [83], [87]. In this study it is proposed to conduct the stability analysis by obtaining the eigenvalues of the monodromy matrix X(T). Note that in the case of a general system of coupled Chapter 5. Dynamics of an engine-mount system 118 B A Figure 5.5: Example of a SDOF system differential equations with periodic coefficients determination of the matrix of fundamental solutions X(t) is not a trivial task. Here it is proposed to conduct a numerical integration on the interval [0,T], for equation (5.25) to obtain X(T) starting from X(0) = I. An explicit 4-stage Runge-Kutta method will be used as a numerical integrator (a computer program was written for this purpose). At first, for illustration purposes, a special case of one-degree-of-freedom system was considered. An example of mechanical system is shown in Fig. 5.5, where there is a crankshaft (body B) rotating about axis b — b and the framework (body A) is allowed to rotate about axis a — a which passes through the mass centre of the framework (point A). The angle of rotation q5 is assumed to be very small. The framework (body A) is attached to the ground by elastic springs which create the rotational stiffness denoted by K55. The equation of motion for this system is obtained from the general equation (5.22) by constraining all degrees of freedom except q5 and assuming that Mp — 0 (no pistons). Assume also that the mass center of body A coincides with point D (see Fig. 5.2). Recall that point D is the projection of the mass centre of body B on the axis of crankshaft rotation. Then the motion of this system (with forcing functions set to zero) can be described by Chapter 5. Dynamics of an engine-mount system 119 the following equation: 3 (1 — esin2tot)q5 — -eu>sin2u>tq5 + \2qs = 0 (5.26) with parameters jBo _ jBo U22 US3 e A = y/KnKJ* + J2B2°) (5.27) «^22 "I" ^ 2 2 ° where J ^ " , JB3° are components of the tensor of moments of inertia defined with respect to axes g i , g 2 , g 3 rigidly connected to the body B. From (5.27) it follows that e < 1, because JJ2 > 0> J£? > ° . J22 > 0-The rigidly connected system g i , g 2 , g 3 is chosen in the following way. The axis gi is taken parallel to ei. As far as, g2,g3 are concerned, they are taken in such a way that the component Jf30 = J32° = 0. Without loss of generality, it also assumed that at the instant t = 0 bases g i , g 2 , g 3 and e i ,e 2 , e 3 coincide. Note that the unbalance radius p which contibutes to the external forcing function does not appear in equation (5.26), because it is not involved in consideration of the homogeneous solutions. One can rewrite the equation (5.26) using a new variable: 7 = u>t, thus t = j/u>, i.e. t can be considered as function of 7, t = t(~f). Then q5(t) — 95(7/^) = £5(7) and 9*6(7) = ft;1; q\(j) = ftA Thus 95 = wg'B(7) 95 = u2q"5{~f) Substituting this into (5.26), one obtains: 3 (1 - esin27)u2q"5(i) - -£U)2sin2<yq'5{~/) + \2q5(j) = 0 Denoting S = X2/UJ2 one can write (1 - e « n 2 7 ) < p ' B ( 7 ) - j |e«n2 7 t f 5 ( 7 ) + £95(7) = 0 (5.28) Chapter 5. Dynamics of an engine-mount system 120 0.10 0.05 S 0.00 -0.05 -0.10 0.85 0.90 0.95 1.00 1.05 1.10 1.15 8 Figure 5.6: Regions of stability and instability for the SDOF system (Fig. 5.5) Note that the coefficients of this equation have a period ir. This equation was numerically integrated and the monodromy matrix .X^Tr) ( 2 x 2 in state-space form) was obtained. Then its eigenvalues (multipliers) were evaluated. Regions of stability and instability (dark area) are presented in the plane of the nondimensional parameters 5, e in Fig. 5.6. For the unstable (dark) area the absolute value of one of the multipliers was > 1. One can see that the greater parameter \s\ which reflects nonsymmetry of the body B, the wider frequency region where resonance can occur. Note that the line S = 1 corresponds to the frequency rotation ui = A, where A is the natural frequency of the system. 5.2.1 Parametric resonance analysis of an engine-mount system The parametric resonance analysis of a 6 degree-of-freedom engine-mount system will be conducted on the basis of equations (5.22), where the right-hand side is assumed zero, because homogeneous solutions are of interest in this section. The equations (5.22) will be represented in state-space form (5.25). Thus there will be a system of 12 coupled differential equations with periodic coefficients. Note again that in equation (5.22) it will be assumed K = K Chapter 5. Dynamics of an engine-mount system 121 No. Real part Imaginary part Absolute value 1 2 3 4 5 6 7 8 9 10 11 12 0.1458943369E+00 0.1458943369E+00 0.4635743596E+00 0.4635743596E+00 -0.3079144813E+00 -0.3079144813E+00 -0.6181632403E+00 -0.6181632403E+00 -0.9896163146E+00 -0.9896163146E+00 -0.8928809929E+00 -0.8928809929E+00 0.9893001778E+00 -0.9893001778E+00 0.8860580191E+00 -0.8860580191E+00 0.9514140383E+00 -0.9514140383E+00 0.7860497493E+00 -0.7860497493E+00 0.1437343028E+00 -0.1437343028E+00 0.4502927186E+00 -0.4502927186E+00 0.1000000000E+01 0.1000000000E+01 0.1000000000E+01 0.1000000000E+01 0.1000000000E+01 0.1000000000E+01 0.1000000000E+01 0.1000000000E+01 0.1000000000E+01 0.1000000000E+01 0.1000000000E+01 0.1000000000E+01 Table 5.1: Eigenvalues of monodromy matrix (multipliers). Example 1, stable case (elastic mount case.) The input data for an example of engine-mount system (a medium-speed diesel type multicylinder engine was considered [62]) are presented in Table D. l (column 1) of Appendix D. It was assumed that all elastic springs ki,k2,k3 (Fig. 5.4) have the same stiffness 0.206e + 07 N/m. A numerical integration on the interval [0,T], for the equation of motion represented in the state-space form (5.25) was produced and X(T) (12 x 12 monodromy matrix) was obtained. An explicit 4-stage Runge-Kutta method was used as a numerical integrator. The eigenvalues i = 1,12, of the monodromy matrix are presented below for two examples. Input data for the 1st example are in Table D. l (column 1, Appendix D). In this example (the crankshaft with Jg° = Jg°) all eigenvalues (Table 5.1) are on the unit circle and not repeated which means a stable case, i.e. the homogeneous solutions of (5.22) are bounded. For the second example just one parameter was changed in input data (column 1, Table Chapter 5. Dynamics of an engine-mount system 122 No. Real part Imaginary part Absolute value 1 2 3 4 5 6 7 8 9 10 11 12 0.1990503196E+00 0.1990503196E+00 0.4635742820E+00 0.4635742820E+00 -0.2761777116E+00 -0.2761777116E+00 -0.6181632403E+00 -0.6181632403E+00 -0.9496800600E+00 -0.9496800600E+00 -0.8930918482E+00 -0.8930918482E+00 0.9799898711E+00 -0.9799898711E+00 0.8860580620E+00 -0.8860580620E+00 0.9611066333E+00 -0.9611066333E+00 0.7860497493E+00 -0.7860497493E+00 0.3132183178E+00 -0.3132183178E+00 0.4498769529E+00 -0.4498769529E+00 0.1000000589E+01 0.1000000589E+01 0.1000000000E+01 0.1000000000E+01 0.1000000044E+01 0.1000000044E+01 0.1000000000E+01 0.1000000000E+01 0.9999989655E+00 0.9999989655E+00 0.1000001161E+01 0.1000001161E+01 Table 5.2: Eigenvalues of monodromy matrix (multipliers). Example 2, unstable case D.l), namely, a crankshaft with Jg" — l.SJg" was considered which gives a greater contri-bution to the time-dependent components of matrices M and D (see (5.17), (5.18)). For this example one can see that absolute values (Table 5.2) of some eigenvalues are slightly greater than 1 (outside the unit circle) which means a weak instability eventually leading to unboundness of the homogeneous solutions. In other words parametric resonance arises. The numerical integrations were produced with different step sizes to verify the obtained values of the monodromy matrix X(T), and consequently its eigenvalues (multipliers). The results (Tables 5.1, 5.2) are presented for the engine's operating frequency 10 Hz (T = 0.1 Several other engine-mount systems operating on different frequencies were numerically tested in this study. Some of them, on certain frequencies, produced such absolute values of the multipliers as 1.2 - 1.4 which are far from the unit circle. As a conclusion one can say that when the relative contribution of time-dependent compo-nents in matrices M, D exceeds a certain level then instability can occur at certain frequencies. Chapter 5. Dynamics of an engine-mount system 123 A full analysis of the input parameters of engine-mount system in terms of their influence (ob-viously different for different parameters) on the multipliers of the system (instability criteria) can be viewed as a subject of future work. 5.3 Steady-state response for the case of periodic matrices M and D Consider a general equation of motion where matrices M, D are T-periodic, and forcing vector-function F(t) is T-periodic as well. A procedure to obtain the T-periodic solution (steady-state response) is proposed using complex Fourier series, and avoiding numerical integration in time starting from some initial conditions. Real periodic matrices M(t), D(t) are represented by complex Fourier series (using con-jugate terms) in the following way: M(t)X(t) + D(t)X(t) + KX{t) = F(t) (5.29) OO oo M(t) = £ ( M , + Mke ) D{t) = J2(Dkeikut + Dke-ikut) where Mk and Dk are complex matrices and u> = 2ir/T. The forcing function is extended into complex form: oo k=—oo where Fk are complex vector-columns. The steady-state response solution is sought in analogous complex form: oo X(t) = £ X, (5.30) k=—oo where Xk are complex vector-columns. Chapter 5. Dynamics of an engine-mount system 124 Substituting (5.30) into (5.29) gives oo Xke fcV 2(M,eaw' + Mte-ilut) + ikwY,{Dieilut + Ae"""*) + E.{ku>)K0 1=0 1=0 = F*eikut ( 5 - 3 1 ) k=—oo In this illustration only complex Young's modulus E^kuj) for the complex stiffness matrix is introduced, although in the general case of isotropic material one can introduce a two-part form (3.3). Note that in the case of numerical integration in time, instead of such parameters as Young's complex moduli given on a set of certain frequencies, the relaxation kernel parameters of the matrix-operator K would be required. The equality in (5.31) is achieved by equating the coefficients of like harmonics e%kwt k = 0, +1, +2, +3, +4,... in the left-hand and right-hand sides of the equation. Therefore a linear system of infinite number of complex equations with respect to the unknown complex vectors Xk, k = 0, +1, +2, +3, +4,... is obtained. In practical applications the number of terms taken in (5.30) will be finite, and conse-quently the system of equations will be of a finite order (contribution of harmonics expwt is assumed to decrease with p -» oo). Note that Re[X(t)] (real part) and Im[X(t)] (imaginary part) of a found complex solution X(t) will be solutions of (5.29) for Re[F(t)] and Im[F(t)} respectively. To illustrate the developed method some numerical results are presented below. Equation (5.22) obtained for an engine-mount system can be considered as a particular example of (5.29). Introduce a new notation of the generalized coordinates for this example xi = qit i = 1,6. In vector form X(t) = q(t) = [x1(t) x2(t) x3(t) x4(t) xs(t) x6(t)]T In the previous section the multipliers (complex eigenvalues of the monodromy matrix) were obtained (Table 5.1,5.2) for some examples of engine-mount system. For these examples Chapter 5. Dynamics of an engine-mount system 125 none of the mulitipliers was equal 1. According to the results from [83] this means that the inhomogeneous equation (5.29) with K — K will have a unique T-periodic solution, which will be referred here as the steady-state response. Note that this result is proved in [83] for the case of differential equations with time-periodic coefficients, which is our case when the stiffness matrix is assumed constant (or time-periodic). In this section a viscoelastic matrix-operator K replaces K in (5.29), however it will be still assumed that equation (5.29) has a unique T-periodic solution (proof is not considered here). The real forcing function (see 5.23) is extended into complex form as follows F(t) = F0 + Fxeiut + F2e2i"\ F^ = F _ 2 = 0 provided that Re[F(i)\ is the given real forcing function (5.23). The quantity u> is the crankshaft rotation frequency [rad/s]. For this example the force contributions arising from non-uniform gas pressure distribution were neglected, i.e., the terms with superscript "gas" in (5.23) were dropped. The time-dependent mass matrix (5.17) and velocity matrix (5.18) were represented in the following form 2 2 M(t) = Y(Mkeikwt + Mke-ik^) D{t) = ^ / V " " + Dke~ik,Jt) k=0 k=0 and a Fourier representation of the complex solution is chosen as follows 4 X(t) = £ Xkeikwt fc=-4 Note that for this example it was enough to keep 9 terms (k =-4,-3,..,3, 4) in (5.30). The numerical results obtained with a larger number of terms were practically identical to ones obtained with 9 terms. Note that Re[X(t)} will be the solution of (5.29) for Re[F(t)}, i.e., for the given real forcing function. Chapter 5. Dynamics of an engine-mount system 126 3.0E-4 •o CO 2.0E-4 \-C -2.0E-4 -< .3 .0E -4 C—'— 1 — 1 — 1 — 1 — 1 —'— 1 — 1 — 1 — 1 —'— 1 — 1 — 1 — 1 — 1 — 1 — 1 — 0.00 0.05 0.10 0.15 0.20 Time, s Figure 5.7: Steady-state responses. Crankshaft with J%2° = Jg", elastic mounts Two models were compared: in the 1st the mass matrix M, and the velocity matrix D were assumed constant, i.e., M = M0, D = D0, where M 0 , D0 are the coefficients in Fourier series. In the 2nd model M, D were time-dependent. The input parameters of the engine-mount system are presented in Table D. l (column 2) of Appendix D. The stiffnesses (complex in general) of the springs (Fig. 5.4) were assumed equal = fc»2 = K3. The steady-state responses are demonstrated by the angular displacements x4,x5,x6 of the engine framework for the duration of one period. There was virtually no difference between the 1st and 2nd model (Fig. 5.7). The solid lines represent both solutions. However when a crankshaft with Jg" ^ Jg" was considered (all the rest input parameters in Table D. l were the same), namely, Jg0 = 9453 kg * TO2, Jg° = 8453 kg * TO2, then time-dependent components in matrices M, D become more noticeable and cause a larger difference in the responses (Fig. 5.8). Note that the mounts were assumed elastic for both models. For the next example only one change comparing with the previous example was done. Chapter 5. Dynamics of an engine-mount system 127 5.0E-4 0.20 Figure 5.8: Steady-state responses: — M,D constant, Crankshaft with J^0 ^ J^° , elastic mounts M,D time-dependent. 3.0E-4 Figure 5.9: Steady-state responses: — M,D constant, Crankshaft with ^ J33, viscoelastic mounts M,D time-dependent. Chapter 5. Dynamics of an engine-mount system 128 Some viscoelastic properties were introduced in the mounts, namely, it was assumed that 4^^ 4 = 0.05 n = 0,1,2,3,4 Kstor{™jj) for the mount material, where Kioss{nw) = the imaginary part, and Kator{nu;) is the real part of the complex stiffness of the springs. One can see that difference between these two models was reduced (Fig. 5.9) in comparison with the elastic mount case (Fig. 5.8). As a conclusion one can say at this stage that in many practical cases when the time-dependent components are small in comparison with the constant components in matrices M and D, the difference in steady-state responses between these two models is negligible, .e., one can use M , D-constant model. However in the case of crankshafts with Jg" ^ Jg", or significant contribution of the reciprocating parts (pistons), the use of M , .D-constant model may yield unreliable results. A full parametric analysis of the effect of different input parameters of engine-mount system on the steady-state response can be viewed as future work. 5.4 Optimization problem for an engine-mount system The steady-state response of the engine-mount system (6-DOF model) yields the information which is necessary to evaluate its performance. The determination of the steady-state response for this model was presented in section 5.3, where the forcing function is represented by a Fourier series with a frequency spectrum w i , ^ , - ,w n . Here the optimization problem of an engine-mount system is posed with respect to prop-erties of the mount material, namely, the values of complex Young's modulus at the frequency spectrum u>i,U2, —,un will be the design variables. Usually n = 2 (1st two harmonics) are es-sential for the type of engine which is of interest in this study (multicylinder, slow, or medium speed diesel engines), all other higher harmonics of the excitation will be neglected. Due to a relatively slow passage from the start (0 frequency) to the operating frequency, a slow (or Chapter 5. Dynamics of an engine-mount system 129 medium speed) engine will be considered as passing through a series of steady-state responses at the intermediate frequencies. Of course, for high speed engines this passage can be quite short, and this assumption will have less relevance. It is assumed that the mounts rest on a rigid foundation. As far as constraints are concerned, the maximum allowable amplitude of dynamic displacements (relative to static equilibrium position) at the mounts are prescribed. For example, at several frequencies f2», i = 1,2,.. (which include intermediate frequencies and the operating one) a system of inequality constraints will be given at each mount as follows ||3i(*,fii)|| <*i I M * , " 0 I I < * 2 IM*.nOII<*s (5.32) where xx,x2, and x 3 are respectively X,Y, and Z displacements of the engine at the mounts. The objective of the optimization procedure is the minimization of the transmitted forces through the mounts. It is assumed for the class of problems considered that the resonant frequencies of the engine-mount system lie well below the operating frequency, which implies that minimum damping is required to minimize the transmitted forces at the operating fre-quency. The introduction of damping is necessary to decrease the dynamic displacements in order to satisfy the inequality constraints (5.32). Therefore the required value of minimum damping depends on limiting values Sx, 82, S3 in (5.32). For intermediate frequencies the minimum value of damping (among those values which yield the satisfaction of (5.32)) will be assigned to the mounts as well. The latest assumption is made not from the considera-tion of transmitted forces, but from the consideration of difficulties in creating of high level damping at low frequencies for the mount material. Therefore the optimization problem posed with the objective to minimize transmitted forces through the mounts (at the operating frequency) and with inequality constraints (5.32) is reduced to the problem of determination of the minimum damping at each frequency in the interval [0 - operating frequency]. Minimum damping means here minimum loss modulus Chapter 5. Dynamics of an engine-mount system 130 (imaginary part of the complex modulus). After the optimization problem is solved one can compute the transmitted forces for the optimum set of design variables, and compare (as a tool for checking) them with ones computed for alternative (not optimum) sets. For each j th mount the norm of transmitted force is computed as H^ 'll = Maxtc[0iT]\\ E.M^Zk^W fe=0 where E*=the complex modulus of the mount material, 7j= stiffness factor, z^ e1"*' is the steady-state response elongation (fcth Fourier harmonic) of the mount. As an objective function, it is proposed to calculate the following quantity T(uj) = Y\\FJ\\ (5-33) i=i which will be called the transmitted force factor. In expression (5.33) N is the number of mounts. Now consider this optimization problem in more detail. For the general case of periodic loading with a period T = ^ the real forcing function can be expressed by a complex Fourier series F(t) = Re[JT Fke^] uk = k^ (5.34) fc=o The steady-state response X(t) (in real form) for the case where M, D are time-dependent matrices can be found in closed form according to the procedure presented in the previous section, i.e., X(t) = Re \k=-l where the number of terms / depends on the number of terms n in expression (5.34). The order of the system to determine the unknown components of vectors Xk will be (21 + 1) * 6, because the system is coupled with respect to the unknown vectors Xk. For computations presented in section 5.3 1 = 4, thus the order of the system was 54. Chapter 5. Dynamics of an engine-mount system 131 With the assumption that M, D matrices are constant (taking the constant term of the Fourier series for these matrices) one can determine the steady-state response as X(t) = Re k=0 where the number of terms coincides with the number of terms in (5.34). Note that it is required to solve a system of order 6 with respect to each unknown component Xk. For the optimization routine where steady-state responses are calculated many times (for different values of design variables) the model with constant M, D matrices would require much less of computer time, than the model with M, D time-dependent. From the other side, it is expedient to use the model with M, D time-dependent when one can expect a large difference in the responses between these two models. It can occur when time-dependent components of matrices M and D are significant, for example, when the engine's crankshaft has Jg" / Jg" (see Fig. 5.8), or when contribution of reciprocating parts is large. A 3-D mount (it can be of cubic, or cylindrical shape etc) is modelled as an appropriate combination of 3 bar elements (Fig. 5.10). The complex modulus of the mount material for these bar elements will be prescribed. Here it is assumed that the mount cross-section is symmetric, thus the lateral stiffnesses are equal and described by one bar element k2. If the cross-section is nonsymmetric, then all the derivations below can be without difficulty generalized by introducing a third parameter k3. Note the stiffnesses of these bar elements are h = £.71 k2 = £.72 where 7 l , 72 are calculated, for example, by the finite element method using an elastic model (with real constant Young's modulus and Poisson's ratio) of the mount. This is done by imparting unit displacements to the mount top plane (Fig. 5.10) in each direction, and calculating the stiffnesses in each direction respectively. Then these stiffnesses divided by the Chapter 5. Dynamics of an engine-mount system 132 77777777777777777P e . 7777 Figure 5.10: Mount model value of elastic modulus of the model yield the stiffness factors'71, 72 which have dimension of length. Coupling between linear lateral displacement and the moment arising is neglected for this model. The values of E* as function of frequency will be considered as the design variables, specifically the values of the loss modulus at the frequencies ui\ and 2u)\. The real part of the complex modulus is computed from the static constraint condition, e.g., if the allowable static displacement of the engine (mass Me) is A, then the storage modulus will be Meg £ 1 (5.35) where N = number of the mounts, g is the gravity acceleration. The value of Ex mentioned above will be the optimum one, because it is clear that any increase in its value would shift the natural frequencies (designed to be below the operating frequency) closer to the operating one, and by that providing an increase in the amplitude of steady-state responses. From the other hand Ei cannot decrease with the frequency (such behavior has not been observed in experiments with polymeric materials). Therefore the value of the storage modulus for the Chapter 5. Dynamics of an engine-mount system 133 optimum material is assumed constant (independent on frequency) and determined by (5.35). 5.4.1 Numerical results for the optimization problem The numerical examples below were obtained for a medium speed multicylinder engine. The input parameters of the engine-mount system are presented in Table D. l , (Appendix D, column 2). One can see that the crankshaft with J22 = J j | ° was considered. For this case the model with constant M, D matrices was utilized as a part of the optimization routine in calculations of steady-state responses. Note that the external forcing vector-function was calculated according to (5.23), where the force contributions arising from nonuniform gas pressure distribution were neglected, i.e., the terms with superscript "gas" were dropped. The minimum values of loss modulus E2(wi) and E2(2u>i) (which provide satisfaction of constraints (5.32)) were obtained by a simple incremental selection (starting from values E2(ui) = E2(2wi) = 0) with the increments A E 2 ( ^ i ) = A £ 2 ( 2 u > i ) = O.OlJEi. Such selection was produced for the operating frequency (which was equal 10 Hz for examples below) and for a set of intermediate frequencies in the interval [0,10] Hz. In such manner the optimum profiles of E2(w) were built and for different values of the unbalance radius (p) (hence for different levels of excitation) are shown in Fig. 5.11. The allowable static displacement of the engine was 0.02 m which determined the value of the storage modulus E\ as shown in Fig. 5.11. The limit values for the dynamic steady-state response displacements in (5.32) were assumed to be Si = S2 = S3 — 0.0001 m. The six natural frequencies of the associated undamped system were 1.715, 3.035, 3.524, 5.456, 5.701, 6.461 Hz. A local peak in the region of 1.5-2.0 (Hz) (Fig. 5.11) is explained by the existence of the 1st natural frequency 1.715 Hz. The stiffness factors of the springs (as was mentioned above they have dimension of Chapter 5. Dynamics of an engine-mount system 134 4.00E6 | i i i i | . i . i | i i i i | i i i i ; i i i i | i i i | i i i i | . i i . | i i i . i 1 1 1 1 i 0 1 2 3 4 5 6 7 8 9 10 11 Frequency, Hz Figure 5.11: Profiles of E2{u)) for different unbalance radii p (m): 1 - p = 0.00014, 2 -p = 0.00016 , 3 - p = 0.00018 , 4 - p = 0.00020 length) were prescribed to be 71 = 72 = 1 (TO). Note that with a change of the assumed stiffness factors from 1, to, say, 71 = 72 = a the storage modulus will change respectively Ei -¥ and the values of the loss modulus will be scaled in the same way E2(u>) —y . However with a change of the ratio of the stiffness factors ^ the optimum profile of E2(u>) will change. This will be shown below. The transmitted force factors are shown respectively in Fig. 5.12 by solid lines. The calculated transmitted force factor for the case of absolutely rigid mounts is shown by dashed lines (it is proportional to UJ2). One can see that introduction of viscoelastic mounts yields a significant reduction in the forces transmitted through the mounts. For this example the transmitted force factor at the operating frequency 10 Hz was 3-5 times less than it would be for the case of rigid mounts. The optimum profiles obtained with different stiffness factor ratios are shown in Fig. 5.13 for two cases: 71 = 1, 72 = 0.5 (TO), and 71 = 1, 72 = 1 (m). The unbalance radius was 0.00018 TO for both cases. The six natural frequencies of the associated undamped system Chapter 5. Dynamics of an engine-mount system 135 Figure 5.12: Transmitted force factor T(u>) for different unbalance radii p (TO): 1 -p = 0.00014 , 2 - p = 0.00016 , 3 - p = 0.00018, 4 - p = 0.00020 for the case = 0.5 were 1.584, 2.302, 3.524, 3.866, 4.356, 5.667 Hz. 7 1 The graphs of transmitted force factors are presented in Fig. 5.14. The solid line represents the previous graph from Fig. 5.12 (72/71 = 1). the dash-dot line represents the (72/71 = 0.5) case, and the dashed line represents the rigid mount case. One can see that (72/71 = 0.5) case yields a lesser transmitted force factor at the operating frequency (10 Hz), than the (72/71 = 1) case in about two times. The stiffness factors of the springs 71 and 72 are determined by the shape of the mount and Poisson's ratio of the mount material (assumed constant for this case), were prescribed here. However the selection of the optimum stiffness factor ratio can be undertaken as an analysis task, i.e., in this case the parameter ^ will be allowed to change. As a conclusion to this section one can notice that the calculated optimum profiles of complex Young's modulus will indicate the ideal (required) properties of a viscoelastic material. Chapter 5. Dynamics of an engine-mount system 136 4.00E6 I i i i i | i w W 2.00E6 "5 T3 O E (/) 1.00E6 00 o I ' ' 1 ' I ' 1 • ' I ' 1 ' 1 I 1 0 1 2 3 4 5 6 7 Frequency, Hz 9 10 11 Figure 5.13: Profiles of ^(w) for different stiffness ratios: solid line - 72/71 = 1, long-dash line - 72/71 = 0.5. Unbalance radius - p = 0.00018 m 1 2 3 4 5 6 7 8 Frequency, Hz Figure 5.14: Transmitted force factor T(u>) for different stiffness ratios: solid line -72/71 = 1) dash-dotted line - 72/71 = 0.5, dashed line - rigid mounts. Unbalance ra-dius - p = 0.00018 m \ Chapter 6 Conclusions and future work In this study the application of exponentials for the relaxation kernel of a viscoelastic operator has been presented and substantiated. Approximation techniques which allow the parameters of the relaxation kernel (as a sum of exponentials) to be determined from experimental data have been developed. Such experimental data which can be used for this purpose are relaxation, or creep curves, or complex moduli (given on a set of frequencies) of elastomeric materials. An explanation of how the finite element procedure is used in treatment of viscoelastic systems has been presented. This procedure can be used to build the mass and elastic stiffness matrices of a system. A hereditary (viscoelastic) stiffness matrix-operator is then obtained by replacement of the elastic constants (Lame's coefficients A and G) in the elastic stiffness matrix by the corresponding viscoelastic operators, or by complex moduli (for steady-state response problems). If Poisson's ratio is assumed constant then only Young's modulus needs to be replaced. A numerical procedure for the determination of complex Young's modulus and complex Poisson's ratio of elastomeric isotropic materials has been developed. This procedure uses as input information certain experimental data concerned with a material specimen. The ob-tained complex Poisson's ratios have shown that for some materials such as CR (chloroprene), NBR (acrylonitrile-butadiene rubber) one can consider Poisson's ratio v„ as constant and real (for the frequency range 10 — 250 Hz). For other materials, e.g., PECH (polyepichloro-hydrin), EAR C — 1002 (Isodamp), and DPNR (natural rubber) the dependence of on 137 Chapter 6. Conclusions and future work 138 frequency is quite noticeable for the mentioned frequency range. Therefore, in general, when a viscoelastic material is considered it is not always possible to assume that vm is constant. A procedure of obtaining steady-state response for discrete viscoelastic systems has been demonstrated and some numerical results have been presented. The effect of complex Pois-son's ratio as a frequency-dependent parameter was significant in some cases, and insignificant in others. The comparison of the numerical results (in terms of steady-state responses) with the ex-perimental ones, obtained for a model of a vibration rig with viscoelastic mounts has demon-strated a good correlation. The theoretical model for this example utilized the complex modulus data calculated through the above mentioned procedure. Finite element formulations for dynamic viscoelastic systems have traditionally been de-scribed in the literature with step-by-step numerical integration (in time) schemes. In this study the closed form solutions in the time-domain for both initial and boundary value prob-lems have been presented for the case when the relaxation kernels of E (Young's modulus operator) and v (Poisson's ratio operator) are represented as a sum of exponentials. Namely, two methods of obtaining closed form solutions in the time-domain by using Laplace transform method and substitution method have been shown. A new approach (substitution method) avoids the difficulties encountered in the use of the Laplace transform approach for multi-degree-of-freedom systems and involves the formulation of an eigenvalue problem. Forced vibration solutions for the case of a periodic forcing function have been also obtained in closed form with the use of the substitution method. An eigenvalue problem for a viscoelastic system has been formulated. The analysis of eigenvalues and eigenvectors has shown that there will be a significant number of overdamped eigenvectors, which correspond to real eigenvalues. The contribution of certain eigenvectors to the homogeneous solution depends essentially on the initial conditions. It is expedient to use the substitution method for the boundary value problems and also Chapter 6. Conclusions and future work 139 in situations when it is advantageous to have the solution in analytical form. It may be noted that the use of numerical integration for the boundary value problem is much more complicated than for the initial value problem. The analysis of conditions of diagonalization of discrete viscoelastic systems has demon-strated that when only Young's modulus operator is introduced (Poisson's ratio is assumed constant) for a system of one homogeneous material and with the proportional viscous damp-ing term, then the equation of motion (in matrix form) is always diagonalizable. If Poisson's ratio is considered as a viscoelastic operator, or the system consists of different viscoelastic materials (even with constant Poisson's ratios), then decoupling of equations of motion, in general, is not possible. An additional investigation is required for this case. An analysis of the free vibration response of a SDOF viscoelastic system (when the relax-ation kernel was assumed to consist of one exponential term) has been presented. In general the response of this system consists of a damped oscillation mode plus an overdamped mode. However there exists a region in the parameter space for which the response consists of three overdamped modes. This region has been analytically defined and expressed in terms of nondimensional parameters characterizing the properties of the material and the inertia property of the system. It has been shown that such a overdamped response is not possible for values - of relaxation time a 1 T a 2ir\/3 where T is the period of the undamped elastic system. In the underdamped region the constitutive parameters can be assigned in such a way that yields prescribed frequencies and decay rates. A formulation of a dynamic model of an engine on viscoelastic mounts has been presented. The equation of motion contains the time-periodic mass and velocity matrices which may lead (under certain conditions) to such phenomena as parametric resonance (dynamic instability). Chapter 6. Conclusions and future work 140 The investigation of the parametric resonance conditions on an example of engine-elastic mount system has shown that in the case of a specific crankshaft (when certain diagonal components of the tensor of moments of inertia are nonequal) and combination of a certain rotation frequency the instability can occur. In general, one can say that when time-dependent components in the mass and velocity matrices become larger then more possibilities exist for instability to take place. A full parametric analysis of the effect of different input parameters on the system's stablitiy, or instablity can be viewed as future work. A steady-state response calculation method has been developed for the case of time-dependent (periodic) matrices in the equation of motion. This method has been applied to an example of engine-mount system and the results have demonstrated that under certain combination of input parameters one can have a significant difference between the responses of two models (one - when mass, velocity matrices are constant, and 2nd - when mass, velocity matrices are time-periodic). Note that upon introduction of damping in the mounts the difference in the responses between these two models is reduced. A full parametric analysis of the effect of different input parameters of the engine-mount system in the creation of such difference can be also viewed as future work. An optimization problem has been posed and solved with the associated constraints and the objective function reflecting the optimum criteria of the performance of the engine-mount system. As the result the optimum parameters of the mount material (storage and loss moduli) have been determined. It was obtained that the transmitted force factor can be reduced several times in comparison with the rigid mount case. It was also discovered that depending on the unbalance radius (level of excitation) there will be different optimum profiles for the loss modulus as function of frequency. A full parametric analysis of the effect of different parameters of an engine-mount system on the optimum profile of the loss modulus can be viewed as future work. Bibliography [1] Turner A . E . , "The use of damping materials for noise reduction of a passenger ship," Journal of sound and vibration, Vol. 10, 1966, pp. 187. [2] Nashif A.D. , Jones D.I.G., Henderson J.P., Vibration damping, McGraw-Hill Book Company, 1993. [3] Tsai C.S., Lee H.H. , "Application of viscoelastic dampers to high-rise buildings," Journal of Structural Engineering, Vol. 119, No. 4, 1993, pp. 1222-1233. [4] McTavish D.J . , Hughes P.C., "Modeling of linear viscoelastic space structures," Journal of Vibration and Acoustics, Vol. 115, 1993, pp. 103-110. [5] Mancuso G. , Sacchi F . , "Main propulsion diesel generator sets with acoustic en-closure and double resilient mounting for low noise application," in "Shipboard Acoustics", Martinus Nijhoff Publishers, 1986, pp.335-350. [6] Christensen R . M . , Theory of viscoelasticity, Academic Press, 1982. [7] Ferry J.D., Viscoelastic properties of polymers, John Wiley and sons, Inc, 1970. [8] Bert C.W., "Material damping: an introductory review of mathematical models, measures and experimental techniques", Journal of Sound and Vibration, Vol. 29, n. 2, 1973, pp. 129-153. [9] Meirovitch L . , Analytical methods in vibrations, The Macmillan company, 1967. [10] Meirovitch, L . , Computational Methods in Structural Dynamics, Sijthoff and No-ordhoff International Publishers, 1980. [11] Weaver, W., Young, D. H . , Timoshenko, S. P., Vibration Problems in Engineering, John Wiley and Sons, 1990. [12] Snowdown J .C. , Vibration and shock in damped mechanical systems, John Wiley and Sons, 1968. [13] Caughey T . K . , M . E . J . O'Kelly, "Classical Normal Modes in Damped Linear Dynamic Systems," Journal of Applied Mechanics, Vol. 12, 1965, pp. 583-588. [14] Thomson, W. T. , Theory of Vibrations with Applications, Prentice Hall Inc., 1993. 141 Bibliography 142 [15] Rao S.S., Mechanical vibrations, Addison-Wesley, 1995. [16] Harris C M . , Crede C . E . , Shock and vibration handbook, McGraw-Hill Book Com-pany, 1961. [17] Hurty, W. C , "Dynamic Analysis of Structural Systems Using Component Modes," AIAA Journal, Vol. 3, April 1965, pp. 678-685. [18] Craig R.R., Ni Z., "Component mode synthesis for model order reduction of non-classically damped systems," Journal of Guidance, Vol. 12, July-August, 1989, pp. 577-584. [19] Muravyov A. , Hutton S.G., "Component mode synthesis for nonclassically damped systems," AIAA Journal, Vol. 34, No. 8, 1996, pp. 1664-1669. [20] Muravyov A. , "Analysis of damped linear dynamic systems and application of component mode synthesis," M.A.Sc. thesis, Mechanical Engineering Department, University of British Columbia, Vancouver, B .C . , 1994. [21] Muravyov A. , Hutton S.G., "A free-interface CMS method for systems with gen-eral damping," Proceedings of the 13th International Modal Analysis Conference, Nashville, USA, 1995, pp. 255-260. [22] Muravyov A. , Hutton S.G.,"A geometrical interpretation of damping for discrete classically damped systems", International Journal of Mechanical Sciences, (in press). [23] Muravyov A. , Hutton S.G., "Free vibration response of discrete systems with gen-eral damping characteristics," Applied Mechanics in the Americas, Vol. 2, Godoy L . A . , Idelsohn S.R., Laura P.A., and Mook D.T. (eds), Santa Fe: A A M and A M C A , 1995, pp.121-126. [24] Rabotnov Yu.,N., Creep problems in structural members, North-Holland Publishing Company, Amsterdam, 1969. [25] Bland D.R., The theory of linear viscoelasticity, Pergamon Press, Oxford, 1960. [26] Creus G.J. , Viscoelasticity - Basic Theory and Applications to concrete structures, Springer-Verlag, 1986. [27] Flugge W., Viscoelasticity, Blaisdell Publishing Company, 1967. [28] Gross B. , Mathematical structure of the theories of viscoelasticity, Hermann & C, 1953. Bibliography 143 [29] Caxini A. , De Donato 0., "Fundamental solutions for linear viscoelastic continua", Journal of Solids and Stuctures, Vol. 29, No. 23, 1992, pp. 2989-3009. [30] Rabotnov Yu.,N., Elements of hereditary solid mechamics, Mir Publishers, Moscow, 1980. [31] Dlubac J .J . , Lee G.F. , Duffy J .V. , Deigan R.J . , Lee J.D., "Comparison of the complex dynamic modulus as measured by three apparatus," in Sound and vibration damping with polymers, American Chemical Society, Washington, D C , 1990, pp. 49-62. [32] Sattinger S.S., "Direct method for measuring the dynamic shear properties of damping polymers," in Sound and vibration damping with polymers, American Chemical Society, Washington, D C , 1990, pp. 79-91. [33] Andrews R.D., Hofman-Bang N. , Tobolsky A . V . , "Elastoviscous Properties of poly-isobutylene. Relaxation of stress in whole Polymer of different molecular weights at elevated temperatures", J. Polymer Sci., Vol. 3, No. 5, 1948, pp. 669-682. [34] Burns J . , Dubbleday P.S., Ting R.Y. , "Dynamic bulk modulus of various elas-tomers," Journal of polymer science, Part B: Polymer Physics Vol. 28, 1990, pp. 1187-1195. [35] Gottenberg W . G . , Christensen R . M . , "An experiment of the mechanical property in shear for a linear, isotropic viscoelastic solid," International Journal of Engineering Sciences, Vol. 2, 1964, pp. 45-57. [36] Oyadiji S.O., Tomlinson G.R., "Relating the complex moduli of viscoelastic ma-terials to the complex stiffness of antivibration mounts," Proceedings of the Smart Structures and Materials conference, February, Florida, USA, 1994, pp.226-237. [37] Plazek D.J . , Vrancken M.N. , John W.B. , "A torsion pendulum for dynamic and creep measurements of soft viscoelastic materials," Trans. Soc. Rheology, Vol. 2, 1951, pp. 39-51. [38] Read B., Dean G.D., The determination of dynamic properties of polymers and composites, Adam Hilger Ltd, 1978. [39] Reid M . C . , Oyadiji S.O., Tomlinson G.R., "Relating material properties and wave effects in vibration isolators," Proceedings of the 59th Shock & Vibration Sympo-sium, Sandia Laboratories, USA, October, 1988, pp. 139-156. [40] Chen C P . , "Characterization of and design with viscoelastic materials," Ph.D the-sis, Mechanical Engineering Department, The University of Iowa, 1989. Bibliography 144 [41] Hartmann B., Lee G.F. , Lee J.D., "Loss factor height and width limits for polymer relaxations," Journal of Acoustic Society, Vol. 95, No. 1, 1994. pp. 226-233. [42] Havriliak Jr., Negami S., "A complex plane representation of dielectric and mechan-ical relaxation processes in some polymers," Polymer, Vol. 8, 1967, pp. 161-210. [43] Bagley R.L. , Torvik P. J . , "Fractional calculus - a different approach to the analysis of viscoelastically damped systems," AIAA Journal, Vol. 21, No. 5, 1983, pp. 741-748. [44] Chern J .T . , "Finite element modeling of viscoelastic materials on the theory of frac-tional calculus," Ph.D thesis, Mechanical Engineering Department, Pennsylvania State University, 1993. [45] Suarez L . E . , Shokooh A. , "On the response of systems with damping materials modeled using fractional calculus," Applied Mechanics in the Americas 2, Godoy L . A . , Idelsohn S.R., Laura P.A., and Mook D.T. (eds), Santa Fe: A A M and A M C A , 1995, pp.147-151. [46] Miller K.S. , Ross B. , An introduction to the fractional calculus and fractional dif-ferential equations, John Wiley and Sons, Inc., 1993. [47] Gradshtein I.S., Ryzhik I.M., Table of integrals, series, and products, Academic Press, 1980. [48] Gromov V . G . , "On the mathematical content of Volterra's principle in the bound-ary value problem of viscoelasticity", Prikl. Matem. i Mekh., Vol. 35, No. 5, 1971, pp. 869-878. [49] "ABAQUS: user's manual", Version 5.3, Hibbit, Earlsson and Sorensen, 1993. [50] Szumski R . G . , "A finite element formulation for the time domain vibration analysis of an elastic-viscoelastic structure," Ph.D thesis, Mechanical Engineering Depart-ment, Georgia Institute of Technology, 1993. [51] Day S.M., Minster J.B. , "Numerical simulation of attenuated wavefields using a Pade approximant method," Geophys. JR. astr. Soc, Vol. 78, 1984, pp. 105-118. [52] Brunner H . , Von der Houwen P.J., The numerical solution of Volterra equations, North-HoUand, 1986. [53] Linz P., Analytical and numerical methods for Volterra equations, SIAM, Philadel-phia, 1985. Bibliography 145 [54] Muravyov A. , "Analytical solutions in the time domain for vibration problems of discrete viscoelastic systems", Journal of Sound and Vibration, Vol. 199, n. 2, 1997, pp. 337-348. [55] Muravyov A. , Hutton S.G., "Closed form solutions and the eigenvalue problem for vibration of discrete viscoelastic systems", Journal of Applied Mechanics, (in press). [56] Muravyov A. , Hutton S.G., "Free vibration of discrete viscoelastic systems and the eigenvalue problem," Proceedings of the 1st South African conference of Applied Mechanics, Eskom Midrand, South Africa, July, 1996, pp. 215-222. [57] Golla D.F . , Hughes P.C., "Dynamics of viscoelastic structures-a time domain, finite element formulation," Journal of Applied Mechanics, Vol. 52, 1985, pp. 897-906. [58] Johnson C D . , Kienholz D.A. , "Finite element prediction of damping in structures with constrained viscoelastic layers," AIAA Journal, Vol. 20, No. 9., 1982, pp. 1284-1290. [59] Booker J.R., Small J . C , "Finite layer analysis of layered viscoelastic materials under three-dimensional loading conditions," International Journal for numerical methods in engineering, Vol. 21, 1985, pp. 1709-1727. [60] Y i S., Fouad Ahmad M . , Hilton H.H. , "Dynamic responses of plates with vis-coelastic free layer damping treatment," Journal of vibration and acoustics, Vol. 118, July, 1996, pp. 362-367. [61] Petrovsky N. , Marine internal combustion engines, Mir Publishers, Moscow, 1960. [62] Henshall S.H., Medium and high speed diesel engines for marine use, St. Stephen's Bristol Press, 1972. [63] Swanson D.A. , Miller L.R., Norris M . A . , "Multidimensional mount effectiveness for vibration isolation," Journal of Aircraft, Vol. 31, No. 1, 1994, pp. 188-196. [64] Swanson D.A. , Wu H.T. , Ashrafiuon H. , "Optimization of aircraft engine suspen-sion systems," Journal of aircraft, Vol. 30, No. 6, 1993, pp. 979-984. [65] Ashrafiuon H . , "Design optimization of aircraft engine-mount systems ," Journal of vibration and acoustics, Vol. 115, Oct., 1993, pp. 463-467. [66] Spiekermann C . E . , Radcliffe C.J . , Goodman E.D. , "Optimal design and simula-tion of vibrational isolation systems," Journal of mechanisms, transmissions, and automation in design, Vol. 107, June, 1985, pp. 271-276. Bibliography 146 [67] Pshenichny B.N. , Danilin Y . M . , Numerical methods in extremal problems, Mir Pub-lishers, Moscow, 1978. [68] Beltzer A.I., Acoustics of solids, Springer-Verlag, 1988. [69] Sturla F . A . , Argento A. , "Free and forced vibrations of a spinning viscoelastic beam," Journal of vibration and acoustics, Vol. 118, July, 1996, pp. 463-468. [70] Struik L . C . , "Free damped vibrations of linear viscoelastic materials," Rhelogica Acta, Vol. 6, No. 2, 1967, pp. 119-129. [71] Vibration and Strength Analysis Program (VAST): User's manual, Version 6.0, Martec Limited, Halifax, Nova Scotia, 1990. [72] Lusternik L . A . , Sobolev V . J . , Elements of functional analysis, Hindustan Publish-ing Corpn., John Wiley and sons, Inc., 1974. [73] Chen G. , Hutton S.G., "An experimental research on the transmissibility of vis-coelastic material mounts," research report (in preparation), Mechanical Engineer-ing Department, University of British Columbia, Vancouver B . C . [74] Bugrov Ya.S., Nikolski S.M., Differential equations, Multiple integrals, Series, Functions of a complex variable, Moscow, Nauka, (in russian), 1989. [75] Nicol T., (editor), UBC Matrix book (A Guide to Solving Matrix Problems), Com-puting Centre, University of British Columbia, Vancouver, B . C . , March, 1982. [76] Noble B. , Applied linear algebra,, Prentice-Hall, Inc., 1969. [77] Foss K . A . , "Coordinates which uncouple the equations of motion of damped linear dynamic systems," Journal of Applied Mechanics," September, 1958, pp. 361-364. [78] Bellman R., Introduction to Matrix Analysis, McGraw-Hill Book Company Inc., 1960. [79] Muravyov A. , Hutton S.G., "Free vibration response characteristics of a simple elasto-hereditary system", Journal of Vibration and Acoustics, (in press). [80] Korn G.A. , Korn T . M . , Mathematical handbook for scientists and engineers, McGraw-Hill Book Company, 1968. [81] Lourie A.,I., Analytical mechanics, Moscow, 1961 (in russian). [82] Muravyov A. , Hutton S.G., "Analysis of an engine-mount system with the time-dependent mass and velocity matrices", submitted to Journal of Sound and Vibra-tion, 1997. Bibliography 147 [83] Iakubovich V . A . , Starzhinskii V . M . , Linear differential equations with periodic co-efficients, John Wiley and sons, 1975. [84] Fomin V . N . , Mathematical theory of parametric resonance in linear consistent sys-tems, Leningrad University, 1972 (in russian). [85] Stevens K . K . , Evan-Iwanowski R . M . , "Parametric resonance of viscoelastic columns," International Journal of Solids and Structures," Vol. 5, pp. 1969, 755-765. [86] Nayfeh A . H . , Mook D.T. , Nonlinear oscillations, John Wiley and sons, 1979. [87] Stoker J .J . , Nonlinear vibrations, New York, 1950. Appendix A Transfer functions, complex moduli The description of the procedure for determination of complex moduli has been presented in section 3.2. Below the complex transfer functions experimentally obtained for each specimen of material in the two tests are represented by the real and imaginary parts. Calculated complex Young's and shear moduli are represented by the real part (N/m2) and the loss factor. Complex Poisson's ratios are represented by the real and imaginary parts. 6 i i i ' ' | i i i i | i i i i | ij 4 h • 6 h 40 80 120 160 200 240 Frequency, Hz Figure A . l : Transfer function in the 1st test for P E C H 148 Appendix A. Transfer functions, complex moduli 149 2.0 i . i i i | 1.5 --1.5 h 40 80 120 i60 200 240 Frequency, Hz Figure A.3: Transfer function in the 1st test for E A R C-1002 Appendix A. Transfer functions, complex moduli 150 40 80 120 160 200 240 Frequency, Hz Figure A.5: Transfer function in the 1st test for DPNR Appendix A. Transfer functions, complex moduli 151 80 120 160 Frequency, Hz Figure A.6: Transfer function in the 2nd test for DPNR Appendix A. Transfer functions, complex moduli 152 c o o c CD CO c CO I magi nary part 40 200 240 80 120 160 Frequency, Hz Figure A.8: Transfer function in the 2nd test for CR Real part Imaginary part 40 200 240 80 120 160 Frequency, Hz Figure A.9: Transfer function in the 1st test for NBR Appendix A. Transfer functions, complex moduli Appendix A. Transfer functions, complex moduli 154 O.OOEO 1 — ' — ' — 1 — 1 — 1 — ' — ' — 1 — ' — 1 — ' — 1 — 1 — ' — 1 — ' — 1 — 1 — ' — 1 — 1 — 1 — 1 — 1 — 1 0 50 100 150 200 250 Frequency, Hz Figure A.11: Complex Young's modulus of P E C H Frequency, Hz Figure A.12: Complex Poisson's ratio of P E C H Appendix A. Transfer functions, complex moduli 155 Appendix A. Transfer functions, complex moduli 156 8.00E6 r — i — i — i — i — | — i — i — i — i — | — i — i — i — i — | — i — i — i — i — | — i — i — < — ' — I 3.5 0 50 100 150 200 250 Frequency, Hz Figure A.16: Complex shear modulus of E A R C-1002 Appendix A. Transfer functions, complex moduli 157 O.OOEO 1 — 1 — 1 — ' — 1 — 1 — ' — ' — 1 — ' — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 0.05 0 50 100 150 200 250 Frequency, Hz Figure A.17: Complex Young's modulus of DPNR Frequency, Hz Figure A.18: Complex Poisson's ratio of DPNR Appendix A. Transfer functions, complex moduli 158 1.50E6 CO 1.00E6 Q . CO 0 DC 5.00E5 0.00E0 Real part Loss factor 50 200 100 150 Frequency, Hz Figure A.19: Complex shear modulus of DPNR I 0.40 0.30 o -10.20 3 CO CO o 0.10 1 0.00 250 1.20E7 1.00E7 l-8.00E6 co 3" 6.00E6 CO CD rr 4.00E6 2.00E6 0.00E0 I 0.6 Real part 50 Loss factor 0.4 o | 0.3 & —-A o I 0.2 0.1 1 0.0 100 200 150 Frequency, Hz Figure A.20: Complex Young's modulus of CR 250 Appendix A. Transfer functions, complex moduli 159 O.OEO Frequency, Hz Figure A.21: Complex Poisson's ratio of CR Appendix A. Transfer functions, complex moduli 160 1.20E7 i — i — i — i — i — | — i — i — i — i — i — i — i — i — i — i — i — . — i — i — | — i — i — i — i — i 0.6 0.00E0 I—'—'—1—'—1—1—'—1—1—1—'—'—'—'—1—'—'—'—1—1—1—1—1—1—' 0.0 0 50 100 150 200 250 Frequency, Hz Figure A.23: Complex Young's modulus of NBR Real part 50 150 200 100 Frequency, Hz Figure A.24: Complex Poisson's ratio of NBR 250 Appendix A. Transfer functions, complex moduli 4.00E6 3.50E6 \-O.OOEO 1 — ' — 1 — ' — ' — 1 — ' — ' — ' — ' — 1 — 1 — 1 — 1 — 1 — 1 — ' — 1 — 1 — 1 — 1 — 1 — ' -0 50 100 150 200 Frequency, Hz Figure A.25: Complex shear modulus of NBR Appendix B Dimensions of the vibration rig -7-7-. / / / / / / / / ',1 / / / / / y / / / 17 0.25 -— 0.45 " 1 • • - Mount 14 12 36 Figure B.l: Vibration rig with viscoelastic mounts, dimensions in [inch] 162 Appendix B. Dimensions of the vibration rig A - A (73 oi 2 X 2 X 1 / 8 angle Figure B.2: Vibration rig, View A - A Appendix C Derivation of the transformation matrix R (for section 5.1) The vectors i x , i 2 , i 3 of the auxiliary rotating basis (Fig. 5.1) are expressed in terms of the ground basis vectors as follows n r 1 0 0 0 cosf sin-y e2 (Cl) 0 —siwy cos~f where 7 = cut + cp0 (angle of crankshaft rotation). Assume that at the moment t (or for the angle 7) the engine framework has the following angular displacements </>i,<j>2,<j>z, then the basis g i , g 2 , g 3 (rigidly fixed with the crankshaft) is actually the basis i i , i 2 , i 3 which is given the angular displacements </>i, (f>2, cfiz, i.e., • i i 12 = 13 ei e2 e3 gfc = U + £ A i f e n * = 1,2,3 (C.2) n = l where AU„ is the change to the i*. basis vector due to its rotation about the basis vector en. These changes are computed as follows (note that angles <j>i,<j)2,4>z are assumed very small): Ain = </>iei x i x = ei e 2 e 3 4>! 0 0 1 0 0 analogously A i i 2 = </>2e2 x ii ei e2 e 3 0 c/>2 0 1 0 0 = -</>2e3 164 Appendix C. Derivation of the transformation matrix R (for section 5.1) Aii3 = 0 3 e 3 x ii = ei e 2 e 3 0 0 0 3 1 0 0 03e2 A i 2 i = 0ie x x i 2 ei e 2 e 3 0i 0 0 0 cos-y siwy = 0 i C O 5 7 e 3 — <f>isin/ye2 A i 2 2 = 0 2 e 2 x i 2 A i 2 3 = 0 3 e 3 x i 2 = ei e 2 e 3 0 0 2 0 0 cos~j siwy ei e 2 e 3 0 0 0 3 0 cosj siwy = <f>2siwyei = -0 3co57ei A i 3 i = 0iei x i 3 = ei e 2 0 —siwy cos-y e 3 0 —(j>\sinyez — 0 i C O 5 7 e 2 ei e 2 e 3 A i 3 2 = 0 2 e 2 x i 3 = 0 02 0 = 0 2 c o 5 7 e i 0 —siwy COSJ ei e 2 e 3 A i 3 3 = 03 e 3 x i 3 = 0 0 03 = (j)3siwye\ 0 —siwy COSJ Appendix C. Derivation of the transformation matrix R (for section 5.1) Therefore substituting the obtained Ai**. in (C.2) one can rewrite it in matrix form 0 03 • gl ll g2 = 12 + . S 3 . is -02 ei (j)\COSr) e2 —<f>isiwy e3 and substituting ( C l ) in the expression above one obtains gl 1 03 -02 ei g2 = fasiwy — cpscosj C057 — (frisinj siwy + <j>iCos~f e2 . S3 . fesin'y + facosj —siwy — 0ico«7 C057 — cpisiwy e3 or in abbreviated form gi = Rikek i = 1,2,3 where matrix R will be called the transformation matrix. Appendix D Input parameters of an engine-mount system The parameters below are mentioned in section 5.1. p = unbalance radius of the crankshaft (the distance between the mass centre of the crankshaft and the axis of rotation) m-A = mass of the engine framework TUB = mass of the crankshaft si, s2, s3 = components of the vector YAD (see Fig. 5.2), where point A = mass centre of the engine framework and point D = projection point of the crankshaft mass centre on the axis of rotation I = length of the crank (see Fig. 5.3, denoted as "/") h — length of the connecting rod (see Fig. 5.3, denoted as "Z2") Ncyl = number of cylinders mp = piston's mass fi = distance between ith cylinder and point D along the axis of rotation (see Fig. 5.3, denoted as "/") 7Poi = phase angle of ith piston «/4 = components of tensor of moments of inertia of the engine framework Jff = components of tensor of moments of inertia of the crankshaft (in the system rigidly connected with the crankshaft) J41, A2 = dimensions of the engine plate which sits on the mounts (see Fig. 5.4) H= vertical distance between e.g. of the engine framework and horizontal plane which is 167 Appendix D. Input parameters of an engine-mount system 168 attached to the mounts k*(nu) = complex stiffness of the springs (Fig. 5.4) at frequency nu>, n = 0,1,2,3,4. Below the values of these parameters are shown in Table D. l , where column 1 contains input parameters for the parametric resonance example, column 2 - input parameters for the steady-state response and optimization problem examples. The values with f are used for the steady-state response example only. For values of p for the optimization problem example see section 5.4.1. For the optimization routine the following additional input parameters are required: A = maximum allowable static displacement (Z direction), Si,S2,S3 = maximum allowable dynamic displacements at the mounts (X,Y,Z directions). Remark. Without loss of generality, at the instant t = 0 it is assumed that unbalance radius has zero phase, i.e., it is oriented as axis Y (vector e2 in Fig. 5.2). The phase angles of the pistons are determined according to the configuration of the cranks and with respect to the unbalance radius direction. The components of J$° are calculated in the coordinate system parallel to the reference system (t = 0) and taken at the crankshaft's e.g. They are so-called components in the fixed with the crankshaft coordinate system. The components of JA are taken at the e.g. (point A) and in the system parallel to the reference system (they are assumed constant). Appendix D. Input parameters of an engine-mount system 169 Parameter 1 2 Parameter 1 2 P - o.ooois* JA 12000. 12000. mA 12181. 12181. J12 10. 10. mB 4116. 4116. JA 13 10. 10. Sl -0.1 -0.1 TA J22 17100. 17100. &2 0.02 0.02 IA J23 20. 20. S3 -0.8 -0.8 TA "33 17100. 17100. I 0.2286 0.2286 TBO J l l 299. 299. h 1.162 1.162 jBo J12 0. 0. Ncyl 6 6 jBo "13 0. o. mp 81. 81. jBo J22 8453. 8453. fi 1.6 1.6 jBo J23 0. 0. h 0.96 0.96 jBo J33 8453. 8453. h 0.32 0.32 Al 0.72 0.72 U -0.32 -0.32 A2 1.8 1.8 h -0.96 -0.96 H 1.2 1.2 U -1.6 -1.6 Mo) - 0.20667+10+ Ipol 0. 0. - 0.20667+10+ lpo2 -60. -120. - 0.20667+10+ lpo3 -120. -240. - 0.206e7+iOf Kpo4 -180. -240. fc,(4u>) - 0.206e7+iOf IpoS -240. -120. - - -lpo6 -300. 0. • - - -Table D . l : Input parameters of the engine-mount system. Units: length - [m], mass -[kg], angle - [degree], Jtj, J?° - [kg * m2], k. - [N/m]
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Discrete dynamic viscoelastic systems and vibration...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Discrete dynamic viscoelastic systems and vibration analysis of an engine supported on viscoelastic mounts Muravyov, Alexander 1997
pdf
Page Metadata
Item Metadata
Title | Discrete dynamic viscoelastic systems and vibration analysis of an engine supported on viscoelastic mounts |
Creator |
Muravyov, Alexander |
Date Issued | 1997 |
Description | An analysis of linear dynamic viscoelastic discrete systems is presented, and its application to the vibration of an engine supported on viscoelastic mounts is discussed. A new procedure is developed which allows exact (closed form) homogeneous solutions in the time domain to be derived for a dynamic system consisting of isotropic viscoelastic components, for which the relaxation kernels are represented as a sum of exponentials. The developed procedure (which is given a name "substitution method") for determination of closed form solutions is extended to the solution of boundary value problems. The application of the substitution method is also extended to the case of periodic loading. Based on this method, a numerical investigation of free and forced vibration responses of some viscoelastic systems is presented. Several approximation techniques are developed in this study which allow the parameters of the relaxation kernels (represented as a sum of exponentials) to be determined from experimental data. Also a numerical procedure for determination of complex moduli of isotropic viscoelastic materials is developed, in which certain experimental data related to the material specimen are required as input information. A hereditary (viscoelastic) stiffness matrix-operator is obtained by replacement of the elastic constants in the elastic stiffness matrix by the corresponding viscoelastic operators, or by complex moduli (for steady-state response problems). Comparison of experimental results (in terms of steady-state responses) with the numerical ones is presented. The sufficient conditions of diagonalization of discrete viscoelastic systems are formulated in this study. Analysis of the conditions of overdamping of a simple viscoelastic single-degree-of-freedom system is conducted, and new results concerned with this analysis are demonstrated. A particular case of a dynamic viscoelastic system (an internal combustion engine on elastomeric mounts) is given special consideration. A new dynamic model of an enginemount system is developed where rotating and reciprocating parts lead to the mass matrix and velocity matrix (matrix-coefficient at the velocity vector) as periodic functions of time. The derivation of the equations of motion on the basis of Lagrange's equations is demonstrated. An analysis of parametric resonance phenomena for some examples of engine-mount systems is conducted. A method for steady-state response calculations for the case of time-dependent (periodic) matrices in the equation of motion is developed and some numerical results are presented. An optimization problem is posed and solved with the associated constraints and the objective function reflecting the optimum criteria of the performance of an engine-mount system. As a result, the optimum parameters of the mount material are determined. |
Extent | 5896884 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-04-20 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0080885 |
URI | http://hdl.handle.net/2429/7435 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mechanical Engineering |
Affiliation |
Applied Science, Faculty of Mechanical Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 1997-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_1997-251241.pdf [ 5.62MB ]
- Metadata
- JSON: 831-1.0080885.json
- JSON-LD: 831-1.0080885-ld.json
- RDF/XML (Pretty): 831-1.0080885-rdf.xml
- RDF/JSON: 831-1.0080885-rdf.json
- Turtle: 831-1.0080885-turtle.txt
- N-Triples: 831-1.0080885-rdf-ntriples.txt
- Original Record: 831-1.0080885-source.json
- Full Text
- 831-1.0080885-fulltext.txt
- Citation
- 831-1.0080885.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0080885/manifest