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 $〈Δr〉2∝tα$ [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<α<1$, that interpolates between Hookean $(α→0)$ and Newtonian $(α→1)$ 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 [2022]. 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. [2527] 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 [4447], 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 vb (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.

Fig. 1
Fig. 1
Close modal
Consider the inertial (x, y, z) and moving $(x′,y′,z′)$ coordinate systems, the latter attached to the base, i.e., $(x0′,y0′,z0′)=(0,vb,0)$. A differential element rotates with angle $ψ(s,t)$ about $z′$ to the coordinate system $(ξ,η,ζ)$, as $[eξ eη eζ]T=Rz(ψ) [ex′ ey′ ez′]T$, where $ei$ is the unit vector along the $i th$ coordinate and $Rz(ψ)$ is an axis rotation matrix. The angular velocity and curvature at any coordinate s and time t can be written, respectively, as $ω(s,t)=(∂ψ/∂t)ez′$ and $ρ(s,t)=(∂ψ/∂s)ez′$. Hence, the total displacement and velocity of any point with respect to (x, y, z) reads
$r=(u−η sin(ψ)) ex+(v+vb+η cos(ψ)) ey∂r∂t=(∂u∂t−η ∂ψ∂t cos(ψ)) ex+(∂v∂t+∂vb∂t−η ∂ψ∂t sin(ψ)) ey$

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*−ds)/ds=(1+∂u/∂s)2+(∂v/∂s)2−1$. Applying the inextensionality assumption e = 0, we have $1+∂u/∂s=(1−(∂v/∂s)2)1/2$. Under the latter equation and negligible vertical shear strains, we have $ψ=tan−1[(∂v/∂s)(1−(∂v/∂s)2)−1/2]$. Using expansion $tan−1(x)=x−13x3+⋯$, the curvature can be approximated up to third-order terms as $ψ≃∂v/∂s+(1/6)(∂v/∂s)3$. Therefore, the angular velocity and curvature of the beam can be approximated as $∂ψ/∂t≃(∂2v/∂t∂s)(1+(1/2)(∂v/∂s)2)$ and the curvature approximated as $∂ψ/∂s≃(∂2v/∂s2)(1+(1/2)(∂v/∂s)2)$. Finally, under assumptions (1) and (2), the strain–curvature relationship takes the form $ε(s,t)=−η∂ψ(s,t)/∂s$.

Fig. 2
Fig. 2
Close modal

2.2 Linear Viscoelasticity: Boltzmann Superposition Principle.

We start with a 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 $δε(t)$, at $t=τ1$, the stress state is given by $σ(t)=G(t−τ1)δε(τ1), t>τ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 τj. Hence, $σ(t)=∑τj, where the limiting case $δτj→0$ yields
$σ(t)=∫−∞tG(t−τ) ε̇(τ) dτ$
(1)

and $ε̇$ denotes the strain rate.

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

The relaxation function G(t) is traditionally expressed as a linear combination exponential functions, which yields the so-called generalized Maxwell form as $G(t)=∑Cie−t/τi$. For the simple case of a single exponential term (a single Maxwell branch), we have $G(t)=Ee−t/τ$. Therefore, in the case of zero initial strain $(ε(0)=0)$, we have
$σ(t)=E∫0te−(t−t̃)/τ ε̇(t̃) dt̃$

which solves the ODE $∂ε/∂t=E−1∂σ/∂t+η−1σ$, where the relaxation time constant $τ=η/E$ is experimentally obtained.

By letting the relaxation function (kernel) in Eq. (1) to have a power-law form $G(t)=Eα g(α)( t−τ )−α$, we obtain
$σ(t)=Eα g(α) ∫−∞tε̇(τ)( t−τ )α dτ$
(2)
where $Eα [Pa⋅sα]$ denotes a pseudo-constant. If we choose $g(α)=1/Γ(1−α)$, the integro-differential operator (Eq. (2)) gives the Liouville-Weyl fractional derivative [51]. Although the lower integration limit of Eq. (2) is $−∞$, assuming causality (the material is quiescent for t <0), Eq. (2) can be rewritten as
$σ(t)=ε(0+) EαΓ(1−α)tα+Eα 0CDtα ε=Eα 0RLDtα ε$
(3)
where $0CDtα$ and $0RLDtα$ denote, respectively, the Caputo and Riemann–Liouville (RL) fractional derivatives [51]
$0CDtαε(t)=1Γ(1−α)∫0tε̇(s)(t−s)α ds, t>0$
(4)
$0RLDtαε(t)=1Γ(1−α)ddt∫0tε(s)(t−s)α ds, t>0$
(5)

which are equivalent when $ε(0)=0$.

Remark 1. Equation (3) represents the SB element [52], which interpolates between the limit cases of a Hookean spring $(α→0)$ and a Newtonian dashpot $(α→1)$.

2.2.2 Multipower Laws Versus Distributed-Order Models.

Generally, materials intrinsically possess a spectrum of power-law relaxations, requiring a distributed-order stress–strain representation. Therefore, the kernel 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 $α∈[αmin,αmax]$ and $β∈[βmin,βmax]$ can be expressed in the general form as
$∫βminβmaxΦ(β;x,t,σ)0*Dtβσ(t) dβ=∫αminαmaxΨ(α;x,t,ε)0*Dtαε(t)dα$
(6)

where $*$ denotes any definition of fractional derivative. The terms $Φ(β;x,t,σ)$ and $Ψ(α;x,t,ε)$ denote distribution functions, where $α↦Ψ(α;x,t,ε)$ and $β↦Φ(β;x,t,σ)$ are continuous mappings in $[αmin,αmax]$ and $[βmin,βmax]$. Furthermore, the dependence of the distributions on the (thermodynamically) conjugate pair $(σ,ε)$ 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, $Φ(β;x,t,σ)$ and $Ψ(α;x,t,ε)$ 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
$(1+∑k=1pσ ak0Dtβk) σ(t)=(c+∑k=1pε bk0Dtαk) ε(t)$
Finally, to obtain the fractional KV model, let $Φ(β)=δ(β)$ and $Ψ(α)=E∞δ(α)+Eαδ(α−α0)$ in Eq. (6). Hence
$σ(t)=E∞ ε(t)+Eα 0RLDtα ε(t), α∈(0,1)$
(7)

2.3 Extended Hamilton's Principle.

We derive the equations of motion by employing the extended Hamilton's principle $∫t1t2(δT−δW) dt=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 $δW=∫Vσ δε dv$ [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.

We refer to Suzuki et al. (Appendix A) for the full derivation of the governing equation presented here. Let $I=∫Aη2 dA, m=ρ/(E∞ I)$ and $Er=Eα/E∞$. We approximate the nonlinear terms up to third order and use the following dimensionless variables: $s*=s/L, v*=v/L, t*=t(mL4)−1/2, Er*=Er(mL4)−α/2, J*=J/(ρL3), M*=M/(ρL), vb*=vb/L$, and derive the strong form of the equation of motion. Therefore, by choosing a proper function space V, the problem reads as: find $v∈V$ such that
$∂2v∂t2+∂2∂s2(∂2v∂s2+∂2v∂s2(∂v∂s)2+Er0RLDtα [∂2v∂s2(1+12(∂v∂s)2)]+12Er(∂v∂s)2 0RLDtα ∂2v∂s2)−∂∂s(∂v∂s(∂2v∂s2)2+Er∂v∂s ∂2v∂s2 0RLDtα ∂2v∂s2)=−vb̈$
(8)
which is subject to the following boundary conditions
$v |s=0=∂v∂s |s=0=0J(∂3v∂t2∂s(1+(∂v∂s)2)+∂v∂s(∂2v∂s∂t)2)+(∂2v∂s2+∂2v∂s2(∂v∂s)2+Er 0RLDtα[∂2v∂s2(1+12(∂v∂s)2)] +12Er (∂v∂s)2 0RLDtα ∂2v∂s2) |s=1=0M(∂2v∂t2+vb̈)−∂v∂s(∂2v∂s2+∂2v∂s2(∂v∂s)2+Er0RLDtα[∂2v∂s2(1+12(∂v∂s)2)]+Er2(∂v∂s)20RLDtα∂2v∂s2)+(∂v∂s(∂2v∂s2)2+Er ∂v∂s∂2v∂s20RLDtα∂2v∂s2)|s=1=0$
(9)

2.4 Weak Formulation.

Common practices in numerical PDEs are mostly well developed for linear equations, but they are still scarce for nonlinear PDEs. Linear theories are usually applied to nonlinear problems assuming sufficiently smooth solutions [56]. We do not intend to develop an analysis for our proposed nonlinear model. Instead, we assume a smooth solution, and employ linear theories. Let $v:ℝ1+1→ℝ$ for $α∈(0,1)$ and $Ω=[0,T]×[0,L]$. We construct a complete solution space $Bα (Ω)$ [42], in order to formulate weak form of Eq. (8). Recalling Eq. (8) as E, then we have $Bα (Ω):={v∈0lHα(Ω) |∫ΩE dΩ<∞}$, where $0lHα(Ω)=0lHα (I;L2 (Ω)∩L2(I;0H2(Ω))$, and $0H2(Ω)={v∈H2(Ω) | v |s=0=∂v/∂s |s=0=0}$. We obtain the weak form of the problem by multiplying Eq. (8) with proper test functions $ṽ(s)∈Bα(Ω)$ that satisfy the boundary conditions $ṽ(0)=∂ṽ/∂s(0)=0$, and integrate over the computational domain $Ωs=[0,1]$. The result is integrated by parts and rearranged, leading to
$∂2∂t2 (∫01v ṽ ds+M v ṽ |s=1+J ∂v∂s ∂ṽ∂s |s=1)+J(∂3v∂t2∂s(∂v∂s)2+∂v∂s(∂2v∂t∂s)2)∂ṽ∂s |s=1+∫01∂2v∂s2 ∂2ṽ∂s2 ds+Er∫010RLDtα [∂2v∂s2] ∂2ṽ∂s2 ds+∫01∂2v∂s2(∂v∂s)2 ∂2ṽ∂s2 ds+∫01∂v∂s(∂2v∂s2) 2 ∂ṽ∂s ds+Er2∫010RLDtα[∂2v∂s2(∂v∂s)2] ∂2ṽ∂s2 ds+Er2∫01(∂v∂s)20RLDtα[∂2v∂s2]∂2ṽ∂s2ds+Er∫01∂v∂s∂2v∂s20RLDtα[∂2v∂s2]∂ṽ∂sds=−vb̈∫01ṽds+Mṽ|s=1$
(10)

2.5 Assumed Mode: Spectral Space Approximation.

We employ the following modal discretization to obtain a reduced-order model of the beam
$v(s,t)≃vN(s,t)=∑n=1Nqn(t) ϕn(s)$
(11)
where the spatial functions $ϕn(s), n=1,2,…,N$ (see Ref. [51], Appendix D) are obtained by solving the corresponding linear eigenvalue problem, and the temporal functions $qn(t), n=1,2,…,N$ are unknown modal coordinates. The assumed modes $ϕn(s)$ in Eq. (11). We then construct the proper finite dimensional spaces of basis/test functions as $VN=ṼN=span{ϕn(x):n=1,2,…,N}$. Since $VN=ṼN⊂V=Ṽ$, problem (Eq. (10)) reads: find $vN∈VN$ such that
$∂2∂t2(∫01vNṽNds+MvNṽN|s=1+J∂vN∂s∂ṽN∂s |s=1)+J(∂3vN∂t2∂s(∂vN∂s)2+∂vN∂s(∂2vN∂t∂s)2)∂ṽN∂s|s=1+∫01∂2vN∂s2 ∂2ṽN∂s2 ds+Er∫010RLDtα [∂2vN∂s2] ∂2ṽN∂s2 ds+∫01∂2vN∂s2(∂vN∂s)2 ∂2ṽN∂s2 ds+∫01∂vN∂s(∂2vN∂s)2 ∂ṽN∂s ds+Er2∫010RLDtα[∂2vN∂s2(∂vN∂s)2] ∂2ṽN∂s2 ds+Er2∫01(∂vN∂s)2 0RLDtα[∂2vN∂s2] ∂2ṽN∂s2 ds+Er∫01∂vN∂s∂2vN∂s20RLDtα[∂2vN∂s2]∂ṽN∂sds=−vḃ̇∫01ṽNds+MṽN|s=1s$
(12)

for all $ṽN∈ṼN$.

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.

Habtour et al. [3] validated the softening response of a nonlinear cantilever beam from local stress-induced, early fatigue damage 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. [6265] 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) ϕ(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),
$Mq¨+J(q¨q2+qq̇2)+Klq+ErCl0RLDtαq+2Knl q3+ErCnl2(0RLDtαq3+3q20RLDtαq)=−Mbvḃ̇$
(13)
in which
$M=∫01ϕ2 ds+M ϕ2(1)+J ϕ′2(1), Mb=∫01ϕ ds+M ϕ(1)Kl=Cl=∫01ϕ″2 ds, Knl=Cnl=∫01ϕ′2 ϕ″2 ds, J=J ϕ′4(1)$
(14)

Remark 4. Any mode of vibration $ϕn(s), n=1,2,…,N$ can be isolated by assuming $ϕn(s)$ is the only active one, and thus yield a similar equation of motion as Eq. (13), where the coefficients in Eq. (14) are obtained based on the active mode $ϕn(s)$.

3 Linearized Equation: Numerical Time Integration.

We linearize the equation of motion (Eq. (13)) by assuming small motions (see Ref. [51], Appendix C), and obtain
$q¨+Er cl 0RLDtαq+kl q=−mbv¨b$
(15)
with initial conditions $q(0)=v(L,0)/ϕ(L), q̇(0)=v̇(L,0)/ϕ(L)$ and the coefficients
$cl=ClM, kl=KlM, mb=MbM$
(16)
The linearized, unimodal form Eq. (15) is equivalent to a fractional KV oscillator. Let a uniform time-grid with N time-steps of size $Δt$, such that $tn=nΔt, n=0, 1, …, N$. We employ the equivalency $0RLDtαq(t)=0CDtαq(t)+q(0)t−α/Γ(1−α)$ 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
$q¨n+1+E*(qn+1−qn+Hαqn+1)+Er cl q0Γ(1−α)tn+1α+klqn+1=−mbv¨b,n+1$
(17)
where $Hk+1α=∑j=0n−1bj(qn−j+1−qn−j)$ represents the discretized history term, with α-dependent convolution coefficients $bj=(j+1)α−jα$ and $E*=(Er cl)/(ΔtαΓ(2−α))$. We approximate the acceleration $q¨n+1$ and velocity $q̇n+1$ through a Newmark-β method
$q¨n+1=a1(qn+1−qn)−a2q̇n−a3q¨nq̇n+1=a4(qn+1−qn)+a5q̇n+a6q¨n$
(18)
with approximation coefficients given in Ref. [66], which are chosen as $β=0.5, γ=0.25$ for unconditional stability. Inserting Eq. (18) into Eq. (17), we obtain the following closed form for $qn+1$:
$qn+1=(a1+E*)qn+a2q̇n+a3q¨n−mbv¨b,n+1−E*[Hn+1α+q0(1−α)(n+1)α]a1+E*+kl$
(19)

Since the Newmark method is $O(Δt2)$-accurate, the overall accuracy is dominated by the L1 scheme, which is $O(Δt2−α)$. We also note that a discretization of a Caputo-variant of Eq. (15) is recovered if we remove the term $q0(1−α)/(n+1)α$ 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 Er = 1 and consider the lumped mass at the tip, with $M=J=1$, i.e., we utilize (D.7) from Ref. [51] for $ϕ(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 sin(ωbt)$ in Eq. (15), where $ωb∈[0.5,3.5]$ and $ab=0.01$ denote, respectively, the base displacement frequency and amplitude. The value $mb=−0.042$ is obtained by Eqs. (16) and (14). We set $q(0)=0, q̇(0)=0$, and time $t∈(0,100]$, with step-size $Δt=10−3$, 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 $ωb=1$ that changes the dissipation nature of the fractional order. Regarding the maximum observed amplitudes, increasing the fractional order in the range $α∈[0.1,0.4]$ decreases and slightly shifts the amplitude peaks to higher (right) frequencies (an anomalous quality). Conversely, as α increases in the range $α∈[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.

Fig. 3
Fig. 3
Close modal

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¨b=0$, and $q(0)=0.01, q̇(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 $α→1$. In addition, an anomalous transient region is observed at the short time-scale $t∈[0.1, 1]$. 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.

Fig. 4
Fig. 4
Close modal
Fig. 5
Fig. 5
Close modal

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.

We investigate the dynamics of system (Eq. (13)) through the method of multiple scales [70,71]. Let the new independent time scales and their derivative operators be defined as $Tm=εm t, Dm=∂/∂Tm, m=0,1,2,…$. It is convenient to utilize another representation of the fractional RL derivative (see Ref. [72], Eq. (5.82)), as the fractional power of conventional time-derivatives, i.e., $0RLDtα=(d/dt)α$. Therefore
$ddt=D0+εD1+⋯, d2dt2=D02+2εD0D1+⋯0RLDtα=(ddt)α=D0α+εαD0α−1D1+⋯$
(20)
The response $q(T0,T1,…)$ can be written as the expansion
$q=q0(T0,T1,…)+εq1(T0,T1,…)+ε2q2(T0,T1,…)+⋯$
(21)
The equation of motion coefficients assume the scalings
$JM=ε mnl, KlM=kl=ω02, ClM=ε cl, KnlM=ε knl, CnlM=ε cnl$
(22)
also consider the base excitation $−(Mb/M)vb̈$ to be a harmonic function in the form $ε F cos(Ω t)$. Thus, employing Eqs. (20)(22), Eq. (13) can be expanded and we collect similar coefficients of zeroth and first orders of ε, which yields
$O(ε0):D02q0+ω02q0=0$
(23)
$O(ε1):D02q1+ω02q1=−2D0D1q0−mnl(q02D02q0+q0(D0q0)2)−ErclD0αq0−2knlq03−Ercnl2D0αq03−3Ercnl2q02D0αq0+F cos(ΩT0)$
(24)
The solution to Eq. (23) is of the form $q0(T0,T1)=A(T1) ei ω0 T0+c.c$, where “c.c” denotes the complex conjugate. Upon substituting this form into Eq. (24), distinct resonance cases are possible, for which we obtain solvability conditions by removing secular terms. Then, we employ the polar form $A=(a/2) ei φ$, where the real valued functions a and $φ$ are the amplitude and phase lag of time response, respectively. Thus, the solution q(t) becomes
$q(t)=a(ε t) cos(ω0 t+φ(ε t))+O(ε)$
(25)

where a and $φ$ 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 $ϕ(s)$, from which we compute $M=1, Kl=Cl=12.3624, Mb=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=−Er ω0α−1 sin(απ2)(12 cl a+38 cnl a3)$
(26)
$dφdT1=clEr2ω0α−1 cos(πα2)+3a24[cnlErω0α−1 cos(πα2)+knlω0−1]$
(27)
Note that the amplitude decays with rate $τd=cl Er ω0α−1 sin(απ/2)$ in Eq. (26), which depends on the pair $(Er, α)$ (Fig. 6). Let the sensitivity index $Sτd,α$, defined as
$Sτd,α=dτddα=π2cl Er ω0α−1 cos(απ2)+cl Er ω0α−1 sin(απ2)log(ω0)$
which is computed and plotted in Fig. 7 for the same set of parameters. There exists a critical value $αcr=−(2/π) tan−1(π/(2 log ω0))$, where $(dSτd,α/dα)=0$. We observe that by increasing α when $α<αcr$ (adding viscosity), the dissipation rate, and thus decay rate, increases; this can be interpreted as a softening (stiffness-decreasing) region. Further increasing α when $α>α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 $α>α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 by
$G*(ω)=E∞+Eαωα(cos(απ2)+i sin(απ2))$
(28)
from which the real part yields the storage modulus $G′(ω)=E∞+Eαωα cos(απ/2)$ and the imaginary part yields the loss modulus $G″(ω)=Eαωα sin(απ/2)$, which represent, respectively, the stored and dissipated energies per loading cycle. Finally, we define the tangent loss $tan δloss=G″(ω)/G′(ω)$, which represents the ratio between the dissipated/stored energies, and hence related to the mechanical damping of the anomalous medium. We set $ω=ω0$ and Er = 1 and demonstrate the results for $tan δloss$ with varying fractional orders α. In Fig. 8 (top), increasing values of α led to increased values of the tangent loss, and the hardening part ($α>α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∞=1, Eα=1$. We also assume the following piecewise strain function: $ε(t)=(1/24)t$, for $0≤t<2.5$ (monotone stress/strain), and $ε(t)=1/10$ for $2.5≤t≤6$ (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, $Ω≈ω0$

Let $Ω=ω0+ε Δ$, where Δ is called the detuning parameter and thus, we write the force function as $12F ei Δ T1 ei ω0 T0+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=−Erω0α−1 sin(πα2)(cl2 a+3cnl8 a3)+12fω0−1 sin(ΔT1−ϕ)$
(29)
$a dϕdT1=12cl Er ω0α−1 cos(πα2) a+34cnl Er ω0α−1 cos(πα2)(x) a3+34 ω0−1 knl a3−12f ω0−1 cos(ΔT1−ϕ)$
(30)
in which the parameters ${α,Er,f,Δ}$ mainly change the frequency response of the system. Equations (29) and (30) can be transformed into an autonomous system by letting $γ=Δ T1−ϕ$. The steady-state solution occurs when $da/dT1=dϕ/dT1=0$, which yields, after algebraic manipulations
$[A1 a+A2 a3]2+[B1 a+B2 a3]2=CA1=cl2Er ω0α−1 sin(πα2), A2=3cnl8Er ω0α−1 sin(πα2)B1=Δ−cl2Er ω0α−1 cos(πα2)B2=−34(cnl Er ω0α−1 cos(πα2)+ω0−1 knl), C=f24 ω02$
(31)
Fig. 6
Fig. 6
Close modal
Fig. 7
Fig. 7
Close modal
Fig. 8
Fig. 8
Close modal
$(A22+B22)a6+(2A1A2+2B1B2)a4+(A12+B12)a2−C=0$
(32)

which is a cubic equation in a2. The discriminant of a cubic equation of the form $ax3+bx2+cx+d=0$ is given as $ϑ=18abcd−4b3d+b2c2−4ac3−27a2d2$. Equation (32) has one real root when $ϑ<0$ and three distinct real roots when $ϑ>0$. The parameters ${α,Er,f,Δ}$ 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 ${α,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 $α=0.4$). The peak amplitudes decrease as α increases. Figure 10 shows the frequency response of the system for different values of ${α,Er}$ when f = 0.5. The amplitude peak decreases with the increase of Er. For both higher values of Er and α, 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].

Fig. 9
Fig. 9
Close modal
Fig. 10
Fig. 10
Close modal

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, J=5008.25, Kl=Cl=98.1058, Mb=−0.648623−2.69692M$, and $Knl=Cnl=2979.66$. Similar to Case 1, we consider the following

• Free vibration, F = 0

In this setting, the amplitude equation preserves its structure, while the phase equation has an extra term
$dadT1=−Er ω0α−1 sin(απ2)(12 cl a+38 cnl a3)$
(33)
$dϕdT1=12cl Er ω0α−1 cos(πα2)+34cnl Er ω0α−1 cos(πα2) a2+34 ω0−1 knl a2−14 mnl ω0 a2$
(34)
• Primary resonance case, $Ω≈ω0$

Similar to free vibration, we obtain the autonomous steady-state solution for amplitude and phase
$[cl2Er ω0α−1 sin(πα2) a+3cnl8Er ω0α−1 sin(πα2) a3]2×[(Δ−cl2Er ω0α−1 cos(πα2))a−34(cnl Er ω0α−1 cos(πα2)+ω0−1 knl+13 mnl ω0)a3]2=f24 ω02$
(35)
which can be written in the same form and coefficients as case 1, except for the coefficient $B2=−(3/4)(cnl Er ω0α−1 cos(πα/2)+ω0−1 knl+(mnl/3)ω0)$, and yields the steady-state amplitude of the system.

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

1.
Sagaut
,
P.
, and
Cambon
,
C.
,
2018
,
Homogeneous Turbulence Dynamics
, Springer International Publishing, Berlin.
2.
Akhavan-Safaei
,
A.
,
Seyedi
,
S.
, and
Zayernouri
,
M.
,
2020
, “
Anomalous Features in Internal Cylinder Flow Instabilities Subject to Uncertain Rotational Effects
,”
Phys. Fluids
,
32
(
9
), p.
094107
.10.1063/5.0021815
3.
Habtour
,
E.
,
Cole
,
D. P.
,
Riddick
,
J. C.
,
Weiss
,
V.
,
Robeson
,
M.
,
Sridharan
,
R.
, and
Dasgupta
,
A.
,
2016
, “
Detection of Fatigue Damage Precursor Using a Nonlinear Vibration Approach
,”
Struct. Control Health Monit.
,
23
(
12
), pp.
1442
1463
.10.1002/stc.1844
4.
Kapnistos
,
M.
,
Lang
,
M.
,
Vlassopoulos
,
D.
,
Pyckhout-Hintzen
,
D.
,
Richter
,
D.
,
Cho
,
D.
,
Chang
,
T.
, and
Rubinstein
,
M.
,
2008
, “
Unexpected Power-Law Stress Relaxation of Entangled Ring Polymers
,”
Nat. Mater.
,
7
(
12
), pp.
997
1002
.10.1038/nmat2292
5.
McKinley
,
G.
, and
Jaishankar
,
A.
,
2013
, “
Critical Gels, Scott Blair and the Fractional Calculus of Soft Squishy Materials
,”
Presentation
.
6.
Wong
,
I.
,
Gardel
,
M.
,
Reichman
,
D.
,
Weeks
,
E.
,
Valentine
,
M.
,
Bausch
,
A.
, and
Weitz
,
D.
,
2004
, “
Anomalous Diffusion Probes Microstructure Dynamics of Entangled F-Actin Networks
,”
Phys. Rev. Lett.
92, p.
178101
.10.1103/PhysRevLett.92.178101
7.
Bonakdar
,
N.
,
Gerum
,
R.
,
Kuhn
,
M.
,
Spörrer
,
M.
,
Lippert
,
A.
,
Schneider
,
W.
,
Aifantis
,
K. E.
, and
Fabry
,
B.
,
2016
, “
Mechanical Plasticity of Cells
,”
Nat. Mater.
,
15
(
10
), pp.
1090
1094
.10.1038/nmat4689
8.
Richeton
,
T.
,
Weiss
,
J.
, and
Louchet
,
F.
,
2005
, “
Breakdown of Avalanche Critical Behaviour in Polycrystalline Plasticity
,”
Nat. Mater.
,
4
(
6
), pp.
465
469
.10.1038/nmat1393
9.
Christensen
,
R.
,
2013
,
Theory of Viscoelasticity: An Introduction
,
Dover Publications, Mineola, NY
.
10.
Aster
,
R.
,
Borchers
,
B.
, and
Thurber
,
C.
,
2018
,
Parameter Estimation and Inverse Problems
, 3rd ed.,
Elsevier, Amsterdam, The
Netherlands.
11.
Bagley
,
R.
,
1989
, “
Power Law and Fractional Calculus Model of Viscoelasticity
,”
AIAA J.
,
27
(
10
), pp.
1412
1417
.10.2514/3.10279
12.
Jaishankar
,
A.
, and
McKinley
,
G.
,
2013
, “
Power-Law Rheology in the Bulk and at the Interface: Quasi-Properties and Fractional Constitutive Equations
,”
Proc. R Soc. A
,
469
(
2149
), p.
20120284
.10.1098/rspa.2012.0284
13.
Nutting
,
P.
,
1921
, “
A New General Law of Deformation
,”
J. Franklin Inst.
,
191
(
5
), pp.
679
685
.10.1016/S0016-0032(21)90171-6
14.
Gemant
,
A.
,
1936
, “
A Method of Analyzing Experimental Results Obtained From Elasto-Viscous Bodies
,”
Physics
,
7
(
8
), pp.
311
317
.10.1063/1.1745400
15.
Bagley
,
R.
, and
Torvik
,
P.
,
1983
, “
A Theoretical Basis for the Application of Fractional Calculus to Viscoelasticity
,”
J. Rheol.
,
27
(
3
), pp.
201
210
.10.1122/1.549724
16.
Naghibolhosseini
,
M.
,
2015
, “
Estimation of Outer-Middle Ear Transmission Using DPOAEs and Fractional-Order Modeling of Human Middle Ear
,” Ph.D. thesis,
City University of New York
,
New York
.
17.
Naghibolhosseini
,
M.
, and
Long
,
G.
,
2018
, “
Fractional-Order Modelling and Simulation of Human Ear
,”
Int. J. Comput. Math.
,
95
(
6–7
), pp.
1257
1273
.10.1080/00207160.2017.1404038
18.
Suzuki
,
J.
,
Zayernouri
,
M.
,
Bittencourt
,
M.
, and
,
G.
,
2016
, “
Fractional-Order Uniaxial Visco-Elasto-Plastic Models for Structural Analysis
,”
Comput. Methods Appl. Mech. Eng.
,
308
, pp.
443
467
.10.1016/j.cma.2016.05.030
19.
Suzuki
,
J.
,
Zhou
,
Y.
,
D'Elia
,
M.
, and
Zayernouri
,
M.
,
2021
, “
A Thermodynamically Consistent Fractional Visco-Elasto-Plastic Model With Memory-Dependent Damage for Anomalous Materials
,”
Comput. Methods Appl. Mech. Eng.
,
373
, p.
113494
.10.1016/j.cma.2020.113494
20.
Shitikova
,
M.
,
Rossikhin
,
Y.
, and
Kandu
,
V.
,
2017
, “
Interaction of Internal and External Resonances During Force Driven Vibrations of a Nonlinear Thin Plate Embedded Into a Fractional Derivative Medium
,”
Procedia Eng.
,
199
, pp.
832
837
.10.1016/j.proeng.2017.09.008
21.
Rossikhin
,
Y.
, and
Shitikova
,
M.
,
1997
, “
Applications of Fractional Calculus to Dynamic Problems of Linear and Nonlinear Hereditary Mechanics of Solids
,”
ASME Appl. Mech. Rev.
,
50
(
1
), pp.
15
67
.10.1115/1.3101682
22.
Samiee
,
M.
,
Akhavan-Safaei
,
A.
, and
Zayernouri
,
M.
,
2020
, “
A Fractional Subgrid-Scale Model for Turbulent Flows: Theoretical Formulation and a Priori Study
,”
Phys. Fluids
,
32
(
5
), p.
055102
.10.1063/1.5128379
23.
Lorenzo
,
C.
, and
Hartley
,
T.
,
2002
, “
Variable Order and Distributed Order Fractional Operators
,”
Nonlinear Dyn.
,
29
(
1/4
), pp.
57
98
.10.1023/A:1016586905654
24.
Atanackovic
,
T.
,
Oparnica
,
L.
, and
Pilipović
,
S.
,
2009
, “
Distributional Framework for Solving Fractional Differential Equations
,”
Integral Transform. Spec. Funct.
,
20
(
3–4
), pp.
215
222
.10.1080/10652460802568069
25.
Caputo
,
M.
,
2001
, “
Distributed Order Differential Equation Modelling Dielectric Induction and Diffusion
,”
Fract. Calc. Appl. Anal.
,
4
(
4
), pp.
421
442
.https://www.researchgate.net/publication/268489043_Distributed_order_differential_equations_modeling_dielectric_induction_and_diffusion
26.
Chechkin
,
A.
,
Gorenflo
,
R.
, and
Sokolov
,
I.
,
2002
, “
Retarding Subdiffusion and Accelerating Superdiffusion Governed by Distributed-Order Fractional Diffusion Equations
,”
Phys. Rev. E
,
66
(
4
), p.
046129
.10.1103/PhysRevE.66.046129
27.
Chechkin
,
A.
,
Gonchar
,
V.
,
Gorenflo
,
R.
,
Korabel
,
N.
, and
Sokolov
,
I.
,
2008
, “
Generalized Fractional Diffusion Equations for Accelerating Subdiffusion and Truncated Lévy Flights
,”
Phys. Rev. E
,
78
(
2
), p.
021111
.10.1103/PhysRevE.78.021111
28.
Li
,
Y.
,
Sheng
,
H.
, and
Chen
,
Y.
,
2011
, “
On Distributed Order Integrator/Differentiator
,”
Signal Process.
,
91
(
5
), pp.
1079
1084
.10.1016/j.sigpro.2010.10.005
29.
Li
,
Y.
, and
Chen
,
Y.
,
2014
, “
Lyapunov Stability of Fractional-Order Nonlinear Systems: A Distributed-Order Approach
,”
ICFDA'14 International Conference on Fractional Differentiation and Its Applications
, IEEE, Catania, Italy, pp.
1
6
.
30.
Duan
,
J.
, and
Baleanu
,
D.
,
2018
, “
Steady Periodic Response for a Vibration System With Distributed Order Derivatives to Periodic Excitation
,”
J. Vib. Control
,
24
(
14
), pp.
3124
3131
.10.1177/1077546317700989
31.
Bagley
,
R.
, and
Torvik
,
P.
,
2000
, “
On the Existence of the Order Domain and the Solution of Distributed Order Equations-Part i
,”
Int. J. Appl. Math.
,
2
(
7
), pp.
865
882
.https://www.researchgate.net/publication/306079170_On_the_Existence_of_the_Order_Domain_and_the_Solution_of_Distributed_Order_Equations_-_Part_I
32.
Kharazmi
,
E.
,
Zayernouri
,
M.
, and
,
G.
,
2017
, “
Petrov–Galerkin and Spectral Collocation Methods for Distributed Order Differential Equations
,”
SIAM J. Sci. Comput.
,
39
(
3
), pp.
A1003
A1037
.10.1137/16M1073121
33.
Kharazmi
,
E.
, and
Zayernouri
,
M.
,
2018
, “
Fractional Pseudo-Spectral Methods for Distributed-Order Fractional PDEs
,”
Int. J. Comput. Math.
,
95
(6–7), pp. 1340–1361.10.1080/00207160.2017.1421949
34.
Ding
,
W.
,
Patnaik
,
S.
,
Sidhardh
,
S.
, and
Semperlotti
,
F.
,
2021
, “
Applications of Distributed-Order Fractional Operators: A Review
,”
Entropy
,
23
(
1
), p.
110
.10.3390/e23010110
35.
Łabȩdzki
,
P.
,
Pawlikowski
,
R.
, and
,
A.
,
2018
, “
Transverse Vibration of a Cantilever Beam Under Base Excitation Using Fractional Rheological Model
,”
AIP Conference Proceedings
,
2029(1), p.
020034
.10.1063/1.5066496
36.
Ansari
,
R.
,
Oskouie
,
M. F.
, and
Gholami
,
R.
,
2016
, “
Size-Dependent Geometrically Nonlinear Free Vibration Analysis of Fractional Viscoelastic Nanobeams Based on the Nonlocal Elasticity Theory
,”
Phys. E Low Dimens. Syst. Nanostruct.
,
75
, pp.
266
271
.10.1016/j.physe.2015.09.022
37.
Faraji Oskouie
,
M.
,
Ansari
,
R.
, and
,
F.
,
2017
, “
Nonlinear Vibration Analysis of Fractional Viscoelastic Euler—Bernoulli Nanobeams Based on the Surface Stress Theory
,”
Acta Mech. Solida Sin.
,
30
(
4
), pp.
416
424
.10.1016/j.camss.2017.07.003
38.
Eyebe
,
G.
,
Betchewe
,
G.
,
,
A.
, and
Kofane
,
T.
,
2018
, “
Nonlinear Vibration of a Nonlocal Nanobeam Resting on Fractional-Order Viscoelastic Pasternak Foundations
,”
Fractal Fract.
,
2
(
3
), p.
21
.10.3390/fractalfract2030021
39.
Lewandowski
,
R.
, and
Wielentejczyk
,
P.
,
2017
, “
Nonlinear Vibration of Viscoelastic Beams Described Using Fractional Order Derivatives
,”
J. Sound Vib.
,
399
, pp.
228
243
.10.1016/j.jsv.2017.03.032
40.
Manley
,
G. A.
,
Gummer
,
A. W.
,
Popper
,
A. N.
, and
Fay
,
R. R.
,
2017
,
Understanding the Cochlea
, Vol.
62
,
Springer International Publishing
, Switzerland.
41.
Samiee
,
M.
,
Zayernouri
,
M.
, and
Meerschaert
,
M.
,
2019
, “
A Unified Spectral Method for FPDEs With Two-Sided Derivatives; Part I: A Fast Solver
,”
J. Comput. Phys.
,
385
, pp.
225
243
.10.1016/j.jcp.2018.02.014
42.
Samiee
,
M.
,
Zayernouri
,
M.
, and
Meerschaert
,
M.
,
2019
, “
A Unified Spectral Method for Fpdes With Two-Sided Derivatives; Part II: Stability, and Error Analysis
,”
J. Comput. Phys.
,
385
, pp.
244
261
.10.1016/j.jcp.2018.07.041
43.
Samiee
,
M.
,
Kharazmi
,
E.
,
Zayernouri
,
M.
, and
Meerschaert
,
M.
,
2018
, “
A Unified Petrov–Galerkin Spectral Method and Fast Solver for Distributed-Order Partial Differential Equations
,”
Comm. App. Math. Com.
, 3(1) pp. 61–90.10.1007/s42967-020-00070-w
44.
Lubich
,
C.
,
1986
, “
Discretized Fractional Calculus
,”
SIAM J. Math. Anal.
,
17
(
3
), pp.
704
719
.10.1137/0517050
45.
Zayernouri
,
M.
,
Ainsworth
,
M.
, and
,
G.
,
2015
, “
Tempered Fractional Sturm–Liouville Eigenproblems
,”
SIAM J. Sci. Comput.
,
37
(
4
), pp.
A1777
A1800
.10.1137/140985536
46.
Suzuki
,
J.
, and
Zayernouri
,
M.
,
2020
, “
A Self-Singularity-Capturing Scheme for Fractional Differential Equations
,”
Int. J. Comput. Math
.,
98
(
5
), pp. 933–960.10.1080/00207160.2020.1792453
47.
Zayernouri
,
M.
, and
Matzavinos
,
A.
,
2016
, “
Fractional Adams–Bashforth/Moulton Methods: An Application to the Fractional Keller–Segel Chemotaxis System
,”
J. Comput. Phys.
,
317
, pp.
1
14
.10.1016/j.jcp.2016.04.041
48.
Lin
,
Y.
, and
Xu
,
C.
,
2007
, “
Finite Difference/Spectral Approximations for the Time-Fractional Diffusion Equation
,”
J. Comput. Phys.
,
225
(
2
), pp.
1533
1552
.10.1016/j.jcp.2007.02.001
49.
Zhou
,
Y.
,
Suzuki
,
J.
,
Zhang
,
C.
, and
Zayernouri
,
M.
,
2020
, “
Implicit-Explicit Time Integration of Nonlinear Fractional Differential Equations
,”
Appl. Numer. Math.
,
156
, pp.
555
583
.10.1016/j.apnum.2020.04.006
50.
Mashayekhi
,
S.
,
Hussaini
,
Y.
, and
Oates
,
W.
,
2019
, “
A Physical Interpretation of Fractional Viscoelasticity Based on the Fractal Structure of Media: Theory and Experimental Validation
,”
J. Mech. Phys. Solids
,
128
, pp.
137
150
.10.1016/j.jmps.2019.04.005
51.
Mainardi
,
F.
,
2010
,
Fractional Calculus and Waves in Linear Viscoelasticity: An Introduction to Mathematical Models
,
World Scientific
, Imperial College Press, UK.
52.
Mainardi
,
F.
, and
,
G.
,
2011
, “
Creep, Relaxation and Viscosity Properties for Basic Fractional Models in Rheology
,”
Eur. Phys. J. Spec. Top.
,
193
(
1
), pp.
133
160
.10.1140/epjst/e2011-01387-1
53.
Meirovitch
,
L.
,
2010
,
Fundamentals of Vibrations
,
Waveland Press, Long Grove, IL.
54.
Bonet
,
J.
, and
Wood
,
R.
,
1997
,
Nonlinear Continuum Mechanics for Finite Element Analysis
,
Cambridge University Press
, New York.
55.
Lion
,
A.
,
1997
, “
On the Thermodynamics of Fractional Damping Elements
,”
Contin. Mech. Thermodyn.
,
9
(
2
), pp.
83
96
.10.1007/s001610050057
56.
,
E.
,
2012
, “
A Review of Numerical Methods for Nonlinear Partial Differential Equations
,”
B. Am. Math. Soc.
,
49
(
4
), pp.
507
554
.10.1090/S0273-0979-2012-01379-4
57.
Azrar
,
L.
,
Benamar
,
R.
, and
White
,
R.
,
1999
, “
Semi-Analytical Approach to the Non-Linear Dynamic Response Problem of s–s and c–c Beams at Large Vibration Amplitudes Part I: General Theory and Application to the Single Mode Approach to Free and Forced Vibration Analysis
,”
J. Sound Vib.
,
224
(
2
), pp.
183
207
.10.1006/jsvi.1998.1893
58.
Tseng
,
W.-Y.
, and
Dugundji
,
J.
,
1971
, “
Nonlinear Vibrations of a Buckled Beam Under Harmonic Excitation
,”
ASME J. Appl. Mech.
,
38
(
2
), pp.
467
476
.10.1115/1.3408799
59.
Loutridis
,
S.
,
Douka
,
E.
, and
,
L.
,
2005
, “
Forced Vibration Behaviour and Crack Detection of Cracked Beams Using Instantaneous Frequency
,”
NDT E Int.
,
38
(
5
), pp.
411
419
.10.1016/j.ndteint.2004.11.004
60.
Hamdan
,
M.
, and
,
M.
,
1997
, “
Large Amplitude Free Vibrations of a Uniform Cantilever Beam Carrying an Intermediate Lumped Mass and Rotary Inertia
,”
J. Sound Vib.
,
206
(
2
), pp.
151
168
.10.1006/jsvi.1997.1081
61.
Lestari
,
W.
, and
Hanagud
,
S.
,
2001
, “
Nonlinear Vibration of Buckled Beams: Some Exact Solutions
,”
Int. J. Solids Struct.
,
38
(
26–27
), pp.
4741
4757
.10.1016/S0020-7683(00)00300-0
62.
Eisley
,
J. G.
,
1964
, “
Nonlinear Vibration of Beams and Rectangular Plates
,”
Z. Angew. Math. Phys.
,
15
(
2
), pp.
167
175
.10.1007/BF01602658
63.
Hsu
,
C.
,
1960
, “
On the Application of Elliptic Functions in Non-Linear Forced Oscillations
,”
Q. Appl. Math.
,
17
(
4
), pp.
393
407
.10.1090/qam/110250
64.
Pillai
,
S.
, and
Rao
,
B. N.
,
1992
, “
On Nonlinear Free Vibrations of Simply Supported Uniform Beams
,”
J. Sound Vib.
,
159
(
3
), pp.
527
531
.10.1016/0022-460X(92)90756-N
65.
Evensen
,
D. A.
,
1968
, “
Nonlinear Vibrations of Beams With Various Boundary Conditions
,”
AIAA J.
,
6
(
2
), pp.
370
372
.10.2514/3.4506
66.
Wriggers
,
P.
,
2008
,
Nonlinear Finite Element Methods
,
Springer
, Berlin, Germany.
67.
Lin
,
R.
,
2018
, “
Comments on “Nonlinear Vibration of Viscoelastic Beams Described Using Fractional Order Derivatives
,”
J. Sound Vib.
,
428
, pp.
195
204
.10.1016/j.jsv.2018.05.015
68.
Svenkeson
,
A.
,
Glaz
,
B.
,
Stanton
,
S.
, and
West
,
B.
,
2016
, “
Spectral Decomposition of Nonlinear Systems With Memory
,”
Phys. Rev. E
,
93
(
2
), p.
022211
.10.1103/PhysRevE.93.022211
69.
Shoshani
,
O.
,
Shaw
,
S.
, and
Dykman
,
M.
,
2017
, “
Anomalous Decay of Nanomechanical Modes Going Through Nonlinear Resonance
,”
Sci. Rep.
,
7
(
1
), p.
18091
.10.1038/s41598-017-17184-6
70.
Nayfeh
,
A.
, and
Mook
,
D.
,
2008
,
Nonlinear Oscillations
,
Wiley
, Hoboken, NJ.
71.
Rossikhin
,
Y. A.
, and
Shitikova
,
M.
,
2010
, “
Application of Fractional Calculus for Dynamic Problems of Solid Mechanics: Novel Trends and Recent Results
,”
ASME Appl. Mech. Rev.
,
63
(
1
), p.
010801
.10.1115/1.4000563
72.
Samko
,
S.
,
Kilbas
,
A.
, and
Marichev
,
O.
,
1993
,
Fractional Integrals and Derivatives: Theory and Applications
,
CRC Press
, Switzerland; Gordon and Breach Science Publishers, Philadelphia, PA.
73.
Shen
,
Y.
,
Yang
,
S.
,
Xing
,
H.
, and
Gao
,
G.
,
2012
, “
Primary Resonance of Duffing Oscillator With Fractional-Order Derivative
,”
Commun. Nonlinear Sci.
,
17
(
7
), pp.
3092
3100
.10.1016/j.cnsns.2011.11.024