Frequency domain analysis of continuous systems with viscous generalized damping

This paper deals with the effects of generalized damping distributions on vibrating linear systems. The attention is focused on continuous linear systems with distributed and possibly non-proportional viscous damping, which are studied in terms of modal analysis, defining and discussing the orthogonality properties of their eigenfunctions. Exact expressions of the frequency response functions obtained by direct integration of the equations of motion are compared with the analogous formulas based on the superposition of modes. In addition, approximate expressions of the frequency response functions of both continuous and discrete (finite element models) systems in terms of their undamped eigenfunctions/eigenvectors are also considered and discussed. The presented methods are explained, compared and validated by means of numerical examples on a clamped-free EulerBernoulli beam.


Introduction
The effects of generalized damping distributions on vibrating linear mechanical systems have been not exhaustively studied in terms of modal analysis, especially as regards to continuous systems, i.e. distributed parameter systems.In fact, continuous systems are seldom modelled considering damping distributions, and when this is the case, the models are almost always based on the proportional damping assumption, i.e. the damping operator can be expressed as a linear combination of the mass and stiffness operators [1].This way of modelling is of course so often adopted since it carries little analytical and computational further effort in addition to the undamped case analysis.But in many real situations the proportional damping assumption is not valid and such a simplified approach does not describe the dynamics of the system with sufficient accuracy.So, in this paper the more general case in which the damping distributions result to be non-proportional is considered.
Since the existing bibliography about modal analysis for distributed parameter systems with generalized damping concerns particular cases [2,3], a complete theoretical statement valid in the general case is herein included.It is nothing but the natural extension of the undamped continuous system theory [1], in which the well-known results for the discrete systems are easily found as a particular case.
After reducing the differential boundary problem to an eigenvalue problem through the separation of variables, the existence of orthogonality relations valid for the general case of non-proportional damping is demonstrated.By means of such relations, the general solution can then be expressed through the so-called expansion theorem.A discussion about the meaning of the modal parameters in case of non-proportional damping is also included.
Since the solution of the above mentioned differential eigenproblem is often a very difficult task, several methods giving approximate solutions have been proposed [1].Among them, this paper includes the description and application of a technique presented in [4], valid when the solution is known for the undamped system only.
For a rather large class of continuous systems, however, analytical solutions (i.e.analytical expressions of the eigenfunctions) can be found.According to [5], a method for the solution of the differential eigenproblem is described, suitable for a class of vibrating continuous systems with non-proportional damping distributions, according to different damping models.
This method, starting from a partition of a continuous system in homogeneous substructures or segments, matches a reduction of the differential equation order with a transfer matrix technique.So it can be easily applied to a large number of continuous vibrating systems with non-proportional damping, provided that closed-form solutions of the undamped case for each substructure are known.Moreover, this approach leads to an easy computer implementation and presents a high computational efficiency, due to the invariance of the matrix dimensions with respect to the number of substructures considered.
The particular case of non-homogeneous Euler-Bernoulli beams with different non-proportional internal and external damping distributions is considered.Actually, only little changes in the coefficients would be required to solve the problem for strings, rods, shafts or Timoshenko beams with viscous or more complicated damping models.The possibility of extending the proposed method also to several kinds of membranes and plates is intrinsic in its formulation.
The attention is then focused on the expression of the frequency response functions (FRFs).This result can now be achieved through the modal analysis approach, since both the modal parameters and modal shapes are available applying the above mentioned analytical methods.The same techniques, however, allow to express the FRFs in a different way, which does not require either eigenvalues or eigenfunctions, i.e. which does not need the solution of any eigenproblem.This result is simply achieved by direct integration of the equations of motion, and since in this case the solution is not expressed by means of a series (indeed it is the sum of the series resulting from the modal approach), it could be very useful for high frequency analysis.
Numerical examples are then included in order to show and to compare both the accuracy and efficiency of the proposed methods.Non-proportional damping distributions are tested on a non-homogeneous Euler-Bernoulli beam in bending vibration and consequently a discussion on the related frequency response diagrams is presented.
Finally, the results are validated by means of a finite element (FE) model, thus showing their reliability in problems involving non-proportional damping distributions.

Modal analysis of continuous systems with viscous generalized damping
In this section some fundamentals of modal analysis for distributed parameter systems are presented.Among them, a statement of the expansion theorem, leading to the general solution, a discussion about the meaning of the modal parameters and finally the expression of the frequency response functions, both in exact and approximate form.

General solution: The expansion theorem
The dynamic behaviour of a continuous system with viscous generalized damping can be described by the following equation of motion where M, C, K are linear homogeneous differential operators and are referred to as mass operator, damping operator (generally non-proportional) and stiffness operator, respectively, f is the external force density, w and x are the displacement and the spatial coordinate in a domain of extension D, respectively (the spatial coordinate as well as the displacement and the external force density can be scalars or vectors, but here the scalar notation is adopted, since in what follows this does not represent a loss of generality) and t is time.
To solve Eq. ( 1), appropriate boundary and initial conditions must be satisfied by w.
If the damping operator can be expressed as a linear combination of the mass and stiffness operators, it is said to be proportional [1], but in this paper the more general case in which it results to be non-proportional is considered.
Recalling that self-adjointness of differential operators in continuous systems corresponds to symmetry of matrices in discrete systems, in the following the operators M, C and K will be supposed to be self-adjoint.This assumption is unnecessary and could easily be removed, as it will be clarified in the end of the present section.Nevertheless it will be adopted since such properties hold in most of the existing models, carrying less cumbersome analytical developments.
With the notation denoting the inner product between two scalar functions w 1 and w 2 over a domain of extension D (if w 1 and w 2 are vector functions, then the integrand in Eq. ( 2) represents their scalar product), a linear differential operator L is said to be self-adjoint if which is a property of symmetry with respect to the inner product.
As usual in modal analysis, the differential boundary problem will be reduced to a differential eigenvalue problem by separating the variables, so the solution will be given by a linear combination of terms in the form w(x, t) = φ(x)q(t).Note that although the global solution will always be real, complex terms φ(x)q(t) are expected: in the case of non-proportional damping the phase of φ(x) will not be constant with respect to x, and consequently the motion will be non-synchronous [5].
Equation (1) can then be rewritten in the state-space form as follows the dot denoting derivation with respect to time, where A and B are linear homogeneous differential operators, w and f are the state vector and the external force density vector, respectively.They can be expressed as so that, if M, C and K are self-adjoint, A and B result to be self-adjoint as well.
The state-space Eq. (4) leads to the differential eigenproblem The solution of this eigenproblem forms an infinite set of pairs of discrete values, each pair characterizing a mode and being related to a pair of eigenvectors (i.e. to a pair of eigenfunctions).In the case of underdamped modes, complex conjugate pairs s n , s * n of eigenvalues are expected, yielding pairs of complex conjugate eigenvectors z n , z * n (i.e.pairs of complex conjugate eigenfunctions φ n , φ * n ).On the other hand, in the case of overdamped modes, pairs of real negative eigenvalues are expected, yielding pairs of real eigenvectors (in the following, the notation for the underdamped case will be adopted, but the results are valid in both cases).
The eigenvectors orthogonality properties can be derived rewriting Eq. ( 6) for the nth and mth mode, respectively, pre-multiplying the first by z T m and the second by z T n , then integrating them both over the spatial domain D, i.e.
which, taking into account the self-adjointness of A and B, yield Due to the orthogonality properties of the eigenvectors z n , any other vector in the same space of functions can be expressed as their linear combination.This statement constitutes what is usually known as the expansion theorem, so that the free response can be written in the form where γ n is a complex scaling factor which depends on the initial conditions.Note that a n and b n , if taken separately, are known unless an indeterminate scaling factor albeit their ratio is fixed, i.e. b n /a n = −s n .
When at least one of the differential operators involved in the model results not to be self-adjoint the expansion theorem still holds, but the orthogonality relations Eq. ( 8) have to be replaced by a set of biorthogonality relations, which require the solution of two eigenproblems for the so-called right and left eigenvectors, respectively [1].

Expression of the frequency response functions
A harmonic excitation force of amplitude f 0 acting with angular frequency ω at a coordinate x f is now considered.Since the system is linear-time-invariant, the response w will still be a harmonic oscillation at the same angular frequency ω.So, taking into account the expansion theorem and dropping the time dependent terms, the state-space equation of motion Eq. ( 4) can be rewritten as where Γ m is a scaling factor and Pre-multiplying by z T n , integrating over the spatial domain D and remembering the orthogonality properties Eq. ( 8), Eq. ( 10) gives where the expression of the modal force f n = (z n , f 0 ) in terms of φ n and f 0 is due to the Dirac distribution properties.
Finally, if the eigenfunctions are normalized with respect to a n , i.e.
taking into account again the expansion theorem and Eqs ( 9) to (11), it is possible to express the system receptance, defined as the ratio of the amplitude of displacement at a coordinate x to the intensity of a single harmonic force acting at a coordinate x f as follows The expressions of other frequency response functions, such as accelerance or mechanical impedance, follow immediately from Eq. (13).

Meaning of the modal parameters
The definitions of the modal parameters which hold in the case of proportional damping, usually referred to as modal mass, modal damping and modal stiffness, can be extended to the non-proportional case according to Despite their dimensions are coincident with those of a modal mass, a modal damping and a modal stiffness, respectively, their properties are not the same.To highlight this concept, it is necessary to put in explicit form their relationships with the eigenvalues s n .
Remembering the previous definitions, the orthogonality relations (8a) yield Since these expressions are formally the same as in the case of proportional damping, then by analogy it is possible to define represents the modulus of the related eigenvalue s n and ζ (np) n defines its phase, but it is very important to underline that ω (np) n is not the modal natural angular frequency ω n and its magnitude depends on the rate of damping, whilst ζ (np) n is not the modal damping ratio ζ n which holds in the case of proportional damping.As a consequence, the n th modal natural frequency cannot be extracted directly from the corresponding eigenvalue s n .
The gap between ω (np) n and ω n gives a measure of the non-proportionality effects, therefore it suggests the definition of the following "modal" index of non-proportionality which is a dimensionless, real, non-negative parameter.For different definitions of non-proportionality indexes the reader is referred to [6].

Approximation of the solution
The solution of the eigenproblem Eq. ( 6) is often a very difficult task (in Section 3, however, an analytical technique suitable for an important class of continuous systems will be described) and several methods have been proposed to by-pass the problem when closed-form expressions for the eigenfunctions are not known.Suffice is to remember the Rayleigh-Ritz approach where the eigenfunctions are defined by a sum of functions, referred to as admissible functions, satisfying the geometric boundary conditions only [1].
In this section a method giving approximate results is presented, valid when the solution is known for the undamped system only.According to this method, as clearly explained in [4], the solution can be approximated by a finite expansion in terms of the undamped system (known) eigenfunctions ϕ, i.e.

w(x, t) ∼
Substituting this expanded form of the solution in the equation of motion Eq. ( 1) and taking into account the orthogonality relations of the eigenfunctions ϕ with respect to M and K, it is possible to rewrite the state-space Eq. (4) as follows where and the N × N matrices M, C and K are built up by means of the following inner products involving the differential operators M, C, K and the eigenfunctions ϕ with i, j = 1, . . .N. It is worth noticing that both M and K are diagonal.The solution of the related algebraic eigenvalue problem, consisting of a set of 2N eigenvalues (say: s (r) n ) and 2N eigenvectors (say: u n ), allows to uncouple the equations of motion Eq. ( 19) introducing the usual coordinate transformation v = Uη (U denoting the eigenvector matrix), which in the frequency domain yields where a (r) n is the nth element of the diagonalization of A and f n is the nth component of the modal force vector U T f .If a harmonic excitation force of amplitude f 0 acting with angular frequency ω at a coordinate xf is considered, according to Eq. ( 20) f n can be expressed by means of the eigenfunctions ϕ of the undamped system as follows Introducing Eqs (23) in ( 22) and taking into account backwards the links among η, v and r, the expansion Eq. (18) yields the system receptance the superscript ∧ denoting normalization with respect to the square root of a (r) n .

Further investigations for a class of continuous systems
In this section an analytical method for the solution of the differential eigenproblem is presented, valid for a class of vibrating continuous systems with non-proportional damping distributions, according to different damping models [5].The results are then applied to the calculation of the FRFs.Such methodology will in particular be implemented for non-homogeneous Euler-Bernoulli beams in bending vibration.However, it could be easily applied also to strings, shafts, rods and Timoshenko beams with any possible boundary conditions.

Solution of the eigenproblem
In the special case of an Euler-Bernoulli beam in bending vibration, the mass, damping and stiffness operators consist of where m(x) is the mass per unit length of beam, c(x) is the external viscous damping distribution, c in (x) is the internal viscous damping distribution (according to the Kelvin-Voigt model, used in conjunction with the assumption that cross-sectional areas remain planar during deformation) and k(x) = EI(x) is the bending stiffness, or flexural rigidity, in which E is the Young's modulus of the material and I is the area moment of inertia [1].
In order to highlight the effects of non-proportional viscous damping, the differential eigenvalue problem Eq. ( 6) will be solved in the special case in which m(x), c(x) (or c in (x)) and k(x) can be considered piecewise constant on D.
Dividing the beam into P segments of length Dx p = x p − x p − 1 (where x 0 = 0, x P = l, length of the beam), and assuming m(x), c(x) (or c in (x)) and k(x) constant on each segment, the differential eigenvalue problem reduces to a set of P fourth-order ordinary differential equations with constant coefficients of the type (the roman number denoting the derivative order with respect to the spatial coordinate) with appropriate boundary conditions, where which holds for external distributed damping and which holds for internal damping.Note that more complicated damping laws, even involving fractional derivatives, could be easily taken into account simply by modifying the definition of a p as a function of s.In any case, s is obviously the same for every segment.
At this stage it is convenient to convert Eq. ( 26) into a set of four first order equations.
According to the state vector definition the solution for each segment can then be expressed as where Φ p is the pth segment eigenvector matrix, Λ p is the pth segment eigenvalue matrix (with eigenvalues λ 1 = a, λ 2 = −a, λ 3 = ia, λ 4 = −ia) and c p is the pth segment constant vector.Moreover, it is possible to show [5] that the solution at any point x p can be written as where the ith segment eigenvectors matrix and its inverse, written as functions of a i , have the form and B i−1 are 4 × 4 matrices obtained by imposing the continuity of displacement, rotation, moment and shear in x = x i − 1.Clearly, these constraints represent the inner boundary conditions between the adjacent beam segments.Note that B 0 = I and that in the absence of external constraints in x i−1 , B i−1 can be written as A more general expression for B, taking into account external constraints of different kinds, is given in [5].
It is now possible to relate the solution y(l) at one end of the beam to the solution y(0) at the other end, which enables to express the boundary conditions at the ends of the beam in the following form where B e are 2 × 4 matrices depending on the kind of constraints and y 1 (0) = φ 1 c 1 .For example, in case of a clamped end, a pinned end or a free end, they simply are B e = 0 0 1 0 0 0 0 1 champed B e = 0 1 0 0 0 0 0 0 pinned B e = 1 0 0 0 0 1 0 0 free Equation (34) form a linear homogeneous system of four algebraic equations in four unknowns (i.e. the constants c 1 ).Thus the solution of the eigenproblem follows directly by setting to zero the determinant of the coefficient matrix of system Eq.( 34), due to the fact that an algebraic system possesses non-trivial solutions if and only if the determinant of its coefficient matrix is zero, and recalling that the elements of the coefficient matrix of system Eq.(34) depend on the (unknown) eigenvalue s.
It should finally be noted that mathematically the eigenfunctions φ result to be classical solutions (i.e.four times continuously differentiable in D) everywhere, except in a finite subset of D (i.e.x = x p , with p = 1, . . ., P − 1): here they result to be weak (in this case at least one time continuously differentiable) as a consequence of the discontinuities introduced in the functions m(x), c(x) and k(x), which have been assumed piecewise constant on D.

Frequency response functions through modal analysis
If the differential eigenproblem has been solved (i.e. both the eigenvalues sn and eigenfunctions s n are available), the FRFs can be calculated according to Eq. ( 13) after the parameters an have been determined.
To this purpose, it is necessary to write in explicit form the relations among the parameters a n , b n , the differential operators M, C, K and the eigenfunctions φ n , i.e.
which are a direct consequence of Eq. ( 8).Thus, since the eigenfunctions φ n are known, it is possible to calculate a n and b n as well as the modal parameters defined in Section 2.3 simply by applying the definition of inner product.
Introducing the notation of Table 1, according to definition Eq. ( 2) and taking into account the spatial domain partition of Section 3.1, the above inner products can be written in quite similar form as Substituting the eigenfunction expressions given by Eqs (30)-(31) into Eq.(37), straightforward but tedious integrations eventually give where H p = Φ −1 p Π 1 p=1 y 1 (0), Π 1 0 = I, (4 × 4 identity matrix) and E (1) p are 4 × 4 matrices whose elements are respectively in which the eigenvalues λ depend from both the modal index n and the spatial domain partition index p, according to the definitions Eqs ( 27) and (28).

Frequency response functions through direct integration
Besides the modal approach just described, the analytical tools developed in Section 3.1 allow to express the FRFs in a different way, which does not require either eigenvalues or eigenfunctions, i.e. which does not need the solution of the eigenproblem.
This result is achieved by direct integration of the equation of motion, and it is nothing but the sum of the series in Eq. ( 13).
The Euler-Bernoulli beam model with piecewise constant distributions described in Section 3.1 is now considered under the effect of a harmonic excitation force of amplitude f 0 acting with angular frequency ω at a coordinate x f .Since the system is linear-time-invariant, the response will still be a harmonic oscillation at the same angular frequency ω.So, dropping the time dependent terms, the equation of motion for each segment of the beam reduces to where the coefficients a ωp and κ ωp , which are constant within each segment, can be defined according to Eqs ( 27) and (28) by substituting s with iω.Equation ( 40) is a non-homogeneous ordinary differential equation with constant coefficients, since the angular frequency is considered as a given parameter.
As in Section 3.1, in order to find the global solution, the four coefficients of c1 have to be determined by imposing four boundary conditions at the ends of the beam.
By assuming, without loss of generality, that the external force acts at a separation point between two segments (say: x f = x p ), and defining the external force vector in the state-space as follows (κ ωf being κ ω evaluated in x f ) then the system yielding the unknown coefficients of c 1 is simply where it is important pointing out the following remarks: -the matrices B e are the same as in system Eq.(34); -the matrices Π, Φ (and Λ) retain their own definitions as in system Eq.( 34), but the subscript ω means that a np has been substituted by a ωp (i.e. in every definition sn has been changed in iω).
As an example, if a homogeneous beam with two different external damping levels forced in x f (with x 1 x f l) is considered, the system Eq.( 42) simply reduces to so the receptance at a coordinate x (with 0 x x 1 ) can be easily written in function of the four coefficients c 1 The coefficients c 1 are generally rather complicated functions of both x f and ω.However, in some particular cases such functions take a very simple form, as for example in the case of a clamped-free homogeneous Euler-Bernoulli beam forced at its free end (x f = l), whose receptance at a coordinate x is given by where As it is expected, if the angular frequency tends to 0, the system receptance Eq. ( 45) tends to the static deflection of the beam, i.e.
and this limit still holds whatever non-proportional damping distribution is added to the beam model.Note that the described method could also be applied to the calculation of the response to distributed harmonic loads, retaining the notation of system Eq.( 42) and introducing a convolution integral.

Numerical example
A numerical example is presented in order to test and to compare the described methods, and eventually to validate the results by means of a FE model, showing their reliability in problems involving non-proportional damping distributions.
As already shown in [5], the proposed approach is characterised by a high computational efficiency, due to the reduced dimensions of the matrices involved in the numerical procedure.The most crucial point of the modal approach is the zero finding routine needed to solve the algebraic system Eq.(34).This problem has been solved applying the secant method to a real function of complex variable [7].All the codes have been written in Matlab and computed by an AMD-Athlon XP1600+ processor.The zero finding routine runs in less than one second and the finite element model presented in Section 5 runs in some tens of seconds.

Analysis of non-proportional damping effects
The selected numerical example concerns a homogeneous Euler-Bernoulli beam clamped in x = 0 and free in x = l with a non-proportional damping distribution consisting of two different levels of external damping according to Fig. 1 (as regards to the effects of internal damping in similar cases, the reader is referred to [5]).
The parameters for each of the two segments in which the beam is divided are as follows: -length l = 0.30 m, l 1 = 0.10 m, l 2 = 0.20 m; -mass density In this example the distributed external damping density on the second segment (l 1 x l) is kept constant, c 2 = 1.675Ns/m 2 , while on the first segment (0 x l 1 ) it varies from (proportional damping case) to infinity (non-proportional damping limit case).So, different levels of non-proportionality can be obtained by increasing the damping on the first segment only.
In the following, the dimensionless parameter will provide a measure of both the non-proportionality and damping levels.
Figure 2 shows the root loci for the first four modes of the beam under the effect of non-proportional external damping.The curves in the proportional damping case can be obtained by varying both c 1 and c 2 keeping χ = 1: for underdamped modes they are a quarter of a circle.For each mode, the two trajectories (proportional and nonproportional case) start from the same point s(prop) corresponding to c 2 = c 1 = 1.675Ns/m 2 .Due to the particular choice of the damped segment lengths, even for the first and second mode relevant differences can be observed between proportional and non-proportional external damping.The curves are nearly coincident with the proportional c(x) Non-proportional external damping distributio n  case only in the neighbourhood of the starting point s (prop) , then they strongly diverge at higher values of damping and never reach the real axis.
The third mode behaves more like the proportional case and becomes overdamped at high values of χ.On the contrary, the fourth mode curve never reaches the real axis but intersects and then tends to for s (lim) 4 for χ → ∞.The asymptotic behaviour of the root loci of the first, second and fourth mode can be explained considering that as χ → ∞, the clamped-free beam under analysis tends to transform into a clamped-free beam of total length l 2 as shown for a similar example in [5], where the same clamped-free beam with a different damping variation gives completely different root loci.
In order to highlight the effects of a non-proportional damping distribution on the frequency response functions, five equally spaced damping levels from χ 1 = 1 to χ 5 = 801 are considered.
In Fig. 3 the modal indexes of non-proportionality N P for the first four modes are depicted versus the five damping levels χ.
Figure 4 shows a FRF corresponding to a displacement measured at a coordinate x = l1 2 due to a single harmonic force acting at a coordinate x = l 1 + l2 2 with χ = 1 (proportional damping).The receptance modulus | x α xf f (ω)| obtained by the modal approach (Section 3.2) with the first four modes is compared with that obtained by direct integration (Section 3.3) and by FE analysis (Section 5.1) with eight undamped modes.The three curves are in a very good agreement, except for the antiresonances, where the modal truncation error becomes important, and away from the natural frequencies of the first four modes (the only terms included in the modal approach).
The influence of the damping level χ on the receptance is highlighted in Fig. 5, where the measurement and forcing points are the same as in the previous case.The curves obtained with direct integration and FE model are perfectly superimposed, while those obtained with the modal approach and four FE modes (not shown) exhibit a modal truncation error of the same order of magnitude as in Fig. 4. Similar results have been obtained for the phase plots as well.

Finite element numerical validation
In this section the results of the analytical methods presented previously are numerically validated by a finite element model, which has been designed with elements based on an assumed displacement field.

The finite element model
A standard beam element has been chosen, i.e. with transverse displacement wi and rotation θ i as degrees of freedom at each node i and cubic interpolation, so that its mass and stiffness matrices can simply be computed or even found in any textbook [8].
As regards the damping matrix, it is written with the same structure of the mass matrix [9], which corresponds to the case of external distributed viscous damping (see Eq. ( 25)).It would be possible to assume the same structure of the stiffness matrix, which would lead to the case of internal damping.
The damping matrix can be non-proportional and the expression of the FRF can be found using the state space approach, i.e. solving a complex eigenproblem.This technique is not commonly implemented in FE procedures because it doubles the dimensions of the matrices, thus significantly increasing the computational effort.
However, it is possible to by-pass this limitation by solving two smaller eigenproblems, i.e. by following a procedure which is very close to that described in Section 3.4.
Two important concepts have to be underlined at this point: -the expansion of the solution is now written in terms of eigenvectors (namely ϕ j , not to be confused with the eigenfunctions ϕ(x)) of the undamped discrete system; -to cut down the computational effort, the order of the system can be reduced by taking into account a subset of only N eigenvectors ϕ j with N M , M being the number of degrees of freedom of the FE model.It should be stressed that the selected sequence of eigenvectors does not necessarily include the first N or even a set of N close to ϕ j , albeit this has been che choice for the numerical examples herewith presented.
Under the hypothesis of a single force of amplitude f 0 acting with angular frequency ω on the physical d.o.f.m, it is then possible to demonstrate that the receptance at a d.o.f.h is The frequency response functions can therefore be expressed as functions of a subset of real eigenvectors ϕ j of the undamped system, and of the complex eigenvalues and eigenvectors s (r) n , u n of the (low order) damped system.

Numerical results
The numerical results achieved with the analytical methods have been compared with a model consisting of 132 elements.As expected, the modal frequencies of this FE model without damping are in a very good accordance with those of a classical Euler-Bernoulli beam, as confirmed also by convergence tests.Note also that the parameters of the system are those of a proper beam with a very high length to thickness ratio.
The root loci represented in Fig. 2 show an almost perfect coincidence between the analytical and numerical results, and also the comparison of the FRFs both in modulus (Figs 4 and 5) and phase is completely satisfying.
In the selected frequency band and taking into account the first eight modes of the undamped system, the FE model receptances are exactly superimposed on the curves obtained through the exact theoretical approach, whilst the absence of higher modes becomes significant in the upper part of the frequency domain (not represented in the figures).Finally, it should be noted that the analytical dotted line in Fig. 4 (modal approach with four modes) is also representative of the effects of using four modes in the FE method.

Conclusions
In this paper two general methods have been proposed to compute the exact frequency response functions of continuous systems with non-proportional damping distributions, focusing the attention on the Euler-Bernoulli beam model.
The first method is based on the modal approach and takes advantage from the orthogonality properties of the eigenfunctions, which have been demonstrated for vibrating continuous systems whose equations of motion are characterized by self-adjoint differential operators.
On the contrary, the second method exploits a direct integration of the equations of motion, thus being not affected by any modal truncation error at high frequencies.
Both methods, starting from a partition of a continuous system in homogeneous substructures, have been developed combining the reduction of the differential equation order with a transfer matrix technique.
As a result, they have shown a high computational efficiency, due in part to the invariance of the dimensions of the matrices involved in the numerical procedures with respect to the number of substructures in which the system has been divided.
The presented methods have then been applied to test the accuracy of a technique based on the approximation of the solution by a finite expansion in terms of the undamped system eigenfunctions, showing its reliability.
Finally, the numerical results have been successfully validated by means of a finite element procedure, in which the computational effort due to the non-proportional damping distributions has been significantly reduced applying again the same technique of approximation based on a selected set of undamped eigenvectors.
The described analytical tools enable a complete frequency domain study of the effects of generalized damping distributions on continuous systems.However, the fundamentals to extend the analysis to the time domain are included as well.
In particular, the introduction of a new modal index of non-proportionality has been proposed, following a discussion about the meaning of the modal parameters in case of non-proportional damping.
Possible applications of these methods could regard the analysis and passive control of vibrating elements consisting of non-homogeneous bars, shafts beams, or more complicated systems, such as for example ducts or pipe-lines, in which the proportional damping assumption could not be valid to describe the dynamics with sufficient accuracy.
Future work will extend the proposed approach to different vibrating continuous systems including more complicated damping laws, even involving fractional derivatives, and possibly the effects of random or/and moving loads.

Table 1
Short notation for the modal parameters