## Abstract

A novel nonlinear dynamics model is developed in this paper to describe the static and dynamic nonlinear behaviors of a rod pendulum partially immersed in still water. The pendulum is hinged above the water level (WL) and subject to nonlinear gravity, hydrostatic, and hydrodynamic loads, all of which are incorporated into the system dynamics. The nonlinear static behavior and stability of the pendulum have been characterized by analyzing the fixed points. It is found that Pitchfork bifurcation governs the relationship between the rod density (the control parameter) and the static equilibrium angle. The pendulum's nonlinear response to external harmonic torque is obtained using harmonic balance method (HBM). The influence of system parameters, including hinge height, rod diameter, and rod density, on the nonlinear frequency response is examined. Upon altering the system parameters, particularly the rod density, it is found that the system exhibits either a softening or a hardening effect.

## 1 Introduction

Pendulums are ideal for investigating nonlinear behavior of dynamical systems because they are simple with just a few significant parameters, such as the mass and length of the rod, as well as the angle of oscillation [1]. These parameters are easily adjustable in order to investigate the impact of nonlinearity on the system's behavior. Consequently, pendulums are frequently studied to gain a deeper understanding of the nonlinear behaviors exhibited by various mechanical systems such as bifurcations, resonance, and chaos [2]. Several pendulum types have been subject to extensive study, exhibiting distinct characteristics and behaviors. These include double pendulums [3,4], inverted pendulum [5], parametric pendulums [6–8], and spring pendulums [9,10]. The latter type consists of a mass attached to a spring, which introduces additional nonlinearities into the system, leading to complex dynamics [1,11].

Extensive research has been undertaken on parametrically excited pendulums (PEP) owing to their wide spectrum of nonlinear phenomena, including both traditional periodic patterns and intricate chaotic behaviors. In their study, Leven and Koch [12] used the Melinkov and averaging approaches to calculate the limits of the subharmonic and homoclinic bifurcations, respectively. Han and Cao [13] carried out an investigation on a PEP. Their numerical results demonstrated the presence of both periodic and chaotic oscillations, as well as the occurrence of smooth and discontinuous behaviors. The dynamic response of the PEP to an external harmonic force was explored by Xu and Wiercigroch [14]. A perturbation approach was used to examine the nonlinear dynamic behavior of the system. The investigation conducted by Szemplinska-Stupmicka et al. [15,16] explored many prerequisites for the manifestation of chaos in PEPs. Furthermore, the investigators conducted an analysis on the features of global bifurcations and fractal basin attractions.

The investigation of the nonlinear dynamic characteristics exhibited by pendulums under the influence of constant and harmonic forces has attracted significant attention in many research endeavors. Han et al. [17–19] conducted a comprehensive investigation on the nonlinear dynamics of linked mass-spring oscillators and rotating pendulums, including both perturbed and unperturbed states. Various forms of bifurcation have been seen during their examination of the undisturbed system in its state of equilibrium. Dabbs and Smith [20] utilized the Melinkov method to examine the nonlinear dynamics of a rotating pendulum system under the influence of a harmonic force. The investigation has also explored the existence of two distinct horseshoe shapes. In a recent study conducted by Al-Solihat and Al Janaideh [1], the nonlinear static and dynamic characteristics of a rotating pendulum with a spinning tip mass were examined. The investigators observed that the gyroscopic moment generated by spinning tip mass could have either a stabilizing or destabilizing effect on the oscillation of the pendulum, depending on the direction and axis of the spin.

The complex geometric nonlinearity arising from the interplay between gravity and spring forces has captured the interest of many investigators on the dynamic behavior of spring pendulums. In their study, Lee and Park [21,22] utilized the multiple scales approach to examine the behavior of a harmonically excited spring-pendulum. The findings of their numerical analysis indicated the presence of Hopf bifurcation and a series of period-doubling bifurcations, finally resulting in the emergence of chaotic dynamics. Gitterman [10] performed a comparable study utilizing the multiple scales technique to estimate a solution for the nonlinear dynamic response of a spring pendulum under vertical oscillations of the hinge point. The dynamics of two planar spring pendulums connected to a horizontally stimulated platform were explored by Kapitaniak et al. [23]. Their numerical findings indicated that the pendulums display synchronized oscillation and rotating motion. Alasty and Shabani [24] conducted an investigation on the forced response of a spring pendulum system in close proximity to its static equilibrium. The effects of internal and external detuning factors, as well as various damping coefficients, on chaotic and quasi-periodic responses were examined. Amer et al. [25] examined the nonlinear dynamic characteristics of a rigid body connected to a damped spring pendulum subjected to harmonic excitation. The angular velocity of the pendulum pivot point remains constant while it traces an oval trajectory. Viet and Nghi [26] developed a nonlinear tuned mass damper utilizing a spring pendulum in order to mitigate horizontal vibration. The proposed tuned mass damper exhibited enhanced performance at elevated vibration levels due to its reduced susceptibility to parameter mistuning. Buoyancy imposes linear and rotational restoring load on a submerged rod pendulum that are analogous to those of traditional linear and torsional springs [27]. However, the nonlinear hydrostatic restoring loads characterization and analyzing their impact on the nonlinear static and dynamic behavior of pendulums have never been investigated.

The dynamics of pendulums swinging in liquids are substantially more intricate than air-swinging pendulums. This is due to the more complex interactions between the pendulum and the surrounding fluid. As a result, little research has been conducted to investigate the dynamics of pendulums subjected to fluid forces. To determine the high-amplitude oscillations of buoyant cylindrical pendulums in water, Mathai et al. [28] developed a simple theoretical model that takes into consideration the nonlinear drag, bearing friction, buoyancy, and added mass forces. The model's numerical predictions have been demonstrated to be in fair agreement with the experiments. Worf et al. [29] analyzed the dynamic behavior of a cylindrical pendulum and the interacting flow dynamics. They endeavored to determine if the two-dimensional model equations adequately represent the three-dimensional dynamics and to evaluate the emergent three-dimensional vortical flow structures. Extensive experiments have been conducted by Lenci et al. [30] to examine the response of a parametric pendulum installed on a floating support and constrained to move vertically under the influence of water waves. The experiments aimed to demonstrate the system's efficacy for harvesting energy from ocean waves. Bos and Wellenes [31] performed an experiment with a pendulum hanging just above the still water level and subjected to monochromatic free surface waves. The added mass and drag forces were included into a simplified dynamic model, which effectively corroborated the experiment findings.

In recent years, there has been a lot of focus on using pendulums to harness wave energy [32–37]. A wave energy converter (WEC) based on a pendulum was developed by Nicola et al. [33] to capture the shorter and higher frequency waves found in the Mediterranean Sea. The numerical results of the dynamic model were validated using a 1:12 scale prototype. Wu et al. [38] investigated a WEC that includes an internal inverted pendulum. The linearized equations of motion have been experimentally validated. Wang et al. [39] developed a mechanism that allows for the adjustment of the length of a pendulum installed on an unmanned marine vehicle, enabling efficient energy harvesting even in the presence of unknown excitation frequencies. A WEC based on a parametrically excited pendulum was investigated by Yurchenko and Alevras [40]. By reducing the influence of gravity, their design enhances the device's efficacy in regular sea stats. For harnessing near-shore water waves, the Oyster WEC was developed [41,42]. The WEC consists of a buoyant flap that is linked to the seabed by a hinge. The flap undergoes oscillatory motion in response to water waves.

According to our exhaustive research, the nonlinear static and dynamic behaviors of a rod pendulum partially immersed and hinged above the water level have never been explored. The primary contributions of this paper can be outlined as follows:

Develop a novel nonlinear dynamic model for a partially immersed rod pendulum in still water. The model considers the exact nonlinear hydrostatic and hydrodynamic loads while the wetted surface of the pendulum varies with time.

Explore the nonlinear static response and stability of the system. The exact hydrostatic force and moment are considered to accurately predict the large static equilibrium angle.

Analyze the nonlinear frequency response of the system to a harmonic external load. Notably, the numerical analysis incorporates the exact nonlinear terms in the equation of motion, without resorting to simplification via Taylor series, to provide precise predictions of the nonlinear dynamic characteristics of the pendulum.

The remaining content of this paper is organized as follows: Section 2 presents the dynamic model development of the pendulum and the formulations of the external loads. In Sec. 3, the nonlinear static equilibrium and stability characteristics are examined. Section 4 analyzed the nonlinear dynamic response to external harmonic torque. In Sec. 5, a parametric study has been conducted to investigate the effects of the system parameters on the nonlinear dynamic response. The major findings and concluding remarks of the paper are summarized in Sec. 6.

## 2 Dynamic Modeling and Equation of Motion

*l*, diameter

*d*, cross-sectional area $A=\pi d2/4$, and density $\rho p$. The rod is suspended from its upper end (

*O*) and undergoes oscillations with an angle of $\theta $ within the plane of the paper. The rod is partially submerged in still water (of density $\rho $). The hinge point is located at a height

*D*above the WL. The rod's center of gravity

*G*is located at its midpoint such that $rG=l/2$. The rod is subject to an external harmonic torque ($\tau e=\tau o\u2009sin(\omega t)$) applied at the hinge (

*O*), where $\tau o$ represents the magnitudes, $\omega $ is the frequency of the forcing, and

*t*denotes the time. The equation of motion for the rotation of the pendulum ($\theta $) can be expressed as

The mass of the rod is $m=\rho pAl$, and $IO=m(13l2+116d2)$ is the moment of inertia of the rod about the hinge axis. The angular velocity and acceleration of a pendulum are denoted by $\theta \u02d9$ and $\theta \xa8$, respectively. The external moments $MG$, $Mhs$, and $Mhd$ are nonlinear functions of $\theta $ and represent the gravitational, hydrostatic, and hydrodynamic moments around (*O*), respectively. The system dynamics incorporate the exact magnitudes of $MG$ and $MB$, while few assumptions are used to calculate the hydrodynamic moment ($Mhd$). These assumptions can be described as follows:

The Morison's equation is used for estimating the hydrodynamic loads exerted on the rod pendulum, which is sufficient for oscillations of slender structures in stagnant fluids [27,43–46].

Higher order hydrodynamic radiation and diffraction loads are considered negligible because of the slender shape of the rod ($d/l$ is small) [43–45]. In addition, other nonlinear hydrodynamic loads can also be ignored.

*XYZ*and the body fixed frame

*xyz*that share the origin

*O*are introduced. The gravitational moment representing the weight force moment about

*O*can be expressed as

*C*) and the rod base center (

*Q*), as shown in Fig. 2, which can be obtained from the geometry as

*B*relative to

*O*. Point

*B*represents the midsection of the submerged length, $rB$ can be then calculated as

where $Iyy=\pi d4/64$ is the second moment of the rod cross section area.

where $\omega na=3g/2l$ is the natural frequency of a rod pendulum freely hanging in air

## 3 Static Equilibrium

When the pendulum is released from specific initial conditions, it will exhibit oscillatory motion around the angle of static equilibrium. However, the amplitude of these oscillations will gradually decrease due to the presence of hydrodynamic damping.

It is essential to investigate the effects of the system parameters, which are inclusively $\alpha $, $\eta D$, and $\eta d$, on the static stability of the system to understand how the rod achieves static equilibrium.

Figure 5 shows the effect of the draft ratio ($\eta D$) on the bifurcation diagrams of the pendulum static equilibrium angle ($\theta e$), where the control parameter is $\alpha $. The results show that $\theta e$ is insensitive to the variations of $\eta D$ when $\alpha >1$ such that the pendulum will stay in vertical postilion ($\theta e=0$) if it is denser than water. However, the bifurcation point $\alpha b$ decreases as $\eta D$ increases. In general, as $\eta D$ decreases, $\theta e$ increases. This implies that a *lighter-than-water* rod will float in a larger ($\theta e$) as it becomes closer to the water level with a larger portion immersed in water (i.e., lower $\eta D$). It is also observable that the rate of variation in $\theta e$ in the vicinity of bifurcation point ($\alpha b$) for bifurcation curves associated with lower $\eta D$ is more pronounced, as depicted in Fig. 5 for $\eta D=0.25$.

The influence of the diameter ratio ($\eta d$) on the bifurcation diagrams of $\theta e$ is illustrated in Fig. 6. Multiple bifurcation diagrams were generated across a range of $\eta d$ values spanning from 0.1 to 0.5. The results indicate that an increase in $\eta d$ is associated with a slight decrease in both $\alpha b$ and $\theta e$. This is due to the fact that an increase in rod diameter results in a greater water plane restoring moment, which opposes the buoyancy moment. As a result, the rod settles at a reduced $|\theta e|$ value. It can be also observed that the impact of $\eta d$ is particularly pronounced in the vicinity of bifurcation points. The three-dimensional graph presented in Fig. 7 illustrates the relationship between the critical density ratio ($\alpha b$) and the parameters $\eta D$ and $\eta d$. As previously illustrated, an increase in $\eta D$ results in a significant decrease in $\alpha b$, whereas an increase in $\eta d$ leads to a relatively lower $\alpha b$.

## 4 Nonlinear Dynamic Behavior

where $A0$, $Ai$, and $Bi$ are the coefficients of the harmonic solution to be determined, and *n* is the number of the harmonic components. A system of nonlinear algebraic equations as a function of the harmonic coefficients can be obtained by substituting Eq. (17) into Eq. (10). The frequency response curves (amplitude of $\theta (\tau )$ versus the nondimensional frequency $\Omega $) are subsequently computed by solving the resulting equations using the arc-length continuation scheme [52,53].

To evaluate the convergence of the solution produced by the HBM (harmonic balance method), the system's frequency response to harmonic torques is calculated using varying numbers of harmonic terms, ranging from $n=2$ to $n=8$, as seen in Fig. 8. The frequency response solution clearly converges to six harmonic terms, with the solutions for $n=6$ and $n=8$ being almost similar. Therefore, 8 harmonic terms are employed for all dynamic simulations.

To verify the accuracy of the solutions obtained using the HBM method, the equation of motion is numerically integrated using the Runge–Kutta (R–K) method. A series of numerical simulations are performed to obtain the time response at different values of $\Omega $ via numerical integration. To attain a precise estimation of the steady-state amplitude, the time duration is set to 3000 cycles. Figure 9 compares the frequency response of $\theta $ using the HBM and numerical integration using Runge–Kutta (R–K) techniques. A significant agreement can be noted between the solutions acquired using the two approaches.

It is important to note that the dynamic model is applicable to pendulum oscillations when the pendulum is partially submerged in water. However, the dynamic model loses its validity when the tip of the pendulum leaves the water, which may occur if there is a significant external torque. In this instant, the hydrostatic and hydrodynamic loads are eliminated, resulting in a system dynamics similar to the traditional pendulum swinging in the air. Furthermore, as previously demonstrated, it is necessary for the pendulum to be slender in order to use Morrison's equation, which is utilized to calculate the hydrodynamic loads.

As was demonstrated earlier, the rod can reach a state of static equilibrium in either the upright posture or when it is floating in an inclined orientation. The system's stiffness is a nonlinear function of the pendulum's equilibrium angle ($\theta e$). Therefore, the nonlinear dynamic behavior of the system is dependent on its nonlinear static properties. In a similar manner, the density ratio ($\alpha $) is the primary factor that drive the nonlinear dynamic behavior of the system. In most cases, we will see three patterns of the system's nonlinear frequency responses when $\alpha >\alpha b$, $\alpha <\alpha b$, and $\alpha \u2248\alpha b$, each of which will be discussed next.

### 4.1 Nonlinear Dynamic Response at $\alpha >\alpha b$.

The nonlinear frequency response of a pendulum that is stimulated by an external harmonic torque of dimensionless amplitude $T0=0.1$ at $\alpha =4$ ($\alpha b=0.7488$) is illustrated in Fig. 10. The Numerical results are obtained by taking into account the dimensionless parameters' values, namely, $\eta D$ = 0.5, $\eta d$ = 0.1, and $Cd$ = 0.3. The frequency-maximum amplitude curve reveals, as shown in Fig. 10, that the effect of hardening is apparent. At $\Omega =0.2$, Fig. 11 depicts the corresponding time response and phase plane portrait. The system exhibits periodic response across the entire frequency spectrum. It is evident that the pendulum oscillates around $\theta e=0$ (vertical position) because the rod is significantly denser than water ($\alpha >1$). It is worth noting that $A0=0$ due to oscillation around the vertical position ($\theta e=0$).

### 4.2 Nonlinear Dynamic Response at $\alpha <\alpha b$.

The nonlinear frequency response of the pendulum, which is subjected to an external harmonic torque of dimensionless amplitude $T0=0.002$ at $\alpha =0.2$ ($\alpha b=0.7488$), is illustrated in Fig. 12. The dimensionless parameter values of $\eta D$ = 0.5, $\eta d$ = 0.1, and $Cd$ = 0.3 are utilized in the numerical analysis. In this particular case, the pendulum undergoes oscillatory motion while being inclined at a nonzero equilibrium angle ($\theta e>0$). Figure 12(a) illustrates the fluctuating response while Fig. 12(b) illustrates the mean response. The fluctuating response exhibits a peak while the mean response demonstrates a dip in the resonance zone. It is noteworthy that the frequency response curve demonstrates a softening effect, in contrast to the prior case when $\alpha >\alpha b$. The hardening effect is commonly ascribed to positive cubic stiffness, whereas the softening effect is induced by negative cubic stiffness [50]. For the system under study, the system's stiffness results from hydrostatic and gravitational restoring loads. According to the equation of motion, gravitation stiffness is positive while buoyancy stiffness is negative. The resulting hardening effect could be an indication that nonlinear stiffness of gravity is greater than that of the buoyancy.

Figure 13 depicts the time response and phase plane portrait at $\Omega =0.2$. Obviously, the pendulum oscillates about the inclined static equilibrium angle ($\theta e=55.95\xb0$). In addition, the time response amplitude is less than $3\u2009deg$, indicating that the system is highly nonlinear when the pendulum oscillates while floating at an inclined angle.

### 4.3 Nonlinear Dynamic Response at $\alpha \u2248\alpha b$.

Figure 14 illustrates the nonlinear frequency response for $\alpha \u2248\alpha b$. On the primary resonance region, the hardening effect is prominent. Notably, there are also a number of super-harmonic resonant peaks, which exhibit a softening effect unlike the primary resonance. This could be attributable to the variations in the nonlinear system stiffness during oscillation while varying the frequency as the density of the rod approaches the critical density ratio ($\alpha b$).

The time responses are depicted in Fig. 15 for various values of $\Omega $, including 0.02, 0.05, 0.1, and 1.5. When the frequency is low, the time responses will have various frequency contents, which is an indication of superharmonic resonances. However, at higher frequencies where $\Omega $ is greater than 0.12, the pendulum oscillates at a frequency identical to the forcing frequency. Figure 16 illustrates the phase plane trajectories that correspond to the time histories that were previously displayed in Fig. 15.

## 5 Parametric Analysis

As previously illustrated, the rod undergoes oscillations about the vertical position ($\theta e=0$) when its density ratio exceeds the bifurcation density ratio ($\alpha b$). However, the rod will oscillate around an inclined configuration if it is sufficiently light such that $\alpha <\alpha b$. The system parameters that affect the nonlinear behavior of the pendulum include the draft ratio ($\eta D$), the diameter ratio ($\eta d$), and the density ratio ($\alpha $), as evident from the equation of motion (Eq. (12)). This section will examine the impacts of these parameters on the nonlinear dynamic response.

The effect of the density ratio ($\alpha $) on the nonlinear frequency response of the pendulum during vertical and inclined oscillations is demonstrated in Figs. 17 and 18, respectively. It is evident that the response decreases as $\alpha $ increases for the vertical configuration, as depicted in Fig. 17. It is also important to note that the bandwidth of the resonance region increases as $\alpha $ decreases. Nevertheless, the fluctuating response of the inclined configuration is influenced by $\alpha $ in an adverse manner, whereby an increase in $\alpha $ results in an increase in the frequency response, as shown in Fig. 18(a). Furthermore, it can be observed that the mean response of the inclined configuration decreases as $\alpha $ increases, as depicted in Fig. 18(b). This finding is consistent with the static equilibrium results previously presented in Fig. 3.

The impact of the draft ratio ($\eta D$) on the nonlinear behavior of the system while oscillating in both the vertical and inclined configurations is depicted in Figs. 19 and 20, respectively. It can be observed that the increase in $\eta D$ leads to a slight decrease in the frequency response of the vertical configurations, as depicted in Fig. 19. Nevertheless, the influence of $\eta D$ is more noticeable on the inclined configuration, resulting in larger fluctuating response for higher values of $\eta D$ [Fig. 20(a)]. Conversely, the mean response decreases as $\eta D$ increases, as depicted in Fig. 20(b). It is also apparent that the effect of $\eta D$ on the mean response aligns with the static equilibrium results previously demonstrated in Fig. 5.

The effect of the diameter ratio ($\eta d$) on the nonlinear response of both the vertical and inclined configurations is illustrated in Figs. 21 and 22, respectively. In contrast to the effects of $\alpha $ and $\eta D$ previously examined, it is evident that larger magnitude of $\eta d$ leads to a proportional amplification in the magnitude of the response for both configurations. This is due to the inverse relationship between $\eta d$ and hydrodynamic damping, in which a higher value of $\eta d$ results in less damping and thus a greater response. This can be demonstrated by examining the equation of motion (Eq. (10)) in which the hydrodynamic damping is proportional to $CD/\eta d$. On the other hand, larger $\eta d$ yields less mean response for the inclined configuration, as depicted in Fig. 22(b). The reduction in the mean response is attributed to the larger water plane restoring moment. This observation is consistent with the static equilibrium findings previously presented in Fig. 6. The effect of the magnitude of the external torque magnitude ($T0$) on the nonlinear frequency response for both configurations is illustrated in Figs. 23 and 24. The results show that larger magnitude of $T0$ yields larger response for both configurations.

## 6 Conclusions

A novel nonlinear dynamic model has been developed for a partially submerged rod pendulum in stagnant water. The model takes into account the nonlinear hydrostatic, hydrodynamic, and gravitational loads. The solutions were obtained by including the exact nonlinear terms in the dimensionless equation of motion, in order to accurately capture the nonlinear static and dynamic properties of the system, particularly when subjected to significant deflections.

The study revealed that the density of the rod has a significant impact on the nonlinear static equilibrium and stability. The relationship between rod density and static equilibrium angle was characterized by a supercritical Pitchfork bifurcation, wherein the control parameter is the density ratio. It was found that rods with greater density settle to the vertical hanging position to attain static equilibrium. In the event that the density of the rod is decreased to a specific critical density ratio, the rod will float in an inclined position to achieve equilibrium. The impacts of the rod diameter and draft ratios on the nonlinear static equilibrium were also examined. The numerical results revealed that an increase in draft ratio leads to an earlier onset of bifurcation, resulting in a lower density ratio at which the bifurcation occurs. Furthermore, an increase in diameter ratio results in a decrease in the magnitude of the equilibrium angle. The study also showed that the static behavior is influenced by the diameter ratio in a manner similar to that of the draft ratio. However, the changes in the equilibrium angle resulting from variations in the diameter ratio are comparatively minor.

The harmonic balance method was used to solve the equation of motion, and the nonlinear frequency response diagrams were obtained through the utilization of the arc-length continuation method. The numerical results revealed that the frequency response exhibited a hardening effect for rods of greater density, specifically those with a density ratio surpassing the critical density ratio. However, the softening effect was noticeable when the rod density ratio was lower than the critical density ratio. In contrast, the system exhibits a pronounced hardening effect when the rod density ratio is in the vicinity of the critical density ratio. Additionally, multiple superharmonic resonance peaks become apparent. The investigation was then focused on the impact of rod's density, diameter, and draft ratios on the nonlinear frequency response of the pendulum when oscillating around the hanging and inclined positions. The parametric study revealed that an increase in the diameter ratio reduces the frequency response of the pendulum due to the reduction in the damping ratio for both configurations. Nevertheless, it has been shown that modifying the density and draft ratios of the rod results in different consequences for both hanging and inclined arrangements. Investigating the effect of the draft ratio revealed that an increase in the draft ratio results in a reduction in the frequency response for the vertical arrangement. Conversely, an adverse influence was observed for the inclined configuration. The investigation of the impact of the density ratio revealed that an increase in the density ratio results in a decrease in the frequency response of the hanging configuration; however, a contrasting effect was found for the inclined configuration.

## Acknowledgment

The author expresses gratitude to King Fahd University of Petroleum and Minerals (KFUPM) for the financial assistance.

## Funding Data

KFUPM Early Career Research Grant (Project No. EC221006; Funder ID: 10.13039/501100004055).

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

The author attests that all data for this study are included in the paper.