Orthogonality of modes of structures when using the exact transcendental stiffness matrix method

This paper presents theory, physical insight and results for mode orthogonality of piecewise continuous structures, including both coincident and non-coincident natural frequencies. The structures are ones for which exact member equations have been obtained by solving the governing differential equations, e.g. as can be done for members of plane frames or prismatic plate assemblies. Such member equations are transcendental functions of the distributed member mass and the frequency. They are used to obtain a transcendental overall stiffness matrix for the structure, from which the natural frequencies are extracted by using the Wittrick-Williams algorithm, prior to using any existing method to find the modes which are examined from the orthogonality viewpoint in this paper. The natural frequencies and modes found are the exact values for the structure in the sense that the usual finite element method approximations are avoided.


Introduction
Structures such as plane or space frames and prismatic plate assemblies, e.g.stiffened panels, are usually composed of continuous members, i.e. beams or plates, which are uniform along their longitudinal direction.Such structures can be called 'piecewise contin-uous', because they are continuous between the joints, at which lumped masses or rotatory inertias may be present.
The natural frequencies and modes of vibration for such structures are often found by using the finite element (FE) method to set up the static overall stiffness matrix K s and the mass matrix M .The natural frequencies and modes are then found by solving the linear eigenvalue problem where ω is the circular frequency and the vector D contains the amplitudes of the N displacement components at the joints of the structure which all vary sinusoidally with time.Numerous methods exist for solving the linear eigenvalue problem of Eq. ( 1) to obtain the eigenvalues ω i (i = 1, 2, . ..) and the associated eigenvectors, i.e. modes, D i with a very high degree of reliability.Such methods can give the ω i and D i of Eq.
(1) to very high accuracy, but they are only approximations to the values for the real problem because of the approximations involved in obtaining the FE member stiffness and mass matrices which are needed to assemble K s and M .These FE approximations can often optionally be avoided when assembling K s by using the exact static member stiffness matrices, e.g. the well known slope-deflection equations when the members are uniform beams of a plane frame.However, the calculation of M in Eq. ( 1) still requires approximations to be made to obtain either lumped or consistent mass matrices.Hence the FE method gives approximations to ω i and D i which approach the correct values for the structure only if the structure is divided into numerous very small finite elements by introducing a very large number of additional joints, i.e. nodes, such that N → ∞.
The present paper deals with a less common, though still quite widely used, approach which gives values of ω i and D i which are potentially (depending upon the convergence methods chosen) exactly correct for the  1) would give if an infinite number of properly formulated finite elements were used, so that N → ∞, and if the computer worked to the impractically large number of significant figures needed to avoid any illconditioning of the results.This is achieved by solving the differential equations (DE's) of each member with the mass per unit length, µ, of the member included in the DE.Insertion of appropriate boundary conditions in the DE gives member equations which are transcendental functions of µ and of ω.The usual FE rules can be used to assemble the overall stiffness matrix K from the member equations.However, the elements of K are now transcendental functions of ω and of the values of µ for the component members.Hence, the eigenvalue problem is a transcendental one and the linear eigensolvers used to solve Eq. ( 1) are no longer applicable, and nor are the rules used to normalise the modes D i and confirm that they are mutually orthogonal.
Methods for finding the natural frequencies of this transcendental eigenvalue problem with complete certainty have been available since the publication of the Wittrick-Williams algorithm [4], which performs a function analogous to that of the Sturm sequence property of linear eigenvalue problems.Such methods can be made to converge on the ω i to very high accuracy.Methods have also been developed [1,3] for finding the associated modes from ω i .These methods have a sound theoretical basis and have proved to be reasonably robust, although D i can be significantly less accurate than ω i .These methods have been extended to include computation of the r modes corresponding to the coincident eigenvalues of multiplicity r, which are simply called 'coincident modes' below.
Very little has been published on the orthogonality of the modes of piecewise continuous structures obtained by the transcendental stiffness matrix method and in particular the authors are unaware of any treatment of the orthogonality of coincident modes.The present paper derives the orthonormalisation conditions for such modes and shows how they can be extended to orthonormalise any set of coincident modes which may be present.The orthonormalised set of coincident modes is not unique and a physical reason is given to show why this is so.
Illustrative examples are solved involving a simple rectangular plane frame with distributed mass used for the members and lumped masses or rotatory inertias optionally present at the nodes.

Derivation of conditions for normalisation and orthogonality of modes
The normalisation and orthogonality conditions for a piecewise continuous structure can be written, on normal Bernoulli-Euler beam type assumptions, as where: summations are over the n c members (i.e.beams for frame structures) and n j joints; integration is over the member length; k, i and j denote the k-th member and the i-th and j-th modes of the structure; µ is the member mass per unit length; u and w are, respectively, along and perpendicular to the member at distance x from one end; U , W and Θ are the joint displacements and rotation in global co-ordinates; M and I are the lumped mass and inertia (if present) at a joint; and δ ij is the Kronecker delta, i.e. δ ij = 0 for i = j and δ ij = 1 for i = j.Equation ( 2) can be derived as follows, from the two starting points that an FE model which uses N lumped masses to represent the distributed mass of each member gives the exact distributed mass solution (barring ill-conditioning) as N → ∞ and that for such an FE model, where D and M are, respectively, the displacement vector and mass matrix of the FE model.Now let N → ∞ and let all the degrees of freedom, D m , of nodes within the lengths of members precede the displacements D c at the joints, giving Eqs. ( 4) and (5), so that, because M is diagonal due to lumped masses being used, Eq. ( 3) can be partitioned to give Eqs.(6) and hence (7).
Because N → ∞, the first triple product of Eq. ( 7) clearly becomes the first summation of Eq. ( 2), since M m contains masses but not rotatory inertias.The second triple product is readily seen, by using Eq. ( 5), to give all the terms in the second summation of Eq.
(2).Hence Eq. ( 2) has been derived from Eq. ( 7).If rotatory inertias are included in the M m of Eq. ( 7), Eq. ( 2) is derived but with an additional summation over the n c members of I k θ ik θ jk dx, where θ and I are the rotation and rotatory inertia of a cross-section, e.g. when the members are beams they are Bernoulli-Euler ones but with the additional effects of rotatory inertia allowed for, i.e. they become Timoshenko members if shear deflection is also included.

Examples for checking orthonormalisation conditions
All correctly calculated non-coincident modes must be mutually orthogonal and so can be used to check the correctness of Eq. ( 2).Because the DE's have been solved the number of modes, n m , that can be computed is unlimited.This can alternatively be deduced as being because the order of Eq. ( 3) approaches infinity.However, in practice n m is usually chosen to be quite small.In the examples, the amplitudes of all 'raw' modes are normalised to satisfy Eq. ( 2) for i = j(i = 1, 2, . . ., n m ).The fact that the orthogonality conditions are satisfied could then be confirmed by showing that the left-hand side of Eq. ( 2) computes to zero for arbitrarily chosen modes i and j(i = j).However, because the computed modes usually have small relative absolute errors (defined in any sensible way) of ε i and ε j , respectively, the left-hand side of Eq. ( 2) usually differs slightly from zero, but should lie in the approximate range The P RT method [1] was used to find all the modes for the examples presented, i.e.KD c = P RT was solved very close to the i-th natural frequency ω i , where P RT was a randomly chosen force vector.Because ω i approximates the natural frequency, the response D c obtained closely approximates the raw mode required.However the mode accuracy ε i is usually poorer than the absolute accuracy ε ωi of the natural frequency.Typically [1,3] −100ε ωi ε i 100ε ωi .To avoid possible ill-conditioning, ε ωi should not be made too small and so the results given in this paper all used ε ωi = 10 −7 .Hence it should almost always be conservative to assume ε i = ε j = 10 −5 , so that Eq. (8) becomes −2 × 10 −5 LHS of Eq. (2) 2 × 10 −5 (9) The test structure is the two-bay single storey rigidly jointed plane frame with identical and inextensible members shown in Fig. 1.Extra joints were introduced at the beam and column mid-lengths to ensure accurate computation of modes such as those of Fig. 3 below, which would otherwise have D = 0. Computations used an adaptation of the simple plane frame program of Howson [2], in which inextensible members were modelled by setting the extensional rigidity to 10 10 EI/l 2 .A higher value was not used to avoid any danger of ill-conditioning and, for security, the absence of ill-conditioning was confirmed by recomputing key results with the extensional rigidity doubled.Results for several pairs of both non-coincident and coincident modes are given in Table 1.The table also shows the effect of adding equal lumped masses m at all three joints, where m is equal to the mass of each member of the frame.Additional results are also shown for when, instead, three lumped inertias of magnitude ml 2 were added at the three joints.
Figure 2 shows the first five non-coincident modes computed in the absence of lumped masses or inertias at the joints.Adding the lumped masses or inertias at the joints altered these modes, but the fifth and sixth modes remained coincident in both cases.(Strictly, they were only very nearly coincident, because the members were not quite inextensible ones).The three modes of Fig. 3 are clearly all possible choices and occur at exactly the clamped/clamped fundamental frequency of each of the five identical members.Note that they are obviously not orthogonal to each other.For example, for the first two modes shown the integrations of Eq. ( 2) are clearly zero for all the members except for the central column, so that they cannot sum to zero.

Orthonormalisation of two coincident modes
Suppose that the Wittrick-Williams algorithm [4] shows that two natural frequencies are coincident, or   appear so because they differ by less than the tolerance ε ωi required during convergence.Suppose also that a small alteration to the structure slightly separates these two natural frequencies.Their correct normalised modes will then be orthogonal to all the non-coincident modes and to each other.Let them be denoted as D i and D j , with j = i + 1.For convenience, and because the two have been shown via Eq.( 7) to be equivalent, Eq. ( 3) is now used to represent what happens when Eq. ( 2) is used.
Because D i and D j are orthogonal to each other, Eq.
(3) gives ) Now the P RT method performs calculations close to, but not at, ω i and ω j .As a consequence, instead of D i and D j , normalised combinations D * i and D * j of D i and D j are found, i.e. in the limit ) where A i , B i , A j and B j are unknown constants but are not unique, because they depend on ε ωi and on the will not usually be orthogonal to one another, because when substituted in the LHS of Eq. ( 3) they do not give zero, but instead give However, D * i , D * j , or any linear combination of them, can be chosen as one of the set of orthogonal modes and then D * i and D * j can be combined appropriately to find a second mode such that all modes, including the coincident ones, are orthogonal.For convenience, and without loss of generality because A i and B i can be chosen arbitrarily, the following proof of the orthogonality of the coincident modes takes D * i as one of the pair of coincident modes and ED * i + F D * j as the other.(The proof that D * i and ED * i + F D * j are orthogonal to all the non-coincident modes is omitted, because it is almost trivially simple, using Eqs.( 3) and ( 12)).
Because the modes are orthogonal, Eq. (3) gives Due to the modes having been normalised, Eq. (3) shows that the first term is equal to E. The second term is equal to F times the LHS of Eq. ( 2), = G say, with G calculated exactly as was done to get the value −0.545 for coincident modes given in the third column of Table 1.Hence Eq. ( 14) reduces to Eq. ( 15) which, since G is known, gives a unique value of −G for E/F and hence, after normalisation, a unique mode.

Conclusions
Three important conclusions follow from Eqs. (13) and ( 15): (1) two coincident modes found by the P RT method are unlikely to be mutually orthogonal despite being orthogonal to all other modes, e.g.see −0.545 in the third column of Table 1; (2) by choosing different values of A i , B i , A j and B j in Eq. ( 12), or different random force vectors for the P RT method, an infinite set of possible pairs of coincident modes can be found and; (3) any one of the pairs can be chosen to generate a pair of mutually orthogonal normalised modes which are also orthogonal to all the non-coincident modes of the structure, e.g. by choosing one of the two modes as one of the orthogonal pair and then using Eq. ( 15) to generate the other.Alternatively, the two modes can first be combined in any desired way to obtain a third mode and then Eq. ( 15) can be used to combine this third mode with either of the other two to obtain a mode orthogonal to it, e.g. to ensure that for a symmetric structure one mode is symmetric and the other antisymmetric, as in Figs.4(a) and (b).
The above arguments can be used recursively when more than two natural frequencies coincide.Thus for r coincident modes: (1) an infinite number of sets of r mutually orthogonal coincident modes can be generated; (2) for each set, one of the modes (or any com-bination of them) can be arbitrarily chosen as one of the orthogonal set; (3) Eq. ( 15) can then be applied to the remaining r − 1 modes in turn to generate a new set of r − 1 modes that are orthogonal to the arbitrarily chosen one, but not to each other; (4) the argument is recursive, i.e. any one of the r − 1 new modes can be chosen arbitrarily, with the remaining r − 2 modes used to generate r − 2 modes that are orthogonal to the first and second modes, but not to each other, etc, and; (5) the completion of this recursive process generates r coincident modes that are orthogonal to each other and to all the non-coincident modes, and involves r − 1 arbitrary choices being made.
For the simple frame of Fig. 1, several possible pairs of orthogonal coincident modes are readily visualised.For example, the four modes of Fig. 4 clearly satisfy moment equilibrium at the joints regardless of whether or not there are lumped masses or inertias there.Hence they are all possible coincident modes at the clamped/clamped member natural frequency.The symmetric mode (a) and anti-symmetric mode (b) clearly satisfy Eq. ( 2) and so are orthogonal to one another, because the products of the left-hand beam and column displacements cancel those of the right-hand beam and column, while products for the central column and any joint lumped masses or rotatory inertias are zero.Similarly (c) and (d) are a possible choice of orthogonal modes (and hence so also are their mirror images) because the products of the displacements sum to zero, having values for the three columns and then the two beams, working from left to right, of −V, +2V, 0, −V and 0, where V is a constant obtained by integration.Physically, mode (a) is a non-coincident mode if the central column is made minutely stiffer than the other members, so that modes (a) and (b) would be the pair of orthogonal coincident modes approached by reducing the stiffness of the central column to that of the other members.Similarly, mode (c) is a non-coincident mode if the right-hand beam is made minutely stiffer than the other members, so that modes (c) and (d) are the pair of orthogonal coincident modes approached as the righthand beam stiffness is reduced to those of the other members.

Fig. 1 .
Fig. 1.Plane frame example.Crosses denote the joints at which lumped masses or rotatory inertias were optionally added.

Fig. 4 .
Fig. 4. Two possible pairs of orthogonal coincident modes, showing the clamped/clamped member amplitudes.random force vector chosen.The modes D * i and D * j

Table 1
(2)ck on the correctness of the orthonormalisation conditions of Eq. (2), as judged by approximate satisfaction (denoted by Yes in the table) of Eq. (9) by the LHS of Eq.(2) * Expected, because two random coincident modes used.