Exact Transfer Function Analysis of Distributed Parameter Systems by Wave Propagation Techniques

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

$$T\frac{{\partial }^{2}W}{\partial X^2}+C_e\frac{\partial W}{\partial t}=m\frac{{\partial }^{2}W}{\partial t^2}$$E1

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

$$w=W/L x=X/L \tau =t/t_0 c_e=C_{e}L/t_0 t_0=L/c_0 c_0=\sqrt{T/m}$$E2

the equation of motion takes the non-dimensional form of

$$w^{″}+c_e\dot{w}=\ddot{w} (0\le x\le 1)$$E3

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

$$\overline{w^{″}}(x;s)+c_{e}s\overline{w}(x;s)=s^2\overline{w}(x;s)$$E4

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

$$\overline{w}(x;s)=C{e}^{i\gamma x}$$E5

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

$${\gamma }^2+s^2-c_{e}s=0\text{or}{\gamma }^2=c_{e}s-s^2$$E6

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

$$\overline{w}(x;s)=C^{+}{e}^{-i\gamma x}+C^{-}{e}^{i\gamma x}$$E7

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 $f(x;s)=e^{-i\gamma x}$ as the field transfer function which relates the wave amplitudes by

$$C^{+}(x+x_0)=f{C}^{+}(x_0) \text{or} C^{-}(x+x_0)=f^{-1}{C}^{-}(x_0)$$E8

where$f^{-1}(x;s)=e^{i\gamma x}$, the wave solution in Eq. (7) can be re-written as

$$\overline{w}(x;s)=f{C}^{+}+f^{-1}{C}^{-}$$E9

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 $\xi =0$ as shown in Fig. 1, where the constraint is a point support consisting of an attached mass (M_c), a transverse spring (K_c), and a viscous damper (C_c).

Figure 1.

Wave reflection and transmission at a discontinuity.

When a positive-traveling wave $C^{+}$ is incident upon the support, it gives rise to a reflected wave $r{C}^{+}$ and a transmitted wave$t{C}^{+}$, where r and t represent the wave reflection and transmission coefficients, respectively. The resulting transverse displacements at the left and right of the constraint are

$${\overline{w}}_l=C^{+}+r{C}^{+}{\overline{w}}_r=t{C}^{+}$$E10

Since ${\overline{w}}_l={\overline{w}}_r$ at =0, one can find

$$1+r=t$$E11

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

$${\overline{w^{\prime }}}_r-{\overline{w^{\prime }}}_l=\phi {\overline{w}}_r$$E12

where represents the kinetic properties of the constraint as

$$\phi (s)=k_c+c_{c}s+m_c{s}^2$$E13

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

$$k_c=K_{c}L/T{c}_c=C_c{c}_0/T{m}_c=M_c{c}_0^2/TL$$E14

Equation (12) leads to

$$\gamma (r+t-1)=i\phi t$$E15

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

$$r=\frac{i\phi }{2\gamma -i\phi }andt=\frac{2\gamma }{2\gamma -i\phi }$$E16

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

$${\overline{w^{\prime }}}_l=-\phi {\overline{w}}_l$$E17

which leads to

$$r=\frac{\gamma +i\phi }{\gamma -i\phi }$$E18

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

$$C_{ir}^{-}=R_{ir}{C}_{ir}^{+}$$E19

Figure 2.

Waves traveling along a multi-span string.

Since$C_{ir}^{-}=f_i{R}_{(i+1)l}{C}_{ir}^{+}$, one can find R_ir in terms of the global wave reflection coefficient on the left side of discontinuity i+1; i.e.,

$$R_{ir}=f_i^2{R}_{(i+1)l}$$E20

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

$$C_{ir}^{+}=t_i{C}_{il}^{+}+r_i{C}_{ir}^{-}and{C}_{il}^{-}=t_i{C}_{ir}^{-}+r_i{C}_{il}^{+}$$E21

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

$$R_{il}=r_i+t_i^2{(R_{ir}^{-1}-r_i)}^{-1}$$E22

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

$$C_1^{+}=r_1{C}_1^{-}{C}_n^{-}=r_n{C}_n^{+}$$E23

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

$$C_{ir}^{+}=T_i{C}_{(i-1)r}^{+}$$E24

Rewriting Eq. (21) by applying $C_{il}^{+}=f_{(i-1)}{C}_{(i-1)r}^{+}$ and$C_{ir}^{-}=R_{ir}{C}_{ir}^{+}$, and then comparing it with Eq. (25), the global wave transmission coefficient at discontinuity i can be found as

$$T_i={(1-r_i{R}_{ir})}^{-1}{t}_i{f}_{(i-1)}$$E25

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

$$C_1^{-}=R_{1r}{C}_1^{+}$$E26

However, recalling $C_1^{+}=r_1{C}_1^{-}$ from Eq. (24), it can be found that

$$(r_1{R}_{1r}-1)C_1^{+}=0$$E27

For nontrivial solutions,

$$F(s)=r_1{R}_{1r}-1=0$$E28

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_nlr_n. This equation simply states that when the string system vibrates at one of its natural frequencies, r₁R_1r1. As a simple example, for a single span uniformly damped string fixed at both ends, r₁r₂1 from Eq. (18) and R_1re²ⁱ from Eq. (20). Therefore, F(s)e²ⁱ10 gives the natural wavenumbers _nn (n1,2,3,...), or in terms of the non-dimensional frequency ${\omega }_n^2+i{c}_e{\omega }_n={(n\pi )}^2$ by Eq. (6) with s replaced by i.

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

$${\overline{w}}_i({\xi }_i)=f_i({\xi }_i)C_i^{+}+f_i^{-1}({\xi }_i)C_i^{-}$$E29

However, due to Eq. (19)

$${\overline{w}}_i({\xi }_i)=[f_i({\xi }_i)+f_i^{-1}({\xi }_i)R_{ir}]C_i^{+}$$E30

Now by applying the global wave transmission coefficient defined in Eq. (25) to $C_i^{+}$ in the above equation, the displacement at any point in span i of the string can be expressed in terms of the wave amplitude in span i1 as

$${\overline{w}}_i({\xi }_i)=[f_i({\xi }_i)+f_i^{-1}({\xi }_i)R_{ir}]T_i{C}_{i-1}^{+}$$E31

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$C_1^{+}$; i.e.,

$${\overline{w}}_i({\xi }_i)=[f_i({\xi }_i)+f_i^{-1}({\xi }_i)R_{ir}]\underset{j=i}{\overset{1}{\Pi }}{T}_j{C}_1^{+} C_1^{+} 0\le {\xi }_i\le l_i$$E32

Figure 3.

Plot of the characteristic equation. The solid and dashed curves represent the real and imaginary parts, respectively.

For example, consider a fixed-fixed undamped string with three supports specified by ₂=5+0.1s², ₃=7+0.1s², and ₄=4+0.1s² according to Eq. (13). l₁0.25, l₂0.3, l₃0.25, and l₄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.

Figure 4.

First three mode shapes obtained from Eq. (33). The solid, dashed, and dotted curves represent the 1^st, 2^nd, and 3^rd modes, respectively.

2.4. Transfer function analysis

Consider a multi-span string subjected to an external point load$\overline{p}(s)$, normalized against tension T, applied at xx₀ 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 $C_1^{\pm }$ and $D_1^{\pm }$ denote the injected waves that travel in the region x<x₀ and xx₀ within span 1, respectively. The transverse displacement of span 1 of the string can be expressed as

Figure 5.

Wave motion in a multi-span string due to a point load.

$${\overline{w}}_1({\xi }_1,x_0;s)=f_1({\xi }_1;s)C_1^{+}+f_1^{-1}({\xi }_1;s)C_1^{-} 0\le {\xi }_1\le l_1$$E33

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

$$C_1^{+}-D_1^{-}=v{D}_1^{+}-C_1^{-}=v$$E34

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

$$v=\frac{\sqrt{-1}}{2\gamma }\overline{p}(s)$$E35

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

$$C_1^{-}=R_{1r}{C}_1^{+}{D}_1^{-}=R_{1r}^{*}{D}_1^{+}$$E36

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

$$C_1^{+}=T_{1}v{T}_1=(1-R_{1r}^{*}{R}_{1r})(1+R_{1r}^{*})$$E37

$$D_1^{+}=T_1^{*}v{T}_1^{*}=(1-R_{1r}{R}_{1r}^{*})(1+R_{1r})$$E38

Note that $T_1$ and $T_1^{*}$ and can be considered as the global wave transmission coefficients that characterize the transmissibility of the wave injected by the external force in the region x<x₀ and xx₀, respectively. It is evident that and $T_1$ and $T_1^{*}$ are different unless the string system is symmetric about the loading point. Applying the results in Eq. (37) and Eq. (38) to Eq. (34), the wave motion in either side of x=x₀ can be found; i.e., in the region x<x₀

$${\overline{w}}_1({\xi }_1,x_0;s)=[f_1({\xi }_1;s)+f_1^{-1}({\xi }_1;s)R_{1r}]T_{1}v0\le {\xi }_1\le l_1$$E39

Now, in the same manner as for the mode shape analysis, since $C_{i+1}^{+}=T_i{C}_i^{+}$ and$C_1^{+}=T_{1}v$, the wave motion in span i on either side of the loading point can be found. For the region x<x₀

$${\overline{w}}_i({\xi }_i,x_0;s)=[f_i({\xi }_i;s)+f_i^{-1}({\xi }_i;s)R_{ir}]\underset{k=i}{\overset{1}{\Pi }}{T}_{k}v0\le {\xi }_i\le l_i$$E40

Note that the ratio

$$G_i({\xi }_i,x_0;s)={\overline{w}}_i({\xi }_i,x_0;s)/\overline{p}(s)$$E41

is the transfer function governing the forced response of any point in span i due to the point loading at xx₀. The Laplace inversion of G_i(_i, x₀; s) is the Green’s function of the problem. Thefore, denoting $L^{-1}$ as the inverse Laplace transform operator, the forced response at any point within any subspan of the multi-span string can be determined from the following convolution integral; e.g., for span i in the region x<x₀

$$w_i({\xi }_i,x_0;\tau )={\displaystyle {\int }_0^{\tau }{G}_i({\xi }_i,x_0;\tau -\overline{\tau })p(\tau )d\overline{\tau }}(i=1,2,\cdots ,m)$$E42

$$G_i({\xi }_i,x_0;\tau )=L^{-1}[\frac{\sqrt{-1}}{2\gamma }[f_i({\xi }_i;s)+f_i^{-1}({\xi }_i;s)R_{ir}]\underset{k=i}{\overset{1}{\Pi }}{T}_1]0\le {\xi }_i\le l_i$$E43

The exact Laplace inversion of G_i(_i, x₀; 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

$$w(x,x_0;\tau )={\displaystyle \sum _{n=1}^N\frac{2\sin n\pi x_0}{n\pi }\sin n\pi x\sin n\pi \tau }$$E44

The corresponding transfer function from Eq. (40) is

$$\overline{w}(x,x_0;s)=\frac{1}{2s(e^{2s}-1)}\{\begin{array}{cc}{e}^{-sx}(e^{2s(1-x_0)}-1)(e^{2sx}-e^{2s{x}_0})& \text{for}0\le x\le x_0\\ e^{-sx}(e^{2s{x}_0}-1)(e^{2sx}+e^{2s(1-x_0)})& \text{for}{x}_0<x\le 1\end{array}$$E45

Figure 6.

Unit impuse response of a single span string near the point of loading at 0.001: (a) M32 and N2,500; (b) M64 and N20,000. The solid and dashed curves represent the solutions from Eq. (43) and Eq. (44), respectively.

Shown in Fig. 6 is the comparison of the response solutions given in Eq. (43) and Eq. (44) when 0.001 near xx₀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$p(\tau )=p_0{e}^{i\omega \tau }$; i.e., a harmonic forcing function, the steady-state response of the problem can be readily found in terms of the complex frequency function defined as

$$H({\xi }_i,x_0;\omega )=G({\xi }_i,x_0;s){|}_{s=i\omega }$$E46

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₀

$$w_i({\xi }_i,x_0;\tau )=|H_i({\xi }_i,x_0;\omega )|p_0{e}^{i(\omega \tau +{\varphi }_i)}0\le {\xi }_i\le l_i$$E47

$$H_i({\xi }_i,x_0;\omega )=\frac{\sqrt{-1}}{2\gamma }[f_i({\xi }_i;\omega )+f_i^{-1}({\xi }_i;\omega )R_{ir}]\underset{k=i}{\overset{1}{\Pi }}{T}_k{\varphi }_i=\angle H_i$$E48

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

$$m\frac{{\partial }^{2}W}{\partial t^2}+C_e\frac{\partial W}{\partial t}+EI\frac{{\partial }^{4}W}{\partial X^4}=P(X,t)$$E49

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

$$w=W/L,x=X/L,\tau =t/t_0,t_0=m{L}^4/EI,c_e=C_e{t}_0/m,p(x,t)=P(X,t)L^3/EI$$E50

the equation of motion takes the non-dimensional form of

$$\ddot{w}+c_e\dot{w}+w^{(4)}=p(x,t)(0\le x\le 1)$$E51

Applying the Laplace transform to Eq. (48) yields

$$s^2\overline{w}(x;s)+c_{e}s\overline{w}(x;s)+{\overline{w}}^{(4)}(s;x)=\overline{p}(x;s)$$E52

Letting$\overline{p}(x;s)=0$, the homogeneous wave solution of Eq. (49) can be assumed as:

$$\overline{w}(x;s)=C{e}^{i\gamma x}$$E53

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

$${\gamma }^4=-(s^2+c_{e}s)$$E54

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

$$\overline{w}(x;s)=C_a^{+}{e}^{-i\gamma x}+C_a^{-}{e}^{i\gamma x}+C_b^{+}{e}^{-\gamma x}+C_b^{-}{e}^{\gamma x}$$E55

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.,

$$C^{+}=\left\{\begin{array}{c}{C}_a^{+}\\ C_b^{+}\end{array}\right\} C^{-}=\left\{\begin{array}{c}{C}_a^{-}\\ C_b^{-}\end{array}\right\}$$E56

and then

$$\overline{w}(x;s)=[\begin{array}{cc}1& 1\end{array}][f(x;s)C^{+}+f^{-1}(x;s)C^{-}]$$E57

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

$$f(x;s)=\left[\begin{array}{cc}{e}^{-i\gamma x}& 0\\ 0& e^{-\gamma x}\end{array}\right]$$E58

which relates the wave amplitudes by

$$C^{+}(x+x_0)=f{C}^{+}(x_0)or{C}^{-}(x+x_0)=f^{-1}{C}^{-}(x_0)$$E59

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

$$\left[\begin{array}{cc}1& 1\\ -i& -1\end{array}\right]C^{+}+\left[\begin{array}{cc}1& 1\\ i& 1\end{array}\right]r{C}^{+}=\left[\begin{array}{cc}1& 1\\ -i& -1\end{array}\right]t{C}^{+}$$E60

$$\left[\begin{array}{cc}1& -1\\ i& -1\end{array}\right]C^{+}+\left[\begin{array}{cc}1& -1\\ -i& 1\end{array}\right]r{C}^{+}=\left[\begin{array}{cc}1-i{k}_r/\gamma & -1-k_r/\gamma \\ i+(k_t+c_{t}s)/{\gamma }^3& -1+(k_t+c_{t}s)/{\gamma }^3\end{array}\right]t{C}^{+}$$E61

where$k_t=K_t{L}^3/EI$, $k_r=K_{r}L/EI$, and$c_t=C_t{t}_0/m$. Solving the above equations gives the local wave reflection and transmission matrices as:

$$r=\frac{1}{2}\left[\begin{array}{cc}(1-i)(\alpha -\beta i)& -(1+i)(\alpha +\beta )\\ -(1-i)(\alpha +\beta )& (1+i)(\alpha +\beta i)\end{array}\right]$$E62

$$t=\frac{1}{2}\left[\begin{array}{cc}2-(1-i)\alpha -(1+i)\beta & (1+i)(\alpha -\beta )\\ (1-i)(\alpha -\beta )& 2-(1+i)\alpha -(1-i)\beta \end{array}\right]$$E63

$$\alpha =\frac{k_r+c_{r}s}{k_r+c_{r}s+(2+2i)\gamma } \beta =\frac{k_t+c_{t}s}{k_t+c_{t}s-(2-2i){\gamma }^3}$$E64

Figure 7.

Wave reflection and transmission at a discontinuity.

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

$$C_{ir}^{-}=R_{ir}{C}_{ir}^{+}$$E65

Since$C_{ir}^{-}=f_i{R}_{(i+1)l}{C}_{ir}^{+}$, one can find R_ir in terms of the global wave reflection matrix on the left side of discontinuity i+1; i.e.,

$$R_{ir}=f_i{R}_{(i+1)l}{f}_i$$E66

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

$$C_{ir}^{+}=t_i{C}_{il}^{+}+r_i{C}_{ir}^{-} C_{il}^{-}=t_i{C}_{ir}^{-}+r_i{C}_{il}^{+}$$E67

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

$$R_{il}=r_i+t_i{(R_{ir}^{-1}-r_i)}^{-1}{t}_i$$E68

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

$$C_1^{+}=r_1{C}_1^{-} C_n^{-}=r_n{C}_n^{+}$$E69

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

$$r={\left[\begin{array}{rr}\hfill \gamma -i{k}_r& \hfill -\gamma -k_r\\ \hfill -i{\gamma }^3-(k_t+c_{t}s)& \hfill {\gamma }^3-(k_t+c_{t}s)\end{array}\right]}^{-1}\left[\begin{array}{rr}\hfill \gamma +i{k}_r& \hfill -\gamma +k_r\\ \hfill i{\gamma }^3-(k_t+c_{t}s)& \hfill -{\gamma }^3-(k_t+c_{t}s)\end{array}\right]$$E70

Now, to determine the global wave transmission matrix T_i, define

$$C_{ir}^{+}=T_i{C}_{(i-1)r}^{+}$$E71

Rewriting Eq. (60) by applying $C_{il}^{+}=f_{(i-1)}{C}_{(i-1)r}^{+}$ and$C_{ir}^{-}=R_{ir}{C}_{ir}^{+}$, and then comparing it with Eq. (68), one can find that

$$T_i={(I_{2\times 2}-r_i{R}_{ir})}^{-1}{t}_i{f}_{(i-1)}$$E72

where I₂₂ 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$C_i^{\pm }$, $R_i$, and $f_i$ are now replaced by$C_i^{\pm }$, $R_i$, and$f_i$, respectively, at the left boundary

$$C_1^{-}=R_{1r}{C}_1^{+}$$E73

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

$$(r_1{R}_{1r}-I_{2\times 2})C_1^{+}=0$$E74

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

$$F(s)=\text{Det}[r_1{R}_{1r}-I_{2\times 2}]=0$$E75

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

$${\overline{w}}_i({\xi }_i)=[\begin{array}{cc}1& 1\end{array}][f_i({\xi }_i)C_i^{+}+f_i^{-1}({\xi }_i)C_i^{-}]$$E76

However, due to Eq. (55)

$${\overline{w}}_i({\xi }_i)=[\begin{array}{cc}1& 1\end{array}][f_i({\xi }_i)+f_i^{-1}({\xi }_i)R_{ir}]C_i^{+}$$E77

Since $C_i^{+}=T_i{C}_{i-1}^{+}$ from Eq. (65),

$${\overline{w}}_i({\xi }_i)=[\begin{array}{cc}1& 1\end{array}][f_i({\xi }_i)+f_i^{-1}({\xi }_i)R_{ir}]T_i{C}_{i-1}^{+}$$E78

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$C_1^{+}$; i.e.,

$${\overline{w}}_i({\xi }_i)=[\begin{array}{cc}1& 1\end{array}][f_i({\xi }_i)+f_i^{-1}({\xi }_i)R_{ir}]\underset{j=i}{\overset{1}{\Pi }}{T}_j{C}_1^{+}{T}_1=I_{2\times 2}0\le {\xi }_i\le l_i$$E79

	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

Table 1.

Nondimensional system parameter used for Fig. 8.

where note that the amplitude ratio between the two wave components $C_a^{+}$ and $C_b^{+}$ can be found from Eq. (59). For example, shown in Fig. 8 are the first three mode shapes of a uniformly damped five-span beam with system parameters specified in Table 1. Once the wave reflection and transmission matrices at each discontinuity and the boundary are determined, one can apply Eq. (60) to find the first three natural wavenumbers ₁=10.294, ₂=12.038, and ₃=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.

Figure 8.

First three mode shapes obtained from Eq. (52). The solid, dashed, and dotted curves represent the 1^st, 2^nd, and 3^rd modes, respectively.

3.3. Transfer function analysis

Consider a multi-span beam subjected to an external force applied at x=x₀ where x₀ is located in subspan m. Let $C^{\pm }$ be the amplitudes of the waves rise in sub-span n as a result of injected waves due to the applied force, and also assume that $C^{\pm }$ satisfy all the continuity conditions at intermediate discontinuities and boundary conditions of the multi-span beam system. The transverse displacement ${\overline{w}}_n(x;s)$ of span n can be expressed in wave form

$${\overline{w}}_n(x;s)=[\begin{array}{cc}1& 1\end{array}][f(x;s)C^{+}+f^{-1}(x;s)C^{-}]$$E80

Now, in order to determine the actual wave amplitudes$$C^{\pm }$$, consider the multi-span beam with arbitrary supports and boundary conditions under a concentrated applied force of$$\overline{p}(s)\delta (x-x_0)$$, where $$\overline{p}(s)$$ is the Laplace transform of p(), as schematically depicted in Fig. 5 with$C_i^{\pm }$, $D_i^{\pm }$, and $R_i$ replaced by$C_i^{\pm }$, $D_i^{\pm }$, and$R_i$, respectively. Since the waves injected at x=x₀ 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 $$C_i^{\pm }$$ and $$D_i^{\pm }$$ be the amplitudes of the waves traveling within subspan i in the region x<x₀ and xx₀, respectively. The discontinuity in the shear force at x=x₀ can be expressed in term of the difference in amplitudes between the incoming and outgoing waves at the discontinuity such that

$$D_1^{+}-C_1^{-}=v C_1^{+}-D_1^{-}=v$$E81

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

$$v=\frac{\overline{p}(s)}{4{\gamma }^3}\left\{\begin{array}{c}\sqrt{-1}\\ 1\end{array}\right\}$$E82

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

$$C_1^{-}=R_{1r}{C}_1^{+}{D}_1^{-}=R_{1r}^{*}{D}_1^{+}$$E83

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

$$C_1^{+}=T_{1}v T_1={(I-R_{1r}^{\ast }{R}_{1r})}^{-1}(I+R_{1r}^{\ast })$$E84

$$D_1^{+}=T_1^{*}v T_1^{*}={(I-R_{1r}{R}_{1r}^{*})}^{-1}(I+R_{1r})$$E85

$T_1$and $T_1^{*}$ in the above equations can be considered as the generalized wave transmission matrices that characterize the transmissibility of the waves injected by the applied force in the region x<x₀ and xx₀, respectively. It is evident that $T_1$ and $T_1^{*}$ are different unless the system is symmetric about x=x₀. 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₀

$${\overline{w}}_1({\xi }_1,x_0;s)=[\begin{array}{cc}1& 1\end{array}][f_1({\xi }_1;s)+f_1^{-1}({\xi }_1;s)R_{1r}]T_{1}v0\le {\xi }_1\le l_1$$E86

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₀, $C_1^{+}=T_{1}v$and$C_{i+1}^{+}=T_i{C}_i^{+}$. Therefore one can progressively expand the solution in Eq. (67) until the expansion terminates at the leftmost boundary. Thus, the transverse displacement of any point within span i in the region x<x₀ due to external loading at x=x₀ can be determined from

$${\overline{w}}_i({\xi }_i,x_0;s)=\frac{\overline{p}(s)}{4{\gamma }^3}[\begin{array}{cc}1& 1\end{array}][f_i({\xi }_i;s)+f_i^{-1}({\xi }_i;s)R_{ir}]{\displaystyle \prod _{k=i}^1{T}_k}\left\{\begin{array}{c}\sqrt{-1}\\ 1\end{array}\right\}0\le {\xi }_i\le l_i$$E87

One can find the same expression in the region xx₀, but in terms of $R_{ir}^{*}$ and $T_k^{*}$ instead. Note that the ratio $$G_i(\xi ,x_0;s)=\overline{w}({\xi }_i,x_0;s)/\overline{p}(s)$$ in Eq. (62) is the transfer function governing the forced response of anypoint in span i due to the point loading at x=x₀. The Laplace inversion of $$G_i(\xi ,x_0;s)$$ is the Green’s function of the problem. Therefore, denoting $L^{-1}$ as the inverse Laplace transform operator, the forced response at any point within any span of the multi-span beam can be found from; e.g., for span i in the region x<x₀,

$$w_i({\xi }_i,x_0;\tau )={\displaystyle {\int }_0^t{G}_i({\xi }_i,x_0;\tau -\overline{\tau })p(\tau )d\overline{\tau }}(i=1,2,\cdots ,m)$$E88

$$G_i({\xi }_i,x_0;\tau )=L^{-1}\left[\frac{1}{4{\gamma }^3}[\begin{array}{cc}1& 1\end{array}][f_i({\xi }_i;s)+f_i^{-1}({\xi }_i;s)R_{ir}]{\displaystyle \prod _{k=i}^1{T}_k}\left\{\begin{array}{c}\sqrt{-1}\\ 1\end{array}\right\}\right]0\le {\xi }_i\le l_i$$E89

The exact Laplace inversion of G_i(, x₀; 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.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

Table 2.

Nondimensional system parameter used for Fig. 9, where c_e=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

Figure 9.

Unit impulse response of the three-span beam sampled at (a) x=0.15 and (b) x=0.825. The solid and dashed curves represent the results from the present analysis and the method of normal mode expansion, respectively.

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$p(\tau )=p_0{e}^{i\omega \tau }$; i.e., a harmonic forcing function, the steady-state response of the problem can be readily found in terms of the complex frequency function defined as

$$H({\xi }_i,x_0;\omega )=G({\xi }_i,x_0;s){|}_{s=i\omega }$$E90

		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

Table 3.

Nondimensional system parameter used for Fig. 8, where c_e=0.

Figure 10.

Frequency response of the six-span beam due to harmonic excitation applied at x=0.39: (a) amplitudes sampled at x=0.325 (solided curve) and x=0.92 (dashed curve); and (b) corresponding phase angle.

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₀

$$w_i({\xi }_i,x_0;\tau )=|H_i({\xi }_i,x_0;\omega )|p_0{e}^{i(\omega \tau +{\varphi }_i)}0\le {\xi }_i\le l_i$$E91

$$H_i({\xi }_i,x_0;\omega )=\frac{1}{4{\gamma }^3}[\begin{array}{cc}1& 1\end{array}][f_i({\xi }_i;\omega )+f_i^{-1}({\xi }_i;\omega )R_{ir}]\underset{k=i}{\overset{1}{\Pi }}{T}_k\left\{\begin{array}{c}\sqrt{-1}\\ 1\end{array}\right\}{\varphi }_i=\angle H_i$$E92

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

$$\frac{EI}{R^3}\frac{{\partial }^3}{\partial {\theta }^3}\left(U-\frac{\partial W}{\partial \theta }\right)-\frac{EA}{R}\left(W+\frac{\partial U}{\partial \theta }\right)=\rho AR\frac{{\partial }^{2}W}{\partial t^2}$$E93

$$\frac{EI}{R^3}\frac{{\partial }^2}{\partial {\theta }^2}\left(U-\frac{\partial W}{\partial \theta }\right)+\frac{EA}{R}\frac{\partial }{\partial \theta }\left(W+\frac{\partial U}{\partial \theta }\right)=\rho AR\frac{{\partial }^{2}U}{\partial t^2}$$E94

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:

$$u=\frac{U}{R}w=\frac{W}{R}\tau =\frac{t}{t_0}{t}_0=R^2\sqrt{\frac{\rho A}{EI}}{k}^2=\frac{I}{A{R}^2}$$E95

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

$$k^2\frac{{\partial }^3}{\partial {\theta }^3}\left(\overline{u}(\theta ;s)-\frac{\partial \overline{w}(\theta ;s)}{\partial \theta }\right)-\left(\overline{w}(\theta ;s)+\frac{\partial u(\theta ;s)}{\partial \theta }\right)=k^2{s}^2\overline{w}(\theta ;s)$$E96

$$k^2\frac{{\partial }^2}{\partial {\theta }^2}\left(\overline{u}(\theta ;s)-\frac{\partial \overline{w}(\theta ;s)}{\partial \theta }\right)+\frac{\partial }{\partial \theta }\left(\overline{w}(\theta ;s)+\frac{\partial \overline{u}(\theta ;s)}{\partial \theta }\right)=k^2{s}^2\overline{u}(\theta ;s)$$E97

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

$$\left\{\begin{array}{c}\overline{w}(\theta ;s)\\ \overline{u}(\theta ;s)\end{array}\right\}=\left\{\begin{array}{c}{C}_w\\ C_u\end{array}\right\}{e}^{i\gamma \theta }$$E98

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

$$\left[\begin{array}{cc}{k}^2({\gamma }^4+s^2)+1& i\gamma (1+{\gamma }^2{k}^2)\\ i\gamma (1+{\gamma }^2{k}^2)& -k^2({\gamma }^2+s^2)-{\gamma }^2\end{array}\right]\left\{\begin{array}{c}{C}_w\\ C_u\end{array}\right\}=0$$E99

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

$${\gamma }^6+(k^2{s}^2-2){\gamma }^4+[1+(1+k^2)s^2]{\gamma }^2+(k^2{s}^2+1)s^2=0$$E100

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

$$\overline{w}(\theta ;s)={\displaystyle \sum _{n=1}^3(C_n^{+}{e}^{-i{\gamma }_n\theta }+C_n^{-}{e}^{i{\gamma }_n\theta }\text{) }}$$E101

$$\overline{u}(\theta ;s)={\displaystyle \sum _{n=1}^3{\alpha }_n(C_n^{+}{e}^{-i{\gamma }_n\theta }-C_n^{-}{e}^{i{\gamma }_n\theta }\text{) }}$$E102

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

$${\alpha }_n=\frac{i{\gamma }_n(k^2{\gamma }_n^2+1)}{(1+k^2){\gamma }_n^2+k^2{s}^2}$$E103

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., $$C^{+}={\left\{\begin{array}{ccc}{C}_1^{+}& C_2^{+}& C_3^{+}\end{array}\right\}}^T$$and$$C^{-}={\left\{\begin{array}{ccc}{C}_1^{-}& C_2^{-}& C_3^{-}\end{array}\right\}}^T$$, then

$$\overline{w}(\theta ;s)=u[f{C}^{+}+f^{-1}{C}^{-}]$$E104

$$\overline{u}(\theta ;s)=\alpha [f{C}^{+}-f^{-1}{C}^{-}]$$E105

where$u=[\begin{array}{ccc}1& 1& 1\end{array}]$, $\alpha =[\begin{array}{ccc}{\alpha }_1& {\alpha }_2& {\alpha }_3]\end{array}$, and f is the field transfer matrix defined by

$$f(\theta ;s)=\left[\begin{array}{ccc}{e}^{-i{\gamma }_1\theta }& 0& 0\\ 0& e^{-i{\gamma }_2\theta }& 0\\ 0& 0& e^{-i{\gamma }_3\theta }\end{array}\right]$$E106

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 ${\gamma }_n$ and ${{\gamma }^{\prime }}_n$ (n=1, 2, 3), the geometric continuity conditions at =0 give

$$\begin{array}{l}\left[\begin{array}{ccc}1& 1& 1\\ {\alpha }_1& {\alpha }_2& {\alpha }_3\\ {\alpha }_1-i{\gamma }_1& {\alpha }_2-i{\gamma }_2& {\alpha }_3-i{\gamma }_3\end{array}\right]C^{+}+\left[\begin{array}{ccc}1& 1& 1\\ -{\alpha }_1& -{\alpha }_2& -{\alpha }_3\\ -{\alpha }_1+i{\gamma }_1& -{\alpha }_2+i{\gamma }_2& -{\alpha }_3+i{\gamma }_3\end{array}\right]r{C}^{+}\\ =\left[\begin{array}{ccc}1& 1& 1\\ {{\alpha }^{\prime }}_1& {{\alpha }^{\prime }}_2& {{\alpha }^{\prime }}_3\\ {{\alpha }^{\prime }}_1-i{{\gamma }^{\prime }}_1& {\alpha }_2-i{{\gamma }^{\prime }}_2& {{\alpha }^{\prime }}_3-i{{\gamma }^{\prime }}_3\end{array}\right]t{C}^{+}\end{array}$$E107

and the kinetic continuity conditions give

$$\begin{array}{l}\left[\begin{array}{ccc}1-i{\gamma }_1{\alpha }_1& 1-i{\gamma }_2{\alpha }_2& 1-i{\gamma }_3{\alpha }_3\\ -{\gamma }_1^2(i{\gamma }_1+{\alpha }_1)& -{\gamma }_2^2(i{\gamma }_2+{\alpha }_2)& -{\gamma }_3^2(i{\gamma }_3+{\alpha }_3)\\ {\gamma }_1({\gamma }_1-i{\alpha }_1)& {\gamma }_2({\gamma }_2-i{\alpha }_2)& {\gamma }_3({\gamma }_3-i{\alpha }_3)\end{array}\right]C^{+}\\ +\left[\begin{array}{ccc}1-i{\gamma }_1{\alpha }_1& 1-i{\gamma }_2{\alpha }_2& 1-i{\gamma }_3{\alpha }_3\\ {\gamma }_1^2(i{\gamma }_1+{\alpha }_1)& {\gamma }_2^2(i{\gamma }_2+{\alpha }_2)& {\gamma }_3^2(i{\gamma }_3+{\alpha }_3)\\ {\gamma }_1({\gamma }_1-i{\alpha }_1)& {\gamma }_2({\gamma }_2-i{\alpha }_2)& {\gamma }_3({\gamma }_3-i{\alpha }_3)\end{array}\right]r{C}^{+}\\ =\left[\begin{array}{ccc}1-i{{\gamma }^{\prime }}_1{{\alpha }^{\prime }}_1& 1-i{{\gamma }^{\prime }}_2{{\alpha }^{\prime }}_2& 1-i{{\gamma }^{\prime }}_3{{\alpha }^{\prime }}_3\\ -{{\gamma }^{\prime }}_1^2(i{{\gamma }^{\prime }}_1+{{\alpha }^{\prime }}_1)& -{{\gamma }^{\prime }}_2^2(i{{\gamma }^{\prime }}_2+{{\alpha }^{\prime }}_2)& -{{\gamma }^{\prime }}_3^2(i{{\gamma }^{\prime }}_3+{{\alpha }^{\prime }}_3)\\ {{\gamma }^{\prime }}_1({{\gamma }^{\prime }}_1-i{{\alpha }^{\prime }}_1)& {{\gamma }^{\prime }}_2({{\gamma }^{\prime }}_2-i{{\alpha }^{\prime }}_2)& {{\gamma }^{\prime }}_3({{\gamma }^{\prime }}_3-i{{\alpha }^{\prime }}_3)\end{array}\right]t{C}^{+}\end{array}$$E108

$${{\alpha }^{\prime }}_n=\frac{i{{\gamma }^{\prime }}_n({k^{\prime }}^2{{\gamma }^{\prime }}_n^2+1)}{(1+{k^{\prime }}^2){{\gamma }^{\prime }}_n^2+{k^{\prime }}^2{s}^2}$$E109

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

$$F(s)=\text{Det}[r_1{R}_{1r}-I_{3\times 3}]=0$$E110

where r₁ is the wave reflection matrix at the clamed boundary, which is

$$r_1={\left[\begin{array}{ccc}-1& -1& -1\\ -{\alpha }_1& -{\alpha }_2& -{\alpha }_3\\ {\alpha }_1-i{\gamma }_1& {\alpha }_2-i{\gamma }_2& {\alpha }_3-i{\gamma }_3\end{array}\right]}^{-1}\left[\begin{array}{ccc}-1& -1& -1\\ -{\alpha }_1& -{\alpha }_2& -{\alpha }_3\\ {\alpha }_1-i{\gamma }_1& {\alpha }_2-i{\gamma }_2& {\alpha }_3-i{\gamma }_3\end{array}\right]$$E111

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.

Figure 11.

Curved beam with three subspan of equal span angle of 60, where k₁=k₃=0.01 and k₂=0.02.

	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

Table 4.

Non-dimensional natural frequencies of the curved beam in Fig. 11.

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.