Linearized equations of motion in multibody dynamics

Abstract

In this paper the linearized equations of motion in multibody dynamics are derived. Explicit expressions for the coefficient matrices are presented and given their physical interpretations. The equations of motion are presented in terms of the mechanical stiffness, its adjoint and the associated differential operators. It is demonstrated how the adjoint matrix may be used to find solutions to the associated algebraic eigenvalue problem. The case of multiple roots of the characteristic equation will result in a generalized eigenvalue problem involving the notion of a Jordan chain. Qualitative properties of the spectrum are derived without explicitly solving the characteristic equation. Finally, the mechanical admittance and its spectral representations are discussed.

Keywords

Multibody dynamics linearized equations vibrations multiple roots generalized eigenvalue problem

1. Introduction

The equations of motion in multibody dynamics constitute a system of ordinary differential equations of the second order. Typically they are non-linear and a detailed analysis of their solutions is usually very complicated. In many cases only qualitative aspects of their properties are accessible and for quantitative information one has to rely on numerical methods. If equilibrium configurations are present for the multibody system, solutions representing its motion in the vicinity of the equilibrium configuration may appear. The existence of these so-called mechanical vibrations is determined by the stability properties of the system at the equilibrium configuration. The analysis of mechanical vibrations is based on the linearized equations of motion. These equations are linear in the configuration coordinates, and their first- and second-order time derivatives. They are usually written in the standard format

M \ddot{x} + A \overset{\cdot}{x} + Bx = f

(1)

where the configuration coordinate $x = x (t) \in R^{n \times 1}$ represents the deviation of the system from the equilibrium configuration, $M \in R^{n \times n}$ is the mass matrix at the equilibrium configuration and $A, B \in R^{n \times n}$ are matrices representing dissipative and stiffness properties, as well as certain aspects of inertial and external forces that may appear as a consequence of, for instance, the choice of configuration coordinates. The column matrix $f = f (t) \in R^{n \times 1}$ , on the right-hand side of the equation, represents the specific time-dependent part of the external force on the system. Equation (1) is one of the classical equations in structural dynamics modelling the vibrations of structures with a finite number of degrees of freedom. Its history goes back to the early days of rational mechanics and the pioneering work by Euler and Lagrange. The status of the subject, as of the late 19th century, is accounted for by Lord Rayleigh [1]. Theories and methods presented in this classical and often-cited work are still of significance in the field. More recent numerical methods, such as the finite element method, have renewed the interest in the equation concerning both formal mathematical aspects as well as numerical techniques (see Hughes [2]).

The formulation of the linear ordinary differential equation (1) requires a calculation of the coefficient matrices $M$ , $A$ and $B$ . Starting from the basic equations of motion in multibody dynamics, one identifies the equilibrium solutions and then performs a formal linearization of the equations at an equilibrium configuration. The mathematical character of the coefficient matrices will then follow from physical properties of the system and from the method underlying the formulation of the basic equations. If the derivation of the basic equations is obtained by using the method of Lagrange, the mass matrix $M$ will always be symmetric (see, for instance, Géradin and Rixen [3]). For many multibody systems this will be the case for the matrices $A$ and $B$ as well. However, there are important systems where this is not so.

In this paper we consider a multibody system consisting of rigid and visco-elastic parts connected by ideal joints. Starting from the basic equations of motion, as presented by Lidström [4], linearized equations at an equilibrium configuration are derived. The coefficient matrices $M, A$ and $B$ are then identified and characterized in terms of their general mathematical properties and physical interpretations.

The differential equation (1) is closely related to the matrix polynomial

S (s) = M s^{2} + As + B, s \in C

(2)

called the dynamic stiffness of the multibody system at the equilibrium configuration. Equation (1) may then be written

S (d_{t}) x = ({Md}_{t}^{2} + A d_{t} + B) x = f

(3)

where $d_{t}$ is the ordinary differential operator with respect to time. It is well known that the associated eigenvalue problem: find $s \in C$ and $w \neq 0_{n \times 1}$ such that

S (s) w = 0_{n \times 1}

(4)

plays an important role for the solution of Equation (3). The existence of a non-trivial solution to the eigenvalue problem requires that $s$ satisfies the characteristic equation $p (s) = 0$ , where $p (s) = det S (s)$ is the characteristic polynomial associated with the dynamic stiffness. The roots of the characteristic equation define the spectrum of the differential operator (3). The adjoint $T (s) = adj S (s)$ associated with the dynamic stiffness is the unique matrix satisfying

S (s) T (s) = T (s) S (s) = p (s) 1_{n \times n}, s \in C

(5)

In this paper we make a systematic use of the adjoint matrix and its related adjoint differential operator $T (d_{t})$ to construct solutions to (3). This is a different approach than usually taken in books on ordinary differential equations such as, for instance, Birkhoff and Rota [5] and Arnold [6]. The standard procedure there is to transform Equation (3) into a system of first-order differential equation, the so-called state space formulation. This is, from the mathematical point of view, a rational and unifying approach. The advantage of the method used here is that it retains the original differential equation (3) with its close connection to the basic principles of mechanics. We suspect that the approach of using the adjoint, if not entirely new, is not too well known in the engineering community.

The objective of this paper is to derive the linearized equations of motion in multibody dynamics, to analyse their properties and characterize their solutions in the cases where the coefficient matrices $A$ and $B$ represent properties of the mechanical system that are not too special. The classical solution procedure of Equation (3) is to try to find a linear coordinate transformation that will uncouple the differential equations (see Foss [7] and Caughey [8]). This procedure works if the coefficient matrices $A$ and $B$ are symmetric and if, in addition, $A M^{- 1} B = B M^{- 1} A$ (see Caughey and O’Kelly [9]). A situation of this kind is, for instance, at hand in the case of so-called proportional damping. If $A$ or $B$ are not symmetric the situation is more complicated. Results covering some of these more general cases are found in Lancaster [10], a work that will be frequently referred to in this paper.

There are many textbooks on mechanical vibrations. One good example is the previously mentioned book by Géradin and Rixen [3], where some attention is given to the derivation of the linearized equation. Concerning the properties of the solutions to the equations, the literature is abundant and a full review will not be undertaken here. Most papers deal with special applications often resulting in equations where the matrices $A$ and $B$ are symmetric and they use solution methods involving changes of coordinates in order to obtain uncoupled or nearly uncoupled equations. The presentation of more general methods, in order to find solutions to the equations, will be one of the main themes of this paper.

The paper is organized as follows. In Section 2 the basic notations and definitions used throughout the paper are reviewed. In Section 3 the multibody system and its configuration coordinates are introduced along with the equations of motion. Section 4 specifies the interaction between parts and constitutive assumptions for rigid and visco-elastic parts. Section 5 defines the equilibrium configuration and derives the linearized equation of motion. The coefficient matrices are calculated and given their appropriate physical interpretations. A linearized version of the power theorem is derived. In Section 6 general properties of the linearized equation are discussed involving, among other things, the existence and uniqueness of solutions. Two elementary examples, illustrating the theory, are presented. In Section 7 free vibrations are analysed using the adjoint differential operator. This will result in the so-called characteristic differential equation. Solutions of this equation, which are well known in the mathematical literature, are accounted for as exponential-polynomials. In Section 8 the important case of simple roots to the characteristic equation is discussed. This will, for instance, result in the eigenvalue problem defined in (4) above. It is demonstrated how the adjoint matrix may be used to find solutions to this eigenvalue problem. In Section 9 the case of multiple roots is discussed. This will result in a generalized eigenvalue problem involving the notion of a Jordan chain. The case of a so-called simple dynamic stiffness matrix is given special attention. In Section 10 the dynamic flexibility is defined as the inverse of the dynamic stiffness. It is evaluated in the case of a simple dynamic stiffness and related to modal properties. In the general case a basic formula for the flexibility is derived, but without an explicit relation to its modal properties. In Section 11 qualitative properties of the spectrum are discussed. This is done without explicitly solving the characteristic equation. The results are illustrated in a simple example. In the final section, Section 12, forced vibrations are characterized using the dynamic flexibility.

To facilitate the reading of this paper some of the mathematical notation, definitions and results used in the main text are filed in an Appendix.

2. Notations

In this paper $N$ denotes the set of natural numbers (not including $0$ ), $R$ denotes the set of real numbers, $R^{+}$ denotes the set of positive real numbers and $C$ the set of complex numbers. If $z \in C$ then $z^{*} \in C$ denotes its complex conjugate. Detailed information on matrix notation and definitions are collected in the Appendix.

Let $E$ denote a three-dimensional Euclidean point space with the corresponding translation vector space $V$ , that is, $\dim (V) = 3$ . Points in $E$ are denoted by $x, y, \dots, X, Y, \dots$ and vectors in $V$ by $a, b, \dots, u, v, \dots$ . The scalar product of two vectors $a and b$ is denoted $a \cdot b$ . The norm $| a |$ of a vector a is defined by $| a | = \sqrt{a \cdot a}$ . The boundary of a set $D \subseteq E$ is denoted $\partial D$ . The space of all second-order tensors A on $V$ , that is, linear mappings $V \to V$ , is denoted $End (V)$ . We write $1_{v}$ for the unit tensor in $End (V)$ and $A^{T}$ for the transpose of $A \in End (V)$ . The space of all symmetric tensors is denoted by $Sym (V) =$ ${A \in End (V) | A^{T} = A}$ and the space of skew-symmetric tensor is denoted $Skew (V) = {A \in End (V) | A^{T} = - A}$ . The scalar product, $A \cdot B$ , of two tensors $A, B \in End (V)$ is defined by $A \cdot B = tr (A^{T} B)$ and the corresponding norm $| A | = \sqrt{A \cdot A}$ .

3. Equations of motion in multibody dynamics

This introductory section follows closely that of Lidström [4]. For more information on details and motivations for the present approach the reader is referred to that paper.

Consider a multibody $B$ consisting of $N$ rigid or elastic parts: $B^{α}, α = 1, \dots, N$ . Let the regular region $B_{t}^{α}$ in $E$ denote the present (spatial) placement of $B^{α}$ , at time t. The function representing the transplacement of the part $B^{α}$ from its reference placement $B_{0}^{α}$ to its present placement $B_{t}^{α}$ is given by the mapping $χ^{α} : B_{0}^{α} \times R \to B_{t}^{α}$ , where

x = χ^{α} (X, t), (X, t) \in B_{0}^{α} \times R

is the present place (at time $t$ ) of the material point $P \in B$ with its referential place denoted by $X \in B_{0}^{α}$ . We assume that it is possible to introduce configuration coordinates

q = (q^{1}, q^{2}, \dots, q^{n}) \in Ω \subseteq R^{n}

where $Ω$ is an open set, so that the set of possible placements (configurations) of the parts of the multibody may (locally) be represented by the mappings

x = χ^{α} (X, t) = χ_{q}^{α} (t, q^{1}, q^{2}, \dots, q^{n}; X) = χ_{q}^{α} (t, q; X), (t, q; X) \in R \times Ω \times B_{0}^{α}

where $χ_{q}^{α} : R \times Ω \times B_{0}^{α} \to E$ is assumed to be a twice continuously differentiable mapping in all its arguments $(t, q; X)$ . The deformation gradient $F^{α}$ and the Green–St. Venant strain tensor $E^{α}$ are defined by $F^{α} = \partial_{X} χ^{α}$ and $E^{α} = \frac{1}{2} ({F^{α}}^{T} F^{α} - 1_{V})$ , respectively.

The motion of the multibody is given by a function $q = q (t), t \in [0, \infty [$ and a mapping

x = χ_{q} (t, q (t); X) \overset{def}{=} χ_{q}^{α} (t, q (t); X), t \in [0, \infty [, X \in B_{0}^{α}

Let $ρ_{0}^{α} = ρ_{0}^{α} (X), X \in B_{0}^{α}$ denote the mass density of part $B^{α}$ in its reference placement. The kinetic energy for the multibody $B$ may be written

T = T (t, q, \overset{\cdot}{q}) = \frac{1}{2} (M_{0} + M_{1} \overset{\cdot}{q} + {\overset{\cdot}{q}}^{T} M_{2} \overset{\cdot}{q})

(6)

where $M_{0} = \sum_{α = 1}^{N} M_{0}^{α}$ , $M_{1} = \sum_{α = 1}^{N} M_{1}^{α}$ , $M_{2} = \sum_{α = 1}^{N} M_{2}^{α}$ and

M_{0}^{α} = M_{0}^{α} (t, q) = m_{00}^{α} \in R, M_{1}^{α} = M_{1}^{α} (t, q) = (\begin{matrix} m_{01}^{α} + m_{10}^{α} & m_{02}^{α} + m_{20}^{α} & \dots & m_{0 n}^{α} + m_{n 0}^{α} \end{matrix}) \in R^{1 \times n},

M_{2}^{α} = M_{2}^{α} (t, q) = (\begin{matrix} m_{11}^{α} & m_{12}^{α} & \dots & m_{1 n}^{α} \\ m_{21}^{α} & m_{22}^{α} & \dots & m_{2 n}^{α} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ m_{n 1}^{α} & m_{n 2}^{α} & \dots & m_{nn}^{α} \end{matrix}) \in R^{n \times n}

where the matrix elements $m_{kl}^{α}$ are defined by

m_{kl}^{α} = m_{kl}^{α} (t, q) = \int_{B_{0}^{α}} \frac{\partial χ_{q}^{α}}{\partial q^{k}} \cdot \frac{\partial χ_{q}^{α}}{\partial q^{l}} ρ_{0}^{α} (X) dv (X), k, l = 0, 1, \dots, n

The matrix $M_{2}$ is symmetric and positive semi-definite. We will here assume that the system of configuration coordinates is regular, which implies that $M_{2}$ is positive definite (see Lidström [4]).

The equations of motion for the multibody may now, according to Lidström [4], be written

{\begin{matrix} {\ddot{q}}^{T} M_{2} = Q^{sum} \\ q (0) = q_{0}, \overset{\cdot}{q} (0) = {\overset{\cdot}{q}}_{0} \end{matrix}

(7)

where the sum of the generalized forces $Q^{sum} = Q^{sum} (t, q, \overset{\cdot}{q})$ represents the sum of external as well as internal forces acting on the multibody system. It is convenient to express the force sum according to

Q^{sum} = Q^{cif} + Q^{i} + Q^{c} + Q^{b}

(8)

where $Q^{cif} = Q^{cif} (t, q, \overset{\cdot}{q})$ , the complementary inertia force, represents inertial forces such as centrifugal and Coriolis forces, $Q^{i} = Q^{i} (t, q, \overset{\cdot}{q})$ is the internal force, $Q^{c} = Q^{c} (t, q, \overset{\cdot}{q})$ is the contact force and $Q^{b} = Q^{b} (t, q)$ is the body force due to gravity. The complementary inertia force is defined by (see the Appendix in Lidström [4] for a definition of derivatives of matrix valued functions)

Q^{cif} = - ({\overset{\cdot}{q}}^{T} (\frac{\partial M_{2}}{\partial q} - \frac{1}{2} {(\frac{\partial M_{2}}{\partial q})}^{t}) \overset{\cdot}{q} + {\overset{\cdot}{q}}^{T} (\frac{\partial M_{2}}{\partial t} + skew (\frac{\partial M_{1}}{\partial q})) + \frac{1}{2} (\frac{\partial M_{1}}{\partial t} - \frac{\partial M_{0}}{\partial q}))

(9)

and the internal force by

Q^{i} = Q^{i} (t, q, \overset{\cdot}{q}) = \sum_{α = 1}^{N} Q^{α, i}

where

Q^{α, i} = - \int_{B_{0}^{α}} S^{α} \cdot \frac{\partial E_{q}^{α}}{\partial q} dv (X)

(10)

is the internal force corresponding to part $B^{α}$ . In Equation (10) $E_{q}^{α} = E_{q}^{α} (t, q; X)$ is the Green–St. Venant strain tensor defined above and $S^{α}$ denotes the second Piola–Kirchhoff stress tensor. The generalized contact force is defined by

Q^{c} = \sum_{α = 1}^{N} Q^{α, c}, Q^{α, c} = \int_{\partial B_{0}^{α}} \frac{\partial χ_{q}^{α}}{\partial q} \cdot t_{0}^{α} da (X)

(11)

where $t_{0}^{α}$ is the traction vector on the boundary surface $\partial B_{0}^{α}$ of part $B^{α}$ . We may employ the decomposition $Q^{α, c} = Q^{α, ic} + Q^{α, ec}$ , where $Q^{α, ic} = Q^{α, ic} (t, q, \overset{\cdot}{q})$ is the internal contact force, that is, the contact force acting on part $B^{α}$ from all other parts of $B$ and $Q^{α, ec} = Q^{α, ec} (t, q, \overset{\cdot}{q})$ is the external contact force on part $B^{α}$ , from the exterior of $B$ . This force is assumed to be prescribed. We may write

Q^{α, ic} = \int_{\partial^{ic} B_{0}^{α}} \frac{\partial χ_{q}^{α}}{\partial q} \cdot t_{0}^{α, ic} da (X), Q^{α, ec} = \int_{\partial^{ec} B_{0}^{α}} \frac{\partial χ_{q}^{α}}{\partial q} \cdot t_{0}^{α, ec} da (X)

where $t_{0}^{α, ec}$ is the prescribed traction vector acting on $B^{α}$ from the exterior of $B$ and $\partial^{ec} B_{0}^{α}$ represents the boundary of $B^{α}$ in contact with the exterior of $B$ , while $\partial^{ic} B_{0}^{α}$ is the contact surfaces between $B^{α}$ and all other parts of $B$ and $t_{0}^{α, ic}$ is the corresponding traction vector. With

Q^{ic} = \sum_{α = 1}^{N} Q^{α, ic} and Q^{ec} = \sum_{α = 1}^{N} Q^{α, ec}

it follows that $Q^{c} = Q^{ic} + Q^{ec}$ . One may write

Q^{ic} = \frac{1}{2} \sum_{α, β = 1}^{N} I^{α β}

(12)

where the interaction $I^{α β}$ , between parts $B^{α}$ and $B^{β}$ , is given by

I_{k}^{β α} = I_{k}^{α β} = \int_{S_{t}^{α β}} x_{k}^{α, β} \cdot t^{α β} da (x), k = 1, \dots, n

(13)

where $S_{t}^{α β}$ is the contact surface between parts $B^{α}$ and $B^{β}$ in their present placements, $x_{k}^{α, β} \overset{def}{=} x_{k}^{α} - x_{k}^{β}$ and $x_{k}^{α} (t, q, x) = \frac{\partial χ_{q}^{α}}{\partial q^{k}} (t, q (t); (χ_{q}^{α})^{- 1} (t, q; x)), α = 1, \dots, N, k = 1, \dots, n$ and $t^{α β}$ is the traction vector acting on $B^{α}$ from $B^{β}$ in their present placement (see Lidström [4]).

The generalized body force $Q^{b} = Q^{b} (t, q)$ is here assumed to be due to gravity. Then

Q^{b} = \sum_{α = 1}^{N} Q^{α, b}, Q_{k}^{α, b} = \int_{B_{0}^{α}} \frac{\partial χ_{q}^{α}}{\partial q^{k}} \cdot g ρ_{0}^{α} dv (X)

where $g$ is the (constant) acceleration of gravity. Thus

Q^{b} = - \frac{\partial V_{g}}{\partial q}, V_{g} = V_{g} (t, q) = \sum_{α = 1}^{N} (x_{c}^{α} - x_{O}) \cdot g m^{α}

(14)

where $m^{α}$ is the mass of part $B^{α}$ and

x_{c}^{α} = x_{c}^{α} (t, q) = x_{O} + \frac{1}{m^{α}} \int_{B_{0}^{α}} (χ_{q}^{α} (t, q; X) - x_{O}) ρ_{0}^{α} dv (X)

is the present place of the centre of mass of part $B^{α}$ . The point $O$ is assumed to have a fixed place $x_{O}$ in Euclidean space.

4. Constitutive assumptions and interactions between parts

The constitutive equation for a visco-elastic part $B^{α}$ is given by

S^{α} (X, t) = C^{α} (X) [E^{α} (X, t)] + D^{α} (X) [{\overset{\cdot}{E}}^{α} (X, t)]

where $C^{α} = C^{α} (X) : Sym (V) \to Sym (V)$ and $D^{α} = D^{α} (X) : Sym (V) \to Sym (V)$ are symmetric, positive semi-definite linear mappings, that is, $A \cdot C^{α} (X) [B] = B \cdot C^{α} (X) [A],$ $A, B \in Sym (V), X \in B_{0}^{α}$ and $A \cdot D^{α} (X) [A] \geq 0, A \in Sym (V), X \in B_{0}^{α}$ . Then

Proposition 4.1

Q^{i} = Q^{i} (t, q, \overset{\cdot}{q}) = - \frac{\partial V_{e} (t, q)}{\partial q} - \frac{\partial D_{e} (t, q, \overset{\cdot}{q})}{\partial \overset{\cdot}{q}}

(15)

where

V_{e} (t, q) = \frac{1}{2} \sum_{α = 1}^{N} \int_{B_{0}^{α}} E_{q}^{α} \cdot C^{α} [E_{q}^{α}] dv (X), D_{e} (t, q, \overset{\cdot}{q}) = \frac{1}{2} \sum_{α = 1}^{N} \int_{B_{0}^{α}} {\overset{\cdot}{E}}_{q}^{α} \cdot D^{α} [{\overset{\cdot}{E}}_{q}^{α}] dv (X)

(16)

are the elastic energy and the Rayleigh dissipation function of the multibody $B$ , respectively.

Proof: We note that

{\overset{\cdot}{E}}_{q}^{α} (t, q, \overset{\cdot}{q}; X) = \frac{\partial E_{q}^{α} (t, q; X)}{\partial t} + \frac{\partial E_{q}^{α} (t, q; X)}{\partial q} \overset{\cdot}{q}

(17)

and consequently $\frac{\partial {\overset{\cdot}{E}}_{q}^{α}}{\partial \overset{\cdot}{q}} = \frac{\partial E_{q}^{α}}{\partial q}$ . Then, considering part $B^{α}$ ,

\begin{matrix} Q^{α, i} (t, q, \overset{\cdot}{q}) = - \int_{B_{0}^{α}} S^{α} \cdot \frac{\partial E_{q}^{α}}{\partial q} dv (X) = - \int_{B_{0}^{α}} (C^{α} [E_{q}^{α}] + D^{α} [{\overset{\cdot}{E}}_{q}^{α}]) \cdot \frac{\partial E_{q}^{α}}{\partial q} dv (X) = \\ - \int_{B_{0}^{α}} C^{α} [E_{q}^{α}] \cdot \frac{\partial E_{q}^{α}}{\partial q} dv (X) - \int_{B_{0}^{α}} D^{α} [{\overset{\cdot}{E}}_{q}^{α}] \cdot \frac{\partial {\overset{\cdot}{E}}_{q}^{α}}{\partial \overset{\cdot}{q}} dv (X) = \\ - \frac{\partial}{\partial q} \frac{1}{2} \int_{B_{0}^{α}} E^{α} \cdot C^{α} [E^{α}] dv (X) - \frac{\partial}{\partial \overset{\cdot}{q}} \frac{1}{2} \int_{B_{0}^{α}} {\overset{\cdot}{E}}^{α} \cdot D^{α} [{\overset{\cdot}{E}}^{α}] dv (X) = - \frac{\partial V_{e}^{α}}{\partial q} - \frac{\partial D_{e}^{α}}{\partial \overset{\cdot}{q}} \end{matrix}

where $V_{e}^{α} = V_{e}^{α} (t, q)$ and $D_{e}^{α} = D_{e}^{α} (t, q, \overset{\cdot}{q})$ are defined in (16). Consequently, the generalized internal force for the multibody may then be written

Q^{i} = Q^{i} (t, q, \overset{\cdot}{q}) = \sum_{α = 1}^{N} Q^{α, i} = \sum_{α = 1}^{N} (- \frac{\partial V_{e}^{α}}{\partial q} - \frac{\partial D_{e}^{α}}{\partial \overset{\cdot}{q}}) = - \frac{\partial V_{e}}{\partial q} - \frac{\partial D_{e}}{\partial \overset{\cdot}{q}}

where

V_{e} = \sum_{α = 1}^{N} V_{e}^{α}, D_{e} = \sum_{α = 1}^{N} D_{e}^{α}

(18)

This proves the proposition. □

Note that if part $B^{α}$ is rigid then $E^{α} = 0$ , ${\overset{\cdot}{E}}^{α} = 0$ and, consequently, in this case, $Q^{α, i} = 0_{1 \times n}$ , which means that rigid parts do not contribute to the generalized internal force of the multibody.

Proposition 4.2 The Rayleigh dissipation function

D_{e} = \frac{1}{2} (D_{e, 0} + D_{e, 1} \overset{\cdot}{q} + {\overset{\cdot}{q}}^{T} D_{e, 2} \overset{\cdot}{q})

(19)

where

D_{e, 0} = D_{e, 0} (t, q) = \sum_{α = 1}^{N} \int_{B_{0}^{α}} \frac{\partial E_{q}^{α}}{\partial t} \cdot D^{α} [\frac{\partial E_{q}^{α}}{\partial t}] dv (X) \in R

D_{e, 1} = D_{e, 1} (t, q) = 2 \sum_{α = 1}^{N} \int_{B_{0}^{α}} \frac{\partial E_{q}^{α}}{\partial t} \cdot D^{α} [\frac{\partial E_{q}^{α}}{\partial q}] dv (X) \in R^{1 \times n}

(20)

D_{e, 2} = D_{e, 2} (t, q) = \sum_{α = 1}^{N} \int_{B_{0}^{α}} {(\frac{\partial E_{q}^{α}}{\partial q})}^{T} \cdot D^{α} [\frac{\partial E_{q}^{α}}{\partial q}] dv (X) \in R^{n \times n}

and $D_{e, 2}$ is a symmetric, positive semi-definite matrix. Consult the Appendix A.1 for a definition of the integrands appearing in (20).

Proof: A demonstration is obtained by inserting (17) into (16)₂. □

Interactions between parts are mediated by the internal contact forces $t^{α β}$ and summarized in $Q^{ic}$ according to Equation (12). Characteristic of multibody systems is the presence of joints that impose constraints on the relative motion between parts of the system. Many joints, in practical applications, may be modelled in terms of so-called simple joints; examples are the cylindrical, prismatic, screw, revolute, spherical and planar joints. In this paper we assume that the q-coordinate system is compatible with the constraints, that is, the constraint conditions are identically satisfied for all $(t, q) \in [0, \infty [\times Ω, \overset{\cdot}{q} \in R^{n}$ . Then, according to Lidström [4]:

I_{k}^{α β} = \int_{S_{t}^{α β}} x_{k}^{α, β} \cdot τ^{α β} da (x), k = 1, \dots, n

where $τ^{α β}$ denotes the tangential composant of the contact force $t^{α β}$ , that is, $τ^{α β} \cdot n^{α β} = 0$ where $n^{α β}$ denotes the orientation of the contact surface $S_{t}^{α β}$ . We now assume that

x_{k}^{α, β} (x, t) \cdot τ^{α β} (x, t) = 0, (x, t) \in S_{t}^{α β} \times [0, \infty [

(21)

The interaction is then said to be ideal and it follows that $I_{k}^{α β} = 0, k = 1, \dots, n$ and according to (12) we then have $Q^{ic} = 0_{1 \times n}$ and consequently

Q^{sum} = Q^{cif} + Q^{i} + Q^{b} + Q^{ec}

(22)

Ideal joints may be equipped with visco-elastic torsion bushing elements that will add to the interactions of the multibody (see Lidström [11]). The interaction may also be mediated by a torque obtained from an electric motor acting between parts $B^{α}$ and $B^{β}$ . These possibilities will, however, not be considered in this paper.

Remark 4.1: The assumption in (21) rules out the influence of Coulomb friction at contact surfaces. The incorporation of Coulomb friction would require an expression for $τ^{α β}$ of the form

τ^{α β} = μ_{kin}^{α β} N^{α β} \frac{v^{α β}}{v^{α β}}, v^{α β} \neq 0

where $μ_{kin}^{α β}$ is the coefficient of kinematic friction at the contact surface $S_{t}^{α β}$ between parts $B^{α}$ and $B^{β}$ , $N^{α β} = N^{α β} (x, t)$ is the normal component of the contact traction, that is, $t^{α β} = n^{α β} N^{α β} + τ^{α β}, N^{α β} \leq 0$ , $v^{α β}$ is the relative tangential material velocity at the contact surface, that is, $v^{α β} = {\overset{\cdot}{x}}^{α, β} - n^{α β} (n^{α β} \cdot {\overset{\cdot}{x}}^{α, β})$ and $v^{α β} = | v^{α β} |$ . If the q-coordinate system is compatible with the constraints then (see Lidström [4])

v^{α β} = v_{0}^{α β} + \sum_{k = 1}^{n} x_{k}^{α, β} {\overset{\cdot}{q}}^{k}, v_{0}^{α β} = x_{0}^{α, β} - n^{α β} (n^{α β} \cdot x_{0}^{α, β})

Thus $\frac{\partial v^{α β}}{\partial {\overset{\cdot}{q}}^{k}} = x_{k}^{α, β}$ and then

I_{k}^{α β} = \int_{S_{t}^{α β}} \frac{\partial v^{α β}}{\partial {\overset{\cdot}{q}}^{k}} \cdot μ_{kin}^{α β} N^{α β} \frac{v^{α β}}{v^{α β}} da (x) = \int_{S_{t}^{α β}} μ_{kin}^{α β} N^{α β} \frac{v^{α β}}{v^{α β}} \cdot \frac{\partial v^{α β}}{\partial {\overset{\cdot}{q}}^{k}} da (x) = \int_{S_{t}^{α β}} μ_{kin}^{α β} N^{α β} \frac{\partial v^{α β}}{\partial {\overset{\cdot}{q}}^{k}} da (x)

(23)

where we have used the identity $v^{α β} \cdot \frac{\partial v^{α β}}{\partial {\overset{\cdot}{q}}^{k}} = v^{α β} \frac{\partial v^{α β}}{\partial {\overset{\cdot}{q}}^{k}}$ . If we introduce the (heuristic) assumption that $μ_{kin}^{α β}$ and $N^{α β}$ do not depend on $\overset{\cdot}{q}$ then (see Géradin and Rixen [3], p. 24–25)

I_{k}^{α β} = \frac{\partial}{\partial {\overset{\cdot}{q}}^{k}} \int_{S_{t}^{α β}} μ_{kin}^{α β} N^{α β} v^{α β} da (x)

and we may write $Q^{ic} = \frac{\partial D_{ic}}{\partial \overset{\cdot}{q}}$ where

D_{ic} = D_{ic} (t, q, \overset{\cdot}{q}) = \frac{1}{2} \sum_{α, β = 1}^{N} \int_{S_{t}^{α β}} μ_{kin}^{α β} N^{α β} v^{α β} da (x)

is the Rayleigh dissipation function corresponding to the internal contact forces due to Coulomb friction. The assumption that $N^{α β}$ , in general, does not depend on $\overset{\cdot}{q}$ is most certainly not correct. Furthermore, $μ_{kin}^{α β}$ may, in more advanced friction models, depend on $N^{α β}$ . A calculation of $I^{α β}$ requires, according to (23), a calculation of the normal force $N^{α β} = N^{α β} (x, t)$ . In order to obtain this one has to introduce a q - coordinate system, which is not compatible with the constraints (see Lidström [12]). We do not consider this much more complicated situation in the present paper.■

The sum of the generalized forces (except $Q^{cif}$ ) may, if $Q^{b}$ and $Q^{i}$ are given by (14) and (15), be written

Q^{i} + Q^{b} + Q^{ec} = - \frac{\partial V_{e}}{\partial q} - \frac{\partial D_{e}}{\partial \overset{\cdot}{q}} - \frac{\partial V_{g}}{\partial q} + Q^{ec} = - \frac{\partial V}{\partial q} - \frac{\partial D_{e}}{\partial \overset{\cdot}{q}} + Q^{ec}

where

V (t, q) = V_{e} (t, q) + V_{g} (t, q)

is the potential energy of the multibody. We have the following division into conservative and non-conservative forces:

Q^{i} + Q^{b} + Q^{ec} = Q^{con} + Q^{non}

where

Q^{con} = - \frac{\partial V}{\partial q}, Q^{non} = - \frac{\partial D_{e}}{\partial \overset{\cdot}{q}} + Q^{ec}

The power expended by these forces is given by

P^{con} = Q^{con} \overset{\cdot}{q} = - \frac{\partial V}{\partial q} \overset{\cdot}{q}, P^{non} = (- \frac{\partial D_{e}}{\partial \overset{\cdot}{q}} + Q^{ec}) \overset{\cdot}{q}

The mechanical energy, E, of the multibody is defined by

E = p \overset{\cdot}{q} - L

where $L = L (t, q, \overset{\cdot}{q}) = T (t, q, \overset{\cdot}{q}) - V (t, q)$ is the Lagrangian, and $p = \frac{\partial T}{\partial \overset{\cdot}{q}}$ is the so-called generalized momentum.

Proposition 4.3 The power theorem for the multibody reads

\overset{\cdot}{E} = - \frac{1}{2} D_{e, 1} \overset{\cdot}{q} - {\overset{\cdot}{q}}^{T} D_{e, 2} \overset{\cdot}{q} + Q^{ec} \overset{\cdot}{q} - \frac{\partial L}{\partial t}

Proof: We have (see Lidström [4], Section 8)

\overset{\cdot}{E} = P^{non} + P^{ic} - \frac{\partial L}{\partial t}

where in this case $P^{ic} = Q^{ic} \overset{\cdot}{q} = 0$ and $P^{non}$ , the power of non-conservative forces, is given by

P^{non} = (- \frac{\partial D_{e}}{\partial \overset{\cdot}{q}} + Q^{ec}) \overset{\cdot}{q} = (- \frac{1}{2} D_{e, 1} - {\overset{\cdot}{q}}^{T} D_{e, 2} + Q^{ec}) \overset{\cdot}{q}

which proves the proposition.■

The power theorem, in the previous proposition, may be reformulated according to the following.

Corollary 4.1

{\overset{\cdot}{E}}^{*} = - \frac{1}{2} D_{e, 1} \overset{\cdot}{q} - {\overset{\cdot}{q}}^{T} D_{e, 2} \overset{\cdot}{q} + Q^{ec} \overset{\cdot}{q} - \frac{\partial L}{\partial t}

(24)

where

E^{*} = \frac{1}{2} {\overset{\cdot}{q}}^{T} M_{2} \overset{\cdot}{q} + V^{*}

is the relative mechanical energy and

V^{*} = V^{*} (t, q) = V (t, q) - \frac{1}{2} M_{0} (t, q)

(25)

is called the modified potential energy of the multibody.

Proof: We have, according to (6)

p = \frac{\partial T}{\partial \overset{\cdot}{q}} = \frac{1}{2} (M_{1} + 2 {\overset{\cdot}{q}}^{T} M_{2})

and then

E = p \overset{\cdot}{q} - L = \frac{1}{2} (M_{1} + 2 {\overset{\cdot}{q}}^{T} M_{2}) \overset{\cdot}{q} - \frac{1}{2} (M_{0} + M_{1} \overset{\cdot}{q} + {\overset{\cdot}{q}}^{T} M_{2} \overset{\cdot}{q}) + V = - \frac{1}{2} M_{0} + \frac{1}{2} {\overset{\cdot}{q}}^{T} M_{2} \overset{\cdot}{q} + V

Invoking Proposition 4.3 one obtains (24). □

5. Linearization of the equations of motion

The equations of motion for the multibody are in this paper given by Equations (7) and (8).

The multibody is said to be non-disturbed if the external force $Q^{ec} = 0_{1 \times n}$ . A fixed configuration $q = q_{eq} \in Ω$ is an equilibrium configuration to a non-disturbed system if $q (t) = q_{eq}, t \in [0, \infty [$ is a solution to (7) with $q_{0} = q_{eq}$ and ${\overset{\cdot}{q}}_{0} = 0_{n \times 1}$ . If we insert $q (t) = q_{eq}, (\overset{\cdot}{q} (t) = 0, \ddot{q} (t) = 0)$ into (7), one obtains the equilibrium equation:

0_{1 \times n} = Q^{sum} (t, q_{eq}, 0_{n \times 1}) = Q^{cif} (t, q_{eq}, 0_{n \times 1}) + Q^{i} (t, q_{eq}, 0_{n \times 1}) + Q^{b} (t, q_{eq})

(26)

Proposition 5.1 If $Q^{b}$ and $Q^{i}$ are given by (14) and (15), respectively, then the equilibrium equation (26) may be written

\frac{\partial V^{*} (t, q_{eq})}{\partial q} + \frac{1}{2} \frac{\partial M_{1} (t, q_{eq})}{\partial t} + \frac{1}{2} D_{e, 1} (t, q_{eq}) = 0_{1 \times n}

(27)

Proof: We have, according to (9):

Q^{cif} (t, q_{eq}, 0_{n \times 1}) = - \frac{1}{2} (\frac{\partial M_{1} (t, q_{eq})}{\partial t} - \frac{\partial M_{0} (t, q_{eq})}{\partial q})

From (15) and (19) it follows that

Q^{i} (t, q_{eq}, 0_{n \times 1}) = - \frac{\partial V_{e} (t, q_{eq})}{\partial q} - \frac{\partial D_{e} (t, q_{eq}, 0_{n \times 1})}{\partial \overset{\cdot}{q}} = - \frac{\partial V_{e} (t, q_{eq})}{\partial q} - \frac{1}{2} D_{e, 1} (t, q_{eq})

Q^{b} (t, q_{eq}) = - \frac{\partial V_{b} (t, q_{eq})}{\partial q}

The equilibrium equation (26) may then be written

\frac{1}{2} (\frac{\partial M_{1} (t, q_{eq})}{\partial t} - \frac{\partial M_{0} (t, q_{eq})}{\partial q}) + \frac{\partial V_{e} (t, q_{eq})}{\partial q} + \frac{1}{2} D_{e, 1} (t, q_{eq}) + \frac{\partial V_{b} (t, q_{eq})}{\partial q} = 0_{1 \times n}

\frac{\partial V^{*} (t, q_{eq})}{\partial q} + \frac{1}{2} \frac{\partial M_{1} (t, q_{eq})}{\partial t} + \frac{1}{2} D_{e, 1} (t, q_{eq}) = 0_{1 \times n}

which proves the proposition. □

Next we introduce the configuration coordinate $x = q - q_{eq}$ representing the deviation from the equilibrium solution. For $q_{eq} \in Ω$ we introduce the open set $U_{q_{eq}} = {x \in R^{n \times 1} | q_{eq} + x \in Ω}$ .

Equation (7)₁ may be written

{\ddot{x}}^{T} M_{2} (t, q_{eq} + x) = Q^{sum} (t, q_{eq} + x, \overset{\cdot}{x}), x \in U_{q_{eq}}

where $Q^{sum}$ now includes the external force $Q^{ec}$ and is given by (22). We have, using a Taylor expansion:

M_{2} (t, q_{eq} + x) = M_{2} (t, q_{eq}) + \frac{\partial M_{2} (t, q_{eq})}{\partial q} x + ε_{M} (t, x) ‖ x ‖

and $ε_{M} : [0, \infty [\times U_{q_{eq}} \to sym (R^{n \times n})$ is continuous at $x = 0_{n \times 1}$ and $ε_{M} (t, 0_{n \times 1}) = 0_{n \times n}$ . Furthermore

\begin{matrix} Q^{sum} (t, q_{eq} + x, \overset{\cdot}{x}) = Q^{sum} (t, q_{eq}, 0_{n \times 1}) + \frac{\partial Q^{sum} (t, q_{eq}, 0_{n \times 1})}{\partial q} x + \frac{\partial Q^{sum} (t, q_{eq}, 0_{n \times 1})}{\partial \overset{\cdot}{q}} \overset{\cdot}{x} + \\ ε_{Q} (t, x, \overset{\cdot}{x}) \sqrt{{‖ x ‖}^{2} + {‖ \overset{\cdot}{x} ‖}^{2}} = Q^{ec} (t, q_{eq}, 0_{n \times 1}) + \frac{\partial Q^{sum} (t, q_{eq}, 0_{n \times 1})}{\partial q} x + \frac{\partial Q^{sum} (t, q_{eq}, 0_{n \times 1})}{\partial \overset{\cdot}{q}} \overset{\cdot}{x} + \\ ε_{Q} (t, x, \overset{\cdot}{x}) \sqrt{{‖ x ‖}^{2} + {‖ \overset{\cdot}{x} ‖}^{2}} \end{matrix}

where $ε_{Q} : [0, \infty [\times U_{q_{eq}} \times R^{n \times 1} \to R^{1 \times n}$ is continuous at $(x, \overset{\cdot}{x}) = (0_{n \times 1}, 0_{n \times 1})$ and $ε_{Q} (t, 0_{n \times 1}, 0_{n \times 1}) = 0_{1 \times n}$ . Thus

\begin{matrix} {\ddot{x}}^{T} M_{2} (t, q_{eq}) - x^{T} \frac{\partial Q^{sum} (t, q_{eq}, 0_{n \times 1})}{\partial q} - {\overset{\cdot}{x}}^{T} \frac{\partial Q^{sum} (t, q_{eq}, 0_{n \times 1})}{\partial \overset{\cdot}{q}} = Q^{sum} (t, q_{eq}, 0_{n \times 1}) + \\ ε_{Q} (x, \overset{\cdot}{x}) \sqrt{{‖ x ‖}^{2} + {‖ \overset{\cdot}{x} ‖}^{2}} - {\ddot{x}}^{T} (\frac{\partial M_{2} (t, q_{eq})}{\partial q} x + ε_{M} (x) ‖ x ‖) \end{matrix}

and from this we extract the linearized equations of motion:

{\ddot{x}}^{T} M_{2} (t, q_{eq}) - x^{T} \frac{\partial Q^{sum} (t, q_{eq}, 0_{n \times 1})}{\partial q} - {\overset{\cdot}{x}}^{T} \frac{\partial Q^{sum} (t, q_{eq}, 0_{n \times 1})}{\partial \overset{\cdot}{q}} = Q^{ec} (t, q_{eq}, 0_{n \times 1})

(28)

where we have used the fact that, according to (26):

Q^{sum} (t, q_{eq}, 0_{n \times 1}) = Q^{cif} (t, q_{eq}, 0_{n \times 1}) + Q^{i} (t, q_{eq}, 0_{n \times 1}) + Q^{b} (t, q_{eq}) + Q^{ec} (t, q_{eq}, 0_{n \times 1}) = Q^{ec} (t, q_{eq}, 0_{n \times 1})

We have

\frac{\partial Q^{sum} (t, q_{eq}, 0_{n \times 1})}{\partial q} = {(\frac{\partial Q^{cif}}{\partial q})}_{eq} + {(\frac{\partial Q^{i}}{\partial q})}_{eq} + {(\frac{\partial Q^{b}}{\partial q})}_{eq} + {(\frac{\partial Q^{ec}}{\partial q})}_{eq}

(29)

\frac{\partial Q^{sum} (t, q_{eq}, 0_{n \times 1})}{\partial \overset{\cdot}{q}} = {(\frac{\partial Q^{cif}}{\partial \overset{\cdot}{q}})}_{eq} + {(\frac{\partial Q^{i}}{\partial \overset{\cdot}{q}})}_{eq} + {(\frac{\partial Q^{b}}{\partial \overset{\cdot}{q}})}_{eq} + {(\frac{\partial Q^{ec}}{\partial \overset{\cdot}{q}})}_{eq}

where we have introduced the simplifying notation

{(\frac{\partial Q^{cif}}{\partial q})}_{eq} = \frac{\partial Q^{cif} (t, q_{eq}, 0_{1 \times n})}{\partial q}, {(\frac{\partial Q^{cif}}{\partial \overset{\cdot}{q}})}_{eq} = \frac{\partial Q^{cif} (t, q_{eq}, 0_{1 \times n})}{\partial \overset{\cdot}{q}}

From (9) it follows that

{(\frac{\partial Q^{c i f}}{\partial q})}_{e q} = - \frac{1}{2} {(\frac{\partial^{2} M_{1}}{\partial q \partial t})}_{e q} + \frac{1}{2} {(\frac{\partial^{2} M_{0}}{\partial q^{2}})}_{e q}, {(\frac{\partial Q^{c i f}}{\partial \dot{q}})}_{e q} = {(\frac{\partial M_{2}}{\partial t})}_{e q} + s k e w {(\frac{\partial M_{1}}{\partial q})}_{e q}

We now have the following.

Proposition 5.2 The linearized equation of motion (28) may be written

{\ddot{x}}^{T} M + {\overset{\cdot}{x}}^{T} I + x^{T} J = F

(30)

where

M = M (t) = M_{2} (t, q_{e})

\begin{matrix} I & = I (t) = - ({(\frac{\partial M_{2}}{\partial t})}_{eq} + skew {(\frac{\partial M_{1}}{\partial q})}_{eq} + {(\frac{\partial Q^{i}}{\partial \overset{\cdot}{q}})}_{eq} + {(\frac{\partial Q^{b}}{\partial \overset{\cdot}{q}})}_{eq} + {(\frac{\partial Q^{ec}}{\partial \overset{\cdot}{q}})}_{eq}) \\ J & = J (t) = - (- \frac{1}{2} {(\frac{\partial^{2} M_{1}}{\partial q \partial t})}_{eq} + \frac{1}{2} {(\frac{\partial^{2} M_{0}}{\partial q^{2}})}_{eq} + {(\frac{\partial Q^{i}}{\partial q})}_{eq} + {(\frac{\partial Q^{b}}{\partial q})}_{eq} + {(\frac{\partial Q^{ec}}{\partial q})}_{eq}) \end{matrix}

(31)

and

F = Q^{ec} (t, q_{eq}, 0_{n \times 1})

Corollary 5.1 If $Q^{b}$ and $Q^{i}$ are given by the assumptions (14) and (15), respectively, then

\begin{matrix} I (t) & = D_{e, 2} (t, q_{eq}) - {(\frac{\partial M_{2}}{\partial t})}_{eq} - skew {(\frac{\partial M_{1}}{\partial q})}_{eq} - {(\frac{\partial Q^{ec}}{\partial \overset{\cdot}{q}})}_{eq} \\ J (t) & = {(\frac{\partial^{2} V^{*}}{\partial q^{2}})}_{eq} + \frac{1}{2} {(\frac{\partial D_{e, 1}}{\partial q})}_{eq} + \frac{1}{2} {(\frac{\partial^{2} M_{1}}{\partial q \partial t})}_{eq} - {(\frac{\partial Q^{ec}}{\partial q})}_{eq} \end{matrix}

(32)

Proof: We have, according to (19):

{(\frac{\partial Q^{i}}{\partial \overset{\cdot}{q}})}_{eq} + {(\frac{\partial Q^{b}}{\partial \overset{\cdot}{q}})}_{eq} = - {(\frac{\partial^{2} D_{e}}{\partial {\overset{\cdot}{q}}^{2}})}_{eq} = - D_{e, 2} (t, q_{e})

since $\frac{\partial^{2} V_{e}}{\partial \overset{\cdot}{q} \partial q} = \frac{\partial^{2} V_{b}}{\partial \overset{\cdot}{q} \partial q} = 0_{n \times n}$ . Furthermore:

{(\frac{\partial Q^{i}}{\partial q})}_{eq} + {(\frac{\partial Q^{b}}{\partial q})}_{eq} = - {(\frac{\partial^{2} V_{e}}{\partial q^{2}})}_{eq} - {(\frac{\partial^{2} D_{e}}{\partial q \partial \overset{\cdot}{q}})}_{eq} - {(\frac{\partial^{2} V_{g}}{\partial q^{2}})}_{eq}

But $\frac{\partial D_{e}}{\partial \overset{\cdot}{q}} = \frac{1}{2} D_{e, 1} + {\overset{\cdot}{q}}^{T} D_{e, 2}$ and then

\frac{\partial^{2} D_{e}}{\partial q \partial \overset{\cdot}{q}} = \frac{1}{2} \frac{\partial D_{e, 1}}{\partial q} + \frac{D_{e, 2}}{\partial q} \overset{\cdot}{q} \Rightarrow {(\frac{\partial^{2} D_{e}}{\partial q \partial \overset{\cdot}{q}})}_{eq} = \frac{1}{2} {(\frac{\partial D_{e, 1}}{\partial q})}_{eq}

and this proves the corollary. □

The matrix $I$ in (32) may be written $I (t) = C (t) + G^{T} (t)$ where

C = sym (I) = D_{e, 2} (t, q_{eq}) - {(\frac{\partial M_{2}}{\partial t})}_{eq} - sym {(\frac{\partial Q^{ec}}{\partial \overset{\cdot}{q}})}_{eq}

(33)

G^{T} = skew (I) = - skew {(\frac{\partial M_{1}}{\partial q})}_{eq} - skew {(\frac{\partial Q^{ec}}{\partial \overset{\cdot}{q}})}_{eq}

Here $C$ is called the damping matrix and $G$ is called the gyroscopic matrix. The matrix $J$ may be written $J (t) = K (t) + H (t)^{T}$ where

\begin{matrix} K & = sym (J) = {(\frac{\partial^{2} V^{*}}{\partial q^{2}})}_{eq} + \frac{1}{2} sym {(\frac{\partial D_{e, 1}}{\partial q})}_{eq} + \frac{1}{2} sym {(\frac{\partial^{2} M_{1}}{\partial q \partial t})}_{eq} - sym {(\frac{\partial Q^{ec}}{\partial q})}_{eq} \\ H^{T} & = skew (J) = \frac{1}{2} skew {(\frac{\partial D_{e, 1}}{\partial q})}_{eq} + \frac{1}{2} skew {(\frac{\partial^{2} M_{1}}{\partial q \partial t})}_{eq} - skew {(\frac{\partial Q^{ec}}{\partial q})}_{eq} \end{matrix}

(34)

$K$ is called the stiffness matrix and $H$ is called the circulatory matrix. Note that

{(\frac{\partial^{2} V^{*}}{\partial q^{2}})}_{eq} = {(\frac{\partial^{2} V_{e}}{\partial q^{2}})}_{eq} + {(\frac{\partial^{2} V_{g}}{\partial q^{2}})}_{eq} - \frac{1}{2} {(\frac{\partial^{2} M_{0}}{\partial q^{2}})}_{eq}

where

\frac{\partial^{2} V_{e}}{\partial q^{2}} (t, q) = \sum_{α = 1}^{N} \int_{B_{0}^{α}} \frac{\partial^{2} E_{q}^{α}}{\partial q^{2}} \cdot C^{α} [E_{q}^{α}] + \frac{\partial E_{q}^{α}}{\partial q} \cdot C^{α} [\frac{\partial E_{q}^{α}}{\partial q}] dv (X), \frac{\partial^{2} V_{g}}{\partial q^{2}} = \sum_{α = 1}^{N} \frac{\partial^{2} x_{c}^{α}}{\partial q^{2}} \cdot g m^{α}

Corollary 5.2 If

{(\frac{\partial D_{e, 1}}{\partial q})}_{eq}, {(\frac{\partial^{2} M_{1}}{\partial q \partial t})}_{eq}, {(\frac{\partial Q^{ec}}{\partial q})}_{eq} \in Sym (R^{n \times n})

then $H = 0_{n \times n}$ and

K = {(\frac{\partial^{2} V^{*}}{\partial q^{2}})}_{eq} + \frac{1}{2} {(\frac{\partial D_{e, 1}}{\partial q})}_{eq} + \frac{1}{2} {(\frac{\partial^{2} M_{1}}{\partial q \partial t})}_{eq} - {(\frac{\partial Q^{ec}}{\partial q})}_{eq}

Proof: This is a direct consequence of (34). □

The coordinate system is called quasi-scleronomic if $M_{0}$ , $M_{1}$ , $M_{2}$ , $D_{e, 1}$ , $D_{e, 2}$ and $V$ are explicit functions of $q$ only, and not of $t$ . The external force is called quasi-scleronomic if

Q^{ec} (t, q, \overset{\cdot}{q}) = Q_{1}^{ec} (t) + Q_{2}^{ec} (q, \overset{\cdot}{q})

The external force then consists of two parts, one prescribed time-dependent part $Q_{1}^{ec}$ and one part $Q_{2}^{ec}$ depending on the state of the multibody. This second part may include so-called follower forces.

Remark 5.1: Traditionally, a coordinate system $x = χ_{q}^{α} (t, q; X)$ is called scleronomic if $\partial_{t} χ_{q} (t, q; X) = 0, (t, q; X) \in [0, \infty [\times Ω \times B_{0}^{α}$ . In this case, $M_{0} = 0$ , $M_{1} = 0_{1 \times n}$ , $\partial_{t} M_{2} = 0_{n \times n}$ , $\partial_{t} E_{q}^{α} = 0$ and, consequently, $D_{0} = 0$ , $D_{1} = 0_{1 \times n}$ .■

Corollary 5.3 If the coordinate system and the external force are quasi-scleronomic then

\begin{matrix} M & = M_{2} (q_{e}), C = D_{e, 2} (q_{eq}) - sym {(\frac{\partial Q_{2}^{ec}}{\partial \overset{\cdot}{q}})}_{eq}, G^{T} = - skew {(\frac{\partial M_{1}}{\partial q})}_{eq} - skew {(\frac{\partial Q_{2}^{ec}}{\partial \overset{\cdot}{q}})}_{eq} \\ K & = {(\frac{\partial^{2} V^{*}}{\partial q^{2}})}_{eq} + \frac{1}{2} sym {(\frac{\partial D_{e, 1}}{\partial q})}_{eq} - sym {(\frac{\partial Q_{2}^{ec}}{\partial q})}_{eq}, H^{T} = \frac{1}{2} skew {(\frac{\partial D_{e, 1}}{\partial q})}_{eq} - skew {(\frac{\partial Q_{2}^{ec}}{\partial q})}_{eq} \end{matrix}

(35)

and

F (t) = Q_{1}^{ec} (t) + Q_{2}^{ec} (q_{eq}, 0_{n \times 1})

where now all matrices $M$ , $C$ , $G$ , $K$ and $H$ are constant matrices.

Proof: The expressions in (35) are direct consequences of (33) and (34). □

Corollary 5.4 If the coordinate system and the external force are quasi-scleronomic and

\frac{\partial E_{q}^{α}}{\partial t} = 0, Q_{2}^{ec} = 0_{1 \times n}

then

K = {(\frac{\partial^{2} V^{*}}{\partial q^{2}})}_{eq}, H = 0_{n \times n}

(36)

Proof: From (20) it follows that

\frac{\partial D_{e, 1}}{\partial q} = 2 \sum_{α = 1}^{N} \int_{B_{0}^{α}} (\frac{\partial^{2} E_{q}^{α}}{\partial q \partial t} \cdot D^{α} [\frac{\partial E_{q}^{α}}{\partial q}] + \frac{\partial E_{q}^{α}}{\partial t} \cdot D^{α} [\frac{\partial^{2} E_{q}^{α}}{\partial q^{2}}] dv (X) = 0_{n \times n}

This, together with (35)_4,5, proves the corollary. □

Introduce the linearized mechanical energy

E = \frac{1}{2} {\overset{\cdot}{x}}^{T} M \overset{\cdot}{x} + \frac{1}{2} x^{T} Kx

(37)

Proposition 5.3 If the coordinate system and the external force are quasi-scleronomic, then the power theorem corresponding to the linearized equations reads

\overset{\cdot}{E} = F \overset{\cdot}{x} - {\overset{\cdot}{x}}^{T} C \overset{\cdot}{x} - {\overset{\cdot}{x}}^{T} Hx

(38)

Proof: Considering the power theorem (24) we have the following Taylor expansion:

\begin{matrix} \frac{1}{2} {\overset{\cdot}{q}}^{T} M_{2} (t, q) \overset{\cdot}{q} = \frac{1}{2} {\overset{\cdot}{x}}^{T} M_{2} (t, q_{eq} + x) \overset{\cdot}{x} = \frac{1}{2} {\overset{\cdot}{x}}^{T} (M_{2} (t, q_{eq}) + \frac{\partial M_{2} (t, q_{eq})}{\partial q} x + ε_{M} (t, x) ‖ x ‖) \overset{\cdot}{x} = \\ \frac{1}{2} {\overset{\cdot}{x}}^{T} M_{2} (t, q_{eq}) \overset{\cdot}{x} + \frac{1}{2} (\frac{\overset{\cdot}{x}}{‖ \overset{\cdot}{x} ‖})^{T} (\frac{\partial M_{2} (t, q_{eq})}{\partial q} x + ε_{M} (t, x) ‖ x ‖) \frac{\overset{\cdot}{x}}{‖ \overset{\cdot}{x} ‖} ‖ \overset{\cdot}{x} ‖^{2} = \\ \frac{1}{2} {\overset{\cdot}{x}}^{T} M \overset{\cdot}{x} + ε_{K} (t, x, \overset{\cdot}{x}) ‖ \overset{\cdot}{x} ‖^{2}, \overset{\cdot}{x} \neq 0_{n \times 1} \end{matrix}

where the mapping $ε_{K} : [0, \infty [\times U_{q_{eq}} \times (R^{n \times 1} \ 0_{n \times 1}) \to R$ is defined by

ε_{K} (t, x, \overset{\cdot}{x}) = \frac{1}{2} (\frac{\overset{\cdot}{x}}{‖ \overset{\cdot}{x} ‖})^{T} (\frac{\partial M_{2} (t, q_{eq})}{\partial q} x + ε_{M} (t, x) ‖ x ‖) \frac{\overset{\cdot}{x}}{‖ \overset{\cdot}{x} ‖}

It satisfies $ε_{K} (t, x, λ \overset{\cdot}{x}) = ε_{K} (t, x, \overset{\cdot}{x}), λ \in R, λ \neq 0$ and the mapping is continuous at $x = 0_{n \times 1}$ and $ε_{K} (t, 0_{n \times 1}, \overset{\cdot}{x}) = 0$ . Furthermore

V^{*} (t, q) = V^{*} (t, q_{eq} + x) = V^{*} (t, q_{eq}) + x^{T} {(\frac{\partial V^{*}}{\partial q})}_{eq} + \frac{1}{2} x^{T} {(\frac{\partial^{2} V^{*}}{\partial q^{2}})}_{eq} x + ε_{V^{*}} (t, x) ‖ x ‖^{2}

where $ε_{V^{*}} : [0, \infty [\times U_{q_{eq}} \to R$ is continuous at $x = 0_{n \times 1}$ and $ε_{V^{*}} (t, 0_{n \times 1}) = 0$ and

D_{e, 1} (t, q) = D_{e, 1} (t, q_{eq}) + x^{T} {(\frac{\partial D_{e, 1}}{\partial q})}_{eq} + ε_{D_{e, 1}} (t, x) ‖ x ‖

D_{e, 2} (t, q) = D_{e, 2} (t, q_{eq}) + x^{T} {(\frac{\partial D_{e, 2}}{\partial q})}_{eq} + ε_{D_{e, 2}} (t, x) ‖ x ‖

where the mappings $ε_{D_{e, 1}} : [0, \infty [\times U_{q_{eq}} \to R$ and $ε_{D_{e, 2}} : [0, \infty [\times U_{q_{eq}} \to R^{1 \times n}$ are continuous at $x = 0_{n \times 1}$ , $ε_{D_{e, 1}} (t, 0_{n \times 1}) = 0$ and $ε_{D_{e, 2}} (t, 0_{n \times 1}) = 0_{1 \times n}$ . Finally

Q_{2}^{ec} (q, \overset{\cdot}{q}) = Q_{2}^{ec} (q_{eq}, 0_{n \times 1}) + x^{T} {(\frac{\partial Q_{2}^{ec}}{\partial q})}_{eq} + {\overset{\cdot}{x}}^{T} {(\frac{\partial Q_{2}^{ec}}{\partial \overset{\cdot}{q}})}_{eq} + ε_{Q_{2}^{ec}} (x, \overset{\cdot}{x}) \sqrt{{‖ x ‖}^{2} + {‖ \overset{\cdot}{x} ‖}^{2}}

and then

Q^{ec} (t, q, \overset{\cdot}{q}) = Q_{1}^{ec} (t) + Q_{2}^{ec} (q, \overset{\cdot}{q}) = F + x^{T} {(\frac{\partial Q_{2}^{ec}}{\partial q})}_{eq} + {\overset{\cdot}{x}}^{T} {(\frac{\partial Q_{2}^{ec}}{\partial \overset{\cdot}{q}})}_{eq} + ε_{Q_{2}^{ec}} (x, \overset{\cdot}{x}) \sqrt{{‖ x ‖}^{2} + {‖ \overset{\cdot}{x} ‖}^{2}}

Inserting this into (24), and omitting higher order terms, one obtains

\begin{matrix} \frac{d}{dt} (\frac{1}{2} {\overset{\cdot}{q}}^{T} M_{2} \overset{\cdot}{q} + V^{*} (q_{eq}) + {(\frac{\partial V^{*}}{\partial q})}_{eq} x + \frac{1}{2} x^{T} {(\frac{\partial^{2} V^{*}}{\partial q^{2}})}_{eq} x) = - \frac{1}{2} (D_{e, 1} (q_{eq}) + x^{T} {(\frac{\partial D_{e, 1}}{\partial q})}_{eq}) \overset{\cdot}{x} - \\ {\overset{\cdot}{x}}^{T} (D_{e, 2} (q_{eq}) + x^{T} {(\frac{\partial D_{e, 2}}{\partial q})}_{eq}) \overset{\cdot}{x} + (F (t) + x^{T} {(\frac{\partial Q_{2}^{ec}}{\partial q})}_{eq} + {\overset{\cdot}{x}}^{T} {(\frac{\partial Q_{2}^{ec}}{\partial \overset{\cdot}{q}})}_{eq}) \overset{\cdot}{x} \end{matrix}

where we have used the fact that, in the present quasi-scleronomic case, $\frac{\partial L}{\partial t} = 0$ . The previous equation is then equivalent to

\begin{matrix} \frac{d}{dt} (\frac{1}{2} {\overset{\cdot}{q}}^{T} M_{2} \overset{\cdot}{q} + \frac{1}{2} x^{T} {(\frac{\partial^{2} V^{*}}{\partial q^{2}})}_{eq} x) = - ({(\frac{\partial V^{*}}{\partial q})}_{eq} x + \frac{1}{2} D_{e, 1} (q_{eq})) \overset{\cdot}{x} - \frac{1}{2} x^{T} {(\frac{\partial D_{e, 1}}{\partial q})}_{eq}) \overset{\cdot}{x} - \\ {\overset{\cdot}{x}}^{T} (D_{e, 2} (q_{eq}) + x^{T} {(\frac{\partial D_{e, 2}}{\partial q})}_{eq}) \overset{\cdot}{x} + (F + x^{T} (\frac{\partial Q_{2}^{ec}}{\partial q}) + {\overset{\cdot}{x}}^{T} (\frac{\partial Q_{2}^{ec}}{\partial \overset{\cdot}{q}})) \overset{\cdot}{x} \end{matrix}

(39)

However, according to the equilibrium condition (27), we have

{(\frac{\partial V^{*}}{\partial q})}_{eq} x + \frac{1}{2} D_{e, 1} (q_{eq}) = 0_{1 \times n}

and then (39) may be written

\begin{matrix} \frac{d}{dt} (\frac{1}{2} {\overset{\cdot}{q}}^{T} M_{2} \overset{\cdot}{q} + \frac{1}{2} x^{T} ({(\frac{\partial^{2} V^{*}}{\partial q^{2}})}_{eq} + \frac{1}{2} sym {(\frac{\partial D_{e, 1}}{\partial q})}_{eq} - sym {(\frac{\partial Q^{ec}}{\partial q})}_{eq}) x) = \\ - {\overset{\cdot}{x}}^{T} (D_{e, 2} (q_{eq}) - sym {(\frac{\partial Q_{2}^{ec}}{\partial \overset{\cdot}{q}})}_{eq}) \overset{\cdot}{x} - x^{T} (\frac{1}{2} skew {(\frac{\partial D_{e, 1}}{\partial q})}_{eq} - skew {(\frac{\partial Q_{2}^{ec}}{\partial q})}_{eq}) \overset{\cdot}{x} + F \overset{\cdot}{x} \end{matrix}

and this proves the proposition. □

Remark 5.2: Multiplying equation (30) by $\overset{\cdot}{x}$ from the right one obtains

{\ddot{x}}^{T} M \overset{\cdot}{x} + {\overset{\cdot}{x}}^{T} I \overset{\cdot}{x} + x^{T} J \overset{\cdot}{x} = {\ddot{x}}^{T} M \overset{\cdot}{x} + {\overset{\cdot}{x}}^{T} (C + G^{T}) \overset{\cdot}{x} + x^{T} (K + H^{T}) \overset{\cdot}{x} = F \overset{\cdot}{x}

But ${\overset{\cdot}{x}}^{T} G^{T} \overset{\cdot}{x} = 0$ and $\overset{\cdot}{E} = {\ddot{x}}^{T} M \overset{\cdot}{x} + x^{T} K \overset{\cdot}{x}$ where $E$ is the linearized mechanical energy defined in (37). This inserted into the equation above results in the power theorem according to (38). If we accept complex valued solutions $z = z (t)$ to (30) then the mechanical energy is defined by

E = \frac{1}{2} {\overset{\cdot}{z}}^{H} M \overset{\cdot}{z} + \frac{1}{2} z^{H} K z

where $z^{H}$ denotes the Hermitian transpose of $z$ (see Appendix A.1). Now assume that $F (t) = 0_{n \times 1}$ and $z (t) = w e^{st}$ , where $w \in C^{n \times 1}$ and $s \in C$ are constants. Then, since $H \in Skew (R^{n \times n})$ :

\overset{\cdot}{E} = - {\overset{\cdot}{z}}^{H} C \overset{\cdot}{z} - {\overset{\cdot}{z}}^{H} H z = - (s^{*} w^{H} e^{s^{*} t}) C (sw e^{st}) - (s^{*} w^{H} e^{s^{*} t}) Hw e^{st} = - {| s |}^{2} e^{2 Re (s) t} w^{H} Cw < 0

if $C$ is positive definite. Thus, the linearized mechanical energy is decreasing. More on this will be discussed in the forthcoming sections.■

6. The linearized equation

If Equation (30) is transposed one obtains the linearized equations of motion in the more familiar format:

M \ddot{x} + A \overset{\cdot}{x} + Bx = f

(40)

where $A = I^{T} = C + G$ , $B = J^{T} = K + H$ and $f = F^{T}$ . We have $M, C, K \in Sym (R^{n \times n})$ and $G, H \in Skew (R^{n \times n})$ . Along with (40) one prescribes the initial conditions

x (0) = x_{0}, \overset{\cdot}{x} (0) = {\overset{\cdot}{x}}_{0}

(41)

where $x_{0}, {\overset{\cdot}{x}}_{0} \in R^{n \times 1}$ are given. A function $x = x (t)$ satisfying (40) and (41) will, under certain conditions, represent a vibrational motion of the multibody in the neighbourhood of the equilibrium configuration. If the external force $f$ is identically zero the equation represents a free vibration of the system. Obviously, in this case and with $x_{0} = 0_{n \times 1}, {\overset{\cdot}{x}}_{0} = 0_{n \times 1}$ , $x (t) = 0_{n \times 1}, t \in [0, \infty [$ is a solution to the initial value problem representing the equilibrium solution.

Example 6.1: A particle P with mass $m$ , which may slide without friction on a horizontal turn-table, is connected to a fixed point A on the table by an ideal linear elastic spring with spring constant $k$ and spring natural (unstressed) length $r_{n}$ (see Figure 1). The table is rotating around a fixed vertical axis through O with the prescribed constant angular velocity $Ω$ and its moment of inertia with respect to the axis is denoted by $I$ . Distance $OA = a$ .

Figure 1.

Two-dimensional multibody.

Introduce configuration coordinates $r, φ$ according to Figure 1. The kinetic energy of the system

T = \frac{1}{2} (M_{0} + M_{1} \overset{\cdot}{q} + {\overset{\cdot}{q}}^{T} M_{2} \overset{\cdot}{q})

where $q = {(\begin{matrix} r & φ \end{matrix})}^{T}, r \in R^{+}, φ \in R$ and

M_{0} = (m (a^{2} + r^{2} + 2 ar \cos φ) + I) Ω^{2}, M_{1} = 2 m Ω (\begin{matrix} a \sin φ & r^{2} + ar \cos φ \end{matrix}), M_{2} = m (\begin{matrix} 1 & 0 \\ 0 & r^{2} \end{matrix})

The elastic potential energy is given by

V_{e} = \frac{1}{2} k (r - r_{n})^{2}

The coordinate system used is thus quasi-scleronomic. We have $Q^{c} = 0_{1 \times 2}$

\begin{matrix} Q^{cif} & = - ({\overset{\cdot}{q}}^{T} (\frac{\partial M_{2}}{\partial q} - \frac{1}{2} {(\frac{\partial M_{2}}{\partial q})}^{t}) \overset{\cdot}{q} + {\overset{\cdot}{q}}^{T} skew (\frac{\partial M_{1}}{\partial q}) - \frac{1}{2} \frac{\partial M_{0}}{\partial q})) = \\ - m ((\begin{matrix} - r {\overset{\cdot}{φ}}^{2} & 2 r \overset{\cdot}{r} \overset{\cdot}{φ} \end{matrix}) + 2 Ω (\begin{matrix} - r \overset{\cdot}{φ} & r \overset{\cdot}{r} \end{matrix}) - Ω^{2} (\begin{matrix} r + a \cos φ & - ar \sin φ \end{matrix})) = \\ - m (\begin{matrix} - r {\overset{\cdot}{φ}}^{2} - 2 Ω r \overset{\cdot}{φ} - Ω^{2} (r + a \cos φ) & 2 r \overset{\cdot}{r} \overset{\cdot}{φ} + 2 Ω r \overset{\cdot}{r} + Ω^{2} ar \sin φ \end{matrix}) \end{matrix}

and

Q^{i} = - \frac{\partial V_{e}}{\partial q} = - (\begin{matrix} \frac{\partial V_{e}}{\partial r} & \frac{\partial V_{e}}{\partial φ} \end{matrix}) = - (\begin{matrix} k (r - r_{n}) & 0 \end{matrix})

The equations of motion (7) and (22) read

\begin{matrix} (\begin{matrix} \ddot{r} & \ddot{φ} \end{matrix}) m (\begin{matrix} 1 & 0 \\ 0 & r^{2} \end{matrix}) = \\ - m (\begin{matrix} - r {\overset{\cdot}{φ}}^{2} - 2 Ω r \overset{\cdot}{φ} - Ω^{2} (r + a \cos φ) - \frac{k}{m} (r - r_{n}) & 2 r \overset{\cdot}{r} \overset{\cdot}{φ} + 2 Ω r \overset{\cdot}{r} + Ω^{2} ar \sin φ \end{matrix}) \end{matrix}

The modified potential energy is given by

V^{*} = V_{e} - \frac{1}{2} M_{0} = \frac{1}{2} k (r - r_{n})^{2} - \frac{1}{2} (m (a^{2} + r^{2} + 2 ar \cos φ) + I) Ω^{2}

and

\frac{\partial V^{*}}{\partial q} = (\begin{matrix} \frac{\partial V^{*}}{\partial r} & \frac{\partial V^{*}}{\partial φ} \end{matrix}) = (\begin{matrix} k (r - r_{n}) + m Ω^{2} (r + a \cos φ) & m Ω^{2} (- ar \sin φ) \end{matrix})

The equilibrium configuration $q_{eq} = {(\begin{matrix} r_{eq} & φ_{eq} \end{matrix})}^{T}$ is determined by Equation (27):

{(\frac{\partial V^{*}}{\partial q})}_{eq} = (\begin{matrix} k (r_{eq} - r_{n}) + m Ω^{2} (r_{eq} + a \cos φ_{eq}) & m Ω^{2} (- a r_{eq} \sin φ_{eq}) \end{matrix}) = 0_{1 \times 2}

which is equivalent to

{\begin{matrix} k (r_{eq} - r_{n}) + m Ω^{2} (r_{eq} + a \cos φ_{eq}) = 0 \\ r_{eq} \sin φ_{eq} = 0 \end{matrix}

with a solution: $r_{eq} = \frac{k r_{n} + m Ω^{2} a}{k - m Ω^{2}}, φ_{eq} = 0 if k > m Ω^{2}$ . This corresponds to a stable equilibrium configuration.

There may be another equilibrium configuration corresponding to $φ_{eq} = π$ , but this is then unstable. We introduce the configuration coordinate $x = q - q_{eq} = {(\begin{matrix} r - r_{eq} & φ - φ_{eq} \end{matrix})}^{T}$ representing the deviation from the equilibrium solution. The linearized equations of motion are given by (40) where

M = M_{2} (q_{e}) = m (\begin{matrix} 1 & 0 \\ 0 & r_{eq}^{2} \end{matrix}), C = 0_{2 \times 2}, G = - skew {(\frac{\partial M_{1}}{\partial q})}_{eq} = 2 m Ω r_{eq} (\begin{matrix} 0 & 1 \\ - 1 & 0 \end{matrix})

K = {(\frac{\partial^{2} V^{*}}{\partial q^{2}})}_{eq} = m (\begin{matrix} m Ω^{2} a r_{eq} & 0 \\ 0 & k - m Ω^{2} \end{matrix}), H = 0_{n \times n}, f = 0_{n \times 1}

In this case we have a non-zero gyroscopic matrix.■

Example 6.2: The planar double pendulum consists of two homogeneous, rigid and slender bars OA and AB with masses $m_{1}$ and $m_{2}$ and lengths $l_{1}$ and $l_{2}$ , respectively. The bars are connected by an ideal revolute joint at point A and bar OA is connected, by an ideal revolute joint, to the fixed point O. The configuration space in this case is the two-dimensional torus. Configuration coordinates $q = {(\begin{matrix} θ & ϕ \end{matrix})}^{T}$ are introduced according to Figure 2. The pendulum is, at the joints at O and A, equipped with elastic torsional springs, with spring constants $κ_{1}$ and $κ_{2}$ , respectively. The pendulum is subjected to a follower point force $f_{0}$ and a point couple $M_{0}$ applied at the end-point B of bar AB and defined according to

f_{0} = f_{0} (ϕ) = (i \sin ϕ - j \cos ϕ) f_{0}, M_{0} = M_{0} (t) = k M_{0} (t)

where $f_{0}$ is a non-zero constant and $M_{0} = M_{0} (t)$ is a prescribed function, $(i j k)$ being an orthonormal basis according to Figure 2. The pendulum may move in a vertical plane. The kinetic energy of the pendulum is given by

T = T (q, \overset{\cdot}{q}) = \frac{1}{2} {\overset{\cdot}{q}}^{T} M_{2} \overset{\cdot}{q}

since $M_{0} = 0$ and $M_{1} = 0_{1 \times 2}$ due to the scleronomic coordinate system.

Figure 2.

The planar double pendulum subjected to a follower force and couple.

The potential energy of the pendulum may be written $V (θ, ϕ) = V_{e} (θ, ϕ) + V_{g} (θ, ϕ)$ where

V_{e} (θ, ϕ) = \frac{1}{2} κ_{1} θ^{2} + \frac{1}{2} κ_{2} (ϕ - θ)^{2}, V_{g} (θ, ϕ) = - (\frac{l_{1}}{2} m_{1} \cos θ + m_{2} (l_{1} \cos θ + \frac{l_{2}}{2} \cos ϕ)) g

are the elastic and gravitational potential energies, respectively. The equations of motion for the pendulum read, according to (7) and (22):

{\ddot{q}}^{T} M_{2} = Q^{cif} + Q^{i} + Q^{ec} + Q^{b}

where in this case

M_{2} = M_{2} (θ, ϕ) = (\begin{matrix} \frac{m_{1} l_{1}^{2}}{3} + m_{2} l_{1}^{2} & \frac{m_{2} l_{1} l_{2}}{2} \cos (ϕ - θ) \\ \frac{m_{2} l_{1} l_{2}}{2} \cos (ϕ - θ) & \frac{m_{2} l_{2}^{2}}{3} \end{matrix})

and

Q^{cif} = - {\overset{\cdot}{q}}^{T} (\frac{\partial M_{2}}{\partial q} - \frac{1}{2} {(\frac{\partial M_{2}}{\partial q})}^{t}) \overset{\cdot}{q} = (\begin{matrix} \frac{m_{2} l_{1} l_{2}}{2} \sin (ϕ - θ) {\overset{\cdot}{ϕ}}^{2} & - \frac{m_{2} l_{1} l_{2}}{2} \sin (ϕ - θ) {\overset{\cdot}{θ}}^{2} \end{matrix})

Q^{i} = (\begin{matrix} θ & ϕ \end{matrix}) (\begin{matrix} - (κ_{1} + κ_{2}) & κ_{2} \\ κ_{2} & - κ_{2} \end{matrix}), Q^{b} = (\begin{matrix} \sin θ & \sin ϕ \end{matrix}) (\begin{matrix} - (\frac{1}{2} m_{1} + m_{2}) l_{1} g & 0 \\ 0 & - \frac{l_{2}}{2} m_{2} g \end{matrix})

Q^{ec} = Q_{1}^{ec} (t) + Q_{2}^{ec} (θ, ϕ)

where

Q_{1}^{ec} (t) = (\begin{matrix} 0 & M_{0} (t) \end{matrix}), Q_{2}^{ec} (θ, ϕ) = (\begin{matrix} l_{1} f_{0} \sin (ϕ - θ) & 0 \end{matrix})

At the equilibrium configuration, given by $θ_{eq} = ϕ_{eq} = 0$ , one has $C = G = 0_{2 \times 2}$ , and

K = {(\frac{\partial^{2} V}{\partial q^{2}})}_{eq} - sym {(\frac{\partial Q_{2}^{ec}}{\partial q})}_{eq} = (\begin{matrix} (κ_{1} + κ_{2}) + l_{1} (\frac{m_{1}}{2} + m_{2}) g + l_{1} f_{0} & - κ_{2} - \frac{l_{1}}{2} f_{0} \\ - κ_{2} - \frac{l_{1}}{2} f_{0} & κ_{2} + \frac{l_{2}}{2} m_{2} g \end{matrix})

H^{T} = - skew {(\frac{\partial Q_{2}^{ec}}{\partial q})}_{eq} = (\begin{matrix} 0 & - 1 \\ 1 & 0 \end{matrix}) \frac{l_{1}}{2} f_{0}, F (t) = Q_{1}^{ec} (t) + Q_{2}^{ec} (0, 0) = (\begin{matrix} 0 & M_{0} (t) \end{matrix})

In this case we thus have a non-zero circulatory matrix. There are three more equilibrium configurations corresponding to $θ_{eq} = 0$ , $ϕ_{eq} = π$ , $θ_{eq} = π$ , $ϕ_{eq} = π$ and $θ_{eq} = π$ , $ϕ_{eq} = 0$ , but these are unstable.■

The matrix polynomial

S (s) = M s^{2} + As + B, s \in C

(42)

is called the dynamic stiffness. The differential equation (40) may then be written ( $d_{t} x = \frac{dx}{dt}$ )

S (d_{t}) x = {Md}_{t}^{2} x + A d_{t} x + Bx = f

(43)

The characteristic polynomial $p = p (s)$ , associated with the dynamic stiffness, is defined by $p (s) = det S (s)$ . It has degree equal to $2 n$ and one may show that

p (s) = s^{2 n} det M + \dots + det B

(44)

where the intermediate terms have more complicated coefficients in terms of $M, A$ and $B$ .

Example 6.1 (continued): The dynamic stiffness matrix of the system is given by

S (s) = M s^{2} + (C + G) s + K + L = (\begin{matrix} {mr}_{eq}^{2} s^{2} + m Ω^{2} a r_{eq} & 2 m Ω r_{eq} s \\ - 2 m Ω r_{eq} s & m s^{2} + k - m Ω^{2} \end{matrix})

and the characteristic polynomial reads

p (s) = m^{2} r_{eq}^{2} s^{4} + m r_{eq} (k r_{eq} + m Ω^{2} (a + 3 r_{eq})) s^{2} + m Ω^{2} a r_{eq} (k - m Ω^{2})

■

We assume, from now on, that $M$ is positive definite. Then $det M \neq 0$ . The spectrum $spec (S)$ , of the dynamic stiffness, is defined by $spec (S) = {s \in C | p (s) = 0}$ . Note that $0 \in spec (S) \Leftrightarrow det B = 0$ . The characteristic polynomial may be factorized according to

p (s) = det M (s - s_{1}) (s - s_{2}) \dots (s - s_{2 n})

(45)

where $s_{i} \in spec (S), i = 1, \dots, 2 n$ . The characteristic polynomial has real coefficients and, consequently, $s \in spec (S)$ if and only if $s^{*} \in spec (S)$ . The roots of the characteristic equation $p (s) = 0$ may then be grouped into two sets $s_{1}, \dots, s_{r} \in R$ and $s_{r + 1}, s_{r + 1}^{*}, \dots, s_{r + c}, s_{r + c}^{*} \in C$ , where $r + 2 c = 2 n$ . Note that the number of real roots must be even and that some of the $s_{i}, i = 1, \dots, 2 n$ may be equal. The number of times the root $s_{i}$ appears in (45) is called its (algebraic) multiplicity and denoted $m_{i} \in N$ ( $1 \leq m_{i} \leq 2 n)$ . We may write

p (s) = det M (s - s_{1})^{m_{1}} (s - s_{2})^{m_{2}} \dots (s - s_{l})^{m_{l}}

(46)

where now all $s_{1}, \dots, s_{l}$ are distinct complex numbers and $\sum_{i = 1}^{l} m_{i} = 2 n$ .

The following well-known result will be frequently used in the forthcoming discussions.

Lemma 6.1 $s_{i} \in C$ is a root of the characteristic equation with multiplicity $m_{i}$ if and only if

p (s_{i}) = p' (s_{i}) = \dots = p^{(m_{i} - 1)} (s_{i}) = 0 and p^{(m_{i})} (s_{i}) \neq 0

Proof: $s_{i} \in C$ is a root with multiplicity $m_{i}$ if and only if $p (s) = det M (s - s_{i})^{m_{i}} r (s)$ , where $r$ is a polynomial with $\deg (r) = 2 n - m_{i}$ and $r (s_{i}) \neq 0$ . Now if $k \leq m_{i}$ then

p^{(k)} (s) = det M (s - s_{i})^{m_{1} - k} (\frac{m_{i}!}{(m_{i} - k)!} r (s) + (s - s_{i}) b_{k} (s))

where $b_{k}$ is a polynomial with $\deg (b_{k}) = 2 n - m_{i} - 1$ . Then $k < m_{i} \Rightarrow p^{(k)} (s_{i}) = 0$ and

p^{(m_{i})} (s_{i}) = det M m_{i}! r (s_{i}) \neq 0

This proves the lemma. □

With $k = 0, 1, \dots$ . we introduce the linear spaces:

$C^{k} ([0, \infty [, R^{n \times 1}) = {q : [0, \infty [\to R^{n \times 1} | q is k times continuously differentiable}$

The kernel (null-space) and range of the linear differential operator $S (d_{t})$ :

$C^{2} ([0, \infty [, R^{n \times 1}) \to$ $C^{0} ([0, \infty [, R^{n \times 1})$ are given by

\ker (S (d_{t})) = {x \in C^{2} ([0, \infty [, R^{n \times 1}) | S (d_{t}) x = 0_{n \times 1}}

(47)

r a n (S (d_{t})) = {g \in C^{0} ([0, \infty [, ℝ^{n \times 1}) | g = S (d_{t}) x, x \in C^{2} (T, ℝ^{n \times 1})}

For $f \in C^{0} ([0, \infty [, R^{n \times 1})$ we define $sol (S (d_{t}), f) = {x \in C^{2} ([0, \infty [, R^{n \times 1}) | S (d_{t}) x = f}$ and we then have the following trivial implications: $sol (S (d_{t}), f) \neq \emptyset \Rightarrow f \in ran (S (d_{t}))$ and $x \in sol (S (d_{t}), f)$ , $x_{h} \in \ker (S (d_{t})) \Rightarrow x + x_{h} \in sol (S (d_{t}), f)$ . Thus we have the following.

Proposition 6.1 Let $x_{p} \in sol (S (d_{t}), f)$ (a particulate solution). Then

sol (S (d_{t}), f) = {x \in C^{2} ([0, \infty [, R^{n \times 1}) | x = x_{p} + x_{h}, x_{h} \in \ker (S (d_{t}))}

Theorem 6.1 (Existence and uniqueness) The initial value problem

S (d_{t}) x = f, x (0) = x_{0}, \overset{\cdot}{x} (0) = {\overset{\cdot}{x}}_{0}

has a unique solution $x \in C^{2} (] - \infty, \infty [, R^{n \times 1})$ .

Proof: By introducing $z = {(\begin{matrix} x & y \end{matrix})}^{T} \in R^{2 n \times 1}$ , where $y = \overset{\cdot}{x}$ , Equation (43) may be written $\overset{\cdot}{z} = g (z)$ , where $g (z) = - A^{- 1} B z + A^{- 1} f$ and

A = (\begin{matrix} A & M \\ M & 0_{n \times n} \end{matrix}), B = (\begin{matrix} B & 0_{n \times n} \\ 0_{n \times n} & - M \end{matrix}), f = (\begin{matrix} f \\ 0_{n \times 1} \end{matrix})

Note that $det M \neq 0 \Rightarrow det A \neq 0$ . The function $g$ is affine and the theorem now follows from general theorems on existence and uniqueness for first-order ordinary differential equations (see, for instance, Arnold [5]). □

Remark 6.1: A formal solution of the initial value problem $\overset{\cdot}{z} = - A^{- 1} B z + A^{- 1} f, z (0) = z_{0}$ is given by (see Arnold [5])

z (t) = e^{- A^{- 1} B t} z_{0} + \int_{0}^{t} e^{- A^{- 1} B (t - s)} A^{- 1} f (s) ds

The character of this solution is determined by the Jordan canonical form of the matrix $A^{- 1} B$ . This approach is mathematically efficient but it lacks some of the physical transparency obtained by directly considering the solution to Equation (43).■

Remark 6.2: A linear change of coordinates $x \to \hat{x} \in R^{n \times 1}$ is obtained by the relation $x = U \hat{x}$ , where $U \in R^{n \times n}$ is a constant, non-singular matrix. Equation (40) is then equivalent to

\hat{M} \ddot{\hat{x}} + \hat{A} \overset{\cdot}{\hat{x}} + \hat{B} \hat{x} = \hat{f}

(48)

where $\hat{M} = U^{T} MU$ , $\hat{A} = U^{T} AU$ , $\hat{B} = U^{T} BU$ and $\hat{f} = U^{T} f$ . Note that ${\hat{M}}^{T} = \hat{M}$ and $det \hat{M} = (det U)^{2} det M$ , $det \hat{B} = (det U)^{2} det B$ . A coordinate change results in the dynamic stiffness $\hat{S} = U^{T} SU$ . The corresponding characteristic polynomial is given by $\hat{p} (s) = (det U)^{2} p (s)$ . Thus $p (s) = 0 \Leftrightarrow \hat{p} (s) = 0$ . A decoupling of the equations in (40) is obtained if one can find a change of coordinates that renders the matrices $\hat{M}$ , $\hat{A}$ and $\hat{B}$ diagonal. In this case the solution of (48) is more or less trivial. It is well known that if $A$ and $B$ are symmetric (non-gyroscopic and non-circulatory system) and if $A M^{- 1} B = B M^{- 1} A$ then there exists a change of coordinates that decouples the equations (Caughey and O’Kelly [9]).■

7. Free vibrations

The function $x \in C^{2} ([0, \infty [, R^{n \times 1})$ represents a free vibration if $x \in \ker (S (d_{t}))$ .

Proposition 7.1 $\dim (\ker (S (d_{t}))) = 2 n$

Proof: See the Appendix A.3.1.

The result in Proposition 7.1 implies that if one has found $2 n$ solutions to the differential equation $S (d_{t}) x = 0_{n \times 1}$ and can show that these are linearly independent then, by combining them linearly, all solutions have been found. When searching for these solutions it is a good strategy to look among the very regular functions, that is, functions belonging to the space of infinitely differentiable functions, that is, assuming $x \in C^{\infty} ([0, \infty [, R)$ .

Introduce the adjoint differential operator $T (d_{t})$ related to $S (d_{t})$ and defined by $T (d_{t}) = adjS (d_{t}) : C^{\infty} ([0, \infty [, R^{n \times 1}) \to$ $C^{\infty} ([0, \infty [, R^{n \times 1})$ . (See the Appendix for the algebraic definition of the adjoint.) We have

S (d_{t}) T (d_{t}) = T (d_{t}) S (d_{t}) = p (d_{t}) 1_{n \times n}

(49)

where

p (d_{t}) = det (S (d_{t})) = det M d_{t}^{2 n} + \dots + det B

(50)

The kernel of the differential operator $p (d_{t}) : C^{\infty} ([0, \infty [, R) \to C^{\infty} ([0, \infty [, R)$ is given by $\ker (p (d_{t})) = {h \in C^{\infty} ([0, \infty [, R) | p (d_{t}) h = 0}$ . The equation $p (d_{t}) h = 0$ is called the characteristic differential equation. A direct consequence of (49) is the following.

Proposition 7.2

x \in \ker (S (d_{t})) \cap C^{\infty} ([0, \infty [, R^{n \times 1}) \Rightarrow x_{i} \in \ker (p (d_{t})), i = 1, \dots, n

In order to answer the question of the dimension of $\ker (p (d_{t}))$ it turns out to be advantageous to complexify the problem and consider

ke r_{c} (p (d_{t})) = {z \in C^{\infty} ([0, \infty [, C) | p (d_{t}) z = 0}

Note that $z \in C^{\infty} ([0, \infty [, C)$ if and only if $z (t) = x (t) + i y (t)$ , where $x, y \in C^{\infty} ([0, \infty [, R)$ , and since the operator $p (d_{t})$ has real coefficients it follows that $p (d_{t}) z = 0 \Leftrightarrow$ $p (d_{t}) x = p (d_{t}) y = 0$ and $p (d_{t}) z = 0 \Leftrightarrow p (d_{t}) z^{*} = 0$ , where $z^{*} (t) = x (t) - i y (t)$ is the complex conjugate of $z$ .

The main theorem, characterizing members of $ke r_{c} (p (d_{t}))$ , is the following.

Theorem 7.1 If the characteristic equation has the distinct roots $s_{i}, i = 1, \dots, l$ with multiplicities $m_{i}, i = 1, \dots, l$ , ( $\sum_{i = 1}^{l} m_{i} = 2 n$ ), then a solution to the characteristic differential equation is given by the exponential-polynomial (expol):

z (t) = \sum_{i = 1}^{l} g_{i} (t) e^{s_{i} t}

(51)

where $g_{i}$ are polynomials with complex coefficients and $\deg (g_{i}) < m_{i}$ .

Proof: A proof of this may be found in a textbook on ordinary differential equations, for instance Birkhoff and Rota [5] (Chapter 3). A proof, complying with the notation used in this paper, is given in Appendix A.3.2.□

The following conclusion is a direct consequence of the previous theorem.

Corollary 7.1 $\dim (ke r_{c} (p (d_{t}))) \geq 2 n$ .

The roots of the characteristic equation may be grouped into two sets $δ_{1}, \dots, δ_{r} \in R$ and $s_{1}, s_{1}^{*}, \dots, s_{k}, s_{k}^{*} \in C$ with multiplicities $μ_{1}, \dots, μ_{r}$ and $m_{1}, \dots, m_{k}$ , where $r + 2 k = l$ and $\sum_{i = 1}^{r} μ_{i} + 2 \sum_{i = 1}^{k} m_{i} = 2 n$ . The expol (51) may then be written

z (t) = \sum_{i = 1}^{r} f_{i} (t) e^{δ_{i} t} + \sum_{i = 1}^{k} ({\hat{g}}_{i} (t) e^{s_{i} t} + {\hat{h}}_{i} (t) e^{s_{i}^{*} t})

(52)

where $f_{i}$ , ${\hat{g}}_{i}$ and ${\hat{h}}_{i}$ are polynomials with complex coefficients and $\deg (f_{i}) < μ_{i}$ , $\deg ({\hat{g}}_{i}), \deg ({\hat{h}}_{i}) < m_{i}$ . Put $s_{i} = σ_{i} + i ω_{i}, σ_{i}, ω_{i} \in R, i^{2} = - 1$ , where $σ_{i}$ is called the damping factor and $ω_{i}$ is called the (damped) natural frequency. The expol (51) may then be written

z (t) = \sum_{i = 1}^{r} f_{i} (t) e^{δ_{j} t} + \sum_{i = 1}^{k} (g_{i} (t) \cos ω_{i} t + h_{i} (t) \sin ω_{i} t) e^{σ_{i} t}

where $\deg (g_{i}), \deg (h_{i}) < m_{i}$ . Obviously, if we are only interested in real valued expols we may choose polynomials $f_{j}, g_{j}, h_{j}$ with real coefficients. It is clear that the functions

t^{i} e^{δ_{j} t}, i = 0, 1, \dots, μ_{j} - 1, j = 1, \dots, r

(53)

t^{i} e^{σ_{j} t} \cos ω_{j} t, t^{i} e^{σ_{j} t} \sin ω_{j} t, i = 0, 1, \dots, m_{j} - 1, j = 1, \dots, k

are linearly independent over $R$ . We have thus demonstrated that $\dim (\ker (p (d_{t}))) \geq 2 n$ and that $\ker (p (d_{t}))$ contains exponential and harmonic functions combined with polynomials. Are there other functions in $\ker (p (d_{t}))$ ? The answer to this question is, no! Before we prove this we return to our original problem, namely, the structure of $\ker (S (d_{t}))$ . We know that $\dim (\ker (S (d_{t}))) = 2 n$ and that if $x \in \ker (S (d_{t}))$ then $x_{i} \in \ker (p (d_{t})), i = 1, \dots, n$ . As candidates for elements in $\ker (S (d_{t}))$ we now take

x (t) = \sum_{j = 1}^{r} F_{j} (t) e^{δ_{j} t} + \sum_{j = 1}^{k} (G_{j} (t) \cos ω_{j} t + H_{j} (t) \sin ω_{j} t) e^{σ_{j} t}

(54)

where

F_{j} (t) = \sum_{i = 1}^{m_{j} - 1} a_{j}^{i} t^{i}, G_{j} (t) = \sum_{i = 1}^{μ_{j} - 1} b_{j}^{i} t^{i}, H_{j} (t) = \sum_{i = 1}^{μ_{j} - 1} c_{j}^{i} t^{i}

and $a_{j}^{i}, b_{j}^{i}, c_{j}^{i} \in R^{n \times 1}$ are constant column matrices. These functions span a linear space of dimension equal to $2 n^{2}$ . If then $x$ , according to (54), belongs to $\ker (S (d_{t}))$ , then there must be a certain relationship between the vectors $a_{j}^{i}$ , $b_{j}^{i}$ and $c_{j}^{i} \in R^{n \times 1}$ . This relationship is established by the condition $S (d_{t}) x = 0_{n \times 1}$ .

8. Simple roots to the characteristic equation

Assume that all roots of the characteristic equation are simple, that is, $μ_{j} = 1, j = 1, \dots, r$ and $m_{j} = 1, j = 1, \dots, k$ . Then $F_{j} (t) = a_{j}$ , $G_{j} (t) = b_{j}$ and $H_{j} (t) = c_{j}$ , where $a_{j}, b_{j}, c_{j} \in R^{n \times 1}$ are constant vectors. Consequently

S (d_{t}) x = S (d_{t}) (\sum_{j = 1}^{r} a_{j} e^{δ_{j} t} + \sum_{j = 1}^{k} (b_{j} \cos ω_{j} t + c_{j} \sin ω_{j} t) e^{σ_{j} t}) =

\sum_{j = 1}^{r} S (d_{t}) (e^{δ_{j} t}) a_{j} + \sum_{j = 1}^{k} (S (d_{t}) (e^{σ_{j} t} \cos ω_{j} t) b_{j} + S (d_{t}) (e^{σ_{j} t} \sin ω_{j} t) c_{j})

where $S (d_{t}) (e^{δ_{j} t}) = S (δ_{j}) e^{δ_{j} t}$ and $S (d_{t}) (e^{σ_{j} t} \cos ω_{j} t) = (M ((σ_{j}^{2} - ω_{j}^{2}) \sin ω_{j} t + 2 σ_{j} ω_{j} \cos ω_{j} t) + A (σ_{j} \sin ω_{j} t + ω_{j} \cos ω_{j} t) + B \sin ω_{j} t) e^{σ_{j} t}$ . Then

S (d_{t}) x = \sum_{j = 1}^{r} S (δ_{j}) a_{j} e^{δ_{j} t} + \sum_{j = 1}^{k} (S (d_{t}) (e^{σ_{j} t} \cos ω_{j} t) b_{j} + S (d_{t}) (e^{σ_{j} t} \sin ω_{j} t) c_{j}) = \sum_{j = 1}^{r} S (δ_{j}) a_{j} e^{δ_{j} t} +

\sum_{j = 1}^{k} ((M ((σ_{j}^{2} - ω_{j}^{2}) \cos ω_{j} t - 2 σ_{j} ω_{j} \sin ω_{j} t) + A (σ_{j} \cos ω_{j} t - ω_{j} \sin ω_{j} t) + B \cos ω_{j} t) e^{σ_{j} t} b_{j} +

(M ((σ_{j}^{2} - ω_{j}^{2}) \sin ω_{j} t + 2 σ_{j} ω_{j} \cos ω_{j} t) + A (σ_{j} \sin ω_{j} t + ω_{j} \cos ω_{j} t) + B \sin ω_{j} t) e^{σ_{j} t} c_{j} =

\sum_{j = 1}^{r} S (δ_{j}) a_{j} e^{δ_{j} t} + \sum_{j = 1}^{k} ((σ_{j}^{2} - ω_{j}^{2}) M b_{j} + 2 σ_{j} ω_{j} M c_{j} + σ_{j} A b_{j} + ω_{j} A c_{j} + B b_{j}) e^{σ_{j} t} \cos ω_{j} t +

\sum_{j = 1}^{k} ((σ_{j}^{2} - ω_{j}^{2}) M c_{j} - 2 σ_{j} ω_{j} M b_{j} + σ_{j} A c_{j} - ω_{j} A b_{j} + B c_{j}) e^{σ_{j} t} \sin ω_{j} t = 0_{n \times 1}, \forall t \in [0, \infty [

which is true if and only if the vectors $a_{j}$ , $b_{j}$ , $c_{j}$ satisfy the equations

{\begin{matrix} S (δ_{j}) a_{j} = 0_{n \times 1}, j = 1, . ., r \\ ((σ_{j}^{2} - ω_{j}^{2}) M + σ_{j} A + B) b_{j} + ω_{j} (2 σ_{j} M + A) c_{j} = 0_{n \times 1} \\ - ω_{j} (2 σ_{j} M + A) b_{j} + ((σ_{j}^{2} - ω_{j}^{2}) M + σ_{j} A + B) c_{j} = 0_{n \times 1} \\ j = 1, \dots, k \end{matrix}

(55)

If instead one uses complex expols

z (t) = \sum_{j = 1}^{r} a_{j} e^{δ_{j} t} + \sum_{j = 1}^{k} (b_{j} e^{s_{j} t} + c_{j} e^{s_{j}^{*} t})

(56)

where now $a_{j}$ , $b_{j}$ , $c_{j} \in C^{n \times 1}$ , then

S (d_{t}) z = \sum_{j = 1}^{r} S (d_{t}) e^{δ_{j} t} a_{j} + \sum_{j = 1}^{k} (S (d_{t}) e^{s_{j} t} b_{j} + S (d_{t}) e^{s_{j}^{*} t} c_{j}) = \sum_{j = 1}^{r} S (δ_{j}) a_{j} e^{δ_{j} t} +

\sum_{j = 1}^{k} (S (s_{j}) b_{j} e^{s_{j} t} + S (s_{j}^{*}) c_{j} e^{s_{j}^{*} t}) = 0_{n \times 1}, \forall t \in [0, \infty [

which is true if and only if the vectors $a_{j}$ , $b_{j}$ , $c_{j}$ satisfy the equations

{\begin{matrix} S (δ_{j}) a_{j} = 0_{n \times 1}, j = 1, . ., r \\ S (s_{j}) b_{j} = 0_{n \times 1} \\ S (s_{j}^{*}) c_{j} = 0_{n \times 1} \\ j = 1, \dots, k \end{matrix}

These conditions appear much simpler than the (equivalent) ones given in (55), this being a beneficial consequence of complexification. Thus for each $s \in spec (S)$ one should determine $w \in C^{n \times 1}$ so that

S (s) w = 0_{n \times 1}

(57)

which means that we are now looking for the kernel $\ker (S (s)) = {w \in C^{n \times 1} | S (s) w = 0_{n \times 1}}$ . Obviously, $s \notin spec (S) \Rightarrow \ker (S (s)) = {0_{n \times 1}}$ . If $s \in spec (S)$ then there exists a vector $w \in C^{n \times 1}, w \neq 0_{n \times 1}$ satisfying (57), that is, $\ker (S (s)) \neq {0_{n \times 1}}$ . The column matrix $w \in C^{n \times 1}, w \neq 0_{n \times 1}$ is called the eigenvector corresponding to the eigenvalue $s \in spec (S)$ . With the notation

w = u + iv, w^{*} = u - iv, u, v \in R^{n \times 1}, i^{2} = - 1

and $s \in spec (S)$ we have $S (s) w = 0_{n \times 1} \Leftrightarrow S (s^{*}) w^{*} = 0_{n \times 1}$ . Note that if $s = δ_{j} \in R$ we may choose $v = 0_{n \times 1}$ and then $w = u \in R^{n \times 1}$ . Now assume that

a_{j} \in R^{n \times 1}, j = 1, \dots, r, w_{j}, w_{j}^{*} \in C^{n \times 1}, j = 1, \dots, k

(58)

are solutions to (57) corresponding to $s = δ_{j} \in R$ , $s = s_{j} \in C$ and $s = s_{j}^{*} \in C$ , respectively. Then (56) may be written

z (t) = \sum_{j = 1}^{r} α_{j} a_{j} e^{δ_{j} t} + \sum_{j = 1}^{k} (β_{j} w_{j} e^{s_{j} t} + γ_{j} w_{j}^{*} e^{s_{j}^{*} t})

(59)

where $α_{j}, β_{j}, γ_{j} \in C$ . The expression (59) contains $r + 2 k = 2 n$ arbitrary complex constants. Note that the column vectors in (59) are linearly dependent. The expression is real valued if and only if $α_{j} \in R$ , $β_{j}^{*} = γ_{j}$ .

Proposition 8.1 Real valued solutions to $S (d_{t}) = 0_{n \times 1}$ of the expol type are given by

x (t) = \sum_{j = 1}^{r} α_{j} a_{j} e^{δ_{j} t} + \sum_{j = 1}^{k} (A_{j} (u_{j} \cos (ω_{j} t + ϕ_{j}) - v_{j} \sin (ω_{j} t + ϕ_{j})) e^{σ_{j} t}

(60)

where $α_{j}, A_{j}, ϕ_{j} \in R$ and $a_{j}, w_{j} = u_{j} + i v_{j}$ satisfy (57).

Proof: We may write $β_{j} = \frac{A_{j}}{2} e^{i ϕ_{j}}, A_{j}, ϕ_{j} \in R$ and $w_{j} = u_{j} + i v_{j}, u_{j}, v_{j} \in R^{n \times 1}$ and $s_{j} = σ_{j} + i ω_{j}, σ_{j}, ω_{j} \in R$ . Then (59) may then be written

\begin{matrix} z (t) = \sum_{j = 1}^{r} α_{j} a_{j} e^{δ_{j} t} + \sum_{j = 1}^{k} (\frac{A_{j}}{2} e^{i ϕ_{j}} (u_{j} + i v_{j}) e^{s_{j} t} + \frac{A_{j}}{2} e^{- i ϕ_{j}} (u_{j} - i v_{j}) e^{s_{j}^{*} t}) = \\ \sum_{j = 1}^{r} α_{j} a_{j} e^{δ_{j} t} + \sum_{j = 1}^{k} A_{j} (u_{j} \cos (ω_{j} t + ϕ_{j}) - v_{j} \sin (ω_{j} t + ϕ_{j})) e^{σ_{j} t} \end{matrix}

Theorem 8.1 If the roots of the characteristic equation are simple then

\ker (S (d_{t})) = {x = x (t) according to (60) | A_{j}, ϕ_{j} \in R, and a_{j}, w_{j} = u_{j} + i v_{j} satisfy (57)}

Proof: Note that $x = x (t)$ , according to (60), contains $2 n$ arbitrary constants: $α_{j}$ , $A_{j}$ and $ϕ_{j}$ . Thus these functions constitute a linear space with dimension $= 2 n$ . However, according to Proposition 7.1, $\dim (\ker (S (d_{t}))) = 2 n$ and this proves the theorem. □

Remark 8.1: We have the following alternative representation of $\ker (S (d_{t}))$ :

\ker (S (d_{t})) = {x = x (t) according to (59) | α_{j} \in R, β_{j}^{*} = γ_{j} \in C and a_{j}, w_{j}, w_{j}^{*} \in C^{n \times 1}

satisfying S (δ_{j}) a_{j} = 0_{n \times 1}, S (s_{j}) w = 0_{n \times 1}, S (s_{j}^{*}) w^{*} = 0_{n \times 1}}

(61)

■

The oscillating terms in (60) are on the form

\tilde{x} (t) = A (u \cos (ω t + ϕ) - v \sin (ω t + ϕ)) e^{σ t}

(62)

This is called the complex natural mode with (damped) natural frequency $ω$ and damping factor $σ$ . With $σ < 0$ the expression (61) represents, for fixed $A, ϕ$ , an oscillating motion with exponentially decreasing amplitude. Component-wise we have

{\tilde{x}}_{i} (t) = A (u_{i} \cos (ω t + ϕ) - v_{i} \sin (ω t + ϕ)) e^{σ t} = A e^{σ t} \sqrt{u_{i}^{2} + v_{i}^{2}} \sin (ω t + ϕ + θ_{i}), i = 1, \dots, n

where $θ_{i} = - \arctan (\frac{v_{i}}{u_{i}})$ . The mode component is equal to zero at times

t_{i, k} = \frac{1}{ω} (- ϕ - θ_{i} + k π), k = \dots - 2, - 1, 0, 1, 2, \dots

These times are in general depending on the specific mode component ${\tilde{x}}_{i}$ . However, if $θ_{i} = θ$ , independent of $i$ then the mode components will be equal to zero simultaneously. In this case

\frac{v_{i}}{u_{i}} = - \tan θ \Leftrightarrow v = - u \tan θ

and then

{\tilde{x}}_{i} = A e^{σ t} \sqrt{1 + \tan^{2} θ} u_{i} \sin (ω t + ϕ + θ)

This is called a real (natural) mode. Note that $w = u + iv = u + i (- u \tan θ) = (1 - i \tan θ) u$ and with $θ = 0$ one obtains $v = 0$ and $w = u$ .

Remark 8.2: Real natural modes are present in many multibody systems. If $G = 0_{n \times n}$ , $H = 0_{n \times n}$ (non-gyroscopic and non-circulatory) and if $C M^{- 1} K = K M^{- 1} C$ then all modes are real. This situation is, for instance, at hand in the case of proportional damping where it is assumed that $C = M \sum_{i = 1}^{k} α_{i} {(M^{- 1} K)}^{i}$ , where $k = 1, 2, \dots$ and $α_{i} \in R, i = 1, \dots, k$ (see Clough and Penzien [13]).■

The following corollary is an immediate consequence of Theorem 8.1.

Corollary 8.1 If the characteristic equation has simple roots then

s \in spec (S) \Rightarrow \dim (\ker (S (s))) = 1

Proof: We have $s \in spec (S) \Rightarrow \dim (\ker (S (s))) \geq 1$ . Now assume that for some $s \in spec (S)$ one has $\dim (\ker (S (s))) \geq 2$ , then it follows by inspection of (60) that $\dim (\ker (S (d_{t}))) \geq 2 n + 1$ , which contradicts the fact that, according to Proposition 7.1, $\dim (\ker (S (d_{t}))) = 2 n$ . Note that given constant, linearly independent vectors $w_{1}, \dots, w_{k} \in C^{n \times 1}$ and $s \in C$ , then the functions $t \to w_{i} e^{st}, i = 1, \dots, k$ are linear independent functions over $C$ in the space $C^{\infty} ([0, \infty [, C^{n \times 1})$ .

This proves the proposition. □

We now summarize the discussion above by giving the following well-known scheme for solving the equation $S (d_{t}) x = 0_{n \times 1}$ .

Set up the dynamic stiffness $S (s)$ .

Calculate the characteristic polynomial $p (s) = det S (s)$ .

Determine the roots $s_{j}, j = 1, \dots, 2 n$ of the characteristic equation $p (s) = 0$ . It is assumed that all roots are simple.

Determine $w_{j} = u_{j} + i v_{j} \neq 0_{n \times 1}$ so that $S (s_{j}) w_{j} = 0_{n \times 1}, j = 1, \dots, 2 n$ .

Insert $δ_{j}, σ_{j}, ω_{j}$ and $u_{j}, v_{j}$ into (60).

In the present case of simple roots we may read off the solution in step (iv) directly from the adjoint matrix $T = T (s) = adj S (s)$ to the dynamic stiffness: We have

S (s) T (s) = T (s) S (s) = p (s) 1_{n \times n}, s \in C

(63)

Proposition 8.2 If $s_{j}$ is a simple root to the characteristic equation then $T (s_{j}) \neq 0_{n \times n}$ and $w_{j}$ (in step (iv) above) may be taken as proportional to any non-zero column vector in the matrix $T (s_{j})$ .

Proof: From relation (63) one obtains the implication $p (s_{j}) = 0 \Rightarrow S (s_{j}) T (s_{j}) = 0_{n \times n}$ . Let $T (s_{j}) = (\begin{matrix} T_{1} (s_{j}) & \dots & T_{n} (s_{j}) \end{matrix})$ , where $T_{i} (s_{j}) \in C^{n \times 1}$ is the ith column in the matrix $T (s_{j})$ . Then $S (s_{j}) T_{i} (s_{j}) = 0_{n \times 1}, j = 1, \dots, n$ . Thus, every $T_{i} (s_{j}) \neq 0_{n \times 1}$ is an eigenvector. The existence of at least one column vector $T_{i} (s_{j}) \neq 0_{n \times 1}$ follows from the assumption that $s_{j}$ is a simple root. By differentiating (63) with respect to $s$ , it follows that

T' (s) S (s) + T (s) S' (s) = p' (s) 1_{n \times n}, s \in C

(64)

and then

T' (s_{j}) S (s_{j}) + T (s_{j}) S' (s_{j}) = p' (s_{j}) 1_{n \times n},

Now if $S (s_{j}) w_{j} = 0_{n \times 1}, w_{j} \neq 0_{n \times 1}$ then

T (s_{j}) S' (s_{j}) w_{j} = p' (s_{j}) w_{j} \neq 0_{n \times 1}

since, according to Lemma 6.1, $p' (s_{j}) \neq 0$ . From this we conclude that $T (s_{j}) \neq 0_{n \times 1}$ . This proves the proposition. □

Remark 8.3: If more than one column in $T (s_{j})$ is non-zero, let us say $T_{i} (s_{j}) \neq 0_{n \times 1}$ and $T_{k} (s_{j}) \neq 0_{n \times 1}$ , then these columns must be linearly dependent.■

If the characteristic equation has simple roots then the function $x \in \ker (S (d_{t}))$ , satisfying the initial conditions $x (0) = x_{0}, \overset{\cdot}{x} (0) = {\overset{\cdot}{x}}_{0}$ , is determined by the conditions

\sum_{j = 1}^{r} α_{j} a_{j} + \sum_{j = 1}^{k} (A_{j} (u_{j} \cos (ϕ_{j}) - v_{j} \sin (ϕ_{j})) = x_{0}

\sum_{j = 1}^{r} d_{j} α_{j} a_{j} + \sum_{j = 1}^{k} A_{j} (\cos φ_{j} (σ_{j} u_{j} - ω_{j} v_{j}) - \sin φ_{j} (ω_{j} u_{j} + σ_{j} v_{j})) = {\overset{\cdot}{x}}_{0}

or in matrix format

A u = (\begin{matrix} x_{0} \\ {\overset{\cdot}{x}}_{0} \end{matrix})

(65)

where

A = (\begin{matrix} a_{1} & \dots & a_{r} & u_{1} & - v_{1} & \dots & u_{k} & - v_{k} \\ d_{1} a_{1} & \dots & d_{r} a_{r} & σ_{1} u_{1} - ω_{1} v_{1} & - (ω_{1} u_{1} + σ_{1} v_{1}) & \dots & σ_{k} u_{k} - ω_{k} v_{k} & - (ω_{k} u_{k} + σ_{k} v_{k}) \end{matrix})

u = (\begin{matrix} α_{1} & \dots & α_{r} & A_{1} \cos φ_{1} & A_{1} \sin φ_{1} & \dots & A_{k} \cos φ_{k} & A_{k} \sin φ_{k} \end{matrix})^{T}

The $2 n$ real constants $α_{1}, \dots, α_{r}$ , $A_{1}, \dots, A_{k}$ and $φ_{1}, \dots, φ_{k}$ are uniquely determined by (65). This follows from the uniqueness theorem, Theorem 6.1. The matrix $A \in R^{2 n \times 2 n}$ is thus non-singular. Note, however, that the vectors $a_{1}, \dots, a_{r}, u_{1}, v_{1}, \dots, u_{k}, v_{k}$ must be linearly dependent. For special properties of the coefficient matrices $A$ and $B$ in (40) the modes will fulfil certain orthogonality properties that will simplify the solution of (65). This is, for instance, the case for the system defined in Remark 8.2 above. In this case we have real modes $u_{i}, i = 1, \dots, n$ and orthogonality properties according to Corollary 8.5 below.

Proposition 8.3 If $s_{i}, s_{j} \in spec (S)$ and $w_{i} \in \ker (S (s_{i}))$ , $w_{j} \in \ker (S (s_{j}))$ then

w_{j}^{T} M w_{i} (s_{i}^{2} - s_{j}^{2}) + w_{j}^{T} (A s_{i} - A^{T} s_{j}) w_{i} + w_{j}^{T} (B - B^{T}) w_{i} = 0

(66)

w_{j}^{T} M w_{i} s_{i} s_{j} (s_{i} - s_{j}) + w_{j}^{T} (A - A^{T}) w_{i} s_{i} s_{j} + w_{j}^{T} (B s_{j} - B^{T} s_{i}) w_{i} = 0

Proof: We have, since $M \in sym (R^{n \times n})$ :

w_{j}^{T} S (s_{i}) w_{i} = w_{j}^{T} M w_{i} s_{i}^{2} + w_{j}^{T} A w_{i} s_{i} + w_{j}^{T} B w_{i} = 0

w_{i}^{T} S (s_{j}) w_{j} = w_{i}^{T} M w_{j} s_{j}^{2} + w_{i}^{T} A w_{j} s_{j} + w_{i}^{T} B w_{j} = w_{j}^{T} M w_{i} s_{j}^{2} + w_{j}^{T} A^{T} w_{i} s_{j} + w_{j}^{T} B^{T} w_{i} = 0

Furthermore

w_{j}^{T} M w_{i} s_{i}^{2} s_{j} + w_{j}^{T} A w_{i} s_{i} s_{j} + w_{j}^{T} B w_{i} s_{j} = 0, w_{j}^{T} M w_{i} s_{j}^{2} s_{i} + w_{j}^{T} A^{T} w_{i} s_{j} s_{i} + w_{j}^{T} B^{T} w_{i} s_{i} = 0

Subtracting these expressions pairwise gives (66). □

The following corollaries are direct consequences of the previous proposition.

Corollary 8.2 If $s_{i} = s_{j}$ and if $w_{i}, w_{j} \in \ker (S (s_{i}))$ then

w_{j}^{T} (A - A^{T}) w_{i} s_{i} + w_{j}^{T} (B - B^{T}) w_{i} = 0

Corollary 8.3 If $A, B \in sym (R^{n \times n})$ then, if $s_{i} \neq s_{j}$ and $w_{i} \in \ker (S (s_{i}))$ , $w_{j} \in \ker (S (s_{j}))$ :

w_{j}^{T} (S (s_{i}) - S (s_{j})) w_{i} = 0

w_{j}^{T} M w_{i} (s_{i} + s_{j}) + w_{j}^{T} A w_{i} = 0, w_{j}^{T} M w_{i} s_{i} s_{j} - w_{j}^{T} B w_{i} = 0

Corollary 8.4 If $A \in skew (R^{n \times n})$ $B \in sym (R^{n \times n})$ then, if $s_{i} \neq s_{j}$ and $w_{i} \in \ker (S (s_{i}))$ , $w_{j} \in \ker (S (s_{j}))$ :

w_{j}^{T} M w_{i} (s_{i} - s_{j}) + w_{j}^{T} A w_{i} = 0

Corollary 8.5 If $A, B \in sym (R^{n \times n})$ and if the modes are real $w_{i} = u_{i} \in R^{n \times 1}$ then, if $s_{i} \neq s_{j}$ , $ω_{i} + ω_{j} \neq 0$ and $w_{i} \in \ker (S (s_{i}))$ , $w_{j} \in \ker (S (s_{j}))$ :

{u_{j}}^{T} M u_{i} = {u_{j}}^{T} A u_{i} = {u_{j}}^{T} B u_{i} = 0

Proof: With $s_{i} = σ_{i} + i ω_{i}$ , $s_{j} = σ_{j} + i ω_{j}$ one obtains from (66)₁

u_{j}^{T} M u_{i} (σ_{i} + σ_{j}) + u_{j}^{T} A u_{i} + i u_{j}^{T} M u_{i} (ω_{i} + ω_{j}) = 0_{n \times 1} \Rightarrow u_{j}^{T} M u_{i} = 0, u_{j}^{T} A u_{i} = 0

and then from (66)₂ $u_{j}^{T} M u_{i} s_{i} s_{j} + u_{j}^{T} B u_{i} = 0_{n \times 1} \Rightarrow u_{j}^{T} B u_{i} = 0$ . □

Remark 8.4: The previous corollary contains the classical orthogonality conditions obtained by assuming Rayleigh damping, that is, $G = 0_{n \times n}$ , $H = 0_{n \times n}$ and $C = α M + β K$ , where $α, β \in R$ (see Lord Rayleigh [1]).

9. Multiple roots

The appearance of multiple roots in a multibody system is often a consequence of some symmetry of the system. Systems with multiple roots are, however, sensitive to disturbances. They may, from the point of view of applications, often be considered as a ‘degenerated’ limit case that may be done away with by arbitrarily small changes of the coefficient matrices $M, A$ and $B$ . Nevertheless it seems important, also from the engineering point of view, to be aware of some of the consequences of the presence of multiple roots. Investigations of multibody systems with general coefficient matrices $A$ and $B$ and with multiple roots are, however, sparse in the engineering literature. One of the few references we have found is Wu and Greif [14], where a vibrating system with null as well as non-null eigenvalues is analysed.

The algebraic multiplicity of a root to the characteristic equation will, in general, give rise to mathematical as well as numerical complications. Let $s$ be a root of the characteristic equation with algebraic multiplicity $m$ . The geometric multiplicity $ν$ of $s$ is defined by $ν = \dim \ker (S (s))$ . It is always the case that $ν \leq m$ . An eigenvalue $s$ is called defect if $ν < m$ . It is the presence of defect eigenvalues that makes the situation for multiple roots more complicated than for simple roots. Note that, in the case of simple roots, it follows from Corollary 8.1 that $ν = m = 1$ .

If $s$ is a double root to the characteristic equation then $x (t) = (w_{0} + w_{1} t) e^{st} \in \ker (S (d_{t}))$ if and only if

S (d_{t}) x = 0_{n \times 1}, t \in [0, \infty [

(67)

where

\begin{matrix} S (d_{t}) (w_{0} e^{st}) = e^{st} S (s) w_{0} and S (d_{t}) (w_{1} t e^{st}) = ({Md}_{t}^{2} + A d_{t} + B) (t e^{st}) w_{1} = \\ e^{st} (2 Ms + A) w_{1} + t e^{st} (M s^{2} + As + B) w_{1} = e^{st} S' (s) w_{1} + t e^{st} S (s) w_{1} \end{matrix}

Thus (67) is equivalent to

e^{st} (S (s) w_{0} + S' (s) w_{1}) + t e^{st} S (s) w_{1} = 0_{n \times 1}, \forall t \in [0, \infty [

which is equivalent to the conditions

{\begin{matrix} S' (s) w_{1} + S (s) w_{0} = 0_{n \times 1} \\ S (s) w_{1} = 0_{n \times 1} \end{matrix}

(68)

Note that since $s \in spec (S)$ there always exists a vector $w_{0} \neq 0_{n \times 1}$ such that $S (s) w_{0} = 0_{n \times 1}$ and then $w_{0}$ , $w_{1} = 0_{n \times 1}$ is a solution to (68), called the trivial solution. A non-trivial solution is obtained if $w_{1} \neq 0_{n \times 1}$ .

If $s$ is a triple root to the characteristic equation then $x (t) = (w_{0} + w_{1} t + w_{2} t^{2}) e^{st} \in$ $\ker (S (d_{t}))$ if and only if

{\begin{matrix} S ″ (s) w_{2} + S' (s) w_{1} + S (s) w_{0} = 0_{n \times 1} \\ 2 S' (s) w_{2} + S (s) w_{0} = 0_{n \times 1} \\ S (s) w_{2} = 0_{n \times 1} \end{matrix}

(69)

This may be generalized according to the following.

Theorem 9.1 If $s$ is a multiple root to the characteristic equation with multiplicity $m$ , then the function

x (t) = \sum_{k = 0}^{m - 1} w_{k} t^{k} e^{st}, w_{k} \in C^{n \times 1}

is a solution to (67) if and only if

{\begin{matrix} S (s) w_{k} + (k + 1) S' (s) w_{k + 1} + \frac{1}{2} (k + 2) (k + 1) S ″ (s) w_{k + 2} = 0_{n \times 1}, 0 \leq k \leq m - 3 \\ S (s) w_{m - 2} + (m - 1) S' (s) w_{m - 1} = 0_{n \times 1} \\ S (s) w_{m - 1} = 0_{n \times 1} \end{matrix}

(70)

For the proof of this theorem we need the following.

Lemma 9.1

{\begin{matrix} (\frac{d}{dt} - s) (t^{k} e^{st}) = k t^{k - 1} e^{st}, if k \geq 1 \\ (\frac{d}{dt} - s) (t^{k} e^{st}) = 0, if k = 0 \end{matrix}, {\begin{matrix} (\frac{d}{dt} - s)^{2} (t^{k} e^{st}) = k (k - 1) t^{k - 2} e^{st}, if k \geq 2 \\ (\frac{d}{dt} - s)^{2} (t^{k} e^{st}) = 0, if k \leq 1 \end{matrix}

Proof: Follows by a straightforward computation. □

Proof: (Theorem 9.1) We have, for $λ, s \in C$ :

S (λ) = S (s) + S' (s) (λ - s) + \frac{1}{2} S ″ (s) (λ - s)^{2}

By replacing $λ$ with $d_{t}$ one obtains

S (d_{t}) = S (s) + S' (s) (d_{t} - s) + \frac{1}{2} S ″ (s) (d_{t} - s)^{2}

and then

S (d_{t}) (t^{k} e^{st}) = S (s) (t^{k} e^{st}) + S' (s) (d_{t} - s) (t^{k} e^{st}) + \frac{1}{2} S ″ (s) (d_{t} - s)^{2} (t^{k} e^{st})

This expression, together with Lemma 9.1, results in

S (d_{t}) (t^{k} e^{st}) = {\begin{matrix} e^{st} S (s), if k = 0 \\ e^{st} (tS (s) + S' (s)), if k = 1 \\ e^{st} (t^{k} S (s) + t^{k - 1} kS' (s) + t^{k - 2} \frac{1}{2} k (k - 1) S ″ (s)), if 2 \leq k \leq m - 1 \end{matrix}

Thus

S (d_{t}) x (t) = \sum_{k = 0}^{m - 1} w_{k} S (d_{t}) (t^{k} e^{st}) = e^{st} S (s) w_{0} + e^{st} (tS (s) + S' (s)) w_{1} +

\sum_{k = 2}^{m - 1} e^{st} (t^{k} S (s) + t^{k - 1} kS' (s) + t^{k - 2} \frac{1}{2} k (k - 1) S ″ (s)) w_{k} = e^{st} S (s) w_{0} + e^{st} (tS (s) + S' (s)) w_{1} +

\sum_{k = 2}^{m - 1} e^{st} t^{k} S (s) w_{k} + \sum_{k = 1}^{m - 2} e^{st} t^{k} (k + 1) S' (s) w_{k + 1} + \sum_{k = 0}^{m - 3} e^{st} t^{k} \frac{1}{2} (k + 2) (k + 1) S ″ (s) w_{k + 2} =

\sum_{k = 0}^{m - 3} e^{st} t^{k} (S (s) w_{k} + (k + 1) S' (s) w_{k + 1} + \frac{1}{2} (k + 2) (k + 1) S ″ (s) w_{k + 3}) +

e^{st} t^{m - 2} (S (s) w_{m - 2} + (m - 1) S' (s) w_{m - 1}) + e^{st} t^{m - 1} S (s) w_{m - 1} = 0_{n \times 1}, \forall t \in [0, \infty [

and this is equivalent to (70). □

If $s \in spec (S)$ there always exists a vector $w_{0} \neq 0_{n \times 1}$ such that $S (s) w_{0} = 0_{n \times 1}$ . A trivial solution to (70) is then given by $w_{0}$ , $w_{1} = \dots = w_{m - 1} = 0_{n \times 1}$ . Equations (70), with the requirement $w_{m - 1} \neq 0_{n \times 1}$ , define a generalized eigenvalue problem where $w_{0}, \dots, w_{m - 1} \in C^{n \times 1}$ are called generalized eigenvectors. A sequence of vectors $w_{0}, \dots, w_{m - 1} \in C^{n \times 1}$ ( $w_{m - 1} \neq 0_{n \times 1}$ ) satisfying (70) is, in the mathematical literature, referred to as a Jordan chain for $S (s)$ of length $m$ (see Gohberg et al. [15]). Note that some of the vectors $w_{0}, \dots, w_{m - 2}$ in the chain might be equal to zero.

If $s_{j}$ is a double root to the characteristic equation then, according to Lemma 6.1, $p (s_{j}) = 0$ , $p' (s_{j}) = 0$ and $p ″ (s_{j}) \neq 0, i = 1, \dots, 2 n$ . It then follows that

S (s_{j}) T (s_{j}) = 0_{n \times n}, S' (s_{j}) T (s_{j}) + S (s_{j}) T' (s_{j}) = 0_{n \times n}

Now let $T (s_{j}) = (\begin{matrix} T_{1} (s_{j}) & \dots & T_{n} (s_{j}) \end{matrix})$ and $T' (s_{j}) = (\begin{matrix} {T_{1}}^{'} (s_{j}) & \dots & {T_{n}}^{'} (s_{j}) \end{matrix})$ , where $T_{i} (s_{j})$ , ${T_{i}}^{'} (s_{j}) \in C^{n \times 1}$ . Assume that $T_{i} (s_{j}) \neq 0_{n \times 1}$ for some $i$ . Then $w_{1} = T_{i} (s_{j})$ , $w_{0} = {T_{i}}^{'} (s_{j})$ constitute a Jordan chain. If $T (s_{j}) = 0_{n \times n}$ then ${T_{i}}^{'} (s_{j}) \neq 0_{n \times 1}$ for some $i$ , since by differentiating (63) twice one obtains

T ″ (s) S (s) + 2 T' (s) S' (s) + T (s) S ″ (s) = p ″ (s) 1_{n \times n}, s \in C

and then if $S (s_{j}) w = 0_{n \times 1}, w \neq 0_{n \times 1}$

2 T' (s_{j}) S' (s_{j}) w = p ″ (s_{j}) w \neq 0_{n \times 1} \Rightarrow T' (s_{j}) \neq 0_{n \times n}

A trivial solution to (70) is then given by $w_{0} = {T_{i}}^{'} (s_{j})$ , $w_{1} = 0_{n \times 1}$ . This procedure may be generalized to a root $s_{j}$ with multiplicity $m_{j}$ .

Proposition 9.1 $s_{j} \in spec (S)$ , with multiplicity $m_{j}$ , if and only if

S (s_{j}) T (s_{j}) = 0_{n \times n}

S' (s_{j}) T (s_{j}) + S (s_{j}) T' (s_{j}) = 0_{n \times n}

S (s_{j}) T^{(k)} (s_{j}) + S' (s_{j}) k T^{(k - 1)} (s_{j}) + \frac{k (k - 1)}{2} S ″ (s_{j}) T^{(k - 2)} (s_{j}) = 0_{n \times n}, 2 \leq k \leq m_{j} - 1

(71)

S (s_{j}) T^{(m_{j})} (s_{j}) + m_{j} S' (s_{j}) T^{(m_{j} - 1)} (s_{j}) + \frac{m_{j} (m_{j} - 1)}{2} S ″ (s_{j}) T^{(m_{j} - 2)} (s_{j}) \neq 0_{n \times n}

Proof: By differentiating the identity (63) $k$ times one obtains

(S (s) T (s))^{(k)} = p^{(k)} (s) 1_{n \times n}, k = 1, 2, \dots

For $k = 1$ this results in

S' (s) T (s) + S (s) T' (s) = p' (s) 1_{n \times n}

For $k = 2, 3, \dots$ , using the Leibnitz formula, one obtains

\begin{matrix} (S (s) T (s))^{(k)} = \sum_{i = 0}^{k} (\begin{matrix} k \\ i \end{matrix}) S^{(i)} (s) T^{(k - i)} (s) = (\begin{matrix} k \\ 0 \end{matrix}) S (s) T^{(k)} (s) + (\begin{matrix} k \\ 1 \end{matrix}) S' (s) T^{(k - 1)} (s) + \\ (\begin{matrix} k \\ 2 \end{matrix}) S ″ (s) T^{(k - 2)} (s) = S (s) T^{(k)} (s) + kS' (s) T^{(k - 1)} (s) + \frac{k (k - 1)}{2} S ″ (s) T^{(k - 2)} (s) \end{matrix}

since $S^{(i)} (s) = 0_{n \times n}$ for $i \geq 3$ . The proposition now follows from Lemma 6.1. □

The structural similarity between Equations (70) and (71) is apparent and suggests a possibility of using the adjoint in the solution of the generalized eigenvalue problem.

Remark 9.1: Note that by changing the order of $S$ and $T$ one obtains

(T (s) S (s))^{(k)} = T^{(k)} (s) S (s) + k T^{(k - 1)} (s) S' (s) + \frac{k (k - 1)}{2} T^{(k - 2)} (s) S ″ (s)

and thus conditions (71)₄ is equivalent to

T^{(m_{j})} (s_{j}) S (s_{j}) + m_{j} S' (s_{j}) T^{(m_{j} - 1)} (s_{j}) + \frac{m_{j} (m_{j} - 1)}{2} T^{(m_{j} - 2)} (s_{j}) S ″ (s_{j}) \neq 0_{n \times n}

■

Remark 9.2: The reduced adjoint $Q (s)$ to the dynamic stiffness $S (s)$ is a matrix polynomial satisfying $S (s) Q (s) = i_{n} (s) 1_{n \times n}$ , where $i_{n} (s)$ is the nth invariant polynomial of $S$ (see Appendix A.1.1). It follows that $p (s_{i}) = 0 \Leftrightarrow i_{n} (s_{i}) = 0$ and then $S (s_{i}) Q (s_{i}) = 0_{n \times n}$ , where $Q (s_{i}) \neq 0_{n \times n}$ . Every non-zero column vector in $Q (s_{i})$ is then an eigenvector to $S$ corresponding to the eigenvalue $s_{i}$ . Now if $p (s) \neq 0$ then, by using the Smith normal form, Theorem A.1.1, one obtains

\begin{matrix} Q (s) = i_{n} (s) S (s)^{- 1} = i_{n} (s) F (s)^{- 1} I (s)^{- 1} E (s)^{- 1} = \\ F (s)^{- 1} diag (\begin{matrix} q_{1} (s) & \dots & q_{n - 1} (s) & 1 \end{matrix}) E (s)^{- 1} \end{matrix}

(72)

where

q_{k} (s) = \frac{i_{n} (s)}{i_{k} (s)} = Π_{j = 1}^{l} {(s - s_{j})}^{μ_{jn} - μ_{jk}}

The invariant polynomials $i_{k} (s), k = 1, \dots, n$ are defined in Appendix A.1. Since the quotient $q_{k}$ is defined for all $s \in C$ , the expression for $Q$ in (72) is valid for all $s \in C$ . Then

Q (s_{i}) = i_{n} (s_{i}) F (s_{i})^{- 1} diag (\begin{matrix} q_{1} (s_{i}) & \dots & q_{n - 1} (s_{i}) & 1 \end{matrix}) E (s_{i})^{- 1}

where

q_{k} (s_{i}) = Π_{j = 1}^{l} {(s_{i} - s_{j})}^{μ_{jn} - μ_{jk}} \neq 0 \Leftrightarrow μ_{jk} = μ_{jn}, k = 1, \dots, n - 1

The number of $q_{k} (s_{i}) \neq 0, k = 1, \dots, n$ will then be equal to the number of elementary divisors of maximal degree. The rank of $Q (s_{i})$ is equal to one plus the number of $q_{k} (s_{i}) \neq 0$ and thus equal to the number of elementary divisors of maximal degree (see Lancaster and Webber [16] (Theorem 1)). Note that if we have simple roots, as in Section 8, then

i_{k} (s) = 1, k = 1, \dots, n - 1, i_{n} (s) = Π_{j = 1}^{2 n} (s - s_{j}) = \frac{1}{det M} p (s)

and

q_{k} (s) = \frac{i_{n} (s)}{i_{k} (s)} = \frac{1}{det M} p (s) \Rightarrow q_{k} (s_{i}) = 0, i = 1, \dots, n - 1 \Rightarrow rankQ (s_{i}) = 1

■

Starting from the identity $S (s) Q (s) = i_{n} (s) 1_{n \times n}$ , in the previous remark, we may now formulate the following.

Proposition 9.2 $s_{j}$ is a zero of the polynomial $i_{n} = i_{n} (s)$ , with multiplicity $α_{j}$ , if and only if

S (s_{j}) Q (s_{j}) = 0_{n \times n}

S' (s_{j}) Q (s_{j}) + S (s_{j}) Q' (s_{j}) = 0_{n \times n}

S (s_{j}) Q^{(k)} (s_{j}) + S' (s_{j}) k Q^{(k - 1)} (s_{j}) + \frac{k (k - 1)}{2} S ″ (s_{j}) Q^{(k - 2)} (s_{j}) = 0_{n \times n}, 2 \leq k \leq α_{j} - 1

S (s_{j}) Q^{(α_{j})} (s_{j}) + α_{j} S' (s_{j}) Q^{(α_{j} - 1)} (s_{j}) + \frac{α_{j} (α_{j} - 1)}{2} S ″ (s_{j}) Q^{(α_{j} - 2)} (s_{j}) \neq 0_{n \times n}

Proof: Similar to the proof of Proposition 9.1. □

Remark 9.3: Not all generalized eigenvalue problems admit a Jordan chain. For instance, considering the case of a double root, we obtain from (68), by multiplying from the left with $w_{0}^{H}$ and $w_{1}^{H}$ ( $H$ means Hermitian transpose, see Appendix A.1), respectively:

{\begin{matrix} w_{1}^{H} S' (s) w_{1} + w_{1}^{H} S (s) w_{0} = 0_{n \times 1} \\ w_{0}^{H} S (s) w_{1} = 0_{n \times 1}, w_{1} \neq 0_{n \times 1} \end{matrix}

(73)

By taking the Hermitian transpose of the last expression one obtains $w_{1}^{H} S (s^{*})^{T} w_{0} = 0_{n \times 1}$ and from this and (73)₁ one gets

w_{1}^{H} S' (s) w_{1} + w_{1}^{H} S (s) w_{0} = w_{1}^{H} S' (s) w_{1} + w_{1}^{H} (S (s) - S (s^{*})^{T}) w_{0} = 0_{n \times 1}

where

S (s) - S (s^{*})^{T} = M (s^{2} - s^{* 2}) + As - A^{T} s^{*} + B - B^{T}

and then

2 w_{1}^{H} M w_{1} s + w_{1}^{H} A w_{1} + w_{1}^{H} M w_{0} (s^{2} - s^{* 2}) + w_{1}^{H} A w_{0} s - w_{1}^{H} A^{T} w_{0} s^{*} + w_{1}^{H} (B - B^{T}) w_{0} = 0_{n \times 1}

Now assume that $A = 0_{nxn}$ and $B = B^{T}$ , then

2 w_{1}^{H} M w_{1} s + w_{1}^{H} M w_{0} (s^{2} - s^{* 2}) = 0

(74)

If $B$ is positive definite then, according to (68)₂:

w_{1}^{H} S (s) w_{1} = w_{1}^{H} M w_{1} s^{2} + w_{1}^{H} B w_{1} = 0_{n \times 1}, w_{1} \neq 0_{n \times 1} \Rightarrow s^{2} = - \frac{w_{1}^{H} B w_{1}}{w_{1}^{H} M w_{1}} \Rightarrow s = \pm i \sqrt{\frac{w_{1}^{H} B w_{1}}{w_{1}^{H} M w_{1}}}

and from this we may conclude that $s^{2} - s^{* 2} = 0$ and then, according to (74), $2 w_{1}^{H} M w_{1} s = 0 \Leftrightarrow w_{1} = 0_{n \times 1}$ , since $det B \neq 0$ and consequently $s \neq 0$ . Thus, a Jordan chain does not exist and Equations (68) are reduced to

S (s) w_{0} = 0_{n \times 1}, w_{0} \neq 0_{n \times 1}

(75)

For $A = 0_{nxn}$ and $B = B^{T}$ positive definite, the same argument will work for roots of any multiplicity $m$ , which means that $w_{1} = \dots = w_{m - 1} = 0_{n \times 1}$ and there exists $w_{0} \neq 0_{n \times 1}$ satisfying (70). In this case one may show that there will be $m$ linearly independent vectors $w_{0}, w_{0, 1}, \dots, w_{0, m - 1}$ satisfying equation (75).■

The dynamic stiffness $S$ is said to be simple if no eigenvalue is in defect (see Lancaster [10]). In this case we have $ν_{i} = m_{i}, i = 1, \dots, l$ . This is equivalent to $rank S (s_{i}) = n - m_{i}$ . It follows that if the dynamic stiffness has simple eigenvalues then it is simple. If there is an eigenvalue $s_{j}$ such that $ν_{j} < m_{j}$ , then $S$ is said to be defective. The dynamic stiffness is defective if and only if there is an eigenvalue $s_{j}$ with an eigenvector $w_{j}$ such that $w^{T} S' (s_{j}) w_{j} = 0$ for all $w$ satisfying $S (s_{j}) w = 0_{n \times 1}$ (see Lancaster [10], (Theorem 4.6)).

Proposition 9.3 For a system with simple dynamic stiffness let $s_{i}, i = 1, \dots, l$ be distinct roots of the characteristic equation with multiplicities $m_{i}, i = 1, \dots, l$ and let $w_{i, 0}, \dots, w_{i, m_{i} - 1}$ be a basis for $\ker (S (s_{i})), i = 1, \dots, l$ . Then the functions $w_{1, 0} e^{s_{1} t}, \dots, w_{1, m_{1} - 1} e^{s_{1} t}, \dots ., w_{l, 0} e^{s_{l} t}, \dots, w_{l, m_{l} - 1} e^{s_{l} t}$ constitute a basis for $ke r_{c} (S (d_{t}))$ and a general member of $ke r_{c} (S (d_{t}))$ may be written

z (t) = \sum_{i = 1}^{l} (\sum_{j = 0}^{m_{i} - 1} α_{ij} w_{i, j}) e^{s_{i} t}, α_{ij} \in C

(76)

Proof: Obviously $w_{i, 0} e^{s_{i} t}, \dots, w_{i, m_{i} - 1} e^{s_{i} t} \in ke r_{c} (S (d_{t})), i = 1, \dots, l$ . From Lemma A.1.1(a) it follows that

\sum_{i = 1}^{l} (\sum_{j = 0}^{m_{i} - 1} λ_{ij} w_{i, j}) e^{s_{i} t} = 0_{n \times 1} \Rightarrow \sum_{j = 0}^{m_{i} - 1} λ_{ij} w_{i, j} = 0_{n \times 1}, i = 1, \dots, l \Rightarrow λ_{ij} = 0

since $w_{i, 0}, \dots, w_{i, m_{i} - 1}$ are linearly independent. Thus, the $2 n$ functions $w_{1, 0} e^{s_{i} t}, \dots, w_{1, m_{1} - 1} e^{s_{1} t}, \dots .$ $, w_{l, 0} e^{s_{l} t}, \dots, w_{l, m_{l} - 1} e^{s_{l} t}$ constitute a basis for $ke r_{c} (S (d_{t}))$ . □

If the roots of the characteristic equation are complex, that is, if $s_{i} = σ_{i} + i ω_{i}, ω_{i} \neq 0,$ $i = 1, \dots, k$ , $l = 2 k$ , then (76) may be written

z (t) = \sum_{i = 1}^{k} ((\sum_{j = 0}^{m_{i} - 1} β_{ij} w_{i, j}) e^{s_{i} t} + (\sum_{j = 0}^{m_{i} - 1} γ_{ij} w_{i, j}^{*}) e^{s_{i}^{*} t})

(77)

The function (77) is real valued if and only if $β_{ij}^{*} = γ_{ij}$ and in this case

z (t) = x (t) = \sum_{i = 1}^{k} ((\sum_{j = 0}^{m_{i} - 1} β_{ij} w_{i, j}) e^{s_{i} t} + (\sum_{j = 0}^{m_{i} - 1} β_{ij} w_{i, j})^{*} e^{s_{i}^{*} t})

(78)

To find the weakest conditions on the coefficient matrices $A$ and $B$ , resulting in a system with simple dynamical stiffness, is an interesting problem. We have looked into the work of Lancaster and co-workers but not found a definitive solution.

From the previous discussion it may be inferred that the solution of the generalized eigenvalue problem, in the general case, is non-trivial. Methods based on the adjoint and reduce adjoint matrices are accounted for in this paper. An alternative approach to the solution of Equation (43) is based on the Smith canonical form for the dynamic stiffness (see Theorem A.1.1). Based on the Smith form the differential equation (43) may be written $S (d_{t}) x =$ $E (d_{t}) I (d_{t}) F (d_{t}) x = E (d_{t}) I (d_{t}) y = f$ , where $y = F (d_{t}) x$ and then $I (d_{t}) y =$ $E (d_{t})^{- 1} f = g$ or explicitly

(\begin{matrix} i_{1} (d_{t}) \\ ⋱ \\ i_{n} (d_{t}) \end{matrix}) (\begin{matrix} y_{1} \\ ⋮ \\ y_{n} \end{matrix}) = (\begin{matrix} g_{1} \\ ⋮ \\ g_{n} \end{matrix})

This system uncouples into $n$ independent scalar differential equations $i_{k} (d_{t}) y_{k} = g_{k}$ , $k = 1, \dots, n$ and the solution of (43) is then given by $x = F (d_{t})^{- 1} y$ (see Gohberg et al. [15]).■

10. The dynamic flexibility

Given a dynamic stiffness $S (s)$ with $p (s) \neq 0$ , the corresponding Dynamic flexibility (Admittance) is defined by $F (s) = S (s)^{- 1}$ . Obviously:

F (s) = \frac{T (s)}{p (s)}

(79)

For the spectral representation of the Dynamic flexibility one needs the following.

Lemma 10.1 Let $s_{i}, i = 1, \dots, l$ be distinct roots of the characteristic equation with multiplicities $m_{i}, i = 1, \dots, l$ . For a simple dynamic stiffness there are linearly independent vectors $w_{i, 1}, \dots, w_{i, m_{i}} \in \ker (S (s_{i})), i = 1, \dots, l$ satisfying the orthogonality conditions

w_{i, j}^{T} S' (s_{i}) w_{i, k} = α_{i} δ_{jk}, j, k = 1, . ., m_{i}

(80)

where $α_{i} = w_{i, k}^{T} S' (s_{i}) w_{i, k} \neq 0, k = 1, \ddot{,} m_{i}, i = 1, \dots, l$ .

Proof: See Lancaster [10] (Theorem 4.5) or [17]. □

Remark 10.1: The vectors $w_{i, 1}, \dots, w_{i, m_{i}}$ may of course be normalized so that $α_{i} = 1, i = 1, \dots, l$ . We here prefer the flexibility offered by allowing $α_{i}$ to be any complex number not equal to zero.■

Remark 10.2: By introducing $x = {(\begin{matrix} x & y \end{matrix})}^{T} \in R^{2 n \times 1}$ , where $y = \overset{\cdot}{x}$ , Equation (43) may be written

A \overset{\cdot}{x} + B x = f

(81)

where

A = (\begin{matrix} A & M \\ M & 0_{n \times n} \end{matrix}), B = (\begin{matrix} B & 0_{n \times n} \\ 0_{n \times n} & - M \end{matrix}), f = (\begin{matrix} f \\ 0_{n \times 1} \end{matrix})

Compare this with the proof of Theorem 6.1. Considering free vibrations we take $f = 0_{2 n \times 1}$ and $x (t) = z e^{st}, s \in C$ , where $z \in R^{2 n \times 1}$ is a constant vector. Then Equation (81) implies $(A s + B) z = 0_{2 n \times 1}$ with a solution $z \neq 0_{2 n \times 1}$ if and only if $det (A s + B) = (- 1)^{n} det M p (s) = 0$ . The matrix $A s + B$ is called the matrix pencil associated to the dynamic stiffness $S (s)$ (see Appendix A.1). If $x (t) = w e^{st}, s \in C$ then $y = \overset{\cdot}{x} = sw e^{st}$ and $z = {(\begin{matrix} w & sw \end{matrix})}^{T} \in R^{2 n \times 1}$ . Let $s_{i}$ be a root of the characteristic equation, that is, $p (s_{i}) = 0$ , with multiplicity $m_{i}$ . Consider vectors $w_{i, 1}, \dots, w_{i, m_{i}} \in \ker (S (s_{i}))$ and the corresponding vectors $z_{i, k} = {(\begin{matrix} w_{i, k} & s_{i} w_{1, k} \end{matrix})}^{T}, k = 1, \dots, m_{i}$ satisfying $(A s_{i} + B) z_{i, k} = 0_{2 n \times 1}$ , $k = 1, \dots, m_{i}$ . Note that the vectors $z_{i, 1}, \dots, z_{i, m_{i}}$ are linearly independent if and only if $w_{i, 1}, \dots, w_{i, m_{i}}$ are linearly independent. Now let $s_{j}$ be a second root with multiplicity $m_{j}$ and vectors $w_{j, 1}, \dots, w_{j, m_{j}} \in \ker (S (s_{j}))$ . Introduce the matrices

W_{i} = (\begin{matrix} w_{i, 1} & \dots & w_{i, m_{i}} \end{matrix}) \in R^{n \times m_{i}}, W_{j} = (\begin{matrix} w_{j, 1} & \dots & w_{j, m_{j}} \end{matrix}) \in R^{n \times m_{j}}

and

Z_{i} = (\begin{matrix} z_{i, 1} & \dots & z_{i, m_{i}} \end{matrix}) = (\begin{matrix} W_{i} \\ s_{i} W_{i} \end{matrix}) \in R^{2 n \times m_{i}}, Z_{j} = (\begin{matrix} z_{j, 1} & \dots & z_{j, m_{j}} \end{matrix}) = (\begin{matrix} W_{j} \\ s_{j} W_{j} \end{matrix}) \in R^{2 n \times m_{j}}

It may be demonstrated that $S (s)$ is simple if and only if the associated pencil, $A s + B$ , is a simple pencil (see Lancaster [10]). Thus, if $S (s)$ is simple then, according to Theorem A.1.2, it is possible to choose vectors $w_{i, 1}, \dots, w_{i, m_{i}} \in \ker (S (s_{i}))$ and $w_{j, 1}, \dots, w_{j, m_{j}} \in \ker (S (s_{j}))$ to obtain

Z_{i}^{T} A Z_{j} = (\begin{matrix} W_{i}^{T} & s_{i} W_{i}^{T} \end{matrix}) (\begin{matrix} A & M \\ M & 0_{n \times n} \end{matrix}) (\begin{matrix} W_{j} \\ s_{j} W_{j} \end{matrix}) = W_{i}^{T} (A + (s_{i} + s_{j}) M) W_{j} = {\begin{matrix} 0_{m_{i} \times m_{j}} if s_{i} \neq s_{j} \\ 1_{m_{i} \times m_{i}} if i = j \end{matrix}

Thus

W_{i}^{T} (A + (s_{i} + s_{j}) M) W_{j} = 0_{m_{i} \times m_{j}}, if s_{i} \neq s_{j}, W_{i}^{T} (A + 2 s_{i} M) W_{i} = 1_{m_{i} \times m_{i}}

The conditions above may be written

W_{i}^{T} (S (s_{i}) - S (s_{j})) W_{j} = 0_{m_{i} \times m_{j}}, W_{i}^{T} S' (s_{i}) W_{i} = 1_{m_{i} \times m_{i}}

(82)

respectively, and (82)₂ is equivalent to the orthogonality conditions (80), taking $α_{i} = 1, i = 1, \dots, l$ . Furthermore

Z_{i}^{T} B Z_{j} = (\begin{matrix} W_{i}^{T} & s_{i} W_{i}^{T} \end{matrix}) (\begin{matrix} B & 0_{n \times n} \\ 0_{n \times n} & - M \end{matrix}) (\begin{matrix} W_{j} \\ s_{j} W_{j} \end{matrix}) = W_{i}^{T} (B - s_{i} s_{j} M) W_{j} = {\begin{matrix} 0_{m_{i} \times m_{j}} if s_{i} \neq s_{i} \\ - s_{i} 1_{m_{i} \times m_{i}} if i = j \end{matrix}

Thus

W_{i}^{T} (B - s_{i} s_{j} M) W_{j} = 0_{m_{i} \times m_{j}}, s_{i} \neq s_{j}, W_{i}^{T} (B - s_{i}^{2} M) W_{i} = - s_{i} 1_{m_{i} \times m_{i}}

(83)

The orthogonality conditions (82)₁ and (83)₁ should be compared with the ones presented in Corollary 8.3 that were obtained for $A, B \in sym (R^{n \times n})$ .■

Now, introduce the matrices

W_{i} = \sum_{j = 1}^{m_{i}} w_{i, j} w_{i, j}^{T}, i = 1, \dots, l

(84)

The matrix $W_{i}$ with $rank (W_{i}) = m_{i}$ is called the residue at $s_{i}$ of the flexibility $F$ . The following theorem expresses the dynamic flexibility in terms of modal properties.

Theorem 10.1 If the dynamic stiffness is simple then

F (s) = \sum_{i = 1}^{l} \frac{W_{i}}{(s - s_{i}) α_{i}}

(85)

Proof: See Lancaster [10] (Theorem 4.3) and [17]. □

Remark 10.3: If the dynamic stiffness is simple then, according to (85), the flexibility has a pole of order one at $s_{i}, i = 1, \dots, l$ .■

Remark 10.4: The dynamic stiffness $S$ is simple for many multibody systems encountered in applications. For instance, if $G = 0_{n \times n}$ and $H = 0_{n \times n}$ (non-gyroscopic and non-circulatory) and if $C M^{- 1} K = K M^{- 1} C$ then $S (s) = M s^{2} + Cs + K$ is simple.■

The following is a direct consequence of the previous theorem.

Corollary 10.1 If the dynamic stiffness is simple and in the case of simple roots we have

F (s) = \sum_{i = 1}^{2 n} \frac{w_{i} w_{i}^{T}}{(s - s_{i}) α_{i}}

(86)

where $α_{i} = w_{i}^{T} S' (s_{i}) w_{i} \neq 0, i = 1, \dots, 2 n$ .

Proof: In this case $W_{i} = w_{i} w_{i}^{T}, i = 1, \dots, 2 n$ . □

The adjoint $T = T (s)$ may be expressed in terms of the modal properties of the dynamic stiffness.

Corollary 10.2 If the dynamic stiffness is simple then

T (s) = det M \sum_{i = 1}^{l} \frac{Π_{j = 1}^{l} {(s - s_{j})}^{m_{j}}}{(s - s_{i}) α_{i}} W_{i}

Proof: We have, according to (79) and Theorem 10.1:

T (s) = p (s) F (s) = \sum_{i = 1}^{l} \frac{p (s)}{(s - s_{i}) α_{i}} W_{i} = det M \sum_{i = 1}^{l} \frac{Π_{j = 1}^{l} {(s - s_{j})}^{m_{j}}}{(s - s_{i}) α_{i}} W_{i}

□

The dynamic flexibility may be expressed in terms of the adjoint according to the following.

Proposition 10.1 If the dynamic stiffness is simple then

\frac{W_{i}}{α_{i}} = \frac{m_{i} T^{(m_{i} - 1)} (s_{i})}{p^{(m_{i})} (s_{i})}, i = 1, \dots, l

(87)

and

F (s) = \sum_{i = 1}^{l} \frac{1}{s - s_{i}} \frac{m_{i} T^{(m_{i} - 1)} (s_{i})}{p^{(m_{i})} (s_{i})}

(88)

Proof: We have

T^{(m_{i} - 1)} (s_{i}) = det M \sum_{j = 1}^{l} {\frac{d^{m_{i} - 1}}{d s^{m_{i} - 1}} (\frac{Π_{j = 1}^{l} {(s - s_{j})}^{m_{j}}}{(s - s_{j}) α_{j}}) |}_{s = s_{i}} W_{j} = \frac{det M}{α_{i}} (m_{i} - 1)! Π_{\binom{j = 1}{j \neq i}}^{k} (s_{i} - s_{j}) W_{i}

and $p^{(m_{i})} (s_{i}) = det M m_{i}! Π_{\binom{j = 1}{j \neq i}}^{l} {(s - s_{j})}^{m_{j}}$ . Combining these two expressions leads to

W_{i} = \frac{T^{(m_{i} - 1)} (s_{i})}{\frac{det M}{α_{i}} (m_{i} - 1)! Π_{\binom{j = 1}{j \neq i}}^{k} (s_{i} - s_{j})} = \frac{α_{i} m_{i} T^{(m_{i} - 1)} (s_{i})}{p^{(m_{i})} (s_{i})}

and this proves the proposition. □

Corollary 10.3 If the dynamic stiffness is simple, and in the case of simple roots, we have

\frac{W_{i}}{α_{i}} = \frac{T (s_{i})}{p' (s_{i})}

If the dynamic stiffness is not simple then the situation is much more complicated. A representation of the flexibility in terms of modal properties, similar to the one in (85), is not known to the author. However, we have the following general result.

Theorem 10.2 Given a dynamic stiffness $S = S (s)$ . If $s_{i}, i = 1, \dots, l$ are roots of the characteristic equation with multiplicities $m_{i}$ then the flexibility may be written

F (s) = \sum_{i = 1}^{l} \sum_{j = 1}^{m_{i}} \frac{R_{ij}}{{(s - s_{i})}^{j}}

(89)

where the residues $R_{ij} \in R^{n \times n}$ are given by

R_{ij} = \frac{1}{(m_{i} - j)!} lim_{s \to s_{i}} \frac{d^{m_{i} - j}}{d s^{m_{i} - j}} ((s - s_{i})^{m_{i}} \frac{T (s)}{p (s)})

(90)

Proof: A partial fraction decomposition of the dynamic flexibility gives

F (s) = \frac{T (s)}{p (s)} = \frac{T (s)}{det M {(s - s_{1})}^{m_{1}} {(s - s_{2})}^{m_{2}} \dots {(s - s_{l})}^{m_{l}}} = \sum_{a = 1}^{l} \sum_{b = 1}^{m_{i}} \frac{R_{ab}}{{(s - s_{a})}^{b}}

where $T (s) = B_{2 (n - 1)} s^{2 (n - 1)} + \dots + B_{k} s^{k} + \dots + B_{0}$ and $R_{ab} \in R^{n \times n}$ are constant matrices. Then

\frac{d^{m_{i} - j}}{d s^{m_{i} - j}} ((s - s_{i})^{m_{i}} \frac{T (s)}{p (s)}) = \sum_{a = 1}^{l} \sum_{b = 1}^{m_{a}} \frac{d^{m_{i} - j}}{d s^{m_{i} - j}} \frac{{(s - s_{i})}^{m_{i}}}{{(s - s_{a})}^{b}} R_{ab} = G (s) + H (s)

where

G (s) = \sum_{\binom{a = 1}{a \neq ì}}^{l} \sum_{b = 1}^{m_{a}} \frac{d^{m_{i} - j}}{d s^{m_{i} - j}} \frac{{(s - s_{i})}^{m_{i}}}{{(s - s_{a})}^{b}} R_{ab}, H (s) = \sum_{b = 1}^{m_{i}} \frac{d^{m_{i} - j} {(s - s_{i})}^{m_{i} - b}}{d s^{m_{i} - j}} R_{ib}

Now $G (s_{i}) = 0_{n \times n}$ . Furthermore

\frac{d^{m_{i} - j} {(s - s_{i})}^{m_{i} - b}}{d s^{m_{i} - j}} = {\begin{matrix} 0 if j < b \\ (m_{i} - b) (m_{i} - b - 1) \dots (j - b + 1) (s - s_{i})^{j - b} if j \geq b \end{matrix}

Thus $H (s_{i}) = (m_{i} - j)! R_{ij}$ and

lim_{s \to s_{i}} \frac{d^{m_{i} - j}}{d s^{m_{i} - j}} ((s - s_{i})^{m_{i}} \frac{T (s)}{p (s)}) = G (s_{i}) + H (s_{i}) = (m_{i} - j)! R_{ij}

and this proves the theorem. □

Remark 10.5: For a simple dynamic stiffness one has $R_{ij} = 0, j = 2, \dots, m_{i}, i = 1, \dots, l$ and $R_{i 1} = \frac{W_{i}}{α_{i}},$ $i = 1, \dots, l$ . Note that in the case with simple roots we have, according to (90):

R_{i 1} = \frac{1}{(1 - 1)!} lim_{s \to s_{i}} ((s - s_{i}) \frac{T (s)}{p (s)}) = \frac{T (s_{i})}{p' (s_{i})}

which is in agreement with Corollary 10.3.■

Remark 10.6: A general formula relating the residues $R_{ij}$ to the spectral properties (of the generalized eigenvalue problem) of a non-simple dynamic stiffness is not known to the author. This problem is, according to the theorem above, related to the representation of the adjoint in terms of modal properties. The reader is referred to Gohberg et al. [15] where some results related to this problem are presented.■

11. Qualitative properties of damping factor and natural frequency

The eigenvalues $s = σ + i ω, σ, ω \in R$ are determined as roots to the characteristic equation. In this section we will derive some qualitative properties of $s$ that follow from physically motivated assumptions on the matrices $A$ and $B$ . This may be done without explicitly solving the characteristic equation. If $p (s) = 0$ , then there exists $w = u + iv, u, v \in R^{n \times 1}$ , $w \neq 0_{n \times 1}$ such that $S (s) w = (M s^{2} + As + B) w = 0_{n \times 1}$ . By multiplying this equation, from the left, with $w^{H}$ one obtains the equation

m s^{2} + as + b = 0

(91)

where

m = w^{H} Mw = u^{T} Mu + v^{T} Mv + i (u^{T} Mv - v^{T} Mu) = u^{T} Mu + v^{T} Mv

a = w^{H} Aw = u^{T} Au + v^{T} Av + i (u^{T} Av - v^{T} Au)

(92)

b = w^{H} Bw = u^{T} Bu + v^{T} Bv + i (u^{T} Bv - v^{T} Bu)

We have the decompositions $A = C + G$ , $B = K + H$ , where $C^{T} = C$ , $G^{T} = - G$ , $K^{T} = K$ and $H^{T} = - H$ . Inserting this into (92)_2,3 results in

a = u^{T} Cu + v^{T} Cv + i 2 u^{T} Gv, b = u^{T} Ku + v^{T} Kv + i 2 u^{T} Hv

If we introduce the real constants

c = u^{T} Cu + v^{T} Cv, g = u^{T} Gv, k = u^{T} Ku + v^{T} Kv, h = u^{T} Hv

(93)

then

a = c + i 2 g, b = k + i 2 h

(94)

The numerical values of the real constants $c, g, k, h$ will of course depend on the actual mode vector $w = u + iv$ . However, some general properties of these constants may be motivated from physical considerations.

Proposition 11.1 If the damping matrix $C$ is positive semi-definite then $c \geq 0$ and if $C$ is positive definite then $c > 0$ . If the stiffness matrix $K$ is positive semi-definite then $k \geq 0$ and if $K$ is positive definite then $k > 0$ .

Proof: This is a direct consequence of (93)₁ and (93)₃, respectively, and the fact that $w = u + iv \neq 0_{n \times 1}$ . □

Lemma 11.1 If the stiffness matrix $K$ is positive definite then $det B \neq 0$ and consequently $p (0) = det B \neq 0$ .

Proof: Assume that $det B = 0$ then there exists $w \neq 0_{n \times 1}$ such that $Bw = 0$ and then $w^{H} Bw =$ $b = u^{T} Ku + v^{T} Kv + i 2 u^{T} Hv = 0 \Rightarrow u^{T} Ku + v^{T} Kv = 0$ , which leads to a contradiction. □

Note that if $p (s) = 0$ then $p (s^{*}) = 0$ and we have $S (s^{*}) w^{*} = (M s^{* 2} + A s^{*} + B) w^{*} = 0_{n \times 1}$ where $w^{*} = u - iv, u, v \in R^{n \times 1}$ and

m s^{* 2} + a^{*} s^{*} + b^{*} = 0

which is equal to the complex conjugation of Equation (91). By multiplying (91) by $s^{*}$ one obtains

m {| s |}^{2} s + a {| s |}^{2} + b s^{*} = m {| s |}^{2} (σ + i ω) + (c + i 2 g) {| s |}^{2} + (k + i 2 h) (σ - i ω) = 0

which is equivalent to

{\begin{matrix} (m {| s |}^{2} + k) σ + c {| s |}^{2} + 2 h ω = 0 \\ (m {| s |}^{2} - k) ω + 2 g {| s |}^{2} + 2 h σ = 0 \end{matrix}

(95)

Proposition 11.2 If $h = 0$ and if $C$ and $K$ are positive definite then $σ < 0$ .

Proof: From (95)₁ it follows that $(m {| s |}^{2} + k) σ + c {| s |}^{2} = 0$ where, according to Lemma 11.1, ${| s |}^{2} > 0$ . Then, since $c > 0$ and $k > 0$ :

σ = - \frac{c {| s |}^{2}}{m {| s |}^{2} + k} < 0

Proposition 11.3 (Gyroscopic conditions) If $h = 0$ , $c = 0$ and $k > 0$ then

σ = 0, ω = 0 or ω = - \frac{g}{m} \pm \sqrt{\frac{k}{m} + {(\frac{g}{m})}^{2}}

(96)

Proof: From (95)₁ it follows that $(m {| s |}^{2} + k) σ = 0 \Rightarrow σ = 0$ and then from (95)₂

(m ω^{2} - k) ω + 2 g ω^{2} = 0 \Leftrightarrow ω (ω^{2} + \frac{2 g}{m} ω - \frac{k}{m}) = 0

and this proves the proposition. □

Remark 11.1: Note that if $K$ is positive definite then, according to Lemma 11.1, $ω = 0$ in (96) is excluded.■

Using (94), Equation (91) may be written

m s^{2} + (c + i 2 g) s + k + i 2 h = 0

(97)

Proposition 11.4 Equation (97) has the solutions $s = σ + i ω$ where

σ = - \frac{c}{2 m} \pm \sqrt{\frac{\sqrt{Δ^{2} + Γ^{2}} + Δ}{2}}, ω = {\begin{matrix} - \frac{g}{m} \pm \sqrt{\frac{\sqrt{Δ^{2} + Γ^{2}} - Δ}{2}}, Γ \geq 0 \\ - \frac{g}{m} \mp \sqrt{\frac{\sqrt{Δ^{2} + Γ^{2}} - Δ}{2}}, Γ < 0 \end{matrix}

(98)

and

Δ = \frac{c^{2}}{4 m^{2}} - \frac{g^{2}}{m^{2}} - \frac{k}{m}, Γ = \frac{cg}{m^{2}} - \frac{2 h}{m}

(99)

Proof: We have

\begin{matrix} s = σ + i ω = - \frac{c + i 2 g}{2 m} \pm \sqrt{{(\frac{c + i 2 g}{2 m})}^{2} - \frac{k + i 2 h}{m}} = \\ - \frac{c}{2 m} - i \frac{g}{m} \pm \sqrt{\frac{c^{2}}{4 m^{2}} - \frac{g^{2}}{m^{2}} - \frac{k}{m} + i (\frac{cg}{m^{2}} - \frac{2 h}{m})} = - \frac{c}{2 m} - i \frac{g}{m} \pm \sqrt{Δ + i Γ} \end{matrix}

where

\sqrt{Δ + i Γ} = {\begin{matrix} \sqrt{\frac{\sqrt{Δ^{2} + Γ^{2}} + Δ}{2}} + i \sqrt{\frac{\sqrt{Δ^{2} + Γ^{2}} - Δ}{2}}, Γ \geq 0 \\ \sqrt{\frac{\sqrt{Δ^{2} + Γ^{2}} + Δ}{2}} - i \sqrt{\frac{\sqrt{Δ^{2} + Γ^{2}} - Δ}{2}}, Γ < 0 \end{matrix}

and this proves the proposition. □

Remark 11.2: Note that we cannot make any definitive decision on the signs in (98). These will depend on the mode shape vector $w$ .■

Proposition 11.5 (Stability) (a) If

c \geq 0, k \leq - \frac{c^{2}}{4 m} - \frac{g^{2}}{m}

(100)

then $σ > 0$ . (b) If

c \geq 0, k > - \frac{c^{2}}{4 m} - \frac{g^{2}}{m}

(101)

then $σ < 0$ if and only if

(kc + 4 gh) c > 4 h^{2} m

(102)

and $σ = 0$ if and only if $(kc + 4 gh) c = 4 h^{2} m$ .

Proof: (a) If $c \geq 0$ then we have the following equivalence:

σ = - \frac{c}{2 m} + \sqrt{\frac{\sqrt{Δ^{2} + Γ^{2}} + Δ}{2}} \geq 0 \Leftrightarrow \frac{\sqrt{Δ^{2} + Γ^{2}} + Δ}{2} \geq \frac{c^{2}}{4 m^{2}}

(103)

From (99)₁ and (100)₁ one obtains

Δ = \frac{c^{2}}{4 m^{2}} - \frac{g^{2}}{m^{2}} - \frac{k}{m} \geq \frac{c^{2}}{4 m^{2}} - \frac{g^{2}}{m^{2}} + \frac{c^{2}}{4 m^{2}} + \frac{g^{2}}{m^{2}} = \frac{c^{2}}{2 m^{2}}

and then

\frac{\sqrt{Δ^{2} + Γ^{2}} + Δ}{2} = \frac{\sqrt{Δ^{2} + Γ^{2}}}{2} + \frac{Δ}{2} > \frac{c^{2}}{4 m^{2}}

and consequently $σ > 0$ . (b) We have the implications

\begin{matrix} σ \leq 0 \Leftrightarrow \frac{\sqrt{Δ^{2} + Γ^{2}} + Δ}{2} \leq \frac{c^{2}}{4 m^{2}} \Leftrightarrow \sqrt{Δ^{2} + Γ^{2}} \leq \frac{c^{2}}{2 m^{2}} - Δ \Rightarrow Δ^{2} + Γ^{2} \leq (\frac{c^{2}}{2 m^{2}} - Δ)^{2} = \\ \frac{c^{4}}{4 m^{4}} + Δ^{2} - \frac{c^{2}}{m^{2}} Δ \Leftrightarrow \frac{c^{4}}{4 m^{4}} - \frac{c^{2}}{m^{2}} Δ \geq Γ^{2} \Leftrightarrow \frac{k}{m^{3}} c^{2} + 4 \frac{gh}{m^{3}} c \geq \frac{4 h^{2} m}{m^{3}} \Leftrightarrow k c^{2} + 4 ghc \geq 4 h^{2} m \end{matrix}

It remains to show the implication

\sqrt{Δ^{2} + Γ^{2}} \leq \frac{c^{2}}{2 m^{2}} - Δ \Leftarrow Δ^{2} + Γ^{2} \leq (\frac{c^{2}}{2 m^{2}} - Δ)^{2}

But

Δ^{2} + Γ^{2} \leq (\frac{c^{2}}{2 m^{2}} - Δ)^{2} \Rightarrow \sqrt{Δ^{2} + Γ^{2}} \leq | \frac{c^{2}}{2 m^{2}} - Δ | = \frac{c^{2}}{2 m^{2}} - Δ if \frac{c^{2}}{2 m^{2}} - Δ > 0

This condition being equivalent to

\frac{c^{2}}{2 m^{2}} - Δ > 0 \Leftrightarrow \frac{c^{2}}{2 m^{2}} - (\frac{c^{2}}{4 m^{2}} - \frac{g^{2}}{m^{2}} - \frac{k}{m}) = \frac{c^{2}}{4 m^{2}} + \frac{g^{2}}{m^{2}} + \frac{k}{m} > 0 \Leftrightarrow k > - \frac{c^{2}}{4 m} - \frac{g^{2}}{m}

σ = - \frac{c}{2 m} + \sqrt{\frac{\sqrt{Δ^{2} + Γ^{2}} + Δ}{2}} > 0

This proves the proposition. □

Remark 11.3: (Instability) From the previous proof it follows that if $c \geq 0$ then

(kc + 4 gh) c < 4 h^{2} m \Rightarrow σ > 0

■

Remark 11.4: It follows from the previous proposition that if $c = 0$ , $k > - \frac{g^{2}}{m}$ and $h \neq 0$ then $σ > 0$ . Thus a gyroscopic system may be destabilized by introducing circulatory forces.■

Corollary 11.1 If inequality (102) is satisfied, $h = 0$ (non-circulatory conditions) and

(a) $c > 0$ then $σ \leq 0$ if and only if $k \geq 0$ ;

(b) $c = 0$ then $σ \geq 0$ .

Proof: (a) From (102) it follows that if $c > 0$ then $k c^{2} \geq 0 \Leftrightarrow k \geq 0$ . (b) If $c = 0$ and $σ < 0$ then inequality (102) leads to a contradiction. Thus $σ \geq 0$ . □

Remark 11.5: Note that (101) is valid if $c \geq 0$ , $g \neq 0$ and $- \frac{g^{2}}{m} < \frac{k}{m} < 0$ . This suggests that the damping term may have a destabilizing influence on a system where the potential energy has a maximum at $x = 0_{n \times 1}$ . A closer analysis will confirm this. See Example 11.1 below and Krechetnikov and Marsden [18]■

Example 11.1: Consider a particle with mass $m_{0}$ connected by linear elastic (massless) springs to a circular ring that is rotating with a given angular velocity $Ω$ (see Figure 3). Let $x^{1}, x^{2}$ denote Cartesian configuration coordinates of the particle relative to the ring and put $x = {(\begin{matrix} x^{1} & x^{2} \end{matrix})}^{T}$ . The elastic energy and the dissipation function of the system are assumed to be given by

V_{e} (x) = \frac{1}{2} m_{0} ω_{0}^{2} x^{T} x, D_{e} (x) = 2 m_{0} ω_{0} ζ {\overset{\cdot}{x}}^{T} \overset{\cdot}{x}

where $m_{0}$ , $ω_{0}$ and $ζ$ are positive constants.

Figure 3.

Two-dimensional multibody system.

The dynamic stiffness matrix of the system $S (s) = M s^{2} + As + B$ where

M = m_{0} 1_{2 \times 2}, A = 2 m_{0} ω_{0} ζ 1_{2 \times 2} + 2 m_{0} Ω (\begin{matrix} 0 & - 1 \\ 1 & 0 \end{matrix}), B = m_{0} (ω_{0}^{2} - Ω^{2}) 1_{2 \times 2}

Thus $C = 2 m_{0} ω_{0} ζ 1_{2 \times 2}$ , $G = 2 m_{0} Ω (\begin{matrix} 0 & - 1 \\ 1 & 0 \end{matrix})$ , $K = m_{0} (ω_{0}^{2} - Ω^{2}) 1_{2 \times 2}$ and $H = 0_{2 \times 2}$ and consequently

\begin{matrix} m = (u^{T} u + v^{T} v) m_{0}, c = (u^{T} u + v^{T} v) 2 m_{0} ω_{0} ζ, g = 2 m_{0} Ω u^{T} (\begin{matrix} 0 & - 1 \\ 1 & 0 \end{matrix}) v = 2 m_{0} Ω (u^{2} v^{1} - u^{1} v^{2}), \\ k = (u^{T} u + v^{T} v) m_{0} (ω_{0}^{2} - Ω^{2}), h = 0 \end{matrix}

The condition (101)₂ may be written

(u^{T} u + v^{T} v) m_{0} (ω_{0}^{2} - Ω^{2}) > - \frac{{(u^{T} u + v^{T} v)}^{2} 4 m_{0}^{2} ω_{0}^{2} ζ^{2}}{4 m_{0}} - \frac{4 m_{0}^{2} Ω^{2} {(u^{2} v^{1} - u^{1} v^{2})}^{2}}{m_{0}}

Without loss of generality we may choose $w^{H} w = u^{T} u + v^{T} v = 1$ . Thus

ω_{0}^{2} - Ω^{2} > - ω_{0}^{2} ζ^{2} - 4 Ω^{2} (u^{2} v^{1} - u^{1} v^{2})^{2} \Leftrightarrow ω_{0}^{2} (1 + ζ^{2}) > Ω^{2} (1 - 4 (u^{2} v^{1} - u^{1} v^{2})^{2})

which is satisfied if

ω_{0}^{2} (1 + ζ^{2}) > Ω^{2}

since $- 1 \leq 1 - 4 (u^{2} v^{1} - u^{1} v^{2})^{2} = 1 - 4 {| u \times v |}^{2} \leq 1$ . The left inequality follows from $1 - 4 {| u \times v |}^{2} \geq$ $1 - 4 {| u |}^{2} {| v |}^{2} \geq 1 - 4 \frac{1}{2} = - 1$ , where we have used the fact that ${| u |}^{2} {| v |}^{2} \leq \frac{{| u |}^{2} + {| v |}^{2}}{2} = \frac{1}{2}$ . The system is thus stable if and only if $k = m_{0} (ω_{0}^{2} - Ω^{2}) > 0 \Leftrightarrow | Ω | < ω_{0}$ .■

The following corollaries are direct consequences of Proposition 11.4.

Corollary 11.2 (Non-gyroscopic and non-circulatory conditions) If $g = 0$ , $h = 0$ then

(a) if $\frac{c^{2}}{4 m^{2}} - \frac{k}{m} > 0$ : $σ = - \frac{c}{2 m} \pm \sqrt{\frac{c^{2}}{4 m^{2}} - \frac{k}{m}}$ , $ω = 0$ (over-damped)

(b) if $\frac{c^{2}}{4 m^{2}} - \frac{k}{m} < 0$ : $σ = - \frac{c}{2 m}$ , $ω = \pm \sqrt{\frac{k}{m} - \frac{c^{2}}{4 m^{2}}}$ (under-damped)

Corollary 11.3 (Gyroscopic conditions) If $c = 0$ , $h = 0$ then

\begin{matrix} (a) if k \geq - \frac{g^{2}}{m} : σ = 0, ω = - \frac{g}{m} \pm \sqrt{\frac{g^{2}}{m^{2}} + \frac{k}{m}} \\ (b) if k < - \frac{g^{2}}{m} : σ = \pm \sqrt{- \frac{g^{2}}{m^{2}} + \frac{k}{m}}, ω = - \frac{g}{m} \end{matrix}

(104)

Proof: If $c = 0$ , $h = 0$ then $Γ = 0$ and consequently, according to (98):

σ = \pm \sqrt{\frac{| Δ | + Δ}{2}}, ω = - \frac{g}{m} \pm \sqrt{\frac{| Δ | - Δ}{2}}

(105)

However, $Δ = - \frac{g^{2}}{m^{2}} - \frac{k}{m}$ and then (104) follows from (105). □

Remark 11.6: Properties of overdamped and gyroscopic systems are discussed by Duffin [19], Rogers [20], Inman [21], Brackwell and Lancaster [22] and Lancaster et al. [23].■

Corollary 11.4 (Non-damped and non-gyroscopic conditions) If $c = 0$ , $g = 0$ then

σ = \pm \sqrt{- \frac{k}{2 m} + \sqrt{\frac{k^{2}}{4 m^{2}} + \frac{h^{2}}{m^{2}}}}, ω = {\begin{matrix} \pm \sqrt{\frac{k}{2 m} + \sqrt{\frac{k^{2}}{4 m^{2}} + \frac{h^{2}}{m^{2}}}}, if h \leq 0 \\ \mp \sqrt{\frac{k}{2 m} + \sqrt{\frac{k^{2}}{4 m^{2}} + \frac{h^{2}}{m^{2}}}}, if h > 0 \end{matrix}

Proof: If $c = 0$ , $g = 0$ then, according to (98):

σ = \pm \sqrt{\frac{\sqrt{Δ^{2} + Γ^{2}} + Δ}{2}}, ω = {\begin{matrix} \pm \sqrt{\frac{\sqrt{Δ^{2} + Γ^{2}} - Δ}{2}}, Γ \geq 0 \\ \mp \sqrt{\frac{\sqrt{Δ^{2} + Γ^{2}} - Δ}{2}}, Γ < 0 \end{matrix}

where $Δ = - \frac{k}{m}$ , $Γ = - \frac{2 h}{m}$ and this proves the corollary. □

Remark 11.7: Let $G = 0_{n \times n}$ , $H = 0_{n \times n}$ and assume that $C$ and $K$ are positive definite. The system is then said to be over-damped if $c^{2} > 4 mk$ for all $w \in C^{n \times 1}, w \neq 0_{n \times 1}$ . Let $C = 0_{n \times n}$ , $H = 0_{n \times n}$ , $G \neq 0_{n \times n}$ and assume that $K$ is non-negative definite. The system is then said to be gyroscopic. It may be shown that if the system is over-damped or gyroscopic then the dynamic stiffness is simple (see Lancaster [10]).■

12. Forced vibrations

Consider the non-homogeneous linearized equation of motion

S (d_{t}) x = f

(106)

where the external force $f = f (t)$ is a prescribe function. The general solution, the response, is given by $x = x_{h} + x_{p}$ , where $x_{h} \in \ker (S (d_{t}))$ and $x_{p}$ is any function satisfying $S (d_{t}) x_{p} = f$ .

Now let $v_{p}$ be a solution to the inhomogeneous characteristic differential equation $p (d_{t}) v_{p} = f$ . If $v_{p}$ is sufficiently differentiable then put $x_{p} = T (d_{t}) v_{p}$ . The function $x_{p} = x_{p} (t)$ is then a solution to (106) since

S (d_{t}) x_{p} = S (d_{t}) T (d_{t}) v_{p} = p (d_{t}) v_{p} = f

Assume that the external force is harmonic $f (t) = f_{0} \sin ω t$ , where $f_{0} \in R^{n \times 1}$ is a constant vector and $ω \in R$ . We may complexify the problem by taking $f (t) = f_{0} e^{st}$ , where $s \in C$ and consider the equation

S (d_{t}) z_{p} = f_{0} e^{st}

(107)

If $z_{p}$ is a solution to (107) for $s = i ω$ then $x_{p} = Im (z_{p})$ is a solution to (106) with $f (t) = f_{0} \sin ω t$ . Assume that $z_{p} (t) = V e^{st}$ is a solution to (107), where $V \in C^{n \times 1}$ is a time-independent vector to be determined. Then

S (d_{t}) z_{p} = S (s) V e^{st} = f_{0} e^{st} \Rightarrow S (s) V = f_{0} \Rightarrow V = F (s) f_{0}

if $p (s) \neq 0$ . Using (89) one obtains

V = V (s) = \sum_{i = 1}^{l} \sum_{j = 1}^{m_{i}} \frac{R_{ij} f_{0}}{{(s - s_{i})}^{j}}

and then

z_{p} (t) = \sum_{i = 1}^{l} \sum_{j = 1}^{m_{i}} \frac{R_{ij} f_{0}}{{(s - s_{i})}^{j}} e^{st}

and, consequently, in the case with harmonic external force then

x_{p} (t) = \sum_{i = 1}^{l} \sum_{j = 1}^{m_{i}} Im (\frac{R_{ij}}{{(i ω - s_{i})}^{j}} e^{i ω t}) f_{0} = \sum_{i = 1}^{l} \sum_{j = 1}^{m_{i}} Im (\frac{R_{ij}}{{(- σ_{i} + i (ω - ω_{i}))}^{j}} e^{i ω t}) f_{0}

(108)

We have the following estimate:

‖ x_{p} (t) ‖ \leq \sum_{i = 1}^{l} \sum_{j = 1}^{m_{i}} \frac{‖ R_{ij} ‖ ‖ f_{0} ‖}{{(\sqrt{σ_{i}^{2} + {(ω - ω_{i})}^{2}})}^{j}}

If $ω = ω_{k}$ (resonance) and $σ_{k} \neq 0$ then the term with $i = k$ and $j = m_{k}$ on the right-hand side of the inequality reads

\frac{‖ R_{k m_{k}} ‖ ‖ f_{0} ‖}{{| σ_{k} |}^{m_{k}}}

(109)

and if $0 < | σ_{k} | << 1$ and $m_{k}$ is large then the denominator is a small number rendering the quotient large if the nominator $‖ R_{k m_{k}} ‖ ‖ f_{0} ‖$ is not too small. Thus, the amplitude of resonances in systems with multiple eigenvalues may become very large.

If the dynamic stiffness is simple then the dynamic flexibility is given in (85) and

z_{p} (t) = \sum_{i = 1}^{l} \frac{W_{i} f_{0}}{(s - s_{i}) α_{i}} e^{st}

In the case of a harmonic external force we then obtain

x_{p} (t) = \sum_{i = 1}^{l} Im (\frac{W_{i}}{(i ω - s_{i}) α_{i}} e^{i ω t}) f_{0}

(110)

We have the following estimate:

‖ x_{p} (t) ‖ \leq \sum_{i = 1}^{l} \frac{‖ W_{i} ‖ ‖ f_{0} ‖}{\sqrt{σ_{i}^{2} + {(ω - ω_{i})}^{2}} | α_{i} |}

where $‖ W_{i} ‖ \leq \sum_{j = 1}^{m_{i}} ‖ w_{i, j} w_{i, j}^{T} ‖ = \sum_{j = 1}^{m_{i}} ‖ w_{i, j} ‖, i = 1, \dots, l$ . If $ω = ω_{k}$ (resonance) and $σ_{k} \neq 0$ then the term with $i = k$ on the right-hand side of the inequality reads

\frac{‖ W_{k} ‖ ‖ f_{0} ‖}{| σ_{k} | | α_{k} |}

(111)

By comparing this expression with (109) it is clear that if the dynamic stiffness is simple then the resonance is not as prominent as it may be in the general case.

If the roots of the characteristic equation are simple then the dynamic flexibility is given by (86) and for a harmonic external force we then obtain the response

x_{p} (t) = \sum_{i = 1}^{2 n} Im (\frac{w_{i} w_{i}^{T}}{(s - s_{i}) α_{i}} e^{i ω t}) f_{0}

From the engineering point of view resonances are important, since they may give rise to a magnification of the response amplitude. In this perspective it is clear that multiple roots may form a greater risk of high amplitudes than simple roots, since the amplitude in (111) compares with

\frac{‖ w_{k} ‖ ‖ f_{0} ‖}{| σ_{k} | | α_{k} |}

in the case of simple roots.

Proposition 12.1 If the roots of the characteristic equation are simple and complex, that is if $s_{i} = σ_{i} + i ω_{i}, ω_{i} \neq 0, i = 1, \dots, n$ , then

F (s) = \sum_{i = 1}^{n} \frac{(α_{i}^{*} w_{i} w_{i}^{T} + α_{i} w_{i}^{*} w_{i}^{* T}) s + w_{i} w_{i}^{T} s_{i}^{*} α_{i}^{*} + w_{i}^{*} w_{i}^{* T} s_{i} α_{i}}{(s - s_{i}) (s - s_{i}^{*}) {| α_{i} |}^{2}}

(112)

and if the external force is harmonic then

x_{p} (t) = \sum_{i = 1}^{n} Im (\frac{(α_{i}^{*} w_{i} w_{i}^{T} + α_{i} w_{i}^{*} w_{i}^{* T}) s + w_{i} w_{i}^{T} s_{i}^{*} α_{i}^{*} + w_{i}^{*} w_{i}^{* T} s_{i} α_{i}}{(i ω - s_{i}) (i ω - s_{i}^{*}) {| α_{i} |}^{2}} e^{i ω t}) f_{0}

(113)

Proof: From (86) it follows that

F (s) = \sum_{i = 1}^{n} (\frac{w_{i} w_{i}^{T}}{(s - s_{i}) α_{i}} + \frac{w_{i}^{*} w_{i}^{* T}}{(s - s_{i}^{*}) α_{i}^{*}}) = \sum_{i = 1}^{n} \frac{(α_{i}^{*} w_{i} w_{i}^{T} + α_{i} w_{i}^{*} w_{i}^{* T}) s + w_{i} w_{i}^{T} s_{i}^{*} α_{i}^{*} + w_{i}^{*} w_{i}^{* T} s_{i} α_{i}}{(s - s_{i}) (s - s_{i}^{*}) {| α_{i} |}^{2}}

The complexified solution is given by

z_{p} (t) = \sum_{i = 1}^{n} \frac{(α_{i}^{*} w_{i} w_{i}^{T} + α_{i} w_{i}^{*} w_{i}^{* T}) s + w_{i} w_{i}^{T} s_{i}^{*} α_{i}^{*} + w_{i}^{*} w_{i}^{* T} s_{i} α_{i}}{(s - s_{i}) (s - s_{i}^{*}) {| α_{i} |}^{2}} f_{0} e^{st}

(114)

and the real solution (113) is obtained by putting $s = i ω$ into (114) and then selecting the imaginary part. □

Remark 12.1: In the case of real modes $w_{i} = u_{i} \in R^{n \times 1}$ and then, according to (112):

F (s) = \sum_{i = 1}^{n} \frac{((α_{i}^{*} + α_{i}) s - s_{i}^{*} α_{i}^{*} - s_{i} α_{i}) u_{i} u_{i}^{T}}{(s - s_{i}) (s - s_{i}^{*}) {| α_{i} |}^{2}}

where

α_{i} = u_{i} S' (s_{i}) u_{i}^{T} = u_{i} (2 M (σ_{i} + i ω_{i}) + A) u_{i}^{T} = 2 σ_{i} μ_{i} + η_{i} + i 2 ω_{i} μ_{i}

and $μ_{i} = u_{i} {Mu}_{i}^{T}$ , $η_{i} = u_{i} {Au}_{i}^{T}$ . If we assume that $2 σ_{i} μ_{i} + η_{i} = 0$ , which is the case for non-gyroscopic and non-circulatory systems, then $α_{i} = i 2 ω_{i} μ_{i}$ , $(α_{i}^{*} + α_{i}) s = 0$ , $s_{i} α_{i} = - 2 ω_{i}^{2} μ_{i} +$ $i 2 σ_{i} ω_{i} μ_{i}$ and $(α_{i}^{*} + α_{i}) s - s_{i}^{*} α_{i}^{*} - s_{i} α_{i} = 4 ω_{i}^{2} μ_{i}$ . Consequently

F (s) = \sum_{i = 1}^{n} \frac{u_{i} u_{i}^{T}}{(s - s_{i}) (s - s_{i}^{*}) μ_{i}}

(115)

■

According to Equation (46) the characteristic polynomial may be written $p (s) =$ $det M (s - s_{1})^{m_{1}} (s - s_{2})^{m_{2}} \dots (s - s_{l})^{m_{l}}$ , where all $s_{1}, \dots, s_{l}$ are distinct complex numbers and $\sum_{i = 1}^{l} m_{i} = 2 n$ . If all roots are complex $s_{i} = σ_{i} + i ω_{i}, ω_{i} \neq 0$ then $l = 2 k$ and we may write

p (s) = det M (s - s_{1})^{m_{1}} (s - s_{1}^{*})^{m_{1}} \dots (s - s_{k})^{m_{k}} (s - s_{k}^{*})^{m_{k}}

where $\sum_{i = 1}^{k} m_{i} = n$ . We end this section with the following theorem.

Theorem 12.1 If the roots of the characteristic equation are complex, that is if $s_{i} = σ_{i} + i ω_{i}, ω_{i} \neq 0$ , $i = 1, \dots, k$ , ( $2 k = l$ ) are distinct roots with multiplicities $m_{i}$ , then for a system with simple dynamic stiffness the solution to the initial value problem

M \ddot{x} + A \overset{\cdot}{x} + Bx = f_{0} \sin ω t

(116)

x (0) = x_{0}, \overset{\cdot}{x} (0) = {\overset{\cdot}{x}}_{0}

is given by

x (t) = \sum_{i = 1}^{k} ((\sum_{j = 0}^{m_{i} - 1} β_{ij} w_{i, j}) e^{s_{i} t} + (\sum_{j = 0}^{m_{i} - 1} β_{ij} w_{i, j})^{*} e^{s_{i}^{*} t}) + \sum_{i = 1}^{k} \frac{1}{2 i} (\frac{W_{i}}{(i ω - s_{i}) α_{i}} e^{i ω t} + \frac{W_{i}^{*}}{(i ω + s_{i}^{*}) α_{i}^{*}} e^{- i ω t}) f_{0}

(117)

where the $n$ complex constants $β_{ij}, i = 1, \dots, k j = 1, \dots, m_{i}$ are uniquely determined by

\sum_{i = 1}^{k} (\sum_{j = 0}^{m_{i} - 1} β_{ij} w_{i, j} + \sum_{j = 0}^{m_{i} - 1} β_{ij}^{*} w_{i, j}^{*}) + \sum_{i = 1}^{k} \frac{1}{2 i} (\frac{W_{i}}{(i ω - s_{i}) α_{i}} + \frac{W_{i}^{*}}{(i ω + s_{i}^{*}) α_{i}^{*}}) f_{0} = x_{0}

(118)

\sum_{i = 1}^{k} (s_{i} \sum_{j = 0}^{m_{i} - 1} β_{ij} w_{i, j} + s_{i}^{*} \sum_{j = 0}^{m_{i} - 1} β_{ij}^{*} w_{i, j}^{*}) + \frac{ω}{2} \sum_{i = 1}^{k} (\frac{W_{i}}{(i ω - s_{i}) α_{i}} - \frac{W_{i}^{*}}{(i ω + s_{i}^{*}) α_{i}^{*}}) f_{0} = {\overset{\cdot}{x}}_{0}

Proof: This is a direct consequence of (78) and (110). □

13. Concluding remarks

In this paper the linearized equations of motion in multibody dynamics have been derived. The multibody is assumed to consist of rigid parts as well as parts consisting of visco-elastic material. Parts are connected by ideal joints and it is assumed that an equilibrium configuration for the multibody system exists. The analysis of the linearized equations of motion is based on a consequent use of the differential operator related to the dynamic stiffness of the system and its adjoint. This approach seems to be novel, at least in the engineering literature. It has the advantage of retaining the original differential equation with its close connection to the basic principles of mechanics. The solution of the linearized equations involves the solution of a (generalized) eigenvalue problem involving a Jordan chain. It is demonstrated how the solution can be obtained from the adjoint of the stiffness matrix. Some seemingly open problems have been identified:

- to find the weakest conditions on the coefficient matrices $A$ and $B$ , resulting in a system with simple dynamical stiffness;

- a general formula relating the residues $R_{ij}$ to the spectral properties (the generalized eigenvalue problem) of a general (non-simple) dynamic stiffness.

A non-trivial extension of the multibody model presented in this paper might involve the incorporation of friction and visco-elastic bushing elements at joints. This will be addressed in a forthcoming paper.

Footnotes

Appendix

Funding

This research received no specific grant from any funding agency in the public, commercial or not-for-profit sectors.

Conflict of interest

None declared.

References

Lord Rayleigh. The theory of sound (two volumes), New York: Dover Publications, 1945.

Hughes

TJR

. The finite element method. Linear static and dynamic finite element analysis. Englewood Cliffs, NJ: Prentice Hall, Inc., 1987.

Géradin

Rixen

. Mechanical vibrations. Theory and applications to structural dynamics. Chichester: John Wiley & Sons, 1997.

Lidström

. On the equations of motion in multibody dynamics. Math Mech Solid 2012; 17: 165–205.

Birkhoff

Rota

G-C

. Ordinary differential equations. 3rd ed.New York: John Wiley and Sons, 1978.

Arnold

. Ordinary differential equations. Cambridge MA: The MIT Press, 1973.

Foss

. Co-ordinates which uncouple the equations of motion of damped linear dynamic systems. J Appl Mech 1958; 25: 361–364.

Caughey

. Classical normal modes in damped linear structures. J Appl Mech 1960; 27: 269–271.

Caughey

O’Kelly

MEJ

. Classical normal modes in damped linear dynamic systems. J Appl Mech 1965; 32: 583–588.

10.

Lancaster

. Lambda-matrices and vibrating systems. Oxford: Pergamon Press Inc., 1966.

11.

Lidström

. On the relative rotation of rigid parts and the visco-elastic torsion bushing element. Math Mech Solid 2013; 18: 788–802.

12.

Lidström

. On the equations of motion in constrained multibody dynamics. Math Mech Solid 2012; 17: 209–242.

13.

Clough

Penzien

. Dynamics of structures. Düsseldorf: McGraw-Hill, 1975.

14.

Greif

. The eigenproblem for damped linear vibrating systems with multiple eigenvalues. J Sound Vib 1989; 131: 197–214.

15.

Gohberg

Lancaster

Rodman

. Matrix polynomials. New York, London: Academic Press, 1982.

16.

Lancaster

Webber

. Jordan chains for lambda matrices. Lin Algebra Appl 1968; 1: 563–569.

17.

Lancaster

. Inversion of lambda-matrices and application to the theory of linear vibrations. Arch Ration Mech Anal 1960; 6: 105–114.

18.

Krechetnikov

Marsden

. Dissipation-induced instabilities in finite dimensions. Rev Mod Phys 2007; 79: 519–553.

19.

Duffin

. A minimax theory for overdamped networks. J Ration Mech Anal 1955; 4: 221–233.

20.

Rogers

. A minimax theory for overdamped systems. Arch Ration Mech Anal 1964; 16: 89–96.

21.

Inman

. A sufficient condition for the stability of conservative gyroscopic systems. J Appl Mech 1988; 55: 895–898.

22.

Brackwell

Lancaster

. Overdamped and gyroscopic vibrating systems. J Appl Mech 1992; 59: 176–181.

23.

Lancaster

Markus

Zhou

. A wider class of stable gyroscopic systems. Lin Algebra Appl 2003; 370: 257–257.