Incorporating Boundary Nonlinearity into Structural Vibration Problems

Abstract

This paper presents a methodology for accurately incorporating the nonlinearity of boundary conditions (BCs) into the mode shapes, natural frequencies, and dynamic behaviour of analytical beam models. Such models have received renewed interest in recent years as a result of their successful implementation in state-of-the-art multiphysics problems. To address the need for this boundary nonlinearity to be more completely captured in the equations of motion, a nonlinear algebra expansion of the classical linear approach for developing solvability conditions for natural frequencies and mode shapes is presented. The method is applicable to any BC that can be accurately represented in polynomial form, either explicitly or through the application of a Taylor expansion; this is the only assumption made in removing the need for the use of analytical approximations of the dynamics themselves. By reducing the BCs of the beam to a system of polynomials, it is possible to utilise the tensor resultant to develop these solvability conditions analogous to the conditions placed on the matrix determinant in linear, classical cases. The approach is first derived for a general set of nonlinear BCs before being applied to two example systems to investigate the importance of including nonlinear tip behaviour in the BCs to accurately predict the system response. In the first, a theoretical, symmetric system, in which a beam is supported by nonlinear springs, is used to explore both the applicability of the methodology and the improvements it can make to the accuracy of the model. Then, the more practical example of a cantilever beam with repulsive magnetic interaction at the tip is used to more explicitly assess the importance of properly incorporating boundary nonlinearity into multiphysics problems.

1. Introduction

The ability to accurately model nonlinear vibrations has become increasingly important as the desire for more efficient and reliable mechanical systems, as well those that interact with external physical fields, has continued to grow. This is largely due to the fact that some nonlinear dynamic behaviour, such as hysteresis and modal interaction, have the potential to be hugely damaging to expensive structures if they are not accurately captured during the design phase. A number of significant system properties relating to dynamic behaviour are influenced by the boundary conditions (BCs). In particular, these conditions have a large impact on the natural frequencies and mode shapes of the underlying linear structure and, hence, the equations of motion. Therefore, it is of the utmost importance that these are correctly captured and incorporated into the equations of motion.

The field of nonlinear dynamics is well-established [1,2,3], with geometrically nonlinear structures demonstrated to exhibit a number of challenging dynamical behaviours. These include, but are not limited to, amplitude-dependent natural frequencies [4,5], bistability of solutions [6,7], interactions between modes [2,8], localisation due to nonlinearity [9], and chaos [10]. This paper addresses the interaction of some of these challenges with the possibility of nonlinear forces defined at the boundary by directly, providing a novel solution that allows such systems to utilise the established and current findings and method from within this varied field. As such, the proposed approach allows systems with nonlinear BCs, such as those in multiphysics problems, to be practically treated in line with those with classical BCs.

For continuous structures, it is common practice to use finite elements to approximate this dynamic behaviour. For complex mechanical systems with multiple subcomponents, it can be beneficial to utilise specialist finite element software to develop an accurate model. However, for simpler structures, such as those at the microscale, it is often possible to develop a thorough understanding of the problem, simply by modelling analytically. There are a number of analytical approaches that can be taken (see [11], for example), but the primary focus of this work is the Ritz–Galerkin method (henceforth referred to simply as the Galerkin method). For the beams considered in this work, this technique assumes that the time- and position-dependent deflection can be discretised as an infinite sum of the product of the projection functions, one of which is solely time-dependent, and the other is solely position-dependent. This approach has been applied to a wide variety of classical beam BCs [12,13] but has more recently been applied to more complicated beam configurations.

Initially, this variation has seen the inclusion of compressive loads [14], as well as the addition of rotational [15] and linear translational [16] springs at an arbitrary point along the beam, though the discussion in these studies is limited to the determination of the natural frequencies and mode shapes of the system, rather than its dynamic behaviour. A similar approach was taken in [17], which considers a beam supported by both a linear translational spring and rotational spring at both ends. In [5], a rotational spring is added to only one end of a pinned-pinned beam, breaking the symmetry of the system and introducing internal resonances between the first and second mode. This highlights the possibility that complex dynamical behaviour can arise in systems with relatively simple BCs.

This discussion has also been expanded to include systems with nonlinear boundary conditions. In [18], a number of such conditions are investigated using Hamilton’s principle, with the discussion including consideration of the free and forced responses of such structures and the variation in frequency and mode shape at higher amplitudes. Perhaps the most significant aspect of this discussion for the present work is that surrounding the variability of the mode shapes depending on the boundary nonlinearity. For two specific beam configurations, characteristic equations that can be solved numerically to find the natural frequencies are developed. Hence, mode shapes of the system, capable of encapsulating the nonlinear tip relationship, can also be created. In the current work, a methodical procedure for finding such an equation for any nonlinear beam BC is developed, providing a nonlinear algebraic parallel to the linear methodology utilised in those systems in the previous paragraph. In doing so, the dependency on an analytical approximation method, such as the harmonic balance method in [18], is removed.

More recently, ref. [19] proposed an analytical method for incorporating the boundary nonlinearity, this time based on a multiple timescale assumption and applying perturbation and modal revision simultaneously. This approach clearly demonstrated the variation in the solutions derived with and without the boundary nonlinearity incorporated, with further validation of an expanded version provided using the finite difference method in [20]. Furthermore, this technique has been applied to an elastic beam with boundary inerter-enhanced nonlinear energy sinks in [21] and expanded to account for pulley support ends in [22]. In recent decades, there has been a marked interest in the more fundamental mathematics of fourth-order, nonlinear boundary value problems, with a series of studies [23,24,25,26] establishing the existence of solutions and generating them through an iterative process, albeit one with a lower rate of convergence that is somewhat prohibitive.

Although there are a number of investigations in which the nonlinearity of the BCs is defined by springs with higher-order stiffnesses [27,28], the methodology proposed in this study can actually be applied to any BC that can be approximated by a Taylor expansion. An example of such a system is presented in [29], who investigate the forced responses of a cantilever beam with magnetic interaction at the tip using a single-mode multiple scale model. Although a hardening behaviour is predicted, consideration is not given to the influence the tip interaction has on the frequencies at which the system vibrates. Recent interest in this area has also included beams with unknown viscoelastic BCs [30], particularly in relation to moving loads, and nonlinear composite beams with moment conditions at the boundary [31]. This continued research in the area reiterates the importance of fully understanding and incorporating the relationship between boundary nonlinearity and nonlinear structural response. Finally, ref. [32] investigated extreme deflections of a cantilever, providing a key reminder of the potential for nonlinearity, even in the absence of the nonlinear BCs discussed here.

It is important to stress that the validity of this method is entirely dependent on the aforementioned ability to represent the true BC using a polynomial expression. If this condition is not met, it is possible that the resultant model will be inaccurate. For example, in the development of this paper, the method was applied to the dynamics of an atomic force microscope (AFMs), approximating the van der Waals force [33] using a third-order Taylor expansion. The approximated equations of motion predicted by the proposed method suggested the system would exhibit a hardening nonlinearity, whereas real AFMs exhibit a softening nonlinearity in the considered region. The examples provided in this paper have been selected based on this ability to represent the BCs using a polynomial expression.

At this point, it is important to highlight the differences and relationship between the boundary nonlinearity—as is the focus of this paper—and the nonlinear terms that arise in the equations of motion. Assuming the widely-applied definition, originally proposed by Rosenberg [34] and extended by Kerschen et al. [35], a nonlinear normal mode (NNM) is defined as a not necessarily synchronous, periodic motion of the system. Note that an alternative, manifold-based definition of NNMs, which extends the concept to damped systems, has been proposed by Shaw and Pierre [36] and expanded by Haller and Ponsioen [37], though this is not the focus of the present study. For the former definition, it is common practice for NNMs to be expressed in terms of the modes of the underlying linear system. Nonlinearity is introduced to the equations of motion through the addition of higher-order strain terms that arise due to the stretching of the structure at larger deformations, often as a result of BCs that restrict the free movement of the system. Increasingly, BCs defined by some nonlinear function of the boundary displacement—i.e., ${\varphi }^n\left(x_{\mathrm{BC}}\right)=F_{\mathrm{NL}}\left(x_{\mathrm{BC}}\right)$, where $x_{\mathrm{BC}}$ is the BC displacement, $F_{\mathrm{NL}}$ is some nonlinear function, and ${\varphi }^n$ is some derivative of the mode shape—have been observed. This is particularly common for multiphysics problems, as discussed below. This work proposes a novel algebraic approach for integrating such nonlinear BCs into analytical vibration models, providing a nonlinear parallel to the treatment of non-classical linear BCs that utilises nonlinear tensor algebra in place of linear algebra. This allows a simple, analytical form of the mode shapes to be found, whilst also allowing the impact of boundary nonlinearity on the natural frequencies to be defined.

2. Analytical Methods

This section outlines the derivation of an analytical model for an Euler–Bernoulli beam with boundary nonlinearity, building on the expansive beam theory literature; an overview of the components of this area that are directly relevant to the present study can be found in [3].

The general equation for the transverse vibration of a beam at small amplitudes can be written as

$${\displaystyle \frac{{\delta }^{2}w}{\delta t^2}}+a^2{\displaystyle \frac{{\delta }^{4}w}{\delta x^4}},\mathrm{where}{a}^2={\displaystyle \frac{EI}{\rho \widehat{A}}},$$

(1)

where $w(x,t)$ is the beam deflection, which is dependent on the position along the beam and time, denoted x and t, respectively; E is the Young’s modulus; I is the second moment of inertia; $\rho$ denotes the density of the beam; and $\widehat{A}$ is its cross-sectional area. The small amplitude assumption is initially included to facilitate the application of the Galerkin method so that a general solution for w can be more readily found.

In applying the Galerkin method, it is assumed that the transverse displacement can be written as a series of functions with distinct spacial and temporal parts. This can be written as

$$w(x,t)=\sum _{j=1}^{\infty }{\varphi }_j\left(x\right)q_j\left(t\right),$$

(2)

where ${\varphi }_j\left(x\right)$ and $q_j\left(t\right)$ denote the mode shapes of the underlying linear system and harmonic functions in time denoting the contribution of the jth mode, respectively.

In this paper, the mode shapes are of particular interest, so the harmonic functions are simply written as

$$q_j\left(t\right)=A\cos \left(p_{j}t\right)+B\sin \left(p_{j}t\right),$$

(3)

where $p_j$ are constants that are defined by the boundary conditions; values for these are found numerically. Applying an arbitrary mode shape, $\varphi \left(x\right)$, and the harmonic form in Equation (3) in Equation (1) allows the general equation to be rewritten as

$${\displaystyle \frac{d^4\varphi }{d{x}^4}}-{\displaystyle \frac{p_j^2}{a^2}}\varphi ={\displaystyle \frac{d^4\varphi }{d{x}^4}}-{\kappa }^4\varphi =0,$$

(4)

where the variable $\kappa$ has been introduced so that $\sin \left(\kappa x\right)$, $\cos \left(\kappa x\right)$, $\sinh \left(\kappa x\right)$, and $\cosh \left(\kappa x\right)$ are particular solutions to Equation (4), with

$$\kappa =\sqrt[4]{{\displaystyle \frac{p_j^2\rho \widehat{A}}{EI}}}.$$

(5)

Now, the general solution to Equation (4) can be written as

$$\varphi \left(y\right)=c_1\cos \left(\kappa x\right)+c_2\sin \left(\kappa x\right)+c_3\cosh \left(\kappa x\right)+c_4\sinh \left(\kappa x\right),$$

(6)

where the $c_k$ are constants defined by the BCs, and $y=x/\ell$ has been introduced to non-dimensionalise the system and remove the dependence on the beam length, ℓ. The specific values of $c_k$ are defined by the boundary conditions of the beam, which has the possibility of being a nonlinear function of displacement in this paper.

If the system is not assumed to have small amplitude vibrations (as in Equation (1)), the length of the beam can cause a geometric nonlinearity to occur. To capture this, it is necessary to consider the tension force that arises in the second derivative of the moment. Full details of the derivation of this additional term are not given here but are readily available in [3], for example. The equations of motion are now given by

$$EI{\displaystyle \frac{{\delta }^{4}w}{\delta x^4}}+\rho \widehat{A}{\displaystyle \frac{{\delta }^{2}w}{\delta t^2}}-{\displaystyle \frac{E\widehat{A}}{2\ell }}{\int }_0^{\ell }{\left({\displaystyle \frac{\delta w}{\delta x}}\right)}^2\mathrm{d}x\left({\displaystyle \frac{{\delta }^{2}w}{\delta x^2}}\right)=0.$$

(7)

By imposing the Galerkin condition (Equation (2)), the equations of motion can be rewritten as

$$EI\sum _{j=1}^{\infty }{\displaystyle \frac{d^4\varphi }{d{x}^4}}{q}_j+\rho \widehat{A}\sum _{j=1}^{\infty }{\varphi }_j{\ddot{q}}_j-{\displaystyle \frac{E\widehat{A}}{2\ell }}\sum _{j=1}^{\infty }{\int }_0^{\ell }\sum _{k=1}^{\infty }{\left({\displaystyle \frac{d{\varphi }_k}{dx}}\right)}^2\mathrm{d}x\left({\displaystyle \frac{{\delta }^2{\varphi }_k}{\delta x^2}}\right)q_j{q}_k^2=0,$$

(8)

where $\ddot{•}$ denotes the second derivative with respect to t. To decouple this equation, we can make use of the orthogonality of the mode shapes. By multiplying Equation (8) by an arbitrary ${\varphi }_n$ and integrating across the length of the beam, the equations of motion can now be expressed as

$$\begin{array}{cc}\hfill EI\sum _{j=1}^{\infty }{\int }_0^{\ell }{\displaystyle \frac{d^4{\varphi }_j}{d{x}^4}}& {\varphi }_n\mathrm{d}x{q}_j+\rho \widehat{A}\sum _{j=1}^{\infty }{\int }_0^{\ell }{\varphi }_j{\varphi }_n\mathrm{d}x{\ddot{q}}_j\hfill \\ \hfill & -{\displaystyle \frac{E\widehat{A}}{2\ell }}\sum _{j=1}^{\infty }{\int }_0^{\ell }\left[\sum _{k=1}^{\infty }\left({\int }_0^{\ell }{\left({\displaystyle \frac{d{\varphi }_k}{dx}}\right)}^2\mathrm{d}x{q}_k^2\right)\left({\displaystyle \frac{d^2{\varphi }_j}{d{x}^2}}\right)q_j\right]\mathrm{d}x=0.\hfill \end{array}$$

(9)

By expanding the various integral terms in this expression (as outlined in more detail in [3]), the nth decoupled equation is given by

$${\ddot{q}}_n+{\displaystyle \frac{EI{\alpha }_{4,n}}{\rho \widehat{A}\ell }}+{\displaystyle \frac{E{\alpha }_{2,n,j}}{2\rho {\ell }^2}}+\sum _{i=1}^N\sum _{j=1}^N\sum _{k=1}^N{\beta }_{i,j}{\alpha }_{2,n,k}{q}_i{q}_j{q}_k=0,$$

(10)

where the $\alpha$ and $\beta$ terms are defined by

$$\begin{array}{cc}\hfill & {\int }_0^{\ell }{\varphi }_n{\varphi }_j=\left\{\begin{array}{cc}\ell ,\hfill & \mathrm{if}n=j\hfill \\ 0,\hfill & \mathrm{otherwise},\hfill \end{array}\right.,\hfill \\ \hfill & {\int }_0^{\ell }{\varphi }_n{\displaystyle \frac{d^2{\varphi }_j}{d{x}^2}}={\alpha }_{2,n,j},\hfill \\ \hfill & {\int }_0^{\ell }{\varphi }_n{\displaystyle \frac{d^4{\varphi }_j}{d{x}^4}}=\left\{\begin{array}{cc}{\alpha }_{4,j},\hfill & \mathrm{if}n=j\hfill \\ 0,\hfill & \mathrm{otherwise},\hfill \end{array}\right.,\hfill \\ \hfill & {\int }_0^{\ell }{\displaystyle \frac{d{\varphi }_i}{dx}}{\displaystyle \frac{d{\varphi }_j}{dx}}={\beta }_{i,j}.\hfill \end{array}$$

(11)

2.1. Linear, Non-Classical Boundary Conditions

Before the inclusion of nonlinear BCs is considered, it is useful to revisit a linear, non-classical beam system, so that direct comparisons may be made when the BCs include nonlinear terms. In particular, we consider the cantilever beam with linear spring support, as investigated in [16]. The left end of the beam is clamped, so the standard BCs are given by

$$\varphi \left(0\right)={\varphi }^{\prime }\left(0\right)=0.$$

(12)

To aid the solution of the normal modes, it is helpful to rewrite Equation (6) in the form

$$\varphi \left(y\right)=c_1\cos \left(\kappa y\right)+c_2\sin \left(\kappa y\right)+(c_3+c_4)\cosh \left(\kappa y\right)+(c_3-c_4)\sinh \left(\kappa y\right),$$

(13)

Then, it can be seen, from the first condition, that $c_1+c_3+c_4=0$ and, from the second condition, that $c_2+c_3-c_4=0$. Solving these equations, the general mode shape can be parametrised—using only two parameters—as

$$\varphi \left(y\right)=c_1(\cos \left(\kappa y\right)-\cosh \left(\kappa y\right)+\sin \left(\kappa y\right)-\sinh \left(\kappa y\right))+c_2(\sin \left(\kappa y\right)-\sinh \left(\kappa y\right)),$$

(14)

For the linear spring support, the BCs are

$${\varphi }^{″}\left(1\right)=0,{\varphi }^{‴}\left(1\right)=f={\widehat{K}}_{\ell }\varphi \left(1\right),\mathrm{where}{\widehat{K}}_{\ell }={\displaystyle \frac{k_{\ell }{\ell }^3}{EI}}.$$

(15)

Here, ${\widehat{K}}_{\ell }$ is the non-dimensionalised linear spring constant. Substituting in the form of $\varphi$ from Equation (14), the non-classical BCs can be written as

$$\begin{array}{cc}\hfill & (\cos \left(\kappa y\right)+\sin \left(\kappa y\right)+\cosh \left(\kappa y\right)+\sinh \left(\kappa y\right))c_1+2(\sin \left(\kappa y\right)+\sinh \left(\kappa y\right))c_2=0,\hfill \\ \hfill & ({\kappa }^3[-\cos \left(\kappa y\right)-\cosh \left(\kappa y\right)-\sin \left(\kappa y\right)-\sinh \left(\kappa y\right)]\hfill \\ \hfill & +{\widehat{K}}_{\ell }[-\cos \left(\kappa y\right)+\cosh \left(\kappa y\right)-\sin \left(\kappa y\right)+\sinh \left(\kappa y\right)])c_1\hfill \\ \hfill & -2({\kappa }^3[\cos \left(\kappa y\right)+\cosh \left(\kappa y\right)]+{\widehat{K}}_{\ell }[\sin \left(\kappa y\right)-\sinh \left(\kappa y\right)])c_2=0\hfill \end{array}$$

(16)

Given that both of these equations are linear expressions in terms $\mathbf{c}={[c_1,c_2]}^T$, it is possible to express this system as

$$\left[A\right]\mathbf{c}=\left[\begin{array}{cc}{a}_{11}& a_{12}\\ a_{21}& a_{22}\end{array}\right]\left[\begin{array}{c}{c}_1\\ c_2\end{array}\right]=\left[\begin{array}{c}0\\ 0\end{array}\right],$$

(17)

where $\left[A\right]$ is defined by the coefficients in Equation (16). The system only has non-trivial solutions if, and only if, its determinant is zero. Therefore, the solvability condition is written as

$$\det \left[A\right]=0.$$

(18)

This expression only has one unknown term, namely, $\kappa$. As such, the solutions of Equation (18) can be used to define the modeshapes in Equation (14).

2.2. Nonlinear Boundary Conditions

In this section, it assumed that both ends of the beam have nonlinear beam supports, so it is necessary to find a nonlinear equivalent to Equation (18). As alluded to in the Introduction, the motivation of this paper is to develop a methodology through which boundary conditions of the form ${\varphi }^n\left(x_{\mathrm{BC}}\right)=F_{\mathrm{NL}}\left(x_{\mathrm{BC}}\right)$ (i.e., containing nonlinear functions of displacement) can be more completely captured in the modal basis of the system and, therefore, the nonlinear strain terms in the equations of motion defining the NNM. Although this derivation considers a generalised structure with nonlinear BCs at both ends, the methodology could easily be simplified to a system in which one BC is nonlinear and the other is linear. As such, this strategy is applicable in a wide variety of multiphysics problems.

The authors note that the complexity of the mathematical strategy and notation has the potential to make this methodology appear more onerous than its implementation. To this end, the following “recipe” for its use is given:

• Express the nonlinear BCs as polynomials, applying a Taylor expansion where necessary.
• Solve the homogeneous BCs, reducing the order of the system to two equations in two unknowns.
• Express the general mode shape in terms of the remaining unknown variables.
• Apply the updated mode shapes in the polynomial BCs.
• Introduce updated polynomial coefficients to simplify the homogeneous equations.
• Assume a general solution for the polynomial equations and introduce a dummy variable to allow them to be expressed in terms of a single variable.
• If necessary, introduce dummy coefficients so that the order of the two polynomials is the same.
• Calculate the tensor resultant of the system, applying Plücker relations, where possible.

These steps are now outlined for a general system with nonlinear BCs. In this case, the BCs are given by

$${\varphi }^{″}\left(0\right)={\varphi }^{″}\left(1\right)=0,{\varphi }^{‴}\left(0\right)=f_0\left(\varphi \left(0\right)\right),{\varphi }^{‴}\left(1\right)=f_{\ell }\left(\varphi \left(1\right)\right),$$

(19)

where $f_0$ and $f_{\ell }$ are nonlinear functions defining the behaviour at the base and tip of the beam, respectively.

• Step 1: Express the nonlinear BCs as polynomials, applying a Taylor expansion where necessary.

The only assumption on the functions defining the BCs are that $f_0,f_{\ell }$ are analytic functions, so that there is coincidence between the functions and their Taylor expansion when the Taylor series is convergent. As such, it can be assumed that the BCs can be approximated as

$${\varphi }^{″}\left(0\right)={\varphi }^{″}\left(1\right)=0,{\varphi }^{‴}\left(0\right)=\sum _{m=0}^M{\zeta }_m\varphi {\left(0\right)}^m,{\varphi }^{‴}\left(1\right)=\sum _{n=0}^N{\eta }_n\varphi {\left(1\right)}^n,$$

(20)

where ${\zeta }_m$ and ${\eta }_n$ are coefficients arising from the respective Taylor expansions of $f_0$ and $f_{\ell }$, and M and N are their respective orders.

• Step 2:  Solve the homogeneous BCs, reducing the order of the system to two equations in two unknowns.

By imposing the general form given in Equation (6), the first two equalities in Equation (20) can be written as

$$\begin{array}{cc}\hfill & {\kappa }^2(c_3-c_1)=0,\hfill \\ \hfill & {\kappa }^2(c_3\cosh \left(\kappa \right)-c_1\cos \left(\kappa \right)+c_4\sinh \left(\kappa \right)-c_2\sin \left(\kappa \right))=0.\hfill \end{array}$$

(21)

These equations can now be solved for $c_3$ and $c_4$:

$$\begin{array}{cc}\hfill & c_3=c_1,\hfill \\ \hfill & c_4={\displaystyle \frac{c_1(\cos \left(\kappa \right)-\cosh \left(\kappa \right))+c_2\sin \left(\kappa \right)}{\sinh \left(\kappa \right)}}.\hfill \end{array}$$

(22)

• Step 3:  Express the general mode shape in terms of the remaining unknown variables.

The definition of these coefficients means that it is now possible to parametrise $\varphi \left(y\right)$ using only $c_1$ and $c_2$:

$$\begin{array}{cc}\hfill & \varphi \left(y\right)=c_1\left[\cos \left(\kappa y\right)+\cosh \left(\kappa y\right)+{\displaystyle \frac{\cos \left(\kappa \right)-\cosh \left(\kappa \right)}{\sinh \left(\kappa \right)}}\sinh \left(\kappa y\right)\right]\hfill \\ \hfill & +c_2\left[\sin \left(\kappa y\right)+{\displaystyle \frac{\sin \left(\kappa \right)}{\sinh \left(\kappa \right)}}\sinh \left(\kappa y\right)\right].\hfill \end{array}$$

(23)

• Step 4:  Apply the updated mode shapes in the polynomial BCs.

By representing the mode shape functions in this way, it is possible to massively simplify the polynomial terms in Equation (20), since the displacements at the boundaries can now be given by

$$\varphi \left(0\right)=2c_1,\varphi \left(1\right)=2(c_1\cos \left(\kappa \right)+c_2\sin \left(\kappa \right)).$$

(24)

Furthermore, it can be noted that the third derivatives at the boundaries can also be more simply expressed as

$$\begin{array}{cc}\hfill & {\varphi }^{‴}\left(0\right)={\kappa }^3\left[c_1{\displaystyle \frac{\cos \left(\kappa \right)-\cosh \left(\kappa \right)}{\sinh \left(\kappa \right)}}+c_2{\displaystyle \frac{\sin \left(\kappa \right)-\sinh \left(\kappa \right)}{\sinh \left(\kappa \right)}}\right],\hfill \\ \hfill & {\varphi }^{‴}\left(1\right)={\kappa }^3[c_1\left({\displaystyle \frac{\cos \left(\kappa \right)-\cos \left(\kappa \right)\cosh \left(\kappa \right)}{\sinh \left(\kappa \right)}}+(\sin \left(\kappa \right)+\sinh \left(\kappa \right))\right)\hfill \\ \hfill & +c_2\left({\displaystyle \frac{\sin \left(\kappa \right)}{tan\left(\kappa \right)}}-\cos \left(\kappa \right)\right)].\hfill \end{array}$$

(25)

Although these terms are slightly more complex, it can be noted that they both take the simple form $A_{\mathrm{BC}}^1{c}_1+A_{\mathrm{BC}}^2{c}_2$ and, hence, can only be linear contributions to the polynomial BCs.

At this point, we now know the specific form of the two unsolved BCs. To aid their solution, it is useful to express these in a more generalised form. As such, the BCs can now be written as

$$\begin{array}{cc}\hfill & x=0:A_0^1{c}_1+A_0^2{c}_2=\sum _{m=0}^M\left[\sum _{k_1+k_2=m}{\mu }_0^{k_1{k}_2}{c}_1^{k_1}{c}_2^{k_2}\right]\hfill \\ \hfill & x=\ell :A_{\ell }^1{c}_1+A_{\ell }^2{c}_2=\sum _{n=0}^N\left[\sum _{k_1+k_2=n}{\mu }_{\ell }^{k_1{k}_2}{c}_1^{k_1}{c}_2^{k_2}\right].\hfill \end{array}$$

(26)

It can be noted that the two BCs are of identical form. Therefore, for the sake of brevity, only a general BC is defined from this point forward, with K denoting the order of equation $f_{B}C$.

• Step 5:  Introduce updated polynomial coefficients to simplify the homogeneous equations.

The right hand side of both equations in Equation (26) gives a more explicit definition of the polynomial BCs in Equation (20). These expressions correspond to the equivalent linear conditions given in Equation (16) but are now given in the form of polynomials with tensor coefficients. Similarly to the linear case, it is useful to rewrite these equalities in the form $f(c_1,c_2)=0$. This is easily carried out for Equation (26), which can now be expressed as

$$\sum _{k=0}^K\left[\sum _{k_1+k_2=k}{\widehat{\mu }}_{\mathrm{BC}}^{k_1{k}_2}{c}_1^{k_1}{c}_2^{k_2}\right]=0,$$

(27)

where

$${\widehat{\mu }}_{\mathrm{BC}}^{k_1{k}_2}=\left\{\begin{array}{cc}{\mu }_{\mathrm{BC}}^{k_1{k}_2}-A_{\mathrm{BC}}^j,\hfill & \mathrm{if}{k}_1+k_2=1\hfill \\ {\mu }_{\mathrm{BC}}^{k_1{k}_2},\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(28)

Now that this system has been defined, it is possible to expand the linear methodology of the previous section using nonlinear algebraic methods, allowing a solution to be found. Further reading and explanation on the steps in this section are thoroughly explained in [38].

For the previously discussed linear system, it was possible to express the equations in the form $\left[A\right]\mathbf{c}=0$, then find the solvability condition by solving $\det \left[A\right]=0$. Given that Equation (27) is now in non-homogeneous tensor form, it is necessary to apply the generalised Craemer rule [38], which allows a solvability condition to be developed by setting the tensor resultant to zero.

• Step 6:  Assume a general solution for the polynomial equations and introduce a universal variable to allow them to be expressed in terms of a single term.

The initial step is to assume that there is some known solution to the system of equations, denoted $(X,Y)$. Here, this would effectively be given by $(X,Y)=(c_1,c_2)$, though this is not immediately implemented to avoid overcomplicating the discussion.

By assuming that this solution is known, it is further possible to rewrite the system in terms of two new variable pairs, either $(z,y)$ or $(x,z)$. The variable z is a dummy variable that allows the system to be written using only one of the solutions from $(X,Y)$. More explicitly, Equation (27) can now be written either as

$$\sum _{k=0}^K\left[\sum _{k_1+k_2=k}{\widehat{\mu }}_{\mathrm{BC}}^{k_1{k}_2}{X}^{k_1}{z}^{k_1}{y}^{k_2}\right]=0,$$

(29)

or

$$\sum _{k=0}^K\left[\sum _{k_1+k_2=k}{\widehat{\mu }}_{\mathrm{BC}}^{k_1{k}_2}{Y}^{k_2}{x}^{k_1}{z}^{k_2}\right]=0.$$

(30)

Now, Equations (29) and (30) are solved by $(z,y)=(1,Y)$ and $(x,z)=(X,1)$, respectively. It can be noted that these systems must necessarily lead to the same overall solution, thus it is only necessary to consider one from this point forward; Equation (29) is chosen here.

• Step 7:  If necessary, introduce dummy coefficients, so that the order of the two polynomials is the same.

Once more, it is possible to make a minor amendment to Equation (27) to aid its solution. It is important to recall that it is not guaranteed that $M=N$. Thus, it is useful to introduce the parameter $\theta =\max \{M,N\}$ and write

$$\sum _{k=0}^{\theta }\left[\sum _{j=0}^k{\widehat{\mu }}_{\mathrm{BC}}^{(m-j),j}{X}^{(m-j)}\right]z^{(\theta -k)}{y}^k=\sum _{k=0}^{\theta }{\nu }_{\mathrm{BC}}^{(\theta -k),k}{z}^{(\theta -k)}{y}^k=0.$$

(31)

Here, we have forced the order of both ${\widehat{\mu }}_{\mathrm{BC}}$ to be the same by expanding the lower-order set with dummy coefficients. Additionally, the assumed constant solution, X, is incorporated into the new tensor coefficients, ${\nu }_{\mathrm{BC}}$, which allows the simplification of the polynomials in z and y. It is worth noting that these coefficients are now defined in terms of $c_1$ and $\kappa$; i.e., ${\nu }_{\mathrm{BC}}={\nu }_{\mathrm{BC}}(c_1,\kappa )$.

• Step 8:  Calculate the tensor resultant of the system, applying Plücker relations, where possible.

As previously mentioned, the linear equivalent of such a system of equations can be solved by ensuring that the determinant is equal to zero. In the nonlinear case, this is mirrored by setting the resultant—a higher-order generalisation of the determinant—to zero. A formal definition of the resultant is given in [38], though is not required here as this work makes use of the fact that the system defined by Equation (31) comprises only two equations in two variables, even though the polynomial order of each equation is $\theta$. When this is the case, one simply needs to find the ordinary resultant [39] of two polynomials in a single variable. This specific case can actually be solved using linear algebra in the following fashion [38]

$$\begin{array}{cc}\hfill & {\mathcal{R}}_{2|\theta }\left(\nu \right)=\hfill \\ & \underset{2\theta \times 2\theta }{\det }\left[\begin{array}{cccccccccc}{\nu }_0^{0,\theta }& {\nu }_0^{1,\theta -1}& {\nu }_0^{2,\theta -2}& \dots & {\nu }_0^{\theta -1,1}& {\nu }_0^{\theta ,0}& 0& 0& \dots & 0\\ 0& {\nu }_0^{1,\theta -1}& {\nu }_0^{2,\theta -2}& \dots & {\nu }_0^{\theta -2,2}& {\nu }_0^{\theta -1,1}& {\nu }_0^{\theta ,0}& 0& \dots & 0\\ & & & & \dots & & & & & \\ 0& 0& 0& \dots & {\nu }_0^{0,\theta }& {\nu }_0^{1,\theta -1}& {\nu }_0^{2,\theta -2}& {\nu }_0^{3,\theta -3}& \dots & {\nu }_0^{\theta ,0}\\ {\nu }_{\ell }^{0,\theta }& {\nu }_{\ell }^{1,\theta -1}& {\nu }_{\ell }^{2,\theta -2}& \dots & {\nu }_{\ell }^{\theta -1,1}& {\nu }_{\ell }^{\theta ,0}& 0& 0& \dots & 0\\ 0& {\nu }_{\ell }^{1,\theta -1}& {\nu }_{\ell }^{2,\theta -2}& \dots & {\nu }_{\ell }^{\theta -2,2}& {\nu }_{\ell }^{\theta -1,1}& {\nu }_{\ell }^{\theta ,0}& 0& \dots & 0\\ & & & & \dots & & & & & \\ 0& 0& 0& \dots & {\nu }_{\ell }^{0,\theta }& {\nu }_{\ell }^{1,\theta -1}& {\nu }_{\ell }^{2,\theta -2}& {\nu }_{\ell }^{3,\theta -3}& \dots & {\nu }_{\ell }^{\theta ,0}\end{array}\right]=0,\hfill \end{array}$$

(32)

where ${\mathcal{R}}_{2|\theta }\left(\nu \right)$ is the tensor resultant of a system, $\nu$, of 2 polynomials of order $\theta$. It can now be recalled that the coefficient tensor ${\nu }_{B}C$ is defined in terms of $c_1$ and $\kappa$. It is possible to gain an understanding of the mode shapes by setting $c_1=1$, since the shape will be scaled as the system vibrates. In doing so, Equation (32) provides a solvability condition that can be numerically solved to define $\kappa$. Although a numerical solution to Equation (32) is possible, as $\theta$ grows, it becomes increasingly difficult to ensure the accuracy of the solution to a sufficient number of decimal places. As such, it can be useful to express this large determinant in terms of elementary, $2\times 2$ determinants, referred to as Plücker relations; this is explored in the following examples and thoroughly outlined in [38].

3. Example: Symmetric, Nonlinear Spring-Supported Beam

To demonstrate the insight that can be found via this method, as well as the practicalities of its implementation, the beam configuration displayed in Figure 1 is introduced. In this model, the beam is supported by identical springs with linear, quadratic, and cubic stiffnesses, given by ${\widehat{K}}_1$, ${\widehat{K}}_2$, and ${\widehat{K}}_3$, respectively.

The BCs for this beam are given by

$$\begin{array}{cc}\hfill & {\varphi }^{″}\left(0\right)={\varphi }^{″}\left(1\right)=0,{\varphi }^{‴}\left(0\right)={\widehat{K}}_1^{\left(0\right)}\varphi \left(0\right)+{\widehat{K}}_2^{\left(0\right)}\varphi {\left(0\right)}^2+{\widehat{K}}_3^{\left(0\right)}\varphi {\left(0\right)}^3,\hfill \\ \hfill & {\varphi }^{‴}\left(1\right)={\widehat{K}}_1^{(\ell )}\varphi \left(1\right)+{\widehat{K}}_2^{(\ell )}\varphi {\left(1\right)}^2+{\widehat{K}}_3^{(\ell )}\varphi {\left(1\right)}^3,\hfill \end{array}$$

(33)

where ${\widehat{K}}_n^{\left(x\right)}$ defines the nth-order spring stiffness at x.

By addressing the homogeneous, second-order BCs, it is possible to rewrite the mode shape equation in Equation (6) as

$$\begin{array}{cc}\hfill \varphi \left(y\right)=& c_1\left(\cos \left(\kappa y\right)+\cosh \left(\kappa y\right)+\left[{\displaystyle \frac{\cos \left(\kappa \right)}{\sinh \left(\kappa \right)}}-\coth \left(\kappa \right)\right]\sinh \left(\kappa y\right)\right)\hfill \\ \hfill & +c_2\left(\sin \left(\kappa y\right)+{\displaystyle \frac{\sin \left(\kappa \right)}{\sinh \left(\kappa \right)}}\sinh \left(\kappa y\right)\right).\hfill \end{array}$$

(34)

Since the BCs defined in Equation (33) take the exact form given in Equation (20), it is not necessary to repeat the detailed algebraic calculations. Instead, it is possible to immediately consider the system as presented in Equation (32), where the coefficients in ${\nu }_0$ and ${\nu }_{\ell }$ are complicated expressions in terms of ${\widehat{K}}_n^{\left(x\right)}$ and $\kappa$. These are explicitly given by

$$\begin{array}{cc}\hfill & {\nu }_0^{0,3}=\left[2{\widehat{K}}_1^{\left(0\right)}+{\kappa }^3\left(\coth \left(\kappa \right)-{\displaystyle \frac{\cos \left(\kappa \right)}{\sinh \left(\kappa \right)}}\right)\right]X+4{\widehat{K}}_2^{\left(0\right)}{X}^2+8{\widehat{K}}_3^{\left(0\right)}{X}^3,\hfill \\ \hfill & {\nu }_0^{1,2}={\kappa }^3\left(1-{\displaystyle \frac{\sin \left(\kappa \right)}{\sinh \left(\kappa \right)}}\right),\hfill \\ \hfill & {\nu }_0^{2,1}=0,\hfill \\ \hfill & {\nu }_0^{3,0}=0,\hfill \\ \hfill & {\nu }_{\ell }^{0,3}=(2{\widehat{K}}_1^{(\ell )}\cos \left(\kappa \right)+{\kappa }^3[(\cosh \left(\kappa \right)-\cos \left(\kappa \right))\coth \left(\kappa \right)-\sin \left(\kappa \right)-\sinh \left(\kappa \right)])X\hfill \\ & +4{\widehat{K}}_2^{(\ell )}{\cos }^2\left(\kappa \right)X^2+8{\widehat{K}}_3^{(\ell )}{\cos }^3\left(\kappa \right)X^3,\hfill \\ \hfill & {\nu }_{\ell }^{1,2}={\kappa }^3(\cos \left(\kappa \right)-\sin \left(\kappa \right)\coth \left(\kappa \right))+2{\widehat{K}}_1^{(\ell )}\sin \left(\kappa \right)\hfill \\ & +8{\widehat{K}}_2^{(\ell )}\sin \left(\kappa \right)\cos \left(\kappa \right)X+24{\widehat{K}}_3^{(\ell )}\sin \left(\kappa \right){\cos }^2\left(\kappa \right)X^2,\hfill \\ \hfill & {\nu }_{\ell }^{2,1}=4{\widehat{K}}_2^{(\ell )}{\sin }^2\left(\kappa \right)+244{\widehat{K}}_3^{(\ell )}{\sin }^2\left(\kappa \right)\cos \left(\kappa \right)X,\hfill \\ \hfill & {\nu }_{\ell }^{3,0}=8{\widehat{K}}_3^{(\ell )}{\sin }^3\left(\kappa \right).\hfill \end{array}$$

(35)

To simplify the investigation into the varying behaviour for this symmetric beam, it is beneficial to use a single parameter, $\gamma$, to define the stiffnesses as

$${\widehat{K}}_1=\gamma ,{\widehat{K}}_2=10\gamma ,{\widehat{K}}_3=100\gamma .$$

(36)

Figure 2, Figure 3 and Figure 4 display the first three mode shapes for $\gamma \in$$[1\times {10}^{-9},1\times {10}^{-6},1\times {10}^{-3},1]$. It can be immediately noted that although the system itself is symmetric, this is distinctly not true for the mode shapes presented in these figures. This behaviour is most pronounced at lower values of $\gamma$. In particular, for $\gamma \le$$1\times {10}^{-6}$, the first mode shapes effectively have a curvature of zero. When the spring stiffness is this low, it is greatly surpassed by the beam stiffness. As such, the first mode shape is equivalent to that of a rigid beam supported by nonlinear springs.

The mode shapes presented in Figure 2, Figure 3 and Figure 4 have been validated using Ansys Mechanical. The system was modelled as a 0.5 m × 0.01 × 0.001 m beam made of structural steel, with a Young’s Modulus of E = 200 GPa and density of $\rho =7850$ kg/${\mathrm{m}}^3$. The beam was modelled using 1D beam elements (BEAM188), each of length 0.01 m, and the nonlinear springs at the boundary were modelled using COMBIN39 elements with user-defined nonlinear stiffness. As can be seen in

Table 1. Natural frequencies (Hz) of the approximations of the cantilever model.

Approximation	${\mathit{\omega }}_{\mathit{n},1}$	${\mathit{\omega }}_{\mathit{n},2}$	${\mathit{\omega }}_{\mathit{n},3}$
Ansys	6.71	21.46	57.57
Free	3.26	20.44	57.23
Linear	3.59	20.49	57.25
Cubic	6.88	21.24	57.48

, the agreement between the Ansys model and predictions from the current work was shown to be good, with a maximum error of 5.26% for the first mode. This higher error is likely to be due to the fact that the user-definition of the nonlinear spring only allows for a small number of discrete data points and the dominance of the spring stiffness in the first mode. The remaining natural frequencies were shown to have errors of less than 1%, with most being significantly below this level. Further validation of the proposed approach has been provided through the dynamic analysis of a full model of the beam created using Ansys. The true behaviour was observed to be close to the true behavior, as presented in Figure 2, which presents a 0–2 Hz sine sweep. However, the way in which nonlinear springs are defined is seen to lead to some minor disagreements.

As $\gamma$ is increased, the mode shapes of the system approach those of the pinned-pinned beam, due to the fact that this increased stiffness limits the extent to which the beam is able to move but without restricting its rotation. That being said, regardless of the value of $\gamma$, the boundaries always have a non-zero displacement. Naturally, as $\gamma \to \infty$, this difference approaches zero. As such, the pinned-pinned beam represents a limit case for the symmetric configuration.

As demonstrated in Equation (10), the nonlinear cross-coupling terms in the ${\alpha }_2$ and $\beta$ matrices are directly derived from the mode shapes. As such, the nature of the curves in Figure 2, Figure 3 and Figure 4 have a large influence on the nonlinear behaviour of the beam. Figure 5 displays the backbone curves for $\gamma \in$$[1\times {10}^{-9},1\times {10}^{-6},1\times {10}^{-3}]$, and it can be seen that as $\gamma$ increases, so does the divergence of the backbone curve from the linear approximation. This is due to the growth in the amount of stretching induced by the BCs, increasing the magnitude of the geometric nonlinearity. Therefore, it is demonstrated that by incorporating this nonlinearity into the BCs, it is possible to predict a wide variety of internal resonance behaviour in beams. While the case $\gamma =$$1\times {10}^{-9}$ presents a backbone curve that is effectively linear, the backbone curves for the cases $\gamma =$$1\times {10}^{-6}$ and $\gamma =$$1\times {10}^{-3}$ present much more complex dynamics. In both cases, there is an inherent hardening effect, in line with the additional stretching that a stiffer spring would facilitate. However, in both cases, there is a specific frequency past which the displacement of the first mode begins to decrease again, indicating that energy has been transferred to a higher-order mode. When $\gamma =$$1\times {10}^{-3}$, the behaviour becomes more complicated still, which is consistent with the complexity of the mode shapes and the apparent interaction between the beam and spring stiffnesses. Further, this highlights the importance of including such nonlinearities in the BCs themselves, as opposed to simply projecting onto a basis mode shapes derived based solely on the linear system. In the following section, this discussion is expanded by considering the use of these polynomial BCs as part of a Taylor expansion for a more complicated boundary interaction, as opposed to simply as a nonlinear spring.

4. Example: Cantilever Beam with Magnetic Tip Interaction

In this section, a more practical example of a cantilever beam, its tip placed equidistantly between two identical magnets, is considered. The magnetic interaction is assumed to occur at a single point at the tip of the beam, and its mass is assumed to be negligible. Note that the clamped end of the beam has a classical boundary condition, so dummy variables are used, as discussed in Section 2.

The repulsive force at the beam tip can be written as

$$F_{\mathrm{tip}}\left(w\right)={\displaystyle \frac{C}{{(D-w)}^2}}+{\displaystyle \frac{C}{{(D+w)}^2}},$$

(37)

where D is the distance between the tip of the undeformed beam and the metal plate below, and C is a constant. As this study is concerned with understanding the mode shapes, and representative responses of the beam itself, D is assumed to be constant. A schematic for this system is presented in Figure 6, in which the tip interaction is found.

Since Equation (37) is not of polynomial form, it is necessary to approximate this using a Taylor expansion about the equilibrium position, as outlined in Step 1. The most appropriate truncation of this expression may not be known a priori or may be determined experimentally. As such, Equation (37) is approximated as both a linear and cubic expression so that the influence of this decision may be assessed. The two truncated expressions are given by

$$\begin{array}{cc}\hfill F_{\mathrm{tip},1}\left(w\right)& ={\displaystyle \frac{2C}{D^3}}w\hfill \\ \hfill F_{\mathrm{tip},3}\left(w\right)& ={\displaystyle \frac{4C}{D^3}}w+{\displaystyle \frac{8C}{D^5}}{w}^3\hfill \end{array}$$

(38)

Now, the complete set of BCs is given by

$$\varphi \left(0\right)={\varphi }^{\prime }\left(0\right)={\varphi }^{\prime }\left(1\right)=0,{\varphi }^{‴}\left(1\right)=F_{\mathrm{tip},n}\left(\varphi \left(1\right)\right).$$

(39)

where $F_{\mathrm{tip},n}$ may represent either truncation.

The procedure defined in Section 2.2 can make use of the consideration given to the classical BC in the linear case. This leads to the mode shape being parametrised as

$$\varphi \left(y\right)=c_1(\cos \left(\kappa y\right)-\cosh \left(\kappa y\right))+c_2(\sin \left(\kappa y\right)-\sinh \left(\kappa y\right)),$$

(40)

Implementing this expression, the BCs at $x=\ell$ can be written as

$$\begin{array}{cc}\hfill & (\cos \left(\kappa \right)+\cosh \left(\kappa \right))c_1+(\sin \left(\kappa \right)+\sinh \left(\kappa \right))c_2=0,\hfill \\ \hfill & {\displaystyle \frac{4C}{D^3}}\varphi \left(1\right)+\left[{\displaystyle \frac{8C}{D^5}}\varphi {\left(1\right)}^3\right]={\kappa }^3((\sin \left(\kappa \right)-\sinh \left(\kappa \right))c_1+(\cos \left(\kappa \right)+\cosh \left(\kappa \right))c_2).\hfill \end{array}$$

(41)

The cubic term is given in square brackets to highlight the fact that a linear or cubic truncation can be used. From this point forward, only the cubic case will be displayed. Expanding this expression, it is possible to collect the $c_n$ terms so that the system can be written as

$$\begin{array}{cc}\hfill 0& ={\sigma }_0^1{c}_1+{\sigma }_0^2{c}_2,\hfill \\ \hfill \sum _{k_1+k_2=3}{\mu }_{\ell }^{k_1{k}_2}{c}_1^{k_1}{c}_2^{k_2}+\sum _{k_1+k_2=2}{\mu }_{\ell }^{k_1{k}_2}{c}_1^{k_1}{c}_2^{k_2}+{\mu }_{\ell }^1{c}_1+{\mu }_{\ell }^2{c}_2+{\mu }_{\ell }^0& ={\sigma }_{\ell }^1{c}_1+{\sigma }_{\ell }^2{c}_2,\hfill \end{array}$$

(42)

where ${\sigma }_B{C}^i$ and $\mu B{C}^{ijk}$ are the collected coefficients from Equation (41). Although it is possible to directly find the resultant of this tensor system, the process is aided by introducing dummy terms to the first of these equations, as in Step 7. Then, the system can be written as

$$\begin{array}{cc}\hfill \sum _{k_1+k_2=3}{\mu }_0^{k_1{k}_2}{c}_1^{k_1}{c}_2^{k_2}+\sum _{k_1+k_2=2}{\mu }_0^{k_1{k}_2}{c}_1^{k_1}{c}_2^{k_2}+{\mu }_0^1{c}_1+{\mu }_0^2{c}_2+{\mu }_0^0& ={\sigma }_0^1{c}_1+{\sigma }_0^2{c}_2,\hfill \\ \hfill \sum _{k_1+k_2=3}{\mu }_{\ell }^{k_1{k}_2}{c}_1^{k_1}{c}_2^{k_2}+\sum _{k_1+k_2=2}{\mu }_{\ell }^{k_1{k}_2}{c}_1^{k_1}{c}_2^{k_2}+{\mu }_{\ell }^1{c}_1+{\mu }_{\ell }^2{c}_2+{\mu }_{\ell }^0& ={\sigma }_{\ell }^1{c}_1+{\sigma }_{\ell }^2{c}_2,\hfill \end{array}$$

(43)

where ${\mu }_0^{ijk}$ are dummy coefficients that allow Equation (42) to be written in the form given in Equation (26). By introducing the change in notation given in Equations (27) and (28), Equation (43) can be written as

$$\begin{array}{cc}\hfill \sum _{k_1+k_2=3}{\widehat{\mu }}_0^{k_1{k}_2}{c}_1^{k_1}{c}_2^{k_2}+\sum _{k_1+k_2=2}{\widehat{\mu }}_0^{k_1{k}_2}{c}_1^{k_1}{c}_2^{k_2}+{\widehat{\mu }}_0^1{c}_1+{\widehat{\mu }}_0^2{c}_2+{\widehat{\mu }}_0^0& =0,\hfill \\ \hfill \sum _{k_1+k_2=3}{\widehat{\mu }}_{\ell }^{k_1{k}_2}{c}_1^{k_1}{c}_2^{k_2}+\sum _{k_1+k_2=2}{\widehat{\mu }}_{\ell }^{k_1{k}_2}{c}_1^{k_1}{c}_2^{k_2}+{\widehat{\mu }}_{\ell }^1{c}_1+{\widehat{\mu }}_{\ell }^2{c}_2+{\widehat{\mu }}_{\ell }^0& =0,\hfill \end{array}$$

(44)

It has been seen that by assuming that $(c_1,c_2)=(X,Y)$ is a solution to this system, it is possible to express it in the following homogeneous form

$$\begin{array}{cc}\hfill Z_3^0{z}^3+Z_2^0{z}^{2}y+Z_1^{0}z{y}^2+Z_0^0{y}^3& =0,\hfill \\ \hfill Z_3^{\ell }{z}^3+Z_2^{\ell }{z}^{2}y+Z_1^{\ell }z{y}^2+Z_0^{\ell }{y}^3& =0,\hfill \end{array}$$

(45)

where

$$\begin{array}{cc}\hfill Z_3^{\mathrm{BC}}& ={\widehat{\mu }}_{\mathrm{BC}}^{30}{X}^3+{\widehat{\mu }}_{\mathrm{BC}}^{20}{X}^2+{\widehat{\mu }}_{\mathrm{BC}}^{10}X+{\widehat{\mu }}_{\mathrm{BC}}^{00},\hfill \\ \hfill Z_2^{\mathrm{BC}}& ={\widehat{\mu }}_{\mathrm{BC}}^{21}{X}^2+{\widehat{\mu }}_{\mathrm{BC}}^{11}X+{\widehat{\mu }}_{\mathrm{BC}}^{01},\hfill \\ \hfill Z_1^{\mathrm{BC}}& ={\widehat{\mu }}_{\mathrm{BC}}^{12}X+{\widehat{\mu }}_{\mathrm{BC}}^{02},\hfill \\ \hfill Z_0^{\mathrm{BC}}& ={\widehat{\mu }}_{\mathrm{BC}}^{03}.\hfill \end{array}$$

(46)

This system is solved by $(z,y)=(1,Y)$. As such, the solvability condition for this system is found by calculating

$${\mathcal{R}}_{2|3}\left(\nu \right)=\det \left[\begin{array}{cccccc}{Z}_3^0& Z_2^0& Z_1^0& Z_0^0& 0& 0\\ 0& Z_3^0& Z_2^0& Z_1^0& Z_0^0& 0\\ 0& 0& Z_3^0& Z_2^0& Z_1^0& Z_0^0\\ Z_3^{\ell }& Z_2^{\ell }& Z_1^{\ell }& Z_0^{\ell }& 0& 0\\ 0& Z_3^{\ell }& Z_2^{\ell }& Z_1^{\ell }& Z_0^{\ell }& 0\\ 0& 0& Z_3^{\ell }& Z_2^{\ell }& Z_1^{\ell }& Z_0^{\ell }\end{array}\right]=0.$$

(47)

Although this form for the resultant may appear somewhat complicated, Equation (47) can be more simply solved through the use of algebraic software. As previously discussed, it is useful to utilise the Plücker relations to simplify this expression. For a system of two cubic polynomials, this is given by [38]

$${\mathcal{R}}_{2|3}\left(\nu \right)=U_3^3-U_2{U}_3{U}_4+U_2^2{U}_5-2U_1{U}_3{U}_5-U_1{V}_3{U}_5,$$

(48)

where

$$\begin{array}{c}\hfill U_1=∥\begin{array}{cc}{Z}_1^0& Z_0^0\\ Z_1^{\ell }& Z_0^{\ell }\end{array}∥,U_2=∥\begin{array}{cc}{Z}_2^0& Z_0^0\\ Z_2^{\ell }& Z_0^{\ell }\end{array}∥,\\ \hfill U_3=∥\begin{array}{cc}{Z}_3^0& Z_0^0\\ Z_3^{\ell }& Z_0^{\ell }\end{array}∥,V_3=∥\begin{array}{cc}{Z}_2^0& Z_1^0\\ Z_2^{\ell }& Z_1^{\ell }\end{array}∥,\\ \hfill U_4=∥\begin{array}{cc}{Z}_3^0& Z_1^0\\ Z_3^{\ell }& Z_1^{\ell }\end{array}∥,U_5=∥\begin{array}{cc}{Z}_3^0& Z_2^0\\ Z_3^{\ell }& Z_2^{\ell }\end{array}∥,\end{array}$$

(49)

Note that the corresponding resultant for the quadratic case would be given by

$${\mathcal{R}}_{2|3}\left(\nu \right)=U_4^2-V_3{U}_5.$$

(50)

First, it is possible to directly investigate the modal characteristics of the various approximations by looking at the mode shapes and natural frequencies. Figure 7 presents the first three mode shapes for the beam in question, showing the linear and cubic approximations in comparison with the free cantilever and pinned cantilever mode shapes. Note that the first of these is often utilised to form the modal basis, whereas the second is included as a reference comparison. The corresponding natural frequencies are compared with the simulated values (again found using Ansys) in Table 1. In Figure 8, the cubic-order model is validated via a linear sine sweep form 0–15 Hz. The predicted solution can be observed to be very close to true solutions generated using Ansys, with some minor deviation caused by the minor difference in ${\omega }_{n,1}$ values.

In both Figure 7 and Table 1, it is apparent that there is a significant difference between the models produced using the linear and cubic truncations. This highlights the importance of testing the accuracy of any Taylor expansion used prior to modelling the system dynamics. In particular, the linear model is close to the free cantilever case, whereas the cubic model is closer to the pinned cantilever model. Furthermore, it can be observed that the linear model natural frequencies become closer to those of the free cantilever as the mode number increases, but there is no similar trend for the linear case. Across all three models, the cubic expansion of the magnetic iteration provides much more accurate natural frequencies, particularly for ${\omega }_{n,1}$.

The effect that these changes have on the system dynamics can be observed in Figure 9, in which the backbone curves for each approximation are displayed. For the linear case, not only is the base of the backbone curve (equal to ${\omega }_{n,1}$) close to that of the free cantilever and not the cubic approximation, it also exhibits a complex softening-to-hardening behaviour that is not present in the more complete cubic approximation, and expected for a cantilever beam with repulsive magnetic interaction at the tip [29]. This model shows that the behaviour should, in fact, be hardening. Interestingly, both models predict a modal interaction with the second mode, though the information provided by the linear model may, again, be of limited use to the user. The discussion of this paragraph demonstrates that the nonlinear BC interaction must be accurately captured to prevent incorrect predictions being made.

5. Conclusions and Discussion

This paper presents a novel method for the incorporation of nonlinear BCs into the modal properties and equations of motion for beams. The proposed technique can be considered as an expansion of the traditional linear methodology, in which the solvability condition arises as a determinant in terms of the coefficients of the general mode shape. Naturally, by allowing the BCs to be nonlinear, this expanded method requires the use of nonlinear algebra, with the determinant equation now replaced by a more general tensor resultant. The technique only assumes that the boundary nonlinearity can be expressed in polynomial form, either directly or through the use of a Taylor expansion.

Two example systems have been used to investigate the applicability of this method. The first of these is a symmetric beam configuration in which both tips are supported by nonlinear springs with linear, quadratic, and cubic stiffness components. The second case is a cantilever beam with a repulsive magnetic interaction at the tip. Through these examples, the importance of accurately including the boundary nonlinearity to the modal characteristics is made clear. The symmetric beam showed increasingly complex backbone curves as the spring stiffness was increased, while both the backbone curves and underlying linear natural frequencies of the cantilever model showed significant error when the BC was not included or incorrectly approximated. As such, the method proposed in this paper provides an innovative method for properly capturing the impact of analytic boundary nonlinearity in dynamical models. In future work, this methodology will be expanded to include more complex geometries and to accommodate non-analytic boundary nonlinearity.