Nondimensional system parameter used for Fig. 8.

## 1. Introduction

The vibrations of elastic structures such as strings, beams, and plates can be described in terms of waves traveling in waveguides (Cremer et al., 1973; Graff, 1975; Fahy, 1987). While the subject of wave propagation has been extensively studied in the fields of acoustics in fluids and solids rather than vibrations of elastic structures, wave analysis techniques have been employed to reveal physical characteristics associated with structural vibrations of elastic media (Argento & Scott, 1995; Kang & Tan, 1998). One of the advantages of the wave analysis technique, when applied to the structural vibration analysis, is its compact and systematic approach to analyze complex structures with discontinuities (Mace, 1984; Yong & Lin, 1989; Kang et al., 2003; Mei & Mace, 2005). Applying the concept of wave reflection and transmission, Mace (1984) obtained the frequency equations of Euler-Bernoulli beams including waves of both propagating and near-field types. By the phase-closure principle, also referred to as the wave-train closure principle (Cremer et al., 1973), Mead (1994) determined natural frequencies of Euler-Bernoulli beams. This principle states that if the phase difference between incident and reflected waves is an integer multiple of *2*, then the waves propagate at a natural frequency and their motions constitute a vibration mode. Based on the same principle, Kang (2007) presented a systematic approach to the free and forced vibration analysis of multi-span beams.

The classical method, known as the normal mode or eigenfunction expansion, of solving the forced vibration problem of a distributed parameter system involves expansion of the forcing function into the eigenfunctions of the associated free vibration problem. While this method is theoretically sound and powerful, the method is difficult to implement when the problem to be solved is a non-self-adjoint system typically due to complicating effects such as damping, discontinuities, or non-classical boundary conditions, for which case obtaining the exact eigensolutions is not often feasible. Although approximate eigensolutions may be used instead of exact ones, the problem still persists in the form of poorly convergent solution and/or significant error in the solution. As an alternative approach to solve forced vibration problems, Yang and Tan (1992) presented a method for evaluating exact closed-form transfer functions for a class of one-dimensional distributed parameter systems. Applying the energy functionals of constrained and combined damped systems, Yang (1996 a, 1996b) presented a method to obtain a closed-form transient response solution in eigenfunction series for a distributed damped system.

The dynamic displacement of any point in an elastic waveguide can be determined by superimposing the amplitudes of the constituent waves traveling along the waveguide, which is a basis of wave propagation. Based on this simple fact, an alternative technique for the free and forced vibration analysis of one-dimensional distributed parameter systems is presented. The method of normal mode expansion is often difficult to implement for nonself-adjoint systems with complicating effects such as mode couplings, non-proportional damping, discontinuities, or arbitrary boundary conditions, since the method requires eigensolutions or assumed normal modes as a priori. However, this alternative analysis technique based on the elastic wave propagation does not pose such a requirement and leads to the exact, closed-form, distributed transfer function of a distributed parameter system. The general wave solution of the equation of motion governing the dynamics of a waveguide is cast into a matrix form in terms of the constituent waves defined in the Laplace domain. The spatial amplitude variation of the traveling wave is represented by the field transfer matrix and the amplitude distortion of the traveling wave incident upon a discontinuity due to geometric or kinetic constraints is described by the local wave reflection and transmission matrices. Combining these matrices in a progressive manner along the waveguide by applying the concepts of global wave reflection and transmission matrices leads to the exact characteristic equation and corresponding mode shapes for the free response analysis and the transfer function of the system for the forced response analysis. The transient response solution for a complex system can be obtained through the Laplace inversion of the transfer function using numerical inversion algorithms. The exact frequency response solution, which includes infinite normal modes of the system, can be obtained in terms of the complex frequency response function from the transfer function. One of the main advantages of this analysis technique is its systematic formulation resulting in a recursive computational algorithm which can be implemented into highly efficient computer codes. This systematic approach also allows modular formulation which can be readily expandable to include additional discontinuities with little alteration to the existing formulation. In addition, it is also computationally advantageous that the technique always results in operating matrices of a fixed size regardless of the number of discontinuities in a waveguide. This analysis technique is applicable to any one-dimensional waveguides (strings, axial rods, torsional bars, beams, and frame structures), in particular systems with multiple point discontinuities such as viscoelastic supports, attached inertias, and geometric/material property changes. The analysis technique is demonstrated using the second order wave equation, fourth order beam equation, and sixth order curved beam equation.

## 2. Second order systems

The free transverse vibration of a taut string, longitudinal vibration of a thin bar, and the torsional vibration of a shaft are governed by the equation of motion in the same form, the wave equation, and thus they are mathematically analogous. Therefore, with no loss of generality, the transverse vibration of the string is taken as a representative problem for the description of the present analysis technique based on wave propagation. The equation governing the transverse motion of a uniformly damped string of span length L is

where W is the transverse displacement, X the spatial variable, t the temporal variable; and T denotes the constant tension, C_{e} the damping coefficient, and m the mass per unit length of the string. With introduction of the following non-dimensional variables and parameters

the equation of motion takes the non-dimensional form of

where the prime (‘) and dot () denote the differentiation with respect to x and , respectively. Applying the Laplace transform to Eq. (3) yields

where s denotes the Laplace variable and zero initial conditions are assumed.

### 2.1. Wave solution

Denoting C as the amplitude of the wave traveling along the string, the solution of Eq. (4) can be assumed in wave form

where is the non-dimensional wavenumber normalized against span length L. Applying the above wave solution to Eq. (4) gives the frequency equation of the problem

from which the general wave solution can be found as the sum of two constituent waves

where the coefficient C represents the amplitude of each wave component with its traveling direction indicated by the plus (+) or () sign. Note that is complex valued for nonzero c_{e}, hence the classification of the wave into propagating and attenuating waves does not apply to this case. Defining

where

### 2.2. Wave reflection and transmission

When a wave traveling along a string is incident upon a discontinuity such as an elastic support, geometric/material property change, or boundary, it is reflected and transmitted at different rates depending on the properties of the discontinuity. The rates of wave reflection and transmission can be determined in terms of wave reflection and transmission coefficients. For example, consider an infinitely long string constrained at a local coordinate _{c}), a transverse spring (K_{c}), and a viscous damper (C_{c}).

When a positive-traveling wave

Since

In addition, the kinetic equilibrium condition at =0 states that

where represents the kinetic properties of the constraint as

Note that k_{c}, c_{c}, and m_{c} in Eq. (13) are the non-dimensional spring constant, damping coefficient, and the attached mass, respectively, defined by

Equation (12) leads to

Combining Eq. (11) and Eq. (15), the wave reflection and transmission coefficients can be found as

When a wave is incident upon a boundary at =0, it is only reflected, therefore

which leads to

In the limiting case where =, it can be seen that r=1 which is the wave reflection coefficient for the classical fixed boundary. When the wave is incident upon a series of discontinuities along its traveling path, it is more computationally efficient to employ the concepts of global wave reflection and transmission coefficients, in particular when the free or forced vibration analysis of a multi-span string is sought. These coefficients relate the amplitudes of incoming and outgoing waves at a discontinuity. Consider wave motion in a multi-span string as illustrated in Fig. 2. Define R_{ir} as the global wave reflection coefficient which relates the amplitudes of negative- and positive-traveling waves on the right side of discontinuity i such that

Since_{ir} in terms of the global wave reflection coefficient on the left side of discontinuity i+1; i.e.,

In addition, by combining the following wave equations at discontinuity i

the relationship between the global wave reflection coefficients on the left and right sides of discontinuity i can be found as

R_{ir} and R_{il} progressively expand to include all the global wave reflection coefficients of discontinuities along the string before terminating its expansion at the boundaries where

While the global wave reflection coefficient relates the amplitudes of waves traveling in the opposite direction of each other within a subspan i, there is a need for another coefficient which relates the amplitudes of waves traveling in the same direction in two adjacent subspans. This inter-span coefficient is particularly useful when the mode shape or forced response of a string with several or more subspans needs to be determined since it allows an intuitive and systematic formulation of the system’s transfer function. Denoting this inter-span wave transfer coefficient as the global wave transmission coefficient T_{i}, define

Rewriting Eq. (21) by applying

The global wave reflection and transmission coefficients are the key elements in determining the exact transfer function of a multi-span string as discussed in the following sections.

### 2.3. Free response analysis

The global reflection and transmission coefficients of waves traveling along a multi-span sting are now combined with the field transfer function to analyze the free response of a multi-span sting. With reference to Fig. 2, consider wave motion in the first span. Based on the definition of the global wave reflection coefficient, at the boundary

However, recalling

For nontrivial solutions,

which is the characteristic equation in terms of the Laplace variable s for the multi-span string with arbitrary discontinuities and boundaries. This remarkably simple expression for the characteristic equation is due to the fact that R_{1r} recursively expands to include all the effects of constraints in the remaining side of the string until its expansion terminates at the rightmost boundary which yields R_{nl}r_{n}. This equation simply states that when the string system vibrates at one of its natural frequencies, r_{1}R_{1r}1. As a simple example, for a single span uniformly damped string fixed at both ends, r_{1}r_{2}1 from Eq. (18) and R_{1r}e^{2i} from Eq. (20). Therefore, F(s)e^{2i}10 gives the natural wavenumbers _{n}n (n1,2,3,...), or in terms of the non-dimensional frequency

The mode shapes of the multi-span string system can be found in a systematic manner by relating wave amplitudes between two adjacent subspans. Let _{i} denote the local coordinate within span i. The transverse displacement at any point in span i can be expressed as

However, due to Eq. (19)

Now by applying the global wave transmission coefficient defined in Eq. (25) to

Assume a disturbance arise in span 1; i.e., a wave originates and starts traveling from the leftmost boundary of the string. Then, by successively applying the global transmission coefficient of each discontinuity on the way up to the first span, the mode shape of span i can be found in terms of wave amplitude

For example, consider a fixed-fixed undamped string with three supports specified by _{2}=5+0.1s^{2}, _{3}=7+0.1s^{2}, and _{4}=4+0.1s^{2} according to Eq. (13). l_{1}0.25, l_{2}0.3, l_{3}0.25, and l_{4}0.2 are assumed. Once the global wave reflection coefficient at each discontinuity has been determined, one can apply Eq. (29) to find the natural frequencies. Shown in Fig. 3 is the plot of the characteristic equation, where the first three natural frequencies are indicated. The mode shapes can be found from Eq. (33) in a systematic way once the global wave transmission coefficient at each discontinuity has been determined. Figure 4 shows the mode shapes for the first three modes.

### 2.4. Transfer function analysis

Consider a multi-span string subjected to an external point load_{0} as shown in Fig. 5. Since the waves injected by the load travel in both direction from the point of loading, a set of local coordinates {, ^{*}} is introduced such that the wave traveling toward each boundary of the string is considered positive as indicated in Fig. 5. Let _{0} and xx_{0} within span 1, respectively. The transverse displacement of span 1 of the string can be expressed as

Defining v as the wave injected by the applied load, the difference in amplitudes between the incoming and outgoing waves at the loading point is

where v can be determined from the geometric and kinetic continuity conditions at =0 as

Applying the global wave reflection coefficient on each side of the loading point gives

where the asterisk (^{*}) signifies that the global wave reflection coefficient is defined in the region xx_{0} to distinguish it from the one defined in the region x<x_{0}. Combining Eq. (35) and Eq. (37), one can determine the ampltitude of the wave that rises at the loading point in each direction

Note that _{0} and xx_{0}, respectively. It is evident that and _{0} can be found; i.e., in the region x<x_{0}

Now, in the same manner as for the mode shape analysis, since _{0}

Note that the ratio

is the transfer function governing the forced response of any point in span i due to the point loading at xx_{0}. The Laplace inversion of G_{i}(_{i}, x_{0}; s) is the Green’s function of the problem. Thefore, denoting _{0}

The exact Laplace inversion of G_{i}(_{i}, x_{0}; s) in close form is not feasible in general, in particular for multi-span string systems. One may have to resort to the numerical inversion of Laplace transforms. It is found that the algorithm known as the fixed Talbot method (Abate & Valko, 2004), which is based on the contour of the Bromwich inversion integral, appears to perform the numerical Laplace inversion in Eq. (42) with a satisfactory accuracy and reasonable computation time. In this method, the accuracy of the results depends only on the number of precision decimal digits (denoted with M in the algorithm) carried out during the inversion. It is found that for well-damped second- and higher order distributed parameter systems, M32 and for lightly damped or undamped systems, M64 gives acceptable results. To demonstrate the effectiveness of the present analysis approach, consider the unit impulse response of a single span undamped string, in particular the response near the point of loading immediately after the loading, which is well known for its deficiency in numerical convergence. The response solution by the method of normal mode expansion is

The corresponding transfer function from Eq. (40) is

Shown in Fig. 6 is the comparison of the response solutions given in Eq. (43) and Eq. (44) when 0.001 near xx_{0}0.5. N2,500 and N20,000 are used for the evaluation of the series solution, while M32 and M64 are used for the numerical Laplace inversion of the transfer function in Eq. (44). It can be seen that the series solution in Eq. (43) fails to accurately represent the actual impulse response behavior with N2,500. This is expected for the series solution since it would take a large number of harmonic terms (N10,000 for this example) to represent such a sharp spike due to the impulse. It can be seen that the result with M32 reasonably represents the actual behavior, and the result with M64 is almost comparable to the series solution with N20,000. However, if one tries to obtain the response at a time very close to the moment of impact, the numerical Laplace inversion becomes extremely strenuous or beyond the machine precision of the computing machine. This is because the expected response would consist of waves that have unrealistically short wavelengths. This is not a unique problem for the present wave approach since the same problem would manifest itself in the series solution given in Eq. (43), requiring an impractically large number of harmonics terms for a convergent solution.

If

Therefore the frequency response at any point within any subspan of the string can be obtained by; e.g., for span i in the region x<x_{0}

One of the main advantages of this approach is its systematic formulation resulting in a recursive computational algorithm which can be implemented into highly efficient computer codes consuming less computer resources. This systematic approach also allows modular formulation which can be easily expandable to include additional subspans with very minor alteration to the existing formulation. Another significant advantage of the present wave-based approach to the forced response analysis of a multi-span string is that the eigensolutions of the system is not required as a priori as in the method of normal mode expansion which assumes the forced response solution in terms of an infinite series of the system eigenfunctions – truncated later for numerical computations. However, exact eigensolutions are often difficult to obtain particularly for non-self-adjoint systems, and also approximated eigensolutions can result in large error. In contrast, the current analysis technique renders closed-form transfer functions from which exact frequency response solutions can be obtained.

## 3. Fourth order systems

The analysis techniques described by using the vibration of a string can be applied to the transverse vibration of a beam of which equation of motion is typically a fourth order partial differential equation. Denoting X and t as the spatial and temporal variables, respectively, the equation of motion governing the transverse displacement W(X,t) of a damped uniform Euler-Bernoulli beam of length L subjected to an external load P(X,t) is

where m denotes the mass per unit length, EI the flexural rigidity, C_{e} the external damping coefficient of the beam. With introduction of the following non-dimensional variables and parameters

the equation of motion takes the non-dimensional form of

Applying the Laplace transform to Eq. (48) yields

Letting

where is the non-dimensional wavenumber normalized against span length L. Applying the above wave solution to Eq. (49) gives the frequency equation of the problem

from which the general wave solution can be found as the sum of four constituent waves

where the coefficient C represents the amplitude of each wave with its traveling direction indicated by the superscript; plus (+) and minus (–) signs denote the wave’s traveling directions with respect to the x-coordinate. The subscripts a and b signify the spatial wave motion of the same type traveling in the opposite direction. Note that is complex valued in general. The general wave solution in Eq. (52) may be re-expressed in vector form by grouping the wave components (wave packet) traveling in the same direction; i.e.,

and then

where f(x;s) is the diagonal field transfer matrix (Mace, 1984) given by

which relates the wave amplitudes by

### 3.1. Wave reflection and transmission matrices

For a wave packet with multiple wave components, the rates of wave reflection and transmission at a point discontinuity can be found in terms of the wave reflection matrix r and wave transmission matrix t, in the same manner described in Section 2.2. When the flexural wave packet in Eq. (52) travels along a beam and is incident upon a kinetic constraint (=0) which consists of, for example, transverse (K_{t}) and rotational (K_{r}) springs and transverse damper (C_{t}), r and t at the discontinuity can be found by applying the geometric continuity kinetic equilibrium conditions at =0; i.e., with reference to Fig. 7, one can establish the following matrix equations

where

However, as previously discussed in Section 2.2, when the wave packet is incident upon a series of discontinuities along its traveling path, it is more computationally efficient to employ the concepts of global wave reflection and transmission matrices. These matrices relate the amplitudes of incoming and outgoing waves at a point discontinuity. When compared to the string problem, the only difference in formulating these matrices for the beam problem is to use vectors and matrices instead of single coefficients. Therefore, with reference to Fig. 2, let R_{ir} as the global wave reflection matrix which relates the amplitudes of negative- and positive-traveling waves on the right side of discontinuity i such that

Since_{ir} in terms of the global wave reflection matrix on the left side of discontinuity i+1; i.e.,

In addition, by combining the following wave equations at discontinuity i

the relationship between the global wave reflection matrices on the left and right sides of discontinuity i can be found as

R_{ir} and R_{il} progressively expand to include all the global wave reflection matrices of intermediate discontinuities along the beam before terminating its expansion at the boundaries. Since incident waves are only reflected at the boundaries, one can find the following wave equations

where r can be found by imposing the force and moment equilibrium conditions at the boundary; e.g., for the same kinetic constraint previously considered

Now, to determine the global wave transmission matrix T_{i}, define

Rewriting Eq. (60) by applying

where I_{22} is the 22 identity matrix.

### 3.2. Free response analysis

The global reflection and transmission matrices of waves traveling along a multi-span beam are now combined with the field transfer matrix to analyze the free vibration of the beam. With reference to Fig. 2, where

However, due to Eq. (57), it can be found that

Applying the condition for non-trivial solutions to the above matrix equation, one can obtain the following characteristic equation in terms of the Laplace variable s

By applying a standard root search technique (e.g., Newton-Raphson method or secant method) to Eq. (71), one can obtain the natural frequencies of the multi-span beam.

The mode shapes of the multi-span beam can be systematically found by relating wave amplitudes between two adjacent subspans, in the same way described in Section 2.3. Defining _{i} as the local coordinate in span i, the transverse displacement of any point in span i can be found as

However, due to Eq. (55)

Since

Assume a wave packet originates and starts traveling from the leftmost boundary of the beam. By successively applying the global transmission matrix of each discontinuity on the way up to the first span, the mode shape of span i can be found in terms of wave amplitude

Discontinuity | ||||||

Constraint | 1 | 2 | 3 | 4 | 5 | 6 |

Location | 0 | 0.15 | 0.35 | 0.6 | 0.82 | 1 |

kt | 3,000 | 2,500 | 1,500 | 3,500 | ||

kr | 0 | 250 | 150 | 100 | 300 | |

ct | 0 | 0 | 0 | 0 | 0 | 0 |

where note that the amplitude ratio between the two wave components _{1}=10.294, _{2}=12.038, and _{3}=14.148, from which the damped natural frequencies of the beam can be determined by using Eq. (52). It should be noted from the computational point of view that the present wave approach always results in operationg matrices of a fixed size regardless of the number of subspans. However if the classical method of separation of variables is applied to solve a multi-span beam problem, the size of matrix that determines the eigensolutions of the problem increases as the number of subspans increases, which may cause strenuous computations associated with large-order matrices.

### 3.3. Transfer function analysis

Consider a multi-span beam subjected to an external force applied at x=x_{0} where x_{0} is located in subspan m. Let

Now, in order to determine the actual wave amplitudes_{0} travel in both directions, a new set of local coordinates {, ^{*}} is defined such that the waves traveling towards each end of the beam are measured positive as indicated in the Fig. 5. Let _{0} and xx_{0}, respectively. The discontinuity in the shear force at x=x_{0} can be expressed in term of the difference in amplitudes between the incoming and outgoing waves at the discontinuity such that

where v is the wave vector injected by the applied force which can be determined by the geometric and kinetic continuity conditions at =0 as

Applying the global wave reflection matrices on each side of discontinuity 1 gives

Combining Eq. (62) and Eq. (63), one can find the wave amplitudes that give rise at the point of loading

_{0} and xx_{0}, respectively. It is evident that _{0}. Applying the results in Eq. (65) and Eq. (66) to Eq. (70), the transverse displacement of any point in span 1 on either side of the point of loading can be found. For example, for span 1 in the region x<x_{0}

In the same manner as for free response analysis, the forced response solution in Eq. (67) can be generalized for a multi-span beam by applying the global wave transmission matrix. For x<x_{0}, _{0} due to external loading at x=x_{0} can be determined from

One can find the same expression in the region xx_{0}, but in terms of _{0}. The Laplace inversion of _{0},

The exact Laplace inversion of G_{i}(, x_{0}; s) in close form is prohibitively difficult for the beam problems considered in this study, in particular multi-span beams with non-classical intermediate supports and boundary conditions. One may have to resort to the numerical inversion of Laplace transforms. Shown in Fig. 9 is the unit impuse response curve of a non-uniformly damped three-span beam, whose system parameters are listed in Table 2, sampled at x=0.15 and x=0.825 for x_{0}= 0.4. The unit impulse response curve obtained from the method of normal mode expansion is also shown for comparison. For the results obtained from the method of normal mode expansion, the mode shape function of sinnx (n=1,..., 10) are used for the expansion and the constraint forces at the intermediate supports are represented by using the Dirac-Delta function. The resulting set of modal equations are numerically integrated using a fixed-step Runge-Kutta algorithm since the beam is nonproportionally damped hence the modal equations cannot be decoupled. It can be seen that both results show an excellent agreement.

Discontinuity | |||||

Constraint | 3 | 2 | 1=1^{*} | 2 | 3 |

Location | 0 | 0.15 | 0.4 | 0.65 | 1 |

kt | 2,000 | 0 | 3,000 | ||

kr | 0 | 0 | 0 | 0 | 0 |

ct | 0 | 15 | 0 | 10 | 0 |

Regarding the numerical Laplace inversion to obtain the transient response of a Euler-Bernoulli beam model due to an impulse, it should be noted that there exists a singular behavior of the response for small values of time; i.e., the initial response. This is not due to the algorithm of the numerical Laplace inversion, but due to the inherent deficiency of the classical Euler-Bernoulli beam theory that neglects the effects of rotary inertia and shear

deformations. According to this simple beam theory, a wave of very high frequency (or very short wavelength) would be predicted instantaneously at remote locations, which is physically unacceptable since this can be interpreted as that such a wave travels at an infinity velocity. This abnormal behavior of the Euler-Bernoulli beam can be readily verified by observing the phase velocity c of the dispersion equation (Eq. (59) with s=i and =c) continues to increase with increasing . Therefore, for a meaningful initial transient response to an impulsive load, one must employ Rayleigh or Timoshenko beam models. Studies on the wave reflection and transmission behavior at various point discontinuities and boundaries using a Timoshenko beam model can be found in Refs. (Tan & Kang, 1999; Mei and Mace, 2005).

If

Discontinuity | ||||||||

Constraint | 4 | 3 | 2 | 1=1^{*} | 2^{*} | 3^{*} | 4^{*} | 5^{*} |

Location | 0 | 0.12 | 0.26 | 0.39 | 0.56 | 0.69 | 0.84 | 1 |

kt | 5,000 | 4,000 | 3,500 | 0 | 3,000 | 5,000 | 0 | |

kr | 100 | 0 | 100 | 0 | 100 | 50 | 0 | 0 |

ct | 0.5 | 0.1 | 0 | 0 | 0 | 0.1 | 0.1 | 0.2 |

Therefore the frequency response at any point within any subspan of the beam can be obtained by; e.g., for span i in the region x<x_{0}

Shown in Fig.10 is the frequency response sampled at x=0.325 and x=0.92 for the non-uniformly damped six-span beam, whose system parameters are listed in Table 3, under harmonic excitation at x=0.39.

## 4. Sixth order systems

The analysis techniques described in the previous sections can be applied to the vibrations of higher order, one-dimensional, distributed parameter systems. For example, the in-plane vibration of a planar curved beam is governed by a sixth order partial differential equation. Consider a thin planar curved bam. Neglecting the effects of rotary inertia and shear deformations, the coupled equations of motion governing the flexural displacement W and the tangential displacement U of the centroidal axis are

where E denotes the Young’s modulus, I the second moment of inertia of the cross-section about the centroid, the angular coordinate, R the constant radius of curvature for the given range of angle , A the cross-sectional area, the mass density, and t the time variable. Details of deriving these equations of motion can be found in Ref. (Graff, 1975). With the following non-dimensional variables and parameters:

where t_{0} is a characteristic time constant and k is the curvature parameter, Eq. (79) takes the following non-dimensional form in the Laplace domain

where s represents the Laplace variable. In order to determine the general wave solutions for the equations of motion in Eq. (81), the following form of harmonic wave solutions are first assumed

where denotes the wavenumber, normalized against R, of the wave traveling along the centroidal axis. Substituting the above wave solutions into Eq. (81) leads to

Since the determinant of the matrix in Eq. (85) must vanish for nontrivial solutions, one can obtain following frequency equation

The six roots (essentially three complex conjugates, namely _{n}(s), n=1,2, 3) of Eq. (86) can be obtained in a closed form by transforming Eq. (86) into a cubic equation. Each conjugate pair is the wavenumber of the constituent wave components corresponding to each of the two wave modes of the curved beam. Therefore, it can be found that

where the coefficient C represents the amplitude of each flexural wave component with its traveling direction indicated by the / signs. The subscript n signifies the spatial wave motion of the same type traveling in the opposite direction. _{n} is the amplitude ratio between the flexural and tangential waves which can be found from Eq. (85) that

The above wave solutions may be recast in vector form by grouping the waves traveling in the same direction in a wave packet; i.e.,

where

One now can see that the analysis techniques described in the previous sections for the multi-span beam can be applied to a multi-span curved beam in the same manner to determine the local wave reflection and transmission matrices, global wave reflection and transmission matrices, and the transfer function (Kang et al., 2003). The only difference is the size of matrix, three by three for the curved beam. It should be also noted that the present analysis techniques are applicable to one-dimensional distributed parameter systems involving other types of discontinuities such as geometric/material changes and cracks, as long as the properties of those discontinuities can be characterized by wave reflection and transmission coefficients or matrices. In what follows, an example of free vibration analysis of a three-span curved beam with curvature changes is presented.

### 4.1. Wave reflection and transmission at a curvature change

Consider a curved beam with two subspans of different curvatures joined at =0. Assuming that the wavenumbers of the waves traveling in each subspan are

and the kinetic continuity conditions give

The local wave reflection and transmission matrices can be found by solving the above equations. The curved beam system shown in Fig. 11 has three subspans of equal span angle of 60 but different radius of curvature. Both ends are clamped. The natural frequencies of the system can be obtained from

where r_{1} is the wave reflection matrix at the clamed boundary, which is

and R_{1r} is the global wave reflection matrix that includes the effects of the curvature changes and the other clamped boundary. The natural frequencies are compared in Table 4 with the results obtained by the finite element method (FEM) (Maurizi et al., 1993; Krishnan and Suresh, 1998) and the cell discretization method (CDM) (Maurizi et al., 1993). Note that the natural frequencies obtained from the present wave analysis are exact since both propagating and attenuating wave components are considered in the formulation.

Present | FEM | CDM | |||

12 curvedelements | 24 straight elements | 100 Degreesof freedom | |||

Mode | Ref. 1 | Ref. 2 | Ref. 1 | Ref. 1 | |

1 | 2.683 | 2.680 | 2.701 | 2.685 | 2.671 |

2 | 4.834 | 4.824 | 4.828 | 4.828 | 4.780 |

3 | 9.565 | 9.536 | 9.543 | 9.543 | 9.413 |

4 | 14.585 | 14.527 | 14.535 | 14.535 | 14.309 |

5 | 21.865 | 21.749 | 21.751 | 21.751 | 21.265 |

## 5. Summary

An alternative approach to the dynamic response analysis of one-dimensional distributed parameter systems by applying the concepts of wave motions in elastic waveguides is presented. The analysis techniques are demonstrated using the vibration of multi-span strings, straight beams, and curved beams with general support and boundary conditions, however other one-dimensional systems such as rods and higher order beam models can be treated in the same manner. The present approach allows a systematic formulation that leads to exact, closed-form, distributed transfer functions from which the transient response and frequency response solutions can be obtained. In addition, the present analysis approach results in recursive computational algorithms that always involve computations of fixed-size matrices regardless of the number of subspans, which can be implemented into highly efficient computer codes. Since neither exact nor approximate eigensolutions are required as a priori, the present wave-based approach is suitable for the forced response analysis of non-self-adjoint systems. There are two limiting cases that may affact the analysis accuracy and numerical efficiency of the present wave approach: (1) when the waveguide contains a very small amount of inertia or flexibilty (such as massless elements or rigid bodies), which results in making the wavelenght of the constituent waves larger than the span length of the waveguide; and (2) when the waveguide is extremely flexible or its response contains a very sharp impulsive spike (such as the impulse response example shown in Section 2.4), which result in making the wavelength of the constituent waves unrealistically small. In paractice, however, most engineering structures have resonable amounts of inertia, flxibility, and also damping such that these two limiting cases will be rarely encountered.