Modiﬁed transfer matrix method for steady-state forced vibration: a system of bar elements

. The Elements by a System of Transfer (EST) method offers exact solutions for various vibration problems of trusses, beams and frames. The method can be regarded as an improved or modiﬁed transfer matrix method where the roundoff errors generated by multiplying transfer arrays are avoided. It is assumed that in a steady state a bar/beam will vibrate with the circular frequency of a harmonic excitation force. The universal equation of elastic displacement (2nd/4th order differential equation) is described as a system of ﬁrst order differential equations in matrix form. For the differential equations the compatibility conditions of a bar/beam element displacements at joint serve as essential boundary conditions. As the natural boundary conditions at joints, the equilibrium equations of elastic forces of bar/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. The magniﬁcation factor at the excitation frequency is determined.


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 [1]. Two phenomena, resonance and dynamic vibration absorption, have been of great interest [2,3]. According to the work-energy theorem valid in structural analysis the sum of the work done by inertial, internal and external forces is zero: where W T is the work done by inertial forces; W i -work done by internal forces; W e -work done by external forces; W f -work done by active forces, e.g. concentrated loads, uniformly distributed loads; W b -work done by constraint forces, e.g. support reactions, internal reactions (a contact force acts at the point of contact between two objects [4]). The work-energy theorem for a frame element [5, p. 60; 6, pp. 683-687; 7, p. 221] can be written as: In the equation above, we consider two load states (see Basic variational principles [8, p.32]) associated with respective deformations and displacements: N x , Q z , M y -internal axial force, shear force, and bending moment of the first load state; λ ,β z ,ψ y -axial, shear, and bending deformations of the second load state; N x | b a , Q y | b a , M y | b a -axial force, shear force and bending moment of the first load state at boundaries a and b; u| b a ,ŵ| b a ,φ y | b a -longitudinal and transverse displacements, and the rotation of the cross section of the second load state at boundaries a and b; q x (x), q z (x) -distributed loads of the first load state; u(x),ŵ(x) -longitudinal and transverse displacements of the second load state; F xi , F zi -force components of the first load state, applied at the point i in xand z-directions, respectively; u i ,ŵ i -longitudinal and transverse displacements of the point i of the second load state; u(x),ẅ(x) -longitudinal and transverse accelerations of the first load state; ρ -mass density; A -cross sectional area of the bar/beam element; D -energy dissipation (entropy production D ≥ 0 [9, p. 8; 10]. In computational structural mechanics, the state-space representation of mechanical systems can be seen as transfer matrix method [11,12,6,7,13]: where Z A , Z L designate the components of the state vectors (displacements, internal forces at the beginning x = 0 and the end of the element x = , respectively); Z p is the element loading vector; U denotes the transfer matrix.
Elastic forces cause elastic axial, shear, and bending deformations:λ =N x /EA,β z =Q z /GA red ,ψ y = M y /EI y . So the work of elastic forces W i can be written as wherê N x ,Q z ,M y -internal axial force, shear force, and bending moment of the second load state; EA -axial cross-sectional stiffness of bar/beam element; GA red -shear cross-sectional stiffness of beam element; EI y -flexural cross-sectional stiffness of beam element.
An external force can be added as an element load (described with a forcing function) or as a nodal load on joints. The general solution to an ordinary differential equation can be obtained by adding to the solution of a homogeneous equation a particular solution obtained with the forcing function [15]. The forcing function for linear time periodic (LTP) system [16] is dealt with in [17, p. 121; 18; 19; 20, p. 248; 13, pp. 26, 93, 94].
With a constraint force (load on joint) on the node, we have a problem with nonhomogeneous boundary conditions that can be converted to an equivalent problem with homogeneous boundary conditions [21, p. 57; 22, p. 43].
Besides, the transfer equations (4) and (6) contain the following boundary conditions: • compatibility equations of the displacements at nodes (geometric/essential boundary conditions); • joint equilibrium equations at nodes (natural boundary conditions); • side conditions (for bending moment, axial and shear force hinges); • support conditions (restrictions on support displacements).
The assembled member-end displacements compatibility conditions at joint node [7, pp. 34, 36] (cf. force (flexibility) method) and the assembled member-end forces equilibrium equations at joint node [7, pp. 40, 41] (cf. displacement (stiffness) method) are included in the nodal transfer equations for natural vibration analysis of tree system [11].
Let us present the system of equations (4) with boundary conditions as follows: where the vector Φ Φ Φ components Φ 1 , Φ 2 , ..., Φ N contain unknown state vectors of element ends (Z A and Z L ) Φ i (i = 1, 2, ..., n), and unknown support reactions (C j ) Φ n+ j (j = 1, 2, ..., m), n + m = N. The spA (ω) is the augmented transfer matrix. The right-hand sideZ (global loading vector) of the equation system contains element loading vectors Z p and nodal loads.
Here, by the improved or modified transfer matrix method, unlike the transfer matrix method [1; 23, p. 236], transfer matrices are not multiplied to find the initial parameters (state vectors). 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 (6) the load vector is set equal to zero [24,Eq. (31)]: For the nontrivial solution Φ Φ Φ i of the homogeneous system (7) we will choose a free variable in accordance with the natural frequency ω i : Here ω i denotes natural (or characteristic, or normal) frequencies that are found numerically by the bisection method. These frequencies are conventionally arranged in sequence from smallest to largest (ω 1 < ω 2 < . . . ω n ).
With all the frequencies that satisfy Eq. (8), the given boundary conditions and transfer equations are met.
The mode shapes are calculated according to Eq. (7). Here 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, the mode shapes are compiled.
To check their correctness, the natural frequencies and mode shapes computed by the EST method were compared with those found by other methods: 1. The natural frequencies of a five-span continuous beam in [13, p. 88 fig. 3.7]. 4. The natural frequencies for a truss structure in [24] gained by the EST method have been compared with true frequencies of natural vibration obtained through constraint equations by Ramsay [28]. The results coincide with the exact ones of the problem in [28, fig. 4]. 5. The displacements and bending moments obtained in [13, p. 95] produced by forced vibration (at frequencies smaller and larger than the first resonance frequency) of the simply supported Euler-Bernoulli beam were compared with those in [17, p. 116] and found to be in close agreement.
To sustain vibration, the energy must be supplied or transferred out. For steady-state forced vibration, resonance and dynamic vibration absorption [29] are of importance. The singular points lying on the frequency axis of steady-state forced vibration response curves are star points and saddle points [29; 30, p. 143].
The total response of the LTP state space model to a complex sinusoid is the sum of the homogeneous and forced responses [31, p. 62]. The resonant frequency identified in [16, p. 4; 32, pp. 2, 3; 33].
At the saddle points of a LTP system, dynamic vibration absorption occurs.

STEADY-STATE FORCED VIBRATIONS OF A BAR
Consider the free body diagram of a differential element of the bar in Fig. 1: . .  We will apply d'Alembert's principle to derive the partial differential equation for longitudinal vibration of the elastic bar: where m = ρA. The constitutive law for the bar relates the displacements u (x) (point strain du/dx) and the normal force N: Let us suppose that the loading is given as p x (x,t) = n x sin ϖt (n x is taken as a constant), then Eq. (11) takes the form The solution u (x,t) is assumed to be in the form of separated functions: where f (x) is only a function of the coordinate x and sin ωt is only a function of the time variable t. Substituting Eq. (13) into Eq. (12) we obtain Let us take a driving frequency ϖ = ω (ϖ = ω n ; ϖ − = ω n − ε, ϖ + = ϖ − + ϖ n ; ω n denotes the natural frequency and ϖ n marks a phase-angle jump (PAJ)). The frequencies ω n at steady-state forced vibration in the amplitude-frequency-plane lying on the abscissa axis are singular (star and saddle) points n = 1, 2, 3, ..., N, N → ∞ [29]. The star points are isolated singular points, and the saddle points are double singular points. Now we get a nonhomogeneous differential equation to find the amplitudes of steady-state output response [31, p. 62]: The homogeneous differential equation below is for finding eigenvalues and eigenvectors The differential equation (16) has a characteristic (or frequency, or secular) equation with repeated roots: here κ is wave number and ω n denotes the natural frequency.
Dimensionless eigenvalue is also used: where is the bar length.
The basic equations (36) of the EST method used in this paper may be considered as an improved transfer matrix method to find the state vectors, eg., Z A and Z L in Eq. (36). Due to the normed fundamental set of solutions (Eq. (21)), the output parameters do not change the zero value of initial parameters at x = 0. Unlike the traditional transfer matrix method [1; 26, 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 (Eqs (6), (36)) 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 Eqs (6) and (36) with the basic equations of the EST method [7, p. 49] that fit the solution of the homogeneous linear ordinary differential equation (16) (see [7, p. 33]). Further we calculate the state vector Z L (x) in Eq. (30) of the nonhomogeneous equation (15). 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.
To solve Eq. (30) we must find the loading vector Z p . The driving frequency ϖ at singular points has the phase-angle jump (PAJ) ϖ k associated with the in-phase/out-of-phase behaviour [39, slide 36]. The dimensionless driving frequency PAJλ k at singular points is also used. For sinusoidal signals a phase shift (PAJ) to the opposite phase is equal to π (out-of-phase) and to the same phase is equal to 2π (in-phase).
A singular point is often associated with a sudden change in the system, solutions in singular points are unstable. In case of undamped harmonic loading, at driving frequency, the amplitudes f (x) at singular points will reach infinity (quality factor Q = ∞), and a phase-angle jump occurs [40, p. 534; 41]. At star points, amplitude changes are significantly larger than at saddle points, where amplitude changes should be determined with low threshold or can be labelled as 'positive' versus 'negative' (see Fig. 8a) [41]. A saddle point phase portrait is shown in [42,p. 181 (LMC 199), fig. 4.7 (c)]. For steady-state forced vibration loading, dimensionless driving frequency phase-angle jumpsλ k appear: • at star points (k = 1, 3, 5, ... odd numbers); • at saddle points (k = 2, 4, 6, ... even numbers).
At singular points, the amplitude f (x) sign changes into reverse (Fig. 8a) associated with the inphase/out-of-phase behaviour [39, slide 36].
The general solution f (x) of the nonhomogeneous differential equation (15) can be expressed as a sum of the general solution f h (x) of the complementary equation (16) and the particular solution f e (x) of the nonhomogeneous differential equation 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 2×2 at x = 0. The normed fundamental set of solutions for the homogeneous differential equation (zero input response) f h (x) is [13, p. 18]: The particular solution of Eq. (15) with zero initial value (zero state response) is obtained using the convolution integral [43, p. 156; 44, p. 279]: or, to be more precise, where G n (x, ξ ) is a normed fundamental set of solutions to the associated homogeneous differential equation (21): In Eq. (23), the load function g n (ξ ) is described: The convolution integrals produce [13, p. 22], where x a = a: the particular solution for F x the particular solution for n x Here it must be taken into account that λ = κ and x − a + is a discontinuity or singularity function (note the use of Macaulay brackets) x − a n where H (x − a) is the Heaviside function. The equation below represents the transfer equations for longitudinal vibration of a bar (sign convention 2 is used): Here, Z L (x) and Z A are the vectors of displacements and forces at the point with x-coordinate (the state vector of output end L = x) and at the beginning of the element x = 0 (the state vector of input end), respectively, Z p is the loading vector of the bar element.
To create the loading vector Z p = Z n + Z F of the transfer equations (30), we use the particular solutions and (27) and (28): The transfer equations (30) of a bar may be expressed as the basic equations of the EST method: or Here, UI 2×4 is an augmented transfer matrix (U 2×2 | −I 2×2 ): Z p is the loading vector; is the length of the bar element; i • is the scaling multiplier for displacements u i . Let us consider a bar of length and axial cross-sectional stiffness EA to determine the natural frequencies. The components of the state vector Z (i) ≡ Φ (i) (i = 1, 2, 3, 4) in index notation are shown in Fig. 2.
In the system (4), the first two equations represent the basic equation (35) of the EST method, the rest being boundary conditions. Boundary conditions for a fixed-free bar: The system of EST-method equations for a fixed-free bar in matrix form:   The determinant of the coefficient matrix of equations (40) equal to zero gives the frequency (or characteristic, or secular) equation where cos κ is equal to zero if Here λ n denotes the dimensionless eigenvalues of the fixed-free bar. Now the normal frequencies ω n of the bar can be calculated by the formula Boundary conditions for a fixed-fixed bar: The system of EST-method equations for a fixed-fixed bar in matrix form: The determinant of the coefficient matrix of equations (45) equal to zero gives the frequency (or characteristic, or secular) equation where sin κ is equal to zero if Here λ n denotes the dimensionless eigenvalues of the fixed-fixed bar. Now the normal frequencies ω n of the fixed-fixed bar can be calculated by the formula ω n = (n − 1) π E ρ , n = 2, 3, 4, . . . .
Example 2.1 (steady-state forced vibration of a fixed-fixed bar). Compose the steady-state frequency response curves of a fixed-fixed bar. Find the displacements and axial forces of the bar in Fig. 3. The system of EST-method equations (6) is where Z is the vector of unknowns: The components of the state vector Z (i) ≡ Φ (i) (i = 1, 2, ..., 8) in index notation are shown in Fig. 4.
In the system of EST-method equations (49), the first four equations represent the basic equation (35) of the method. The following two equations are displacement compatibility and joint equilibrium at node b (cf. relations between the status quantities in [45, p. 37]).   Now we apply the restrictions on support displacements (support conditions): Sparsity pattern of the coefficient matrix spA of the system of equations (49) is shown in Fig. 5. Figure 6 shows the dependence of the determinant of the coefficient matrix of equations (49) on angular frequency ω of the bar (0.6 · ω 1 = 9743. The relationship between dimensionless and natural eigenvalues (λ i and ω i , respectively) of the fixedfixed bar is given by the formula With this formula, we convert the natural eigenvalues ω i to dimensionless eigenvalues λ i : λ 1 = 3.141593,    The exact roots of frequency equation (dimensionless eigenvalues λ n and natural frequencies ω n ) of the fixed-free bar can be calculated using formulas (47) and (48). The natural frequencies of one-and two-element fixed-fixed bars proved to be equal. The first mode shapes are found with a GNU Octave script. We divide the bar of length into two elements of lengths 1 and 2 (see program script excerpt 2.1). The third column of the complementary equations of system (49) is shifted to the right-hand side and the equations are solved with the least-squares method.   (27)) can be represented by the response curves with the eigenvalues (saddle and star points) lying on the abscissa axis and the axial displacement/force amplitude on the ordinate axis [29; 30, p. 143] (Fig. 8a,b).
To construct response curves, we use the element load F x . The particular solution for Z p = Z F of the transfer equations (30) is that is, with equations (49), the basic equation (35) of the EST method is used (Z p = Z F ): We perform the calculations of axial displacement amplitudes and elastic strain energy starting from circular frequency ω = 0.1 up to ω = 1.13670×10 5 s −1 . It can be seen in Fig. 8a,b that star points abscissas on the frequency-axis are 16240, 48716, 81193, 113670, and saddle points abscissas are 32477, 64954, 97432. As the star points are related to resonance [40, p. 521], the amplitude becomes larger approaching infinity. In this case (friction not considered) we have nonequilibrium and are unable to make use of the theory of minimum energy dissipation rate [9, p. 8].
Integrating the first term in Eq.   For steady-state forced vibrations the extended framework of Hamilton's principle can be used [46]. The concepts of Hamilton's principle are: elastic strain energy (U ); kinetic energy (T ); Rayleigh's dissipation function (R); initial and boundary conditions.
The improved or modified transfer matrix method uses the strong form of boundary conditions (see p. 3). For compatible initial value, the normed fundamental set of solutions (21) is used.
At dynamic equilibrium (u = 0), the second term in equation (58) is For the real/actual work W (a) i of internal forces, the kinematically admissible displacements of the second load state are equal to the real/actual displacements of the first load state at time instants t 1 and t 2 , i.e., where U is the elastic strain energy of bar. From the last term of Eq. (58) we get the expression for kinetic energy Extended Hamilton's principle [46]: where V is the work done by external loads. Figure 8a,b show the steady-state frequency response curves of the fixed-fixed bar in the following intervals of frequency ω: On the frequency axis, the singular points 1, 3, 5, 7 are star points and 2, 4, 6 are saddle points.
The elastic energy U of the bar Eq. (60) is calculated with Simpson's rule:   In Fig.10, amplitudes of axial displacements and forces at 60% of the first natural frequency (ω = 0.6 ω 1 ) and the excitation force amplitude F b = 100 N are shown. The amplitudes are found with element and nodal loadings (cf. [45, p. 37]).
The calculation results of both loadings cases are identical.
The magnification factors for the dynamic displacement k s d and dynamic normal force k N d at the excitation frequency ω = 9743.2 s −1 are different: where u st is the static displacement and N st is the static normal force. The saddle points are associated with dynamic vibration absorption (described in [47; 48, p. 373] for the linear time invariant (LTI) systems). Figure 9 shows the dependence of the axial displacement and force amplitudes on the saddle point angular frequencies for fixed-fixed bar. The driving force at x = 0.5 m.
In Fig. 10, amplitudes of axial displacements and forces at 60% of the first natural frequency (ω = 0.6 ω 1 ) and the excitation force amplitude F b = 100 N are shown. The amplitudes are found with element and nodal loadings (cf. [45, p. 37]).
The calculation results of both loadings cases are identical. The magnification factors for the dynamic displacement k s d and dynamic normal force k N d at the excitation frequency ω = 9743.2 s −1 are different:

CONCLUSIONS
A modified transfer matrix method has been developed for solving systems of first order differential equations for vibration with a set of initial value and boundary conditions. The steady-state frequency response curves of a bar/beam are composed. Singular points (star and saddle points) lie on the frequency axis of the   response curves. At these points, the frequencies coincide with the frequencies determined by the homogeneous differential equation. Star points represent resonance frequencies. In the case of undamped harmonic vibration, dynamic vibration absorption takes place at the saddle points. When the force application point coincides with the node location of a standing wave, the standing wave occurs at the respective saddle point frequency.