Anomalous Nonlinear Dynamics Behavior of Fractional Viscoelastic Beams

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 [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.

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

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.

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$⁠.

and $ε̇$ denotes the strain rate.

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

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)$⁠.

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.

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.

for all $ṽN∈ṼN$⁠.

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.

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.

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)$⁠.

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 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 $ϕ(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.

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.

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.

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.

where a and $φ$ are obtained by real/imaginary parts separation.

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 ω₀ 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 E_r = 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)