Stability,bifurcation,and large-amplitude vibration analysis of a symmetric magnetic spherical pendulum

Abstract

Existing studies on the symmetric spherical pendulum are limited to small- and moderate-amplitude vibrations. This study was conducted to obtain accurate solutions for analysis of the large-amplitude vibration of a symmetric magnetic spherical pendulum using the continuous piecewise linearization method (CPLM). The stability conditions and bifurcation of the pendulum were derived based on the critical points, while the CPLM was used to estimate the frequency response and vibration histories to less than 0.1% and 1.0% relative error respectively when compared to numerical solutions. The CPLM was found to be significantly more accurate than the Laplace transform homotopy perturbation method and predicted the large-amplitude bi-stable vibrations accurately. The stability analysis that was conducted enabled the characterization of all bounded symmetric vibrations based on the relationship between the cyclotron frequency and azimuthal velocity, whereas the bifurcation analysis confirmed that the symmetric vibrations can undergo pitchfork bifurcation that results in transition from single-well to double-well (or bi-stable) vibrations and vice versa. Finally, a parametric analysis was conducted to study the effect of the cyclotron frequency and uniform azimuthal velocity on the frequency–amplitude response and vibration histories The parametric analysis showed that the frequency–amplitude response has a strong dependence on the cyclotron frequency and azimuthal velocity for all amplitudes. On the other hand, the oscillation profile only depends on the cyclotron frequency and azimuthal velocity for some amplitudes. The results of this study can be applied in the design of energy harvesters and elliptic tanks for liquid transport.

Introduction

Background

Pendulum systems are one of the most investigated vibration systems. They can be used for pedagogical purposes,^1–5 understanding nonlinear phenomena,^6–8 analyzing the vibration of actual mechanical systems,^9–11 and energy harvesting applications.^12,13 The vibration of a spherical pendulum is a classical mechanics problem.^14,15 Its motion is characterized by simultaneous planar oscillations in the vertical plane and rotational motion in the horizontal plane while moving in a gravitational field.

In recent times, the vibration of the spherical pendulum in both gravitational and electromagnetic fields have been of interest due to energy harvesting potential.^12,13 The resulting vibration model derived by considering both gravitational and electromagnetic fields is somewhat complicated and difficult to solve. The reason is because the vibration model is characterized by trigonometric nonlinearities which are not directly solvable by asymptotic or series-based approximate analytical methods. The latter methods require that the trigonometric nonlinearities are expanded in Taylor's series and approximated by a few terms of the series (see, e.g. Moatimid and Amer ¹⁶ ). Since Taylor's series has a slow convergence, using a few terms in the series introduces significant truncation errors which become obvious under conditions of strong nonlinearity or large-amplitude vibration. This makes the accurate prediction of the large-amplitude vibration response challenging. There is therefore a need to derive approximate analytical solutions that are capable of estimating the large-amplitude vibration response of the spherical pendulum accurately. This article is aimed at deriving an accurate periodic solution for the large-amplitude vibration analysis of a magnetic spherical pendulum. We also determine the conditions for stable vibrations and examined the bi-stable vibrations. The results of the present study can be used in the design of energy harvesters and elliptic tanks for liquid transport.

Literature survey

Predicting the vibration response and stability of the spherical pendulum has been subject of many investigations. Landau and Lifshitz ¹⁵ formulated the complete vibration model of the spherical pendulum in polar coordinates and derived the first integral of the planar motion. They concluded that the exact periodic solution can be expressed in terms of the incomplete elliptic integral of the second kind but did not solve for it. Olsson ¹⁷ investigated the precession of a spherical pendulum restricted to small-amplitude vibrations. Based on results of perturbation analysis, he concluded that the precession rate of the elliptical orbit of the spherical pendulum is proportion to the area of the orbit. Olsson ¹⁸ conducted a follow-up study on the precession of the spherical pendulum undergoing small-amplitude vibrations. The vibration model was derived in both polar and Cartesian coordinates. It was shown that the theoretical solution for the precession orbit of the small-amplitude vibration agreed with the measurements of a simple laboratory experiment. Miles^19,20 investigated the stability and weakly nonlinear vibration of a forced damped spherical pendulum. The studies showed that the resonant response of the spherical pendulum could produce stable non-planar vibrations when the planar motion is unstable. Also, the pendulum could produce periodic doubling and chaotic response under certain parametric conditions. Miles and Zou ²¹ studied the detuning of a weakly nonlinear spherical pendulum subject to parametric excitations. They examined the bifurcation of the spherical pendulum and discovered the existence of saddle node, hopf, subcritical pitchfork, and supercritical pitchfork bifurcations.

The common feature of the older studies is that they involve a spherical pendulum moving in a gravitational field only. Recent attention has now shifted to the motion of a spherical pendulum in both gravitational and electromagnetic fields due to the potential for energy harvesting.^12,13 To this end, the magnetic spherical pendulum has become a system of interest. Cushman and Bates ²² derived a 2-DOF vibration model for a magnetic spherical pendulum with electrically charged mass moving in both magnetic and gravitational fields. They showed that the magnetic spherical pendulum can be used to study the Lagrange top with reduced symmetry about the body axis. However, they did not investigate the vibration response of the magnetic spherical pendulum. Porta and Montiel ²³ obtained analytical solutions for the free vibration of the magnetic spherical pendulum by first transforming the vibration models from polar coordinates to Cartesian coordinates and then considering only linear to weakly nonlinear conditions. Moatimid and Amer ¹⁶ applied the Laplace transform homotopy perturbation method (LTHPM) to derive approximate analytical solutions for a freely vibrating magnetic spherical pendulum with uniform horizontal velocity. They derived a 1-DOF model for the planar motions of the pendulum which they simplified further by using truncated Taylor series approximations for the trigonometric functions. They validated their approximate LTHPM solution using numerical results for small- and moderate-amplitude vibrations. They also performed a linear analysis to examine the stability of the equilibrium points. An obvious gap in the previous studies on the magnetic spherical pendulum is that they did not investigate the large-amplitude vibrations or the bi-stable vibrations of the magnetic spherical pendulum. In the present study, we attempt to close this gap and examine the conditions for the existence of symmetric vibrations through a comprehensive qualitative analysis of the magnetic spherical pendulum.

There are no exact solutions in terms of elementary functions to the nonlinear vibration models of pendulum systems. Consequently, approximate solutions are normally applied to solve nonlinear pendulum models. Currently, there are many approximate analytical methods for solving nonlinear vibration models, e.g. modified Lindstedt Poincare method, ²⁴ harmonic balance method, ²⁵ Newton harmonic balance method, ²⁶ homotopy analysis method, ²⁷ homotopy perturbation method, ²⁸ variational iteration method, ²⁹ Adomian decomposition method, ³⁰ cubication method, ³¹ quintication method, ³² energy balance method, ³³ amplitude frequency formulation, ³⁴ Mickens’ iteration method, ³⁵ global residue harmonic balance method, ³⁶ and method of multiple scales. ³⁷ Big-Alabo and Ossia ³⁸ reviewed many of these methods and highlighted a range of limitations, which include that most approximate methods:

require complex algebra for more accurate higher-order approximations,^{24,25,27–29,30,36}

produce significant errors for large-amplitude or strong nonlinear conditions,^24,25,28

are limited to symmetric oscillators,^31,32

only work with polynomial functions in the restoring force (Cheung et al., 1991; Nayfeh and Mook, 1995; Liao, 2004; He, 1999a; Adomian, 1988; Lu, 2022),^{24,25,27,28,30,36} and

are not capable of capturing essential anharmonic profile that characterizes the vibration histories of strong nonlinear vibrations.^33,34

An approximate method for the solution of nonlinear conservative oscillators that has proven to overcome the above challenges and combines accuracy with simplicity is the continuous piecewise linearization method (CPLM).^38,39 A major advantage of the CPLM is that the solution remains simple irrespective of the complexity of the restoring force. Additionally, the CPLM does not require higher-order terms, accepts arbitrary initial conditions, and works with polynomial and non-polynomial functions in the restoring force. The CPLM has been applied to solve pendulum systems undergoing strong nonlinear vibrations, ⁴⁰ slider-crank mechanism with complex non-polynomial restoring force, ⁴¹ and ship roll motion with strong asymmetric nonlinearity. ⁴² Given the above advantages of the CPLM, it was applied in the present study to investigate the frequency–amplitude response and the vibration histories of the symmetric magnetic spherical pendulum.

The literature survey shows that previous studies on the vibration of the magnetic spherical pendulum^12,16,21,22 are limited to small- to moderate-amplitude vibrations and have not investigated its bi-stable vibration response. Therefore, this study focuses on the stability, bifurcation, and large-amplitude vibration of a freely vibrating symmetric magnetic spherical pendulum using the CPLM. The conditions for stability were explored more comprehensively compared to the study conducted by Moatimid and Amer, ¹⁶ and a bifurcation analysis was conducted based on the critical points. In particular, the conditions for double-well or bi-stable vibrations, not considered in other studies, were established, and an accurate periodic solution for the bi-stable vibration histories was obtained.

Scope and contribution of study

The scope of the present study is limited to the theoretical analysis of the stability, bifurcation, and large-amplitude vibration analysis of a magnetic spherical pendulum constrained to a uniform azimuthal velocity. The stability and bifurcation of the pendulum were examined through qualitative analysis of the critical points. The equation to determine the bifurcation points was derived based on the parameters of the pendulum, and all bounded symmetric vibrations were characterized. The latter showed the conditions under which bi-stable symmetric vibrations can occur. For the periodic solution, the CPLM algorithm was applied to derive accurate solutions for the frequency response and vibration histories of the large-amplitude vibration under uni- and bi-stable conditions. The CPLM was shown to be more accurate than published approximate results. Finally, parametric analysis was conducted to examine the effect of the azimuthal velocity and cyclotron frequency on the vibration response of the pendulum.

Formulation of vibration model

The magnetic spherical pendulum consists of an electrically charged bob of mass, m, carrying a positive charge, e, and suspended from a fixed point by a weightless inextensible string of length, l. The motion of the bob is influenced by a constant gravitational field, g, and a constant magnetic field, B. The magnetic field is induced by a magnetic monopole at the equilibrium point of the bob. A diagrammatic representation of the physical model of the magnetic spherical pendulum is shown in Figure 1(a). During motion the string of the pendulum forms an angle $\theta$ with the $z$ -axis, while the pendulum bob sweeps through an horizontal angle $\varphi$ in the $x-y$ plane, also known as the azimuthal angle. The equation of motion for the magnetic spherical pendulum in Figure 1(a) has been derived by Moatimid and Amer, ¹⁶ but we present the derivation here because our total potential energy was derived differently and for completeness.

Figure 1.

(a) Physical model of the magnetic spherical pendulum. (b). Travel path of the planar motion of the pendulum from its equilibrium position.

The pendulum possesses kinetic and potential energies due to the gravitational and magnetic fields. First, we consider the kinetic and potential energies due to the gravitational field. Given that the motion of the pendulum is in 3-D, its position at any instance can be defined by the position vector:

$$\stackrel{\rightharpoonup }{r}=ai+bj-ck$$

(1)where i, j, and k are the unit vectors in the x, y, and z directions, respectively, and a, b, and c are geometric distances derived from Figure 1(a) as $a=l\sin \theta \cos \varphi$ , $b=l\sin \theta \sin \varphi$ , and $c=l\cos \theta$ . Therefore,

$$\stackrel{\rightharpoonup }{r}=(l\sin \theta \cos \varphi )i+(l\sin \theta \sin \varphi )j-(l\cos \theta )k$$

(2)The velocity vector is:

$$\begin{array}{r}\stackrel{\rightharpoonup }{v}=\frac{d\stackrel{\rightharpoonup }{r}}{dt}=l(\dot{\theta }\cos \theta \cos \varphi -\dot{\varphi }\sin \theta \sin \varphi )i+l(\dot{\theta }\cos \theta \sin \varphi +\dot{\varphi }\sin \theta \cos \varphi )j+(l\dot{\theta }\sin \theta )k\end{array}$$

(3)from which the magnitude of the velocity was obtained as follows:

$$v^2=|\stackrel{\rightharpoonup }{v}{|}^2=l^2[{\dot{\theta }}^2+{\dot{\varphi }}^2{\sin }^2\theta ]$$

(4)Then, the kinetic energy due to the motion of the pendulum is:

$$K{E}_g={\displaystyle \frac{1}{2}m{l}^2[{\dot{\theta }}^2+{\dot{\varphi }}^2{\sin }^2\theta ]}$$

(5)The gravitational potential energy of the pendulum depends on the elevation of the bob from its equilibrium position (see Figure 1(b)) and is given as:

$$P{E}_g=mgh=mgl(1-\cos \theta )$$

(6)Next, we determine the energy of the pendulum due to the magnetic field. The only energy introduced by the influence of the magnetic field is the magnetic potential energy. This energy is a scalar product of the Lorentz force and the position vector of the pendulum. In the absence of an electric field, the Lorentz force on the charged pendulum bob due to the influence of a magnetic field is $e(\stackrel{\rightharpoonup }{v}\times \stackrel{\rightharpoonup }{B})$ . Therefore, the magnetic potential energy can be expressed as:

$$P{E}_m=e(\stackrel{\rightharpoonup }{v}\times \stackrel{\rightharpoonup }{B})\cdot \stackrel{\rightharpoonup }{r}$$

(7)where $\stackrel{\rightharpoonup }{B}$ is a uniform magnetic field vector. Assuming the magnetic field acts only in the positive z direction, then $\stackrel{\rightharpoonup }{B}=Bk$ where B is the magnitude of the uniform magnetic field. After evaluating the vector and scalar products we get:

$$P{E}_m=eB{l}^2\dot{\varphi }{\sin }^2\theta$$

(8)Hence, the total potential and kinetic energies of the pendulum are respectively:

$$U=P{E}_g+P{E}_m=mgl(1-\cos \theta )+eB{l}^2\dot{\varphi }{\sin }^2\theta$$

(9)and

$$T=K{E}_g={\displaystyle \frac{1}{2}m{l}^2[{\dot{\theta }}^2+{\dot{\varphi }}^2{\sin }^2\theta ]}$$

(10)Consequently, the Lagrangian of the magnetic spherical pendulum is:

$$L=T-U={\displaystyle \frac{1}{2}m{l}^2[{\dot{\theta }}^2+{\dot{\varphi }}^2{\sin }^2\theta ]-mgl(1-\cos \theta )-eB{l}^2\dot{\varphi }{\sin }^2\theta }$$

(11)which can be rewritten as:

$$L=m{l}^2\left[{\displaystyle \frac{{\dot{\theta }}^2}{2}+{\sin }^2\theta \left({\displaystyle \frac{{\dot{\varphi }}^2}{2}-\dot{\varphi }{\omega }_c}\right)}\right]-mgl(1-\cos \theta )$$

(12)where ${\omega }_c=eB/m$ is the cyclotron frequency or gyrofrequency in SI unit. Note that in the Gauss unit system the cyclotron frequency must be multiplied by the Lorentz factor of $1/c$ where c is the velocity of light.

Now, the equations of motion of the freely vibrating undamped spherical magnetic pendulum can be derived from the Euler–Lagrangian equation given as:

$${\displaystyle \frac{d}{dt}\left({\displaystyle \frac{\mathrm{\partial }L}{\mathrm{\partial }\dot{q}}}\right)-{\displaystyle \frac{\mathrm{\partial }L}{\mathrm{\partial }q}=0}}$$

(13)where q is a generalized coordinate representing any one of the two coordinates, $\theta$ and $\varphi$ , that define the position of the pendulum.

For the $\theta$ -coordinate,

$$\ddot{\theta }+{\displaystyle \frac{g}{l}\sin \theta -({\dot{\varphi }}^2-2\dot{\varphi }{\omega }_c)\cos \theta \sin \theta =0}$$

(14)For the $\varphi$ -coordinate,

$${\displaystyle \frac{d}{dt}(m{l}^2{\sin }^2\theta (\dot{\varphi }-{\omega }_c))=0}$$

(15)which implies that $m{l}^2{\sin }^2\theta (\dot{\varphi }-{\omega }_c)$ is equal to a constant, say $\gamma$ , that represents the horizontal component of the angular momentum (i.e. the azimuthal momentum). Thus,

$$\dot{\varphi }={\displaystyle \frac{\gamma }{m{l}^2{\sin }^2\theta }+{\omega }_c}$$

(16)From equation (16), it is obvious that the azimuthal velocity generally changes with time in a manner that balances the time variation of the vertical angle such that the azimuthal momentum is conserved. Since the azimuthal momentum is constant, one can express it as the product of a uniform azimuthal velocity and a fixed inertia. The assumption of a uniform azimuthal velocity introduces symmetry in the vibrations of the spherical pendulum. ²² Therefore, by taking the uniform azimuthal velocity to be $\dot{\varphi }=\mu$ , we can establish from equation (14) that the vibration of the symmetric magnetic spherical pendulum is governed by:

$$\ddot{\theta }+{\omega }_0^2\sin \theta -{\displaystyle \frac{1}{2}\mathrm{\Gamma }\sin 2\theta =0}$$

(17)where ${\omega }_0=\sqrt{g/l}$ and $\mathrm{\Gamma }={\mu }^2-2{\omega }_c\mu$ with $\{\mu ,{\omega }_c\}\ge 0$ . If $\mu =0$ then equation (17) reduces to the classical simple pendulum oscillating in a gravitational field, whereas $\mu \ne 0$ and ${\omega }_c=0$ represents a conical pendulum oscillating in a gravitational field or a pendulum with rotating support where $\mu$ becomes the dimensionless rotational velocity. Hence, equation (17) can be used to investigate the nonlinear vibrations of the Watt governor. ⁴³ In a situation where $\mu =2{\omega }_c$ , the horizontal momentum of the spherical pendulum cancels the effect of the magnetic field so that the system oscillates as a simple pendulum in a gravitational field.

Equation (17) is the model for the planar symmetric vibration of the unforced undamped magnetic spherical pendulum, and it can be written as:

$$\ddot{\theta }+f(\theta )=0$$

(18)where $f(\theta )$ is the conservative restoring force of the spherical magnetic pendulum subject to uniform azimuthal velocity and is given as:

$$f(\theta )=\sin \theta -{\displaystyle \frac{1}{2}\mathrm{\Gamma }\sin 2\theta }$$

(19)Note that in equation (19) the assumption, ${\omega }_0=1.0$ , was applied without loss of generality of the analysis.

Integrating equation (17) and applying the initial conditions, $\theta (0)=0$ and $\dot{\theta }(0)=v_0$ , we get:

$$2({\dot{\theta }}^2-v_0^2)=4(\cos \theta -1)+\mathrm{\Gamma }(1-\cos 2\theta )$$

(20)In order to implement the CPLM, knowledge of the restoring force as well as the range of oscillation is required. The restoring force can be obtained from equation (19) whereas the oscillation range is from $-A$ to A, where A is the oscillation amplitude which occurs at the point where $\dot{\theta }=0$ . Hence, the amplitude of oscillation for a given initial velocity and zero initial displacement was obtained from equation (20) as:

$$A=\{\begin{array}{c}{\tan }^{-1}\left({\displaystyle \frac{\sqrt{(2-v_0^2)\mathrm{\Gamma }-2+2\sqrt{1+(v_0^2-2)\mathrm{\Gamma }+{\mathrm{\Gamma }}^2}}}{1-\sqrt{1+(v_0^2-2)\mathrm{\Gamma }+{\mathrm{\Gamma }}^2}}}\right)\mathrm{if}\mathrm{\Gamma }\le 2-v_0^2\\ \pi +{\tan }^{-1}\left({\displaystyle \frac{\sqrt{(2-v_0^2)\mathrm{\Gamma }-2+2\sqrt{1+(v_0^2-2)\mathrm{\Gamma }+{\mathrm{\Gamma }}^2}}}{1-\sqrt{1+(v_0^2-2)\mathrm{\Gamma }+{\mathrm{\Gamma }}^2}}}\right)\mathrm{if}\mathrm{\Gamma }>2-v_0^2\end{array}$$

where $\mathrm{\Gamma }>0$ . Note that in arriving at equation (21), the trigonometric identity $\cos A=1/\sqrt{1+{\tan }^{2}A}$ was applied to equation (20) and the resulting expression in terms of $\tan A$ was solved to get equation (21). This rather algebraically complex approach was used because ${\tan }^{-1}$ has an unrestricted input domain, i.e. $-\mathrm{\infty }\le x\le \mathrm{\infty }$ , unlike ${\cos }^{-1}$ and ${\sin }^{-1}$ that are limited to input domain of $-1<x<1$ . If $\mathrm{\Gamma }=2-v_0^2$ , then the amplitude of vibration is exactly $90\circ$ or $\pi /2$ radians. This implies that the equation $\mu ={\omega }_c+\sqrt{{\omega }_c^2+2-v_0^2}$ maps the initial and parametric conditions for which the vibration amplitude of the spherical magnetic pendulum is exactly $90\circ$ . Figure 2 illustrates a map of the conditions that produce a $90\circ$ vibration amplitude.

Figure 2.

Mapping of initial and parametric conditions for $90\circ$ vibration amplitude. Left, 3-D map, and right, 2-D map.

Qualitative analysis of the vibration response

Potential function

The potential function is an analytical expression for the potential energy (energy possessed by a body due to its state or position) of a system. The potential function represents the energy stored in the system when it is displaced from its equilibrium position. Understanding the potential function is essential for analyzing the behavior of oscillating systems as it provides insight into the stability of the system and the energy changes within the system as it oscillates. The potential function, $V(\theta )$ , is defined as the integral of the restoring force. Therefore,

$$V(\theta )={\displaystyle \frac{\mathrm{\Gamma }}{4}\cos 2\theta -\cos \theta }$$

(22)The only parameter on which the potential function depends on is $\mathrm{\Gamma }$ . To demonstrate the effect of $\mathrm{\Gamma }$ on the potential function, we investigated three cases as presented in Figures 3 to 5. The figures show the potential function plot as well as the corresponding phase plots obtained using equation (20). On the phase plots, the label “ $c$ ” represents the center of oscillation while “ $s$ ” represents a saddle node where periodic oscillations cease to exist. The magnetic spherical pendulum exhibits a single-well oscillation with a center at the origin when $\mathrm{\Gamma }<0$ and $\mathrm{\Gamma }=0$ , whereas it can produce a double-well oscillation with two centers that are equidistant from the origin when $\mathrm{\Gamma }>0$ . The pendulum does not produce any symmetric single-well oscillation with the center outside the origin. For the cases of $\mathrm{\Gamma }=0$ and $\mathrm{\Gamma }>0$ , the vibration amplitude falls in the range $0\circ <A<180\circ$ , and for the case of $\mathrm{\Gamma }<0$ , the vibrations may be restricted to the range $0\circ <A<A_L$ where $A_L=180\circ -{\tan }^{-1}\left(\sqrt{{\mathrm{\Gamma }}^2-1}\right)>90\circ$ is the limiting amplitude at which point and beyond which there is no vibration. The formula for $A_L$ points to the fact that $\mathrm{\Gamma }<-1.0$ . This implies that the following inequalities must be satisfied: ${\omega }_c-\sqrt{{\omega }_c^2-1}<\mu <{\omega }_c+\sqrt{{\omega }_c^2-1}$ and ${\omega }_c\ge 1.0$ .

Figure 3.

Potential function and phase plots for $\mathrm{\Gamma }=0$ ( $\mu =2.0\mathrm{rad}/\mathrm{s}$ , and ${\omega }_c=1.0\mathrm{rad}/\mathrm{s}$ ).

Figure 4.

Potential function and phase plots for $\mathrm{\Gamma }<0$ ( $\mu =2.0\mathrm{rad}/\mathrm{s}$ , and ${\omega }_c=2.0\mathrm{rad}/\mathrm{s}$ ).

Figure 5.

Potential function and phase plots for $\mathrm{\Gamma }>0$ ( $\mu =3.0\mathrm{rad}/\mathrm{s}$ , and ${\omega }_c=1.0\mathrm{rad}/\mathrm{s}$ ).

Stability analysis

Two system parameters, the azimuthal velocity and the cyclotron frequency, were varied to investigate the system's qualitative response. If the system returns to the equilibrium position with time, the system is said to possess a stable equilibrium. On the other hand, if the system deviates from its equilibrium position or produces an unbounded response after a small displacement, it is said to possess unstable equilibrium. The stable and unstable points can be obtained from the potential function. These points occur when the potential gradient is zero and are known as critical points ( ${\theta }_c$ ). Therefore,

$${\displaystyle \frac{dV}{d\theta }=\sin {\theta }_c(1-\mathrm{\Gamma }\cos {\theta }_c)=0}$$

(23)We deduce from equation (23) that the critical points occur at ${\theta }_c=0$ and ${\theta }_c=\pm {\tan }^{-1}\left(\sqrt{{\mathrm{\Gamma }}^2-1}\right)$ . The latter case suggests that non-zero real solutions of the critical points can only be obtained when $|\mathrm{\Gamma }|>1.0$ . The implication is that the condition for the existence of bi-stable oscillations is $\mu >{\omega }_c+\sqrt{{\omega }_c^2+1}>2{\omega }_c$ (see Figure 5). Furthermore, the first solution of the critical points shows that a center of oscillation exists at the origin while the second solution shows the existence of two centers of oscillation equidistant from the origin. A possible inference would be the existence of bifurcations from single-well to double-well oscillations and vice versa. This possibility was examined in section “Bifurcation analysis.”

A vibrating system is stable at minimum potential energy and unstable at maximum potential energy. This implies that we can find the stable and unstable points of the critical points by taking the second derivative of the potential function. Therefore, the conditions for stable and unstable points are, respectively,

$${{\displaystyle \frac{d^{2}V}{d{\theta }^2}}|}_{\theta ={\theta }_c}>0\mathrm{and}{{\displaystyle \frac{d^{2}V}{d{\theta }^2}}|}_{\theta ={\theta }_c}<0$$

(24)These conditions lead to:

$$\mathrm{\Gamma }>{\displaystyle \frac{\cos {\theta }_c}{2{\cos }^2{\theta }_c-1}\mathrm{and}\mathrm{\Gamma }}>{\displaystyle \frac{\cos {\theta }_c}{2{\cos }^2{\theta }_c-1}}$$

for stable and unstable points, respectively.

The stability of the pendulum based on the critical points has been illustrated graphically in Figure 6. The plots show the stable points in blue color and the unstable points in red color. The two points, X ( $v_0=1.85rad/s,{\omega }_c=0.5rad/s$ and $\mu =0.5rad/s$ ) and Y ( $v_0=1rad/s,{\omega }_c=2.3rad/s$ and $\mu =5rad/s$ ), shown on the graph of stable points represent distinct qualitative responses of the magnetic spherical pendulum. Point X represents one of the centers of a double-well oscillation, while point Y represents the center of a single-well oscillation. On the other hand, the graph of the unstable points shows that vibrations are only possible when $A<180\circ$ . Furthermore, it can be observed that a limiting amplitude in the range $90\circ <A_L<180$ exist under certain conditions, whereas there are no periodic oscillations around the origin under other conditions. These observations correspond to the cases $\mathrm{\Gamma }<-1.0$ and $\mathrm{\Gamma }>1.0$ , respectively. Instances of other parametric conditions of $\mathrm{\Gamma }$ that produce bounded symmetric oscillations and the corresponding stability conditions are summarized in Table 1.

Figure 6.

3-D stability plots showing stable points (left) and unstable points (right).

Table 1.

Characterization of the vibration of the magnetic spherical pendulum.

	Cases	Parametric condition	Stability condition
1	$\mu >{\omega }_c+\sqrt{{\omega }_c^2+1}$ for ${\omega }_c\ge 0$	$\mathrm{\Gamma }>1.0$	Double-well with two centers equidistant from the origin
2	$\mu ={\omega }_c+\sqrt{{\omega }_c^2+1}$ for ${\omega }_c\ge 0$	$\mathrm{\Gamma }=1.0$	Single-well with a center at the origin
3	${\omega }_c+\sqrt{{\omega }_c^2-1}<\mu <{\omega }_c+\sqrt{{\omega }_c^2+1}$ for ${\omega }_c>1.0$	$-1.0<\mathrm{\Gamma }<1.0$ ,	Single-well with a center at the origin
4	$\mu ={\omega }_c+\sqrt{{\omega }_c^2-1}$ for ${\omega }_c>1.0$	$\mathrm{\Gamma }=-1.0$	Single-well with a center at the origin
5	${\omega }_c-\sqrt{{\omega }_c^2-1}<\mu <{\omega }_c+\sqrt{{\omega }_c^2-1}$ for ${\omega }_c>1.0$	$\mathrm{\Gamma }<-1.0$	Single-well with a center at the origin and limiting amplitude, $90\circ <A_L<180$
6	$\mu ={\omega }_c-\sqrt{{\omega }_c^2-1}$ for ${\omega }_c>1.0$	$\mathrm{\Gamma }=-1.0$	Single-well with a center at the origin
7	$0\le \mu <{\omega }_c-\sqrt{{\omega }_c^2-1}$ for ${\omega }_c>1.0$	$-1.0<\mathrm{\Gamma }\le 0$	Single-well with a center at the origin

Bifurcation analysis

Bifurcation analysis involves systematically varying parameters in a dynamical system to explore how its response or trajectories change qualitatively. It offers insight into the stability, patterns, and responses of the dynamical system to varying conditions. The numerical values of a set of system parameters that produces a qualitative change in the system's response is known as the bifurcation point. In the present study, a bifurcation analysis was conducted by considering the qualitative change in the stability conditions of the system. The two system parameters, i.e., ${\omega }_c$ and $\mu$ , were considered in determining the bifurcation points. The condition for the occurrence of bifurcation is:

$${{\displaystyle \frac{d^{2}V}{d{\theta }^2}}|}_{\theta ={\theta }_c}=\cos {\theta }_c-\mathrm{\Gamma }({\cos }^2{\theta }_c-{\sin }^2{\theta }_c)=0$$

(25)The bifurcation points occur on the ${\theta }_c=0$ line. Therefore, equation (25) simplifies to $\mathrm{\Gamma }=1$ , from which we derive the bifurcation equation:

$$\mu ={\omega }_c+\sqrt{{\omega }_c^2+1}$$

(26)Bifurcation diagrams showing the effects of the uniform azimuthal velocity and cyclotron frequency on the stability of the magnetic spherical pendulum are plotted in Figures 7 and 8. The stable points are plotted with solid red lines, and the unstable points are plotted with dashed red lines. The magnetic spherical pendulum exhibits a pitchfork bifurcation for any given value of ${\omega }_c\ge 0$ . The bifurcation is called a pitchfork bifurcation because the diagram of the critical points looks like a pitchfork. For ${\omega }_c\ge 0$ , the bifurcation results in transition from single-well oscillations to double-well oscillations. In contrast, there is no bifurcation when $\mu =1.0$ . However, for values of $\mu >1.0$ , a pitchfork bifurcation occurs such that the system transits from a bi-stable or double-well oscillations to a single-well oscillation. The bifurcation points are determined from equation (26) and are plotted in Figure 9. Above the curve of the bifurcation points, the pendulum experiences double-well oscillations. Below the curve, single-well oscillations occur with a center at the origin. For the single-well oscillations, a limiting amplitude that is less than 180° occurs if ${\omega }_c-\sqrt{{\omega }_c^2-1}\le \mu \le {\omega }_c+\sqrt{{\omega }_c^2-1}$ ; otherwise the vibration amplitude is limited by the geometrical constraint imposed by the pendulum, i.e. $0<A<180\circ$ . The line $\mu =0$ represents the vibrations of a magnetic pendulum without rotation, whereas the origin of the plot ( $\mu =0$ , ${\omega }_c=0$ ), which is also below the curve, represents the vibrations of a simple pendulum. Both of these conditions produce single-well oscillations that are only limited by the geometrical constraint of the pendulum.

Figure 7.

Bifurcation plots showing the effect of the uniform azimuthal velocity on the stability of the spherical magnetic pendulum.

Figure 8.

Bifurcation plots showing the effect of the cyclotron frequency on the stability of the magnetic spherical pendulum.

Figure 9.

Bifurcation points of the magnetic spherical pendulum.

Periodic solution

Continuous piecewise linearization method

The CPLM is based on piecewise discretization of the force-displacement or compliance response of a vibrating system. Each discretization of the nonlinear compliance response is approximated using a linear compliance from which a linearized version of the equation of motion is obtained. The solution of the linearized equation of motion provides an approximate solution for the nonlinear equation of motion over the range of discretization. Hence, the CPLM was formulated such that it updates the solution of the linearized equation as we move from one discretization to the next. This way, an approximate solution for any nonlinear conservative system can be obtained in an unconditionally stable scheme. ³⁸ Apparently, an extension of the CPLM to solve nonlinear vibration models of non-conservative systems requires numerical solution to the nonlinear algebraic equation resulting from determining the time. ⁴⁴ Hence, it is most effective for conservative systems where the CPLM uses only closed-form solutions.

For conservative systems with symmetric vibrations, the range of discretization is a quarter cycle (i.e. $0<\theta <A$ ), and for asymmetric vibrations, the discretization is done over a half cycle, i.e. ( $B<\theta <A$ ) where A and B are the amplitudes on either side of the vibration cycle and $A\ne B$ . Traditionally, the discretization procedure is based on n equal divisions of the solution domain. In this case, the step size of each discretization is constant and given as $\mathrm{\Delta }\theta =A/n$ for symmetric oscillators and $\mathrm{\Delta }\theta =|A-B|/n$ for asymmetric oscillators. A constant step size converges fast when the nonlinearity of the compliance is generally uniform within the solution domain. Under these circumstances, a maximum error less than 0.1% could be achieved with $n\le 100$ . On the other hand, in situations where there are wide variations in the nonlinearity of the compliance response, as in very stiff systems, the convergence of the CPLM with constant step size is very slow and require a very large n, e.g. 1000, to achieve the needed accuracy. Consequently, the computational time to achieve the needed accuracy rises considerably.

In order to accelerate the convergence and reduce the computational time, unequal divisions can be applied to discretize the solution domain so that the regions of the compliance response that are characterized by strong nonlinear response can use small step sizes, while the regions with weak to moderate nonlinearity can use large step sizes (see Figure 10). One such technique that utilizes the concept of unequal divisions has been proposed by Big-Alabo and Ezekwem, ⁴⁵ who applied the CPLM to investigate the nonlinear transient response of the temperature of a combustible iron particle. The discretization technique is such that the step size reduces with successive discretization as we move away from the origin to the limit of the solution domain. This is because the nonlinearity becomes stronger as we approach the limit of the solution domain. Based on this technique, the step size for the $i$ th discretization of a symmetric oscillator can be expressed as ⁴⁵ :

$$\mathrm{\Delta }\theta (i)=\mathrm{\Delta }{\theta }_i={\displaystyle \frac{2(n-i+1)}{n(n+1)}A,1\le i\le n}$$

(27a)when the discretization is performed from zero displacement to the amplitude. In the opposite direction, the step size becomes:

$$\mathrm{\Delta }{\theta }_i={\displaystyle \frac{2i}{n(n+1)}A,1\le i\le n}$$

(27b)

Figure 10.

Discretization of nonlinear compliance of a symmetric conservative oscillator (a) equal (constant) step size (b) unequal step size.

According to the CPLM algorithm, the linearized force within each discretization of the compliance response could be expressed as ³⁸ :

$$f_{rs}(\theta )=\pm |K_{rs}|({\theta }_{}-{\theta }_r)+F_r$$

(28)where r and $s=r+1$ denote the initial and the end points of each discretization respectively, $F_r=f({\theta }_r)$ and $K_{rs}$ is the linearized stiffness given as:

$$K_{rs}={\displaystyle \frac{f({\theta }_s)-f({\theta }_r)}{{\theta }_s-{\theta }_r}}$$

(29)Substituting equation (27) in equation (18) gives the linearized equation of motion for each discretization as:

$$\ddot{\theta }\pm |K_{rs}|\theta =\pm |K_{rs}|{\theta }_r-F_r$$

(30)Equation (30) is a nonhomogeneous second-order linear differential equation and its solution, which is straightforward, depends on whether the sign of $|K_{rs}|$ is positive or negative.

Positive linearized stiffness, ${\mathit{K}}_{\mathit{r}\mathit{s}}>0$

The solution of equation (30) if $K_{rs}>0$ is given as:

$$\theta (t)=R_{rs}\sin ({\omega }_{rs}t+{\phi }_{rs})+C_{rs}$$

(31)where ${\omega }_{rs}=\sqrt{|K_{rs}|}$ ; $C_{rs}={\theta }_r-F_r/|K_{rs}|$ ; and $R_{rs}=\sqrt{{({\theta }_r-C_{rs})}^2+{({\dot{\theta }}_r/{\omega }_{rs})}^2}$ . The initial conditions and other parameters are determined based on the oscillation stage. For the half cycle that moves from $+A$ to $\text{–}A$ the initial conditions are ${\theta }_r={\theta }_{r-1}-\mathrm{\Delta }\theta$ and ${\dot{\theta }}_r=-\sqrt{\left|2\underset{A}{\overset{{\theta }_r}{\int }}-f(\theta )\right|}$ , and the other solution parameters for this half cycle are:

$${\phi }_{rs}=\{\begin{array}{c}\pi +{\tan }^{-1}\left[{\displaystyle \frac{{\omega }_{rs}({\theta }_r-C_{rs})}{{\dot{\theta }}_r}}\right],if{\dot{\theta }}_r<0\\ \\ \\ \pi /2,if{\dot{\theta }}_r=0\end{array}$$

(32)

$$\mathrm{\Delta }t=\{\begin{array}{c}{\displaystyle \frac{1}{{\omega }_{rs}}\left[0.5\pi +{\cos }^{-1}\left({\displaystyle \frac{{\theta }_s-C_{rs}}{R_{rs}}}\right)-{\phi }_{rs}\right],if({\theta }_s-C_{rs})<R_{rs}}\\ \\ \\ {\displaystyle \frac{1}{{\omega }_{rs}}(0.5\pi -{\phi }_{rs}),if({\theta }_s-C_{rs})\ge R_{rs}}\end{array}$$

(33)where $\mathrm{\Delta }t$ is the time range for each discretization. The time at the end of each discretization is $t_s=t_r+\mathrm{\Delta }t$ and the end conditions ${\theta }_s$ and ${\dot{\theta }}_s$ are obtained by replacing r with s in the initial condition formulas.

Essentially, the CPLM uses the discretization and linearization of the nonlinear compliance response to derive corresponding time solutions for points on the exact velocity-displacement or phase relationship. Therefore, the CPLM generates the solution ( $t_s$ , ${\theta }_s$ , ${\dot{\theta }}_s$ ) at the endpoint of each discretization and this can be plotted to produce the displacement and velocity histories. In this sense, the CPLM can be implemented as a point-wise solution. However, it should be noted that the CPLM solution for each discretization is applicable over a time range and not just a point. More details about the point-wise and range-wise features of the CPLM and its unique attributes that distinguishes it from traditional numerical schemes have been discussed in ref. ³⁸

On the other hand, the initial conditions for the half cycle that moves from $-A$ to $+A$ are ${\theta }_r={\theta }_{r-1}+\mathrm{\Delta }\theta$ and ${\dot{\theta }}_r=\sqrt{\left|2\underset{A}{\overset{{\theta }_r}{\int }}-f(\theta )\right|}$ , and the other solution parameters for this half cycle are:

$${\phi }_{rs}=\{\begin{array}{c}{\tan }^{-1}\left[{\displaystyle \frac{{\omega }_{rs}({\theta }_r-C_{rs})}{{\dot{\theta }}_r}}\right],if{\dot{\theta }}_r>0\\ \\ \\ -\pi /2,if{\dot{\theta }}_r=0\end{array}$$

(34)

$$\mathrm{\Delta }t=\{\begin{array}{c}{\displaystyle \frac{1}{{\omega }_{rs}}\left[0.5\pi -{\cos }^{-1}\left({\displaystyle \frac{{\theta }_s-C_{rs}}{R_{rs}}}\right)-{\phi }_{rs}\right],if({\theta }_s-C_{rs})<R_{rs}}\\ \\ \\ {\displaystyle \frac{1}{{\omega }_{rs}}(0.5\pi -{\phi }_{rs}),if({\theta }_s-C_{rs})\ge R_{rs}}\end{array}$$

(35)Negative linearized stiffness, ${\mathit{K}}_{\mathit{r}\mathit{s}}<0$

The solution of equation (30) if $K_{rs}<0$ is given as:

$$\theta (t)=A_{rs}{e}^{{\omega }_{rs}t}+B_{rs}{e}^{-{\omega }_{rs}t}+C_{rs}$$

(36)where ${\omega }_{rs}=\sqrt{|K_{rs}|}$ ; $C_{rs}={\theta }_r+F_r/|K_{rs}|$ ; $A_{rs}$ and $B_{rs}$ are integration constants that were derived based on the initial conditions (same as in $K_{rs}>0$ ) as:

$$A_{rs}={\displaystyle \frac{1}{2}\left({\theta }_r+{\displaystyle \frac{{\dot{\theta }}_r}{{\omega }_{rs}}-C_{rs}}\right)}$$

(37a)

$$B_{rs}={\displaystyle \frac{1}{2}\left({\theta }_r-{\displaystyle \frac{{\dot{\theta }}_r}{{\omega }_{rs}}-C_{rs}}\right)}$$

(37b)The time range for each discretization was obtained by applying the end conditions to the solution. Thus,

$$\mathrm{\Delta }t=\{\begin{array}{c}{\displaystyle \frac{1}{{\omega }_{rs}}In\left[{\displaystyle \frac{({\theta }_s-C_{rs})+\mathrm{sgn}({\dot{\theta }}_r)\sqrt{{({\theta }_s-C_{rs})}^2-4A_{rs}{B}_{rs}}}{2A_{rs}}}\right],{({\theta }_s-C_{rs})}^2>4A_{rs}{B}_{rs}}\\ \\ \\ {\displaystyle \frac{1}{{\omega }_{rs}}In\left[{\displaystyle \frac{({\theta }_s-C_{rs})}{2A_{rs}}}\right],{({\theta }_s-C_{rs})}^2\le 4A_{rs}{B}_{rs}}\end{array}$$

(38)where $\mathrm{sgn}(x)$ is the signum function. Lastly, the nonlinear frequency can be determined as:

$${\omega }_{CPLM}={\displaystyle \frac{\pi }{2{\sum }_{i=1}^n\mathrm{\Delta }{t}_i}}$$

(39)

Application of CPLM to the magnetic spherical pendulum

During implementation of the CPLM, we need to know the conservative restoring force and the solution domain. The conservative restoring force can be obtained from the equation of motion (see equation (19)), and the solution domain for the magnetic spherical pendulum is that of a symmetric oscillator. The solution domain is used for discretization while the restoring force is used to determine the main CPLM parameters. For the magnetic spherical pendulum, the main CPLM parameters are:

$$F_r=f({\theta }_r)=\sin {\theta }_r-{\displaystyle \frac{1}{2}\mathrm{\Gamma }\sin 2{\theta }_r}$$

(40)

$$K_{rs}={\displaystyle \frac{2(\sin {\theta }_s-\sin {\theta }_r)-\mathrm{\Gamma }(\sin 2{\theta }_s-\sin 2{\theta }_r)}{2({\theta }_s-{\theta }_r)}}$$

(41)

$${\dot{\theta }}_r=\pm \sqrt{2(\cos {\theta }_r-\cos A)-{\displaystyle \frac{1}{2}\mathrm{\Gamma }(\cos 2{\theta }_r-\cos 2A)}}$$

(42)The remaining CPLM parameters are calculated from these three parameters.

The magnetic spherical pendulum has the following initial conditions: $\theta (0)=0$ and $\dot{\theta }(0)=v_0$ . This means that the CPLM will be initially applied for the solution domain $0<\theta <A$ where A is given in equation (21). Afterward, it can be applied to $A>\theta >-A$ and vice versa. The CPLM procedure as explained in section “Continuous piecewise linearization method” and eqs. (40)–(42) were used to develop a bespoke computer program in Mathematica™ that was used to implement and simulate the CPLM solution for the magnetic spherical pendulum.

Convergence and accuracy of the CPLM

In studying the vibration of the magnetic spherical pendulum, we consider the following response categories: small-amplitude ( $0\circ <A<15\circ$ ), moderate-amplitude ( $15\circ \le A<60\circ$ ), large-amplitude ( $60\circ \le A<150\circ$ ), and very large-amplitude ( $150\circ \le A<180\circ$ ) vibrations. The reference solution used for the convergence and accuracy tests is the numerical solution of equation (17) obtained using the NDSolve function in Mathematica™. The convergence test was conducted to ascertain the number of discretization required to achieve a pre-set relative error (RE) of 0.1% in estimating the time period. The relative error was calculated as $ϵ=100\text{\%}\times |1-T_{CPLM}/T_{num}|$ .

The convergence test was conducted using equal and unequal discretization for large-amplitude ( $A=110\circ$ , $\mu =0.25$ rad/s and ${\omega }_c=5.0$ rad/s) and very large-amplitude ( $A=175\circ$ , $\mu =0.12$ rad/s and ${\omega }_c=5.0$ rad/s) vibrations. Note that for each amplitude that was considered, the values of $\mu$ and ${\omega }_c$ used here and hereafter were selected to satisfy the conditions in Table 1. The numerical result for the time period of the large-amplitude vibration is $T_{num}=14.7264\mathrm{s}$ and that of the very large-amplitude vibration is $T_{num}=10.9361\mathrm{s}$ . The results of the convergence test give the RE of the CPLM in estimating $T_{num}$ for each n and are shown in Table 2 and Figure 11. It is obvious that the CPLM converges unconditionally as n increases, and the convergence accelerates rapidly when unequal discretization was used. Typically, the CPLM was found to converge to less than 0.1% RE for the unequal discretization when $n=20$ , whereas the same convergence was achieved for the equal discretization when $n=160$ . Accelerating the CPLM convergence using unequal discretization saves computational time and significantly improves the accuracy of the solution for a given number of discretization. Hence, $n\le 20$ and the unequal discretization were used to generate all other CPLM results.

Figure 11.

Comparison of convergence of equal and unequal discretization: (a) large-amplitude vibration ( $A=110\circ$ ) and (b) very large-amplitude vibration ( $A=175\circ$ ).

Table 2.

Convergence test for equal and unequal discretization.

Large-amplitude			Very large-amplitude
$\mathit{n}$	Equal	Unequal	Equal	Unequal
10	8.07158	0.482816	10.2196	0.498242
12	5.90034	0.28061	7.07642	0.288204
15	4.0154	0.145166	4.60942	0.148586
18	2.94126	0.0849147	3.29525	0.0867523
20	2.46202	0.06233257	2.72854	0.0636112
25	1.69893	0.0324217	1.84964	0.0330065
35	0.983906	0.0121893	1.05268	0.012312
50	0.559683	0.00443733	0.591656	0.00438039
70	0.332006	0.00182934	0.348239	0.00170969
90	0.225782	0.00103012	0.235768	0.000890577
110	0.166288	0.00070711	0.173127	0.000559348
130	0.129053	0.000552361	0.134067	0.000400586
150	0.103934	0.000469141	0.107787	0.000315172
160	0.0942848	0.000441753	0.0977075	0.000287074

The accuracy test was conducted to determine the error of the CPLM solution in predicting the displacement and velocity histories using the number of discretization obtained from the convergence test. Also, the CPLM solution was compared with the solution of the second-order LTHPM derived by Moatimid and Amer. ¹⁶ The second-order LTHPM solution was derived for the initial conditions: $\theta (0)=0$ , and $\dot{\theta }(0)=1.0rad/s$ . The results of the accuracy test are shown in Figures 12 to 14. The CPLM solution performs significantly better than the second-order LTHPM as shown in Figure 12. The LTHPM solution was found to deviate particularly as the displacement approaches the amplitude where the restoring force is increasingly nonlinear. The maximum RE of the CPLM was 0.329542% whereas the maximum RE of the second-order LTHPM was 75.8537%. Furthermore, the numerical solution of the frequency was 0.999445 rad/s, while the frequency estimates of the CPLM and LTHPM were 0.999493 rad/s (i.e. 0.00488336% RE) and 1.008790 rad/s (i.e. 0.935059% RE), respectively. The vibration histories at points Y and X of Figure 6 and the error of the CPLM were plotted in Figures 13 and 14. The conditions at these points represent large-amplitude vibrations with strong nonlinear response. The CPLM result matches the numerical result and captures the anharmonic vibration profile excellently. In Figure 14, the bi-stable vibration produces a velocity-history with two global maximum points and two global minimum points where the kinetic energy of the pendulum is a maximum. It can be observed that the displacement at the maximum kinetic energy is not zero due to the initial conditions used. The maximum RE of the CPLM in Figures 13 and 14 was found to be 0.747266%. The accuracy test shows that CPLM solution with unequal discretization and $n=20$ produces very accurate results for all possible bounded periodic oscillations of the magnetic spherical pendulum.

Figure 12.

Comparison of CPLM ( $n=20$ ) and LTHPM solutions of the displacement and velocity histories to corresponding numerical solutions ( $\mu =0.10\mathrm{rad}/\mathrm{s}$ and ${\omega }_c=1.0\mathrm{rad}/\mathrm{s}$ ).

Figure 13.

Comparison of CPLM ( $n=20$ ) to numerical solutions for conditions at point Y of Figure 6.

Figure 14.

Comparison of CPLM ( $n=20$ ) to numerical solutions for bi-stable conditions at point X of Figure 6.

Frequency–amplitude response and vibration histories

The periodic response, in terms of the frequency–amplitude response and vibration histories, of the magnetic spherical pendulum was investigated to ascertain the effects of the cyclotron frequency and uniform azimuthal velocity. The results of the investigation are plotted in Figures 15 to 30 for small-amplitude, moderate-amplitude, large-amplitude, and very large-amplitude vibrations. The plots were generate for $n\le 20$ and the maximum relative error of the CPLM in predicting the frequency–amplitude response was less than $0.1\text{\%}.$

Figure 15.

(a) Effect of uniform azimuthal velocity on the frequency–amplitude response for small-amplitude vibrations: $1\circ \le \theta \le 15\circ$ (b) CPLM error plot ( ${\omega }_c=5.0\mathrm{rad}/\mathrm{s}$ ).

Figure 16.

Effect of uniform azimuthal velocity ( $\mu$ ) on the displacement-history (left plot) and velocity-history (right plot) of a small-amplitude vibration; $A=10\circ$ ( ${\omega }_c=5.0\mathrm{rad}/\mathrm{s}$ ).

Figure 17.

(a) Effect of $\mu$ on the frequency–amplitude response at moderate-amplitude vibrations: $15\circ \le \theta \le 60\circ$ (b) CPLM error plot ( ${\omega }_c=5.0\mathrm{rad}/\mathrm{s}$ ).

Figure 18.

Effect of uniform azimuthal velocity ( $\mu$ ) on the displacement-history (left plot) and velocity-history (right plot) of a moderate-amplitude vibration; $A=40\circ$ ( ${\omega }_c=5.0\mathrm{rad}/\mathrm{s}$ ).

Figure 19.

(a) Effect of $\mu$ on the frequency–amplitude response at large-amplitude vibrations: $60\circ \le \theta \le 150\circ$ (b) CPLM error plot ( ${\omega }_c=5.0\mathrm{rad}/\mathrm{s}$ ).

Figure 20.

Effect of uniform azimuthal velocity ( $\mu$ ) on the displacement-history (left plot) and velocity-history (right plot) of a large-amplitude vibration; $A=110\circ$ ( ${\omega }_c=5.0\mathrm{rad}/\mathrm{s}$ ).

Figure 21.

(a) Effect of $\mu$ on the frequency–amplitude response at very large-amplitude vibrations: $150\circ \le \theta \le 179\circ$ (b) CPLM error plot ( ${\omega }_c=5.0rad/s$ ).

Figure 22.

Effect of uniform azimuthal velocity ( $\mu$ ) on the displacement-history (left plot) and velocity-history (right plot) of a very large-amplitude vibration; $A=172\circ$ ( ${\omega }_c=5.0\mathrm{rad}/\mathrm{s}$ ).

Figure 23.

(a) Effect of ${\omega }_c$ on the on the frequency–amplitude response for small-amplitude vibrations: $1\circ \le \theta \le 15\circ$ (b) corresponding error plot ( $\mu =0.3\mathrm{rad}/\mathrm{s}$ ).

Figure 24.

(a) Effect of ${\omega }_c$ on the displacement-history (left plot) and velocity-history (right plot) of a small-amplitude vibration; $A=10\circ$ ( $\mu =0.3rad/s$ ).

Figure 25.

(a) Effect of ${\omega }_c$ on the frequency–amplitude response at moderate-amplitude vibrations: $15\circ \le \theta \le 60\circ$ (b) corresponding error plot ( $\mu =0.3\mathrm{rad}/\mathrm{s}$ ).

Figure 26.

(a) Effect of ${\omega }_c$ on the displacement-history (left plot) and velocity-history (right plot) of a moderate-amplitude vibration; $A=40\circ$ ( $\mu =0.3\mathrm{rad}/\mathrm{s}$ ).

Figure 27.

(a) Effect of ${\omega }_c$ on the frequency–amplitude response at large-amplitude vibrations: $60\circ \le \theta \le 150\circ$ (b) corresponding error plot ( $\mu =0.3\mathrm{rad}/\mathrm{s}$ ).

Figure 28.

(a) Effect of ${\omega }_c$ on the displacement-history (left plot) and velocity-history (right plot) of a large-amplitude vibration; $A=110\circ$ ( $\mu =0.3\mathrm{rad}/\mathrm{s}$ ).

Figure 29.

(a) Effect of ${\omega }_c$ on the frequency–amplitude response at very large-amplitude vibrations: $150\circ \le \theta \le 179\circ$ (b) corresponding error plot ( $\mu =0.3\mathrm{rad}/\mathrm{s}$ ).

Figure 30.

(a) Effect of ${\omega }_c$ on the displacement-history (left plot) and velocity-history (right plot) of a very large-amplitude vibration; $\mathrm{A}=172\circ$ ( $\mu =0.3\mathrm{rad}/\mathrm{s}$ ).

For the small-amplitude vibration, an increase in the cyclotron frequency or azimuthal velocity results in a stronger nonlinear stiffness with increased nonlinear frequency. However, the amplitude has no significant effect on the nonlinear frequency. During moderate-amplitude vibration, the nonlinear frequency increases with the cyclotron frequency or azimuthal velocity and decreases with increase in the amplitude. For the large-amplitude vibration, the nonlinear frequency increases with the cyclotron frequency or azimuthal velocity, but it decreases with an increase in the amplitude up to a point beyond which the nonlinear frequency increases with amplitude. This qualitative change in the frequency–amplitude response is indicative of the existence of bifurcation due to variation in the amplitude. For the very large-amplitude vibration, the nonlinear frequency increases with the cyclotron frequency or azimuthal velocity. Again, the nonlinear frequency increases slightly with amplitude for changing azimuthal velocity (see Figure 21). In the case of very large-amplitude vibrations with changing cyclotron frequency (see Figure 29), the nonlinear frequency decreases slightly with the increase in amplitude up to certain value of the cyclotron frequency beyond which it increases slightly with increase in amplitude. This behavior is indicative of the existence of a bifurcation due to changes in the cyclotron frequency. Bifurcation driven by the cyclotron frequency was illustrated in Figures 8 and 9 and has been discussed in section “Convergence and accuracy of the CPLM.”

The plots of the vibration histories (see Figures 16, 18, 20, 22, 24, 26, 28, and 30) show the effects of the azimuthal velocity and cyclotron frequency on the displacement and velocity profiles of the magnetic spherical pendulum for small-amplitude, moderate-amplitude, large-amplitude, and very large-amplitude vibrations. In all cases, except for the large-amplitude vibration, the azimuthal velocity and cyclotron frequency showed no significant effect on the displacement and velocity profiles. The effects of varying the azimuthal velocity and cyclotron frequency only show up on the response period (or frequency). However, for the large-amplitude vibration, it was observed that an increase in either of the azimuthal velocity or cyclotron frequency leads to stronger anharmonic profile, thereby confirming a stronger nonlinear response. We can then conclude that the azimuthal velocity and cyclotron frequency can be used as tuning parameters to effectively adjust the nonlinear response of the magnetic spherical pendulum.

Conclusion

The oscillation of a spherical pendulum is both gravitational and electromagnetic field is of interest because of its potential for energy harvesting. Previous studies on the symmetric magnetic spherical pendulum have been limited to small- and moderate-amplitude vibrations or weakly nonlinear conditions. In this article, we derived accurate periodic solutions to investigate the large-amplitude vibrations of a symmetric magnetic spherical pendulum. The vibration model was derived as a single degree-of-freedom symmetric oscillator in terms of the planar motion, and it reveals that the symmetric vibration of the magnetic spherical pendulum depends on the uniform azimuthal velocity and the cyclotron frequency, the latter representing the effect of the magnetic field on the charged pendulum mass. The restoring force of the magnetic spherical pendulum is a trigonometric function of sine and cosine terms which makes its dynamics complex and difficult to solve accurately with approximate analytical techniques; especially during large-amplitude vibrations when the pendulum possess very strong nonlinearity. In this study, the CPLM was applied to obtain accurate solutions for the large-amplitude vibrations of the magnetic spherical pendulum. Using an improved discretization technique that accelerates the convergence of the CPLM, it was observed that frequency solutions with less than 0.1% relative error and vibration histories with less than 1.0% relative error can be achieved with a few discretization of the solution domain. For the magnetic spherical pendulum, $n\le 20$ , was sufficient to achieve the required level of accuracy for all possible vibrations, including large-amplitude and very large-amplitude vibrations. However, the CPLM results can be significantly improved by using $n>20$ (see Table 2). The CPLM solution based on $n=20$ was several orders more accurate than the LTHPM.

Stability analysis of the magnetic spherical pendulum revealed the existence of single-well and double-well oscillations. The center of vibration for all symmetric single-well oscillations is at the origin or equilibrium point of the pendulum, while the double-well oscillations have their centers at two points equidistant from the origin. For the single-well oscillations, the stability analysis showed the existence of vibrations constrained by a limiting amplitude that is less than the kinematic constraint imposed by the pendulum. On the other hand, the bifurcation analysis showed that the pendulum exhibits a pitchfork bifurcation from single- to double-well oscillations and vice versa. The type of symmetric vibration experienced by the pendulum depends on the relative magnitudes of the uniform azimuthal velocity and cyclotron frequency. Therefore, the conditions for the different types of symmetric vibration were obtained and summarized in Table 1. Furthermore, an equation to determine the bifurcation points was derived and plotted to characterize the regions for double-well and single-well oscillations.

Finally, a parametric analysis was conducted to study the effect of the cyclotron frequency and uniform azimuthal velocity on the frequency–amplitude response and vibration histories. The results of the parametric studies (Figures 15 to 30) revealed the following:

The frequency increases with the cyclotron frequency and uniform azimuthal velocity for small-amplitude, moderate-amplitude, large-amplitude, and very large-amplitude vibrations.

The frequency shows no significant change with amplitude during small-amplitude vibrations, whereas it decreases as amplitude increases during moderate-amplitude vibrations.

For large-amplitude vibrations, the frequency decrease with increasing amplitude until it gets to an amplitude beyond which the frequency increases with amplitude. This implies a bifurcation due to change in amplitude. For the very large-amplitude vibration, the frequency increases only slightly with the increase in amplitude.

The oscillation profiles showed no observable change with increase in cyclotron frequency or uniform azimuthal velocity during small-amplitude, moderate-amplitude, and very large-amplitude vibrations. However, a stronger anharmonic profile was observed in the velocity-history of the large-amplitude vibrations when the cyclotron frequency or uniform azimuthal velocity was increased.

The results of the present study can be applied in the design of energy harvesting systems and elliptic tanks for liquid transport. The present study was limited to symmetric vibrations of the magnetic spherical pendulum. However, the pendulum can also undergo asymmetric vibrations when the azimuthal velocity is not constant. This can be the subject of future investigations to deepen understanding of the dynamics of the magnetic spherical pendulum. Furthermore, experimental investigations can be conducted to provide validation for the theoretical analysis and deepen understanding of the system's energy harvesting potential.