Open access peer-reviewed chapter

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

By Bongsu Kang

Submitted: November 18th 2010Reviewed: April 21st 2011Published: September 9th 2011

DOI: 10.5772/22379

Downloaded: 2323

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

T2WX2+CeWt=m2Wt2E1

where W is the transverse displacement, X the spatial variable, t the temporal variable; and T denotes the constant tension, Ce 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  τ=t/t0  ce=CeL/t0  t0=L/c0  c0=T/mE2

the equation of motion takes the non-dimensional form of

w+cew˙=w¨  (0x1)E3

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

w¯(x;s)+cesw¯(x;s)=s2w¯(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

w¯(x;s)=CeiγxE5

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

γ2+s2ces=0orγ2=cess2E6

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

w¯(x;s)=C+eiγx+CeiγxE7

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 ce, hence the classification of the wave into propagating and attenuating waves does not apply to this case. Defining f(x;s)=eiγxas the field transfer function which relates the wave amplitudes by

C+(x+x0)=fC+(x0) or  C(x+x0)=f1C(x0)E8

wheref1(x;s)=eiγx, the wave solution in Eq. (7) can be re-written as

w¯(x;s)=fC++f1CE9

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 ξ=0as shown in Fig. 1, where the constraint is a point support consisting of an attached mass (Mc), a transverse spring (Kc), and a viscous damper (Cc).

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 rC+and a transmitted wavetC+, 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

w¯l=C++rC+w¯r=tC+E10

Since w¯l=w¯rat =0, one can find

1+r=tE11

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

w¯rw¯l=φw¯rE12

where represents the kinetic properties of the constraint as

φ(s)=kc+ccs+mcs2E13

Note that kc, cc, and mc in Eq. (13) are the non-dimensional spring constant, damping coefficient, and the attached mass, respectively, defined by

kc=KcL/Tcc=Ccc0/Tmc=Mcc02/TLE14

Equation (12) leads to

γ(r+t1)=iφtE15

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

r=iφ2γiφandt=2γ2γiφE16

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

w¯l=φw¯lE17

which leads to

r=γ+iφγiφ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 Rir 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

Cir=RirCir+E19

Figure 2.

Waves traveling along a multi-span string.

SinceCir=fiR(i+1)lCir+, one can find Rir in terms of the global wave reflection coefficient on the left side of discontinuity i+1; i.e.,

Rir=fi2R(i+1)lE20

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

Cir+=tiCil++riCirandCil=tiCir+riCil+E21

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

Ril=ri+ti2(Rir1ri)1E22

Rir and Ril progressively expand to include all the global wave reflection coefficients of discontinuities along the string before terminating its expansion at the boundaries where

C1+=r1C1Cn=rnCn+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 Ti, define

Cir+=TiC(i1)r+E24

Rewriting Eq. (21) by applying Cil+=f(i1)C(i1)r+andCir=RirCir+, and then comparing it with Eq. (25), the global wave transmission coefficient at discontinuity i can be found as

Ti=(1riRir)1tif(i1)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

C1=R1rC1+E26

However, recalling C1+=r1C1from Eq. (24), it can be found that

(r1R1r1)C1+=0E27

For nontrivial solutions,

F(s)=r1R1r1=0E28

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 R1r 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 Rnlrn. This equation simply states that when the string system vibrates at one of its natural frequencies, r1R1r1. As a simple example, for a single span uniformly damped string fixed at both ends, r1r21 from Eq. (18) and R1re2i from Eq. (20). Therefore, F(s)e2i10 gives the natural wavenumbers nn (n1,2,3,...), or in terms of the non-dimensional frequency ωn2+iceωn=(nπ)2by 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

w¯i(ξi)=fi(ξi)Ci++fi1(ξi)CiE29

However, due to Eq. (19)

w¯i(ξi)=[fi(ξi)+fi1(ξi)Rir]Ci+E30

Now by applying the global wave transmission coefficient defined in Eq. (25) to Ci+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

w¯i(ξi)=[fi(ξi)+fi1(ξi)Rir]TiCi1+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 amplitudeC1+; i.e.,

w¯i(ξi)=[fi(ξi)+fi1(ξi)Rir]Πj=i1TjC1+  C1+  0ξiliE32

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 2=5+0.1s2, 3=7+0.1s2, and 4=4+0.1s2 according to Eq. (13). l10.25, l20.3, l30.25, and l40.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 1st, 2nd, and 3rd modes, respectively.

2.4. Transfer function analysis

Consider a multi-span string subjected to an external point loadp¯(s), normalized against tension T, applied at xx0 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 C1±and D1±denote the injected waves that travel in the region x<x0 and xx0 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.

w¯1(ξ1,x0;s)=f1(ξ1;s)C1++f11(ξ1;s)C1  0ξ1l1E33

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

C1+D1=vD1+C1=vE34

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

v=12γp¯(s)E35

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

C1=R1rC1+D1=R1r*D1+E36

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

C1+=T1vT1=(1R1r*R1r)(1+R1r*)E37
D1+=T1*vT1*=(1R1rR1r*)(1+R1r)E38

Note that T1and T1*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<x0 and xx0, respectively. It is evident that and T1and T1*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=x0 can be found; i.e., in the region x<x0

w¯1(ξ1,x0;s)=[f1(ξ1;s)+f11(ξ1;s)R1r]T1v0ξ1l1E39

Now, in the same manner as for the mode shape analysis, since Ci+1+=TiCi+andC1+=T1v, the wave motion in span i on either side of the loading point can be found. For the region x<x0

w¯i(ξi,x0;s)=[fi(ξi;s)+fi1(ξi;s)Rir]Πk=i1Tkv0ξiliE40

Note that the ratio

Gi(ξi,x0;s)=w¯i(ξi,x0;s)/p¯(s)E41

is the transfer function governing the forced response of any point in span i due to the point loading at xx0. The Laplace inversion of Gi(i, x0; s) is the Green’s function of the problem. Thefore, denoting L1as 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<x0

wi(ξi,x0;τ)=0τGi(ξi,x0;ττ¯)p(τ)dτ¯(i=1,2,,m)E42
Gi(ξi,x0;τ)=L1[12γ[fi(ξi;s)+fi1(ξi;s)Rir]Πk=i1T1]0ξiliE43

The exact Laplace inversion of Gi(i, x0; 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,x0;τ)=n=1N2sinnπx0nπsinnπxsinnπτE44

The corresponding transfer function from Eq. (40) is

w¯(x,x0;s)=12s(e2s1){esx(e2s(1x0)1)(e2sxe2sx0)for0xx0esx(e2sx01)(e2sx+e2s(1x0))forx0<x1E45

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

Ifp(τ)=p0eiωτ; 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(ξi,x0;ω)=G(ξi,x0;s)|s=iω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<x0

wi(ξi,x0;τ)=|Hi(ξi,x0;ω)|p0ei(ωτ+ϕi)0ξiliE47
Hi(ξi,x0;ω)=12γ[fi(ξi;ω)+fi1(ξi;ω)Rir]Πk=i1Tkϕi=HiE48

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

m2Wt2+CeWt+EI4WX4=P(X,t)E49

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

w=W/L,x=X/L,τ=t/t0,t0=mL4/EI,ce=Cet0/m,p(x,t)=P(X,t)L3/EIE50

the equation of motion takes the non-dimensional form of

w¨+cew˙+w(4)=p(x,t)(0x1)E51

Applying the Laplace transform to Eq. (48) yields

s2w¯(x;s)+cesw¯(x;s)+w¯(4)(s;x)=p¯(x;s)E52

Lettingp¯(x;s)=0, the homogeneous wave solution of Eq. (49) can be assumed as:

w¯(x;s)=CeiγxE53

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

γ4=(s2+ces)E54

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

w¯(x;s)=Ca+eiγx+Caeiγx+Cb+eγx+CbeγxE55

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+={Ca+Cb+}  C={CaCb}E56

and then

w¯(x;s)=[11][f(x;s)C++f1(x;s)C]E57

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

f(x;s)=[eiγx00eγx]E58

which relates the wave amplitudes by

C+(x+x0)=fC+(x0)orC(x+x0)=f-1C(x0)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 (Kt) and rotational (Kr) springs and transverse damper (Ct), 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

[11i1]C++[11i1]rC+=[11i1]tC+E60
[11i1]C++[11i1]rC+=[1ikr/γ1kr/γi+(kt+cts)/γ31+(kt+cts)/γ3]tC+E61

wherekt=KtL3/EI, kr=KrL/EI, andct=Ctt0/m. Solving the above equations gives the local wave reflection and transmission matrices as:

r=12[(1i)(αβi)(1+i)(α+β)(1i)(α+β)(1+i)(α+βi)]E62
t=12[2(1i)α(1+i)β(1+i)(αβ)(1i)(αβ)2(1+i)α(1i)β]E63
α=kr+crskr+crs+(2+2i)γ  β=kt+ctskt+cts(22i)γ3E64

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 Rir 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

Cir=RirCir+E65

SinceCir=fiR(i+1)lCir+, one can find Rir in terms of the global wave reflection matrix on the left side of discontinuity i+1; i.e.,

Rir=fiR(i+1)lfiE66

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

Cir+=tiCil++riCir  Cil=tiCir+riCil+E67

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

Ril=ri+ti(Rir1ri)1tiE68

Rir and Ril 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

C1+=r1C1  Cn=rnCn+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=[γikrγkriγ3(kt+cts)γ3(kt+cts)]1[γ+ikrγ+kriγ3(kt+cts)γ3(kt+cts)]E70

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

Cir+=TiC(i1)r+E71

Rewriting Eq. (60) by applying Cil+=f(i1)C(i1)r+andCir=RirCir+, and then comparing it with Eq. (68), one can find that

Ti=(I2×2riRir)1tif(i1)E72

where I22 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, whereCi±, Ri, and fiare now replaced byCi±, Ri, andfi, respectively, at the left boundary

C1=R1rC1+E73

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

(r1R1rI2×2)C1+=0E74

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)=Det[r1R1rI2×2]=0E75

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

w¯i(ξi)=[11][fi(ξi)Ci++fi1(ξi)Ci]E76

However, due to Eq. (55)

w¯i(ξi)=[11][fi(ξi)+fi1(ξi)Rir]Ci+E77

Since Ci+=TiCi1+from Eq. (65),

w¯i(ξi)=[11][fi(ξi)+fi1(ξi)Rir]TiCi1+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 amplitudeC1+; i.e.,

w¯i(ξi)=[11][fi(ξi)+fi1(ξi)Rir]Πj=i1TjC1+T1=I2×20ξiliE79
Discontinuity
Constraint123456
Location00.150.350.60.821
kt3,0002,5001,5003,500
kr0250150100300
ct000000

Table 1.

Nondimensional system parameter used for Fig. 8.

where note that the amplitude ratio between the two wave components Ca+and Cb+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 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.

Figure 8.

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

3.3. Transfer function analysis

Consider a multi-span beam subjected to an external force applied at x=x0 where x0 is located in subspan m. Let C±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±satisfy all the continuity conditions at intermediate discontinuities and boundary conditions of the multi-span beam system. The transverse displacement w¯n(x;s)of span n can be expressed in wave form

w¯n(x;s)=[11][f(x;s)C++f1(x;s)C]E80

Now, in order to determine the actual wave amplitudesC±, consider the multi-span beam with arbitrary supports and boundary conditions under a concentrated applied force ofp¯(s)δ(xx0), where p¯(s)is the Laplace transform of p(), as schematically depicted in Fig. 5 withCi±, Di±, and Rireplaced byCi±, Di±, andRi, respectively. Since the waves injected at x=x0 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 Ci±and Di±be the amplitudes of the waves traveling within subspan i in the region x<x0 and xx0, respectively. The discontinuity in the shear force at x=x0 can be expressed in term of the difference in amplitudes between the incoming and outgoing waves at the discontinuity such that

D1+C1=v    C1+D1=vE81

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=p¯(s)4γ3{11}E82

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

C1=R1rC1+D1=R1r*D1+E83

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

C1+=T1v    T1=(IR1rR1r)1(I+R1r)E84
D1+=T1*v    T1*=(IR1rR1r*)1(I+R1r)E85

T1and T1*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<x0 and xx0, respectively. It is evident that T1and T1*are different unless the system is symmetric about x=x0. 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<x0

w¯1(ξ1,x0;s)=[11][f1(ξ1;s)+f11(ξ1;s)R1r]T1v0ξ1l1E86

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<x0, C1+=T1vandCi+1+=TiCi+. 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<x0 due to external loading at x=x0 can be determined from

w¯i(ξi,x0;s)=p¯(s)4γ3[11][fi(ξi;s)+fi1(ξi;s)Rir]k=i1Tk{11}0ξiliE87

One can find the same expression in the region xx0, but in terms of Rir*and Tk*instead. Note that the ratio Gi(ξ,x0;s)=w¯(ξi,x0;s)/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=x0. The Laplace inversion of Gi(ξ,x0;s)is the Green’s function of the problem. Therefore, denoting L1as 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<x0,

wi(ξi,x0;τ)=0tGi(ξi,x0;ττ¯)p(τ)dτ¯(i=1,2,,m)E88
Gi(ξi,x0;τ)=L1[14γ3[11][fi(ξi;s)+fi1(ξi;s)Rir]k=i1Tk{11}]0ξiliE89

The exact Laplace inversion of Gi(, x0; 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 x0= 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
Constraint321=1*23
Location00.150.40.651
kt2,00003,000
kr00000
ct0150100

Table 2.

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

Ifp(τ)=p0eiωτ; 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(ξi,x0;ω)=G(ξi,x0;s)|s=iωE90
Discontinuity
Constraint4321=1*2*3*4*5*
Location00.120.260.390.560.690.841
kt5,0004,0003,50003,0005,0000
kr100010001005000
ct0.50.10000.10.10.2

Table 3.

Nondimensional system parameter used for Fig. 8, where ce=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<x0

wi(ξi,x0;τ)=|Hi(ξi,x0;ω)|p0ei(ωτ+ϕi)0ξiliE91
Hi(ξi,x0;ω)=14γ3[11][fi(ξi;ω)+fi1(ξi;ω)Rir]Πk=i1Tk{11}ϕi=HiE92

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

EIR33θ3(UWθ)EAR(W+Uθ)=ρAR2Wt2E93
EIR32θ2(UWθ)+EARθ(W+Uθ)=ρAR2Ut2E94

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=URw=WRτ=tt0t0=R2ρAEIk2=IAR2E95

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

k23θ3(u¯(θ;s)w¯(θ;s)θ)(w¯(θ;s)+u(θ;s)θ)=k2s2w¯(θ;s)E96
k22θ2(u¯(θ;s)w¯(θ;s)θ)+θ(w¯(θ;s)+u¯(θ;s)θ)=k2s2u¯(θ;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

{w¯(θ;s)u¯(θ;s)}={CwCu}eiγθ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

[k2(γ4+s2)+1iγ(1+γ2k2)iγ(1+γ2k2)k2(γ2+s2)γ2]{CwCu}=0E99

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

γ6+(k2s22)γ4+[1+(1+k2)s2]γ2+(k2s2+1)s2=0E100

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

w¯(θ;s)=n=13(Cn+eiγnθ+CneiγnθE101
u¯(θ;s)=n=13αn(Cn+eiγnθCneiγnθ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

αn=iγn(k2γn2+1)(1+k2)γn2+k2s2E103

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+={C1+C2+C3+}TandC={C1C2C3}T, then

w¯(θ;s)=u[fC++f1C]E104
u¯(θ;s)=α[fC+f1C]E105

whereu=[111], α=[α1α2α3], and f is the field transfer matrix defined by

f(θ;s)=[eiγ1θ000eiγ2θ000eiγ3θ]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 γnand γn(n=1, 2, 3), the geometric continuity conditions at =0 give

[111α1α2α3α1iγ1α2iγ2α3iγ3]C++[111α1α2α3α1+iγ1α2+iγ2α3+iγ3]rC+=[111α1α2α3α1iγ1α2iγ2α3iγ3]tC+E107

and the kinetic continuity conditions give

[1iγ1α11iγ2α21iγ3α3γ12(iγ1+α1)γ22(iγ2+α2)γ32(iγ3+α3)γ1(γ1iα1)γ2(γ2iα2)γ3(γ3iα3)]C++[1iγ1α11iγ2α21iγ3α3γ12(iγ1+α1)γ22(iγ2+α2)γ32(iγ3+α3)γ1(γ1iα1)γ2(γ2iα2)γ3(γ3iα3)]rC+=[1iγ1α11iγ2α21iγ3α3γ12(iγ1+α1)γ22(iγ2+α2)γ32(iγ3+α3)γ1(γ1iα1)γ2(γ2iα2)γ3(γ3iα3)]tC+E108
αn=iγn(k2γn2+1)(1+k2)γn2+k2s2E109

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)=Det[r1R1rI3×3]=0E110

where r1 is the wave reflection matrix at the clamed boundary, which is

r1=[111α1α2α3α1iγ1α2iγ2α3iγ3]1[111α1α2α3α1iγ1α2iγ2α3iγ3]E111

and R1r 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 k1=k3=0.01 and k2=0.02.

PresentFEMCDM
12 curvedelements24 straight elements100 Degreesof freedom
ModeRef. 1Ref. 2Ref. 1Ref. 1
12.6832.6802.7012.6852.671
24.8344.8244.8284.8284.780
39.5659.5369.5439.5439.413
414.58514.52714.53514.53514.309
521.86521.74921.75121.75121.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.

© 2011 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution-NonCommercial-ShareAlike-3.0 License, which permits use, distribution and reproduction for non-commercial purposes, provided the original is properly cited and derivative works building on this content are distributed under the same license.

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

Bongsu Kang (September 9th 2011). Exact Transfer Function Analysis of Distributed Parameter Systems by Wave Propagation Techniques, Recent Advances in Vibrations Analysis, Natalie Baddour, IntechOpen, DOI: 10.5772/22379. Available from:

chapter statistics

2323total chapter downloads

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

Related Content

This Book

Next chapter

Phase Diagram Analysis for Predicting Nonlinearities and Transient Responses

By Juan Carlos Jáuregui

Related Book

First chapter

Theory of Tribo-Systems

By Xie You-Bai

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More About Us