Frequency response of a mass grounded by linear and nonlinear springs in series: An exact analysis

Abstract

An exact analytical solution for the frequency response of a system consisting of a mass grounded by linear and nonlinear springs in series was derived. The system is applicable in the study of beams carrying intermediate rigid mass. The exact solution was derived by naturally transforming the integral of the governing nonlinear ODE into a form that can be expressed in terms of elliptic integrals. Hence, the exact frequency–amplitude solution was derived in terms of the complete elliptic integral of the first and third kinds. Periodic solutions were found to exist for all real values of $\epsilon$ and for $\eta >0$ . On the other hand, linear or weakly nonlinear frequency–amplitude responses were found to occur during small-amplitude vibrations ( $A\le 0.1$ ), very large-amplitude vibrations with strong hardening nonlinearity ( $A\ge 10$ and $\epsilon \gg 1.0$ ), and for all amplitudes when $\eta \cong 1$ . Simulations showed that the system’s periodic response is significantly influenced by the nonlinearity parameter ( $\epsilon$ ), linearity parameter ( $\eta$ ), and amplitude of vibration ( $A$ ). Besides, the existence of bifurcation points at $\epsilon =0$ and different values of $\eta$ was confirmed. Lastly, an approximate frequency solution obtained using the He’s frequency–amplitude formulation was found to produce errors less than 1.0% for a wide range of input parameters. In conclusion, the present study provides a benchmark solution for verification of other approximate solutions, and the He’s frequency–amplitude formulation can be used to obtain fast and accurate solutions for complex nonlinear systems.

Introduction

Nonlinear vibration is a subject of immense scientific and engineering importance. Many physical phenomena such as chaos, ¹ bifurcations,^2,3 jumps, ⁴ sub-harmonic resonance, ⁵ softening and hardening backbone response,^5,6 parametric response, ⁷ stability analysis,^3,5 and fractal oscillations^8–10 can only be understood through nonlinear vibration analysis. The nonlinearity in a vibrating system may arise from static and/or dynamic inertia effects, damping, size-dependent mechanics, fractal mechanics, elastic or spring stiffness, and external excitation. Of these, the most common source of nonlinearity encountered in vibrating systems is the elastic or spring stiffness. The governing model of a nonlinear vibrating system always manifests itself in the form of nonlinear PDEs and ODEs. The latter are exactly solvable only in a few cases. For most cases, approximate analytical or numerical solutions are normally employed. Many approximate analytical methods have been developed^11–16 and are capable of providing physical insight into the response of a nonlinear vibrating system.

Springs or elastic members are an integral part of machines and in some cases a combination of these are required to design machine components. The spring combination can be in series and/or parallel connection. The vibration analysis of a rigid component attached to a number of springs is easier if an equivalent spring can be determined for the spring combination. In the case of linear springs, the exact equivalent spring is simply determined from Newton's third law as can be seen from any introductory vibration textbook, for example, ‘Mechanical Vibrations’ by Rao. ¹⁷ However, in some real life applications one or more springs in the combination are nonlinear and cannot be approximated by a linear spring if an accurate analysis of system’s dynamic response is required. Sanmiguel-Rojas et al ¹⁸ applied the equivalent rigidity method to investigate the free nonlinear vibration of a mass connected to two cubic nonlinear springs in parallel and series. They obtained an approximate solution of the equivalent spring stiffness for both parallel and series connections. They generalized the solution of the parallel connection for any degree of polynomial nonlinearity but could not generalize the solution of the series connection, which was more complicated. In a recent study, ¹⁹ the principle of Newton’s third law was extended to determine the exact equivalent spring for the series connection of two nonlinear springs with any degree of polynomial nonlinearity. The analysis works for two nonlinear springs with the same linear and nonlinear polynomial terms but different coefficients (i.e. stiffness constants). The solution breaks down when the two springs have different number of polynomial terms. Therefore, the general problem of finding the exact equivalent spring for the series connection of any two nonlinear springs with polynomial-type nonlinearities is still unsolved. The implication is that in many cases, a rigid mass attached to two nonlinear springs connected in series cannot be simplified to a system with an exact equivalent spring (at least for now).

A fundamental vibrating system with a simple spring configuration but without an exact equivalent spring is that of a rigid mass attached to a linear and a nonlinear spring connected in series (see Figure 1). This system can be used to study the vibration of a cantilever beam carrying an intermediate lumped mass. ²⁰ Bayat et al ²¹ applied an averaging method to derive an approximate equivalent spring that was applied to study the free vibration of the system in Figure 1. However, an exact analysis of the system in Figure 1 would require developing the vibration model in terms of the motion of the connection point of the springs and the motion of the mass. ²² This gives rise to differential-algebraic equations that can be reduced to a strong nonlinear single degree-of-freedom motion describing the relative displacement of the connection point and the mass. Then, the motions of the mass ( $y$ ) and connection point ( $x$ ) are determined from the solution to the relative motion.

Figure 1.

Schematic of a mass attached to linear and nonlinear springs connected in series.

Perhaps, the first analytical study of Figure 1 was conducted by Telli and Kopmaz. ²² They obtained two approximate solutions using a sixth-order harmonic balance method and a modified Lindstedt–Poincare method. Their analysis showed that the approximate analytical solutions were only accurate for a limited range of inputs, mostly those that gave rise to a weakly nonlinear response. This is not surprising because they employed asymptotic expansion methods which have been shown to be generally inaccurate for strong nonlinear responses. ²³ Lai and Lim ²⁴ improved the classical harmonic balance method in Ref. 22 by applying a linearization to the nonlinear ODE. The improved method called linearized harmonic balance method was used to obtain high-order approximations of the solution that were shown to be algebraically simpler than its classical counterpart. Third-order approximations of the linearized harmonic balance method were derived and shown to predict weakly and strongly nonlinear responses with reasonable accuracy. Bayat et al ²¹ provided an approximate analytical solution to the free vibration of the system in Figure 1 using the Hamiltonian approach. They also derived an equivalent stiffness for the system that gave the same periodic solution as the Hamiltonian approach. The main drawback to the Hamiltonian approach is that it is based on the assumption of a purely harmonic solution and therefore, it is not able to predict the anharmonic response of the relative motion as found in Figures 8 and 11 of Ref. 22. However, it was shown that the Hamiltonian approach could provide accurate solutions for inputs resulting in weakly nonlinear response of the system. A more accurate approximate solution capable of predicting the strong nonlinear response of the system was provided in Ref. 20. The approximate solution was based on a homotopy Pade technique, which is a combination of the homotopy analysis method and the Pade approximation technique. Whereas the results of the homotopy Pade technique were generally found to be in very good agreement with corresponding results of the numerical solution, the process of deriving the solution is complex and would require symbolic computational packages to aid the algebraic manipulation.

In other studies on the free vibration of a mass grounded by linear and nonlinear springs in series, qualitative analysis was conducted. Big-Alabo ²⁵ developed phase portraits and demonstrated the presence of strong nonlinear vibrations around the origin. Furthermore, an equation that can be used to perform a detailed qualitative analysis was derived and very accurate periodic solutions were provided using the continuous piecewise linearization method. A more comprehensive qualitative analysis was presented by Sun and Wu. ²⁶ They demonstrated the existence of various centres of vibration and established some conditions for the existence of periodic solutions. They further demonstrated the existence of asymmetric vibrations in the system, but their analysis did not consider the conditions under which weak or strong nonlinear vibrations can occur. They also derived an approximate periodic solution based on the Newton harmonic balance method, but they did not investigate the accuracy of the method using numerical examples.

The approximate studies that have been conducted on the vibration of the mass grounded by linear and nonlinear springs in series highlight the importance and complex nature of the system. The approximate studies^{18,20–22,24–26} are applicable only for a limited set of inputs, with the possible exception of Ref. 25, and the accuracy of the approximate analytical solutions was verified using numerical solutions because an exact analytical solution has not been derived yet. Furthermore, the conditions for which periodic solutions exist and the general conditions under which weak and strong nonlinear vibrations can occur have not been fully investigated. In this paper, the exact frequency–amplitude response for the free vibration of a mass grounded by linear and nonlinear springs in series was derived from the first integral of the differential equation in terms of the complete elliptic integral of the first and third kinds. Although the complete elliptic integrals are special functions which do not have exact analytical solutions, they can be evaluated to very high accuracy using efficient numerical algorithms found in software packages such as MATLAB and Mathematica. In this study, the EllipticK and EllipticPi functions in Mathematica were used to evaluate the elliptic integral of the first and third kinds to machine epsilon precision. Hence, deriving the exact solution in terms of elliptic integral posed no computational challenge. The derived exact solution enabled the determination of some conditions for which periodic solutions can be obtained. Likewise, some general conditions under which weak and strong nonlinear vibrations occur were established. Furthermore, qualitative analysis was conducted using phase plots and an approximate frequency–amplitude formula was derived.

Nonlinear vibration model

The vibration system consists of a mass ( $m$ ) attached to two springs connected in series as shown in Figure 1. Spring 1 is a spring with linear stiffness $k_1$ for which one end is attached to a support and the other end is connected to spring 2. The latter is a cubic nonlinear spring also connected to the mass and with linear stiffness $k_2$ and cubic nonlinear stiffness $k_3$ . This means that the connection point of the springs can be displaced as the mass is accelerated. If the connection point of the springs is given an initial displacement while the mass is stationary, the forces generated in the springs give rise to the following equation:

$k_{1}x-k_2\left(y-x\right)-k_3{\left(y-x\right)}^3=0$

(1)

which is a mathematical statement of the constrain on springs 1 and 2. On the other hand, if the mass is given an initial displacement, then both the mass and connection point are in motion.

The variational principle for a conservative system of the form $\ddot{u}+f\left(u\right)=0$ is defined as follows ²¹ :

$J={\int }_0^{T/4}\left(-\frac{1}{2}{\dot{u}}^2+F\left(u\right)\right)dt$

(2a)

where $F\left(u\right)=\int f\left(u\right)du$ . The Lagrangian can be obtained from the variational principle as follows:

$L=-\frac{dJ}{dt}=\frac{1}{2}{\dot{u}}^2-F\left(u\right)$

(2b)

in equation (2b), the first term is the kinetic energy per unit mass while the second term is the potential energy per unit mass. From Figure 1,

$L=\frac{1}{2}{\dot{y}}^2-\frac{k_2}{2m}{\left(y-x\right)}^2-\frac{k_3}{4m}{\left(y-x\right)}^4=\frac{1}{2}{\left(\dot{u}+\dot{x}\right)}^2-\frac{k_2}{2m}{u}^2-\frac{k_3}{4m}{u}^4$

(3a)

where $u=y-x$ . The vibration model for the mass can be derived from the Lagrange equation:

$\frac{d}{dt}\left(\frac{\partial L}{\partial \dot{u}}\right)=\frac{\partial L}{\partial u}$

(3b)

Hence, substituting equation (3a) in (3b) gives the following:

$m\left(\frac{d^{2}u}{d{t}^2}+\frac{d^{2}x}{d{t}^2}\right)+k_{2}u+k_3{u}^3=0$

(4)

Note that equation (4) can be derived also by considering the balance of forces. Equations (1) and (4) are a pair of differential-algebraic equations that can be solved simultaneously using numerical methods. Alternatively, we can reduce the differential-algebraic equations to a single degree-of-freedom in terms of the relative motion as follows. First, equation (1) can be written as follows:

$x=\xi u+\epsilon \xi u^3$

(5)

where $\xi =k_2/k_1$ and $\epsilon =k_3/k_2$ . Then, substituting equation (5) in (4) gives the following equation:

$m\left(1+\xi +3\epsilon \xi u^2\right)\frac{d^{2}u}{d{t}^2}+6m\epsilon \xi u{\left(\frac{du}{dt}\right)}^2+k_{2}u+k_3{u}^3=0$

(6)

Now, dividing equation (6) by $m\left(1+\xi \right)$ we get the following:

$\left(1+3\eta \epsilon u^2\right)\frac{d^{2}u}{d{t}^2}+6\eta \epsilon u{\left(\frac{du}{dt}\right)}^2+{\omega }_0^2\left(u+\epsilon u^3\right)=0$

(7)

where $\eta =\xi /(1+\xi )=k_2/\left(k_1+k_2\right)$ , ${\omega }_0^2=k_2/\left[m\left(1+\xi \right)\right]$ , and the initial conditions are: $u\left(0\right)=\mathrm{A}$ and $du\left(0\right)/dt=0$ . In non-dimensional form, equation (7) can be expressed as follows:

$\left(1+3\eta \epsilon A^2{v}^2\right)\frac{d^{2}v}{d{\tau }^2}+6\eta \epsilon A^{2}v{\left(\frac{dv}{d\tau }\right)}^2+v+\epsilon A^2{v}^3=0$

(8)

where $\tau ={\omega }_{0}t$ and $v=u/A$ are the dimensionless time and relative displacement, respectively. The initial conditions for equation (8) are as follows: $v\left(0\right)=1$ and $dv\left(0\right)/d\tau =0$ . Equation (8) shows that the vibration model for the mass grounded by linear and nonlinear springs in series is characterized by inertia and static nonlinearities. ²⁰ The static nonlinearity is a cubic nonlinearity while the inertia nonlinearity consists of a quadratic nonlinearity and a dynamic nonlinearity, which takes the form of a quadratic damping. It should be noted that the quadratic damping is not a dissipative term but rather it adds to the overall restoring force of the system.

Exact analytical solution

Derivation of the exact frequency–amplitude response

For a nonlinear oscillator of the form:

$I\left(v\right)\frac{d^{2}v}{d{\tau }^2}+k\left(\frac{dI}{dv}\right){\left(\frac{dv}{d\tau }\right)}^2=Q\left(v\right)$

(9)

where $v\left(0\right)=A$ and $dv\left(0\right)/d\tau =0$ , the state-space relationship can be expressed as follows ²⁶ :

$\frac{dv}{d\tau }=\pm \sqrt{\frac{2\left[h\left(A\right)-h\left(v\right)\right]}{{\left[I\left(v\right)\right]}^{2k}}}$

(10)

where $h\left(v\right)=-\int {\left[I\left(v\right)\right]}^{2k-1}Q\left(v\right)dv$ . For equation (2), $I\left(v\right)=1+3\eta \epsilon A^2{v}^2$ , $Q\left(v\right)=-\left(v+\epsilon A^2{v}^3\right)$ , and $k=1$ . Therefore,

$\frac{dv}{d\tau }=\pm \frac{\sqrt{2\left(1-v^2\right)+\left(1+3\eta \right)\epsilon A^2\left(1-v^4\right)+2\eta {\epsilon }^2{A}^4\left(1-v^6\right)}}{\sqrt{2}\left(1+3\eta \epsilon A^2{v}^2\right)}$

(11)

Considering the first quarter of the vibration cycle (the velocity is negative) gives:

${\int }_0^{\tau }\tau d\tau =\sqrt{2}{\int }_{v_{\tau }}^1\frac{1+3\eta \epsilon A^2{v}^2}{\sqrt{2\left(1-v^2\right)+\left(1+3\eta \right)\epsilon A^2\left(1-v^4\right)+2\eta {\epsilon }^2{A}^4\left(1-v^6\right)}}dv$

(12)

where $1\ge v_{\tau }\ge 0$ is the dimensionless relative displacement at time $\tau$ . At the end of the first quarter of the vibration cycle, $\tau =T/4$ and $v_{\tau }=0$ . Then, the exact time period can be obtained from equation (12) as follows:

$T=4\sqrt{2}{\int }_0^1\frac{1+3\eta \epsilon A^2{v}^2}{\sqrt{2\left(1-v^2\right)+\left(1+3\eta \right)\epsilon A^2\left(1-v^4\right)+2\eta {\epsilon }^2{A}^4\left(1-v^6\right)}}dv$

(13)

The main goal of the present article is to derive the exact analytical solution to equation (13).

The integral for the exact period contains the square root of a rational polynomial function. In order to derive the exact period, equation (13) has to be transformed into a standard definite integral that is expressible in terms of a complete elliptic integral solution. Hence, by algebraic simplification we have:

$T=4\sqrt{\frac{1}{\eta {\epsilon }^2{A}^4}}{\int }_0^1\frac{1+3\eta \epsilon A^2{v}^2}{\sqrt{\left(a_0+a_2{v}^2+v^4\right)\left(1-v^2\right)}}dv$

(14)

where $a_0=a_2+1/\left(\eta {\epsilon }^2{A}^4\right)$ and $a_2=1+\left(1+3\eta \right)/2\eta {\epsilon A}^2$ . Applying a change of variable with $z=v^2$ gives:

$T=2\sqrt{\frac{1}{\eta {\epsilon }^2{A}^4}}{\int }_0^1\frac{1+3\eta \epsilon A^{2}z}{\sqrt{\left(1-z\right)z\left(a_0+a_{2}z+z^2\right)}}dz$

(15)

The discriminant of the quadratic expression in the square root of equation (15) is negative (since $4a_0>a_2^2$ ), which implies that it has complex factors as follows:

$a_0+a_{2}z+z^2=\left(z+\frac{a_2}{2}+\frac{i\sqrt{4a_0-a_2^2}}{2}\right)\left(z+\frac{a_2}{2}-\frac{i\sqrt{4a_0-a_2^2}}{2}\right)$

(16)

Hence, putting equation (16) in (15) gives:

$T=2\sqrt{\frac{1}{\eta {\epsilon }^2{A}^4}}\left[{\int }_0^1\frac{dz}{\sqrt{\left(1-z\right)z\left(z-{\sigma }_1\right)\left(z-{\sigma }_2\right)}}+3\eta \epsilon A^2{\int }_0^1\frac{zdz}{\sqrt{\left(1-z\right)z\left(z-{\sigma }_1\right)\left(z-{\sigma }_2\right)}}\right]$

(17)

where ${\sigma }_1=-\left(a_2+i\sqrt{4a_0-a_2^2}\right)/2$ and ${\sigma }_2=-\left(a_2-i\sqrt{4a_0-a_2^2}\right)/2$ are complex conjugates. To solve equation (17), we refer to the table of elliptic integrals by Byrd and Friedman. ²⁷ Specifically, the result of formula 259.03 on page 133 of Ref. 27 that is applicable to the integral of a quartic polynomial where two factors are complex conjugates is of interest and is given as follows:

${\int }_b^y\frac{z^{r}dz}{\sqrt{\left(a-z\right)\left(z-b\right)\left(z-c\right)\left(z-\overline{c}\right)}}=\frac{G{\left(aB+bA\right)}^r}{{\left(A-B\right)}^r}{\displaystyle \sum }_{j=0}^r\frac{{{\alpha }_2^{r-j}\left(\alpha -{\alpha }_2\right)}^{j}r!}{\left(r-j\right)!j!}{\int }_0^{\psi }\frac{d\theta }{{\left(1+\alpha \hspace{0.17em}\text{cn}\hspace{0.17em}\left(\theta \right)\right)}^j}$

(18)

where $r\ge 0$ , $a\ge y>b$ , $\{a,b\}\mathbb{\in }\mathbb{R}$ , and $\{c,\overline{c}\}\mathbb{\in }\mathbb{C}$ represents a pair of complex conjugates. The parameters in equation (18) are defined as follows:

$\begin{array}{l}\alpha =\left(A-B\right)/\left(A+B\right)\hspace{0.17em}\\ {\alpha }_2=\left(bA-aB\right)/\left(bA+aB\right)\hspace{0.17em}\\ A=\sqrt{{\left(a-b_1\right)}^2+b_2^2}\hspace{0.17em}\\ B=\sqrt{{\left(b-b_1\right)}^2+b_2^2}\hspace{0.17em}\\ G={\left(AB\right)}^{-1/2}\hspace{0.17em}\\ b_1=\left(c+\overline{c}\right)/2\hspace{0.17em}\\ b_2^2=-{\left(c-\overline{c}\right)}^2/4\hspace{0.17em}\\ \text{cn}\hspace{0.17em}\psi =\cos \hspace{0.17em}\phi =\frac{\left(a-y\right)B-\left(y-b\right)A}{\left(a-y\right)B+\left(y-b\right)A}\hspace{0.17em}\end{array}\hspace{0.17em}\}$

(19)

where $\phi =\text{am}\hspace{0.17em}\left(\psi \right)$ is the amplitude of $\psi$ . Equation (18) transforms the definite integral of an irrational square root of a fourth-degree polynomial with a pair of conjugate complex zeros to the definite integral of a function of the Jacobi cosine function. As demonstrated afterwards, the latter integral is expressible in terms of one or more of the three standard elliptic integrals depending on the value of $r$ .

The left-hand side of equation (18) shows that $r=0$ and $r=1$ give the same form as the first and second integrals of the exact time period (see equation (17)). Hence, for $r=0$ we have

${\int }_b^y\frac{dz}{\sqrt{\left(a-z\right)\left(z-b\right)\left(z-c\right)\left(z-\overline{c}\right)}}=G{\int }_0^{\psi }d\theta =G\times F\left(\phi ,m\right)$

(20)

where $F\left(\phi ,m\right)$ is the incomplete elliptic integral of the first kind defined as follows:

$F\left(\phi ,m\right)={\int }_0^{\phi }{\left(1-m\hspace{0.17em}{\sin }^2\hspace{0.17em}\theta \right)}^{-1/2}d\theta$

(21)

and the arguments of $F\left(\phi ,m\right)$ are defined as follows:

$\phi ={\cos }^{-1}\left[\frac{\left(a-y\right)B-\left(y-b\right)A}{\left(a-y\right)B+\left(y-b\right)A}\right]$

(22a)

$m=\frac{{\left(a-b\right)}^2-{\left(A-B\right)}^2}{4AB}$

(22b)

Comparing the first integral in equation (17) with equation (20), we get $a=1$ , $b=0$ , $c={\sigma }_1$ , $\overline{c}={\sigma }_2$ , and $y=1$ . Substituting these values into equations (19) and (22) gives $A=\sqrt{1+a_0+a_2}$ , $B=\sqrt{a_0}$ , $G={\left[a_0\left(1+a_0+a_2\right)\right]}^{-1/4}$ , $\alpha =\left(\sqrt{1+a_0+a_2}-\sqrt{a_0}\right)/\left(\sqrt{1+a_0+a_2}+\sqrt{a_0}\right)$ , ${\alpha }_2=-1$ , $\phi =\pi$ , and

$m=\frac{1}{2}-\frac{2a_0+a_2}{4\sqrt{a_0\left(1+a_0+a_2\right)}}$

Therefore,

${\int }_0^1\frac{dz}{\sqrt{\left(1-z\right)z\left(z-{\sigma }_1\right)\left(z-{\sigma }_2\right)}}=\frac{2}{\sqrt[4]{a_0\left(1+a_0+a_2\right)}}K\left(\frac{1}{2}-\frac{2a_0+a_2}{4\sqrt{a_0\left(1+a_0+a_2\right)}}\right)$

(23)

Note that the identity $F\left(\pi ,m\right)=2K\left(m\right)$ was taken into account to arrive at equation (23).

Next, we treat the second integral in equation (17) by substituting $r=1$ in equation (18) to get:

${\int }_b^y\frac{zdz}{\sqrt{\left(a-z\right)\left(z-b\right)\left(z-c\right)\left(z-\overline{c}\right)}}=\frac{G\left(aB+bA\right)}{\left(A-B\right)}\left[{\alpha }_{2}F\left(\phi ,m\right)+\left(\alpha -{\alpha }_2\right){\int }_0^{\psi }\frac{d\theta }{1+\alpha \hspace{0.17em}\text{cn}\hspace{0.17em}\left(\theta \right)}\right]$

(24)

Comparing the second integral in equation (17) with equation (24), $a=1$ , $b=0$ , $c={\sigma }_1$ , $\overline{c}={\sigma }_2$ , and $y=1$ . These results are the same as those obtained for the first integral in equation (17). So, the previous results derived for the parameters in equations (19) and (22) are applicable here. Consequently, we have:

${\int }_0^1\frac{zdz}{\sqrt{\left(1-z\right)z\left(z-{\sigma }_1\right)\left(z-{\sigma }_2\right)}}=\frac{2}{\left(\sqrt{1+a_0+a_2}-\sqrt{a_0}\right)}{\left(\frac{a_0}{1+a_0+a_2}\right)}^{1/4}\left[\left(\frac{\sqrt{1+a_0+a_2}}{\sqrt{1+a_0+a_2}+\sqrt{a_0}}\right){\int }_0^{\psi }\frac{d\theta }{1+\alpha \hspace{0.17em}\text{cn}\hspace{0.17em}\left(\theta \right)}-K\left(\frac{1}{2}-\frac{2a_0+a_2}{4\sqrt{a_0\left(1+a_0+a_2\right)}}\right)\right]$

(25)

The integral in equation (25) can be evaluated based on formula 341.03 on page 206 of Ref 27 as follows:

${\int }_0^{\psi }\frac{d\theta }{1+\alpha \hspace{0.17em}\text{cn}\hspace{0.17em}\left(\theta \right)}=\frac{1}{1-{\alpha }^2}\left[\mathrm{\Pi }\left(\phi ,\frac{{\alpha }^2}{{\alpha }^2-1},m\right)-\alpha f_1\right]$

(26)

where $\mathrm{\Pi }\left(\phi ,n,m\right)$ is the incomplete elliptic integral of the third kind defined as follows:

$\mathrm{\Pi }\left(\phi ,n,m\right)={\int }_0^{\phi }\frac{d\varphi }{\left(1-n\hspace{0.17em}{\sin }^2\hspace{0.17em}\varphi \right)\sqrt{1-m\hspace{0.17em}{\sin }^2\hspace{0.17em}\varphi }}$

(27)

and

$f_1=\{\begin{array}{l}{\left[\frac{1-{\alpha }^2}{m+m^{\prime }{\alpha }^2}\right]}^{1/2}\hspace{0.17em}{\tan }^{-1}\left({\left[\frac{m+m^{\prime }{\alpha }^2}{1-{\alpha }^2}\right]}^{1/2}\text{sd}\psi \right)\hspace{0.17em}\text{if}\hspace{0.17em}\frac{{\alpha }^2}{{\alpha }^2-1}<m\hspace{0.17em}\\ \text{sd}\psi \hspace{0.17em}\text{if}\hspace{0.17em}\frac{{\alpha }^2}{{\alpha }^2-1}=m\hspace{0.17em}\\ \frac{1}{2}{\left[\frac{{\alpha }^2-1}{m+m^{\prime }{\alpha }^2}\right]}^{1/2}\hspace{0.17em}{\log }_e\left[\frac{\sqrt{m+m^{\prime }{\alpha }^2}\text{dn}\psi -\sqrt{{\alpha }^2-1}\text{sn}\psi }{\sqrt{m+m^{\prime }{\alpha }^2}\text{dn}\psi +\sqrt{{\alpha }^2-1}\text{sn}\psi }\right]\hspace{0.17em}\text{if}\hspace{0.17em}\frac{{\alpha }^2}{{\alpha }^2-1}>m\hspace{0.17em}\end{array}$

(28)

In equations (28a)–(28c), $m^{\prime }=1-m$ , $\text{sn}$ is the Jacobi sine function, $\text{dn}$ is the Jacobi delta function, and $\text{sd}=\text{sn}/\text{dn}$ . Since $A>B>0$ , then $\alpha <1$ . This implies that ${\alpha }^2/\left({\alpha }^2-1\right)<0<m$ and equation (28a) is applicable to equation (26). Hence,

${\int }_0^{\psi }\frac{d\theta }{1+\alpha \hspace{0.17em}\text{cn}\hspace{0.17em}\left(\theta \right)}=\frac{1}{1-{\alpha }^2}\left[\mathrm{\Pi }\left(\phi ,\frac{{\alpha }^2}{{\alpha }^2-1},m\right)-\alpha {\left[\frac{1-{\alpha }^2}{m+m^{\prime }{\alpha }^2}\right]}^{1/2}\hspace{0.17em}{\tan }^{-1}\left({\left[\frac{m+m^{\prime }{\alpha }^2}{1-{\alpha }^2}\right]}^{1/2}\text{sd}\psi \right)\right]$

(29)

Using the identities ${\text{sn}}^2+{\text{cn}}^2=1$ and ${\text{dn}}^2+m\hspace{0.17em}{\text{sn}}^2=1$ , we get $\text{sn}\hspace{0.17em}\left(\psi \right)=0$ , $\text{dn}\hspace{0.17em}\psi =1$ , and $\text{sd}\hspace{0.17em}\psi =0$ . Taking these results into account and substituting the expressions for the parameters $\phi$ , $\alpha$ , and $m$ gives:

${\int }_0^{\psi }\frac{d\theta }{1+\alpha \hspace{0.17em}\text{cn}\hspace{0.17em}\left(\theta \right)}=\frac{{\left(\sqrt{1+a_0+a_2}+\sqrt{a_0}\right)}^2}{2\sqrt{a_0\left(1+a_0+a_2\right)}}\mathrm{\Pi }\left(-\frac{{\left(\sqrt{1+a_0+a_2}-\sqrt{a_0}\right)}^2}{4\sqrt{a_0\left(1+a_0+a_2\right)}},\frac{1}{2}-\frac{2a_0+a_2}{4\sqrt{a_0\left(1+a_0+a_2\right)}}\right)$

(30)

where the identity $\mathrm{\Pi }\left(\pi ,n,m\right)=2\mathrm{\Pi }\left(n,m\right)$ has been applied. Substituting equation (30) in (25) and simplifying gives:

${\int }_0^1\frac{zdz}{\sqrt{\left(1-z\right)z\left(z-{\sigma }_1\right)\left(z-{\sigma }_2\right)}}=\frac{2}{\left(\sqrt{1+a_0+a_2}-\sqrt{a_0}\right)}{\left(\frac{a_0}{1+a_0+a_2}\right)}^{1/4}\hspace{0.17em}\left[\left(\frac{\sqrt{1+a_0+a_2}+\sqrt{a_0}}{2\sqrt{a_0}}\right)\mathrm{\Pi }\left(-\frac{{\left(\sqrt{1+a_0+a_2}-\sqrt{a_0}\right)}^2}{4\sqrt{a_0\left(1+a_0+a_2\right)}},\frac{1}{2}-\frac{2a_0+a_2}{4\sqrt{a_0\left(1+a_0+a_2\right)}}\right)-K\left(\frac{1}{2}-\frac{2a_0+a_2}{4\sqrt{a_0\left(1+a_0+a_2\right)}}\right)\right]$

(31)

Then, from equations (17), (23), and (31) the exact period was derived as follows:

$T=4\sqrt{\frac{1}{\eta {\epsilon }^2{A}^4}}\left[\left(\frac{1}{\sqrt[4]{q_1{q}_2}}-p\right)K\left(m\right)+p\left(\frac{\sqrt{q_2}+\sqrt{q_1}}{2\sqrt{q_1}}\right)\hspace{0.17em}\mathrm{\Pi }\left(n,m\right)\right]$

(32)

where

$m=\frac{1}{2}-\frac{q_2+q_1-1}{4\sqrt{q_1{q}_2}}\hspace{0.17em}$

(33)

$n=-\frac{{\left(\sqrt{q_2}-\sqrt{q_1}\right)}^2}{4\sqrt{q_1{q}_2}}$

(34)

$p=\frac{3\eta \epsilon A^2}{\left(\sqrt{q_2}-\sqrt{q_1}\right)}{\left(\frac{q_1}{q_2}\right)}^{1/4}$

(35)

and the $q_i$ (for $i=1,2$ ) parameters are given as follows:

$q_1=a_0=1+\frac{1+3\eta }{2\eta {\epsilon A}^2}+\frac{1}{\eta {\epsilon }^2{A}^4}$

(36)

$q_2=1+a_0+a_2=3+\frac{1+3\eta }{\eta \epsilon A^2}+\frac{1}{\eta {\epsilon }^2{A}^4}$

(37)

Hence, the exact frequency–amplitude solution for the periodic response can be expressed as follows:

${\omega }_{ex}=\frac{2\pi }{T}=\frac{\pi \sqrt{\eta {\epsilon }^2{A}^4}}{2\left[\left(\frac{1}{\sqrt[4]{q_1{q}_2}}-p\right)K\left(m\right)+p\left(\frac{\sqrt{q_2}+\sqrt{q_1}}{2\sqrt{q_1}}\right)\hspace{0.17em}\mathrm{\Pi }\left(n,m\right)\right]}$

(38)

It should be noted that equation (38) gives the exact non-dimensional frequency, that is, exact frequency divided by the constant ${\omega }_0$ , where ${\omega }_0$ is the frequency of the system when $\epsilon =0$ .

Conditions for periodic solutions

Having derived the exact frequency–amplitude solution in a closed-form, we now investigate the conditions under which periodic solutions exist. The square root in the numerator of equation (38) suggests that real periodic solutions exist for this system for positive and negative values of $\epsilon$ , and for positive values of $\eta$ only. Physically speaking, the parameter $\epsilon$ determines the nonlinear strength of the nonlinear spring while the parameter $\eta$ is a dimensionless factor that compares the linearity of spring 2 to that of spring 1. As $\epsilon$ increases, the nonlinear strength of spring 2 increases, and vice versa. Also, $\epsilon >0$ represents a spring with hardening nonlinearity while $\epsilon <0$ is for softening nonlinearity. For the $\eta$ parameter and if $k_1,k_2>0$ , then $k_1>k_2$ implies that $0<\eta <1/2$ , $k_1=k_2$ implies that $\eta =1/2$ , and $k_1<k_2$ implies that $1/2<\eta <1$ . The case of $\eta >1$ occurs when $\xi =k_2/k_1<-1$ and the condition $\xi <-1$ is satisfied when $k_1<-k_2$ . Considering that ${\omega }_0=\sqrt{k_2/\left[m\left(1+\xi \right)\right]}$ , then $k_1$ must be positive and $k_2$ must be negative for ${\omega }_0$ to have real values. This implies that spring 2 is capable of introducing a bifurcation in the vibration of the system when the linearity factor is greater than unity. These conditions define when periodic oscillations can occur and were taken into consideration while simulating the frequency–amplitude response of the system.

Approximate analytical solution for the frequency–amplitude response

Many approximate methods exist for deriving the frequency–amplitude response of nonlinear conservative oscillators. In this section, the frequency–amplitude formulation (FAF) by He ²⁸ has been applied. This approach provides a very fast means of obtaining the frequency–amplitude solution and produces the same results with many other established approximate analytical methods such as Hamiltonian approach and homotopy perturbation method. Since its formulation, a number of improvements have been proposed by several authors.^29–34 The improvements are normally carried out on the residual function, ³⁵ averaging technique for error minimization,^32,35 initial trial functions,^30,33,34 and/or the form of the expression for the FAF. ³¹ Here, we apply the classical approach of the FAF with weighted-average residuals ²⁸ to get a fast approximation of the frequency–amplitude solution to equation (8).

Since equation (8) is in a non-dimensional form, we assume two initial trial solutions of the form:

$v_1\left(\tau \right)=\cos \hspace{0.17em}\tau \hspace{0.17em}\text{and}\hspace{0.17em}{v}_2\left(\tau \right)=\cos {\omega }_a\tau$

(39)

where ${\omega }_a$ is the approximate frequency of the FAF, which is given as follows:

${\omega }_a=\sqrt{\frac{{\overline{R}}_2-{\omega }_a^2{\overline{R}}_1}{{\overline{R}}_2-{\overline{R}}_1}}$

(40)

and ${\overline{R}}_1$ and ${\overline{R}}_2$ are the residual averages defined as follows:

${\overline{R}}_1=\frac{2}{\pi }{\int }_0^{2/\pi }\mathfrak{D}\left(v_1\left(\tau \right)\right)v_1\left(\tau \right)d\tau$

(41a)

${\overline{R}}_2=\frac{2{\omega }_a}{\pi }{\int }_0^{2{\omega }_a/\pi }\mathfrak{D}\left(v_2\left(\tau \right)\right)v_2\left(\tau \right)d\tau$

(41b)

In equations (41a), (41b) $\mathfrak{D}$ is a differential function obtained by making the right-hand side of the vibration (or oscillator) equation zero. Hence,

$\mathfrak{D}\left(v\right)=\left(1+3\eta \epsilon A^2{v}^2\right)\frac{d^{2}v}{d{\tau }^2}+6\eta \epsilon A^{2}v{\left(\frac{dv}{d\tau }\right)}^2+v+\epsilon A^2{v}^3$

(42)

Substituting equations (42) and (39) in (41a), (41b) we get:

${\overline{R}}_1=\frac{3\epsilon A^2\left(1-\eta \right)}{8}$

(43a)

${\overline{R}}_2=\frac{3\epsilon A^2\left(1-\eta {\omega }_a^2\right)+4\left(1-{\omega }_a^2\right)}{8}$

(43b)

Then, from equations (40) and (43a) and (43b) the approximate frequency based on the FAF was derived as follows:

${\omega }_a=\sqrt{\frac{4+3\epsilon A^2}{4+3\eta \epsilon A^2}}$

(44)

Equation (44) is the same result derived if we solve equation (8) using the Hamiltonian approach. ²¹

In contrast to the conditions for periodic oscillations deduced from the exact solution, equation (44) implies the possibility of periodic solutions for $\eta <0$ and predicts that periodic solutions do not exist when $\epsilon =-4/\left(3A^2\right)$ and $\epsilon =-4/\left(3\eta A^2\right)$ . In the former case of $\epsilon$ , the approximate solution predicts the occurrence of buckling failure, whereas in the latter case of $\epsilon$ it predicts a singularity in the frequency response. These predictions do not agree with the observed behaviour of the system as demonstrated in the next Section. In many circumstances, approximate analytical solutions are used to predict nonlinear vibration response and gain physical insight into the behaviour of the system. However, a comparison of the predicted physics of the exact and approximate analytical solutions of the present system shows that the latter can give misleading conclusions. This makes it necessary to validate physical behaviour deduced from approximate analytical solutions using alternative methods.

Results and discussions

Exact frequency–amplitude response

The complete elliptic integrals of the first and third kind can be evaluated to machine precision using the EllipticK and EllipticPi functions in Mathematica™, respectively. Hence, equation (38) was simulated using these Mathematica™ functions to demonstrate the effect of the system parameters ( $\epsilon$ , $\eta$ , and $A$ ) on the frequency–amplitude response. The input values used for the simulations are $m=1.0$ , $k_1=20.0$ , $k_2=20.0$ , and $k_3=40.0$ except otherwise stated. For each simulation, the absolute error of the numerical solution (i.e. $ϵ=\left|{\omega }_{ex}-{\omega }_{num}\right|$ ) was calculated and plotted. The numerical solution was obtained by using the NIntegrate function in Mathematica™ to evaluate equation (13).

Figures 2–5 show the dependence of the frequency–amplitude response on the nonlinear factor, $\epsilon$ . For weak (Figure 2) and moderate (Figure 3) hardening nonlinearity, the frequency depends strongly on $\epsilon$ and $A$ , while for strong hardening nonlinearity (Figure 4), the dependence of frequency on $\epsilon$ is appreciable only during small to moderate amplitudes. For large amplitude vibrations under strong hardening nonlinearity, the frequency is approximately constant at ${\omega }_{ex}=\sqrt{2}$ . On the other hand, a strong dependence of the frequency on $\epsilon$ and $A$ was observed for a softening nonlinearity (Figure 5). The effect of an increase in $\epsilon$ was observed to be an increase in the frequency for a hardening nonlinearity and a decrease in the frequency for a softening nonlinearity. These effects were more pronounced at higher amplitudes.

Figure 2.

Frequency–amplitude response and error in numerical solution for $0.01\le \epsilon \le 0.1$ .

Figure 3.

Frequency–amplitude response and error in numerical solution for  $0.1\le \epsilon \le 1.0$ .

Figure 4.

Frequency–amplitude response and error in numerical solution for  $1.0\le \epsilon \le 10.0$ .

Figure 5.

Frequency–amplitude response and error in numerical solution for  $-10.0\le \epsilon <0.0$ .

Figures 6–8 illustrate the dependence of the frequency–amplitude response on the linearity factor, $\eta$ . A strong dependence of the frequency on $\eta$ and $A$ was observed for small, moderate, and large values of $\eta$ . The effect of an increase in $\eta$ is a decrease in the frequency, and this effect is more evident at larger amplitudes. It was observed that when $\eta =1.0$ , the frequency response is constant for all amplitudes. The implication is that a linear response is produced when $\eta =1.0$ . This observation can also be explained by noting that as $\eta \to 1.0$ , $k_2\gg k_1$ and $\epsilon =k_3/k_2\to 0$ . Hence, substituting $\epsilon =0$ in equation (8) reduces the nonlinear equation of motion to a linear one with a frequency of ${\omega }_{ex}=1.0$ , which is what we see in Figures 7 and 8 when $\eta =1.0$ .

Figure 6.

Frequency–amplitude response and error in numerical solution for  $0.01\le \eta \le 0.1$ .

Figure 7.

Frequency–amplitude response and error in numerical solution for $0.1\le \eta \le 1.0$ .

Figure 8.

Frequency–amplitude response and error in numerical solution for  $1.0\le \eta \le 10.0$ .

Furthermore, it was observed that the frequency increased when $\epsilon >0$ and reduced when $\epsilon <0$ for a given amplitude (see Figures 2–5). This indicates that a bifurcation point exists at $\epsilon =0$ . Bifurcation is the change in the qualitative behaviour of a dynamic system as a parameter of the system is varied. The point at which the change occurs is known as the bifurcation point. As the parameter $\epsilon$ moves from positive to negative values, the system changes from hardening to softening nonlinearity. At $\epsilon =0$ , the system behaves as a linear oscillator with purely harmonic response. Similarly, it was observed that the frequency increased with amplitude when $\eta <1.0$ and decreased with amplitude when $\eta >1.0$ (see Figures 6–8). Therefore, a bifurcation point exists at $\eta =1.0$ . As $\eta$ changes from $\eta <1.0$ to $\eta >1.0$ , spring 2 causes the system’s response to change from a hardening response to a softening response. The change in the response characteristics accounts for the bifurcation observed with the change in the $\eta$ parameter.

Finally, Figures 2–8 show the absolute error obtained by solving equation (13) numerically. These plots show that the numerical error was oscillatory in nature and generally negligible (close to machine epsilon precision) for most cases. However, there were instances where the numerical solution produced significant errors and the frequency estimate was accurate to less than five decimal places (see Figures 4, 7 and 8). Therefore, the numerical method was not found to be reliable for all possible frequency solutions of the system under consideration.

Phase response

The phase response of the relative motion can be obtained from equation (11), and phase diagrams for different conditions of the present system are shown in Figures 9–12. The ratio of the maximum velocity to the maximum displacement of the system, that is $r_{\mathit{max}}=\sqrt{1+\left(1+3\eta \right)\epsilon A^2/2+\eta {\epsilon }^2{A}^4}$ , can be used to characterize its qualitative response. If $r_{\mathit{max}}=1.0$ , which happens when $\epsilon =0$ , then we have a purely harmonic response and the frequency is independent of the amplitude. The implication is that for $r_{\mathit{max}}$ close to unity, the system exhibits a weakly nonlinear response. For ratios far away from unity, that is $r_{\mathit{max}}\ll 1.0$ or $r_{\mathit{max}}\gg 1.0$ , the system exhibits a strong nonlinear response. In other words, one can apply the rule of thumb that if $0.8\le r_{\mathit{max}}\le 1.2$ the system’s response is weakly nonlinear and if otherwise it is moderately or strongly nonlinear.

Figure 9.

Phase portrait. Left: $\eta =0.75$ ; $\epsilon =0.5$ ; $A=1.0$ and right: $\eta =0.75$ ; $\epsilon =2.0$ ; $A=1.0$ .

Figure 10.

Phase portrait. Left: $\eta =1.0$ ; $\epsilon =0.0$ ; $A=1.0$ ; middle: $\eta =1.5$ ; $\epsilon =0.5$ ; $A=1.0$ and right: $\eta =1.5$ ; $\epsilon =2.0$ ; $A=1.0$ .

Figure 11.

Phase portrait. Left: $\eta =0.25$ ; $\epsilon =-1.0$ ; $A=1.0$ and right: $\eta =0.75$ ; $\epsilon =-1.0$ ; $A=1.0$ .

Figure 12.

Phase portrait. Left: $\eta =0.75$ ; $\epsilon =-1.0$ ; $A=0.50$ and right: $\eta =1.5$ ; $\epsilon =-1.0$ ; $A=0.50$ .

In Figure 9, the nonlinearity parameter ( $\epsilon$ ) was increased for $\eta =0.75<1.0$ and $A=1.0$ while in Figure 10 it was increased for $\eta =1.5>1.0$ and $A=1.0$ . The phase diagrams show a stronger nonlinear response with increase in $\epsilon$ . The ratio of the maximum velocity to the maximum displacement as $v\to 1.0$ is significantly greater than one which confirms the presence of a strong nonlinear response. Also, in Figure 10, the phase diagram for $\eta =1.0$ was shown and reveals a purely harmonic response of the system. This confirms the existence of a bifurcation point at $\eta =1.0$ because there is a change from nonlinear response to linear response and back to nonlinear response as $\eta$ increases from less than one to greater than one. Again this observation corroborates the results in Figures 6–8.

Figures 11 and 12 show the phase diagram of the system when $\epsilon <0$ , which implies a softening nonlinearity. The phase diagram shows a change from a single centre oscillation about the origin to a two centre oscillation as $\eta$ increases. This observation illustrates the existence of bifurcation points in the parameter, $\eta$ . For Figure 11, the bifurcation occurs at $\eta =0.41165$ while in Figure 12 it occurs at $\eta =1.3334$ . These bifurcation points are different from the ones identified in Figures 6–8 and demonstrate the complex nature of the present system. It can be concluded from the phase diagrams that the bifurcation of the system with respect to any of its three parameters (i.e. $\epsilon$ , $\eta$ , and $A$ ) also depends on the values of the other two parameters. A complete bifurcation analysis of the present oscillator is beyond the present scope and is left for future work.

Approximate frequency–amplitude response

In a previous section an approximate analytical solution for the nonlinear frequency–amplitude response was derived using He’s FAF. In this section, the accuracy of the approximate frequency–amplitude solution was investigated by comparing with the exact solution as shown in Tables 1 and 2. The results in Table 1 are for the case of hardening nonlinearity ( $\epsilon >0$ ) while the results in Table 2 are for the case of softening nonlinearity ( $\epsilon <0$ ). In both cases, the approximate frequency–amplitude solution gives a maximum relative error of less that 1.0%, which is remarkable considering the simplicity of the method used and the frequency–amplitude formula derived. Therefore, we concluded that the frequency–amplitude formula is capable of providing quick and simple estimates of the nonlinear frequency response of complex oscillators.

Table 1.

Exact ( ${\omega }_{ex}$ ) and approximate ( ${\omega }_a$ ) frequency: $\eta =0.75$ and $\epsilon =1.5$ .

A	${\omega }_a$	${\omega }_{ex}$	% Error FAF
0.1	1.00139	1.00140	0.00100
0.2	1.00543	1.00548	0.00497
0.3	1.01169	1.01192	0.02273
0.4	1.01963	1.02020	0.05587
0.5	1.02862	1.02968	0.10294
0.6	1.03810	1.03974	0.15773
0.7	1.04762	1.04984	0.21146
0.8	1.05683	1.05958	0.25954
0.9	1.06552	1.06870	0.29756
1.0	1.07357	1.07706	0.32403
2.0	1.12122	1.12433	0.27661
3.0	1.13778	1.13966	0.16496
4.0	1.14470	1.14589	0.10385
5.0	1.14815	1.14895	0.06963
6.0	1.15009	1.15067	0.05041
7.0	1.15129	1.15172	0.03734
8.0	1.15207	1.15241	0.02950
9.0	1.15262	1.15288	0.02255
10.0	1.15301	1.15323	0.01908

Table 2.

Exact ( ${\omega }_{ex}$ ) and approximate ( ${\omega }_a$ ) frequency: $\eta =0.25$ and $\epsilon =-1.0$ .

A	${\omega }_a$	${\omega }_{ex}$	% Error FAF
0.1	0.997178	0.997179	0.00010
0.2	0.988600	0.988618	0.00182
0.3	0.973913	0.974007	0.00965
0.4	0.952479	0.952789	0.03254
0.5	0.923287	0.924081	0.08592
0.6	0.884783	0.886503	0.19402
0.7	0.834560	0.837788	0.38530
0.8	0.768706	0.773603	0.63301
0.9	0.680283	0.682166	0.27603

Table 3 shows the frequency estimates of the exact, numerical and approximate analytical solutions for different combinations of the system parameters under small-, moderate-, and large-amplitude vibrations. All frequency estimates are rounded to six significant figures. In all cases considered, the numerical solution produces a near machine precision compared to the exact solution. On the other hand, the approximate analytical solution based on the frequency–amplitude formula generally produced errors that were less than 1.0% except for some few cases of moderate- and large-amplitude vibrations with $\epsilon <0$ and $\eta \gg 1.0$ . The exact results in Table 3 can be used as benchmark to determine the accuracy of other approximate analytical solutions of the system under investigation.

Table 3.

Exact ( ${\omega }_{ex}$ ), numerical ( ${\omega }_{num}$ ), and approximate ( ${\omega }_a$ ) frequency for various combinations of the system parameters for small-, moderate-, and large-amplitude vibrations.

	A	$\epsilon$	$\eta$	${\omega }_a$	${\omega }_{num}$	${\omega }_{ex}$	${ϵ}_a$  (%)	${ϵ}_{num}$  (%)
Small amplitude	0.10	0.1	0.5	1.00019	1.00019	1.00019	0.0000020469	4.44089 × 10−13
	0.10	1.0	0.5	1.00187	1.00187	1.00187	0.000201245	5.32907 × 10−13
	0.10	10.0	0.5	1.01791	1.01809	1.01809	0.0170737	4.32987 × 10−13
	0.10	−1.0	0.5	0.998116	0.998118	0.998118	0.00020901	4.996 × 10−13
	0.10	−10.0	0.5	0.980326	0.980571	0.980571	0.0249468	5.77316 × 10−13
	0.10	2.0	0.1	1.00672	1.00672	1.00672	0.0000358567	5.55112 × 10−13
	0.10	2.0	1.0	1.00000	1.00000	1.00000	0	0
	0.10	2.0	5.0	0.971692	0.971068	0.971068	0.0642532	9.21485 × 10−13
Moderate amplitude	1.0	0.1	0.5	1.01791	1.01809	1.01809	0.0170737	4.55191 × 10−13
	1.0	1.0	0.5	1.12815	1.13356	1.13356	0.477223	6.21725 × 10−13
	1.0	10.0	0.5	1.33771	1.34510	1.34510	0.548951	6.88338 × 10−13
	1.0	−0.05	0.5	0.990400	0.990456	0.990456	0.00564502	4.88498 × 10−13
	1.0	−0.5	0.5	0.877058	0.891806	0.891806	1.6537	5.32907 × 10−13
	1.0	2.0	0.1	1.47442	1.48154	1.48154	0.480714	6.21725 × 10−13
	1.0	2.0	1.0	1.00000	1.00000	1.00000	0	0
	1.0	2.0	5.0	0.542326	0.53239	0.53239	1.86633	8.21565 × 10−13
Large amplitude	10.0	0.1	0.5	1.33771	1.34510	1.34510	0.548951	7.21645 × 10−13
	10.0	1.0	0.5	1.40500	1.40616	1.40616	0.0826793	1.35447 × 10−12
	10.0	10.0	0.5	1.41327	1.41340	1.41340	0.00864372	1.06581 × 10−12
	10.0	−0.001	0.5	0.980326	0.980571	0.980571	0.0249468	4.21885 × 10−13
	10.0	−0.0065	0.5	0.823217	0.864768	0.864768	4.80488	8.41549 × 10−12
	10.0	2.0	0.1	3.07205	3.08295	3.08295	0.353585	8.99281 × 10−13
	10.0	2.0	1.0	1.00000	1.00000	1.00000	0	0
	10.0	2.0	5.0	0.448403	0.448248	0.448248	0.0344897	1.39888 × 10−12

Conclusions

The exact frequency–amplitude solution of a mass grounded by linear and nonlinear springs in series was derived in a closed-form based on the complete elliptic integral of the first and third kinds. This was achieved by transforming the exact integral of the nonlinear ODE governing the relative motion into a form that can be evaluated using standard tables of elliptic integrals. The closed-form solution showed that the response of the system depended strongly on the static nonlinearity of spring 2 ( $\epsilon$ ), the linearity factor of springs 1 and 2 ( $\eta$ ), and the vibration amplitude ( $A$ ). Periodic solutions were found to exist for all real $\epsilon$ and for $\eta >0$ . No periodic solution exists when $\eta <0$ .

The simulated results revealed the following:

a) The frequency increased with an increase in the hardening nonlinearity, that is, when $\epsilon >0$ .

b) The frequency decreased with an increase in the softening nonlinearity, that is, when $\epsilon <0$ .

c) The frequency decreased with an increase in the linearity factor, that is, when $\eta >0$ .

d) The frequency responses observed in (a) to (c) above are more pronounced at higher amplitudes.

e) The frequency is approximately constant for small amplitudes ( $A\le 0.1$ ), very large amplitudes with strong hardening nonlinearity ( $A\ge 10$ and $\epsilon \gg 1.0$ ), and for all amplitudes when $\eta \cong 1$ . Under these conditions, the system exhibits a linear or weakly nonlinear response.

f) The ratio of the maximum velocity to the maximum displacement of the system, that is, $r_{\mathit{max}}=\sqrt{1+\left(1+3\eta \right)\epsilon A^2/2+\eta {\epsilon }^2{A}^4}$ , can be used to determine the nonlinear strength of the system. For $r_{\mathit{max}}\cong 1.0$ , the system exhibits a weakly nonlinear response. For ratios far away from unity, that is, $r_{\mathit{max}}\ll 1.0$ or $r_{\mathit{max}}\gg 1.0$ , the system exhibits a strong nonlinear response. As a rule of thumb, $0.8\le r_{\mathit{max}}\le 1.2$ implies that the system’s response is weakly nonlinear and can be accurately predicted using expansion methods such as the perturbation methods or equivalent linearization methods. Otherwise, the system’s response is moderately or strongly nonlinear.

g) The system has bifurcation points at $\epsilon =0$ and $\eta =1$ . Other bifurcation points were identified and showed that the bifurcation of the system with respect of any one of the three parameters ( $\epsilon$ , $\eta$ , and $A$ ) depends on the values of the other two parameters.

h) The traditional numerical method can produce significant errors in some cases and cannot be relied upon in all cases. Therefore, the present exact solution should be used to verify the accuracy of other approximate solutions instead of numerical solutions.

i) He’s FAF is a simple method that can be used to derive fast and reasonably accurate approximate solution to the nonlinear frequency response of complex vibrating systems. However, it was observed that the approximate analytical solution produced misleading predictions of the conditions for periodic oscillation. This implies that such predictions should be subjected to further verification using more rigorous alternatives such as experimental results, exact analysis, or finite element solutions.

Finally, the present study reveals the complex nature of the oscillations of the mass grounded by linear and nonlinear springs in series. A thorough qualitative analysis is required to provide a comprehensive set of conditions under which the system experiences different responses. The qualitative analysis should include bifurcation and stability analysis. Also, the extension of the present differential-algebraic equation modelling approach to investigate other types of nonlinear springs in series will provide a basis for analysing the vibration of machine components with different combinations of elastic members.