# Normal Modes of Vibration for Non-Linear Continuous Systems

Steven Shaw, Christophe Pierre

<figure><img src="image_1.png" /></figure>

## To cite this version:

Steven Shaw, Christophe Pierre. Normal Modes of Vibration for Non-Linear Continuous Systems. Journal of Sound and Vibration, 1994, 169 (3), pp.319-347. 10.1006/jsvi.1994.1021. hal-01471098

HAL Id: hal-01471098

https://hal.science/hal-01471098

Submitted on 28 Feb 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
NORMAL MODES OF VIBRATION FOR
NON-LINEAR CONTINUOUS SYSTEMS

S. W. SHAW AND C. PIERRE

Department of Mechanical Engineering and Applied Mechanics,
The University of Michigan, Ann Arbor, Michigan 48109-2125, U.S.A.

A definition and a constructive methodology for normal nodes of motion are developed for a class of vibratory systems the dynamics of which are governed by non-linear partial differential equations. The definition for normal modes is given in terms of the dynamics on two-dimensional invariant manifolds in the system phase space. These manifolds are a continuation of the planes which represent the well-known normal modes of the linearized system. A local asymptotic approximation of the geometric structure of these manifolds can be obtained by using an approach which follows that used for generating center manifolds. The procedure also provides the non-linear ordinary differential equations which govern individual modal dynamics and a physical description of the system configuration when it is undergoing a modal motion. In this paper, the general theory is described for the application of vibrations of continuous media, but it can easily be extended to other situations. In order to demonstrate the power of the approach and its unique procedural aspects, three examples involving beam vibrations are worked out in detail. The examples are conservative, simply supported beams. The first demonstrates the methodology for a linear beam model, the second is a beam on a non-linear elastic foundation, and the third example is a beam with non-linear torsional springs attached at each end. In these examples, the simplicity of the mode shapes of the linearized model yields relatively simple calculations that do not obscure the important features of the procedure.

# 1. INTRODUCTION

The concept of normal modes of motion is well developed for a wide class of dynamical systems, typically defined by linear ordinary or partial differential with constant coefficients, which are obtained by linearizing equations of motion about an equilibrium configuration (see, for example, Weaver *et al.* [1]). These normal modes represent a special set of motions of the system in which it behaves like a system of lower order (specifically, of second order in most vibrations applications). For linear, conservative, non-gyroscopic systems that are spatially distributed, each normal mode has associated with it a mode shape, given by the corresponding eigenfunction, and a natural frequency, determined by the corresponding eigenvalue. When damping and/or gyroscopic terms are present, normal modes also exist, but the situation is not so simple, as normal modes are typically travelling wave motions. However, for virtually all linear systems, the dynamics of each individual mode is governed by a second order linear modal oscillator which is uncoupled from all other modal oscillators.

The present work demonstrates the underlying geometrical structure which exists in the phase space for typical non-linear distributed parameter systems near an equilibrium (namely, the structure which is the source of the special normal mode motions). This geometrical structure is comprised of a set of invariant, normal mode manifolds. In this
paper, it is shown how these manifolds provide a means of constructing mode shapes and
modal oscillators for motions of such non-linear systems about equilibrium.

It is important to note that a complete analogy with linear modal analysis is not possible
for non-linear systems, simply due to the fundamental fact that superposition does not
hold. However, the existence of normal modes of motion for non-linear systems gives
significant insight into the dynamics of these systems and provides a systematic method
for carrying out model reduction.

The method presented herein is an extension of that proposed and developed recently
by the authors for finite-dimensional dynamical systems [2]. The motivating ideas are
from the theory of invariant manifolds, an approach which, at least in the engineering
community, has been completely overlooked as a means of defining normal modes. The
method employed is quite simple in concept: it determines a set of two-dimensional
invariant manifolds in the phase space which represent normal mode motions for the
non-linear system, and then provides the equations of motion which dictate the dynamics
on these manifolds. Utilizing the existence of the linear normal modes, one can locally
construct formal asymptotic series expansions for these normal mode manifolds and the
attendant modal oscillators. These series represent a form of non-linear separation of
variables which exactly recovers the linear dynamics when non-linearities are neglected and
systematically produces the non-linear corrections to the normal modes. Such an approach
is very adaptable, and may find use in many areas where one is interested in the non-linear
dynamics of spatially distributed systems.

The oscillators which describe the modal dynamics are completely uncoupled on the
individual modal manifolds (as they must be from invariance), except in cases where
internal resonance occurs. In such situations, there exists a coupling between the resonant
modes which cannot be eliminated. The present method signals such resonances in the
usual manner—by the appearance of small divisors.

Virtually all previous works on non-linear normal modes, for example those of
Rosenberg [3], Rand [4], and Vakakis [5], have been restricted to finite-dimensional,
conservative systems. Efforts to generate normal modes for non-linear continuous systems
include the works of Bennouna et al. [6, 7] and Szemplinska [8], which are also restricted
to conservative systems and, in addition, assume a separation of temporal and spatial
behavior as the first step in the analysis. Such methods are essentially a combination of
harmonic balance and eigenfunction expansions. The approaches found in references
[3–5] are more similar in spirit to the approach taken here, however, in that the form
of the system’s time behavior in a given mode is not prescribed *a priori*, but is
determined by differential equations of motion which are derived for these modal motions.
In reference [2], the authors introduced such a methodology, based on the concept of
invariant manifolds, which generalized the normal mode results for finite-dimensional
systems to include cases in which damping or gyroscopic terms are present. The extension
presented herein further generalizes the formulation to the important class of systems
modelled by partial differential equations. It differs from previous methods for continuous
systems, e.g., those found in references [6–8], in that it may be applied to non-conservative
and/or gyroscopic systems, and, more importantly, in that the behavior of the modal
dynamics is not specified *a priori*, but is determined from differential equations
which govern the modal dynamics. Furthermore, the approach presented herein does
not require a series representation of the mode shapes at the outset of the solution
procedure, as do the methods found in references [6–8]. Rather, the present method
generates, in a systematic manner, the differential equations which the mode shapes
must satisfy. It turns out that these equations often (but not always) are conveniently
solved by a series solutions, but the assumption is in no way required at the outset of the
procedure. Also, the present method is founded on basic principles of dynamical systems; namely, the representation of normal modes by invariant manifolds in the phase space, which provides a general theoretical formulation that is not restricted by the amplitude of motion. Also, while the general case is not solvable in closed form except in special situations, it naturally admits a constructive solution procedure for weakly non-linear systems. In addition to the above advantages, the present approach offers a new way to attack linear problems and provides a very clean alternative to the use of complex normal modes for systems with damping or gyroscopic effects. This will be described in detail in a forthcoming paper [9].

The paper is arranged as follows. Section 2 provides a description of the general approach for the class of problems considered. The examples are presented in section 3, and the paper closes with a discussion in section 4.

## 2. THE GENERAL METHOD

The motivation for the present work lies in applications from structural vibrations, and thus equations of motion which are of second order in time derivatives and of arbitrary order in spatial derivatives are considered. Furthermore, the class of systems is restricted to those which, when linearized, have simple normal modes which arise from a discrete spectrum, i.e., no systems with continuous spectral components are considered. The general method is developed in detail for one-dimensional motions of one-dimensional continua, e.g., planar, transverse vibrations of beams. The means of generalizing the method to other cases is described at the end of this section.

Consider the vibrations of a spatially extended system which occupies a one-dimensional region $\Omega$. Material points in the body in its equilibrium state are labelled by a position scalar $s$. When undergoing a general motion, the point $s$ is displaced such that the displacement for that point is $u(s, t)$, which is a scalar function of $s$ and time, $t$. The velocity of the material point at $s$ is given by $v(s, t) = \partial u(s, t)/\partial t$. It is assumed that $(u(s, t), v(s, t)) = (0, 0)$ is the equilibrium configuration of the system, this is easily achieved by a suitable choice of $u$. Under quite general circumstances, the equations of motion for excursions about equilibrium can be written in first order form as

$$ \begin{gathered} \partial u(s, t)/\partial t = v(s, t), \\ \partial v(s, t)/\partial t = F(u(s, t), v(s, t)), \quad s \in \Omega - \partial\Omega, \end{gathered} \quad (1) $$

with boundary conditions of

$$ B(u(s, t), v(s, t)) = 0, \quad s \in \partial\Omega, \quad (2) $$

where $F$ denotes some non-linear operator, typically of integro-differential type, in the spatial variable $s$, $\partial\Omega$ denotes the boundary of the region $\Omega$, and $B$ denotes the spatial boundary condition operators. Often one must invert linear or non-linear inertia operators in order to achieve this form (see, for example, Hsieh et al. [10]).

The definition of normal modes based on invariant manifolds is as follows: a normal mode of motion for a non-linear, autonomous system is a motion which takes place on a two-dimensional invariant manifold in the system's phase space. This manifold has the following properties: it passes through the stable equilibrium point $(u, v) = (0, 0)$ of the system and at $(u, v) = (0, 0)$ it is tangent to a plane which is an eigenspace of the system linearized about $(u, v) = (0, 0)$.

The key observation in the development of these normal modes is the following: in a normal mode motion, if the displacement $u_0(t) = u(s_0, t)$ and velocity $v_0(t) = v(s_0, t)$ of a
single point, $s = s_0$, are known, then the entire displacement and velocity fields can be determined by the dynamics of that single point. This can be stated mathematically, at least in some neighborhood of equilibrium, as follows:

$$u(s, t) = U(u_0(t), v_0(t), s, s_0), \quad v(s, t) = V(u_0(t), u_0(t), s, s_0), \quad (3)$$

where $U$ and $V$ represent functional relationships which relate the entire $(u, v)$ field to $(u_0, v_0)$ and which satisfy the boundary conditions. These relations are essentially a type of constraint, such that a motion which satisfies them is completely determined by the state $(u_0, v_0)$. The concomitant dynamics of the system are governed by a second order non-linear ordinary differential equation, that is, it behaves like a simple non-linear oscillator. A more geometrical way of viewing equations (3) is that they represent a two-dimensional manifold in the infinite-dimensional phase space of the system. There will be one such manifold for each normal mode.

The following obvious identities on $U$ and $V$, evaluated at $s = s_0$, are noteworthy:

$$U(u_0(t), v_0(t), s_0, s_0) = u_0(t), \quad V(u_0(t), v_0(t), s_0, s_0) = v_0(t). \quad (4)$$

Conditions dictate that these manifolds be invariant for the equations of motion are now derived. The approach followed here is similar to that used in the construction of center manifolds (see Carr [11]). The aim is to eliminate all time derivatives from the equations of motion by using the constraint conditions, after which one is left with equations which describe the geometry of the invariant manifolds; these are given by $(U, V)$ in terms of $u_0$ and $v_0$. Once the manifolds are obtained, the dynamics on them are obtained directly by evaluating the equations of motion on the normal mode manifold.

The process begins by taking a time derivative of equation (3), which yields

$$\frac{dU}{dt} = \frac{\partial U}{\partial u_0} u_{0,t} + \frac{\partial U}{\partial v_0} v_{0,t}, \quad \frac{dV}{dt} = \frac{\partial V}{\partial u_0} u_{0,t} + \frac{\partial V}{\partial v_0} v_{0,t}. \quad (5)$$

(For notational simplicity, the explicit dependence of $u, v, u_0, v_0, U$ and $V$ on $s, s_0$, and $t$ is dropped wherever no confusion is possible. Also, note that time derivatives are denoted here by $(\cdot)_t \equiv \partial(\cdot)/\partial t$ or by $(\cdot)_v \equiv \partial(\cdot)/\partial t$, as convenient.) The equations of motion are then used to replace $u$, and $v, v_t$ with $F(u, v)$, $u_{0,t}$ by $v_0$, and $v_{0,t}$ by $F(u, v)|_{s=s_0}$ (that is, $F(u, v)$ evaluated at $s = s_0$). Then the constraint is enforced everywhere by replacing $u$ by $U$ and $v$ by $V$, including in the boundary conditions. This results in the following semi-linear, hyperbolic partial differential equations, which are to be solved for $U$ and $V$ as functions of $u_0, v_0, s$ and $s_0$:

$$V = \frac{\partial U}{\partial u_0} v_0 + \frac{\partial U}{\partial v_0} [F(U, V)]_{s=s_0},$$

$$F(U, V) = \frac{\partial V}{\partial u_0} v_0 + \frac{\partial V}{\partial v_0} [F(U, V)]_{s=s_0}, \quad s \in \Omega - \partial\Omega, \quad (6)$$

with boundary conditions

$$B(U, V) = 0, \quad s \in \partial\Omega. \quad (7)$$

As is typical in such procedures, these equations are at least as difficult to solve as the original differential equations, but they do allow for solutions to be obtained in the form of power series, here in terms of $u_0$ and $v_0$. At this point, it should be noted that the solution procedure given below for these equations restricts one to considering local non-linear effects up to some order in an asymptotic sense, but that equations (6) and (7) are valid in a more global sense. In fact, if other means are used for solving them, for example, if special symmetries allow for global solutions or if numerical solutions can be obtained,
then the range of validity of the modal manifolds and the modal dynamics can be extended. One obvious limitation of this formulation is that the functional form given in equation (3) breaks down at any point where a manifold bends back on itself, such that it becomes multi-valued in terms of $u_0$ and/or $v_0$. However, this is of no real concern in applications in which first order non-linear effects are of interest.

The dynamics of the system, restricted to a normal mode manifold, is then easily captured by substituting into the equations of motion the corresponding solution for $U$ and $V$, performing all spatial operations, and then evaluating everything at $s = s_0$. This results in the following two-dimensional dynamical system for $u_0$ and $v_0$:

$$du_0/dt = v_0, \quad dv_0/dt = F(u(s, t), v(s, t))|_{s=s_0}, \quad (8)$$

or, in second order form,

$$\frac{d^2 u_0}{dt^2} - F\left(u(s, t), \frac{\partial u(s, t)}{\partial t}\right)|_{s=s_0} = 0, \quad (9)$$

which is the modal oscillator associated with the manifold defined by $U$ and $V$. It describes the dynamics of the point $s_0$ when the system is undergoing a purely single mode motion, while equation (3) provides the physical manner in which the system is configured during the modal motion. There is one such modal oscillator for each non-linear normal mode, and each normal mode for the linearized system has a corresponding non-linear normal mode (except in cases of internal resonance).

The form of equation (8), which trivially converts to a second order differential equation in $u_0$, is denoted as the *oscillator form* of the modal dynamics and emerges directly from the manner in which the problem is formulated. Specifically, it is due to the use of $u_0$ and $v_0$ as the two independent variables which span the two-dimensional invariant manifolds. This may not be the most convenient representation for all applications, but it has many appealing features for mechanics problems.

Local solutions for these manifolds can be obtained by a procedure which is similar in spirit to that used in constructing center manifolds (Carr [11]). One proceeds by assuming an approximate, asymptotic series solution for equation (6) and the boundary conditions (7), of the form

$$\begin{align}
U(u_0, v_0, s, s_0) &= a_1(s, s_0)u_0(t) + a_2(s, s_0)v_0(t) \nonumber \\
&\quad + a_3(s, s_0)u_0(t)^2 + a_4(s, s_0)u_0(t)v_0(t) \nonumber \\
&\quad + a_5(s, s_0)v_0(t)^2 + a_6(s, s_0)u_0(t)^3 + \dots, \nonumber \\
V(u_0, v_0, s, s_0) &= b_1(s, s_0)u_0(t) + b_2(s, s_0)v_0(t) \nonumber \\
&\quad + b_3(s, s_0)u_0(t)^2 + b_4(s, s_0)u_0(t)v_0(t) \nonumber \\
&\quad + b_5(s, s_0)v_0(t)^2 + b_6(s, s_0)u_0(t)^3 + \dots, \tag{10}
\end{align}$$

which is a type of non-linear separation of variables. The $a_i$'s and $b_i$'s contain information about the spatial distribution for a given mode, and $u_0$ and $v_0$ represent the attendant modal displacement and velocity, i.e., the time behavior. Note that the mode shape depends on the amplitudes of $u_0$ and $v_0$, since the relative contributions of the various $a_i$'s and $b_i$'s vary depending on the magnitudes of $u_0$ and $v_0$. In fact, the mode shape depends not only on the peak amplitude of motion, but the spatial configuration of the system changes as a function of time *during a given motion*. This is in contrast to the methods which initially assume spatial and temporal separation [6-8].

The $a_i$'s and $b_i$'s are solved for by substituting the series into equations (6) and boundary conditions (7), expanding in terms of $u_0$ and $v_0$, gathering terms of like powers in $u_0$ and
$v_0$, and matching coefficients of $u_0^m v_0^n$ with $m, n = 0, 1, 2, \dots$ and $m + n \ge 1$. This results in a sequence of boundary value problems for the $a_i$'s and $b_i$'s in the $s$ variable. These problems are uncoupled in sequential order (as is typical in such series solutions), allowing for a direct solution strategy. It can be shown that the solutions obtained in this way for the coefficients of the linear terms in $u_0$ and $v_0$; that is, $a_1, a_2, b_1$ and $b_2$ recover, in a rather unusual manner, the linear eigensolution with a peculiar normalization. (See section 3.1 for an example; the general proof is to be published [9].) This implies that the manifolds represented by $U$ and $V$ are two-dimensional invariant manifolds which are tangent to the invariant planes which represent the eigenspaces for the linearized problem.

It is also worth noting that the identities given in equation (4) above and the independence of $u_0$ and $v_0$ require the following conditions to hold for the $a_i$'s and the $b_i$'s:

$$ \begin{aligned} a_1(s_0, s_0) &= 1, & a_2(s_0, s_0) &= 0, & a_j(s_0, s_0) &= 0, & j &= 3, 4, 5, \dots, \\ b_1(s_0, s_0) &= 0, & b_2(s_0, s_0) &= 1, & b_j(s_0, s_0) &= 0, & j &= 3, 4, 5, \dots \end{aligned} \quad (11) $$

By considering $s$ evaluated at another point, say $s = \bar{s}_0$, other identities which the $a_i$'s and $b_i$'s must satisfy can be obtained. Some of these are derived and presented in Appendix A. Since $s_0$ is a free parameter, and not an independent variable, the explicit dependence of the $a_i$'s and the $b_i$'s on $s_0$ is dropped from this point on. This freedom for selecting $s_0$ will have some interesting consequences in the example problems presented below.

This section is closed with some remarks. First, although the above solution strategy is based on series expansions about equilibrium, which infers that the linear results are simply being extended, the normal mode manifolds are more globally defined and should be thought of as the *source* of the linear eigenmodes, and not as mere extensions of them. Also, the special features of linear eigenmodes—that they are generally uncoupled and can be directly superimposed—arise from the fact that all non-linear coupling terms in the equations of motion expressed in terms of non-linear modal co-ordinates vanish *by definition* when the system is linearized.

Second, it may be possible to represent the normal mode manifolds with co-ordinates other than $(u_0, v_0)$. For example, for pointwise measurement, control or actuation, it might be more suitable to use the displacements (or velocities) of two points on the structure for each mode. The formulation will not be as convenient, but is feasible in principle.

Finally, although the above formulation has been derived for one-dimensional vibrations of one-dimensional continua, the basic ideas can be generalized to other situations, including $m$-dimensional vibrations of $n$-dimensional continua. In that case, $\Omega$ is an $n$-dimensional region and $s$ is a $n$-dimensional vector, while $u, v, U, V$ and $F$ will be $m$-dimensional (see Figure 1). For example, for transverse vibrations of a plate, $n = 2$ and $m = 1$, and for non-planar vibrations of a slender beam, $n = 1$ and $m = 2$. In the

<figure><img src="image_3.png" /><figcaption>Figure 1. Schematic diagram for the general theory.</figcaption></figure>
general case, the formulation differs from the one given above in that $u_0$ cannot simply be $u(s, t)$ evaluated at $s_0$, but must be taken to be a particular element of $u(s, t)$ evaluated at $s_0$. This is so because $u_0$ must be a scalar in order for the $(u_0, v_0)$ dynamics to be two-dimensional; that is, those of a normal mode. The generation of equations (6) and the construction of the modal dynamics from equation (8) require the use of the element of $F$ which corresponds to the direction chosen for $(u_0, v_0)$. Care must be taken in selecting this element, since modes which have no motion in the chosen direction will be completely overlooked. In this general setting, the $a_i$'s and $b_i$'s are $m$-dimensional vectors, and these provide the required distribution of $u_0$ amongst all the elements of the vectors $u$ and $v$. Also, note that in cases for which $n > 1$, the identities given in equations (4) and (11) are not valid and must be modified.

### 3. EXAMPLES

The examples are selected to illustrate the method and the types of results which can be obtained. The example systems are chosen from transverse vibrations of beams in one spatial dimension; this simplifies the calculations and allows for closed form solutions to be obtained in many cases. In each example, the equation of motion considered has been non-dimensionalized, so that the minimum number of parameters appears. Also, the standard notation is employed in which subscripts $s$ and/or $t$ denote partial derivatives with respect to $s$ and/or $t$, respectively. Whenever it is clear, overdots are used for derivatives with respect to $t$ and primes are used for derivatives with respect to $s$.

Three examples are presented, each of which is a variation on the central theme of transverse vibrations of simply supported Euler-Bernoulli beams. The case of a cantilever beam with geometric curvature and inertial non-linearities is considered in reference [10]. The first example considers the simplest case—transverse vibrations of a linear, simply supported beam. Although this example may seem to be too trivial to occupy journal pages, the features of the present approach are quite unusual and require demonstration. Furthermore, this example provides results which are used in the following examples. The second system considered is a simply supported beam which lies on a non-linear elastic foundation. The third example is a linear simply supported beam with non-linear torsional springs at each end. This demonstrates how the method can be used to handle discrete non-linear elements, such as those which often appear in boundary conditions.

The calculations were carried out with the aid of the computer assisted symbolic manipulation program Mathematica, running on a NeXT computer.

#### 3.1. EXAMPLE 1: TRANSVERSE VIBRATIONS OF A LINEAR, SIMPLY SUPPORTED BEAM

This example considers, from a point of view different from that found in any textbook, the classical problem of transverse vibrations of an Euler-Bernoulli beam. The first order form of the equation of motion is given by

$$u_t = v, \quad v_t = -u_{ssss}, \quad s \in (0, 1), \qquad (12)$$

where $u = u(s, t)$ is the transverse displacement of the beam, $v = v(s, t)$ is the transverse velocity of the beam, $s$ is the independent spatial variable, and $t$ is time. The beam is pinned at the left end and simply supported at the right end. This allows for non-zero axial displacement at the right end and avoids membrane strain effects. The associated boundary conditions are

$$u(0, t) = u(1, t) = 0, \quad u_{ss}(0, t) = u_{ss}(1, t) = 0. \qquad (13)$$
A search for normal mode motions begins by assuming that there exists at least one motion for which the entire beam behaves like a second order oscillator. The implementation of this assumption is carried out by writing the displacement and velocity fields for the beam as being dependent on the displacement and velocity of a single point on the beam, the point $s = s_0$. Writing $u(s, t) = u_0(t)$ and $v(s, t) = v_0(t)$ (or simply $u_0$ and $v_0$ when convenient), this dependence is expressed as given in equation (3). In this linear case, the relationship relating $U$ and $V$ to $u_0$ and $v_0$ is linear, and is given by the leading order terms of the series expansions, equation (10), as

$$ 
\begin{aligned}
u(s, t) &= U(u_0, v_0, s, s_0) = a_1(s, s_0)u_0(t) + a_2(s, s_0)v_0(t), \\
v(s, t) &= V(u_0, v_0, s, s_0) = b_1(s, s_0)u_0(t) + b_2(s, s_0)v_0(t).
\end{aligned} 
\quad (14) $$

The procedure is described in section 2 is now continued by taking a time derivative of these two equations to obtain

$$ 
\begin{aligned}
u_t &= \frac{\partial U}{\partial u_0} u_{0,t} + \frac{\partial U}{\partial v_0} v_{0,t} = a_1 u_{0,t} + a_2 v_{0,t}, \\
v_t &= \frac{\partial U}{\partial v_0} u_{0,t} + \frac{\partial U}{\partial u_0} v_{0,t} = b_1 u_{0,t} + b_2 v_{0,t}.
\end{aligned} 
\quad (15) $$

In order to remove all time derivatives, the following substitutions are now made by utilizing the equation of motion: $v(s, t)$ for $u_t(s, t)$, $-u_{ssss}(s, t)$ for $v_t(s, t)$, $v_0(t)$ for $u_{0,t}(t)$ and $-u_{ssss}(s_0, t)$ for $v_{0,t}(t)$. This yields the two equations

$$ 
\begin{aligned}
v(s, t) &= a_1(s)v_0(t) - a_2(s)u_{ssss}(s_0, t), \\
-u_{ssss}(s, t) &= b_1(s)v_0(t) - b_2(s)u_{ssss}(s_0, t).
\end{aligned} 
\quad (16) $$

The final step in setting up the equations which are to be solved for the mode shapes is to substitute everywhere into the above for $u$, $v$, and their derivatives with respect to $s$, the expressions given in equation (14) for $U$ and $V$. Note that some of these expressions are evaluated at $s = s_0$ after derivatives in $s$ have been carried out. This results in the equations

$$ 
\begin{aligned}
b_1(s)u_0(t) + b_2(s)v_0(t) &= a_1(s)v_0(t) - a_2(s)(a_1'''(s_0)u_0(t) + a_2'''(s_0)v_0(t)), \\
-(a_1'''(s)u_0(t) + a_2'''(s)v_0(t)) &= b_1(s)v_0(t) - b_2(s)(a_1'''(s_0)u_0(t) + a_2'''(s_0)v_0(t)),
\end{aligned} 
\quad (17) $$

where primes denote derivatives with respect to $s$. These are the specific versions of equation (6) for this example. At this stage, the independence of $u_0$ and $v_0$ is brought to bear by requiring that these equations be valid for all values of $u_0$ and $v_0$. In particular, setting $u_0$ and $v_0$ equal to zero in turn yields four independent equations, exactly as needed in order to solve for $a_1$, $a_2$, $b_1$ and $b_2$. Equivalently, these are the equations obtained by simply equating coefficients of $u_0$ and $v_0$. The resulting four equations are

$$ 
\begin{aligned}
b_1(s) + a_2(s)a_1'''(s_0) &= 0, & a_1(s) - b_2(s) - a_2(s)a_2'''(s_0) &= 0, \\
a_1'''(s) - b_2(s)a_1'''(s_0) &= 0, & a_2'''(s) + b_1(s) - b_2(s)a_2'''(s_0) &= 0.
\end{aligned} 
\quad (18) $$

These four coupled, algebraic-differential, homogeneous equations have a very unusual structure. Note that the identities $a_1(s_0) = 1$, $a_2(s_0) = 0$, $b_1(s_0) = 0$ and $b_2(s_0) = 1$ are satisfied by equation (18). At this point it should also be noted that the $b_i$'s appear in a linear, non-differential manner; this is due to the fact that no spatial derivatives of $v(s, t)$ appear in the equations of motion. This allows for expression of the $b_i$'s in terms of the $a_i$'s in the first two equations, and thus for a simple elimination of the $b_i$'s from the last
two equations, yielding two coupled, fourth order, homogeneous, differential equations in the $a_i$'s.

Before discussing the solution of these equations, the boundary conditions need to be considered. First, note that boundary conditions on the velocity are required in this formulation; these are obtained by simply taking derivatives with respect to $t$ of the boundary conditions on $u$ and using the definition $v = u_i$. A direct substitution of the forms for $U$ and $V$ in terms of the $a_i$'s and $b_i$'s, and the use of the independence of $u_0$ and $v_0$ yield the following boundary conditions on the $a_i$'s and $b_i$'s:

$$ 
\begin{aligned}
a_1(0) &= a_1(1) = a_2(0) = a_2(1) = 0, & b_1(0) &= b_1(1) = b_2(0) = b_2(1) = 0, \\
a_1''(0) &= a_1''(1) = a_2''(0) = a_2''(1) = 0, & b_1''(0) &= b_1''(1) = b_2''(0) = b_2''(1) = 0.
\end{aligned} 
\quad (19) $$

Note that there are 16 boundary conditions: the eight conditions on the $a_i$'s are for the two fourth order differential equations which the $a_i$'s must satisfy, while the eight conditions for the $b_i$'s must also be satisfied, as the $b_i$'s are linearly related to the $a_i$'s.

In equations (18), the terms which depend on $s$ appear in a linear manner, but some have coefficients which depend on the solution evaluated at $s_0$. This results in a non-linear structure, but one which is amenable to solution strategies used for linear differential equations. Note that if one attempts to solve these equations by considering those terms which involve the solution evaluated at $s = s_0$ as constants, the systems eigenvalues will depend on the solution! This bothersome structure is not as difficult as it initially appears—the solution of a homogeneous problem typically has a free constant, and in this case this constant is not free but is fixed in such a manner that the equations are satisfied. This is a direct consequence of the formulation, and is in fact required if the identities in equation (11) are to be satisfied. This is clearly demonstrated in what follows.

The following solution procedure is used for the present. It is first observed that for a conservative, non-gyroscopic system such as the one considered here, the displacement and velocity are not coupled (at linear order) and therefore $a_2$ and $b_1$ are identically zero; this automatically satisfies the first and fourth equations in equation (18). The remaining two equations in equation (18) simplify directly to the identity

$$ a_1(a) = b_2(s), \quad (20) $$

and the differential equation

$$ a_1'''(a) = a_1(s)a_1'''(s_0) = 0, \quad (21) $$

a much simpler situation. Note that by evaluating this equation at $s = s_0$, the condition $a_1(s_0) = 1$ is returned.

At this stage, the $a_1'''(s_0)$ term is treated as a constant and the usual solution is assumed in terms of a linear combination of $\sin(\lambda s)$, $\cos(\lambda s)$, $\sinh(\lambda s)$ and $\cosh(\lambda s)$. The boundary conditions on $a_1$ nullify all terms except the $\sin(\lambda s)$ term, and furthermore require the interesting identity

$$ \lambda^4 = a_1'''(s_0) = (n\pi)^4, \quad n = 1, 2, 3, \dots, \quad (22) $$

for a non-trivial solution to exist. The solution is thus of the form

$$ a_1(s) = A \sin(\lambda s) = A \sin(n\pi s), \quad n = 1, 2, 3, \dots \quad (23) $$

This solution must satisfy the condition $a_1(s_0) = 1$. This condition, or equation (22), directly gives the following condition on the constant $A$:

$$ A = 1/\sin(n\pi s_0), \quad n = 1, 2, 3, \dots \quad (24) $$
Hence, the solution for $a_1$ and $b_2$ is

$$a_1(s) = b_2(s) = \frac{\sin(n\pi s)}{\sin(n\pi s_0)}, \quad n = 1, 2, 3, \dots \quad (25)$$

which has the desired (and expected) spatial dependence and which satisfies the required identities; namely, that $U$ and $V$ evaluated at $s = s_0$ yield $u_0$ and $v_0$, respectively. Note that there are a countable infinity of modes and that these agree with the well-known linear eigensolution. This is most easily seen by considering a motion in the first mode with a peak mid-span deflection of $A$, in which case the mode shape is given by $A \sin(n\pi s)$. For the present formulation, such a motion has, for $s_0 \in (0, 1)$, $u_0 = A \sin(n\pi s_0)$, verifying the equivalence. Also note that the restriction that $s_0$ should not be a nodal point for the $n$th linear mode appears here. This is reasonable, since the approach will obviously break down if $s_0$ is chosen at a node. This restriction causes no difficulty, since one is free to choose a different $s_0$ for each mode; and a choice of, say $s_0 = 1/2n$ yields nice results for the linear modes. Finally, it should be remarked that since there exists a countable infinity of solutions for $a_1$, it might be preferable to denote each modal solution by a subscript $n$. For notational simplicity, this is done for $U$ and $V$ but not for the $a_i$'s or the $b_i$'s.

The normal modes are thus given by the functions

$$U_n(u_0(t), v_0(t), s, s_0) = \frac{\sin(n\pi s)}{\sin(n\pi s_0)} u_0(t), \quad V_n(u_0(t), v_0(t), s, s_0) = \frac{\sin(n\pi s)}{\sin(n\pi s_0)} v_0(t), \quad (26)$$

which is valid as long as $s_0$ is not a node. The equations which govern the dynamics for a normal mode are determined by substituting equations (26) into the equations of motion (12) and evaluating these at $s = s_0$. This leads to ordinary differential equations which govern the motion of the point $s_0$, from which the motion of the entire beam is dictated by $U_n$ and $V_n$. The general form for this case is

$$\frac{\partial U_n}{\partial t}(u_0(t), v_0(t), s_0, s_0) = V_n(u_0(t), v_0(t), s_0, s_0),$$

$$\frac{\partial V_n}{\partial t}(u_0(t), v_0(t), s_0, s_0) = -\frac{\partial^4 U_n}{\partial s^4}(u_0(t), v_0(t), s_0, s_0). \quad (27)$$

Using the identities at $s = s_0$ this simplifies to

$$\dot{u}_0 = v_0, \quad \dot{v}_0 = -(n\pi)^4 u_0, \quad (28)$$

where a dot represents a time derivative. In second order form, these equations simply become

$$\ddot{u}_0 + (n\pi)^4 u_0 = 0, \quad (29)$$

from which the modal natural frequencies of $\omega_n = (n\pi)^2$ are directly obtained.

It is important to note that in this formulation one first obtains the mode shapes, after which the modal dynamics, including the natural frequencies, are determined from the differential equations which govern the modal dynamics. This property is a consequence of the method. In this approach, one does not need to know the form of the time behavior (e.g., exponential or sinusoidal in the linear case) of the system in a modal motion in order to construct its normal modes. The feature allows the approach to be extended to non-linear systems.
## 3.2. EXAMPLE 2: TRANSVERSE VIBRATIONS OF A SIMPLY SUPPORTED BEAM ON A NON-LINEAR ELASTIC FOUNDATION

This example considers the vibrations of a simply supported beam which is attached to an elastic foundation. The foundation has a non-dimensional linear stiffness coefficient of $k$ and a non-dimensional cubic stiffness coefficient of $\gamma$. The physical system is depicted in Figure 2. The equations of motion for this case are

$$u_t = v, \quad v_t = -u_{ssss} - ku - \gamma u^3, \quad s \in (0, 1), \qquad (30)$$

and the boundary conditions are the same as those in the first example (equations (13)). The normal mode shapes and their attendant dynamics for this non-linear problem will be generated out to third order. To that end, the series expansions for $U$ and $V$ (equation (10)) are taken out to third order in $u_0$ and $v_0$, thus requiring that solutions be obtained for all $a_j$ and $b_j$ for $j = 1, 2, ..., 9$. The boundary conditions on the $a_i$'s and $b_i$'s are determined by direct substitution to be ($j = 1, 2, 3, ...$)

$$a_j(0) = a_j(1) = b_j(0) = b_j(1) = 0, \quad a_j''(0) = a_j''(1) = b_j''(0) = b_j''(1) = 0. \qquad (31)$$

The solution procedure can be carried out by extending the steps from the first example to include non-linear terms, or one can simply tackle the equivalent problem of obtaining a series solution to equation (6). At the linear order either approach leads to the following equations for $a_1, a_2, b_1$ and $b_2$:

$$\begin{align}
b_1(s) + a_2(s)(a_1'''(s_0) + ka_1(s_0)) &= 0, & a_1(s) - b_2(s) - a_2(s)(a_2'''(s_0) + ka_2(s_0)) &= 0, \nonumber \\
a_1'''(s) + ka_1(s) - b_2(s)(a_1'''(s_0) + ka_1(s_0)) &= 0, \nonumber \\
a_2'''(s) + ka_2(s) + b_1(s) - b_2(s)(a_2'''(s_0) + ka_2(s_0)) &= 0. \tag{32}
\end{align}$$

Again, note that $a_1(s_0) = 1$, $a_2(s_0) = 0$, $b_1(s_0) = 0$ and $b_2(s_0) = 1$ satisfy these equations. Also, these equations reduce to those from the first example when $k = 0$. The boundary conditions are given in equation (31). The solution again has $a_2$ and $b_1$ identically zero, along with the identity

$$a_1(s) = b_2(s).$$

After some cancellation, this yields the differential equation

$$a_1'''(s) - a_1(s)a_1'''(s_0) = 0. \qquad (33)$$

Again, as in Example 1, the $a_1'''(s_0)$ term is treated as a constant and the solution is found to contain only a $\sin(\lambda s)$ term. For non-trivial solutions to exist, the following identity is also required:

$$\lambda^4 = a_1'''(s_0) = (n\pi)^4, \quad n = 1, 2, 3, ... \qquad (34)$$

<figure><img src="image_4.png" /><figcaption>Figure 2. The physical system for Example 2: a simply supported beam attached to an elastic foundation with linear stiffness coefficient $k$ and cubic stiffness coefficient $\gamma$.</figcaption></figure>
Following the first example, this again yields the linear mode shapes

$$a_1(s) = b_2(s) = \frac{\sin(n\pi s)}{\sin(n\pi s_0)}, \quad n = 1, 2, 3, \dots,$$

with the following modal oscillators at the linear mode

$$\ddot{u}_0 + ((n\pi)^4 + k)u_0 = 0, \quad n = 1, 2, 3, \dots \quad (35)$$

At this stage the known result is recovered, in which a linear elastic foundation shifts the natural frequencies while leaving the mode shapes unaffected (cf., reference [1], page 459). Note that for this non-linear problem, the requirement appears that $s_0$ should not be a node of the linear mode.

The ability of the method to handle non-linear problems is now demonstrated by proceeding to quadratic and cubic order terms.

At the quadratic order there are six equations, obtained by gathering coefficients of $u_0^2$, $u_0v_0$ and $v_0^2$ from each of the two equations (6). These are given below for the case in which simplifications have been made using the results from the linear order; namely, $a_2(s) = b_1(s) = 0$, $b_2(s) = a_1(s)$, $a_1(s_0) = 1$ and $a_1'''(s_0) = (n\pi)^4$.

From the first equation, we have the $u_0^2$ coefficient,

$$(k + (n\pi)^4)a_4(s) + b_3(s) = 0;$$

the $u_0v_0$ coefficient

$$-2a_3(s) + 2(k + (n\pi)^4)a_5(s) + b_4(s) = 0;$$

and the $v_0^2$ coefficient

$$-a_4(s) + b_5(s) = 0.$$

From the second equation, we have the $u_0^2$ coefficient, \hfill (36)

$$-ka_3(s) - a_3'''(s) + a_1(s)(ka_3(s_0) + a_3'''(s_0)) + (k + (n\pi)^4)b_4(s) = 0;$$

the $u_0v_0$ coefficient,

$$-ka_4(s) - a_4'''(s) + a_1(s)(ka_4(s_0) + a_4'''(s_0))$$

$$-2b_3(s) + 2(k + (n\pi)^4)b_5(s) = 0;$$

and the $v_0^2$ coefficient,

$$-ka_5(s) - a_5'''(s) + a_1(s)(ka_5(s_0) + a_5'''(s_0)) - b_4(s) = 0.$$

These equations have some similarities to those encountered at the linear level. In particular, the equations are of algebraic-differential type, are homogeneous and admit the trivial solution. Note that terms such as $a_j(s_0)$ have been left in to show the structure of the equations, but that due to identities (11) all such terms must be zero at all non-linear orders; that is, for all $j \ge 3$. It is also interesting to note that each individual term involves only one spatially dependent function. In some cases this is $a_1(s)$, in which case the coefficient is one or more of the unknowns or its derivatives evaluated at $s = s_0$. Finally, note that the system non-linearity, characterized by the cubic stiffness coefficient $\gamma$, does not appear in these equations.

The solution at this order must be approached in a manner different to that used at the linear order, since the quadratic terms must be zero in the absence of the non-linearity. This leads to the conclusion that the trivial solution is the required one at this order. This can be seen by considering the alternative: any non-trivial solution of equation (36) would not vanish in the case $\gamma = 0$, leaving arbitrary constants and
non-linear terms in the modes of the linear system. These terms must be zero in order to satisfy the identities given in equation (11). Therefore, the trivial solution is taken at the quadratic order:

$$a_j(s) = b_j(s) = 0, \quad j = 3, 4, 5. \tag{37}$$

In fact, since the problem contains only odd (linear and cubic) order terms, only odd order terms are to be expected in the modes. This is verified by the cubic order calculations which are considered next.

At the cubic order there are eight equations, obtained by gathering coefficients of $u_0^3$, $u_0^2 v_0$, $u_0 v_0^2$ and $v_0^3$ from each of the two equations (6). These are given below for the case in which simplifications have been made using the results from the linear and quadratic orders.

From the first equations, we have the $u_0^3$ coefficient,

$$(k + (n\pi)^4)a_7(s) + b_6(s) = 0;$$

the $u_0^2 v_0$ coefficient,

$$-3a_6(s) + 2(k + (n\pi)^4)a_8(s) + b_7(s) = 0;$$

the $u_0 v_0^2$ coefficient,

$$-2a_7(s) + 3(k + (n\pi)^4)a_9(s) + b_8(s) = 0;$$

and the $v_0^3$ coefficient,

$$-a_8(s) + b_9(s) = 0.$$

From the second equation, we have the $u_0^3$ coefficient,

$$\gamma(a_1(s) - a_1(s)^3) - ka_6(s) - a_6''''(s) + a_1(s)(ka_6(s_0) + a_6''''(s_0)) + (k + (n\pi)^4)b_7(s) = 0; \tag{38}$$

the $u_0^2 v_0$ coefficient,

$$-ka_7(s) - a_7''''(s) + a_1(s)(ka_7(s_0) + a_7''''(s_0)) - 3b_6(s) + 2(k + (n\pi)^4)b_8(s) = 0;$$

the $u_0 v_0^2$ coefficient,

$$-ka_8(s) - a_8''''(s) + a_1(s)(ka_8(s_0) + a_8''''(s_0)) - 2b_7(s) + 3(k + (n\pi)^4)b_9(s) = 0;$$

and the $v_0^3$ coefficient,

$$-ka_9(s) - a_9''''(s) + a_1(s)(ka_9(s_0) + a_9''''(s_0)) - b_8(s) = 0.$$

Note again that all terms of the form $a_j(s_0)$ for $j > 2$ are actually zero. These eight algebraic-differential equations are first reduced by solving the first four equations for the $b_i$'s in terms of the $a_i$'s (this step is trivial since only one $b_i$ appears in each equation). This result is then substituted into the remaining four differential equations, leading to a set of four equations which are to be solved for the $a_i$'s. These equations have a convenient structure in that they constitute two coupled pairs of equations, one in $a_6$ and $a_8$ and the other in $a_7$ and $a_9$. The two equations involving $a_7$ and $a_9$ are homogeneous, similar to the quadratic equation (36), and, using the same reasoning—the system non-linearity does not appear in them—the trivial solution must be taken from $a_7$ and $a_9$.
From the linear-algebraic equations, this immediately implies that $b_6$ and $b_8$ must be also zero. Thus,

$$a_7(s) = a_9(s) = b_6(s) = b_8(s) = 0. \qquad (39)$$

The fact that these terms are zero is due to the conservative, non-gyroscopic nature of the system. The remaining terms correspond to terms in the mode shapes and modal oscillators which represent standing wave normal modes, as expected for this system. This will become evident when the mode shapes and modal oscillators are constructed.

The remaining two equations are coupled, fourth order, non-homogeneous, differential equations in terms of $a_6$ and $a_8$. These are given by (with the $a_j(s_0)$ terms set to zero)

$$ \begin{aligned} & (2k + 3(n\pi)^4)a_6(s) - 2(k^2 + 2k(n\pi)^4 + (n\pi)^8)a_8(s) - a_6'''(s) + a_1(s)a_6''''(s_0) \\ & = \gamma(-a_1(s) + a_1(s)^3), \\ & - 6a_6(s) + 6ka_8(s) + 7(n\pi)^4a_8(s) - a_8'''(s) + a_1(s)a_8''''(s_0) = 0, \end{aligned} \qquad (40) $$

where the boundary conditions are given in equation (31). Note that the non-homogeneous term is proportional to the non-linear stiffness coefficient $\gamma$. It is the term which generates the distortions of the mode shapes which arise from the non-linearity. Once the solutions for $a_6$ and $a_8$ are determined, the other non-trivial terms, $b_7$ and $b_9$, are easily obtained from equation (38) above.

Only the particular solution of equations (40) is sought at this (or any) non-linear order, since the non-homogeneous terms arise directly from the non-linearity. Therefore, the trivial solution for $a_6$, $a_8$, $b_7$ and $b_9$ will, as required, be returned when the system is linearized, i.e., when $\gamma = 0$. Due to the nature of the boundary conditions, the particular solution can be written as a Fourier sine series. In fact, some simple trigonometric identities show that the non-homogeneous terms, due to the simple character of the non-linearity, are composed of only the first and third spatial harmonics, i.e., $\sin(n\pi s)$ and $\sin(3n\pi s)$ (recalling the solution for $a_1$). Thus, the desired solution will be a finite series of the form

$$a_6(s) = \beta_1 \sin(n\pi s) + \beta_3 \sin(3n\pi s), \quad a_8(s) = \eta_1 \sin(n\pi s) + \eta_3 \sin(3n\pi s). \qquad (41)$$

Substitution of these into the differential equations and projection onto $\sin(n\pi s)$ and $\sin(3n\pi s)$ results in four linear equations for the $\beta$'s and $\eta$'s. Using the vector of unknowns defined as

$$\mathbf{w} = (\beta_1, \beta_3, \eta_1, \eta_3)^T,$$

these linear equations can be written in matrix form as $A\mathbf{w} = \mathbf{b}$, where $A$ is the $4 \times 4$ matrix which consists of the four row vectors

$$ \begin{aligned} A_1 &= \left( \frac{3(k + (n\pi)^4)}{2}, \frac{3k}{2} - 2kS_0^2 + \frac{243(n\pi)^4}{2} - 162(n\pi)^4 S_0^2, -k^2 - 2k(n\pi)^4 - (n\pi)^8, 0 \right), \\ A_2 &= (0, k - 39(n\pi)^4, 0, -k^2 - 2k(n\pi)^4 - (n\pi)^8), \\ A_3 &= \left( -3, 0, \frac{7(k + (n\pi)^4)}{2}, \frac{3k}{2} - 2kS_0^2 + \frac{243(n\pi)^4}{2} - 162(n\pi)^4 S_0^2 \right), \\ A_4 &= (0, -3, 0, 3k - 37(n\pi)^4), \end{aligned} \qquad (42) $$

where $S_0 = \sin(n\pi s_0)$. This matrix has a determinant of

$$360\pi^4 n^4 (k + (n\pi)^4)^2 (-k + 9(n\pi)^4), \qquad (43)$$
which is non-zero for $k > 0$, except in the case $k = 9(n\pi)^4$. The source of this singularity is discussed below. The vector **b**, which results from the non-homogeneous terms, is given by

$$ \mathbf{b} = \left( \frac{\gamma(3 - 4S_0^2)}{8S_0^3}, \frac{-\gamma}{8S_0^3}, 0, 0 \right)^T \quad (44) $$

Note that this vector is zero when the non-linearity is absent, i.e., when $\gamma = 0$. The solution of this linear system of equations is

$$ \begin{aligned} \beta_1 &= \frac{\gamma(-3 + 4S_0^2)(-3k + 37(n\pi)^4)}{1280(n\pi)^4 S_0^3 (-k + 9(n\pi)^4)}, & \beta_3 &= \frac{\gamma(-3k + 37(n\pi)^4)}{1280(n\pi)^4 S_0^3 (-k + 9(n\pi)^4)}, \\ \eta_1 &= \frac{-3\gamma(-3 + 4S_0^2)}{1280(n\pi)^4 S_0^3 (-k + 9(n\pi)^4)}, & \eta_3 &= \frac{-3\gamma}{1280(n\pi)^4 S_0^3 (-k + 9(n\pi)^4)}. \end{aligned} \quad (45) $$

The series solution of the normal modes is now available up to third order in $u_0$ and $v_0$.

The non-linear normal modes feature amplitude dependent corrections which consist of the first and third spatial harmonics of the associated linear mode shapes. They are determined by substituting the solutions for the $\beta$'s and the $\eta$'s into equation (41), yielding $a_6$ and $a_8$, which are then used to obtain $b_7$ and $b_9$ from equations (38). The final result is given by $(n = 1, 2, 3, \dots)$.

$$ \begin{aligned} U_n &= u_0 a_1(s) \left( 1 - \gamma (a_1(s)^2 - 1) \frac{(-3k u_0^2 + 37(n\pi)^4 u_0^2 - 3v_0^2)}{320(n\pi)^4 (-k + 9(n\pi)^4)} \right) + \dots, \\ V_n &= v_0 a_1(s) \left( 1 - 3\gamma (a_1(s)^2 - 1) \frac{-k u_0^2 + 39(n\pi)^4 u_0^2 - v_0^2}{320(n\pi)^4 (-k + 9(n\pi)^4)} \right) + \dots. \end{aligned} \quad (46) $$

Several remarks are now in order. First, note that since $a_1(s_0) = 1$, the identities given in equation (4) are exactly satisfied. Also, by recalling the sin $(nns)$ nature of $a_1(s)$, it is seen that these modes contain spatial harmonics of first and third order, and that these contributions enter in an amplitude dependent manner—leading to amplitude dependent mode shapes. It is also worthy of note that the effect of the non-linearity on the shapes decreases for higher mode number. This is consistent with the observation that the elastic foundation has a lesser influence on higher modes. Finally, for $k = 9(n\pi)^4$, a singularity occurs in the modes; this is precisely where a 3:1 internal resonance occurs between the $n$th and the $3n$th modes. In such a case, an unremovable coupling exists between those modes, and the present formulation fails. This situation can be handled by constructing the four-dimensional invariant manifold on which the coupled two-mode dynamics take place. This same situation arises in finite-dimensional systems (see reference [2]), and the details of its resolution are left for future work.

Note also that the non-linear normal modes for this example contain only those terms in $u_0$ and $v_0$ which allow for synchronous modal motions, as expected for a conservative, non-gyroscopic system. The synchronicity is easily seen by observing that whenever $u_0$ is zero, that is, at any instant at which the displacement at $s_0$ is zero, the entire displacement field is zero. Similarly, whenever the velocity at $s_0$ is zero, $v_0 = 0$, the entire velocity field is zero. Thus, the non-linear normal mode motions are standing waves for the beam, and this makes for simple graphical representations of the non-linear normal modes. The presence of non-conservative and/or gyroscopic type terms will generally destroy this synchronicity, resulting in travelling wave motions for the non-linear normal modes, just
as in the linear case. The general formulation given in section 2 can accommodate these effects, but in practice obtaining solutions will be formidable, simply because the equations for the $a_i$'s and $b_i$'s will be more complicated.

It is also worth noting that the linear parameter $k$ influences the non-linear mode shapes, even though it has no effect on the linear mode shapes. Its influence becomes especially important when the parameters approach the internal resonance condition.

Mode shapes are most naturally observed by considering an instant at which the beam reaches its peak displacement; that is, when $v_0 = 0$. In order to plot a mode shape, one can simply set $v_0 = 0$ in the expression for $U_n$ (equation (46)), choose $n$ for the desired mode, select a value for $s_0$, select parameter values, and then plot the resulting function vs. $s$ for a set of amplitudes which are varied by changing $u_0$. For the present example we have taken $s_0 = 1/2n$ for the $n$th non-linear normal mode, which is an anti-node for each linear mode. In Figures 3(a)-(c) are shown non-linear mode shapes for the first three modes for $k=0$ and $\gamma = 1 \times 10^6$ for three different values of $u_0$. In these figures, the linear modes shapes, obtained by simply taking $\gamma = 0$, are shown as dashed lines. Also shown, in Figure 3(d), is the second mode for the case $k = 13\,800$ (which is close to the case of a 3:1 internal resonance between the second and sixth modes; that is, $k \approx 9 (2\pi)^4$) and $\gamma = 1 \times 10^6$. This shows the large modal distortion which occurs near the singularity in the second mode. In this case, the first and third modes are virtually unchanged from the $k=0$ case. (The large values of $\gamma$ and $k$ result from rescaling.) The upper amplitude ranges clearly demonstrate the effects of the non-linearity on the first mode. Also, as expected, the foundation has a less dramatic effect on the higher modes. In this example, the non-linearity is hardening ($\gamma > 0$), and note that, as expected, the non-linearity "flattens" the modes for large amplitudes. A softening non-linearity would accentuate the modal peaks.

<figure><img src="image_5.png" /><figcaption>Figure 3. The mode shapes for Example 2: non-linear normal mode shapes are shown as solid lines and the linear mode shapes are shown as dashed lines for comparison. The base point is at $s_0 = 1/2n$ for the $n$th mode. Peak amplitudes at $s_0$ are $u_{0,max} = 0.01666...$, $0.0333...$ and $0.0500$. In each non-linear case $\gamma = 1 \times 10^6$. The linear mode shapes are independent of the value of $k$: (a) the first mode for $k=0$; (b) the second mode for $k=0$; (c) the third mode for $k=0$; (d) the second mode for $k = 13\,800$, near the 3:1 internal resonance condition.</figcaption></figure>
The equations which govern the dynamics of the modes are obtained by substituting the solutions for $U_n$ and $V_n$ into the equations of motion and then evaluating them at $s = s_0$. The result in second order form in terms of $u_0$ defines the following *modal oscillators*:

$$ \ddot{u}_0 + ((n\pi)^4 + k)u_0 - \frac{\gamma}{16(-k + 9(n\pi)^4)} (u_0^3(9k - 111(n\pi)^4 + 4(k + (n\pi)^4)\sin(n\pi s_0)^2) \\ + u_0\dot{u}_0^2(-9 + 12\sin(n\pi s_0)^2)) + \cdots = 0. \quad (47) $$

It is initially puzzling to find $s_0$ appearing in the coefficients of this differential equation, since the characteristics of the motion seem to depend on the choice of $s_0$—clearly not a feature of a normal mode motion. In particular, the frequency correction from the non-linearity cannot depend on $s_0$ if a normal mode motion is to be synchronous. This perplexing aspect of the formulation is settled as one proceeds in the analysis. At this stage, the amplitude dependent frequencies of oscillations for the normal modes are computed *via* Linstedt's method. This procedure requires the amplitude of oscillation to be specified by a parameter—here the assumption is that, to first order, the motion is given by $u_0(t) \approx U_0 \sin(\omega_0 t)$, where $\omega_0 = \sqrt{((n\pi)^4 + k)}$ is the linear natural frequency of the nth mode. This leads to the non-linear frequency of oscillation of

$$ \omega_n = \sqrt{k + (n\pi)^4} \left( 1 + \frac{9U_0^2\gamma}{32(k + (n\pi)^4)\sin(n\pi s_0)^2} \right) + \cdots \quad (48) $$

Here it seems as if all hope is lost, since the frequency of oscillation appears to depend on $s_0$! However, note that it also depends on the amplitude $U_0$, and this salvages things, as can be seen by the following argument. Consider a purely modal motion for an odd numbered mode, which has a maximum mid-point displacement of magnitude A. For such a motion the amplitude of $U_0$ depends on $s_0$ as follows to leading order: $U_0 \approx A \sin(n\pi s_0)$. Therefore, the frequency of oscillation is given by

$$ \omega_n = \sqrt{k + (n\pi)^4} \left( 1 + \frac{9A^2\gamma}{32(k + (n\pi)^4)} \right) + \cdots \quad (49) $$

for any $s_0$, which clearly shows that the net frequency is independent on $s_0$. A similar argument holds for the even numbered modes. Note that this frequency is exactly that which is predicted by the procedure of using the nth *linear* mode shape as an assumed mode, projecting the non-linear equations of motion onto that mode, and then computing the frequency using a perturbation method. However, higher order frequency corrections will differ. Also, the resulting non-linear modal oscillators are not the same. This is an artifact of the co-ordinates used for the linear part of the modes. However, the present approach systematically generates the modal distribution due to non-linear effects, and these will play a role in higher order terms in the modal oscillators.

Finally, note that the terms in the modal oscillator are all conservative in nature. If the system were non-conservative at the linear order and had cubic non-linearities, the $a_2, a_7, a_9, b_1, b_6$ and $b_8$ terms would generally be non-zero, and would lead to $v_0, v_0^3$ and $v_0 u_0^2$ terms in the modal oscillators. The paper by Shaw and Pierre [2] provides an example which demonstrates this in the finite-dimensional case.
### 3.3. EXAMPLE 3: TRANSVERSE VIBRATIONS OF A SIMPLY SUPPORTED BEAM WITH A NON-LINEAR TORSIONAL SPRING AT EACH END

This example demonstrates the ability of the method to handle discrete non-linear elements and non-linearities in boundary conditions. Many of the details in the previous examples are skipped over and only the essential points are given in detail.

The system is depicted in Figure 4. It consists of a simply supported, linear beam with identical torsional springs attached at each end. These springs are taken to be inertia-free and purely non-linear; that is, they provide no moment at the linear approximation. This assumption preserves the linear mode shapes and may be thought of as an approximation for boundary condition imperfections. For the present case the springs are assumed to provide a restoring moment at each end which is proportional to the cube of the slope at that end.

The equation of motion for $s \in (0, 1)$ is linear and is the same as in Example 1. In this case, however, the boundary conditions on the displacement are given by

$$u(0, t) = u(1, t) = 0, \quad u_{xx}(0, t) = \alpha(u_s(0, t))^3, \quad u_{ss}(1, t) = -\alpha(u_s(1, t))^3, \quad (50)$$

where $\alpha$ is a positive constant. The boundary conditions on the velocity are determined by taking a time derivative of the above, which yields

$$\begin{align}
v(0, t) &= v(1, t) = 0, & v_{xx}(0, t) &= 3\alpha(u_s(0, t))^2 v_s(0, t), \nonumber \\
v_{ss}(1, t) &= -3\alpha(u_s(1, t))^2 v_s(1, t). \tag{51}
\end{align}$$

The corresponding boundary conditions for the $a_i$'s and $b_i$'s are obtained by assuming a modal motion, $u = U$ and $v = V$, directly substituting the series expansions (10) for $U$ and $V$ into the above, expanding in terms of $u_0$ and $v_0$, and gathering terms of like powers in $u_0$ and $v_0$. This procedure yields

$$a_j(0) = b_j(0) = 0, \quad a_j(1) = b_j(1) = 0, \quad j = 1, 2, 3, 4, \dots, \quad (52)$$

$$a_i''(0) = b_i''(0) = 0, \quad a_i''(1) = b_i''(1) = 0, \quad i = 1, 2, 3, 4, 5, \quad (53)$$

$$a_6''(0) - \alpha(a_1'(0))^3 = 0, \quad a_7''(0) - 3\alpha(a_1'(0))^2 a_2'(0) = 0,$$

$$a_8''(0) - 3\alpha(a_2'(0))^2 a_1'(0) = 0, \quad a_9''(0) - \alpha(a_2'(0))^3 = 0,$$

$$b_6''(0) - 3\alpha(a_1'(0))^2 b_1'(0) = 0,$$

$$b_7''(0) - 6\alpha a_1'(0) a_2'(0) b_1'(0) - 3\alpha (a_1'(0))^2 b_2'(0) = 0,$$

$$b_8''(0) - 3\alpha(a_2'(0))^2 b_1'(0) - 6\alpha a_1'(0) a_2'(0) b_2'(0) = 0,$$

$$b_9''(0) - 3\alpha(a_2'(0))^2 b_2'(0) = 0, \quad a_6''(1) + \alpha(a_1'(1))^3 = 0,$$

$$a_7''(1) + 3\alpha(a_1'(1))^2 a_2'(1) = 0, \quad a_8''(1) + 3\alpha(a_2'(1))^2 a_1'(1) = 0,$$

$$a_9''(1) + \alpha(a_2'(1))^3 = 0, \quad b_6''(1) + 3\alpha(a_1'(1))^2 b_1'(1) = 0,$$

$$b_7''(1) + 6\alpha a_1'(1) a_2'(1) b_1'(1) + 3\alpha (a_1'(1))^2 b_2'(1) = 0,$$

$$b_8''(1) + 3\alpha(a_2'(1))^2 b_1'(1) + 6\alpha a_1'(1) a_2'(1) b_2'(1) = 0,$$

$$b_9''(1) + 3\alpha(a_2'(1))^2 b_1'(1) = 0. \quad (54)$$

Since the differential equations for the $a_i$'s and $b_i$'s are identical to those in the previous example with $\gamma = 0$ and $k = 0$, they are not presented again. The linear part of this problem is identical to the first example and the same solutions are obtained for $a_1$, $a_2$, $b_1$ and $b_2$.
<figure><img src="image_6.png" /><figcaption>Figure 4. The physical system for Example 3: a simply supported beam with non-linear torsional springs attached at each end. The springs are purely cubic, with stiffness coefficient $\alpha$.</figcaption></figure>

Furthermore, the quadratic part of the problem is identical to the quadratic part of the second example and yields the trivial solution for $a_3, a_4, a_5, b_3, b_4$ and $b_5$.

The interesting part of this problem lies in the cubic order. The algebraic-differential equations for $a_6, a_7, a_8, a_9, b_6, b_7, b_8$ and $b_9$ are the same as those of the second example (equation (38)) with $\gamma$ and $k$ set equal to zero. With $\gamma = 0$, the equations are homogeneous. However, in this case, a non-trivial solution exists due to the non-homogeneous nature of the boundary conditions. The solution procedure in this case can be streamlined by using directly the conservative, non-gyroscopic nature of the problem and observing that the trivial solution for $a_7, a_8, b_6$ and $b_8$ satisfies four of the algebraic-differential equations and their attendant boundary conditions. As in the previous problem, the remaining four equations are easily reduced to two, since the $b_i$'s can be solved for in terms of the $a_i$'s directly from the non-differential equations. The result of this procedure is the following pair of homogeneous equations for $a_6$ and $a_8$:

$$ 
\begin{aligned}
3(n\pi)^4 a_6(s) - 2(n\pi)^6 a_8(s) - a_6'''(s) + a_1(s)a_6''''(s_0) &= 0, \\
-6a_6(s) + 7(n\pi)^4 a_8(s) - a_8'''(s) + a_1(s)a_8''''(s_0) &= 0,
\end{aligned} 
\quad (55) 
$$

which are the same as equation (40) from the second example with $\gamma = 0$ and $k = 0$. Using the known solutions at the linear level, the boundary conditions for $a_6$ and $a_8$ simplify, from equation (54),

$$ 
\begin{aligned}
a_6(0) = a_6(1) = a_8(0) = a_8(1) &= 0, & a_6''(0) &= \mu_n, \\
a_6''(1) = (-1)^{n+1}\mu_n, & a_8''(0) = a_8''(1) &= 0,
\end{aligned} 
\quad (56) 
$$

where

$$ \mu_n = \alpha \left( \frac{n\pi}{\sin(n\pi s_0)} \right)^3. $$

After $a_6$ and $a_8$ are obtained, $b_7$ and $b_9$ can be determined directly from equations (38).

In order to solve for $a_6$ and $a_8$, the problem is transformed to a non-homogeneous system of equations with homogeneous boundary conditions using a standard procedure [12, 13]. Since the non-homogeneity occurs for $a_6$, a new function $\bar{a}_6(s)$ is defined as follows:

$$ \bar{a}_6(s) = a_6(s) - \mu_n h(s), \quad (57) $$

with $h(s)$ chosen such that the system of equations in terms of $\bar{a}_6(s)$ and $a_8(s)$ has homogeneous boundary conditions. The choice for $h(s)$ is not unique: it must simply satisfy the boundary conditions

$$ h(0) = h(1) = 0, \quad h''(0) = 1, \quad h''(1) = (-1)^{n+1}, \quad (58) $$

in order to render the boundary conditions on $\bar{a}_6$ homogeneous. Here the static solution for the beam with concentrated moments at the ends is chosen for $h(s)$:

$$ h(s) = \left( \frac{-2 + (-1)^n}{6} \right) s + \frac{s^2}{2} - \left( \frac{1 + (-1)^n}{6} \right) s^3. \quad (59) $$
This change in functions yields a new system in terms of $\bar{a}_6(s)$ and $a_8(s)$, as follows:

$$3(n\pi)^4\bar{a}_6(s) - 2(n\pi)^8a_8(s) - \bar{a}_6'''(s) + a_1(s)a_6'''(s_0) = -3(n\pi)^4\mu_n h(s), \\ -6\bar{a}_6(s) + 7(n\pi)^4a_8(s) - a_8'''(s) + a_1(s)a_8'''(s_0) = 6\mu_n h(s), \quad (60)$$

where $h'''(s) = 0$ has been used. The boundary conditions on $\bar{a}_6(s)$ are

$$\bar{a}_6(0) = \bar{a}_6(1) = 0, \quad \bar{a}_6''(0) = \bar{a}_6''(1) = 0. \quad (61)$$

The solutions for $\bar{a}_6(s)$ and $a_8(s)$ can now be obtained by Fourier sine series, as follows:

$$\bar{a}_6(s) = \sum_{i=1}^{\infty} \beta_i \sin(i\pi s), \quad a_8(s) = \sum_{i=1}^{\infty} \eta_i \sin(i\pi s), \quad (62)$$

Note that in this case an infinite series is required for an "exact" solution at this order. However, since the present approach is asymptotic, only a finite number of terms are retained in practice. The equations for $\beta_i$ and $\eta_i$ are linear, and are obtained by projecting onto the functions $\sin(j\pi s)$. This procedure has a twist that needs to be pointed out. The term $a_1(s)a_6'''(s_0)$ in the first equation and its counterpart $a_1(s)a_8'''(s_0)$ in the second equation yield and interesting result when they are projected onto $\sin(j\pi s)$. First, note that the $a_6'''(s_0)$ term contains all the $\beta_i$'s and that these are not directly eliminated by the projection onto $\sin(j\pi s)$, since this term does not depend on $s$ (it is evaluated at $s_0$). However, the $a_1(s)$ term is proportional to $\sin(n\pi s)$ and is therefore orthogonal to all $\sin(j\pi s)$ functions, except in the case $j=n$. Thus, the equations for $\beta_j$ and $\eta_j$ are pairwise uncoupled for all $j \neq n$, but are completely coupled for $j=n$. Since only the $j=n$ pair of equations contains all the $\beta_i$'s and $\eta_i$'s, the solution can be easily obtained as follows. First, for all $j \neq n$ the equations for $\beta_j$ and $\eta_j$ are pairwise uncoupled and easily solved. Their solution can then be substituted into the $j=n$ equation in order to determine $\beta_n$ and $\eta_n$. This procedure is now described in more detail.

The resulting equations to be solved for $\beta_j$ and $\eta_j$ are determined to be

$$-(j\pi)^4\beta_j + 3(n\pi)^4\beta_j - 2(n\pi)^8\eta_j + \frac{\delta_{nj}}{\sin(n\pi s_0)} \sum_{i=1}^{\infty} \beta_i(i\pi)^4 \sin(i\pi s_0) = \frac{6\mu_n \pi n^4}{j^3} (1 + (-1)^{j+n}), \\ -(j\pi)^4\eta_j + 7(n\pi)^4\eta_j - 6\beta_j + \frac{\delta_{nj}}{\sin(n\pi s_0)} \sum_{i=1}^{\infty} \eta_i(i\pi)^4 \sin(i\pi s_0) = \frac{-12\mu_n}{\pi^3 j^3} (1 + (-1)^{j+n}), \quad (63)$$

where $\delta_{nj}$ is the Kronecker delta. For $j \neq n$ these equations are pairwise uncoupled and have the following solution:

$$\beta_j = \frac{-6\mu_n n^4 (j^4 - 3n^4) (1 + (-1)^{j+n})}{\pi^3 j^3 (j^4 - n^4) (j^4 - 9n^4)}, \quad \eta_j = \frac{12\mu_n j (1 + (-1)^{j+n})}{\pi^7 (j^4 - n^4) (j^4 - 9n^4)}. \quad (64)$$

Note that the rate of convergence for the infinite series for $\bar{a}_6$ and $a_8$ is rapid since the terms decay at the rate $j^{-7}$. Also, note that each odd (respectively, even) numbered mode contains only odd (respectively, even) order spatial harmonics.

The solution for the $j=n$ case is obtained by setting $j=n$ and solving the resulting equations for $\beta_n$ and $\eta_n$ in terms of the $j \neq n$ solutions. The result is

$$\beta_n = \frac{1}{9(n\pi)^4} \left( 36\pi n \mu_n - \frac{1}{\sin(n\pi s_0)} \sum_{i=1, i \neq n}^{\infty} (7\beta_i + 2(n\pi)^4\eta_i)(i\pi)^4 \sin(i\pi s_0) \right), \\ \eta_n = \frac{-1}{3(n\pi)^8 \sin(n\pi s_0)} \sum_{i=1, i \neq n}^{\infty} (2\beta_i + (n\pi)^4\eta_i)(i\pi)^4 \sin(i\pi s_0). \quad (65)$$
The non-linear normal modes are now available and are given by

$$ 
\begin{aligned}
U_n &= \frac{\sin(n\pi s)}{\sin(n\pi s_0)} u_0 + \left[ \sum_{i=1}^{\infty} \beta_i \sin(i\pi s) + \mu_n h(s) \right] u_0^3 + \left[ \sum_{i=1}^{\infty} \eta_i \sin(i\pi s) \right] u_0 v_0^2 + \dots, \\
V_n &= \frac{\sin(n\pi s)}{\sin(n\pi s_0)} v_0 + b_7 u_0^2 v_0 + b_9 v_0^3 + \dots,
\end{aligned} 
\quad (66) 
$$

where the expressions for $h$, $\mu_n$, $\beta_i$ and $\eta_i$ are given above. It is important to remember that the $i=n$ term in the infinite series is special. The expressions for $b_7$ and $b_9$ are given in terms of $a_6$ and $a_8$ in equation (38) and are not presented in detail here.

Since, in practice, one must truncate the infinite series involved in the solution, it is useful to consider their convergence properties. It has been noted above that the terms in the series for $\tilde{a}_6$ and $a_8$ are proportional to $j^{-7}$; this implies good convergence properties for the non-linear correction terms in the mode shapes. However, the terms in the series for $\beta_n$ and $\eta_n$ are proportional to $i^{-3}$ due to the $i^4$ term multiplying $\beta_i$ and $\eta_i$ for $i \neq n$. Thus, in order to achieve a desired accuracy for a given mode, several terms must be taken in the infinite series appearing in the solution for $\beta_n$ and $\eta_n$. However, consistent accuracy can be maintained by retaining fewer terms in the series for $\tilde{a}_6$ and $a_8$ than the number of terms taken in the series for $\beta_n$ and $\eta_n$. This fact was not exploited in the calculations presented, but it may prove to be useful in some applications.

Plots of the non-linear mode shapes are determined by the process described for the previous example, except that in the present case the Fourier series must be truncated. Again $s_0 = 1/2n$ is chosen. In the calculations, $N$ terms have been used for each series in general. The first three modes are generated with the following number of terms retained: for the first mode $N=8$, for the second mode $N=12$, and for the third mode $N=16$. The expressions for $U_n$ for $n=1, 2, 3$ are given in Appendix B for the stated number of terms. In Figure 5 are shown each of the first three linear (dashed lines) and non-linear (solid lines) mode shapes for two different amplitudes of $u_0$ with $\alpha = 200$. Smaller amplitudes are used for higher modes, since the non-linearity has a more dramatic effect on higher modes. The results presented were compared with those obtained using twice the number of terms in Fourier series for each case and there was no discernible difference in the mode shapes at the amplitudes shown in Figure 5, thereby confirming satisfactory convergence of the series solutions.

Note that the non-linear spring has exactly the expected effect on the first and second modes: while, for small amplitudes, the beam behaves in a linear manner, the beam becomes more restrained in its end rotations for large amplitudes. The effect on the third mode is different—the primary effect of the springs is to decrease the mid-span deflection. This effect is the opposite to that initially expected, and the explanation for it is unknown. (In fact, linear torsional springs have the same effect on the third linear mode.) In each case the amount of distortion is exactly that required to maintain synchronous motion for a mode, and it has been determined in a systematic manner based on first principles.

The modal oscillators are determined in the usual manner, in fact, from equation (27) from the first example (since the equations of motion for $s \in (0, 1)$ are the same) with the $U_n$ determined from the present example. The result is

$$ \ddot{u}_0 + (n\pi)^4 u_0 + \left[ \sum_{i=1}^{\infty} \beta_i (i\pi)^4 \sin(i\pi s_0) \right] u_0^3 + \left[ \sum_{i=1}^{\infty} \eta_i (i\pi)^4 \sin(i\pi s_0) \right] u_0 v_0^2 + \dots = 0, \quad (67) $$
<figure><img src="image_7.png" /><figcaption>Figure 5. The mode shapes for Example 3: non-linear normal mode shapes are shown as solid lines and the linear mode shapes are shown as dashed lines for comparison. The base point is at $s_0 = 1/2n$ for the nth mode. In each non-linear case $\alpha = 200$. (a) The first mode for peak amplitudes at $s_0$ of $u_{0,max} = 0.025$ and $0.050$; (b) the second mode for peak amplitudes at $s_0$ of $u_{0,max} = 0.015$ and $0.030$; (c) the third mode for peak amplitudes at $s_0$ of $u_{0,max} = 0.0075$ and $0.0150$.</figcaption></figure>

where $h'''(s_0) = 0$ has been used. The modal oscillators and the frequency corrections for a specific mode are most easily presented by evaluating the modal oscillators at a specific $s_0$, which is again taken to be $1/2n$. The resulting oscillators and their amplitude dependent frequencies are as follows:

for the first mode (with $N = 8$ terms),

$$ \ddot{u}_0 + \pi^4 u_0 + \alpha \left( \frac{296\,431\pi^4}{69\,069} u_0^3 - \frac{20\,155}{23\,023} u_0 \dot{u}_0^2 \right) + \cdots = 0, \qquad \omega_1 = \pi^2 + \frac{3A^2\pi^2\alpha}{2} + \cdots; \tag{68} $$

for the second mode (with $N = 12$ terms),

$$ \ddot{u}_0 + 16\pi^4 u_0 + \alpha \left( \frac{15\,776\pi^4}{231} u_0^3 - \frac{62}{77} u_0 \dot{u}_0^2 \right) + \cdots = 0, \qquad \omega_2 = 4\pi^2 + 6A^2\pi^2\alpha + \cdots; \tag{69} $$
for the third mode (with $N = 16$ terms),

$$ 
\begin{aligned}
& \ddot{u}_0 + 81\pi^4 u_0 + \alpha \left( \frac{2\,444\,117\,969\,493\pi^4}{2\,988\,690\,614} u_0^3 - \frac{54\,658\,600\,391}{2\,988\,690\,614} u_0 \dot{u}_0^2 \right) + \cdots = 0, \\
& \omega_3 = 9\pi^2 + \frac{27A^2\pi^2\alpha}{2} + \cdots.
\end{aligned}
\quad (70) 
$$

$A$ is the amplitude of oscillation at $s = s_0$.

This non-linearity is hardening for $\alpha > 0$ and has the expected effect on the frequency. Note that the non-linearity has a greater influence on higher modes for a given amplitude, although smaller amplitudes will generally occur for higher modes. Again these modal oscillators are different from those obtained by projecting the equations of motion onto the linear modes, but the non-linear frequencies obtained from these two approaches are identical, at least to third order in amplitude.

### 4. CONCLUSIONS AND SOME DIRECTIONS FOR FUTURE WORK

This work represents, to the authors' knowledge, the first systematic approach to the definition and generation of normal modes of motion for non-linear, continuous vibratory systems. The methodology is based on the concept of invariant manifolds for dynamical systems and it is constructive for weakly non-linear systems. Using asymptotic series expansions, it provides the physical nature of the non-linear mode shapes and the associated modal dynamics. The formulation recovers the well-known linear eigenmodes when non-linearities are absent, although the approach is completely different from the standard eigenvalue problem formulation. In fact, the non-linear normal modes should be thought of as the source of the linear modes, not merely as extensions of them.

An important distinction between the present approach and those which use harmonic balance and eigenfunction expansions is that, while both approaches lead to mode shapes which depend on the peak amplitude of motion, the formulation presented herein allows for changes in the beam shape *during a given motion*. This is evident by considering that the solution for $u_0$ will, for the systems considered, be a periodic function of time, and therefore the relative contributions of the various linear mode shapes to $U$ and $V$ from $u_0$ and $v_0$ will also be time-periodic. This is not accounted for in the harmonic balance approach, and the two methods lead to different expressions for the non-linear mode shapes. Since the invariant manifold technique is less restrictive and is based on first principles of dynamical systems, the shapes predicted by this approach will more accurately represent those of non-linear systems undergoing single mode motions.

Although a limited class of problems has been considered in this paper (those which are governed by equations of the form given by equation (1)), it should be clear that the basic approach is very general and can be applied to a wider class of problems. The examples presented here demonstrate the power of the approach in handling a variety of non-linearities in a very systematic manner. It should be noted that the calculations can become quite involved, especially if the boundary conditions are less simple than those considered in the examples, and/or if non-conservative and/or gyroscopic terms are present. However, with widely available computer-assisted symolic manipulators, many problems of practical interest can be solved by this method. Currently, the authors are determining the normal modes for finite deformations of a cantilever beams (Hsieh et al. [10]), a problem that involves another complication: a non-linear inertia operator. While the calculations are cumbersome, they are manageable and solutions can be obtained.
It is also important to point out that the range of validity of the asymptotic solutions is limited, and therefore there exists an intermediate range of ($u_0, v_0$) amplitudes over which non-linear effects are correctly captured. For very small amplitudes, the linear and non-linear mode shapes are indistinguishable, and as amplitudes are increased, the non-linear effects start to become evident. However, at some point the non-linear mode shapes predicted by the series solutions begin to differ from the shape required for synchronous motion. Estimates for this range of validity can be obtained by comparing the mode shapes generated by selecting different values for $s_0$ for a given beam amplitude. These shapes should, of course, be nearly the same, and are so over a range of amplitudes. However, as the amplitude is increased, these shapes begin to diverge from one another, signalling the breakdown of the approximation. In all of the examples presented here, the mode shapes were compared for at least three different values of $s_0$ at the maximum amplitude shown in order to verify the validity of the approximation.

In the applied mathematics literature there are some results which are pertinent to the present work, two of which are the following. First, the invariant manifolds which are used here to represent normal modes of motion are known as the "standard foliation" of the unstable, center and stable manifolds of a dynamical system near an equilibrium point. Mathematical proofs of their existence for the finite-dimensional case can be found in the paper by Fenichel [14], and certain non-degeneracy conditions on the eigenvalues are required. These conditions correspond precisely to the internal resonance conditions encountered in Example 2. Second, there exists a recently developed class of computational methods, known as *non-linear Galerkin techniques*, which also make use of invariance properties of modes. That approach is substantially different in spirit to the present one; it is more computationally oriented and is used for the simulation of large-scale fluid mechanics problems [15].

Several extensions, generalizations and applications of the methodology described in this paper remain to be tackled, and the authors are currently engaged in work on several of them. The most significant ones are mentioned below.

Although in all example problems the normal mode invariant manifolds have been observed to reduce to the linear eigenmodes in the linearized case, a general proof of the equivalence of the two approaches for linear systems needs to be formulated. This proof was recently completed by the authors, and will be presented in a paper that emphasizes the features of an invariant manifold approach to normal modes and modal analysis for linear systems [9].

A generalization of the present method may offer substantial benefits in the area of model reduction for non-linear continuous systems. Reduced order models are typically obtained by performing a linear modal analysis of the non-linear system (via orthogonal projection) and subsequently ignoring modes that are deemed non-essential to the non-linear dynamics. This truncation procedure, however, inherently neglects the two-way exchange of energy, or *contamination*, between the modes of interest and the ignored modes, which is automatically generated by the projection of the *non-linear* equations of motion onto the *linear* modes. The neglect of these interactions causes deteriorating convergence and accuracy of the modal analysis procedure. The non-linear normal modes defined in this paper, however, allow for the definition of uncontaminated single mode non-linear models. By generalizing the procedure in such a manner that invariant manifolds of dimension 2N are generated, one could obtain equations in which contamination between the N modes of interest and the ignored modes is eliminated, while the modal coupling within the modelled modes is retained. This procedure would therefore allow for the generation of very "clean," uncontaminated, reduced order models in a
systematic way. This may have important implications in structural dynamics and in the
design and implementation of control systems.

As it stands now, the method cannot handle internal resonances. In this situation, there
exists unremovable coupling between the interacting modes, which renders the approach
useless for those modes. In the case of *M* resonant modes, a 2*M*-dimensional manifold
is required in order to capture the dynamics of the interacting modes. Such manifolds
can be obtained by expressing the response of the system as a function of the displace-
ment and velocities of *M* points in the system and generalizing the procedures presented
above. This is similar in spirit to the generalization described immediately above, and
will no doubt involve substantial calculations, even for the simplest cases of internal
resonance.

In reference [2], the authors presented a method for non-linear, finite-dimensional
systems which is analogous to the one presented herein. In that work, several simulations
were carried out which demonstrated the invariant nature of non-linear normal
modes developed from an invariant manifold approach. Furthermore, a non-linear
co-ordinate transformation which relates the physical system co-ordinates to the non-
linear modal co-ordinates was defined and used in the examples. This co-ordinate
transformation and demonstration of invariance are not so straightforward for continuous
systems. These topics are currently under investigation and will be presented in a future
paper.

Another class of non-linear systems of interest are those with piecewise linear forces, such as systems with dry friction or clearances. The asymptotic series approximate obviously fails in this situation, because the invariant manifolds corresponding to the normal modes are not smooth surfaces, but pieces of planes and/or surfaces joined together along the boundaries of the various linear regions. The solution of the equations defining the geometry of such manifolds may be attainable by employing matching methods.

In some situations, such as those involving systems with special symmetries, the
equations for the invariant manifold geometry may possess global solutions, which
would alleviate the need for local series expansions. Much of the past work on normal
modes for non-linear discrete systems has been centered on such systems [3, 5]. It is
expected that some types of non-linearities for continuous systems will also yield such
global solutions. For example, the equations of motion for a fixed-fixed beam with
a mid-line stretching non-linearity are known to yield invariant dynamics in terms of
modal amplitudes when a Galerkin procedure is performed using the linear mode shapes
as trial functions (see references [7] or [16]). This is a good indication that the linear mode
shapes will provide a global solution for the normal mode manifold equations for that
system.

The stability of non-linear modal motions is also of interest. It is well known that when
amplitudes become large, the normal modes can become unstable and bifurcate [5]. This
issue is important in the consideration of reduced order models since one needs to include
a sufficient number of modes so as to ensure the overall stability of the reduced model:
that is, the system should be stable to small disturbances in the unmodeled modes. Such
stability considerations can be examined using methods such as averaging.

Finally, the authors have recently developed on alternative approach to this class
of problems which first utilizes a usual Galerkin projection of the non-linear partial
differential equation (1) onto the linear modes shapes [17]. The method described in
reference [2] is then applied to the infinite set of coupled non-linear ordinary differential
equations, leading to results which are equivalent to those presented here. That formu-
lation does, however, avoid some of the subtleties encountered in the present procedure.
# ACKNOWLEDGMENTS

This work was supported in part by grants from NSF and ARO. The authors are grateful to Steve Wiggins of Caltech for pointing out the results of Fenichel [14] and to the work on non-linear Galerkin techniques [15].

# REFERENCES

1. W. WEAVER, S. P. TIMOSHENKO and D. H. YOUNG 1990 *Vibration Problems in Engineering*. New York: John Wiley.

2. S. W. SHAW and C. PIERRE 1993 *Journal of Sound and Vibration* **164**, 85-124. Normal modes for non-linear vibratory systems.

3. R. M. ROSENBERG 1966 *Advances in Applied Mechanics* **9**, 155-242. On non-linear vibrations of systems with many degrees of freedom.

4. R. H. RAND 1974 *International Journal of Non-linear Mechanics* **9**, 363-368. A direct method for non-linear normal modes.

5. A. VAKAKIS 1990 *Ph.D. dissertation, California Institute of Technology*. Analysis and identification of linear and non-linear normal modes in vibrating systems.

6. M. M. K. BENNOUNA and R. G. WHITE 1984 *Journal of Sound and Vibration* **96**, 309-331. The effects of large vibration amplitudes on the fundamental mode shape of a clamped-clamped uniform beam.

7. R. BENAMAR, M. M. K. BENNOUNA and R. G. WHITE 1991 *Journal of Sound and Vibration* **149**, 179-195. The effect of large vibration amplitudes on the mode shapes and natural frequencies of thin elastic structures, part I: simply supported and clamped-clamped beams.

8. W. SZEMPLINSKA 1990 *The Behaviour of Non-linear Vibrating Systems, Vol.* 2. Dordrecht: Kluwer.

9. C. PIERRE and S. W. SHAW 1993 An invariant manifold approach to modal analysis of linear systems (in preparation).

10. S. R. HSIEH, S. W. SHAW and C. PIERRE 1993 *International Journal of Solids and Structures* (to appear). Normal modes for large amplitude vibrations of a cantilever beam.

11. J. CARR 1981 *Applications of Centre Manifold Theory*. New York: Springer-Verlag.

12. R. L. BISPLINGHOFF and H. ASHLEY 1962 *Principles of Aeroelasticity*. New York: John Wiley.

13. K. GU and B. H. TONGUE 1987 *American Society of Mechanical Engineers, Journal of Applied Mechanics* **54**, 904-909. A method to improve the modal convergence for structures with external forcing.

14. N. FENICHEL 1971 *Indiana University Mathematical Journal* **21**, 193-255. Persistence and smoothness of invariant manifolds for flows.

15. M. MARION and R. TEMAM 1989 *SIAM Journal on Numerical Analysis* **26**, 1139-1157. Non-linear Galerkin methods.

16. S. W. SHAW 1988 *Journal of Sound and Vibration* **24**, 329-343. Chaotic dynamics of a slender beam rotating about its longitudinal axis.

17. S. W. SHAW and C. PIERRE 1992 American Society of Mechanical Engineers Volume DE-50, *Nonlinear Vibrations*, 1-5. On non-linear normal modes.

# APPENDIX A: ADDITIONAL IDENTITIES FOR THE $a_i$'S AND $b_i$'S

Consider the motion of a body undergoing a purely modal motion and let two base points, $\bar{s}_0$ and $\bar{s}_0$, be selected. These will play the roles of two values for $s_0$ in the formulation. Using these two points, the modal motion can be written either as

$$ 
\begin{aligned}
u(s, t) &= \bar{U}(\bar{u}_0(t), \bar{v}_0(t), s, \bar{s}_0) \\
&= a_1(s, \bar{s}_0)\bar{u}_0(t) + a_2(s, \bar{s}_0)\bar{v}_0(t) + \dots, \\
v(s, t) &= \bar{V}(\bar{u}_0(t), \bar{v}_0(t), s, \bar{s}_0) \\
&= b_1(s, \bar{s}_0)\bar{u}_0(t) + b_2(s, \bar{s}_0)\bar{v}_0(t) + \dots,
\end{aligned} 
\quad (A1) 
$$
or as

$$
\begin{align*}
u(s, t) &= \hat{U}(\hat{u}_0(t), \hat{v}_0(t), s, \hat{s}_0) \\
&= a_1(s, \hat{s}_0)\hat{u}_0(t) + a_2(s, \hat{s}_0)\hat{v}_0(t) + \dots, \\
v(s, t) &= \hat{V}(\hat{u}_0(t), \hat{v}_0(t), s, \hat{s}_0) \\
&= b_1(s, \hat{s}_0)\hat{u}_0(t) + b_2(s, \hat{s}_0)\hat{v}_0(t) + \dots, \tag{A2}
\end{align*}
$$

where $\bar{u}_0(t) = u(\bar{s}_0, t)$, etc.

Now $\bar{U}$ and $\bar{V}$ are evaluated at $s = \bar{s}_0$, and $\hat{U}$ and $\hat{V}$ are evaluated at $s = \bar{s}_0$. This yields
the following:

$$
\begin{align*}
u(\bar{s}_0, t) &= \bar{U}(\bar{v}_0(t), \bar{v}_0(t), \bar{s}_0, \bar{s}_0) = \bar{u}_0(t) \\
&= a_1(\bar{s}_0, \bar{s}_0)\bar{u}_0(t) + a_2(\bar{s}_0, \bar{s}_0)\bar{v}_0(t) + \dots, \tag{A3}
\end{align*}
$$

$$
\begin{align*}
v(\bar{s}_0, t) &= \bar{V}(\bar{u}_0(t), \bar{v}_0(t), \bar{s}_0, \bar{s}_0) = \bar{v}_0(t) \\
&= b_1(\bar{s}_0, \bar{s}_0)\bar{u}_0(t) + b_2(\bar{s}_0, \bar{s}_0)\bar{v}_0(t) + \dots, \tag{A4}
\end{align*}
$$

$$
\begin{align*}
u(\bar{s}_0, t) &= \hat{U}(\hat{u}_0(t), \hat{v}_0(t), \bar{s}_0, \hat{s}_0) = \bar{u}_0(t) \\
&= a_1(\bar{s}_0, \hat{s}_0)\hat{u}_0(t) + a_2(\bar{s}_0, \hat{s}_0)\hat{v}_0(t) + \dots, \tag{A5}
\end{align*}
$$

$$
\begin{align*}
v(\bar{s}_0, t) &= \hat{V}(\hat{u}_0(t), \hat{v}_0(t), \bar{s}_0, \hat{s}_0) = \bar{v}_0(t) \\
&= b_1(\bar{s}_0, \hat{s}_0)\hat{u}_0(t) + b_2(\bar{s}_0, \hat{s}_0)\hat{v}_0(t) + \dots. \tag{A6}
\end{align*}
$$

Next, the series defined for $\hat{u}_0$ and $\hat{v}_0$ in terms of $\bar{u}_0$ and $\bar{v}_0$ in equations (A3) and (A4), respectively, are substituted for $\hat{u}_0$ and $\hat{v}_0$ in the series defined in equations (A5) and (A6), and the results are expanded out to third order in $\bar{u}_0$ and $\bar{v}_0$. Coefficients of like powers in $\bar{u}_0^k \bar{v}_0^m$ are then collected and equated. Including all terms to third order, the result is a set of 18 non-linear algebraic identities relating the $a_i(\bar{s}_0, \hat{s}_0)$'s, $a_j(\bar{s}_0, \bar{s}_0)$'s, $b_j(\bar{s}_0, \hat{s}_0)$'s and $b_j(\hat{s}_0, \bar{s}_0)$'s. Since these are rather lengthy expression in the general case, only the simplified versions for the conservative, non-gyroscopic case with cubic non-linearities are presented here. In this case, $a_i = 0$ for $i = 2, 3, 4, 5, 7$ and $b_j = 0$ for $j = 1, 3, 4, 5, 6, 8$, and $b_2 = a_1$, which leaves the following non-trivial identities:

$$
\begin{align*}
1 &= a_1(\bar{s}_0, \hat{s}_0)a_1(\hat{s}_0, \bar{s}_0), \\
0 &= a_6(\bar{s}_0, \hat{s}_0)a_1^3(\hat{s}_0, \bar{s}_0) + a_1(\bar{s}_0, \hat{s}_0)a_6(\hat{s}_0, \bar{s}_0), \\
0 &= a_8(\bar{s}_0, \hat{s}_0)a_1^3(\hat{s}_0, \bar{s}_0) + a_1(\bar{s}_0, \hat{s}_0)a_8(\hat{s}_0, \bar{s}_0), \\
0 &= b_7(\bar{s}_0, \hat{s}_0)a_1^3(\hat{s}_0, \bar{s}_0) + a_1(\bar{s}_0, \hat{s}_0)b_7(\hat{s}_0, \bar{s}_0), \\
0 &= b_9(\bar{s}_0, \hat{s}_0)a_1^3(\hat{s}_0, \bar{s}_0) + a_1(\bar{s}_0, \hat{s}_0)b_9(\hat{s}_0, \bar{s}_0). \tag{A7}
\end{align*}
$$

Note that since these hold for all $\bar{s}_0$ and $\hat{s}_0$, they are true for the general case of $\bar{s}_0 = s$
and $\hat{s}_0 = s_0$. Therefore, these identities can be used for checking the validity of approximate
solutions for the $a(s, s_0)$'s and $b(s, s_0)$'s. In this general case, the first identity states simply
that

$$
1 = a_1(s, s_0)a_1(s_0, s), \tag{A8}
$$

which implies that if one simply switches $s$ and $s_0$ in $a_1(s, s_0)$, the inverse of $a_1(s, s_0)$ is
obtained. Note also that for $\bar{s}_0 = \hat{s}_0 = s_0$, these identities are consistent with those given in
equation (11).
More identities can be obtained by extending this process to include more $s_0$ points and/or more terms in the expansions.

APPENDIX B: NORMAL MODES FOR EXAMPLE 3

The first mode approximation using eight terms is

$$
\begin{equation}
\begin{aligned}
U_1 = u_0 \sin(\pi s) &+ \alpha \left[ u_0^3 \left( \frac{\pi^3 s}{2} (s-1) + \frac{213\,633\,857}{55\,255\,200} \sin(\pi s) \right. \right. \\
& \left. - \frac{13}{2160} \sin(3\pi s) - \frac{311}{2\,002\,000} \sin(5\pi s) - \frac{1199}{82\,045\,600} \sin(7\pi s) + \dots \right) \\
& + u_0 v_0^2 \left( \frac{255\,019}{18\,418\,400\pi^4} \sin(\pi s) + \frac{1}{80\pi^4} \sin(3\pi s) \right. \\
& \left. + \frac{5}{16\,016\pi^4} \sin(5\pi s) + \frac{7}{239\,200\pi^4} \sin(7\pi s) + \dots \right) \right] + \dots
\end{aligned}
\tag{B1}
\end{equation}
$$

The second mode approximation using 12 terms is

$$
\begin{equation}
\begin{aligned}
U_2 = u_0 \sin(2\pi s) &+ \alpha \left[ u_0^3 \left( \frac{4\pi^3}{3}(-s+3s^2-2s^3) \right. \right. \\
& \left. + \frac{58\,228}{15\,015} \sin(2\pi s) - \frac{13}{70} \sin(4\pi s) - \frac{13}{2160} \sin(6\pi s) - \frac{253}{335\,920} \sin(8\pi s) \right. \\
& \left. - \frac{311}{2\,002\,000} \sin(10\pi s) - \frac{431}{9\,999\,990} \sin(12\pi s) + \dots \right) \\
& + u_0 v_0^2 \left( \frac{61}{80\,080\pi^4} \sin(2\pi s) + \frac{1}{35\pi^4} \sin(4\pi s) + \frac{1}{1280\pi^4} \sin(6\pi s) \right. \\
& \left. + \frac{2}{20\,995\pi^4} \sin(8\pi s) + \frac{5}{256\,256\pi^4} \sin(10\pi s) \right. \\
& \left. + \frac{1}{185\,185\pi^4} \sin(12\pi s) + \dots \right) \right] + \dots
\end{aligned}
\tag{B2}
\end{equation}
$$

The third mode approximation using 16 terms is

$$
\begin{equation}
\begin{aligned}
U_3 = u_0 \sin(3\pi s) &+ \alpha \left[ u_0^3 \left( \frac{27\pi^3 s}{2} (s-1) + \frac{793\,881}{7280} \sin(\pi s) \right. \right. \\
& + \frac{362\,712\,763\,007\,451\,569}{125\,889\,386\,947\,658\,880} \sin(3\pi s) \\
& + \frac{1\,253\,151}{884\,000} \sin(5\pi s) - \frac{7\,079\,319}{166\,313\,840} \sin(7\pi s) - \frac{13}{2160} \sin(9\pi s) \\
& \left. - \frac{47\,232\,639}{33\,700\,707\,040} \sin(11\pi s) - \frac{92\,897\,199}{217\,682\,978\,240} \sin(13\pi s) \right]
\end{aligned}
\tag{B3}
\end{equation}
$$
$$ -\frac{311}{2002000}\sin(15\pi s) + \cdots) + u_0 v_0^2 \left( \frac{81}{7280\pi^4}\sin(\pi s) \right. \\ + \frac{80958693571253723}{3399013447586789760\pi^4}\sin(3\pi s) - \frac{405}{7072\pi^4}\sin(5\pi s) \\ + \frac{567}{484880\pi^4}\sin(7\pi s) + \frac{1}{6480\pi^4}\sin(9\pi s) + \frac{891}{25319840\pi^4}\sin(11\pi s) \\ \left. + \frac{1053}{99081920\pi^4}\sin(13\pi s) + \frac{5}{1297296\pi^4}\sin(15\pi s) + \cdots \right) \right] + \cdots \quad (B3) $$

For a check on the accuracy of the approximation, the values at $s = s_0 = 1/2n$ are evaluated and are found to be

$$ \begin{aligned} U_1(s_0 = \frac{1}{2}) &= u_0 + \left( \frac{4482544}{1157625} - \frac{\pi^3}{8} \right) \alpha u_0^3 + \cdots \\ &= u_0 - 0.0036 \alpha u_0^3 + \cdots, \\ U_2(s_0 = \frac{1}{4}) &= u_0 + \left( \frac{13108}{3375} - \frac{\pi^3}{8} \right) \alpha u_0^3 + \cdots \\ &= u_0 + 0.0081 \alpha u_0^3 + \cdots, \\ U_3(s_0 = \frac{1}{6}) &= u_0 + \left( \frac{196819966743566}{3385135128375} - \frac{15\pi^3}{8} \right) \alpha u_0^3 + \cdots \\ &= u_0 + 0.0057 \alpha u_0^3 + \cdots, \end{aligned} \quad (B4) $$

which are sufficiently accurate for reasonable values of $u_0$. As more terms in the Fourier series are taken, these all converge to $u_0$; that is $a_6(s_0)$ approaches zero.