## Abstract

Fractional models and their parameters are sensitive to intrinsic microstructural changes in anomalous materials. We investigate how such physics-informed models propagate the evolving anomalous rheology to the nonlinear dynamics of mechanical systems. In particular, we study the vibration of a fractional, geometrically nonlinear viscoelastic cantilever beam, under base excitation and free vibration, where the viscoelasticity is described by a distributed-order fractional model. We employ Hamilton's principle to obtain the equation of motion with the choice of specific material distribution functions that recover a fractional Kelvin–Voigt viscoelastic model of order *α*. Through spectral decomposition in space, the resulting time-fractional partial differential equation reduces to a nonlinear time-fractional ordinary differential equation, where the linear counterpart is numerically integrated through a direct L1-difference scheme. We further develop a semi-analytical scheme to solve the nonlinear system through a method of multiple scales, yielding a cubic algebraic equation in terms of the frequency. Our numerical results suggest a set of *α*-dependent anomalous dynamic qualities, such as far-from-equilibrium power-law decay rates, amplitude super-sensitivity at free vibration, and bifurcation in steady-state amplitude at primary resonance.

## 1 Introduction

Nonlinearities are intrinsic in real physical systems, arising from multiple sources, such as changes in geometry, material response (e.g., ageing), and boundary effects (e.g., emergence of boundary-layers/shock). We focus on the analysis of nonlinear systems where anomalous dynamics arise from nonlocal/history effects. Despite the existing “nearly-pure” systems, where standard features evolve to anomalous qualities, e.g., laminar-to-turbulent flows [1,2] and dislocation pile-up in plasticity [3], in our study, the anomaly source is from the employment of extraordinary materials.

Power-law rheology is a characteristic of a wide range of anomalous materials. This complex response exhibits macroscopic *memory-effects* through single-to-multiple power-law relaxation/creep [4] and dynamic storage/dissipation in visco-elasticity [5]. Such power laws are multiscale fingerprints of spatio-temporal subdiffusive processes in fractal-like microstructures, where the mean squared displacement of constituents/defects scales nonlinearly in time as $\u2329\Delta r\u232a2\u221dt\alpha $ [6]. When subject to mechanical loads, such materials undergo microstructural changes, e.g., rearrangement/unfolding of polymer networks/chains [4], plastic stretching/buckling of microfibers [7], formation, arresting, relaxation of dislocations [8], which propagate the evolving microrheology to the larger scales. Classical (integer-order) viscoelastic models accurately fit exponential data with a limited number of relaxation times [9]. However, a large number of relaxation time is usually required to estimate complex hereditary behaviors observed for a broad class of anomalous (nonexponential) materials. This yields high-dimensional parameter spaces, worsening the conditioning of already ill-posed parameter estimations [10]. In addition, multi-exponential approximations are mere truncations of power-law relaxation [11], only providing accurate results for short times, thus lacking predictability and requiring recalibration for multiple time-scales [12].

Fractional differential equations (FDEs) are predictive tools for anomalous materials across multiple time-scales. Nutting [13] and Gemant [14] showed that power law kernels are more descriptive of creep/relaxation. Bagley and Torvik [15] proposed a link between fractional viscoelasticity and polymer dynamics through dynamic moduli. The building block of fractional viscoelasticity is the so-called *Scott-Blair* (SB) element with fractional order $0<\alpha <1$, that interpolates between Hookean $(\alpha \u21920)$ and Newtonian $(\alpha \u21921)$ elements. Distinct arrangements of SB elements model multiple experimentally observed power laws through multiterm FDEs. Such flexible and compact mathematical tool led to multidisciplinary developments, e.g., in bio-engineering [16,17], visco-elasto-plasticity [18,19], among others [20–22]. The most general forms of viscoelasticity are described by distributed order differential equations (DODEs) [23,24], where fractional-order distributions (distributions of SB elements) code evolving, heterogeneous multiscale material properties. DODEs were employed for anomalous diffusion in Refs. [25–27] with applications also in control theory and signal processing [28,29], vibration [30], frequency domain analysis [31], and uncertainty quantification [32,33]. We refer the reader to Ref. [34] for a thorough review on applications of DODEs.

Regarding the dynamics of anomalous beams, Łabȩdzki et al. [35] investigated the resonance of Euler–Bernoulli piezo-electric beams by introducing fractional derivatives in the equation of motion, and solved the strong form of the system using a Rayleigh–Ritz method. Ansari et al. [36] studied the free vibration of a nonlocal, fractional Kelvin–Voigt (KV) Euler–Bernoulli nanobeam. A direct Ritz method in space and a fractional Adams–Moulton scheme in time were employed, observing increased damping for larger fractional orders. Under the same model, Faraji Oskouie et al. [37] incorporated surface stress effects through the Gurtin–Murdoch theory. Recently, Eyebe et al. [38] studied the nonlinear vibration of a nanobeam over a fractional Winkler–Pasternak foundation, utilizing D'Alembert principle to obtain a nonlinear system of equations, solved by a method of multiple scales. Lewandowski and Wielentejczyk [39] analyzed the nonlinear, steady-state vibration of a fractional Zener beam, obtaining amplitude equations with explicit finite element tangent matrices. The authors also studied the stability and parametric influence of their model.

An important application of interest is the dynamics of human's outer *hair cells* inside the fluid-filled cochlea of the inner ear. The corresponding structures are indeed “beam-shaped” whose dysfunction leads to sensorineural hearing loss. Direct measurements from the cochlea are not possible, which makes computational modeling a valuable (noninvasive) tool to study the health and function of the hair cells. The outer hair cells sit on the basilar membrane (i.e., the beam base) and are embedded into the tectorial membrane at the other end. The sound transmission in the cochlea leads to a pressure difference across the basilar membrane leading to the vibration of hair cells [40]. The cochlea has been modeled previously based on a fractional-modeling approach [16,17]. The proposed distributed-order fractional modeling in the current work can be readily employed in studying the fluid-induced vibrations of the hair cells, leading to new experiment setups and sensor developments.

The sophistication of numerical methods allowed numerous applications of fractional models in the last two decades, such as spectral methods for spatio-temporal discretization of FDEs [41,42] and DODEs [43]. Among many schemes for time-fractional integration [44–47], for simplicity, we outline the L1 finite difference scheme by Lin and Xu [48] and refer to Ref. [49] for a brief review of numerical methods for time-fractional ODEs. Despite the existing works on nonlinear vibration of fractional viscoelastic beams, they employed direct Ritz discretizations in the strong forms of the equations of motion, requiring more smoothness on basis functions. Spectral methods for nonlinear fractional beam models with proper finite dimensional function spaces suitable for fractional operators are still lacking. Furthermore, from the rheology standpoint, studying the emergence of anomalous dynamics from evolving extraordinary material properties and their sensitivity is fundamental for physics- and mathematically informed learning of constitutive laws from data, and also requires more attention.

In this work, we study how the evolving power-law rheology, in the language of fractional constitutive laws lead to (counterintuitive) anomalous dynamics in mechanical systems. Our representation of choice is a geometrically nonlinear, fractional KV Euler–Bernoulli cantilever beam under free and forced vibration, where:

The fractional KV model with order

*α*is obtained from both the Boltzmann superposition principle and general distributed-order viscoelasticity.Motivated by the effect of evolving fractal microstructures on macroscopic material dynamics, we study the effects of fractional orders on the continuum response.

We employ Hamilton's principle to avoid a non-trivial decomposition of conservative (elastic) and non-conservative (viscous) parts of SB elements.

The weak form of the governing equation is obtained, followed by a single-mode approximation in space.

We solve the resulting nonlinear system through a method of multiple scales, followed by a sensitivity analysis of amplitude decay rates with respect to

*α*.

Our numerical and semi-analytical experiments demonstrate several anomalous responses linked to far-from-equilibrium dynamics, such as *α*-dependent hardening-like drifts in linear amplitude–frequency behavior and long-term power law response. An interesting new result indicating a super sensitivity of amplitude response with respect to *α* was obtained, which could be potentially related to the identification of early damage precursors, before the onset of macroscopic plasticity and cracks [3]. Specifically, a softening-like behavior is observed until a critical *α*-value, followed by a hardening-like response, both justified from the constitutive standpoint. This motivates the notion of evolving anomalies, where the changing fractal material microstructure drives the fractional operator form through *α* [50]. Finally, we observe the usual bifurcation behavior under steady-state amplitude at primary resonance in the presence of geometrical nonlinearity, which is also driven by material nonlinearity through the fractional order.

This work is organized as follows: In Sec. 2, we introduce the model assumptions, fractional viscoelasticity definitions, the extended Hamilton's principle, and obtain the strong/weak forms of equation of motion under base excitation. We employ assumed modes in space to reduce the problem to a system of fractional ODEs. In Sec. 3, we obtain the linearized equation of motion. Perturbation analysis in carried out in Sec. 4 to solve the resulting nonlinear fractional ODE. Numerical results are presented, followed by conclusions in Sec. 5. We refer to the arXiv version of this work for more detailed formulation steps and proofs.

## 2 Mathematical Formulation

We formulate our anomalous physical system and discuss the main assumptions utilized in our model derivations.

### 2.1 Kinematics of the Anomalous Cantilever Beam.

Let the geometrically nonlinear, viscoelastic, Euler–Bernoulli cantilever beam with length *L*, subject to harmonic vertical base excitation *v _{b}* (Fig. 1) with the following assumptions: (

**1)**The stretch over the neutral axis, warping, and shear strains are negligible. (

**2)**The beam is slender with symmetric cross section, and undergoes pure bending. (

**3)**The length

*L*, cross section area

*A*, linear density

*ρ*, and rotational inertia

*J*of the lumped mass

*M*are constant. (

**4)**The axial and lateral displacements along the length are denoted, respectively, by

*u*(

*s*,

*t*) and

*v*(

*s*,

*t*). (

**5)**We consider the in-plane vertical vibration of the beam and reduce the problem to 1-dimension.

*x*,

*y*,

*z*) and moving $(x\u2032,y\u2032,z\u2032)$ coordinate systems, the latter attached to the base, i.e., $(x0\u2032,y0\u2032,z0\u2032)=(0,vb,0)$. A differential element rotates with angle $\psi (s,t)$ about $z\u2032$ to the coordinate system $(\xi ,\eta ,\zeta )$, as $[e\xi \u2009\u2009e\eta \u2009\u2009e\zeta ]T=Rz(\psi )\u2009[ex\u2032\u2009\u2009ey\u2032\u2009\u2009ez\u2032]T$, where $ei$ is the unit vector along the $i\u200ath$ coordinate and $Rz(\psi )$ is an axis rotation matrix. The angular velocity and curvature at any coordinate

*s*and time

*t*can be written, respectively, as $\omega (s,t)=(\u2202\psi /\u2202t)ez\u2032$ and $\rho (s,t)=(\u2202\psi /\u2202s)ez\u2032$. Hence, the total displacement and velocity of any point with respect to (

*x*,

*y*,

*z*) reads

Let an arbitrary neutral axis element *CD* with initial length *ds* and material coordinate *s* deform to configuration $C*D*$ (see Fig. 2). The displacement components of points *C* and *D* are denoted by (*u*, *v*) and $(u+du,v+dv)$, respectively. The axial strain *e*(*s*, *t*) at *C* is given by $e=(ds*\u2212ds)/ds=(1+\u2202u/\u2202s)2+(\u2202v/\u2202s)2\u22121$. Applying the inextensionality assumption *e* = 0, we have $1+\u2202u/\u2202s=(1\u2212(\u2202v/\u2202s)2)1/2$. Under the latter equation and negligible vertical shear strains, we have $\psi =tan\u22121[(\u2202v/\u2202s)(1\u2212(\u2202v/\u2202s)2)\u22121/2]$. Using expansion $\u2009tan\u22121(x)=x\u221213x3+\cdots $, the curvature can be approximated up to third-order terms as $\psi \u2243\u2202v/\u2202s+(1/6)(\u2202v/\u2202s)3$. Therefore, the angular velocity and curvature of the beam can be approximated as $\u2202\psi /\u2202t\u2243(\u22022v/\u2202t\u2202s)(1+(1/2)(\u2202v/\u2202s)2)$ and the curvature approximated as $\u2202\psi /\u2202s\u2243(\u22022v/\u2202s2)(1+(1/2)(\u2202v/\u2202s)2)$. Finally, under assumptions (**1)** and (**2)**, the strain–curvature relationship takes the form $\epsilon (s,t)=\u2212\eta \u2202\psi (s,t)/\u2202s$.

### 2.2 Linear Viscoelasticity: Boltzmann Superposition Principle.

*bottom–up*derivation of our SB building block model through the Boltzmann superposition principle. Then, in a

*top–bottom*approach, we demonstrate the recovery of the fractional KV model from a general distributed-order form. Under linear viscoelasticity, applying a small step strain increase, denoted by $\delta \epsilon (t)$, at $t=\tau 1$, the stress state is given by $\sigma (t)=G(t\u2212\tau 1)\delta \epsilon (\tau 1),\u2009\u2009t>\tau 1$, where

*G*(

*t*) denotes the relaxation function. The Boltzmann superposition principle states that the total stress at time

*t*is obtained from the linear superposition of all infinitesimal step strains applied at prior times

*τ*. Hence, $\sigma (t)=\u2211\tau j<tG(t\u2212\tau j)(\delta \epsilon (\tau j)/\delta \tau j)\delta \tau j$, where the limiting case $\delta \tau j\u21920$ yields

_{j}and $\epsilon \u0307$ denotes the strain rate.

#### 2.2.1 Exponential Relaxation (Classical Models) Versus Power-Law Relaxation (Fractional Models).

*G*(

*t*) is traditionally expressed as a linear combination exponential functions, which yields the

*so-called*generalized Maxwell form as $G(t)=\u2211Cie\u2212t/\tau i$. For the simple case of a single exponential term (a single Maxwell branch), we have $G(t)=Ee\u2212t/\tau $. Therefore, in the case of zero initial strain $(\epsilon (0)=0)$, we have

which solves the ODE $\u2202\epsilon /\u2202t=E\u22121\u2202\sigma /\u2202t+\eta \u22121\sigma $, where the relaxation time constant $\tau =\eta /E$ is experimentally obtained.

*pseudo-constant*. If we choose $g(\alpha )=1/\Gamma (1\u2212\alpha )$, the integro-differential operator (Eq. (2)) gives the Liouville-Weyl fractional derivative [51]. Although the lower integration limit of Eq. (2) is $\u2212\u221e$, assuming causality (the material is quiescent for

*t*<

*0), Eq. (2) can be rewritten as*

which are equivalent when $\epsilon (0)=0$.

#### 2.2.2 Multipower Laws Versus Distributed-Order Models.

*G*(

*t*) in Eq. (1) would contain a power-law distribution instead of a single power-law in Eq. (2). Considering nonlinear viscoelasticity with material heterogeneities, the distributed order constitutive equations over

*t*>

*0 with orders $\alpha \u2208[\alpha min,\alpha max]$ and $\beta \u2208[\beta min,\beta max]$ can be expressed in the general form as*

where $*$ denotes any definition of fractional derivative. The terms $\Phi (\beta ;x,t,\sigma )$ and $\Psi (\alpha ;x,t,\epsilon )$ denote distribution functions, where $\alpha \u21a6\Psi (\alpha ;x,t,\epsilon )$ and $\beta \u21a6\Phi (\beta ;x,t,\sigma )$ are continuous mappings in $[\alpha min,\alpha max]$ and $[\beta min,\beta max]$. Furthermore, the dependence of the distributions on the (thermodynamically) conjugate pair $(\sigma ,\epsilon )$ introduces the notion of nonlinear viscoelasticity, and the dependence on a material coordinate *x* induces material heterogeneties in space.

*Remark 2*. The pairs (

*α*

_{min},

*α*

_{max}) and (

*β*

_{min},

*β*

_{max}) are theoretical lower/upper bounds in definition (Eq. (6)). In general, $\Phi (\beta ;x,t,\sigma )$ and $\Psi (\alpha ;x,t,\epsilon )$ can arbitrarily confine the integration domain in each realization of practical rheological problems. If we let the distributions be linear combinations of Dirac-delta functions, Eq. (6) recovers a multiterm FDE

### 2.3 Extended Hamilton's Principle.

We derive the equations of motion by employing the extended Hamilton's principle $\u222bt1t2(\delta T\u2212\delta W)\u2009dt=0$, where *δT* and *δW* denote the variations of kinetic energy and total work [53]. The only external source of our system with volume $V$ is the base excitation, which is linearly superposed to the vertical displacement *v*(*t*), and therefore contributes to the kinetic energy in the inertial (Lagrangian) coordinate system. Hence, the total work only contains the internal work done by the stresses, with variation expressed as $\delta W=\u222bV\sigma \u2009\delta \epsilon \u2009dv$ [54].

*Remark 3*. Although employing Lagrange's equations would be more straightforward from the variational perspective in standard rheology, it would require an elastic/viscous stress separability, corresponding to conservative/nonconservative energies, which, as remarked in Eq. (1), is nontrivial in the time-domain. Existing studies have successfully provided such decomposition for SB elements at the energy (functional) level [55], opening interesting grounds to employ more complex fractional visco-elastic models (with fractional derivatives of stress) into variational frameworks. Nevertheless, due to the simplicity (regarding the relaxation representation) of the fractional KV model, we choose the extended Hamilton's principle, as it allows us to work with the total virtual work, avoiding further mathematical complexity.

*V*, the problem reads as: find $v\u2208V$ such that

### 2.4 Weak Formulation.

*E*, then we have $B\alpha \u2009(\Omega ):={v\u22080lH\alpha (\Omega )\u2009|\u222b\Omega E\u2009d\Omega <\u221e}$, where $0lH\alpha (\Omega )=0lH\alpha \u2009(I;L2\u2009(\Omega )\u2229L2(I;0H2(\Omega ))$, and $0H2(\Omega )={v\u2208H2(\Omega )\u2009|\u2009\u2009v\u2009|s=0=\u2202v/\u2202s\u2009|s=0=0}$. We obtain the weak form of the problem by multiplying Eq. (8) with proper test functions $v\u0303(s)\u2208B\alpha (\Omega )$ that satisfy the boundary conditions $v\u0303(0)=\u2202v\u0303/\u2202s(0)=0$, and integrate over the computational domain $\Omega s=[0,1]$. The result is integrated by parts and rearranged, leading to

### 2.5 Assumed Mode: Spectral Space Approximation.

for all $v\u0303N\u2208V\u0303N$.

### 2.6 Single-Mode Approximation.

In general, the modal discretization Eq. (11) in Eq. (12) leads to a coupled nonlinear system of fractional ODEs. We note that fractional operators already impose numerical challenges, which are exacerbated under nonlinearities, leading to failure of existing numerical schemes. However, without loss of generality, we can assume that only the primary mode of motion is involved in the dynamics of the system of interest.

#### 2.6.1 Why is Single-Mode Approximation Useful?.

Although single-mode approximations are simplistic in nature, they encapsulate the most fundamental dynamics and the highest energy mode in the motion of nonlinear systems, making them capable of capturing complex structural behaviors. Azrar et al. [57] showed sufficient approximations of single- and multimodal representations for the nonlinear vibration of a simply supported beam under a harmonic distributed force. Tseng and Dugundji [58] obtained similar results between single- and two-mode approximations for nonlinear vibrations of clamped–clamped beams far from the crossover region. Loutridis et al. [59] implemented a crack detection method for beams using a single DOF system with time-varying stiffness. In Ref. [60], the effects of base stiffness and attached mass on the nonlinear, flexural free vibrations of beams were studied. Lestari and Hanagud [61] studied the nonlinear free vibrations of buckled beams with elastic end constraints, where the single-mode assumption led to a closed-form solution in terms of elliptic functions.

*prior*to crack formation. Their findings demonstrate that the pragmatism of a single-mode approximation yields sufficient sensitivity of the amplitude response with respect to the nonlinear stiffness, making their framework an effective practical tool for early fatigue detection. We refer to Refs. [62–65] for additional applications. Therefore, we let the anomalous dynamics of our system be driven by the fractional-order

*α*, and following the aforementioned studies, we replace Eq. (11) with the single-mode discretization $vN=q(t)\u2009\varphi (s)$ (where we let

*N*=

*1 and drop the subscript 1 for simplicity). Upon substituting in Eq. (12), we obtain the unimodal equation of motion as (see Ref. [51], Appendix B),*

## 3 Linearized Equation: Numerical Time Integration.

*N*time-steps of size $\Delta t$, such that $tn=n\Delta t,\u2009n=0,\u20091,\u2009\u2026,\u2009N$. We employ the equivalency $0RLDt\alpha q(t)=0CDt\alpha q(t)+q(0)t\u2212\alpha /\Gamma (1\u2212\alpha )$ between the RL and Caputo definitions into Eq. (15) and evaluate the resulting from implicitly at $t=tn+1$. The fractional Caputo derivative is discretized through an L1-difference scheme [48], which yields

*α*-dependent convolution coefficients $bj=(j+1)\alpha \u2212j\alpha $ and $E*=(Er\u2009cl)/(\Delta t\alpha \Gamma (2\u2212\alpha ))$. We approximate the acceleration $q\xa8n+1$ and velocity $q\u0307n+1$ through a Newmark-

*β*method

Since the Newmark method is $O(\Delta t2)$-accurate, the overall accuracy is dominated by the L1 scheme, which is $O(\Delta t2\u2212\alpha )$. We also note that a discretization of a Caputo-variant of Eq. (15) is recovered if we remove the term $q0(1\u2212\alpha )/(n+1)\alpha $ from Eq. (19). Consider two numerical tests: The first under harmonic base excitation, and the second under a free-vibration response. For both, we set *E _{r}* = 1 and consider the lumped mass at the tip, with $M=J=1$, i.e., we utilize (D.7) from Ref. [51] for $\varphi (s)$, which yields the coefficients $cl=kl=1.24$. Our developed software was implemented in python and run in a 6-core, 2.6 GHz Intel machine with 32 GB RAM.

### 3.1 Harmonic Base Excitation.

We let the harmonic form $vb=ab\u2009sin(\omega bt)$ in Eq. (15), where $\omega b\u2208[0.5,3.5]$ and $ab=0.01$ denote, respectively, the base displacement frequency and amplitude. The value $mb=\u22120.042$ is obtained by Eqs. (16) and (14). We set $q(0)=0,\u2009q\u0307(0)=0$, and time $t\u2208(0,100]$, with step-size $\Delta t=10\u22123$, and evaluate the maximum displacement amplitude after attaining a steady-state response. Figure 3 illustrates the amplitude *versus* base frequency response under varying values of *α*. We observe the existence of a critical point at $\omega b=1$ that changes the dissipation nature of the fractional order. Regarding the maximum observed amplitudes, increasing the fractional order in the range $\alpha \u2208[0.1,0.4]$ decreases and slightly shifts the amplitude peaks to higher (right) frequencies (an anomalous quality). Conversely, as *α* increases in the range $\alpha \u2208[0.5,0.6]$, the peak amplitudes slightly shift toward the lower (left) frequencies, which is also observed in standard systems with the increase of modal damping. We note that although there are discussions in the literature about the nonexistence of steady-state solutions for fractional oscillators utilizing the Caputo definition under nonzero initial conditions [67], it is not the case in this work. Here, both definitions are equivalent in Eq. (15) due to the employment of homogeneous initial conditions.

### 3.2 Free Vibration.

Following the emergence of a critical point nearby the standard natural frequency in Fig. 3, we solve Eq. (15) under free-vibration employing RL and Caputo definitions, with $v\xa8b=0$, and $q(0)=0.01,\u2009q\u0307(0)=0$. Figure 4 (left) illustrates the obtained results for *q*(*t*) for varying *α* using a RL definition. We observe an *α*-dependent amplitude decay that recovers a classical oscillator as $\alpha \u21921$. In addition, an anomalous transient region is observed at the short time-scale $t\u2208[0.1,\u20091]$. Conversely, in Fig. 4 (right), anomalies are present at large time-scales through a (far-from-equilibrium) power-law relaxation, while the short-time behavior is “standard-like.” Such contrast provides interesting insights on modeling desired anomalous ranges in such power-law materials. When replacing the SB element with a dashpot (see Fig. 5), neither anomalous dynamics are present. These results are in agreement with the power-law and exponential relaxation kernels described in Sec. 2.2. Since the SB element yields a constitutive interpolation between springs and dashpots, it contributes to both effective stiffness and damping ratio of the system, hence increasing values of *α* (decreasing stiffness), yield a reduced frequency response. The observed anomalous amplitude decay when employing the Caputo definition is also in agreement with long time dynamics of fractional linear oscillators reported in Ref. [68]. Finally, anomalous decay rates were also investigated in Ref. [69] for an extended theory of decay of classical vibration models under nonlinear resonance, where sharper “nonexponential” decays rates in dynamical system variables were observed nearby resonance.

## 4 Perturbation Analysis of the Nonlinear Equation

We use perturbation analysis to solve our nonlinear FDE, which is reduced to an algebraic form to be solved for steady-state amplitude and vibration phase.

### 4.1 Method of Multiple Scales.

*ε*, which yields

*a*and $\phi $ are the amplitude and phase lag of time response, respectively. Thus, the solution

*q*(

*t*) becomes

where *a* and $\phi $ are obtained by real/imaginary parts separation.

#### 4.1.1 Case 1: No Lumped Mass at the Tip.

We let $M=J=0$, and thus, we employ (D.8) from Ref. for the eigenfunctions $\varphi (s)$, from which we compute $M=1,\u2009Kl=Cl=12.3624,\u2009Mb=0.782992$, and $Knl=Cnl=20.2203$. We consider the following cases:

Free vibration,

*F*= 0: Super Sensitivity to*α*The governing equations for amplitude and phase are given by$dadT1=\u2212Er\u2009\omega 0\alpha \u22121\u2009sin(\alpha \pi 2)(12\u2009cl\u2009a+38\u2009cnl\u2009a3)$(26)Note that the amplitude decays with rate $\tau d=cl\u2009Er\u2009\omega 0\alpha \u22121\u2009sin(\alpha \pi /2)$ in Eq. (26), which depends on the pair $(Er,\u2009\alpha )$ (Fig. 6). Let the sensitivity index $S\tau d,\alpha $, defined as$d\phi dT1=clEr2\omega 0\alpha \u22121\u2009cos(\pi \alpha 2)+3a24[cnlEr\omega 0\alpha \u22121\u2009cos(\pi \alpha 2)+knl\omega 0\u22121]$(27)which is computed and plotted in Fig. 7 for the same set of parameters. There exists a critical value $\alpha cr=\u2212(2/\pi )\u2009tan\u22121(\pi /(2\u2009log\u2009\omega 0))$, where $(dS\tau d,\alpha /d\alpha )=0$. We observe that by increasing$S\tau d,\alpha =d\tau dd\alpha =\pi 2cl\u2009Er\u2009\omega 0\alpha \u22121\u2009\u2009cos(\alpha \pi 2)+cl\u2009Er\u2009\omega 0\alpha \u22121\u2009sin(\alpha \pi 2)log(\omega 0)$*α*when $\alpha <\alpha cr$ (adding viscosity), the dissipation rate, and thus decay rate, increases; this can be interpreted as a softening (stiffness-decreasing) region. Further increasing*α*when $\alpha >\alpha cr$ conversely leads to decreased decay rates, which can be interpreted as a hardening (more stiffening) region. We also note that*α*_{cr}, dictating the super-sensitivity region where anomalous softening/hardening transitions take place, solely depending on the standard natural frequency*ω*_{0}given in Eq. (22). Although the observed hardening response when $\alpha >\alpha cr$ might at first seem counter-intuitive, we remark that here the notions of softening and hardening have mixed anomalous nature regarding energy dissipation and time-scale dependent material stress response. Similar anomalous dynamics were also observed in ballistic yield stress responses in fractional visco-elasto-plasticity [18]. In the following, we perform two numerical tests at the constitutive level of the fractional KV model Eq. (7) to justify the observed behavior in Fig. 6. We employ the*tangent loss*and the stress–strain response under monotone loads/relaxation.*Dissipation via tangent loss*: Taking the Fourier transform of Eq. (7), we obtain the dynamic modulus $G*$ [51], given byfrom which the real part yields the storage modulus $G\u2032(\omega )=E\u221e+E\alpha \omega \alpha \u2009cos(\alpha \pi /2)$ and the imaginary part yields the loss modulus $G\u2033(\omega )=E\alpha \omega \alpha \u2009sin(\alpha \pi /2)$, which represent, respectively, the stored and dissipated energies per loading cycle. Finally, we define the tangent loss $tan\u2009\delta loss=G\u2033(\omega )/G\u2032(\omega )$, which represents the ratio between the dissipated/stored energies, and hence related to the mechanical damping of the anomalous medium. We set $\omega =\omega 0$ and$G*(\omega )=E\u221e+E\alpha \omega \alpha (cos(\alpha \pi 2)+i\u2009sin(\alpha \pi 2))$(28)*E*= 1 and demonstrate the results for $tan\u2009\delta loss$ with varying fractional orders_{r}*α*. In Fig. 8 (top), increasing values of*α*led to increased values of the tangent loss, and the hardening part ($\alpha >\alpha cr$) is not associated with higher storage in the material. Instead, the increasing dissipation with*α*suggests an increasing damping of the mechanical structure.*Stress–time response for monotone loads/relaxation*: We demonstrate how increasing fractional orders for the model lead to increased hardening for sufficiently high strain rates. We discretize Eq. (7) through the L1-scheme [48] and set $E\u221e=1,\u2009E\alpha =1$. We also assume the following piecewise strain function: $\epsilon (t)=(1/24)t$, for $0\u2264t<2.5$ (monotone stress/strain), and $\epsilon (t)=1/10$ for $2.5\u2264t\u22646$ (stress relaxation). Figure 8 (bottom) illustrates that even for relatively low strain rates, there is a ballistic region nearby*t*= 0 where higher fractional orders yield higher stresses, characterizing a rate-dependent stress-hardening response. However, due to the dissipative nature of SB elements, the initially ballistic material softens after a critical time, due to its faster relaxation nature.Primary Resonance Case, $\Omega \u2248\omega 0$

Let $\Omega =\omega 0+\epsilon \u2009\Delta $, where Δ is called the detuning parameter and thus, we write the force function as $12F\u2009ei\u200a\Delta \u200aT1\u2009ei\u200a\omega 0\u200aT0+c.c$. In this case, the force function also contributes to the secular terms. Therefore, we find the governing equations of solution amplitude and phase as$dadT1=\u2212Er\omega 0\alpha \u22121\u2009sin(\pi \alpha 2)(cl2\u2009a+3cnl8\u2009a3)+12f\omega 0\u22121\u2009sin(\Delta T1\u2212\varphi )$(29)in which the parameters ${\alpha ,Er,f,\Delta}$ mainly change the frequency response of the system. Equations (29) and (30) can be transformed into an autonomous system by letting $\gamma =\Delta \u2009T1\u2212\varphi $. The steady-state solution occurs when $da/dT1=d\varphi /dT1=0$, which yields, after algebraic manipulations$a\u2009d\varphi dT1=12cl\u2009Er\u2009\omega 0\alpha \u22121\u2009\u2009cos(\pi \alpha 2)\u2009a+34cnl\u2009Er\u2009\omega 0\alpha \u22121\u2009\u2009cos(\pi \alpha 2)(x)\u2009a3+34\u2009\omega 0\u22121\u2009knl\u2009a3\u221212f\u2009\omega 0\u22121\u2009\u2009cos(\Delta T1\u2212\varphi )$(30)$[A1\u2009a+A2\u2009a3]2+[B1\u2009a+B2\u2009a3]2=CA1=cl2Er\u2009\omega 0\alpha \u22121\u2009sin(\pi \alpha 2),\u2003A2=3cnl8Er\u2009\omega 0\alpha \u22121\u2009sin(\pi \alpha 2)B1=\Delta \u2212cl2Er\u2009\omega 0\alpha \u22121\u2009\u2009cos(\pi \alpha 2)B2=\u221234(cnl\u2009Er\u2009\omega 0\alpha \u22121\u2009\u2009cos(\pi \alpha 2)+\omega 0\u22121\u2009knl),\u2003C=f24\u2009\omega 02$(31)

which is a cubic equation in *a*^{2}. The discriminant of a cubic equation of the form $ax3+bx2+cx+d=0$ is given as $\u03d1=18abcd\u22124b3d+b2c2\u22124ac3\u221227a2d2$. Equation (32) has one real root when $\u03d1<0$ and three distinct real roots when $\u03d1>0$. The parameters ${\alpha ,Er,f,\Delta}$ dictate the values of coefficients ${A1,A2,B1,B2,C}$, the value of *ϑ*, and thus the number of admissible steady-state amplitudes. We see that for fixed values of ${\alpha ,Er,f}$, by sweeping the detuning parameter Δ from lower to higher excitation frequency, the stable steady-state amplitude bifurcates into two stable branches and one unstable branch, where they converge back to a stable amplitude by further increasing Δ. Figure 9 (top) shows the bifurcation diagram by sweeping Δ under varying *α* when $Er=0.3$ and *f* = 1. The solid and dashed black lines are the stable and unstable amplitudes, respectively. Note that the bifurcation points (red) are strongly related to the value of *α*, where adding extra viscosity to the system leads to a faster recovery from unstable states. Figure 9 (bottom) shows the magnitude of steady-state amplitudes versus excitation frequency. As a forward sweep is imposed on the excitation frequency, the steady-state amplitude increases, reaches a peak value, and then jumps down to a stable path (red dashed line for $\alpha =0.4$). The peak amplitudes decrease as *α* increases. Figure 10 shows the frequency response of the system for different values of ${\alpha ,Er}$ when *f* = 0.5. The amplitude peak decreases with the increase of *E _{r}*. For both higher values of

*E*and

_{r}*α*, the amplitude peaks drift toward lower frequencies, displaying a softening behavior, while eliminating the bifurcation instability. The

*α*-dependent observed softening behavior is in qualitative agreement with other nonlinear systems reported in the literature [67], including Duffing oscillators [73] and damage-induced softening [3].

#### 4.1.2 Case 2: Lumped Mass at the Tip.

We let $M=J=1$, and employ eigenfunctions (D.7) from Ref. [51], which yields $M=1+70.769J+7.2734M,\u2009J=5008.25,\u2009Kl=Cl=98.1058,\u2009Mb=\u22120.648623\u22122.69692M$, and $Knl=Cnl=2979.66$. Similar to Case 1, we consider the following

Free vibration,

*F*= 0In this setting, the amplitude equation preserves its structure, while the phase equation has an extra term$dadT1=\u2212Er\u2009\omega 0\alpha \u22121\u2009sin(\alpha \pi 2)(12\u2009cl\u2009a+38\u2009cnl\u2009a3)$(33)$d\varphi dT1=12cl\u2009Er\u2009\omega 0\alpha \u22121\u2009\u2009cos(\pi \alpha 2)+34cnl\u2009Er\u2009\omega 0\alpha \u22121\u2009\u2009cos(\pi \alpha 2)\u2009a2+34\u2009\omega 0\u22121\u2009knl\u2009a2\u221214\u2009mnl\u2009\omega 0\u2009a2$(34)Primary resonance case, $\Omega \u2248\omega 0$

Similar to free vibration, we obtain the autonomous steady-state solution for amplitude and phasewhich can be written in the same form and coefficients as case 1, except for the coefficient $B2=\u2212(3/4)(cnl\u2009Er\u2009\omega 0\alpha \u22121\u2009\u2009cos(\pi \alpha /2)+\omega 0\u22121\u2009knl+(mnl/3)\omega 0)$, and yields the steady-state amplitude of the system.$[cl2Er\u2009\omega 0\alpha \u22121\u2009sin(\pi \alpha 2)\u2009a+3cnl8Er\u2009\omega 0\alpha \u22121\u2009sin(\pi \alpha 2)\u2009a3]2\xd7[(\Delta \u2212cl2Er\u2009\omega 0\alpha \u22121\u2009\u2009cos(\pi \alpha 2))a\u221234(cnl\u2009Er\u2009\omega 0\alpha \u22121\u2009\u2009cos(\pi \alpha 2)+\omega 0\u22121\u2009knl+13\u2009mnl\u2009\omega 0)a3]2=f24\u2009\omega 02$(35)

## 5 Summary and Discussion

We studied anomalous nonlinear dynamics driven by the application of extraordinary materials, modeled as a nonlinear fractional KV viscoelastic cantilever beam. A spectral method was employed for spatial discretization of the equation of motion, reducing it to a set of nonlinear fractional ODEs. The corresponding system was linearized and time-fractional integration was carried out through a L1-difference scheme and a Newmark method. A semi-analytical method of multiple scales was employed to the nonlinear system solution. We performed numerical experiments on the system response with varying fractional orders that represent different stages of material evolution, where we observed:

Anomalous,

*α*-dependent vibration amplitude drift, with a low-frequency critical point in linear forced vibration.Short- and long-time anomalous behaviors in linear free vibration, respectively, for RL and Caputo definitions.

Super sensitivity of the amplitude response with respect to the fractional model parameters at free vibration.

Critical decay rate sensitivity behavior with respect to

*α*, with higher/lower decay rate threshold*α*_{cr}separating softening/hardening regimes.A bifurcation behavior under steady-state amplitude at primary resonance case.

The choice of a fractional KV model represents materials in the intersection of anomalous/standard behaviors, where the SB element yields a power-law response, and the Hookean spring models the instantaneous response of many engineering materials. In addition, the shifts in amplitude–frequency response with respect to *α* motivate future studies on the downscaling of fractional operators to the far-from-equilibrium dynamics (polymer caging/reptation, dislocation avalanches) in evolving microstructures [50]. Regarding modifications of the current model, different material distribution functions can be chosen, leading to application-based material design for a wide range of complex materials/systems, including micro-electromechanical systems. Finally, regarding numerical discretizations, additional active vibration modes can be employed, as well as faster time-fractional numerical methods, to better approximate the fundamental dynamics of the presented system.

The developed model will be used in future to simulate the excitation of the outer hair cells inside the cochlea. The individual hair cells will be modeled as “cantilever beams” with variable sizes along the basilar membrane. In addition to the nonlocal and history effects in hair cell biomechanics, the *mechanoelectrical transduction* of the hair cells will be incorporated into the model.

## Funding Data

ARO Young Investigator Program Award (W911NF-19-1-0444; Funder ID: 10.13039/100000183).

National Science Foundation Award (DMS-1923201; Funder ID: 10.13039/100000183).

MURI/ARO (W911NF-15-1-0562; Funder ID: 10.13039/100000183).

AFOSR Young Investigator Program (Award FA9550-17-1-0150).

NIH NIDCD (Grant No. K01DC017751; Funder ID: 10.13039/100000183).

## References

**3**(1) pp. 61–90.10.1007/s42967-020-00070-w