Modified transfer matrix method for steady-state forced vibration: a system of beam elements∗

Abstract. The EST (Elements by a System of Transfer equations) method offers exact solutions for various vibration problems of trusses, beams and frames. The method can be regarded as an improved or modified transfer matrix method where the roundoff errors generated by multiplying transfer arrays are avoided. It is assumed that in a steady state a beam will vibrate with the circular frequency of an excitation force. The universal equation of elastic displacement (4th order differential equation) is described as a system of first order differential equations in matrix form. For the differential equations, the compatibility conditions of a beam element displacements at joint serve as essential boundary conditions. As the natural boundary conditions at joints, the equilibrium equations of elastic forces of beam elements are considered. At the supports, restrictions to displacements (support conditions) have been applied. For steady-state forced vibration, the phenomena of dynamic vibration absorption near the saddle points are observed, and the response curves for displacement amplitude and elastic energy are calculated.


INTRODUCTION
One of the problems in structural engineering has been predicting the response of a structure or mechanical system to external steady-state forced vibration [2].Two phenomena, resonance and dynamic vibration absorption, have been of great interest [3,4].
In computational structural mechanics the state-space representation of mechanical systems can be seen as an application of the transfer matrix method [5][6][7][8][9]: where Z A , Z L designate the components (displacements, internal forces) of the state vectors at the beginning (x = 0) and end (x = ) of the element; Z p is the element loading vector; U denotes the transfer matrix.
An external force can be added as an element load (described by a forcing function) or as a nodal load on joints.The general solution to an ordinary differential equation can be obtained by adding a particular solution gained by the forcing function [12] to the solution of a homogeneous equation.
If there is a constraint force (load on joint) on the node, we have a problem with non-homogeneous boundary conditions that can be converted to an equivalent problem with homogeneous boundary conditions [21, p. 57; 22, p. 43].
Let us regard a set of transfer equations (similar to Eq. ( 1)) interconnected through the boundary conditions to a complete system: where the vector Φ Φ Φ components Φ k (k = 1, 2, ..., N) are unknown state vectors of element ends and support reactions, where state vectors Z A and Z L components are Φ i (i = 1, 2, ..., n) and dynamic support reactions are C j ≡ Φ n+ j (j = 1, 2, ..., m), n + m = N.The term spA (ω) is an augmented transfer matrix.The righthand side Z (global loading vector) of the equation system contains element loading vectors Z p and nodal loads.Boundary conditions play an important role for transfer equations [23, p. 1115].The beam elements in Eq. ( 2) are interconnected through boundary conditions1 [8, pp. 34-48]: • compatibility equations of the displacements at nodes (geometric/essential boundary condition); • joint equilibrium equations at nodes (natural boundary condition); • side conditions (for bending moment, axial and shear force hinges); • support conditions (restrictions on support displacements).
Here, by the improved or modified transfer matrix method, unlike the transfer matrix method (TMM) [2; 25, p. 236], transfer matrices are not multiplied to find the initial parameters (state vectors) [8, p. 49].Hence the roundoff errors generated by multiplying transfer arrays are avoided.We will scale up (multiply) the displacements by the scaling multiplier.After solving the system of linear equations, we scale down (unscale) the initial parameter vectors of the elements dividing each of the displacements found by the scaling multiplier.
In a modal analysis, for the system of equations (2) the load vector is set to zero [26, eq. ( 31)]: For the nontrivial solution Φ Φ Φ i of the homogeneous system (3), we will choose a free variable in accordance with the natural frequency ω i : Here ω i denotes different natural (or characteristic, or normal) frequencies that are found numerically by the bisection method.These values are conventionally arranged in sequence from smallest to largest (ω 1 < ω 2 < ...ω n ).
For all the frequencies picked out from Eq. ( 4), the given boundary conditions and transfer equations are met.
The mode shapes are calculated according to Eq. ( 3), where the column of free variables is shifted to the right-hand side, and the system of equations obtained is solved with the least-squares method.After finding the initial parameters, we compile the mode shapes.
Free vibrations of beams with boundary and initial conditions are dealt with in [27].
To sustain vibration, energy must be supplied or tranferred out.Two phenomena, resonance and dynamic vibration absorption, are of importance in steady-state forced vibrations [28].On the frequency axis of steady-state forced vibration response curves, singular points -star and saddle points -lie [28; 29, p. 143].
In the saddle points of a LTP system, dynamic vibration absorption occurs.

STEADY-STATE FORCED VIBRATION OF BEAMS
Consider the free body diagram of a differential element of the beam in Fig. 1, where ẅ = ∂ 2 w/∂t 2 .We will apply d'Alembert's principle to extract the partial differential equation for transverse vibration of an elastic beam: where m = ρA with ρ as the mass density and A as the cross sectional area of the beam element.Now the Euler-Bernoulli hypotheses are used.The constitutive law for the beam relating the displacements w (x) (curvature d 2 w/dx 2 ) and the bending moment M y is where EI y is the bending stiffness of the beam.Combining Eqs ( 5), ( 6) and ( 7), we get (cf.[31, p. 1029]) Let us suppose that for the solution w (x,t), space and time given as separated functions (cf.[32, pp.9-10; 33, eq. ( 8); 34, p. 21; 27]): Fig. 1.Inertial force acting on a differential beam element.Here f (x) is a function of the independent variable x, and in the general case, θ is a complex frequency: where ω = Re (θ ) and ζ ω = Im (θ ) (see double imaginary characteristic roots on a complex plane [35, p. 29; 36; 37]).At steady-state vibration, Im (θ ) describes: 1) in time-domain, the absence of damping, and 2) in requency-domain, a phase-angle jump (PAJ) at natural frequencies (cf.[38]).The couplings between the frequencies of input and output signals can be taken into account [20,39,40; 41, pp. 13-14; 42].
Here the stability of a system is held together with its equilibrium state: a) at negative values of ζ ω < 0, the slightly disturbed equilibrium state remains stable; b) at positive values of ζ ω > 0, the slightly disturbed equilibrium state becomes unstable; c) at zero values of ζ ω = 0 (see sign principle [43, p. 444; 44]), the equilibrium state remains neutral/indifferent, the system is marginally stable.
We start counting the coordinate at x = x • = 0 and time at t = t • = 0 when the steady state frequency is same as the driving frequency (cf.differential equation being described as a system of first order differential equations [34, p. 29]).The natural exponential function is equal to 1 if exp (i (ϖ − θ )t = 0); it means that there are two possibilities: t = t • = 0 (cf.initial conditions in [27, p. 11]), or at 0 < t ≤ nT , where T is a period of vibration, the condition must be satisfied (cf.boundary conditions [27, p. 11]).
Let us start our investigation with a neutral/indifferent equilibrium (ζ ω = 0, see the sign principle).The zero values of ζ ω = Im (ω) = 0 divide the parameter space into regions by the stability.
If the area A and moment of inertia I y of the beam cross-section are constant, then using the assumption of Eq. ( 13), we get from the differential equation ( 12) a non-homogeneous 4th order differential equation to find the amplitudes of steady-state output response [20, p. 62]: Equation ( 14) matches [33, eqs (9), (12)].At steady-state forced vibration, for systems of periodically intermittent time [47] with (n + 1)T ≥ t ≥ nT , the frequencies ω n lying on the abscissa axis of the amplitude-frequency-plane are given as singular points with n = 1, 2, 3, ..., N and N → ∞ [28,49].The isolated singular points are star points, and the double singular points are saddle points.
The homogeneous differential equation below serves to find eigenvalues and eigenvectors: Here we use an auxiliary variable κ where κ represents the repeated roots of a characteristic (or frequency, or secular) equation of the linear differential equation for a beam: We have double real and double imaginary characteristic roots (κ 1,2 and κ 3,4 , respectively) for dynamical systems [36,37].Dimensionless eigenvalues where is the beam length.
The basic equation (2) of the EST method used in this paper may be considered as an improved transfer matrix method to find the state vectors, e.g., Z A and Z L in Eq. (1).Due to the normed fundamental set of solutions, the output parameters do not change the zero value of initial parameters at x = 0. Unlike the traditional transfer matrix method [2; 25, p. 236], here the transfer matrices are not multiplied to find the initial parameters.The novelty of this approach lies in the initial parameter vectors found by compiling sparse linear systems of equations incorporating transfer equations and boundary conditions (Eq. ( 2)) that are solved directly.Thus, the roundoff errors generated by multiplying transfer arrays are avoided.
First, we determine the state vectors Z A and Z L in Eq. ( 2) with the basic equations of the EST method [8, p. 49] that fit the solution of the homogeneous linear ordinary differential equation (15) (see [8, p. 33]).Further we calculate the state vector Z L (x) in Eq. ( 2) which is consistent with the non-homogeneous equation (14).The EST method makes use of the variation of parameters to solve problems of steady-state forced vibrations as well as statics of structural systems with interconnected elements.
In order to solve Eq. ( 37), we need to find the loading vector Z p .The frequency ω at singular points has a phase-angle jump Δω n associated with in-phase/out-of-phase behaviour [53, slide 36].The dimensionless frequency phase-angle jump Δλ k at singular points is also used.For sinusoidal response, a phase shift to the opposite phase is equal to π (out-of-phase) and a shift to the same phase is equal to 2π (in-phase).
At singular points, the amplitude f (x) sign changes into reverse (Fig. 10a) associated with the inphase/out-of-phase behaviour [53, slide 36].
The general solution f (x) of the non-homogeneous differential equation ( 14) can be expressed as a sum of the general solution f h (x) of the complementary equation ( 15) and the particular solution f e (x) of the non-homogeneous differential equation (14).
The fundamental set of solutions to the differential equation ( 16) has the form We norm the fundamental set of solutions ( 20) so that the Wronskian W (x) (normalized fundamental matrix) is the determinant of the identity matrix I 4×4 at x = 0.The normed fundamental set of solutions for the homogeneous differential equation is given below [9, p. 67; 33, eq. ( 15)]: The functions K i (κx) are also called Krylov-Duncan functions [57, p. 192 1  and 2).
Table 1.A cyclic order of Krylov-Duncan functions derivatives Derivatives First Second Third Fourth Table 2. Initial values: derivatives of normed fundamental solutions at x • = 0 Derivatives First Second Third Fourth Krylov Duncan functions and their derivatives satisfy permutations [57, p. 192; 17, p. 247] (see Tables 1  nd 2).Table 1.A cyclic order of Krylov-Duncan functions derivatives Derivatives First Second Third Fourth Table 2. Initial values: derivatives of normed fundamental solutions at x • = 0 Derivatives First Second Third Fourth The particular solution of Eq. ( 14), with zero initial value (zero state response), is obtained using the convolution integral [58, p. 156; 59, p. 279; 9, p. 92] or, to be more precise, Here G n (x, ξ ) is the normed fundamental set of solutions to the associated homogeneous differential equation, given by Eqs ( 21) -( 24): The load function g n (ξ ) from Eq. ( 25) is described with (cf.Eq. ( 14)).We get the following particular solutions, where the relationship λ = κ is taken into account and x a = a.For q z f 4 e (x) = q z E I y This particular solution compares well with solutions in [57, p. 197] and [2, p. 143]. For The present particular solution compares well with the solution in [14, p. 120].
To create the loading vector Z p = Z q + Z F + Z M of the transfer equations, we use the particular solutions (31), (32) and (33).
Let us present the transfer equations for vibration of a Euler-Bernoulli beam (sign convention 2 is used): Fig. 3. Indices of state vector components for a beam element.
Components of the state vector are displacements, rotations, shear forces and bending moments at the ends of the element shown in Fig. 3.
Here, to find the initial parameter vectors Z A , we improve or modify the transfer matrix method.In the traditional transfer matrix method, to find the initial parameters (state vectors) the transfer matrices are multiplied.To avoid the roundoff errors generated by multiplying transfer arrays, we will compile sparse linear systems of equations containing transfer equations and boundary conditions.With these sparse equations, designated by us as the basic equations of the EST method, we find the initial parameter vector Z A .The sparse equations ( 2) can be expressed as the basic equations of the EST method for a beam: hence Here Z p is the loading vector, and UI 4×8 is the augmented transfer matrix (U 4×4 | −I 4×4 ): where is the length of the beam element, and i • is the scaling multiplier for displacements and rotations (e.g., i 0 = 1.0, i 0 = EI basic / basic ).In system (2), the first four equations represent the basic equation ( 39) of the EST method, the rest being boundary conditions.
The boundary conditions of a cantilever beam: Equated to zero, the determinant of the coefficient matrix of equations (39) (cf.Eq. ( 4)) -the boundary conditions of which are expressed by Eq. ( 43) (see [9, p. 72]) -allows us to form the frequency (or characteristic, or secular) equation for a cantilever beam.After reducing the determinant and substituting κ = λ , we obtain the frequency equation (cf.[14, p. 110; 60, p. 527]) The boundary conditions of a fixed-fixed beam: Equated to zero, the determinant of the coefficient matrix of equation ( 39) (cf.Eq. ( 4)) -the boundary conditions of which are expressed by Eq. ( 45) (see [9, p. 73]) -makes it possible to write the frequency (or characteristic, or secular) equation for a fixed-fixed beam.We obtain the frequency equation after reducing the determinant and substituting κ = λ (cf.[14, p. 109; 60, p. 527]): The boundary conditions of a simply supported beam (with a pin connection on one end and a roller support on the other): Equated to zero, the determinant of the coefficient matrix of equation (39) (cf.Eq. ( 4)) -the boundary conditions of which are expressed by Eq. ( 47) (see [9, p. 74]) -makes it possible to write the frequency (or characteristic, or secular) equation for a simply supported beam.We obtain the frequency equation after reducing the determinant and substituting κ = λ (cf.[14, p. 108; 60, p. 527]): shλ sin λ = 0 (if λ = 0, then sin λ = 0), ( 48) The boundary conditions of a propped cantilever beam (fixed on one end, the free end resting on a roller support): Equated to zero, the determinant of the coefficient matrix of equations (39) (cf.Eq. ( 4)) -the boundary conditions of which are expressed by Eq. ( 49) (see [9,  Example 2.1 (steady-state forced vibration of a fixed-fixed beam).Compose the steady-state frequency response curves of a fixed-fixed beam.Find the amplitudes of displacements, angles of rotation, bending moments and shear forces of the beam in Fig. 4 under the excitation force F b e iϖt .
The fixed-fixed beam is of length = 0.229 m; the distance a from the left end of the beam to point b is 0.1145 m; the cross-sectional height is 0.79 mm and width 1.27 cm; elastic or Young's modulus E = 205 GN/m 2 ; shear modulus G = 76.53228GN/m 2 ; mass density ρ = 7.870×10 3 kg/m 3 ; A red = A/1.2(A is the cross-sectional area).The forcing amplitude F b = 0.889644 N with frequencies ϖ b = 376.99s −1 ∼ = f b = 60 Hz and ϖ b = 628.32s −1 ∼ = f b = 100 Hz (cf.[62]).
The relationship between the angular frequency ω and dimensionless frequency λ is given by the formula (cf.[9, p. 77]) The system of EST-method equations (40) (cf.Eq. ( 2)) is where Z is the vector of unknowns: The state vector input Z a and output Z c components are displacements and forces at the ends of the beam ac in Fig. 5:  The components of the state vector Z (i) (i = 1, 2, 3, .., 8) (cf.Φ (i) in Eq. ( 2)) in index notation are shown in Fig. 5.
In system (52), the first four equations represent the basic equation ( 39) of the EST method, the rest are boundary conditions: For nontrivial solutions of the homogeneous system (52), natural frequencies ω i of the beam are found.Figure 6 illustrates the dependence of the determinant of the coefficient matrix of Eqs ( 52 With this formula, we convert the natural eigenvalues ω i to dimensionless eigenvalues λ i : To apply a nodal load, we divide the beam into two elements (Fig. 7).The natural frequencies of the fixed-fixed beam with two elements are equal to the frequencies found with one element (Fig. 5).Dependence of the determinant on angular frequency Fig. 6.Natural frequencies of fixed-fixed beam.The components of the state vector Z (i) (i = 1, 2, 3, ..., 16) (cf.Φ (i) in Eq. ( 2)) in index notation are shown in Fig. 7.In the system of EST-method equations (52), the first eight equations represent the basic equation ( 39) of the method.The following four equations represent compatibility of the displacements and joint equilibrium at node b: Now we apply the boundary conditions as restrictions on support displacements: The sparsity pattern of the coefficient matrix spA of Eqs (52) of the fixed-fixed beam with two elements (Fig. 7) is shown in Fig. 8.
For the nontrivial solution of the homogeneous system (52), we will choose a free variable in accordance with the natural frequency ω i (see program script excerpt 2.1).
The first mode shapes are found with a GNU Octave script.We divide the beam of length into two elements of lengths 1 and 2 (see program script excerpt 2.1).The fifth column is shifted to the right-hand side and the equations are solved with the least-squares method.After finding the initial parameters, we compile the mode shapes.Steady-state forced vibration can be represented by response curves [28] with the eigenvalues (saddle and star points) lying on the abscissa axis and the transversal displacement and elastic strain energy amplitudes on the ordinate axis [29, p. 143].
where U is the elastic strain energy; GA red is the shear stiffness of a beam, and EI y is the flexural rigidity of a beam.The elastic energy U of a fixed-fixed beam calculated with Simpson's rule: where Δ = /16 and f 1 To construct response curves, we use the particular solution Z p = Z F (Eq. ( 35)).After finding the initial parameters, the displacements and forces are found with the transfer equations (37).
The star points are related to resonance [54, p. 521], and that is why the amplitude increases steadily, approaching infinity.
On the frequency axis (ω-axis) (Fig. 10), the location of the natural frequencies ω k matches that of the singular points lying on the same axis.Here the natural frequencies ω k with k = 1, 3, 5, ... are double real auxiliary roots (Eq.( 18)).Calculations show that near the singular points of odd numbers, the equilibrium state is changing rapidly and signs of the amplitudes become reverse.The singular points are star points; the equilibrium state is unstable.The elastic energy changes rapidly towards the maximum (see Fig. 10b).
Also located on the frequency axis (Fig. 10) are the natural frequencies ω k (k = 2, 4, 6, ...), which are double imaginary auxiliary roots (Eq.( 18)).Calculations show that near the singular points of even numbers, the equilibrium state is changing slowly and signs of the amplitudes become reverse.The singular points are saddle points.The equilibrium state is marginally stable (neutral/indifferent) (see Fig. 10b).
The neighbourhood of double imaginary characteristic roots in saddle points needs further investigation.In saddle points, the pumping frequencies periodically shift "harmonics" [41, p. 18].It is important to perceive if there exist regions of "undetermined" states of equilibrium with the so-called "parasitic selfexcitation" (cf.gain of the pumping frequencies in saddle points).This phenomenon appears and disappears, often caused by tiny variations of the parameters in these regions of equilibrium states [64, p. 658].
The saddle points are associated with dynamic vibration absorption described in [65] for LTI systems.
In Figs 11 and 12, the deflection, slope, shear force and bending moment amplitudes at saddle point frequencies (the LTP system theory) of a fixed-fixed beam under steady-state forced vibrations are shown.

CONCLUSIONS
A modified transfer matrix method has been developed for solving the systems of first order differential equations for vibration with a set of initial values and boundary conditions.Steady-state frequency response curves of a beam are composed with singular points (star and saddle points) lying on the frequency axis of the response curves.At these points, the frequencies coincide with these determined by the homogeneous differential equation.Star points represent resonance frequencies.In the case of undamped vibration, dynamic vibration absorption takes place at saddle points.

ACKNOWLEDGEMENT
The publication costs of this article were covered by the Estonian Academy of Sciences.

Fig. 3 .
Fig. 3. Indices of state vector components for a beam element.

Compatibility equations of displacements 9 Fig. 8 .
Fig. 8. Sparsity pattern of matrix spA of the fixed-fixed beam.

Figure 9
Figure9depicts the four displacement mode shapes of the fixed-fixed beam.Steady-state forced vibration can be represented by response curves[28] with the eigenvalues (saddle and star points) lying on the abscissa axis and the transversal displacement and elastic strain energy amplitudes on the ordinate axis[29, p. 143].
Fig. 10.Steady-state frequency response curves of fixed-fixed beam.
Fig. 10.Steady-state frequency response curves of fixed-fixed beam.

Fig. 12 .
Fig. 12. Fixedfixed beam slope and bending moment amplitudes at saddle point frequencies.(Continued on the next page.)

Fig. 12 .
Fig. 12. Fixed-fixed beam slope and bending moment amplitudes at saddle point frequencies.

Figure 13 Fig. 13 .
Figure 13  illustrates the amplitudes of displacements and forces of the fixed-fixed beam at 60 Hz frequency.

Figure 14
Figure 14  illustrates the amplitudes of displacements and forces of the fixed-fixed beam at 100 Hz frequency.

Table 1 .
The cyclic order of Krylov-Duncan functions derivatives