A theoretical and computational model has been developed for the nonlinear motion of an inextensible beam undergoing large deflections for cantilevered and free–free boundary conditions. The inextensibility condition was enforced through a Lagrange multiplier which acted as a constraint force. The Rayleigh–Ritz method was used by expanding the deflections and the constraint force in modal series. Lagrange's equations were used to derive the equations of motion of the system, and a fourth-order Runge–Kutta solver was used to solve them. Comparisons for the cantilevered beam were drawn to experimental and computational results previously published and show good agreement for responses to both static and dynamic point forces. Some physical insights into the cantilevered beam response at the first and second resonant modes were obtained. The free–free beam condition was investigated at the first and third resonant modes, and the nonlinearity (primarily inertia) was shown to shift the resonant frequency significantly from the linear natural frequency and lead to hysteresis in both modes.

## Introduction

Nonlinear responses of the inextensible beam have been studied in various contexts and configurations in the literature. Lacarbonara [1] provides a contemporary overview of and includes several methods for enforcing inextensibility for several boundary conditions. While the cantilevered beam configuration has been more widely studied [2–10], there is relatively little in the literature on the free–free configuration. The theoretical solution to the free–free case is of great interest since it is difficult to model experimentally. With the advent of lighter composite materials used for wings and rotors, and the increasing need for inflatable or deployable spacecraft heat shields and panels, both the cantilever and free–free beam configurations are fundamental to understanding these new system dynamics. This work compares the results of the current model to that of other published cantilever models and extends the model to the free–free configuration.

Several studies have used Hamilton's principle and Lagrange's equation to derive the equations of motion for a beam. Three such studies utilize a change of coordinates to a Lagrangian coordinate system such that the axes are fixed to the beam throughout its deflection [2,5,11]. They also use analytical approaches to solve the equations of motion; Crespo da Silva and Glynn [3] as well as Lacarbonara and Yabuno [11] use perturbation methods, while Hamdan and Dado [5] use the two-term harmonic balance method.

The present work instead uses an Eulerian coordinate system to describe the beam deflection, a formulation developed by Novozhilov [12] which has been used previously in the literature [4,8,9]. This coordinate system is intuitive for understanding the beam deflection in both the longitudinal and transverse directions. In addition, the present work does not seek an analytical solution to these equations of motion, but rather employs a fourth-order Runge–Kutta numerical solution method to solve for the beam's motion in time.

Following the work of Tang et al. [9], this work uses the Rayleigh–Ritz method and the Lagrange equations are used to derive the equations of motion. As demonstrated in the work of Dowell and McHugh [4] as well as Luongo and Di Edigio [6], a Lagrange multiplier, *λ*, is added to the Lagrange equations to enforce the inextensibility condition, defined as zero axial strain along the beam [12]. This differs from the method previously used [8,9] in which the constraint equation was directly incorporated into both the potential and kinetic energy terms. This allows for an alternative and perhaps clearer understanding of how the inextensibility affects the system.

The results presented here for the cantilevered beam reiterate results published [3,9,10] which observed that the first bending mode acts nearly linearly. The previously published results also indicate that the cantilever's second bending mode exhibits a softening nonlinearity due to nonlinear inertia [3,10]. Here, we show agreement with this assessment and also extend it to free–free beams. The free–free beam shows some softening at the first bending resonant peak and shows a substantial softening at the third bending resonance.

## Methods

### A Derivation of the Equations of Motion.

Figure 1 illustrates the cantilevered and free–free geometries and the coordinate system used herein.

*f*, respectively, are as follows [4,9]. Note that the potential energy expression is from Tang et al. [9] in which the strain term from Novozhilov's formulation [12] is expanded in a Taylor series in $\u2202w/\u2202x$

The constraint *f* is derived from Novozhilov's [12] strain relations for an inextensible beam. This relation between the longitudinal and the transverse slope ensures that the strain along the beam is always zero, so *f* = 0. This constraint has also already been applied to the potential energy term, *V*, to express it only as a function of the transverse deflection *w* [4,8,9]. Unlike the previous derivations [8,9], however, the effect of the constraint on the kinetic energy term will naturally be incorporated later in the derivation. It may be noted that the constraint *f* is dimensionless.

*i*correlates with

*u*modes,

*j*with

*w*, and

*k*with

*λ*

The modal functions Ψ are chosen to satisfy the geometric constraints, i.e., $u=w=(\u2202w/\u2202x)=0$ for a clamped end and *λ* = 0 for a free end. These modal functions are all normalized such that the length variable *ξ* runs from 0 to 1 and the endpoints are at 1,0, or −1 depending on the geometric constraints.

In the examples in this paper, rigid body translation for *w* is included for the free–free beam, whereas rigid body rotation is excluded.

### Kinetic Energy.

*w*after going through the constraint equation

The rest of the derivatives of kinetic energy are zero, and therefore, only one kinetic energy term is entered into the Lagrange equation for each of the spatial coordinates.

### Potential Energy.

The other derivatives of the potential energy are zero.

### Lagrange Multiplier.

*f*in Eq. (5) is

### System of Equations.

*u*and

_{i}*w*and normalizing the integral terms such that the integration is from 0 to 1 with respect to $x/L\u2261\xi $ produce several useful definitions. Index notation for the one-dimensional vectors and the two-dimensional matrices can be rewritten in boldface matrix notation as indicated after the right arrows. However, the three- and four-dimensional tensors will be left in index form to avoid the introduction of tensor notation

_{j}A discussion on vector/matrix size and characteristics may be appropriate at this point. Vectors **u**, **w**, and $\lambda $ are length *I*, *J*, and *K*, respectively. Matrices $Mu$ and $Mw$ are diagonal matrices of size *I*^{2} and *J*^{2}, respectively. Matrix **A** is of size *I* × *K*. Tensor *B* is three-dimensional of size *K* × *J* × *J*, and *P* is four-dimensional of size *J*^{4}. Matrices $\zeta ,\u2009\omega $, and $\omega 2$ are diagonal matrices of length *J* with each nonzero entry as the corresponding damping coefficient, modal natural frequency, or modal natural frequency squared, respectively.

*J*, and the final term in Eq. (42) is a tensor summation which simplifies to a vector of length

*K*. These may be represented in their final vector form as follows:

### Solution.

To solve this system of equations, differentiate Eq. (45) in time twice and solve for $u\u0308$. Then, substitute $u\u0308$ from Eq. (43) into this new equation to get an expression for *λ _{k}*. Finally, this expression for

*λ*is substituted into Eq. (44) to get one equation of motion for $wj\xaf$. From this, both $u$ and

_{k}**λ**may then be solved once $w$ is known.

*λ*in terms of

*w*

*K*and has units of force. Finally, we substitute this expression for $\lambda $ into Eq. (44) to solve for $w\u0308$, which provides the equation of motion for the transverse deflection

This single equation of motion (49) is marched through time using a fourth-order Runge–Kutta solver, Matlab's ODE45. The physical transverse deflection is then transformed from the modal coefficients using Eq. (7), which allows the longitudinal deflection to be solved from the constraint Eq. (5).

This method easily allows investigations of different boundary conditions. By altering the dimensionless mode shapes outside of the ODE solver, one may compute coefficient matrices for different boundary conditions one time only. This “bank” of dimensionless modal coefficients for different boundary conditions decreases computational time and up-front effort.

To understand the influence of the stiffness nonlinearity versus the inertia nonlinearity, each of these terms is isolated and explored. Four cases are shown for each loading condition: stiffness only, inertia only, fully nonlinear (FNL), and linear. The linear case considers only the first three terms of Eq. (44): linear inertia, linear damping, and linear stiffness. The stiffness only case considers the linear terms plus the fourth term which represents nonlinear stiffness; the inertia only case considers the linear terms plus the final term as nonlinear inertia; and the FNL case considers all terms.

### Physical Model Properties.

To validate our model, static and dynamic comparisons are drawn to Tang et al. [9], and therefore, we use their model parameters. Material properties used are shown in Table 1. Width, thickness, and length are identical to reported values, whereas the modulus *E* and the density *ρ* were chosen such that natural frequencies reported by Tang et al. were matched. For each case, the modal damping is $\zeta =0.01$ for all modes. This quantity, again, was chosen to match published methods and results [9].

## Results

### Cantilevered Beam Results.

For cantilevered static tests, varying mass loads acting under gravity are placed at the tip of the beam, and transverse and longitudinal deflection data are calculated. The *u* deflection is calculated via the constraint equation. Modal convergence occurs after four modes for the static cases, but here five modes are used in each coordinate, as reported by Tang et al. Figure 2 demonstrates agreement between published results—including experimental results—and the current results for the deflection normalized by *L* in the *w* and *u* directions, respectively. The agreement shown here between the current model and the previously published results is encouraging and validating. It is thought that any disagreement between computational and experimental data is due to a difference of the manufacturer's reported versus actual material properties. Future studies will validate this idea by testing material properties in house. It may be noted that the stiffness of the beam is inversely related to the slope of these plots, indicating that the beam stiffens under higher loads.

Comparing the present *dynamic* results to those published [9] is useful in identifying new insights into the cantilevered beam problem. Forward and backward frequency sweeps are conducted to demonstrate nonlinear effects of the model. The forcing amplitude is set at *F* = 0.147 *N* and acts at $xF=0.7*L$ in the *w* direction to match the published configuration [9]. The frequency is swept from 4.5 to 6.5 Hz, capturing the behavior around the linear natural frequency at 5.5 Hz. A convergence study was performed and it was determined that three modes in *w* and six modes in *u* and *λ* are necessary to ensure modal convergence at all time steps. The beam vibration reaches a steady-state after about 35 s, and the final second at each frequency is used to calculate the root-mean-square (RMS) response of each frequency step.

In Fig. 3, these frequency sweeps are shown normalized to the linear natural frequency, Ω. Each case isolates the stiffness and inertia nonlinearities and compares to the full nonlinear case which includes both nonlinearities. For the isolated stiffness or inertia nonlinearity cases, hysteresis occurs. By plotting forward and backward frequency sweeps, both stable branches of the nonlinear response curves are shown. However, the unstable branch which connects the two is not illustrated because the numerical solution always tends to one of these stable branches.

It is interesting to note that by solving these cases at higher frequency resolutions, the fully nonlinear behavior is somewhat different from previously reported by Tang et al. [9]. The present results agree that the stiffness nonlinearity increases the resonant frequency, while the inertia nonlinearity decreases it, as well that the fully nonlinear case shows a balance between the two nonlinear effects. However, the maximum RMS amplitude of the fully nonlinear case is not significantly less than those of either the isolated stiffness or inertia nonlinearity cases, as previously reported. Its steep peak has important implications for experimental verification of such a case, as the frequency sweep must have very high resolution to avoid the risk of passing over the maximum resonant amplitude. Also note that the full nonlinear case shows essentially no hysteresis, and as such, the forward and backward frequency sweeps are essentially identical.

The offsetting contribution of the inertia and stiffness nonlinearities is in agreement with results previously published [3,9,10] which show that the first resonant mode acts nearly linearly.

Figure 4 shows the second resonant response of the cantilever beam. Here, 3 *w* modes and 6 *u* and *λ* modes were used again. Note that the fully nonlinear beam result is significantly softened, even though the stiffness only case shows a strong stiffening effect. Hysteresis is strong and the resonant frequency is shifted by about 6.5% from the linear resonant frequency.

To illustrate the relationship between forcing amplitude and the nonlinear effects, Fig. 5 plots the linear and nonlinear responses for different loads with a backward frequency sweep across the second natural frequency. Even at low loads and small deflections (under 5% of the beam length), the inertia nonlinearity is still seen to shift the resonance frequency downward. At higher loads, as expected, the beam exhibits stronger hysteresis as the resonance frequency continues to decrease.

### Free–Free Results.

The free–free beam model was subjected to similar tests as the cantilevered beam model, but the free–free beam exhibits additional subtleties.

To assess the free–free model, static computations were again performed, using ten modes, and compared to a finite element method (FEM) solution. The beam was held in static rigid body equilibrium by applying an upward force load *F* at *x* = *L*/2 and downward forces *F*/2 at *x* = 0 and *x* = *L*. Total deflection in the *w* direction is computed as the difference in *w* between *x* = *L*/2 and *x* = *L*. As a comparison, the same configuration was generated and computed in the commercial FEM software ANSYS. Figure 6 compares the results of different forces with different solving methods. It illustrates the difference between the nonlinear and linear cases and shows that our model reasonably agrees with the FEM code.

Upon assessment of the model via static computations, frequency sweeps showed the nature of the nonlinear effects on resonance and hysteresis. A sinusoidal force with amplitude *F* = 1.47 *N* was applied in the *w* direction at *x _{F}* =

*L*/2. 4

*w*modes and 9

*u*and

*λ*modes were computed as determined by a convergence study. Notably, since the rigid body translational mode was included, the beam could react to the forces exerted on it by accelerating in the

*w*direction. To inhibit this motion, if the beam's rigid body translation was greater or less than

*L*, a weak uniformly distributed force of ± 0.02

*N*was applied to contain the beam. Deflection is measured as the difference in maximum

*w*from the resting position of the beam after this rigid-body motion has been removed. The beam motion reaches steady-state after about 7 s, and the final 0.1 s is used to calculate the RMS response to each frequency.

Figures 7(a) and 7(b) plot the RMS deflection data for the forward and backward frequency sweeps, respectively, around the first linear natural frequency. Unlike the cantilevered beam, the full nonlinear case for the free–free beam shows more influence from the inertia nonlinearity, and indeed the resonant peak is shifted downward by about 3.3% from the linear natural frequency. As previously noted with the cantilevered beam, the frequency resolution is increased near the resonance peaks for each case.

Figure 8 plots the third resonant RMS deflection data for the forward and backward frequency sweeps, respectively. The third resonant frequency was chosen for study instead of the second because exciting the second would induce rigid body rotation of the beam, which is not yet fully developed in this model. Here, 4 *w* modes and 14 *u* and *λ* modes are computed. Note that more modes in *u* and *λ* are necessary to reach modal convergence for the third mode excitation than for the first mode excitation. Like the first resonant mode response, the fully nonlinear beam is strongly softened, with the resonant frequency shifting by about 8.2% of the linear natural frequency.

Figure 9(a) plots the linear and nonlinear responses for different loads with a backward frequency sweep. As seen in the case of the cantilevered beam, the nonlinear inertia is prevalent in shifting the resonance frequency downward at all loads.

## Conclusions and Future Work

A new formulation and solution method of the equations of motion for inextensible beams has been developed and assessed. In the first resonant mode, it has been shown that the nonlinear cantilevered beam behaves more linearly than previously thought. Hysteresis is not shown, and the nonlinear resonant response is similar to the linear response for a representative case. In the second resonant mode, the cantilevered beam response notably softens due to nonlinear effects at small deflections of 5% of the beam length.

In addition, the free–free beam configuration has been explored, and it has been determined that the nonlinear inertia in this configuration is very important. Hysteresis is clearly apparent in the first and third resonant modes of free–free case, as the nonlinearity substantially decreases the resonance frequency at deflections of only 5% and 1.5% of the beam length, respectively.

The trend toward more dominant nonlinear inertia effects for higher resonant modes is perhaps due to the increased velocity and the decreased amplitude with which the beam vibrates at higher modes.

Future work will investigate other representative cases, and different types of forcing will be included, such as oscillations upon a statically deflected beam. An exploration of the rotational rigid body mode is also needed to better understand the limitations and utility of this formulation. The formulation will be expanded to nonlinear plates and shells, and the work will be supplemented with experimental validation studies.

## Funding Data

Air Force Research Laboratory SMART Scholarship.

## Nomenclature

- EI =
beam bending stiffness

*f*=constraint function

*F*=force applied

*L*=beam length

*m*=mass per unit length of beam

*T*=kinetic energy

*u*=longitudinal deflection

*V*=potential energy

*w*=transverse deflection

*x*=spatial coordinate

*x*=_{F}location of applied force

*λ*=Lagrange multiplier

- $L$ =
Lagrangian

- $\u0307$ =
denotes time derivative