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.
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 , 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  and dynamic storage/dissipation in visco-elasticity . 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 . When subject to mechanical loads, such materials undergo microstructural changes, e.g., rearrangement/unfolding of polymer networks/chains , plastic stretching/buckling of microfibers , formation, arresting, relaxation of dislocations , 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 . 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 . In addition, multi-exponential approximations are mere truncations of power-law relaxation , only providing accurate results for short times, thus lacking predictability and requiring recalibration for multiple time-scales .
Fractional differential equations (FDEs) are predictive tools for anomalous materials across multiple time-scales. Nutting  and Gemant  showed that power law kernels are more descriptive of creep/relaxation. Bagley and Torvik  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 , that interpolates between Hookean and Newtonian 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 , frequency domain analysis , and uncertainty quantification [32,33]. We refer the reader to Ref.  for a thorough review on applications of DODEs.
Regarding the dynamics of anomalous beams, Łabȩdzki et al.  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.  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.  incorporated surface stress effects through the Gurtin–Murdoch theory. Recently, Eyebe et al.  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  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 . 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 . Among many schemes for time-fractional integration [44–47], for simplicity, we outline the L1 finite difference scheme by Lin and Xu  and refer to Ref.  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 . 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 α . 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.
Let an arbitrary neutral axis element CD with initial length ds and material coordinate s deform to configuration (see Fig. 2). The displacement components of points C and D are denoted by (u, v) and , respectively. The axial strain e(s, t) at C is given by . Applying the inextensionality assumption e = 0, we have . Under the latter equation and negligible vertical shear strains, we have . Using expansion , the curvature can be approximated up to third-order terms as . Therefore, the angular velocity and curvature of the beam can be approximated as and the curvature approximated as . Finally, under assumptions (1) and (2), the strain–curvature relationship takes the form .
2.2 Linear Viscoelasticity: Boltzmann Superposition Principle.
and denotes the strain rate.
2.2.1 Exponential Relaxation (Classical Models) Versus Power-Law Relaxation (Fractional Models).
which solves the ODE , where the relaxation time constant is experimentally obtained.
which are equivalent when .
2.2.2 Multipower Laws Versus Distributed-Order Models.
where denotes any definition of fractional derivative. The terms and denote distribution functions, where and are continuous mappings in and . 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.
2.3 Extended Hamilton's Principle.
We derive the equations of motion by employing the extended Hamilton's principle , where δT and δW denote the variations of kinetic energy and total work . The only external source of our system with volume 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 .
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 , 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.
2.4 Weak Formulation.
2.5 Assumed Mode: Spectral Space Approximation.
for all .
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.  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  obtained similar results between single- and two-mode approximations for nonlinear vibrations of clamped–clamped beams far from the crossover region. Loutridis et al.  implemented a crack detection method for beams using a single DOF system with time-varying stiffness. In Ref. , the effects of base stiffness and attached mass on the nonlinear, flexural free vibrations of beams were studied. Lestari and Hanagud  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.
3 Linearized Equation: Numerical Time Integration.
Since the Newmark method is -accurate, the overall accuracy is dominated by the L1 scheme, which is . We also note that a discretization of a Caputo-variant of Eq. (15) is recovered if we remove the term 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 , i.e., we utilize (D.7) from Ref.  for , which yields the coefficients . 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 in Eq. (15), where and denote, respectively, the base displacement frequency and amplitude. The value is obtained by Eqs. (16) and (14). We set , and time , with step-size , 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 that changes the dissipation nature of the fractional order. Regarding the maximum observed amplitudes, increasing the fractional order in the range decreases and slightly shifts the amplitude peaks to higher (right) frequencies (an anomalous quality). Conversely, as α increases in the range , 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 , 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 , and . 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 . In addition, an anomalous transient region is observed at the short time-scale . 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. . Finally, anomalous decay rates were also investigated in Ref.  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.
where a and are obtained by real/imaginary parts separation.
4.1.1 Case 1: No Lumped Mass at the Tip.
We let , and thus, we employ (D.8) from Ref. for the eigenfunctions , from which we compute , and . We consider the following cases:
Free vibration, F = 0: Super Sensitivity to αThe governing equations for amplitude and phase are given by(26)Note that the amplitude decays with rate in Eq. (26), which depends on the pair (Fig. 6). Let the sensitivity index , defined as(27)which is computed and plotted in Fig. 7 for the same set of parameters. There exists a critical value , where . We observe that by increasing α when (adding viscosity), the dissipation rate, and thus decay rate, increases; this can be interpreted as a softening (stiffness-decreasing) region. Further increasing α when 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 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 . 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 , given byfrom which the real part yields the storage modulus and the imaginary part yields the loss modulus , which represent, respectively, the stored and dissipated energies per loading cycle. Finally, we define the tangent loss , which represents the ratio between the dissipated/stored energies, and hence related to the mechanical damping of the anomalous medium. We set and Er = 1 and demonstrate the results for with varying fractional orders α. In Fig. 8 (top), increasing values of α led to increased values of the tangent loss, and the hardening part () is not associated with higher storage in the material. Instead, the increasing dissipation with α suggests an increasing damping of the mechanical structure.(28)
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  and set . We also assume the following piecewise strain function: , for (monotone stress/strain), and for (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,Let , where Δ is called the detuning parameter and thus, we write the force function as . In this case, the force function also contributes to the secular terms. Therefore, we find the governing equations of solution amplitude and phase as(29)in which the parameters mainly change the frequency response of the system. Equations (29) and (30) can be transformed into an autonomous system by letting . The steady-state solution occurs when , which yields, after algebraic manipulations(30)(31)
which is a cubic equation in a2. The discriminant of a cubic equation of the form is given as . Equation (32) has one real root when and three distinct real roots when . The parameters dictate the values of coefficients , the value of ϑ, and thus the number of admissible steady-state amplitudes. We see that for fixed values of , 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 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 ). The peak amplitudes decrease as α increases. Figure 10 shows the frequency response of the system for different values of 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 , including Duffing oscillators  and damage-induced softening .
4.1.2 Case 2: Lumped Mass at the Tip.
We let , and employ eigenfunctions (D.7) from Ref. , which yields , and . 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(33)(34)
Primary resonance case,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 , and yields the steady-state amplitude of the system.(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 . 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.
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).