Dynamics of heavy and buoyant underwater pendulums

Abstract

The humble pendulum is often invoked as the archetype of a simple, gravity driven, oscillator. Under ideal circumstances, the oscillation frequency of the pendulum is independent of its mass and swing amplitude. However, in most real-world situations, the dynamics of pendulums is not quite so simple, particularly with additional interactions between the pendulum and a surrounding fluid. Here we extend the realm of pendulum studies to include large amplitude oscillations of heavy and buoyant pendulums in a fluid. We performed experiments with massive and hollow cylindrical pendulums in water, and constructed a simple model that takes the buoyancy, added mass, fluid (nonlinear) drag and bearing friction into account. To first order, the model predicts the oscillation frequencies, peak decelerations and damping rate well. An interesting effect of the nonlinear drag captured well by the model is that, for heavy pendulums, the damping time shows a non-monotonic dependence on pendulum mass, reaching a minimum when the pendulum mass density is nearly twice that of the fluid. Small deviations from the model’s predictions are seen, particularly in the second and subsequent maxima of oscillations. Using time-resolved particle image velocimetry (TR-PIV), we reveal that these deviations likely arise due to the disturbed flow created by the pendulum at earlier times. The mean wake velocity obtained from PIV is used to model an extra drag term due to incoming wake flow. The revised model significantly improves the predictions for the second and subsequent oscillations.

1 Introduction

The first known study on pendulums dates back to Galileo Galilei in 1605. The story goes that Galileo observed the lamplighter pushing one of the swaying chandeliers in the Pisa cathedral. Galileo timed the swings with his pulse and concluded that, although the amplitude decreased, the time of each swing remained constant. Half a century later this observation inspired Christiaan Huygens to invent the pendulum clock, which until 1930 has set the standard for accurate timekeeping. Based on Newton’s laws of motion Huygens could derive that an ideal, frictionless pendulum has a period $T=2\unicode[STIX]{x03C0}\sqrt{L/g}$ for small amplitudes, where $L$ is its length, and $g$ the acceleration due to gravity (Huygens 1986). However, as is well known to any clockmaker, real pendulums are often far from ideal. Even in vacuum, and when friction is negligible, one has to take into account the finite amplitude and mass distribution of the pendulum in order to correctly predict its period. When a pendulum swings in a fluid such as air or water, additional fluid forces have to be taken into account. In fact, it was through careful observations of deviations from ideal pendulum motion that the nature and strength of many of these fluid forces were brought to light. Concepts like buoyancy, added mass and fluid friction (viscosity) emerged concurrently with the study of pendulums swinging in different fluid environments (Stokes 1851).

In recent years there has been a renewed interest in the use of simple pendulums to probe fluid–structure interactions. Neill, Livelybrooks & Donnelly (2007) conducted experiments on steel and brass pendulums oscillating in different fluids. Their results agreed fairly well with the motion predicted by taking the buoyancy and added mass into account, implying that viscous corrections played only a minor role in the oscillation frequency. In an extension to this, Bolster, Hershberger & Donnelly (2010) conducted experiments on a pendulum oscillating in a viscous fluid. In all of the above studies, the researchers focused on small amplitude vibrations of dense pendulums in viscous fluids, where the dynamics could be reasonably modelled using the ideal flow approximation, with the inclusion of a linear drag in the equation of pendulum motion. Many studies have focused on the flow disturbances induced by an oscillating body in a fluid. For one-dimensional oscillations at low to moderate Reynolds numbers, the effect of a body’s oscillation on the wake around it has been studied in some detail (Tatsuno & Bearman 1990; Tatsuno & Karasudani 1993). Various types of wakes have been quantified depending on the body oscillation frequency, velocity, fluid viscosity and the characteristic length scale of the problem. More recent studies have addressed the influence of flow disturbances on pendula stability. Obligado, Puy & Bourgoin (2013) found that the equilibrium of a pendulum facing an incoming flow displays both bi-stability and hysteresis. Further, flow disturbances due to incoming turbulence can modify the stability map by changing the drag acting on the pendulum.

Figure 1. Schematics of the experimental set-up for (a) heavy and (b) buoyant pendulums. (c) Cylinder mass was varied by inserting brass disks into the hollow cylinder. (d) Experimental set-up used for performing particle image velocimetry (PIV) measurements.

Here we extend the scope of pendulum–fluid interaction studies to the regime of large amplitude motions of heavy and buoyant cylindrical pendulums in a dense fluid. We have performed controlled experiments in which we simultaneously track the motion of the pendulum and the flow field around it, in a regime where the Reynolds number (based on the cylinder diameter) $Re\sim \mathit{O}(10^{4})$ . The control parameter in our experiment is the cylinder to fluid density ratio $m^{\ast }\equiv \unicode[STIX]{x1D70C}_{p}/\unicode[STIX]{x1D70C}_{f}$ . The ratio of the cylinder diameter to the pendulum length $D/L$ , the fluid properties ( $\unicode[STIX]{x1D70C}_{f}$ and $\unicode[STIX]{x1D708}$ ) and the release angle ( $\unicode[STIX]{x1D703}_{0}=90^{\circ }$ ) are fixed. We vary $m^{\ast }$ , covering both buoyant ( $m^{\ast }<1$ ) and heavy ( $m^{\ast }>1$ ) pendulums. We build on the knowledge of buoyancy, added mass and inertial drag to construct a basic model for the pendulum motion, which we compare in detail to the experiments. We find that a reduced added mass as compared to the ideal flow case leads to better agreement between experiment and model predictions. We use insights from particle image velocimetry (PIV) measurements to model a higher-order fluid–structure interaction, by which initially induced fluid motions affect the force on the pendulum at later instants.

2 Experiments

The experiments were conducted using a cylindrical pendulum, with length $L=203~\text{mm}$ , diameter $D=32~\text{mm}$ and span $W=300~\text{mm}$ . This resulted in a cylinder span to diameter ratio $W/D\approx 10$ . The cylinder attached to the pendulum arm was constructed from polyvinyl chloride (PVC) pipe sealed with two plastic caps. Masses could be inserted inside the pipe (see figure 1 c), which let us study both heavy and buoyant pendulums in water.

The pendulum was placed in a water tank ( $800~\text{mm}\times 400~\text{mm}\times 500~\text{mm}$ ), and pivoted on a metal rod (radius $R=3~\text{mm}$ ), which was fixed to the walls using suction cups. The wire connecting the cylinder to the pivot axis had a diameter of 1.2 mm. This was thick enough to not bend during oscillations, while also providing low resistance to the flow. The cylinder ends remained at least $1.3D$ away from the side walls. The lowest position of the heavy cylinder during its motion was $5D$ above the base of the water tank, and similarly, the highest position of the buoyant cylinder during its motion was $5D$ below the water level. Schematics of the basic experimental arrangement for the heavy and buoyant cases are shown in figures 1(a) and 1(b), respectively.

The release mechanism consisted of a semicircular cavity with two curved metal wire segments. The cylinder was held against the semicircular cavity, and locked in place by the metal wires. The release was done by pulling the wire, which ensured a gentle release of the pendulum without creating disturbances in the water. By inverting these parts, the buoyant pendulums could also be released without disturbing the water. The release position was at a $90^{\circ }$ angle from the vertical. A waiting interval of 15 min ensured that any residual flow from a previous run had damped out.

A high-speed camera (Photron Fastcam 1024PCI) was used to record the pendulum motion. The cylinder movement was detected using a circle detection method (Mathai et al. 2015). The standard deviation of error in centre detection was approximately 0.27 mm for all cases, i.e. ${<}1\,\%$ of the pendulum diameter. Additionally, we performed time-resolved particle image velocimetry (TR-PIV) to measure the flow field surrounding the pendulum. The flow was seeded with fluorescent polystyrene tracer particles (diameter ${\approx}125~\unicode[STIX]{x03BC}\text{m}$ ). The seeding particle size was small enough to have both a small Stokes number and a small Stokes/Froude ratio (Mathai et al. 2016). A high-speed Laser (Litron LDY-303HE) along with cylindrical optics was used to create a light sheet through the mid-section of the cylinder. Mirrors were placed as shown in figure 1(d), which ensured that the shadows cast by the cylinder were removed by reflected light. This enabled us to have well-lit tracer particles all around the cylinder. A double-frame camera (Imager sCMOS) was used at a frame rate of 25 frames per second. The particle seeding density for PIV was such that an interrogation window contained around 6 particles, which is considered good for accuracy. The inter-pulse time $\unicode[STIX]{x0394}t$ was varied from 1 to 3 ms. $\unicode[STIX]{x0394}t$ was optimised by ensuring that most of the tracer particles remained within the interrogation window, with a maximum movement ${\approx}1/4$ th of the window width. The PIV analysis was done with two stage processing and 50 % overlap using LaVision software.

3 Model equation of motion

The equation of motion of a pendulum oscillating in a fluid can be expressed by Newton’s law in the angular form:

(3.1) $$\begin{eqnarray}I\frac{\text{d}^{2}\unicode[STIX]{x1D703}}{\text{d}t^{2}}=\sum \unicode[STIX]{x1D70F}.\end{eqnarray}$$

Here, $\unicode[STIX]{x1D703}$ is the angular position of the pendulum with respect to the equilibrium position, $I$ the moment of inertia of the system and $\unicode[STIX]{x1D70F}$ the net torque acting on the system (see figure 2 a). The pendulum mass and moment of inertia come mostly from the cylinder at a distance $L$ from the rotation axis. Since $2L/D>10$ , the cylinder moment of inertia about its mass centre $I_{cm}=mD^{2}/8$ is negligible compared to the pendulum moment of inertia $mL^{2}$ . The effect of the pendulum accelerating through the fluid is modelled by means of an added mass $m_{a}$ . Therefore, the effective mass of the system can be written as $m_{eff}=m+m_{a}$ . With these approximations, we model the system as a point mass with a moment of inertia: $I_{p}=m_{eff}L^{2}$ .

The torques acting on the system can be written as $\sum \unicode[STIX]{x1D70F}=\sum r\times F$ , where $F$ are the individual forces, and $r$ the corresponding arm lengths. The forces acting on the system are gravity $F_{g}\equiv \unicode[STIX]{x1D70C}_{p}{\mathcal{V}}g$ , buoyancy $F_{B}\equiv \unicode[STIX]{x1D70C}_{f}{\mathcal{V}}g$ , fluid drag $F_{D}$ and a bearing friction $F_{f}\approx \unicode[STIX]{x1D707}_{f}F_{N}$ (see figure 2 a). Here, $\unicode[STIX]{x1D70C}_{p}$ and $\unicode[STIX]{x1D70C}_{f}$ are the particle and fluid densities, respectively, ${\mathcal{V}}$ is the cylinder volume, $\unicode[STIX]{x1D707}_{f}$ is the friction coefficient and $F_{N}\approx (F_{B}-F_{g})\cos \unicode[STIX]{x1D703}$ is the normal force acting at the bearing. The tension due to the centrifugal force is negligible. The Galileo number $Ga\equiv \sqrt{gD^{3}|m^{\ast }-1|}/\unicode[STIX]{x1D708}$ . The velocity scale in $Ga$ is a gravitational velocity $v_{g}=\sqrt{gD|1-m^{\ast }|}$ . The Reynolds number is defined as $Re=v_{max}D/\unicode[STIX]{x1D708}$ , where $v_{max}$ is the measured maximum velocity of the cylinder in the first swing. We note that $Ga\sim \mathit{O}(10^{3}{-}10^{4})$ gives a predictive estimate of the $Re$ range in the experiment (see table 1). The drag force may be written as $F_{D}=(1/2)\unicode[STIX]{x1D70C}_{f}A_{p}C_{D}v_{p}^{2}$ , where $C_{D}$ is the drag coefficient, $A_{p}$ is the projected area of the cylinder and $v_{p}\equiv L\text{d}\unicode[STIX]{x1D703}/\text{d}t$ is the instantaneous cylinder velocity (Batchelor 2000). By defining $m^{\ast }\equiv \unicode[STIX]{x1D70C}_{p}/\unicode[STIX]{x1D70C}_{f}$ and $\tilde{t}\equiv t\sqrt{g/L}$ , the equation of motion can be written in non-dimensional form as

(3.2) $$\begin{eqnarray}m_{eff}^{\ast }\frac{\text{d}^{2}\unicode[STIX]{x1D703}}{\text{d}\tilde{t}^{2}}=-k\sin \unicode[STIX]{x1D703}-c\left|\frac{\text{d}\unicode[STIX]{x1D703}}{\text{d}\tilde{t}}\right|\frac{\text{d}\unicode[STIX]{x1D703}}{\text{d}\tilde{t}}-h|\text{cos}\,\unicode[STIX]{x1D703}|\text{sgn}\left(\frac{\text{d}\unicode[STIX]{x1D703}}{\text{d}\tilde{t}}\right),\end{eqnarray}$$

where $m_{eff}^{\ast }=(m^{\ast }+m_{a}^{\ast })$ ; $k=|m^{\ast }-1|$ ; $c=(1/2)C_{D}A_{p}L/{\mathcal{V}}$ ; and $h=\unicode[STIX]{x1D707}_{f}|m^{\ast }-1|R/L$ .

Figure 2. (a) Schematic of the cylindrical pendulum at an instant when it swings to the right. Blue arrows indicate the forces on the pendulum: buoyancy $F_{B}$ , weight $F_{g}$ , drag $F_{D}$ , normal force due to the pendulum arm $F_{N}$ and a bearing friction $F_{f}$ . The dashed arrow and $v_{p}$ denote the instantaneous velocity of the cylinder. (b) Angular position ( $\unicode[STIX]{x1D703}$ ) versus time for $m^{\ast }=4.98$ and $m^{\ast }=1.16$ cases. (p1) and (p2) – peak deceleration points.

Table 1. $Ga$ and $Re$ for a selected number of mass-density ratios $m^{\ast }$ used in the experiment.

In order to solve (3.2), $m_{a}^{\ast }$ , $C_{D}$ and $\unicode[STIX]{x1D707}_{f}$ have to be set. It is well known that added mass coefficients $m_{a}^{\ast }$ for oscillating cylinders depend on a variety of parameters, including the frequency of oscillation, distance to boundaries, free surfaces, etc. Most studies have focused on the case of cylinders near a free surface and for small oscillations (Dong 1978; Tatsuno & Bearman 1990; Konstantinidis 2013; Koo & Kim 2015; Mathai et al. 2017, 2018). However, in complex situations with relative motions, curved trajectories and unsteady three-dimensional wakes with flow separation, $m_{a}^{\ast }$ can deviate significantly from the two-dimensional potential flow added mass coefficient. Moreover, the fact that the cylinder span is finite will induce a three-dimensional (3-D) flow near the cylinder ends, allowing some fluid to move to the sides. Therefore the cylinder motion is expected to accelerate less fluid in the 3-D case as compared to 2-D potential flow. Nevertheless, to begin with, we choose the 2-D potential flow added mass coefficient of an infinite cylinder, $m_{a}^{\ast }=1.0$ .

The drag coefficient range in our experiments could be estimated from the range of $Re$ (see table 1). We begin with a simplified mean drag coefficient $C_{D}\approx 1.2$ based on the Reynolds number range in our experiments (Lienhard 1966). Note that the influence of Reynolds number and vortex shedding on $C_{D}$ are not included at this stage, and will be discussed later. $\unicode[STIX]{x1D707}_{f}$ is expected to lie in the range [0.2–0.3] for lubricated steel-on-steel contact. We estimated $\unicode[STIX]{x1D707}_{f}\approx 0.2$ , based on tests performed using very heavy pendulums in air, since in these cases the drag and added mass forces were at a minimum.

Figure 2(b) shows $\unicode[STIX]{x1D703}$ versus $t$ predicted by the model for heavy pendulums with $m^{\ast }=1.16$ and $m^{\ast }=4.98$ . Clearly, the amplitude decay and oscillation frequency depend on the mass density ratio. With this in mind, we explore the predictions for a range of $m^{\ast }$ , covering heavy and buoyant underwater pendulums.

Figure 3. (a) Contour plot of $\unicode[STIX]{x1D703}$ versus time $t$ and mass-density ratio $m^{\ast }$ . (b) Phase portraits from numerical solution of (3.2) showing the evolution of angular velocity $\unicode[STIX]{x1D714}_{pend}$ versus angular position $\unicode[STIX]{x1D703}$ for $m^{\ast }=0.02$ and $m^{\ast }=1.98$ cases. Note that these cases were chosen with identical driving: $|F_{B}-F_{g}|$ . Phase portraits obtained for available experimental data (not shown here) show similar behaviour.

4 Results

4.1 Model predictions

We first present some typical pendulum motions predicted by the model equation of motion. In figure 3(a), we show a contour plot of the angular position $\unicode[STIX]{x1D703}$ versus time $t$ for continuous variation of $m^{\ast }$ . We notice a clear asymmetry about $m^{\ast }=1$ line, which separates the heavy cases from the buoyant cases. For instance, the oscillation frequency is slightly higher for a buoyant cylinder as compared to a heavy cylinder that has identical driving $|F_{B}-F_{g}|$ . On the buoyant side ( $m^{\ast }<1$ ), the oscillations damp faster with decreasing $m^{\ast }$ , owing to the decreasing inertia. Other aspects of the pendulum motion are altered in going from $m^{\ast }>1$ to $m^{\ast }<1$ . In figure 3(b), we show the phase portraits of two cases, one buoyant ( $m^{\ast }=0.02$ ) and other heavy ( $m^{\ast }=1.98$ ), with identical driving $|F_{B}-F_{g}|$ . The low mass-density ratio of $m^{\ast }=0.02$ chosen here corresponds to a cylinder made from a very light material such as expanded polystyrene (Mathai et al. 2015). The buoyant pendulum with ( $m^{\ast }=0.02$ ) accelerates quickly, reaching a high angular speed in a short period of time. At the same time, owing to its low inertia, which in this case comes entirely from the fluid (or $m_{a}^{\ast }$ ), the motion damps out quickly. Therefore, the peak angular displacement at first minimum $\unicode[STIX]{x1D703}_{0}$ is lower for the buoyant pendulum as compared to the heavy pendulum ( $m^{\ast }=1.98$ ), while the peak angular velocity reached $\unicode[STIX]{x1D714}_{peak}$ is higher for the buoyant cylinder.

4.2 Model versus experiment

We now present a one-to-one comparison between the model and experiment. In figure 4(a), we show the normalised frequency of oscillation. As expected, the frequency of oscillation has a clear dependence on the mass ratio $m^{\ast }$ . When the pendulum mass density is close to the fluid density (or $m^{\ast }\rightarrow 1$ ), the frequencies are low. For heavy cylinders, the frequency increases with $m^{\ast }$ , and asymptotically approaches the frequency of a large amplitude simple pendulum $f_{\unicode[STIX]{x03C0}/2}^{\ast }$ . The predictions of the model (3.2) are given by the blue and red dashed curves. The blue curve corresponds to the prediction when the potential flow added mass $m_{a}^{\ast }=1.0$ is used. This slightly underpredicts the oscillation frequency. A best fit to the data is obtained for $m_{a}^{\ast }=0.53$ . Thus, for $C_{D}=1.2$ , $m_{a}^{\ast }=0.53$ and $\unicode[STIX]{x1D707}_{f}=0.2$ , the frequency predictions of the model are in good agreement with the experiments for both buoyant and heavy cylinders.

Figure 4. (a) Normalised oscillation frequency $f^{\ast }$ from experiment and model (3.2). Here, $f$ is averaged over the first four oscillations, and normalised by $f_{sp}=(1/2\unicode[STIX]{x03C0})\sqrt{g/L}$ . For the full equation of motion, using potential flow added mass $m_{a}^{\ast }=1.0$ leads to an underprediction of $f^{\ast }$ . For a reduced $m_{a}^{\ast }=0.53$ , the predictions of (3.2) match well with the experiments. Remarkably, an undamped model for a small amplitude simple pendulum, and with the same $m_{a}^{\ast }=0.53$ , reproduces the curve, providing evidence that the nonlinear drag plays only a weak role in the oscillation frequency $f^{\ast }$ . The green curve shows that an undamped model with large initial amplitude $\unicode[STIX]{x1D703}_{0}=\unicode[STIX]{x03C0}/2$ underpredicts the oscillation frequency. (b) Peak pendulum deceleration, $a_{max}$ versus $m^{\ast }$ .

It is interesting to compare the predictions of the model (3.2) when the drag and the bearing friction are ignored altogether. The solid green curve in figure 4(a) shows the frequency prediction for large amplitude oscillations when the potential flow added mass alone is accounted for. This underpredicts the oscillation frequency. However, when the same undamped model is used for small amplitude oscillations, the predictions improve significantly (see solid blue curve). This is because the nonlinear drag is highly effective in quickly reducing the oscillation amplitudes to modest values. A further improvement is obtained upon using a reduced added mass of $m_{a}^{\ast }=0.53$ . Neill et al. (2007) had noted that linear drag had only a small effect on the oscillation frequency. Here we found that the same holds for the nonlinear drag term.

The frequency of oscillation represents an averaged quantity for the pendulum motion considered here. A more sensitive quantity would be the peak deceleration  $a_{max}$ of the oscillating pendulum, which signifies the instant when the largest force acts on the pendulum. Following the initial instants of release of the cylinder, a peak in the deceleration is experienced when the cylinder reaches the end of the first swing, i.e. points such as (p1) and (p2) in figure 2(b). In figure 4(b), we compare the experimentally measured peak deceleration with the predictions of the model (3.2). The peak deceleration changes significantly with $m^{\ast }$ , and the model predictions are in reasonable agreement with the experimental measurements. Interestingly, the value of the added mass coefficient has only a minor role in $a_{max}$ . Instead, the fluid drag has a leading role in the peak deceleration of the pendulum. Thus, the frequency curve (figure 4 a) provides a probe that singles out the added mass effect, while the deceleration curve (figure 4 b) probes mainly the nonlinear drag’s effect.

Figure 3(a) also provides insight into an interesting aspect of the oscillation decay for heavy pendulums. For $m^{\ast }>1$ cases, the oscillations decay quicker with increasing $m^{\ast }$ . This occurred despite the higher inertia of heavier pendulums, and seems counter-intuitive based on our common knowledge of pendulums oscillating in air. However, this trend cannot continue to large $m^{\ast }$ , since in the limit of very large inertia ( $m^{\ast }\gg 1$ ) the decay rate should reduce with $m^{\ast }$ . This implies the existence of an optimal heavy pendulum ( $m_{opt}^{\ast }>1$ ) for which one observes the quickest damping. We will provide an approximate model to predict $m_{opt}^{\ast }$ , corresponding to the quickest damped pendulum.

Figure 5. (a) Time to decay 60 % of the initial amplitude $\tilde{\unicode[STIX]{x1D70F}}_{60\,\%}$ versus mass-density ratio $m^{\ast }$ . (b) Amplitude envelope $\unicode[STIX]{x1D703}_{m}$ versus $m^{\ast }$ after a time $\unicode[STIX]{x1D70F}_{ref}=3\unicode[STIX]{x03C0}\sqrt{L/g}$ . For the heavy pendulums, an optimal $m^{\ast }$ and optimal $f^{\ast }$ are visible in both experiment and simulations. (c) $\tilde{\unicode[STIX]{x1D70F}}_{60\,\%}$ versus normalised oscillation frequency $f^{\ast }=f/f_{sp}$ , where $f_{sp}=(1/2\unicode[STIX]{x03C0})\sqrt{g/L}$ . (d) $\unicode[STIX]{x1D703}_{m}$ versus $f^{\ast }$ at $\unicode[STIX]{x1D70F}_{ref}=3\unicode[STIX]{x03C0}\sqrt{L/g}$ . Note that (c,d) show heavy cases only. The basic model overpredicts the decay time and amplitude envelope for all cases.

4.3 Amplitude decay

The initial damping is determined mainly by the nonlinear drag; therefore, we ignore the bearing friction force in the equation of motion (3.2). This yields: $m_{eff}^{\ast }\text{d}^{2}\unicode[STIX]{x1D703}/\text{d}\tilde{t}^{2}=-k\sin \unicode[STIX]{x1D703}-c|\text{d}\unicode[STIX]{x1D703}/\text{d}\tilde{t}|\text{d}\unicode[STIX]{x1D703}/\text{d}\tilde{t}$ . In the absence of fluid drag ( $c=0$ ), the solution is given by $\unicode[STIX]{x1D703}(\tilde{t})=\unicode[STIX]{x1D703}_{0}\cos \unicode[STIX]{x1D714}_{0}\tilde{t}$ , with $\unicode[STIX]{x1D714}_{0}=\sqrt{k/m_{eff}^{\ast }}$ , and $\unicode[STIX]{x1D703}_{0}$ the initial angle. If the damping constant $c$ is small such that the pendulum amplitude is approximately constant over one period, we can assume that the solution is still approximately the same, but with a maximum amplitude $\unicode[STIX]{x1D703}_{m}(\tilde{t})$ which slowly decays in time. Consequently, the energy in the oscillations evolves as $E(\tilde{t})\approx (1/2)k\unicode[STIX]{x1D703}_{m}(\tilde{t})^{2}$ . The overall decay rate of the energy can be equated to the average work done by the drag force over one period $T=2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D714}_{0}$ :

(4.1) $$\begin{eqnarray}\frac{\text{d}E}{\text{d}\tilde{t}}=\frac{\unicode[STIX]{x1D714}_{0}}{2\unicode[STIX]{x03C0}}\int _{0}^{2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D714}_{0}}F_{D}(\tilde{t})\dot{\unicode[STIX]{x1D703}}(\tilde{t})\,\text{d}\tilde{t}.\end{eqnarray}$$

Using $E(\tilde{t})\approx (1/2)k\unicode[STIX]{x1D703}_{m}(\tilde{t})^{2}$ and $F_{D}=-c\dot{\unicode[STIX]{x1D703}}^{2}$ , we obtain a differential equation for the slowly varying envelope $\unicode[STIX]{x1D703}_{m}(t)$ of the form: $\text{d}\unicode[STIX]{x1D703}_{m}/\text{d}\tilde{t}\approx -{\mathcal{B}}(m^{\ast })\unicode[STIX]{x1D703}_{m}^{2}$ , where ${\mathcal{B}}(m^{\ast })=(4c/3\unicode[STIX]{x03C0})\sqrt{|m^{\ast }-1|/(m^{\ast }+m_{a}^{\ast })^{3}}$ . This yields the time to decay to a certain fraction ( $1-\unicode[STIX]{x1D6FC}$ ) of the initial amplitude $\unicode[STIX]{x1D703}_{0}$ as $\tilde{\unicode[STIX]{x1D70F}}_{(1-\unicode[STIX]{x1D6FC})}=(\unicode[STIX]{x1D6FC}/(1-\unicode[STIX]{x1D6FC}))/\unicode[STIX]{x1D703}_{0}{\mathcal{B}}(m^{\ast })$ .

For $m^{\ast }>1$ , and using the best fit $m_{a}^{\ast }=0.53$ , ${\mathcal{B}}$ has a maxima at $m_{opt}^{\ast }\approx 1.75$ . In other words, for a given $\unicode[STIX]{x1D703}_{0}$ the damping is the quickest when the pendulum is nearly twice the density of the fluid. Alternately, one could plot the amplitude envelope $\unicode[STIX]{x1D703}_{m}$ after a dimensionless time $\tilde{\unicode[STIX]{x1D70F}}_{ref}\equiv \unicode[STIX]{x1D70F}_{ref}/\sqrt{L/g}$ for different $m^{\ast }$ . This yields an expression $\unicode[STIX]{x1D703}_{m}=1/{\mathcal{B}}(m)\tilde{\unicode[STIX]{x1D70F}}_{ref}+(1/\unicode[STIX]{x1D703}_{0})$ . In figure 5(a,b) we show comparisons between the model and experiment for $\tilde{\unicode[STIX]{x1D70F}}_{60\,\%}$ and $\unicode[STIX]{x1D703}_{m}$ , respectively, versus mass-density ratio $m^{\ast }$ . The black symbols denote the model predictions by solving the full equation of motion (3.2). The black curve shows the prediction by the approximate model described above. Note that we first fit an envelope to the amplitude decay, using the maxima (peaks) of $\unicode[STIX]{x1D703}$ versus $t$ curve, i.e. a curve that grazes through the maxima of the amplitude curve. The optimal damping at $m_{opt}^{\ast }\approx 1.75$ is clearly visible in both the simulations and the reduced model. The blue data points are from experiment. A similar behaviour with a minima at an optimal $m^{\ast }$ is visible in the experiments. At the same time we note that both $\tilde{\unicode[STIX]{x1D70F}}_{60\,\%}$ and $\unicode[STIX]{x1D703}_{m}$ are lower in the experiment. Further, the optimal $m^{\ast }$ value in experiment is slightly lower as compared to the model predictions. For very heavy pendulums ( $m^{\ast }\gg 1$ ), the damping ${\mathcal{B}}\propto 1/m^{\ast }$ , which leads to the linear relation $\tilde{\unicode[STIX]{x1D70F}}_{60\,\%}\propto m^{\ast }$ . This is also the commonly encountered situation for pendulums oscillating in air.

The above discussed non-monotonic damping behaviour may be understood as the interplay between the higher speeds achieved by the heavier pendulums, and the nonlinear growth of the drag in proportion to the square of the speed. The phenomenon observed here may hold some analogies to the added damping for oscillating cylinders (Dong 1978), which is usually expressed in terms of the oscillation frequency. Therefore, in figure 5(c,d) we plot $\tilde{\unicode[STIX]{x1D70F}}_{60\,\%}$ and $\unicode[STIX]{x1D703}_{m}$ , respectively, as a function of the normalised oscillation frequency $f^{\ast }$ . The $\unicode[STIX]{x1D703}_{m}$ plot shows a clear minima at $f^{\ast }\approx 0.6$ . Similar to the $m^{\ast }$ plots (figure 5 a,b), the model overpredicts both $\tilde{\unicode[STIX]{x1D70F}}_{60\,\%}$ and $\unicode[STIX]{x1D703}_{m}$ , which occurs despite using the optimal added mass coefficient $m_{a}^{\ast }=0.53$ determined earlier. Clearly, the origin of these must lie in the inadequacy of the simplified drag model used ( $C_{D}=1.2$ ). Nevertheless, the constant $C_{D}$ case can be viewed as the most basic model, which captures many essential features of the damping. Using a higher value of $C_{D}$ seems the logical choice. In the following we will show that accounting for a few phenomenological effects can improve the predictions. We will introduce these modifications to the model and equation of motion.

4.4 Drag corrections

Firstly, the influence of the instantaneous $Re$ on $C_{D}$ can be taken into account. This leads to a faster decay at small amplitudes (or small $v_{p}$ ), since $C_{D}$ is higher at lower Reynolds numbers (Hoerner 1965). Secondly, the forces due to vortex shedding behind the cylinder could be modelled into the drag term. The drag coefficient in the presence of vortex shedding may be approximated as $C_{Dvs}=C_{D}(1+k\sin \unicode[STIX]{x1D714}_{vs}t)$ , where $k\sim 0.1$ (Govardhan & Williamson 2000), and $\unicode[STIX]{x1D714}_{vs}=2\unicode[STIX]{x03C0}Stv_{p}/D$ is the vortex shedding frequency. Here $St$ is the Strouhal number (Govardhan & Williamson 2005), $v_{p}$ the cylinder velocity and $D$ the cylinder diameter. This leads to a marginal reduction in the predicted peaks of the oscillations.

Figure 6(a,b) compares the experiment and model predictions for the angular position $\unicode[STIX]{x1D703}$ versus time $t$ for a heavy pendulum with $m^{\ast }=4.98$ . The oscillation frequency matches well between the experiment and the model. However, deviations are seen in the second peak of oscillation, marked as (r) in figure 6(b). The deviation persists for all of the subsequent oscillations. This suggests that the major factor causing the deviations between the model and experiment is still missing in the equation of motion.

Figure 6. (a) Angular position $\unicode[STIX]{x1D703}$ versus time $t$ for $m^{\ast }=4.98$ . The red curve shows the model predictions. (b) Shows zoom-in of the second swing. The model overpredicts the maximum amplitude (point marked as (r) in a,b). (c–e) Normalised horizontal velocity $U_{x}/v_{0}$ (left) and normalised vorticity $\unicode[STIX]{x1D714}/(v_{0}/D)$ (right) at the three instants marked in (a) as (p), (q) and (r), respectively. Here, $v_{0}$ is the measured maximum speed of the pendulum.

4.5 PIV and wake corrections

To understand the origin of the deviation in the second and subsequent peaks of the oscillations, we employed high-speed PIV to quantify the flow around the cylinder during its motion. Figure 6(c) shows the normalised velocity and vorticity fields at a time instant when the cylinder is moving towards the right during its first swing (point (p) in figure 6 a). The wake behind the cylinder is clearly visible, while the flow ahead of the cylinder has negligible vorticity i.e. nearly irrotational (figure 6 c (right)). The wake is unsteady and is expected to induce unsteady forces on the cylinder during its swing.

One can expect that the mean drag on the cylinder at this instant is fairly predicted by the nonlinear drag term $F_{D}$ in the equation of motion (3.2). Consequently, the first swing of the pendulum is captured well by our model. Next we focus at a later instant in time (point (q) in figure 6 a,b), when the cylinder has just completed its first swing and is returning to the left. Figure 6(d) shows the velocity and vorticity fields at this instant. In this case, however, the cylinder is moving towards a disturbed background flow. The horizontal velocity plot (see left-side plot of figure 6 d) indicates that the cylinder faces an incoming flow to the right. The effect of this disturbed flow is not taken into account in our model. Therefore, as the cylinder travels further through the disturbed flow, it experiences a greater resistance than what is predicted by the drag term in our model. Figure 6(e) shows the flow field at the instant when the cylinder completed its reverse swing (point (r) in figure 6 a,b). By this time, the original wake in front of the cylinder appears to have dissipated.

Figure 7. (a) Schematic showing the estimation of mean wake velocity $U_{f}$ from a sample PIV flow field. A sample video of the PIV flow field is provided as supplemental material available at https://doi.org/10.1017/jfm.2018.867. $U_{f}$ is computed inside a $2D\times 2D$ window located at a selected angular position $\unicode[STIX]{x1D703}_{sel}$ . (b) Decay of wake velocity $U_{f}$ in time at various angular positions $\unicode[STIX]{x1D703}_{sel}$ . $t_{pass}$ is the time when the cylinder passed $\unicode[STIX]{x1D703}_{sel}$ . The initial oscillation of the wake shows clear indication of vortex shedding. However, after this phase the wake decay is smooth. Inset shows the same plot on log–log scale. While at short times $U_{f}$ shows oscillations, the long-time decay shows a power law decay.

One of the approximations of our model is that the cylinder velocity $v_{p}$ nearly equals the relative velocity $v_{rel}$ between the cylinder and the flow. However, the PIV snapshots reveal a strong mean flow after the first swing of the pendulum. This disturbed flow in the wake could result in $v_{rel}$ deviating significantly from $v_{p}$ , which is the likely reason for the overprediction of the maxima of the second and subsequent oscillations. We look at the mean fluid velocity $U_{f}$ at different angular positions in the cylinder wake. Figure 7(a) shows details of the $2D\times 2D$ window at a selected angular position $\unicode[STIX]{x1D703}_{sel}$ used in the estimation of $U_{f}$ . The evolution of $U_{f}$ in time $t-t_{pass}$ is shown in figure 7(b). Here, $t_{pass}$ is the time when the cylinder passed $\unicode[STIX]{x1D703}_{sel}$ . The wake decay shows some interesting characteristics. Initially, $U_{f}$ oscillates, which is a clear indication of vortex shedding in the cylinder wake during the initial moments after the cylinder has passed $\unicode[STIX]{x1D703}_{sel}$ . We also see variations in the peak for different angular positions $\unicode[STIX]{x1D703}_{sel}$ . This part of the curve is highly unsteady and difficult to model. However, later in time $U_{f}$ decays in a gradual and monotonic way. The inset to the figure shows the same plot on log–log scale. The long-time decay appears to follow a power law decay.

Figure 8. (a) Decay of the normalised wake velocity $\tilde{U} _{f}$ versus normalised angular position $\unicode[STIX]{x1D703}^{\ast }$ . The wake velocity is normalised by the pendulum velocity at equilibrium position during the first swing $v_{p0}$ , and the angular position is normalised by the peak angular position at the end of the first swing $\unicode[STIX]{x1D703}_{max}$ . Both $v_{p0}$ and $\unicode[STIX]{x1D703}_{max}$ can be found by solving (3.2), and do not require any input from PIV. The normalisation leads to a reasonable collapse of the data despite the wide variation in Reynolds number from $Re\sim 4000$ (for $m^{\ast }=1.16$ ) to $Re\sim 34\,000$ (for $m^{\ast }=4.98$ ). Note that the arrow of time is pointing opposite to the arrow of $\unicode[STIX]{x1D703}$ , as indicated in (a). (b) Error (%) in the peak amplitude after including the history force due to wake flow. The coloured bands show the error ranges for the original and revised model. The dashed lines show the mean errors for the two models, which decreases from 14.8 % (for the original model) to 3.8 % (for the revised model).

Several studies have addressed the decay of wakes for rising and falling particles (Wu & Faeth 1993; Bagchi & Balachandar 2004). These revealed many aspects of the temporal and spatial decays of wakes behind spheres in quiescent fluids, and in turbulent flows. Alméras et al. (2017) disentangled the wake decay behind rising bubbles into near wake and far fields, with their characteristic decay rates. The wake decay observed in figure 7 also represents the combined effect of the local decay of the wake and the advection of the wake with a velocity $U_{f}$ , i.e. $\text{d}U_{f}/\text{d}t=(\unicode[STIX]{x2202}U_{f}/\unicode[STIX]{x2202}t)+(U_{f}/L)(\unicode[STIX]{x2202}U_{f}/\unicode[STIX]{x2202}\unicode[STIX]{x1D703})$ . These terms can be very difficult to disentangle for the wide range of $m^{\ast }$ (or $Re$ ) in our experiments. Alternatively, one can look at the wake velocity at locations just upstream of the cylinder during its return swing. Figure 8(a) shows the normalised wake velocity $\tilde{U} _{f}$ as a function of normalised angular position $\unicode[STIX]{x1D703}^{\ast }$ for four different mass ratios. Here $U_{f}$ is normalised by the pendulum speed at equilibrium position of each $m^{\ast }$ case, thus representing a characteristic speed for that $m^{\ast }$ . For the heavier pendulums, $v_{p0}$ is comparable to the gravitational velocity scale $v_{g}$ , but it deviates as $m^{\ast }$ is reduced. Similarly, the maximum angle reached by the pendulum at the end of its first swing $\unicode[STIX]{x1D703}_{max}$ (point (q) in figure 6 a,b) is used in the normalisation for $\unicode[STIX]{x1D703}$ . For clarity, we also show the arrow of increasing time, indicating that the return swing starts at $\unicode[STIX]{x1D703}^{\ast }\rightarrow 1$ . With these normalisations the data show a reasonable collapse. It is remarkable that the data collapse despite the large variation in $m^{\ast }$ : [1.15, 4.98], which corresponds to a Reynolds number variation of almost a decade ( $Re\in$  [4000, 34 000]). At the initial instant, i.e. when $\unicode[STIX]{x1D703}^{\ast }\rightarrow 1$ , the cylinder sees a large incoming flow velocity. As the pendulum swings back to $\unicode[STIX]{x1D703}^{\ast }\sim 0$ , the wake velocity has almost completely decayed. Beyond this we observe a surprising increase in the incoming flow velocity. This increase arises from the cylinder wake during the first swing, where the pendulum had its highest swing velocity. Since the cylinder velocity at this position was large, the wake has not decayed completely.

We obtain an approximate fit for the variation of $\tilde{U} _{f}$ versus $\unicode[STIX]{x1D703}^{\ast }$ (see figure 8 a). With the wake flow modelled using $\tilde{U} _{f}(\unicode[STIX]{x1D703}^{\ast })$ shown above, we can include the effect of the incoming flow through a modified drag term: $F_{Df}\approx (1/2)\unicode[STIX]{x1D70C}_{f}A_{p}C_{D}(v_{p}-U_{f})^{2}$ , where $U_{f}=\tilde{U} _{f}v_{p0}$ . This can be implemented as a modified torque that replaces the original drag term in (3.2). The modified equation of motion reads:

(4.2) $$\begin{eqnarray}m_{eff}^{\ast }\frac{\text{d}^{2}\unicode[STIX]{x1D703}}{\text{d}\tilde{t}^{2}}=-k\sin \unicode[STIX]{x1D703}-c\left|\frac{\text{d}\unicode[STIX]{x1D703}_{rel}}{\text{d}\tilde{t}}\right|\frac{\text{d}\unicode[STIX]{x1D703}_{rel}}{\text{d}\tilde{t}}-h|\text{cos}\,\unicode[STIX]{x1D703}|\text{sgn}\left(\frac{\text{d}\unicode[STIX]{x1D703}}{\text{d}\tilde{t}}\right),\end{eqnarray}$$

where the relative angular velocity $\text{d}\unicode[STIX]{x1D703}_{rel}/\text{d}\tilde{t}=(v_{p}-U_{f})/\sqrt{Lg}$ . Note that (4.2) requires no new input parameters, and works for the full range of Reynolds number in our experiments, i.e. $Re\in [4000,34\,000]$ . The improvements in the predictions of the amplitude decay are presented in figure 8(b), which shows the error in the predicted amplitude of the second peak of oscillation for all mass ratios ( $m^{\ast }>1$ ) studied. The predictions of the revised model are shown along with those of the original model. Accounting for the wake flow results in significant improvement in the second peak for all $m^{\ast }$ cases. The mean error reduces from ${\sim}14.8$ to 3.8 % after including the influence of the wake flow into the equation of motion. These results demonstrate the importance of fluid–structure coupling for pendulums oscillating at large amplitudes in viscous fluids.

5 Conclusions and discussion

In this work, we have studied the large amplitude oscillations of heavy and buoyant cylindrical pendulums in water. The oscillation frequency and peak deceleration of the pendulum depend on the mass-density ratio $m^{\ast }$ . Accounting for buoyancy and added mass alone, as done in previous investigations (Neill et al. 2007), does not provide a complete picture of the dynamics. We have developed a basic theoretical model, which accounts for the nonlinear drag force and the bearing friction in addition to the buoyancy and the added mass. With the inclusion of these terms, the model predictions are in reasonable agreement with the experimental observations for a wide range of mass ratios. The added mass coefficient in experiment is found to be $m_{a}^{\ast }=0.53$ , i.e. lower that the two-dimensional (2-D) potential flow value for a cylinder. A reduced $m_{a}^{\ast }$ as compared to the 2-D potential flow added mass can be expected, given that the flow is viscous and three-dimensional, and also due to the finite span of the pendulum. This result is consistent with existing literature on oscillating bodies in viscous fluids (Dong 1978; Konstantinidis 2013; Koo & Kim 2015).

The theoretical model presented here provides fair predictions of the oscillation frequency and peak decelerations for a wide range of $m^{\ast }$ . However, the model overpredicts the peak oscillation amplitudes for the second and subsequent oscillations. Introducing a Reynolds number dependent drag, along with vortex shedding forces, can provide marginal improvements to the model predictions. While modelling $C_{D}$ as a function of $Re$ improves the later oscillations, the vortex shedding term mainly influences the first oscillation.

Particle image velocimetry (PIV) measurements have revealed that the major factor causing the deviations between experiment and the simplified model is the disturbed flow surrounding the cylinder, which was not accounted for in the original model. We have used insights from PIV measurements to obtain a simple model for the incoming wake flow for a wide range of mass ratios. With this history force modelled, the revised model predictions are significantly improved; the mean error reducing from 14.8 to 3.8 %. The predictions are also improved for other mass ratios for which PIV experiments were not performed.

Even with the wake history force included, the current model is still quite basic. In reality, the dynamics is highly nonlinear, with changes in direction, curvilinear trajectories and wide variations in instantaneous $Re$ . Under such conditions, exact analytical expressions for the drag, added mass and history forces are not available. Fully resolved direct numerical simulations (immersed boundary (Mittal & Iaccarino 2005), or Physalis (Naso & Prosperetti 2010)) can provide better insights into the flow-induced forces. On a different note, the present experiments have shown the importance of fluid drag for a bluff body oscillating in a fluid. An interesting extension to the study would be to use a more streamlined body, for which the nonlinear drag would be less dominant. In the absence of major flow separations, we can expect the frequency and dynamics to conform slightly better with a potential flow based added mass.