Fractional Kelvin-Voigt Model for Beam Vibrations: Numerical Simulations and Approximation by a Classical Model

1. Introduction

A cantilever beam with piezoelectric layers and a tip mass under base excitation is a widely studied model in piezoelectric energy harvesting [1]. Damping in such beams has traditionally been modeled using classical derivative-based approaches [2]. Alternatively, fractional-order models have been employed to describe damping effects more accurately, capturing memory-dependent and nonlocal characteristics of viscoelastic materials [3,4,5].

The Galerkin and Rayleigh-Ritz methods are commonly employed to solve problems involving vibrating beams and other mechanical structures such as rods and plates [3,6]. These methods enable the decomposition of the solution of the governing equation into two components: spatial and temporal. Typically, the temporal solution is predetermined, as it often corresponds to the natural vibration modes of the mechanical structure under consideration or is represented by appropriately chosen, well-established functions.

In this study, the vibration of a beam with a tip mass under base excitation is analyzed. The damping in the system is modeled using a fractional Kelvin-Voigt constitutive equation for a viscoelastic material [7]. The application of the Rayleigh-Ritz method reduces the problem to a set of independent fractional differential equations, each corresponding to a distinct natural vibration mode of the beam. These equations take the following form:

$$\ddot{u}\left(t\right)+{\mu \cdot D}^{\alpha }u\left(t\right)+{\omega }^{2}u\left(t\right)=f\left(t\right),$$

(1)

where $u\left(t\right)$ represents the temporal component of the solution associated with a specific natural vibration mode, $\mu >0$ and $\omega >0$ are fixed real parameters, $0<\alpha <1$ is the order of the fractional derivative $D^{\alpha }$, and $f\left(t\right)$ represents the excitation term corresponding to the given vibration mode.

In the literature, the model described by equation (1) is referred to as the fractional Kelvin-Voigt (FKV) oscillator [8]. It attracted the attention of many researchers [9,10,11,12,13]. Some studies employ the Riemann-Liouville definition of the fractional derivative [9], while others use the Caputo definition [12]. Although analytical solutions to Eq. (1) have been derived [11,14,15], numerical methods are often preferred for solving this equation due to the complexity of the analytical solutions and the computational challenges associated with their practical application. In contrast, classical models with integer-order derivatives do not exhibit such difficulties.

This naturally raises the question: under what conditions can fractional models be replaced by classical models, and how should the parameters of classical models be selected to best approximate the behavior of fractional models? A similar problem has been investigated in [16,17,18,19] using different approaches. The Fourier transform and the response of oscillators to a unit impulse were analyzed in [16]. In [17] the total mechanical energy of fractional operators was used to determine the damping coefficient of a classical model corresponding to a fractional model. The relationship between fractional and classical models was established in [18] by matching their dispersion relations. Finally, in [19] a measure of the incompatibility between fractional and classical models, referred to as the divergence coefficient, was introduced, and optimization methods were employed to refine the parameter selection.

This article presents a method for determining the solution of Eq. (1) using closed-form analytical formulas, complemented by a numerical approach to accelerate the convergence of slowly converging series, specifically the Wynn-epsilon method [20]. Typically, equations of this type are solved using purely numerical techniques [21,22] due to the complexity of their analytical solutions. However, obtaining solutions in a closed form is of significant importance, as they provide a benchmark for validating numerical methods and assessing their accuracy.

To determine the classical model—a damped harmonic oscillator, corresponding to a given FKV oscillator—two approaches are employed. In the first approach, the parameters of the classical model—namely, the undamped angular frequency and the damping ratio—are identified by analyzing the shape of the response graph. This includes examining the distribution of zeros and the positions of local maxima, with the logarithmic decrement method applied to the latter [23]. The second approach utilizes the divergence coefficient introduced in [19]. The parameters of the classical model are estimated by minimizing the divergence coefficient using a genetic algorithm. The results obtained from both approaches are compared and discussed.

The paper is structured as follows. The introduction outlines the current state of research and presents the scope of the study in relation to previous work. In the next section, the fractional Kelvin–Voigt (FKV) oscillator equation is obtained as the temporal part of the motion of a cantilever beam with a tip mass subjected to base excitation. The following section presents analytical solutions to both the homogeneous and nonhomogeneous forms of the FKV equation for Riemann–Liouville and Caputo derivatives, and discusses two methods for identifying the parameters of an equivalent harmonic oscillator. The results section includes validation of the numerical method based on analytical formulas, evaluation of the identified harmonic models, and their performance in the nonhomogeneous case. The final section presents the conclusions of the study.

2. Formulation, Solution, and Computational Methods

2.1. Governing Equation and Solution Using the Rayleigh-Ritz Method

Consider a cantilever beam, clamped at $x=0$, with a tip mass $M_{tip}$​ at its free end $x=L$, subjected to base excitation, as shown in Figure 1. Assume that the relationship between stress and strain in the beam’s material is governed by a constitutive equation known as the fractional Kelvin-Voigt model [7]:

$$\sigma \left(x,z,t\right)=E_0\cdot \epsilon \left(x,z,t\right)+{E_0{\tau }^{\alpha }\cdot D}^{\alpha }\epsilon \left(x,z,t\right),$$

(2)

where $E_0$ is a material constant, $\tau$ is retardation time, and $\alpha$ is the order of the fractional derivative $D^{\alpha }$, with $0<\alpha <1$, which is taken with respect to $t$ only.

According to Euler-Bernoulli’s beam bending theory, the relationship between strain $\epsilon (x,z,t)$ and the transverse displacement $w(x,t)$ of the beam is given by:

$$\epsilon \left(x,z,t\right)=-z\cdot {\displaystyle \frac{{\partial }^{2}w\left(x,t\right)}{\partial x^2}},$$

(3)

The bnding moment $M\left(x,t\right)$ in the beam cross section $C$ is determined as:

$$M\left(x,t\right)={\iint }_{C}z\cdot \sigma \left(x,z,t\right)dydz,$$

(4)

By substituting Eqs. (2)-(4) into the equilibrium equation:

$$\rho A\cdot {\displaystyle \frac{{\partial }^{2}w\left(x,t\right)}{\partial t^2}}={\displaystyle \frac{{\partial }^{2}M\left(x,t\right)}{\partial x^2}},$$

(5)

one obtains the governing equation for the transverse displacement $w\left(x,t\right)$:

$$\rho A\cdot {\displaystyle \frac{{\partial }^{2}w\left(x,t\right)}{\partial t^2}}+E_{0}I{\displaystyle \frac{{\partial }^{4}w\left(x,t\right)}{\partial x^4}}+E_0{\tau }^{\alpha }I{D}^{\alpha }{\displaystyle \frac{{\partial }^{4}w\left(x,t\right)}{\partial x^4}}=0$$

(6)

where $I$ is the moment of inertia of the cross-section with respect to the bending axis $y$, $A=W\cdot H$ is the cross-sectional area, $\rho$ is a density of beam’s material. The absolute transverse displacement $w\left(x,t\right)$ can be expressed as:

$$w\left(x,t\right)=w_{rel}\left(x,t\right)+w_b\left(x,t\right),$$

(7)

where $w_{rel}\left(x,t\right)$ is the displacement relative to a moving coordinate system with its origin at the clamped end of the beam, and $w_b\left(x,t\right)$ is the base displacement. It is assumed that:

$$w_b\left(x,t\right)=A_0\mathrm{s}\mathrm{i}\mathrm{n}\left(\mathsf{\Omega }t\right),$$

(8)

where $A_0$ and $\mathsf{\Omega }$ are amplitude and angular frequency of the base excitation respectivelly. Substituting Eqs. (7) and (8) into Eq. (6) yields:

$$\rho A\cdot {\displaystyle \frac{{\partial }^2{w}_{rel}\left(x,t\right)}{\partial t^2}}+E_{0}I{\displaystyle \frac{{\partial }^4{w}_{rel}\left(x,t\right)}{\partial x^4}}+E_0{\tau }^{\alpha }I{D}^{\alpha }{\displaystyle \frac{{\partial }^4{w}_{rel}\left(x,t\right)}{\partial x^4}}=\rho A{A}_0{\mathsf{\Omega }}^2\mathrm{s}\mathrm{i}\mathrm{n}\left(\mathsf{\Omega }t\right),$$

(9)

The boundary conditions for displacement and angle at clamped end $x=0$ are:

$$w_{rel}\left(0,t\right)=0, {\displaystyle \frac{\partial w_{rel}}{\partial x}}\left(0,t\right)=0.$$

(10)

The boundary condition for the bending moment at the free end $x=L$ is:

$$E_{0}I {\displaystyle \frac{{\partial }^2{w}_{rel}}{\partial x^2}}\left(L,t\right)+E_0{\tau }^{\alpha }{D}^{\alpha }{\displaystyle \frac{{\partial }^2{w}_{rel}}{\partial x^2}}\left(x,L\right)+I_{tip}{\displaystyle \frac{{\partial }^2{w}_{rel}}{\partial t^2}}\left(L,t\right)=0,$$

(11)

The boundary condition for the transverse force at the free end $x=L$ is the following:

$$E_{0}I {\displaystyle \frac{{\partial }^3{w}_{rel}}{\partial x^3}}\left(L,t\right)+E_0{\tau }^{\alpha }{D}^{\alpha }{\displaystyle \frac{{\partial }^3{w}_{rel}}{\partial x^3}}\left(x,L\right)+M_{tip}{\displaystyle \frac{{\partial }^2{w}_{rel}}{\partial t^2}}\left(L,t\right)=0,$$

(12)

where $M_{tip}$ is the mass of the body attached to the free end of the beam and $I_{tip}$ is its moment of inertia about $x=0$ (for point mass $M_{tip}$ it would be equal to $I_{tip}=L^2{\cdot M}_{tip}$).

To obtain an approximate solution to Eq. (9), the Rayleigh-Ritz method is employed. In this approach, the transverse displacement $w_{rel}(x,t)$ is assumed to be a series expansion in terms of trial functions that satisfy the boundary conditions:

$$w_{rel}\left(x,t\right)={\displaystyle \sum _{i=1}^{\infty }}{\phi }_i\left(x\right)\cdot u_i\left(t\right)$$

(13)

where ${\phi }_i\left(x\right)$ are the trial functions representing the mode shapes of the beam, and $u_i\left(t\right)$ are the generalized coordinates that describe the time-dependent motion of the system. The mode shapes of a cantilever beam are given by [1,24]:

$${\phi }_i\left(x\right)={\displaystyle \frac{1}{\sqrt{\rho A L}}}\left(\mathrm{c}\mathrm{o}\mathrm{s}\left({\displaystyle \frac{{\lambda }_i}{L}}x\right)-\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}\left({\displaystyle \frac{{\lambda }_i}{L}}x\right)+{\sigma }_i\left(\mathrm{s}\mathrm{i}\mathrm{n}\left({\displaystyle \frac{{\lambda }_i}{L}}x\right)-\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}\left({\displaystyle \frac{{\lambda }_i}{L}}x\right)\right)\right)$$

(14)

where:

$${\sigma }_i={\displaystyle \frac{\sin {\lambda }_i-\sinh {\lambda }_i}{\cos {\lambda }_i-\cosh {\lambda }_i}}$$

(15)

and dimensionless eigenvalues ${\lambda }_i$ are determined from the characteristic equation:

$$1+\cos \lambda \cosh \lambda =0$$

(16)

Eigenfunctions defined in Eq. (14) are orthogonal, which means that for $i\ne j$, the following conditions hold:

$${\displaystyle \underset{0}{\overset{L}{\int }}}{\phi }_i\left(x\right)\cdot {\phi }_j\left(x\right)dx=0, {\displaystyle \underset{0}{\overset{L}{\int }}}{\phi }_i\left(x\right)\cdot {\displaystyle \frac{d^4{\phi }_j\left(x\right)}{d{x}^4}}dx=0$$

(17)

for $i=j$, these integrals are nonzero.

Multiplying Eq. (9) by ​${\phi }_k\left(x\right)$ and integrating over the range $[0,L]$, we obtain, based on the orthogonality conditions of the natural vibration modes (17), an equation of the form (1), as previously introduced, where the function $f\left(t\right)$ is defined as:

$$f\left(t\right)=F_0\sin \left(\mathsf{\Omega }t\right),$$

(18)

For some $F_0$ associated with the $k$-th vibration mode of the beam.

2.2. Solution of the Temporal Equation and Numerical Methods

In this section, we present closed-form solutions for Eq. (1) and discuss numerical methods for their efficient computation.

We begin with the homogeneous equation for $f\left(t\right)=0$:

$$\ddot{u}\left(t\right)+{\mu \cdot D}^{\alpha }u\left(t\right)+{\omega }^{2}u\left(t\right)=0,$$

(19)

We consider the above equation with initial conditions: $x\left(0\right)=x_0$, $\dot{x}\left(0\right)=v_0$.

Based on Theorem 5.2 and Proposition 5.2 from [14], and after applying the necessary transformations, the solution of Eq. (19) for the Riemann–Liouville fractional derivative takes the form::

$$x_{RL}\left(t\right)=x_0{x}_1^{\left(RL\right)}\left(t\right)+v_0{x}_2^{\left(RL\right)}\left(t\right),$$

(20)

where:

$$x_1^{\left(RL\right)}\left(t\right)=\sum _{n=0}^{\infty }\sum _{k=0}^{\infty }{\displaystyle \frac{{\left(-1\right)}^{n+k}\left(\begin{array}{c}n+k\\ n\end{array}\right){\omega }^{2n}{\mu }^k}{\mathsf{\Gamma }(2\left(n+k\right)-\alpha k+1)}}{t}^{2\left(n+k\right)-\alpha k}$$

(21)

$$x_2^{\left(RL\right)}\left(t\right)=t\sum _{n=0}^{\infty }\sum _{k=0}^{\infty }{\displaystyle \frac{{\left(-1\right)}^{n+k}\left(\begin{array}{c}n+k\\ n\end{array}\right){\omega }^{2n}{\mu }^k}{\mathsf{\Gamma }(2\left(n+k\right)-\alpha k+2)}}{t}^{2\left(n+k\right)-\alpha k}$$

(22)

Similarly, applying Theorem 5.13 and Proposition 5.8 from [14], the solution for Eq. (19) in the case of the Caputo fractional derivative takes the form:

$$x_C\left(t\right)=x_0{x}_1^{\left(C\right)}\left(t\right)+v_0{x}_2^{\left(RL\right)}\left(t\right),$$

(23)

where:

$$x_1^{\left(C\right)}\left(t\right)=x_1^{\left(RL\right)}\left(t\right)+\mu t^{2-\alpha }\sum _{n=0}^{\infty }\sum _{k=0}^{\infty }{\displaystyle \frac{{\left(-1\right)}^{n+k}\left(\begin{array}{c}n+k\\ n\end{array}\right){\omega }^{2n}{\mu }^k}{\mathsf{\Gamma }(2\left(n+k\right)-\alpha (k+1)+3)}}{t}^{2\left(n+k\right)-\alpha k},$$

(24)

For the nonhomogeneous Eq. (1), solutions can be obtained for both Riemann-Liouville and Caputo fractional derivatives using Proposition 5.5 and Proposition 5.11, respectively, in the form:

$$x\left(t\right)=x_{hom}\left(t\right)+{\displaystyle \underset{0}{\overset{t}{\int }}}\left(t-\tau \right)\cdot G\left(t-\tau \right)\cdot f\left(\tau \right)d\tau ,$$

(25)

where the function $G\left(t\right)$ is defined as:

$$G\left(t\right)=\sum _{n=0}^{\infty }\sum _{k=0}^{\infty }{\displaystyle \frac{{\left(-1\right)}^{n+k}\left(\begin{array}{c}n+k\\ n\end{array}\right){\omega }^{2n}{\mu }^k}{\mathsf{\Gamma }(2\left(n+k\right)-\alpha k+2)}}{t}^{2\left(n+k\right)-\alpha k}$$

(26)

The function $x_{hom}\left(t\right)$ represents the solution of the homogeneous Eq. (19) for either the Riemann-Liouville or Caputo derivative, as given in Eqs. (20) and (23), respectively.

Direct computation of $x_1^{\left(RL\right)}\left(t\right), x_2^{\left(RL\right)}\left(t\right), x_1^{\left(C\right)}\left(t\right)$ and $G\left(t\right)$ by truncating the series in Eqs. (21), (22), (24), and (26) is only effective for small values of t. For larger $t$, naive truncation leads to significant numerical errors due to cancellation and round-off in the summation of highly oscillatory terms. To overcome this limitation and extend the computation to a broader range of $t$, we employed the Wynn-$\epsilon$ algorithm [20], as implemented in Mathematica’s built-in function SequenceLimit[]. The integral in Eq. (25) was evaluated using the trapezoidal rule. Our calculations were performed using Mathematica v. 12.1.

Applying the Wynn-$\epsilon$ method requires transforming the double series into a single series. Since the original series are absolutely convergent, their terms can be grouped and summed in an arbitrary order. Here, summation was performed over $N=n+k$, varying from zero to infinity. The computations required high floating-point precision, and all model parameters were converted into symbolic fractions to facilitate exact symbolic calculations. While this approach eliminated rounding errors, it significantly increased computation time.

The computational domain $0\le t\le T$ was chosen to approximately span ten fundamental periods of the homogeneous solution $x_{hom}\left(t\right)$. To determine $T$, the fundamental frequency ${\omega }_0$ of the corresponding harmonic oscillator was estimated using the relation derived in [19]:

$${\omega }_0=\sqrt{{\omega }^2+{\mu }^{{\displaystyle \frac{2}{2-\alpha }}}},$$

(27)

The final time was then set as $T=10\cdot T_0$, with $T_0={\displaystyle \frac{2\pi }{{\omega }_0}}$.

The time step for evaluating the integral in Eq. (25) using the trapezoidal rule was chosen as $\mathsf{\Delta }t=\min \{0.01\cdot T_0, 0.01\cdot {\displaystyle \frac{2\pi }{\mathsf{\Omega }}}\}$.

2.3. Identification of Equivalent Harmonic Oscillator Parameters: Geometric and Optimization Approaches

The classical damped harmonic oscillator is well known to be governed by equation:

$$\ddot{x}\left(t\right)+2\xi {\omega }_0\cdot \dot{x}\left(\mathrm{t}\right)+{\omega }_0^{2}x\left(t\right)=0$$

(28)

subject to the initial conditions $x\left(0\right)=x_0$, $\dot{x}\left(0\right)=v_0$. Here, ${\omega }_0$​ denotes the undamped natural angular frequency of the system, and $\xi$ is the damping ratio. The solution to Eq. (28) can be written in the form:

$$x\left(t\right)=x_0\cdot x_1^{\left(HO\right)}\left(t\right)+v_0\cdot x_2^{\left(HO\right)}\left(t\right),$$

(29)

where the functions $x_1^{\left(HO\right)}\left(t\right)$ and $x_2^{\left(HO\right)}\left(t\right)$ are expressed in terms of exponential and trigonometric functions.

Our aim is to determine an equivalent harmonic oscillator represented by Eq. (28) that approximates the behavior of the FKV oscillator, defined by Eq. (19). Solutions of the Eq. (19) for the case of the Riemann-Liouville and Caputo fractional derivatives, given by formulas (20) and (23) respectively, share a common component corresponding to the initial condition $\dot{x}\left(0\right)=v_0$​. Consequently, if a harmonic oscillator can be identified for the case of the Caputo fractional derivative, it remains the same for the case of Riemann-Liouville fractional derivative. In this work, we choose to focus primarily on the Caputo fractional derivative. This is motivated by the fact that, unlike the Riemann–Liouville derivative, the Caputo formulation allows for the use of physically meaningful initial conditions expressed in terms of the function and its integer-order derivatives, such as position and velocity.

2.3.1. Geometric Approach to Parameter Identification

As a result of the numerical method of solving equation (19), we have the values of the components $x_1^{\left(C\right)}\left(t\right)$ and $x_2^{\left(RL\right)}\left(t\right)$ defined in Eqs. (24) and (22), respectively, at discrete moments of time $t_i=i\cdot \mathsf{\Delta }t$, for some $\mathsf{\Delta }t>0$, $i=\mathrm{0,1},\dots ,N$. Here we discuss how, based on these data, to obtain approximations of the parameters ${\omega }_0$ and $\xi$ of a harmonic oscillator that approximates a given FKV oscillator.

Let the sequences $z_j^{\left(1\right)}, j=\mathrm{1,2},\dots ,M_1$ and $m_k^{\left(1\right)}, k=\mathrm{1,2},\dots ,K_1$ represent, respectively, the approximate locations of the zeros and the arguments corresponding to local maxima of the function $x_1^{\left(C\right)}\left(t\right)$. Similarly, the sequences $z_j^{\left(2\right)}, j=\mathrm{1,2},\dots ,M_2$ and $m_k^{\left(2\right)}, k=\mathrm{1,2},\dots ,K_2$ are defined for the function $x_2^{\left(RL\right)}\left(t\right)$. The approximate parameters of the harmonic oscillator corresponding to the given FKV oscillator are determined from the following formulas:

$${\omega }_0={\displaystyle \frac{\pi }{\frac{1}{M_1-1}\sum _{i=1}^{M_1-1}(z_{i+1}^{\left(1\right)}-z_i^{\left(1\right)})+\frac{1}{M_2-1}\sum _{i=1}^{M_2-1}(z_{i+1}^{\left(2\right)}-z_i^{\left(2\right)})}},$$

(30)

$$\xi ={\displaystyle \frac{1}{2\pi }}\left({\displaystyle \frac{1}{K_1-1}}\sum _{i=1}^{K_1-1}\mathit{ln}{\displaystyle \frac{m_i^{\left(1\right)}}{m_{i+1}^{\left(1\right)}}}+{\displaystyle \frac{1}{K_2-1}}\sum _{i=1}^{K_2-1}\mathit{ln}{\displaystyle \frac{m_i^{\left(2\right)}}{m_{i+1}^{\left(2\right)}}} \right),$$

(31)

In the denominator of equation (30), we approximate half the fundamental period for the harmonic oscillator we are looking for, so dividing $\pi$ by half the period give an approximation of the circular frequency ${\omega }_0$. In equation (31), we use the concept of logarithmic damping decrement to determine $\xi$.

2.3.2. Genetic Algorithm Approach to Parameter Identification

To assess the quality of fit between a given FKV oscillator and a classical harmonic oscillator, we adopted a measure of model similarity defined in [19] called the divergence coefficient:

$$I\left(x_{FKV}\left(t\right),x_{HO}\left(t\right)\right)≔\max \left\{{\displaystyle \frac{{‖x_1^{\left(FKV\right)}\left(t\right)-x_1^{\left(HO\right)}\left(t\right)‖}_{\infty }}{{‖x_1^{\left(FKV\right)}\left(t\right)‖}_{\infty }}},{\displaystyle \frac{{‖x_2^{\left(FKV\right)}\left(t\right)-x_2^{\left(HO\right)}\left(t\right)‖}_{\infty }}{{‖x_2^{\left(FKV\right)}\left(t\right)‖}_{\infty }}} \right\},$$

(32)

where ${‖\cdot ‖}_{\infty }$ is supremum norm on interval$0\le t\le T$for some $T>0$ [25], $x_{FKV}\left(t\right)$ and $x_{HO}\left(t\right)$ are defined in (23) and (28) respectively. Divergence coefficient has an interpretation as the maximum relative difference between components of the functions $x_{FKV}\left(t\right)$ and $x_{HO}\left(t\right)$. Function $x_{FKV}\left(t\right)$ is fixed and we look for the function $x_{HO}\left(t\right)$ that minimizes the divergence coefficient (32). In this case we write $I\left(x_{HO}\left(t\right)\right)$.

Let the FKV oscillator with parameters $\mu ,\alpha ,\omega$ be fixed. We use the genetic algorithm to find the minimum value of the divergence coefficient and the parameters ${\omega }_0$ and $\xi$ of harmonic oscillator, which correspond to the given FKV oscillator. The calculations were carried out using MATLAB's ‘Global Optimisation’ toolbox. The objective function, the divergence coefficient, was minimized using MATLAB's built-in function ga. Although the ga function allows adjusting many parameters, after experimenting with various settings, we decided to use its default options. Changes in the parameters typically did not significantly improve the result's quality, but often substantially increased computation time.

3. Results

3.1. Verification of the Wynn-Epsilon Method for Solving the FKV Oscillator Equation

To validate the numerical method used to evaluate the series solutions of Eq. (19), we consider the special case $\alpha =1$. In this case, Eq. (19) reduces to the classical damped harmonic oscillator equation (28), with parameters given by $\xi =\mu /2\omega$ and ${\omega }_0=\omega$. Accordingly, we substitute $\alpha =1$ into Eqs. (21), (22), and (24), and compare the results obtained using the Wynn-epsilon method with the known analytical functions $x_1^{\left(HO\right)}\left(t\right)$ and $x_2^{\left(HO\right)}\left(t\right)$ appearing in Eq. (29). This comparison allows us to verify the correctness and numerical accuracy of the proposed approach.

Figure 2, Figure 3 and Figure 4 present a comparison between the analytical solution components of the classical harmonic oscillator, $x_1^{\left(HO\right)}\left(t\right)$ and $x_2^{\left(HO\right)}\left(t\right)$ (solid black lines), and the corresponding approximations obtained using the Wynn-epsilon method for the fractional model with Riemann–Liouville (red dashed lines) and Caputo (blue dashed lines) derivatives.

Subfigures labeled (a) show the components associated with the initial displacement condition $x\left(0\right)=x_0$​, while subfigures labeled (b) correspond to the initial velocity condition $\dot{x}\left(0\right)=v_0$​. It should be noted that in the subfigures (b), only two curves are plotted, since the solution components related to the initial velocity are identical for both the Riemann–Liouville and Caputo cases.

A very good agreement can be observed between the Wynn-$\epsilon$ approximation of the $x_1^{\left(C\right)}\left(t\right)$ and the function $x_1^{\left(HO\right)}\left(t\right)$ in Figure 2, Figure 3 and Figure 4(a), with the maximum relative absolute difference being on the order of ${10}^{-16}$. A similar observation holds for the Wynn-$\epsilon$ approximation $x_2^{\left(RL\right)}\left(t\right)$ and the function $x_2^{\left(HO\right)}\left(t\right)$ shown in Figure 2, Figure 3 and Figure 4(b).

In contrast, noticeable discrepancies occur between $x_1^{\left(RL\right)}\left(t\right)$ (red dashed lines in Figure 2, Figure 3 and Figure 4(a)) and $x_1^{\left(HO\right)}\left(t\right)$, with maximum relative differences reaching approximately 0.36 in Figure 2(a), 0.55 in Figure 3(a), and 0.09 in Figure 4(a). These deviations are related to the incompatibility of initial conditions in the Riemann–Liouville formulation. However, it can be shown that the difference $x_1^{\left(C\right)}\left(t\right)-x_1^{\left(RL\right)}\left(t\right)$ tends to zero as $t\to \infty$, indicating that the solutions converge asymptotically.

To validate the numerical solutions of the nonhomogeneous fractional Kelvin–Voigt oscillator equation with zero initial conditions, we once again consider the limiting case $\alpha =1$. In this case, the equation reduces to the classical nonhomogeneous damped harmonic oscillator:

$$\ddot{x}\left(t\right)+2\xi {\omega }_0\cdot \dot{x}\left(\mathrm{t}\right)+{\omega }_0^{2}x\left(t\right)=\mathrm{f}\left(\mathrm{t}\right)$$

(33)

where ${\omega }_0=\omega$, $\xi =\mu /2\omega$, and the forcing function $f\left(t\right)$ is defined in Eq. (18). We denote the solution of Eq. (33) under zero initial conditions by $x_{HO}\left(t\right)$. This solution serves as a reference for assessing the accuracy of the numerical results obtained for the fractional model with $\alpha =1$.

Let $x_{FKV}\left(t\right)$ denote the solution of Eq. (1) with zero initial conditions for $\alpha =1$. In this case, the homogeneous component vanishes, and the solution is given by the convolution integral appearing in Eq. (25).

To illustrate the accuracy of the proposed numerical approach, Figure 5 presents a comparison between the exact solution $x_{HO}\left(t\right)$ of the classical oscillator equation (33) and the approximation $x_{FKV}\left(t\right)$ obtained using the Wynn-epsilon method. The comparison is carried out for two selected fractional Kelvin–Voigt oscillators subjected to the external excitation defined in Eq. (18), with parameters fixed as $F_0=1$ and $\mathsf{\Omega }=5$.

The agreement between the exact solution and the Wynn-epsilon approximation is very good. The maximum relative difference between the two curves is on the order of ${10}^{-6}$.

3.2. Identification of Equivalent Harmonic Oscillator Parameters: Geometric and Optimization Approaches

In this section, we compare two methods for identifying the parameters of a classical harmonic oscillator that best approximates the response of a given fractional Kelvin–Voigt (FKV) oscillator. For each method, denoted by the superscript $\left(m\right)$, we define the identified optimal parameters of the equivalent harmonic oscillator as ${\omega }_0^{\left(m\right)}$ and ${\xi }^{\left(m\right)}$, corresponding to the undamped natural frequency and the damping ratio, respectively.

The resulting optimal approximation is denoted by $x_{HO}^{\left(m\right)}\left(t\right)$, and its accuracy is evaluated using formula ${\epsilon }_m=I\left(x_{HO}^{\left(m\right)}\right(t\left)\right)$, where $I(\cdot )$ denotes the divergence coefficient defined in Eq. (32), which quantifies the discrepancy between the fractional and harmonic oscillator responses.

Table 1 presents a comparison of parameter identification results for several fractional Kelvin–Voigt (FKV) with parameters $\left(\mu ,\alpha ,\omega \right)$ shown in the first column. For each case, the identified parameters of the equivalent harmonic oscillator obtained using the geometric method (second column) and the genetic algorithm (third column) are listed, along with the corresponding values of the divergence coefficient ${\epsilon }_m$, which quantifies the discrepancy between the fractional and classical models.

Analysis of the results reveals that the genetic algorithm generally yields lower values of the divergence coefficient ${\epsilon }_m$, indicating better agreement between the models. However, this is not always the case, as illustrated by rows 1 and 7 of Table 1, where the geometric method performs slightly better. In certain cases (rows 7 and 10), the divergence coefficient reaches relatively high values, suggesting that the classical harmonic model may not always provide a satisfactory approximation of the fractional dynamics.

Figure 6, Figure 7, Figure 8 and Figure 9 6–9 present a comparison between the solution components of selected FKV oscillators and those of classical harmonic oscillators identified using two different methods: a geometric approach and a genetic algorithm. Each figure is composed of two parts (a) and (b). Figures (a) show comparison of the displacement-related components $x_1^{\left(C\right)}\left(t\right)$ and $x_1^{\left(HO\right)}\left(t\right)$, while figures (b) comparison of the velocity-related components $x_2^{\left(RL\right)}\left(t\right)$ and $x_2^{\left(HO\right)}\left(t\right)$. Here, $x_1^{\left(HO\right)}\left(t\right)$ and $x_2^{\left(HO\right)}\left(t\right)$ correspond to components of the harmonic oscillator solution ​$x_{HO}^{\left(m\right)}\left(t\right)$ computed using the optimal parameters obtained via the respective identification method.

The results illustrated in Figure 6, Figure 7, Figure 8 and Figure 9 provide additional insight into the effectiveness of both identification methods in approximating the behavior of various FKV oscillators. The observations are consistent with the divergence values ${\epsilon }_m$​ reported in Table 1.

In Figure 6, corresponding to parameters $\left(\mu ,\alpha ,\omega \right)=\left(\mathrm{1,0.1,5}\right)$, both methods yield harmonic oscillators with relatively high divergence values (around 0.08), despite the overall good visual agreement of the solution components. Deviations are mainly noticeable near the extrema of the curves, as seen in part (a) of the figure.

Figure 7 shows results for $\left(\mu ,\alpha ,\omega \right)=\left(\mathrm{1,0.5,10}\right)$, where the divergence coefficients are somewhat different—0.037 for the geometric method and 0.026 for the genetic algorithm—yet in both cases the visual agreement between the fractional and classical components is excellent.

In Figure 8, with parameters $\left(\mu ,\alpha ,\omega \right)=\left(\mathrm{10,0.5,10}\right)$, a more pronounced discrepancy is observed: the geometric approach produces a relatively large value of the ${\epsilon }_m$, while the genetic algorithm yields a smaller value. This case corresponds to a system with strong damping, which may contribute to the increased difficulty in achieving an accurate harmonic approximation.

Finally, Figure 9, representing $\left(\mu ,\alpha ,\omega \right)=\left(\mathrm{12,0.7,10}\right)$, shows that neither method succeeds in producing a good match. The divergence coefficient exceeds 0.1 in both cases, reaching as high as 0.25 for the geometric approach. This poor agreement, clearly visible in both components, suggests that for this configuration, a classical harmonic oscillator may not adequately capture the dynamics of the fractional system.

3.3. Testing Harmonic Models in the Nonhomogeneous Case

In this subsection, we investigate the performance of the identified harmonic oscillator models in the context of the nonhomogeneous fractional Kelvin–Voigt equation (Eq. (1)) with zero initial conditions. The solution of this equation, given by Eq. (25), reduces to a convolution integral, as the homogeneous part vanishes. For comparison, we consider the classical damped harmonic oscillator (Eq. (33)) with zero initial conditions and identical sinusoidal forcing, whose solution can be expressed in closed form using elementary functions such as sine, cosine, and exponential terms.

The harmonic oscillator parameters used in this comparison were obtained using the identification methods described in the previous subsection, for fixed values of the FKV model parameters $(\mu ,\alpha ,\omega )$. Two representative FKV oscillators were selected for this analysis: $\left(\mathrm{1,0.1,5}\right)$ and $\left(\mathrm{1,0.5,10}\right)$. The corresponding harmonic oscillator parameters are listed in Table 1 (rows 1 and 5, respectively).

To quantitatively assess the agreement between the solutions of the fractional model $x_{FKV}\left(t\right)$ and the classical harmonic oscillator $x_{HO}\left(t\right)$, we use a modified version of the divergence coefficient (32) introduced earlier:

$$I\left(x_{FKV}\left(t\right),x_{HO}\left(t\right)\right)≔{\displaystyle \frac{{‖x_{FKV}\left(t\right)-x_{HO}\left(t\right)‖}_{\infty }}{{‖x_{FKV}\left(t\right)‖}_{\infty }}}$$

(34)

Specifically, when the indicator defined in Eq. (34) is evaluated using the solution $x_{HO}\left(t\right)$ obtained with parameters identified in the previous subsection (using either the geometric method or the genetic algorithm), we denote its value by ${\epsilon }_m$​.

Figure 10 and Figure 11 present a comparison between the time-domain responses of the fractional oscillator $x_{FKV}\left(t\right)$ and the corresponding harmonic oscillator approximations $x_{HO}\left(t\right)$, obtained using the geometric method (red dashed lines) and the genetic algorithm (blue dashed lines). The fractional response is plotted in black as a reference.

For each selected FKV oscillator, two scenarios are considered: one in which the forcing frequency $\mathsf{\Omega }$ is far from the resonance frequency of the harmonic model, and another in which $\mathsf{\Omega }$ is chosen to match the resonance frequency of the harmonic model (we set $\mathsf{\Omega }=\omega$). This setup allows us to examine how the accuracy of the harmonic approximation depends on the proximity to resonance.

The results shown in Figure 10 and Figure 11 highlight the influence of the forcing frequency on the accuracy of harmonic oscillator approximations in the nonhomogeneous case. Each figure compares the solution of the FKV equation with those of two harmonic oscillators identified using the geometric method and the genetic algorithm, respectively.

For the FKV oscillator with parameters $\left(\mu ,\alpha ,\omega \right)=\left(\mathrm{1,0.1,5}\right)$, subfigure 10(a) illustrates the resonance case with $\mathsf{\Omega }=5$, where substantial discrepancies are observed: the divergence coefficients reach ${\epsilon }_m=0.3875$ for the geometric method and ${\epsilon }_m=0.3408$ for the genetic algorithm. In contrast, subfigure 10(b), corresponding to the off-resonance case with $\mathsf{\Omega }=10$, shows very good agreement, with ${\epsilon }_m$ dropping to $0.0168$ and $0.0036$, respectively.

A similar trend is visible for the oscillator defined by $\left(\mu ,\alpha ,\omega \right)=\left(\mathrm{1,0.5,10}\right)$. In this case, the off-resonance response (subfigure 11(a)), $\mathsf{\Omega }=10$ yields relatively low divergence values: ${\epsilon }_m=0.0293$ for the geometric method and ${\epsilon }_m=0.0284$ for the genetic algorithm. However, the resonance case shown in subfigure 11(b), with $\mathsf{\Omega }=10$, leads to a complete breakdown of the harmonic approximation, with divergence values near to 1 for both methods.

These results demonstrate that while harmonic oscillator models can accurately reproduce the fractional response under off-resonance excitation, they may fail dramatically near resonance. This is due to the high sensitivity of the system’s steady-state behavior in the resonant regime, where subtle differences in the underlying dynamics—especially those introduced by the fractional-order operator—are strongly amplified and no longer well captured by classical models.

To further investigate the sensitivity of the harmonic approximation with respect to the forcing frequency, we compute the divergence coefficient ${\epsilon }_m$​ for a range of excitation frequencies. Specifically, we vary the ratio $\omega /\mathsf{\Omega }$ from $0.05$ to $1.95$ in steps of $0.05$, while keeping the parameters of the FKV oscillator fixed. This analysis is performed for both identification methods and for the two representative FKV oscillators considered earlier.

The plots in Figure 12 reveal that the divergence coefficient ${\epsilon }_m$ reaches a sharp peak in the vicinity of $\omega /\mathsf{\Omega }=1$, confirming that the accuracy of the harmonic approximation deteriorates significantly near resonance.

For the first FKV oscillator (Figure 12(a)), the peak is less pronounced and slightly shifted away from $\omega /\mathsf{\Omega }=1$. This shift reflects the fact that the actual resonance frequency of the corresponding harmonic oscillator does not exactly match the parameter $\omega$ of the FKV model.

In contrast, for the second oscillator (Figure 12(b)), the resonance frequency of the harmonic model is close to $\omega =10$, resulting in a stronger and more centered peak. These observations highlight the limitations of harmonic approximations in capturing resonant behavior in fractional-order systems.

4. Conclusions

This paper considered the dynamic response of a cantilever beam with a tip mass subjected to base excitation. The beam material was assumed to obey a fractional Kelvin–Voigt constitutive law. By applying modal decomposition, a time-domain equation governing the selected mode shape was derived. It was shown that the resulting model corresponds to the FKV oscillator equation.

A method for computing accurate time-domain solutions based on analytical series expressions was presented. The series, evaluated using the Wynn-epsilon algorithm, provide high accuracy at the cost of significant computational effort. While purely numerical methods may be more efficient in practical applications, the ability to compute highly accurate reference solutions based on analytical formulas is of great importance, as it enables reliable validation and benchmarking of approximate numerical schemes. The method was validated against exact solutions for the case $\alpha =1$. It was also observed that, for the Riemann–Liouville derivative, the component associated with the initial displacement condition deviates from the exact solution. This discrepancy is not due to numerical inaccuracy but rather reflects the known issues related to initial conditions for this type of derivative.

For a given FKV oscillator, two methods were proposed for identifying the parameters of an equivalent classical harmonic oscillator: a geometric approach and a genetic algorithm. Both methods yielded comparable accuracy in most cases, although the genetic algorithm generally produced better results.

Finally, the identified harmonic oscillators were applied to the nonhomogeneous version of the problem with sinusoidal forcing. The harmonic approximation was shown to be accurate when the forcing frequency was far from resonance, but its performance degraded significantly near the resonance frequency of the harmonic model.