Exploring the non-linear oscillation of a rigid sphere on an elastic half-space

Nonlinear oscillations can be found in a wide range of physical systems [1]. Extended research results have been published regarding the mathematical analysis of the complex behavior of nonlinear oscillators [1]. It is significant to mention that the scientific interest regarding the nonlinear behavior of physical systems remains undiminished up to date for research [2–13] or educational purposes [14–19]. However, the procedure to find the expressions that describe the behavior of nonlinear oscillators with respect to time is a challenging task. The basic goal when a nonlinear oscillation is explored is to find approximate solutions of the related differential equation in terms of elementary functions. Using this approach, the solution can be easily applied to practical applications.

In this paper an interesting nonlinear oscillator will be explored, the rigid sphere—elastic half-space system. An elastic half-space is an elastic material that extends infinitely in all directions including depth with the surface at the top to be considered as boundary [20]. The behavior of an elastic half-space is described by its Young's modulus and Poisson's ratio. Despite the fact that an elastic half-space is an ideal mathematical concept, it has been used in a wide range of applications from nanotechnology [21–24] and contact mechanics [25, 26] to seismics [27–30]. In fact, is a useful tool that can be used to represent the behavior of materials that have a linear elastic response [31]. Thus, it can be used as a mean to provide a useful insight in a significant number of physical phenomena in which elastic objects interact.

The problem to be studied in this paper is the nonlinear oscillation of a rigid sphere on an elastic half-space. The aforementioned motion is described by a second-order nonlinear ordinary differential equation that presents significant complexity [32]. In a previous publication, it was shown that for small displacements, the sphere's motion can be considered as harmonic [32]. However, in this paper, cases with significant displacement range will be explored. The basic goal of this paper is to find accurate solutions using elementary functions that describe the behavior of the aforementioned physical system. In addition, the relation of the motion's period with respect to the range of displacement will be explored. Towards this purpose, a simple method will be used that relies on the energy conservation principle. The approach presented in this paper is simple and can be used in post-graduate physics lectures to provide a methodology for the analysis of complicated non-linear phenomena.

Assume a rigid sphere with mass m that is placed at the surface of an elastic half-space. The system elastic half-space—sphere is placed at a constant gravitational field, $\overrightarrow {g}$. At an arbitrary position the applied net force on the sphere equals to:

$$\begin{equation}{\Sigma} \overrightarrow {F}= \overrightarrow {\mathrm{w}}+ \overrightarrow {F}.\end{equation} \tag{ 1 }$$

In equation (1), $\overrightarrow {F}$ is the applied force on the sphere by the elastic half-space and $\overrightarrow {w}$ is the sphere's weight. According to the Hertz contact theory [32, 33],

$$\begin{equation}F=\frac{4}{3}{E}^{{\ast}}{R}^{1/2}{y}^{3/2}.\end{equation} \tag{ 2 }$$

In equation (2), F is the applied force, E* is the reduced modulus of the elastic half-space, R is the sphere's radius and y is the deformation of the surface from its initial position. The reduced modulus and the Young's modulus are connected by the following equation [32, 33]:

$$\begin{equation}{E}^{{\ast}}=\frac{E}{1-{\mathrm{v}}^{2}},\end{equation} \tag{ 3 }$$

where E is the Young's modulus and v is the Poisson's ratio. At this point it must be noted that equation (2) is valid under the assumption that y ≪ R (i.e. the deformation of the surface is significantly smaller compared to the sphere's radius). A typical limit presented in the literature is y ⩽ 0.1R [21, 32, 33]. In addition, at the equilibrium position the deformation of the surface is y = y_eq and the net force on the sphere is zero,

$$\begin{equation}{\Sigma} \overrightarrow {F}= \overrightarrow {0}{\Rightarrow}w-F=0{\Rightarrow}mg=\frac{4}{3}{E}^{{\ast}}{R}^{1/2}{y}_{\mathrm{e}\mathrm{q}}^{3/2}.\end{equation} \tag{ 4 }$$

At an arbitrary position y = y_eq + Δy, the differential equation that describes the sphere's motion is the following,

$$\begin{equation*}m\frac{{\mathrm{d}}^{2} \overrightarrow {y}}{\mathrm{d}{t}^{2}}= \overrightarrow {w}+ \overrightarrow {F}{\Rightarrow}m\frac{{\mathrm{d}}^{2}y}{\mathrm{d}{t}^{2}}=w-\frac{4}{3}{E}^{{\ast}}{R}^{1/2}{y}^{3/2}{\Rightarrow}\end{equation*}$$

$$\begin{equation}m\frac{{\mathrm{d}}^{2}y}{\mathrm{d}{t}^{2}}=\frac{4}{3}{E}^{{\ast}}{R}^{1/2}{y}_{\mathrm{e}\mathrm{q}}^{3/2}-\frac{4}{3}{E}^{{\ast}}{R}^{1/2}{y}^{3/2}.\end{equation} \tag{ 5 }$$

Equation (5) is a second-order nonlinear ordinary differential equation and as a result the sphere's motion is not harmonic. Free body diagrams for the equilibrium and an arbitrary position described respectively by equations (4) and (5) are also presented in figure 1.

Zoom In Zoom Out Reset image size

To derive the y = f(t) function, the energy conservation principle will be applied. The range of the sphere's displacements is assumed to be 0 ⩽ y ⩽ y₀, where y₀ = 0.1R (i.e. according to the restriction that was presented in the introduction theory). The force applied on the sphere by the elastic half-space is conservative since $\mathrm{r}\mathrm{o}\mathrm{t}\left( \overrightarrow {F}\right)= \overrightarrow {0}$. As a result, F = −dU_el/dy. In addition, using the work-energy theorem,

$$\begin{equation*}\mathrm{d}{U}_{\mathrm{e}\mathrm{l}}=-\mathrm{d}{W}_{\mathrm{e}\mathrm{l}}{\Rightarrow}\mathrm{d}{U}_{\mathrm{e}\mathrm{l}}=- \overrightarrow {F}\cdot \mathrm{d} \overrightarrow {y}{\Rightarrow}{U}_{\mathrm{e}\mathrm{l}}\left(y\right)=-{\int }_{0}^{y} \overrightarrow {F}\cdot \mathrm{d} \overrightarrow {y}.\end{equation*}$$

Hence, the potential energy of the elastic deformation (using also equation (2)) of the half-space can be calculated as below,

$$\begin{equation}{U}_{\mathrm{e}\mathrm{l}}\left(y\right)=-{\int }_{0}^{y}-\frac{4}{3}{E}^{{\ast}}{R}^{\frac{1}{2}}{y}^{\frac{3}{2}}\mathrm{d}y=\frac{8}{15}{E}^{{\ast}}{R}^{\frac{1}{2}}{y}^{\frac{5}{2}}.\end{equation} \tag{ 6 }$$

Thus, the energy conservation principle between the two positions as presented in figure 2 can be written as follows (assuming that at the initial position the sphere's velocity is zero, y-coordinate origin is at the undeformed surface and the y-axis is pointing down):

$$\begin{align}\hfill {E}_{k\left(1\right)}+{U}_{\mathrm{e}\mathrm{l}\left(1\right)}+{U}_{\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{v}\left(1\right)}& ={E}_{k\left(2\right)}+{U}_{\mathrm{e}\mathrm{l}\left(2\right)}+{U}_{\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{v}\left(2\right)}{\Rightarrow}\hfill \\ \hfill mg{y}_{0}& =\frac{1}{2}m{\left(\frac{\mathrm{d}y}{\mathrm{d}t}\right)}^{2}+\frac{8}{15}{E}^{{\ast}}{R}^{1/2}{y}^{5/2}+mg\left({y}_{0}-y\right){\Rightarrow}\hfill \\ \hfill 0& =\frac{1}{2}m{\left(\frac{\mathrm{d}y}{\mathrm{d}t}\right)}^{2}+\frac{8}{15}{E}^{{\ast}}{R}^{1/2}{y}^{5/2}-mgy.\hfill \end{align} \tag{ 7 }$$

Zoom In Zoom Out Reset image size

Using equation (7) it is derived that,

if y = 0, then $\frac{\mathrm{d}y}{\mathrm{d}t}=0$ (as it was already mentioned).

In addition, $\frac{\mathrm{d}y}{\mathrm{d}t}=0$ for y = y₀. Thus, (using also equation (4)), it is concluded:

$$\begin{equation}\frac{8}{15}{E}^{{\ast}}{R}^{1/2}{y}_{0}^{5/2}=mg{y}_{0}{\Rightarrow}{y}_{0}={\left(\frac{15mg}{8{E}^{{\ast}}{R}^{1/2}}\right)}^{2/3}{\Rightarrow}{y}_{\mathrm{e}\mathrm{q}}={\left(\frac{5}{2}\right)}^{-2/3}{y}_{0}.\end{equation} \tag{ 8 }$$

As a result, the position at which the net force is zero (equilibrium position) is not at y = y₀/2. In addition, using equation (4), equation (7) can be written in the form,

$$\begin{equation}0=\frac{1}{2}m{\left(\frac{\mathrm{d}y}{\mathrm{d}t}\right)}^{2}+\frac{8}{15}{E}^{{\ast}}{R}^{1/2}{y}^{5/2}-\frac{4}{3}{E}^{{\ast}}{R}^{1/2}{y}_{\mathrm{e}\mathrm{q}}^{3/2}y.\end{equation} \tag{ 9 }$$

By substituting equation (8) to (9), it can be derived,

$$\begin{equation}0=\frac{1}{2}m{\left(\frac{\mathrm{d}y}{\mathrm{d}t}\right)}^{2}+\frac{8}{15}{E}^{{\ast}}{R}^{1/2}{y}^{5/2}-\frac{8}{15}{E}^{{\ast}}{R}^{1/2}{y}_{0}^{3/2}y.\end{equation} \tag{ 10 }$$

To facilitate the analysis the dimensionless parameter x will be used, where x = y/R. Thus,

$$\begin{equation*}0{\leqslant}y{\leqslant}{y}_{0}{\Rightarrow}0{\leqslant}xR{\leqslant}0.1R{\Rightarrow}0{\leqslant}x{\leqslant}0.1.\end{equation*}$$

Using, the transformation y = xR, equation (10) can be written in the form,

$$\begin{align}\hfill {\left(\frac{\mathrm{d}x}{\mathrm{d}t}\right)}^{2}& =\frac{16}{15}\frac{{E}^{{\ast}}R}{m}{x}_{0}^{3/2}x-\frac{16}{15}\frac{{E}^{{\ast}}R}{m}{x}^{5/2}{\Rightarrow}\frac{\mathrm{d}x}{\mathrm{d}t}=\sqrt{\frac{{E}^{{\ast}}R}{m}}\sqrt{\left(\frac{16}{15}{x}_{0}^{3/2}x-\left(\frac{16}{15}\right){x}^{\frac{5}{2}}\right)}{\Rightarrow}\hfill \\ \hfill {\Delta}t& =\sqrt{\frac{m}{{E}^{{\ast}}R}}{\int }_{0}^{x}\sqrt{\frac{1}{\frac{16}{15}\left({x}_{0}^{\frac{3}{2}}x-{x}^{\frac{5}{2}}\right)}}\mathrm{d}x,\hfill \end{align}$$

where

$$\begin{equation}0{\leqslant}x{\leqslant}0.1.\end{equation} \tag{ 11 }$$

Or, assuming that, t₀ = 0, Δt = t − t₀ = t,

$$\begin{equation}t=\sqrt{\frac{m}{{E}^{{\ast}}R}}{t}^{\prime }\end{equation} \tag{ 12 }$$

where,

$$\begin{equation}{t}^{\prime }={\int }_{0}^{x}\sqrt{\frac{1}{\frac{16}{15}\left({x}_{0}^{\frac{3}{2}}x-{x}^{\frac{5}{2}}\right)}}\mathrm{d}x,\quad 0{\leqslant}x{\leqslant}0.1.\end{equation} \tag{ 13 }$$

The magnitude t' calculated by equation (13) is a dimensionless parameter that if it is multiplied with the factor $\sqrt{\frac{m}{{E}^{{\ast}}R}}$, it will provide the moment at which the sphere is at the arbitrary position y = xR. The solution of the integral provided by equation (13) is presented as follows,

$$\begin{equation*}{t}^{\prime }\left(x\right)=\frac{1}{2}\sqrt{15}\sqrt{\frac{x}{{x}_{0}^{\frac{3}{2}}}}{}_{2}F_{1}\left(\frac{1}{3},\frac{1}{2};\frac{4}{3};\frac{{x}^{\frac{3}{2}}}{{x}_{0}^{\frac{3}{2}}}\right),\end{equation*}$$

where, ₂ ${F}_{1}\left(a,b;c;z\right)$ is the p = 2 and q = 1 generalised hypergeometric function. It must be also noted that the abovementioned function is only valid in the regime ${\left(\frac{{x}_{0}}{x}\right)}^{3/2}{ >}1$ which is almost the complete domain 0 ⩽ x ⩽ x₀ = 0.1. Although the aforementioned equation cannot be inverted, by tabulating $\left({t}^{\prime }\left(x\right),x\right)$, with x₀ ∈ (0, 0.1) and plotting such list, the graph of the inverted function, i.e. of x = f(t') plot is presented in figure 3(a). Many functions were tested to fit the data. It was found that a Fourier series function, perfectly fits the data. More specifically, a function in the following form was used for data fitting,

$$\begin{equation*}x\left({t}^{\prime }\right)={a}_{0}+\sum\limits _{n=1}^{k}\enspace {a}_{n}\enspace \mathrm{cos}\left(n{\omega }^{\prime }{t}^{\prime }\right)+\sum\limits _{n=1}^{k}\enspace {b}_{n}\enspace \mathrm{sin}\left(n{\omega }^{\prime }{t}^{\prime }\right).\end{equation*}$$

Zoom In Zoom Out Reset image size

Using k = 2, the fit was perfect since the R-squared coefficient resulted equal to 1 ($\left.{R}_{\mathrm{s}.\mathrm{c}.}^{2}=1.0000\right)$. Thus, the x = f(t') is presented as follows:

$$\begin{equation}x\left({t}^{\prime }\right)={a}_{0}+{a}_{1}\enspace \mathrm{cos}\left({\omega }^{\prime }{t}^{\prime }\right)+{b}_{1}\enspace \mathrm{sin}\left({\omega }^{\prime }{t}^{\prime }\right)+{a}_{2}\enspace \mathrm{cos}\left(2{\omega }^{\prime }{t}^{\prime }\right)+{b}_{2}\enspace \mathrm{sin}\left(2{\omega }^{\prime }{t}^{\prime }\right).\end{equation} \tag{ 14 }$$

In equation (14), a₀ = 0.047 60, a₁ = −0.049 99, b₁ = −8.776 × 10⁻⁵, a₂ = 0.002 408, b₂ = −2.918 × 10⁻⁵, ω' = 0.6506.

The function that represents the position (y) with respect to time (t), y = f(t) can be then found (figure 3(b)),

$$\begin{equation}y\left(t\right)=R\left({a}_{0}+{a}_{1}\enspace \mathrm{cos}\left(\omega t\right)+{b}_{1}\enspace \mathrm{sin}\left(\omega t\right)+{a}_{2}\enspace \mathrm{cos}\left(2\omega t\right)+{b}_{2}\enspace \mathrm{sin}\left(2\omega t\right)\right)\end{equation} \tag{ 15 }$$

where, 0 ⩽ t < ∞.

In equation (15), t is given by equation (12), and $\omega ={\omega }^{\prime }\sqrt{\frac{{E}^{{\ast}}R}{m}}$, (i.e. ω't' = ωt).

In addition, according to figure 3, at x = 0, t' = 0, while at x = 0.1, t' = 4.8288. The half-period of the motion (duration of the motion from the position y = 0 to the position y = y₀) is equal to T/2, where $\frac{T}{2}=\sqrt{\frac{m}{{E}^{{\ast}}R}}{t}^{\prime }=4.8288\sqrt{\frac{m}{{E}^{{\ast}}R}}$, thus $T=9.6575\sqrt{\frac{m}{{E}^{{\ast}}R}}={T}^{\prime }\sqrt{\frac{m}{{E}^{{\ast}}R}}$ (where, T' = 9.6575).

However, in order to test if equation (15), accurately represents the position of the sphere with respect to time, it should be investigated whether the solution satisfies the second-order nonlinear ordinary differential equation, presented at section 2, (equation (5)). For this purpose, equation (5) was written as follows,

$$\begin{align}\hfill \frac{{\mathrm{d}}^{2}y}{\mathrm{d}{t}^{2}}& =\frac{4{E}^{{\ast}}{R}^{1/2}{y}_{\mathrm{e}\mathrm{q}}^{3/2}}{3m}-\frac{4{E}^{{\ast}}{R}^{1/2}}{3m}{y}^{3/2}{\Rightarrow}\frac{{E}^{{\ast}}{R}^{2}}{m}\frac{{\mathrm{d}}^{2}x}{\mathrm{d}{t}^{\prime 2}}\hfill \\ \hfill & =\frac{4{E}^{{\ast}}{R}^{2}{x}_{\mathrm{e}\mathrm{q}}^{3/2}}{3m}-\frac{4{E}^{{\ast}}{R}^{2}}{3m}{x}^{3/2}{\Rightarrow}\frac{{\mathrm{d}}^{2}x}{\mathrm{d}{t}^{{\prime }^{2}}}=\frac{8}{15}{x}_{0}^{3/2}-\frac{4}{3}{x}^{\frac{3}{2}}.\hfill \end{align}$$

Thus, the functions,

$$\begin{equation*}z\left({t}^{\prime }\right)=\frac{8}{15}{x}_{0}^{3/2}-\frac{4}{3}{\left[x\left({t}^{\prime }\right)\right]}^{3/2}\quad \mathrm{a}\mathrm{n}\mathrm{d}\quad q\left({t}^{\prime }\right)=\frac{{\mathrm{d}}^{2}x\left({t}^{\prime }\right)}{\mathrm{d}{t}^{\prime 2}}\end{equation*}$$

were plotted and compared (figure 4(a)). In particular, $q\left({t}^{\prime }\right)$ function, was fitted to $z\left({t}^{\prime }\right)$ data in the domain 0 ⩽ t' ⩽ 5T'. The R-squared coefficient resulted 1, ($\left.{R}_{\mathrm{s}.\mathrm{c}.}^{2}=1.0000\right)$. The same accuracy was presented for bigger domains as well (e.g. 0 ⩽ t' ⩽ 100T'). As a result, equation (15) is an accurate solution of equation (5). The position of the sphere with respect to time over the interval 0 ⩽ t ⩽ 5T is also presented in figure 4(b).

Zoom In Zoom Out Reset image size

In section 3 it was assumed that the initial position of the sphere (i.e. the position at t = 0) was y = 0. In this section, several other cases will be examined. If the elastic half-space is deformed when the sphere is at its initial position (i.e. y = y_in at t = 0), then ${U}_{\mathrm{e}\mathrm{l}\left(1\right)}=\frac{8}{15}{E}^{{\ast}}{R}^{1/2}{y}_{\mathrm{i}\mathrm{n}}^{5/2}$ and U_grav(1) = mgy₀. In addition, at an arbitrary position 2, ${U}_{\mathrm{e}\mathrm{l}\left(2\right)}=\frac{8}{15}{E}^{{\ast}}{R}^{1/2}{\left({y}^{\prime }+{y}_{\mathrm{i}\mathrm{n}}\right)}^{5/2}$ and ${U}_{\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{v}\left(2\right)}=mg\left({y}_{0}-{y}^{\prime }\right)$, where y' is the sphere's distance from its initial position y_in and y₀ the oscillation's range of displacements. In addition, it is assumed that E_k(1) = 0 and ${E}_{k\left(2\right)}=\frac{1}{2}m{\left(\frac{\mathrm{d}{y}^{\prime }}{\mathrm{d}t}\right)}^{2}$ (figure 5). Thus, the conservation of energy principle can be written as follows:

$$\begin{equation*}{E}_{k\left(1\right)}+{U}_{\mathrm{e}\mathrm{l}\left(1\right)}+{U}_{\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{v}\left(1\right)}={E}_{k\left(2\right)}+{U}_{\mathrm{e}\mathrm{l}\left(2\right)}+{U}_{\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{v}\left(2\right)}{\Rightarrow}\end{equation*}$$

$$\begin{align*}\hfill \frac{8}{15}{E}^{{\ast}}{R}^{1/2}{y}_{\mathrm{i}\mathrm{n}}^{5/2}+mg{y}_{0}& =\frac{1}{2}m{\left(\frac{\mathrm{d}{y}^{\prime }}{\mathrm{d}t}\right)}^{2}+\frac{8}{15}{E}^{{\ast}}{R}^{1/2}{\left({y}^{\prime }+{y}_{\mathrm{i}\mathrm{n}}\right)}^{5/2}+mg\left({y}_{0}-{y}^{\prime }\right)\hfill \end{align*}$$

$$\begin{equation}\frac{1}{2}m{\left(\frac{\mathrm{d}{y}^{\prime }}{\mathrm{d}t}\right)}^{2}=\frac{8}{15}{E}^{{\ast}}{R}^{1/2}{y}_{\mathrm{i}\mathrm{n}}^{5/2}-\frac{8}{15}{E}^{{\ast}}{R}^{\frac{1}{2}}{\left({y}^{\prime }+{y}_{\mathrm{i}\mathrm{n}}\right)}^{\frac{5}{2}}+mg{y}^{\prime }.\end{equation} \tag{ 16 }$$

Thus, $t=\sqrt{\frac{m}{{E}^{{\ast}}R}}{t}^{\prime }$, where

$$\begin{equation}{t}^{\prime }={\int }_{0}^{{x}^{\prime }}\sqrt{\frac{1}{\frac{16}{15}{x}_{0}^{\frac{3}{2}}{x}^{\prime }+\left(\frac{16}{15}\right)\left[{x}_{\mathrm{i}\mathrm{n}}^{\frac{5}{2}}-{\left({x}_{\mathrm{i}\mathrm{n}}+{x}^{\prime }\right)}^{\frac{5}{2}}\right]}}\mathrm{d}{x}^{\prime }\qquad \mathrm{a}\mathrm{n}\mathrm{d}\enspace 0{\leqslant}{x}^{\prime }{\leqslant}{x}_{0},\enspace {x}^{\prime }={y}^{\prime }/R.\end{equation} \tag{ 17 }$$

Zoom In Zoom Out Reset image size

A numerical methodology was followed to derive the x' = f(t') function in each case. Briefly, multiple calculations of the integral as provided from equation (17) were performed as follows:

$$\begin{equation*}0{\leqslant}{x}^{\prime }{\leqslant}{x}_{1}^{\prime }\to \mathrm{c}\mathrm{a}\mathrm{l}\mathrm{c}\mathrm{u}\mathrm{l}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{n}\mathrm{g}\enspace {t}_{1}^{\prime }\end{equation*}$$

$$\begin{equation*}0{\leqslant}{x}^{\prime }{\leqslant}{x}_{1}^{\prime }+{\Delta}{x}^{\prime }\to \mathrm{c}\mathrm{a}\mathrm{l}\mathrm{c}\mathrm{u}\mathrm{l}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{n}\mathrm{g}\enspace {t}_{1}^{\prime }+{\Delta}{t}^{\prime }\end{equation*}$$

$$\begin{equation*}0{\leqslant}{x}^{\prime }{\leqslant}{x}_{1}^{\prime }+2{\Delta}{x}^{\prime }\to \mathrm{c}\mathrm{a}\mathrm{l}\mathrm{c}\mathrm{u}\mathrm{l}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{n}\mathrm{g}\enspace {t}_{1}^{\prime }+2{\Delta}{t}^{\prime }\end{equation*}$$

$$\begin{equation*}\cdot \end{equation*}$$

$$\begin{equation*}\cdot \end{equation*}$$

$$\begin{equation*}\cdot \end{equation*}$$

$$\begin{equation*}0{\leqslant}{x}^{\prime }{\leqslant}{x}_{0}\to \mathrm{c}\mathrm{a}\mathrm{l}\mathrm{c}\mathrm{u}\mathrm{l}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{n}\mathrm{g}\enspace {t}_{0}^{\prime }\end{equation*}$$

It was selected, ${x}_{1}^{\prime }={\Delta}{x}^{\prime }=0.001$. Using the aforementioned integral (equation (17)) the x' = f(t') and the y = f(t) functions were determined (where y = y' + y_in). In particular, the y = f(t) function was determined for 5 additional cases (for initial position equal to 0.01R,  0.02R,  0.03R,  0.04R and 0.05R, without initial velocity). In every case, the solution was in the form of equation (15). The fitting factors for each case are presented in table 1. In addition, the y = f(t) graphs for each case (included the case studied in section 3) are presented in figure 6. In every case ${y}_{\mathrm{e}\mathrm{q}}=0.1{\left(\frac{5}{2}\right)}^{-2/3}R$ (in the cases of section 4 there is also a correlation of m, E*, R in order to achieve that: ${y}_{\mathrm{e}\mathrm{q}}={\left(\frac{5}{2}\right)}^{-2/3}0.1R$). In other words, we use 'the exact same sphere' on the 'exact same half-space'. The only difference is that the kinetic energy of the sphere at the equilibrium position in each case is smaller compared to the case of section 3. Thus, the range of displacements will be also smaller.

Zoom In Zoom Out Reset image size

Using figure 6(a), it can be concluded that the period of the motion depends slightly on the range of displacements. In particular, there is a slight decrease in the period's motion as the range of displacements decreases. Additionally, in figure 6(b) the phase diagram for each case presented in figure 6(a) is shown.

Several calculations of the motion's period with respect to the range of displacements were performed and presented in figure 7. It was assumed for every case that ${y}_{\mathrm{e}\mathrm{q}}={\left(\frac{5}{2}\right)}^{-2/3}{y}_{0}$. An exponential fit showed that the function T = f(y₀) follows the exponential law:

$$\begin{equation}T=\sqrt{\frac{m}{{E}^{{\ast}}R}}\left(9.132\enspace 00{e}^{-0.1149{y}_{0}}+0.070\enspace 56{e}^{21.8600{y}_{0}}\right).\end{equation} \tag{ 18 }$$

Zoom In Zoom Out Reset image size

For very small amplitudes around the equilibrium position y = y_eq (i.e. for very small range of displacements y₀) it is concluded that $T\approx 9.20\sqrt{\frac{m}{{E}^{{\ast}}R}}$. However, this is the case of harmonic motion. The aforementioned statement can be easily proved using a Taylor's series expansion (keeping only the linear term):

$$\begin{equation*}F=\frac{4}{3}{E}^{{\ast}}{R}^{1/2}{y}^{3/2}\cong F\left({y}_{\mathrm{e}\mathrm{q}}\right)+\frac{\mathrm{d}F\left({y}_{\mathrm{e}\mathrm{q}}\right)}{\mathrm{d}y}\left(y-{y}_{\mathrm{e}\mathrm{q}}\right){\Rightarrow}\end{equation*}$$

$$\begin{equation}F=2{E}^{{\ast}}{R}^{\frac{1}{2}}{y}_{\mathrm{e}\mathrm{q}}^{\frac{1}{2}}y-\frac{2}{3}{E}^{{\ast}}{R}^{\frac{1}{2}}{y}_{\mathrm{e}\mathrm{q}}^{\frac{3}{2}}.\end{equation} \tag{ 19 }$$

At an arbitrary position y, the differential equation that describes the sphere's motion is the following,

$$\begin{equation*}m\frac{{\mathrm{d}}^{2} \overrightarrow {y}}{\mathrm{d}{t}^{2}}= \overrightarrow {w}+ \overrightarrow {F}{\Rightarrow}m\frac{{\mathrm{d}}^{2}y}{\mathrm{d}{t}^{2}}=w-2{E}^{{\ast}}{R}^{\frac{1}{2}}{y}_{\mathrm{e}\mathrm{q}}^{\frac{1}{2}}y+\frac{2}{3}{E}^{{\ast}}{R}^{\frac{1}{2}}{y}_{\mathrm{e}\mathrm{q}}^{\frac{3}{2}}.\end{equation*}$$

Using also equation (4) we conclude in:

$$\begin{equation}\frac{{\mathrm{d}}^{2}y}{\mathrm{d}{t}^{2}}+\frac{2{E}^{{\ast}}{R}^{\frac{1}{2}}{y}_{\mathrm{e}\mathrm{q}}^{\frac{1}{2}}}{m}\left(y-{y}_{\mathrm{e}\mathrm{q}}\right)=0.\end{equation} \tag{ 20 }$$

In addition, ${y}_{\mathrm{e}\mathrm{q}}={\left(\frac{5}{2}\right)}^{-2/3}{y}_{0}={\left(\frac{5}{2}\right)}^{-\frac{2}{3}}0.1R$. Thus, the angular frequency in this case is $\omega =\sqrt{\frac{2{E}^{{\ast}}{R}^{\frac{1}{2}}{y}_{\mathrm{e}\mathrm{q}}^{\frac{1}{2}}}{m}}={y}_{\mathrm{e}\mathrm{q}}^{1/4}\sqrt{\frac{2{E}^{{\ast}}{R}^{\frac{1}{2}}}{m}}=\sqrt{2}{\left[0.1{\left(\frac{5}{2}\right)}^{-\frac{2}{3}}\right]}^{1/4}\sqrt{\frac{{E}^{{\ast}}R}{m}}=0.68\sqrt{\frac{{E}^{{\ast}}R}{m}}$, so the period of the motion is $T=\frac{2\pi }{\omega }\approx 9.20\sqrt{\frac{m}{{E}^{{\ast}}R}}$.

In addition, it has been also previously shown that for small displacements the sphere's motion can be considered as harmonic with period [32],

$$\begin{equation}T=2\pi \sqrt{\frac{m}{2r{E}^{{\ast}}}},\end{equation} \tag{ 21 }$$

where r is the radius of the circle at contact depth (which is considered to be constant for very small displacements). Hence, assuming a very small displacement from the equilibrium position ${y}_{\mathrm{e}\mathrm{q}}=0.1{\left(\frac{5}{2}\right)}^{-2/3}R$, the radius of the circle at contact depth is approximately,

$$\begin{equation}r=\sqrt{R{y}_{\mathrm{e}\mathrm{q}}}=0.233R.\end{equation} \tag{ 22 }$$

By substituting, equation (22) to (21), it can be easily concluded that $T\approx 9.20\sqrt{\frac{m}{{E}^{{\ast}}R}}$.

Equation (21) can be approximately used for any case, if the average radius at contact depth $\left(\bar{r}\right)$ is used,

$$\begin{equation}T=2\pi \sqrt{\frac{m}{2\bar{r}{E}^{{\ast}}}}.\end{equation} \tag{ 23 }$$

The average radius at contact depth is determined as follows:

$$\begin{equation}\bar{r}=\sqrt{R\bar{y}}\quad \mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}\quad \bar{y}=\frac{1}{T/2}{\int }_{0}^{T/2}y\left(t\right)\mathrm{d}t,\end{equation} \tag{ 24 }$$

Equations (23) and (24) provide an accurate approach in the domain 0 ⩽ y₀ ⩽ 0.025R. If y₀ = 0.1R the aforementioned approach leads to 1.49% lower value regarding the motion's period. However, equations (23) and (24) provide the physical significance of the exponential law presented by equation (18). In particular, the term $2\bar{r}{E}^{{\ast}}$ can be considered as the 'average stiffness' of the system. For small oscillations over the position ${y}_{\mathrm{e}\mathrm{q}}={\left(\frac{5}{2}\right)}^{-\frac{2}{3}}0.1R$, the radius at contact depth is $\bar{r}\approx \sqrt{R{y}_{\mathrm{e}\mathrm{q}}}={\left(0.1\right)}^{1/2}{\left(\frac{5}{2}\right)}^{-1/3}R$. If the oscillation's range of displacements increases it can be easily shown that $\bar{y}{< }{y}_{\mathrm{e}\mathrm{q}}$, and as a result the period is greater compared to the previous value. This fact explains the period's increase with the increase of the range of displacements.

In addition, an interesting question also arises; how small the range of displacements should be in order to consider the motion as harmonic? According to equation (19),

$$\begin{align}\hfill \frac{4}{3}{E}^{{\ast}}{R}^{1/2}{y}^{3/2}& \cong 2{E}^{{\ast}}{R}^{\frac{1}{2}}{y}_{\mathrm{e}\mathrm{q}}^{\frac{1}{2}}y-\frac{2}{3}{E}^{{\ast}}{R}^{\frac{1}{2}}{y}_{\mathrm{e}\mathrm{q}}^{\frac{3}{2}}{\Rightarrow}{y}^{3/2}\hfill \\ \hfill & \cong \frac{3}{2}{y}_{\mathrm{e}\mathrm{q}}^{1/2}y-\frac{1}{2}{y}_{\mathrm{e}\mathrm{q}}^{\frac{3}{2}}{\Rightarrow}{x}^{3/2}\hfill \\ \hfill & \cong \frac{3}{2}{x}_{\mathrm{e}\mathrm{q}}^{\frac{1}{2}}x-\frac{1}{2}{x}_{\mathrm{e}\mathrm{q}}^{\frac{3}{2}}.\hfill \end{align} \tag{ 25 }$$

In equation (25), ${x}_{\mathrm{e}\mathrm{q}}=0.1{\left(\frac{5}{2}\right)}^{-2/3}$ and 0 ⩽ x ⩽ 0.1. Equation (25) is valid for a small range of displacements around x_eq. To find 'how small' the range of displacements should be a graphical representation of the functions $f\left(x\right)={x}^{3/2}$ and $g\left(x\right)=\frac{3}{2}{x}_{\mathrm{e}\mathrm{q}}^{1/2}x-\frac{1}{2}{x}_{\mathrm{e}\mathrm{q}}^{\frac{3}{2}}$ is presented in figure 8(a). In addition, the function $h\left(x\right)=f\left(x\right)-g\left(x\right)$ is also presented in figure 8(b). The harmonic motion approximation is valid for $h\left(x\right)\to 0$. According to figure 8(c), a rational approximation should be:

$$\begin{equation*}0.0539{\leqslant}x{\leqslant}0.0547{\Rightarrow}0.0539R{\leqslant}y{\leqslant}0.0547R.\end{equation*}$$

Zoom In Zoom Out Reset image size

In other words, the sphere's initial position should be y_in = 0.0539R. In this case the sphere's motion can be accurately considered as harmonic.

As an alternative approximate solution of equation (5) (in terms of the range of sphere's displacements) the following function can be used,

$$\begin{equation}y={y}_{0}{\left\vert \mathrm{sin}\left(\frac{\omega t}{2}\right)\right\vert }^{n}+{y}_{\mathrm{i}\mathrm{n}}.\end{equation} \tag{ 26 }$$

In particular, equation (26) was fitted to the y versus t data for the case of y₀ = 0.1R and y_in = 0 and n resulted in 2.271 (in addition, ω' = 0.6506, since also in this case ωt = ω't'). Furthermore, the R-squared coefficient resulted in ${R}_{\mathrm{s}.\mathrm{c}.}^{2}=0.9997$. Thus, equation (26) can be considered as an acceptable approximate solution. Equation (26) does not provide the accuracy of equation (15) (see figure 9), however it has the significant advantage that it is expressed in terms of the range of displacements (y₀) in each case. In addition, it reveals that y ⩾ 0 in every case. Thus, the only unknown parameter that must be determined is the exponential factor n. Equation (26) was also tested for cases with smaller range of displacement and the results are presented in figure 10. In every case, ωt = ω't' (see also table 1). As it can be clearly seen, equation (26) can be better applied for a range of displacements within the domain 0.28R ⩽ y ⩽ 0.65R. However, it also provides acceptable solutions for the other cases. An interesting fact to be noted is that the range of values of n parameter was in a relatively small range 2.090 ⩽ n ⩽ 2.271 for all the examined cases.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this paper, the nonlinear oscillation of a rigid sphere on an elastic half-space was explored. It was proven that a Fourier series expansion as provided by equation (15) is an accurate solution of differential equation (5) and can be used to provide the sphere's position with respect to time for different range of displacements around the equilibrium position. In addition, an approximate solution (equation (26)) that is expressed in terms of the range of displacements can be also used. A significant result that was proven by the presented analysis is that the period of the motion depends on the range of displacements (figures 6(a) and 7(a)). However, for small displacements around the equilibrium position the motion's period becomes constant,

$$\begin{equation}T\approx 9.20\sqrt{\frac{m}{{E}^{{\ast}}R}}\end{equation} \tag{ 27 }$$

and as a result, the motion can be considered as harmonic. At this point it should be noted that the most typical result arising from the harmonic oscillator is that the period is independent of the amplitude of the oscillation. The aforementioned nonlinear behavior, that becomes linear for small displacements range, is similar to the simple pendulum case. More specifically, when dealing with large oscillations, the complete period of a pendulum is given by an expansion in power series of the initial angle at which the pendulum is dropped [34],

$$\begin{equation}T=2\pi \sqrt{\frac{l}{g}}\left(1+\frac{1}{16}{\theta }_{0}^{2}+\frac{11}{3072}{\theta }_{0}^{4}+\frac{173}{737\enspace 280}{\theta }_{0}^{6}+\cdots \right).\end{equation} \tag{ 28 }$$

In addition, in the presented by this paper case, using a Taylor's series expansion of equation (18) and assuming that the equilibrium position is at y = 0 it is concluded:

$$\begin{equation}T=9.202\enspace 56\sqrt{\frac{m}{{E}^{{\ast}}R}}\left(1+0.0536{y}_{0}+1.8385{y}_{0}^{2}+13.3488{y}_{0}^{3}+\cdots \right).\end{equation} \tag{ 29 }$$

Thus, in this case the complete period of the sphere's motion on the elastic half-space is given by an expansion in power series of the motion's range of displacement. In addition, as it is obvious from the aforementioned Taylor's series expansion, for small values of y₀, $T\approx 9.20\sqrt{\frac{m}{{E}^{{\ast}}R}}$ (the case of harmonic motion as it has been already mentioned).

It must be also noted that the analysis presented in this paper can be also applied in other cases of applying forces that are not proportional to the body's displacement (forces of the several form F = cy^u ). In the above-mentioned cases the resulting differential equation to be solved has the form y'' = mg − cy^u . Those are the cases of forces applied by non-linear springs, by liquids (when a rigid body is partially skunked into the liquid), by half-spaces on rigid bodies (with different geometries), by elastic membranes etc. Thus, future work will focus on the equations of motion that describe the aforementioned cases.

In this paper the nonlinear oscillation of a rigid sphere on an elastic half-space was explored. It was found that a Fourier series function (using two sinusoidal and two cosinusoidal terms terms) perfectly fits to the data of the position of the sphere with respect to time. Moreover, as an additional approximate solution, a function in terms of the range of displacement was also provided. Special attention was given to find solutions in terms of elementary functions in order for the solutions to be easily applied in relative applications. It was also shown that the oscillation's period depends on the range of displacements with an exponential law. In particular, if the range of displacements increases, the period of the motion increases since the average stiffness of the system decreases.

The presented by this paper results are similar to the ones of the oscillation of the simple pendulum. In particular, in both cases, the period of the motion can be described by a power series in terms of the range of displacement and for small oscillation's amplitudes the motion can be considered as harmonic.