A Study on the Increase of Numerical Stability and Accuracy of the Transfer Matrix Method

Problem Statement: The transfer matrix method is a very useful tool for the static and dynamic analysis of structures. There are a number of issues though that worsens the numerical stability and the accuracy of this method. Approach: This study proposed a simple technique that can be used to handle these numerical difficulties and overcome the problems they give. Its main idea was to apply the method twice starting from two far ends of the structure. Results: An example from the calculation of the sensitivity function between two points in a dynamic system is presented. The results presented the very big potential of the proposed method. The improvement of the stability was clear in the graphs of the results. An initial study on the limitations of the proposed technique was also briefly given, together with some initial thoughts on how to overcome them. Finally, an idea of a possible use of the method for the maintenance studies of a high-speed rotor is presented, showing the very big variety of applications this methodology can be applied into. Conclusions: The proposed technique was very simple and effective, and hence it should be applied whenever the transfer matrix method was used.


INTRODUCTION
The transfer matrix method is a very efficient method for the solution of the differential equations describing the dynamic behavior of engineering systems. According to [1] , it can be used for the analysis of continuous beams, plates and shells, turbinegenerator shafts, crank shafts, and many other structures. It is a method widely accepted in practice and many books have been written that adapt it to different applications, as for example rotor dynamics [2] , or panels and aircraft structures [3] .
Despite the fact that this method can be used in such a big variety of applications, it has a very big drawback that finally limits its use. When the frequency of the examined system is high, there appear a number of numerical instabilities, leading to inaccurate results or even complete inability of the method to reach a solution. This problem was pointed out since the initial steps of the method development [1] and the limitations imposed were highlighted by almost everyone who dealt with it thereafter.
The aim of this study is to propose a new idea for solving the inherent numerical problems of the transfer matrix method. A detailed theoretical analysis of the proposed idea is given by using a representative example, the analytical solution of which is known and can be used for the presentation of the very promising results and potential of the proposed method.

MATERIALS AND METHODS
The idea behind the transfer matrix method is the split of the modeled structure into a number of elements, each being described by a matrix quantity, which relates the state vector at one side of the element to the state vector at the other. The state vector describes the state of the internal forces and displacements at the points where the divisions of the system have been made. The matrices connecting the state vectors can be either so called point matrices that describe discontinuities along the structure (e.g. inertia forces, shape discontinuities, reactions), or field matrices that describe the elastic behavior of the structure. More information on the formulation and theory of the method are given in [1] and are also briefly presented in the following paragraphs. Figure 1 presents schematically all the notation used in the transfer matrix method: Z is the state vector, V is the vertical load, M is the moment, is the rotational frequency, Fig. 1: Transfer matrix notation m is the mass, y is the deflection, is the density of the shaft, E is the Young's modulus, L (or l) is the length and I is the moment of inertia. Subscript i refers to the "station" where we define a state vector and superscripts L and R refer to the left and right of these "stations".
If we consider the problem to be one dimensional, the state vector for this case will be: The differential equation describing the dynamic behavior of the "elastic part" between two "stations" is given by: = The area of the cross section K = A factor dependant of the cross section shape that takes into account the effect of shear on the deflection.
Next we use the following notation: Hence, the differential equation can be written: Solving this equation and writing the result in matrix form and using the following equations: we get: Where, F is the field matrix, given by Eq. 12. In this equation it is: The combination of hyperbolic and basic trigonometric functions in the above equations is one reason for the numerical instabilities of this method.
If at a point along the "elastic part" there is a discontinuity of any of the components of the state vector due to a concentrated mass, a gyroscopic load or a support of the structure, then a point matrix P is introduced to account for this discontinuity. For example, if we have a concentrated mass, it introduces a discontinuity in shear forces as follows: Or in a matrix form: In a similar way, the loading on a station can be added. Multiplying the above defined transfer matrices one after the other as we move along the structure and applying the boundary conditions, we can calculate the state vector at any position we want, thus calculating all the dynamic characteristics. A more detailed explanation of the whole procedure is out of the scope of this paper, but more details are given in [1,3] In order to demonstrate the method and also describe the present idea for overcoming its numerical difficulties, we will consider the example of a flat panel in 2D that is under a forcing pressure profile of given spatial and temporal characteristics: In Fig. 2 Zn(0) and Zn(N) are the state vectors describing the dynamic characteristics of the two ends of the panel along y. We have also assumed that the description of the panel's deflection along x follows a predefined function (e.g., Fourier series expansion). According to [3] , the sensitivity function of the structure represents the sensitivity of it, as manifested at a point r to a harmonic excitation having a frequency and a sinusoidal spacewise variation defined by a wave number vector k. The closed form equation giving this sensitivity function for the studied example will be: In Eq. 22 H is the frequency response function, X and Y are the functions giving the deformed shape of the structure for given boundary conditions and P mn is the shape of the forcing pressure. For simplicity, we consider the panel at the present example to be simply supported in all its four edges and hence we have: Let's also consider the forcing pressure to have a sinusoidal spatial variation with wavenumbers: This will result in the three-dimensional spatial distribution of forcing pressure shown in Fig. 3.
If we were to use the transfer matrices method to calculate the sensitivity function, we would have to use just a field matrix to describe the problem, since there are no discontinuities in loading or shape in this case. Using also as loading in this case a surface grid of weighted point loads, where the weighting shape function is a sinusoidal function with the known k x and k y , the result of the transfer matrix method will be the sensitivity function.
In the literature (e.g. [1] ), various numerical difficulties of the application of the transfer matrix method are pointed out. Among them, there is the application of it at high frequencies, the use of intermediate stiff supports in the structure, or the inclusion of structures of big length. Pestel, E.C. [1] proposed a method called "Delta-matrix method" that could deal with some of the problems and improve the accuracy of calculations of critical frequencies, but it could not be used for mode shape or other calculations that involve displacement characteristics. In the proposed method, the transfer matrix method is applied twice, starting from both the two ends of the structure, Zn(0) and Zn(N), and progressing towards the right and left respectively. The results of each case are used up to the middle of the structure, leading to the usage of smaller lengths in transfer matrix calculations. That way, the numerical instabilities that are present at points far from the initial state vector are diminished. A more detailed description of the procedure and its results are given next.

RESULTS
Using the unmodified transfer matrix methodology and the closed form solution equation, we calculated the sensitivity function at 1500Hz for a panel of a = 0.768m, b = 0.328m and thickness of 1.6mm. The material of the panel is assumed to be aluminum with a density of 2700 Kg m− 3 , Young's modulus 70GPa and Poisson's ratio 0.33. The results are shown in the Fig. 4. More specifically, Fig. 4 shows the comparison of the analytic results (solid line) with the results of the transfer matrix method without the use of the methodology of this paper (doted line).
It is obvious that starting from y=0 there is a very good agreement between the two results, but as we move to bigger y, the numerical stability problems of the transfer matrix method become clear. Based on the form and the terms of the field matrix, we come to the conclusion that as y is getting bigger, so do the elements of the field matrix and hence we end up in a situation where the successive multiplications of matrices or determinant calculations lead to classical round off problems, like for example additions of numbers some orders of magnitude different, or subtraction of almost equal numbers.
The present idea in order to overcome this problem was to split the solution of the problem to two different problems. One deals with the calculations up to It is clear from the Fig. 5 that using this approach, the numerical stability improves considerably. Nevertheless, we should point out that in case the size of the structure is very big, it is expected to have some instabilities or discontinuities in the middle of it, where the results of the two solutions merge. This effect can be somehow alleviated by averaging the two results, but anyway results should be interpreted carefully in such cases.

DISCUSSION
The transfer matrices method is a very efficient way of solving dynamic analysis problems of several structures. Its inherent numerical accuracy problems have lead to its limited application in high frequency analyses and analysis of big structures. The method proposed in this paper alleviates these numerical difficulties and extend the application range of transfer matrices. It can also be the basis for further research on the subject, probably targeted to specific applications.
In order to demonstrate the range of different applications the proposed method can be applied, it was used by the author in the development of a computer program for the calculation of the frequency response function of rotors. The aim of this was to be able to calculate the attenuation of the signal measured by a vibration transducer at one point of the rotor due to a fault present at another point. This result is very useful during the assessment of the operational integrity of rotating systems, given that it is possible to derive conclusions regarding measurements of faulty components that are not close to the measuring point. Since some of the faults appear at high frequencies or have high characteristic frequencies themselves, it would not be possible to use transfer matrices in this problem without the application of the method proposed in this study.

CONCLUSION
This study proposes a very simple but effective method for increasing the numerical stability and accuracy of the transfer matrix method. It is based on the simultaneous solution of the problem in two directions and a combination of the two results. The effectiveness of the method is presented using as an example the calculation of the sensitivity function of a flat panel loaded by a given pressure profile, the analytical solution of which is also calculated.
The proposed method can be used in all transfer matrices applications, as for example in the calculation of the dynamic characteristics of rotors (critical speeds, mode shapes, response to unbalance), or for the calculation of frequency response function of such systems.