## Abstract

The modeling of nonlinear dynamical systems subject to strong and evolving nonsmooth nonlinearities is typically approached via integer-order differential equations. In this study, we present the possible application of variable-order (VO) fractional operators to a class of nonlinear lumped parameter models that have great practical relevance in mechanics and dynamics. Fractional operators are intrinsically multiscale operators that can act on both space- and time-dependent variables. Contrarily to their integer-order counterpart, fractional operators can have either fixed or VO. In the latter case, the order can be function of either independent or state variables. We show that when using VO equations to describe the response of dynamical systems, the order can evolve as a function of the response itself; therefore, allowing a natural and seamless transition between widely dissimilar dynamics. Such an intriguing characteristic allows defining governing equations for dynamical systems that are evolutionary in nature. Within this context, we present a physics-driven strategy to define VO operators capable of capturing complex and evolutionary phenomena. Specific examples include hysteresis in discrete oscillators and contact problems. Despite using simplified models to illustrate the applications of VO operators, we show numerical evidence of their unique modeling capabilities as well as their connection to more complex dynamical systems.

## 1 Introduction

Fractional calculus (FC) is the study of differential and integral operators of noninteger-order. Although the branch of fractional calculus was created almost simultaneously to its integer-order counterpart, the mathematics and its applications are considerably less developed. Due to their differ-integral nature, fractional operators are intrinsically multiscale. While time fractional operators enable memory effects (i.e., the response of a system is a function of its past history), space fractional operators can account for nonlocality and scale effects. In recent years, there has been a surge of interest in fractional-order operators and in their applications to the simulation of physical problems. Among the areas that have seen the largest number of applications, we mention transport processes in complex media [14], formulation of constitutive equations for viscoelastic materials [57], nonlocal elasticity [811], and model-order reduction of lumped parameter systems [12,13]. In all these studies, fractional operators with constant-order (CO) have typically been used. Comprehensive reviews on constant-order fractional calculus can be found in Refs. [14] and [15].

As fractional-order operators can be seen as an analytical continuation of integer-order operators, variable-order (VO) operators can be seen as the natural extension of CO fractional operators. In other terms, unlike the integer-order that can only vary in discrete steps, the fractional-order allows a continuous variation with any arbitrary step. Furthermore, in VO operators, the order can vary either as a function of an independent variable of integration (or differentiation) or as a function of some other dependent variable. These properties of VO calculus provide the basis for the development of a computational framework with potentially unlimited possibilities in terms of modeling complex physical phenomena. As an example, the reaction kinetics of proteins has been found to experience fractional-order relaxation mechanisms. Glöckle and Nonnenmacher [16] established that the fractional-order has a temperature dependence and defined a variable-order governing equation where the order was defined to be a function of the temperature.

Although the extension from CO to VO operators may seem somewhat natural, the first comprehensive discussion on these operators was only recently given by Samko and Ross [17,18]. In VO operators, the order can be function of time, space, or even of an independent external variable (e.g., temperature or applied loads). The formalism introduced in Refs. [1719] has led to research on the application of VO operators to the modeling of anomalous diffusion in complex structures [20,21]. Coimbra [22] has used VO operators to model oscillators under nonlinear viscoelastic forces. Diaz and Coimbra [23] have investigated the dynamics and control of a VO Van der Pol oscillator. All these studies recognized and took advantage of the intrinsic memory capability of VO operators (see Ref. [19]), and of the way, this property could be leveraged to describe more accurately the dynamics of nonlinear systems. The interested reader can find a comprehensive review of the applications of VO operators in Ref. [24].

In this study, we show how one of the most remarkable properties of VO-based physical models consists in their evolutionary nature. We will also show how this property can play a critical role in the simulation of nonlinear dynamical systems. More specifically, as the VO formalism allows updating the system's order depending on its instantaneous response (and, potentially, its historical response), the same theoretical model can evolve seamlessly to describe widely dissimilar dynamics without the need for changing the underlying governing equation. In the following, we will provide examples and applications of this remarkable property to nonlinear dynamics with particular attention to contact problems and hysteretic systems. Further, as pointed out in Refs. [22] and [23], when using VO fractional calculus (VO-FC), the nonlinear behavior of constant-order differential equations can often be modeled via linear operators. The most immediate effect of this property is the possibility to extend, under certain conditions, many of the analysis tools available for linear systems to the nonlinear ones.

An important step to promote the use of fractional-order models for the simulation of complex systems is to establish the connection between the physical (e.g., material parameters) and the mathematical (e.g., the law of variation of the order) properties. Despite the increasing amount of research dedicated to exploring this specific aspect, a comprehensive approach is still lacking. In this study, we start filling this gap by presenting physics-driven constitutive laws for order variations for specific types of nonlinear problems that are of great interest in mechanics.

The remainder of the paper is structured as follows. First, we introduce briefly the VO operators based on the Riemann–Liouville (RL) definition. Next, we discuss how VO operators can be formulated in order to enable evolutionary governing equations. Then, we present the application of this approach to the modeling of contact dynamics and nonlinear hysteretic response. We will discuss how the order of the VO operator can capture transitions in both the operating regime (e.g., from linear to nonlinear) and the underlying physical phenomena (e.g., contacts and hysteresis).

## 2 Variable-Order Fractional Operators

In their seminal study on VO-FC, Lorenzo and Hartley [19] presented several time-domain definitions of VO-RL fractional operators. The discriminant factor between the different definitions consisted in the memory behavior of the fractional operator. Here, we report only the definition that will be used in this work.

Definition 1: If g(t) and$β(t)$are continuous real-valued functions on$(a0,t)$, the left-handed VO-RL differential operator to the order$β(t)>0$with the lower bound a is defined as
$aRLDtβ(t)g(t)=dmdtm[1Γ(m−β(t))∫atg(τ)(t−τ)1+β(t)−mdτ]+dmdtmψ(g,−q,a0,a,t)+ψ(h,m,a0,a,t)$
(1)
where$a0≤ais the upper integer bound on$β(t), q=m−β(t)$, and$Γ(·)$is the Gamma function. The VO-RL operator is initialized at a0 such that$g(t)=0 ∀ t≤a0≤a$. $ψ(g,−q,a0,a,t)$is the initialization function defined as
$ψ(g,−q,a0,a,t)=1Γ(q)∫a0a(t−τ)q−1g(τ)dτ$
(2)
The term h is defined as
$h=a0RLDt−qg(t)=1Γ(q)∫a0t(t−τ)q−1g(τ)dτ$
(3)

It is well known that the use of RL operators in fractional-order differential equations requires fractional-order boundary conditions whose physical interpretation is more elusive than their integer-order counterpart [14 ]. The initialization function$ψ(t)$shifts the initial time instant (from a to a0) of the interval over which the fractional derivative is defined. As$g(t)=0 ∀ t≤a0≤a$, all the fractional-order derivatives of g(t) at a0 are zero. Thus, we do not require fractional-order boundary conditions while solving fractional-order differential equations with the initialized RL derivative. In other terms, the function$ψ(t)$accounts for the effect of the history of g(t) and its derivatives up to the order m over the interval$[a0,a]$hence allowing bypassing the use of fractional boundary conditions. Mathematically,$ψ(t)$brings to the definition of the fractional derivative the effect of fractionally differentiating g(t) from a0 to a [25 ]. Further details on the initialization procedure and its importance can be found inRef. [25 ].

## 3 Evolutionary Governing Equations: The Role of Variable-Order-Riemann–Liouville Operators Applied to Constants

It is well-known that different definitions of a fractional derivatives are not all equivalent to each other. Differences are particularly pronounced when they are calculated at the bounds of the domain of integration or when the argument is a constant. The only strict requirement for every definition of a fractional derivative is that its value should coincide with the value of the corresponding integer-order derivative when $β(t)=m$, where $m∈N$ (the set of integers). Within this context, an example of particular interest for the following study consists of the fixed-order RL derivative of a constant which is not equal to zero, unless $β(t)=m$. Although, at first, this characteristic might appear unsettling, particularly in light of classical integer-order calculus, we will show that it is a very convenient property to model several nonlinear and discontinuous behaviors in dynamics.

In order to provide the necessary mathematical context, we first present the RL derivative of a constant function. Consider the function g(t) such that $g(t)=0 ∀ t≤0$, and $g(t)=C0 ∀ t>0$, where C0 is an arbitrary constant. Note that we take $g(t)=0 ∀ t≤0$ for proper initialization. The RL derivative of g(t) over the interval $(0,t)$ to a constant fractional-order β0 is given by [14]
$0RLDtβ0C0=C0 t−β0Γ(1−β0)$
(4)

where $m−1≤β0.

Now, consider a piecewise continuous function $β(t)$ defined via a continuous real-valued function s(t) on the domain $(0,t)$ as
$β(t)=exp(−s0 s(t))$
(5)

where the function s(t) is designed to capture the physical mechanism producing the order variation, and $s0∈ℝ+$ is a scaling factor that allows calibrating the order variation based on the characteristic response of the physical system. A detailed discussion of the procedure followed to select both the function s(t) and the numerical value of s0 for the different physical problems is reported in each specific section.

With s0 being a properly chosen constant, the limiting behavior of the order $β(t)$ is found to be

$β(t)→∞ ∀ s(t)<0$
(6a)
$β(t)→0 ∀ s(t)>0$
(6b)

The function $β(t)$ is now constant on intervals determined according to Eq. (6). We apply a VO-RL operator, denoted by $0RLDtβ(t)(·)$, to the constant function g(t) over the interval $(0,t)$. Equations (4)(6) lead to
$0RLDtβ(t)g(t)={limβ(t)→∞C0 t−β(t)Γ(1−β(t)) s(t)<0limβ(t)→0C0 t−β(t)Γ(1−β(t)) s(t)>0$
(7)
which under the limiting conditions simplifies to
$0RLDtβ(t)C0={0 s(t)≤0C0 s(t)>0$
(8)

Note that when s(t)=0, that is when the function changes sign, the VO-RL derivative of the constant function is exactly equal to zero as noted by the equality $s(t)≤0$ in Eq. (7).

Before applying the above concepts to mechanical systems, we present a purely illustrative example in order to provide some context on the methodology discussed above. Consider the VO $β(t)$ with $s(t)=t(t−3)(t−3.5)(t−4)(t−4.5)(t−5)$, that is a function that changes sign at specific time instants. The VO-RL operator $0RLDtβ(t)C0$ with $C0=2$ is given in Fig. 1. In this example, we assumed $s0=105$. This formulation allows the term $0RLDtβ(t)C0$ to switch between 0 and 2 when s(t) changes its sign.

Fig. 1
Fig. 1
Close modal

It is exactly this switching behavior that can be exploited to simulate certain nonlinear dynamical properties of mechanical systems. More specifically, consider defining VO operators as a part of a governing equation such that its variation can capture changes in the properties of the dynamical systems such as, for example, a change in stiffness. The change in stiffness could be abrupt, such as the type associated with closing contacts between multiple bodies, or it could be smooth, such as the one associated with the transition from a linear elastic to a nonlinear plastic material behavior due to large deformations. In all these cases, the response of the system changes from initially linear to, potentially, highly nonlinear. This change in the underlying dynamics of the system is captured via the VO $β(t)$, through the function s(t). It appears that the change in the VO $β(t)$ results in an implicit reformulation of the underlying equations of motion following a change in the underlying mechanisms dominating the response of the system. Note that this approach marks a significant departure from the traditional modeling of nonlinear systems where the functional (nonlinear) form of the equation as well as the existence of nonlinearities must be assumed a priori.

## 4 Variable-Order-Fractional Calculus for the Simulation of Nonlinear Dynamic Systems

In order to illustrate the unique capabilities of VO operators and of the corresponding governing equations, we introduce two case studies that are particularly relevant for the simulation of many dynamical systems. More specifically, the two applications focus on either modeling the dynamics of contacts between two impacting bodies or nonlinear irreversible (hysteretic) behavior.

The analysis of contacts between physical components is important in several real world applications such as, for example, structural analysis of couplings between the rolling stock of trains [26,27], dynamic analysis of loosely jointed structures [28,29] and breathing cracks [3032], and detection of loose bolted joints [33,34] or delamination in composite materials [35]. These systems are characterized by a bilinear change in the local stiffness which typically varies as a function of time, depending on the nature of the excitation and of the dynamic response of the system. In a similar way, there are many systems in which the material properties evolve from linear elastic to nonlinear plastic. Polymeric foams used extensively for impact attenuation in the automobile industry, and steel, used in structural members, exhibit hysteretic behavior [3638]. Hysteresis is also found, as an example, in systems subject to random vibrations [39,40], and in metal cutting [41]. Depending on the level of maximum and accumulated strain, material's behavior can evolve from linear elastic to hardening plastic, or linear elastic to nonlinear plastic.

Currently, available approaches to these classes of problems are based on nonlinear integer-order differential models which are typically solved by finite element techniques. Although integer-order models are invaluable tools for analysis, they also exhibit an important limitation which consist of their inability to evolve between different governing equations. Their ability to account for nonlinearities must be integrated in the model a priori often requiring somewhat arbitrary considerations on the elements that will experience the nonlinear behavior. On the contrary, VO-FC provides a natural platform to approach these problems due to its ability to formulate governing equations capable of evolving according to the system response.

In the following, we show how the characteristic behavior of the VO-RL operators when applied to constants can be used to simulate challenging contact problems. More specifically, we consider two different configurations: (1) dynamic contact between two lumped masses, and (2) dynamic contact between a vibrating beam and a finite stiffness boundary. Further, we present the application of the nonlinear behavior of VO-RL operator to the modeling of nonlinear hysteresis. We show that the order variation can be crafted to represent any linear or nonlinear material law such as those describing linear or nonlinear hysteretic behavior.

### 4.1 Contact Dynamics for Discrete Particles.

Consider a system of two oscillators separated by a dead zone of width d as shown in Fig. 2. Given the system of oscillators, the absolute displacement of the individual masses is indicated by $ξ1(t)$ and $ξ2(t)$, respectively. The two oscillators are connected to a spring–damper system fixed at one end and are driven by externally applied forces $F1(t)$ and $F2(t)$. The oscillators are separated by a dead-zone (i.e., a gap) of initial amplitude d such that the spring (K0) and damper (C0) are engaged only when $ξ1(t)−ξ2(t)>d$.

Fig. 2
Fig. 2
Close modal
Following the application of external loads, the dynamics of the two-oscillator system could be either linear or nonlinear depending on the status of the contact (i.e., whether the gap between them closes or not). More specifically, when the amplitude of oscillation is such that the gap is cleared, the two oscillators enter in contact with each other. Before the gap closes, the two oscillators are decoupled and exhibit a linear dynamic response. When the gap closes, the two oscillators couple in a nonlinear fashion. In the following, we show that the order of the VO-RL operator can be crafted appropriately so that the stiffness term Keq (see Fig. 2) can capture exactly the effect of the closing gap and, hence, the occurrence of the contact. The equations of motion of the system written in classical integer-order form are
$[M100M2]{ξ1˙˙ξ2˙˙}+[C100C2]{ξ1˙ξ2˙}+[K100K2]{ξ1ξ2}={F1F2}$
(9a)
when $ξ1−ξ2≤d$, and
$[M100M2]{ξ1˙˙ξ2˙˙}+[C0+C1−C0−C0C0+C2]{ξ1˙ξ2˙}+[K0+K1−K0−K0K0+K2]{ξ1ξ2}={F1F2}$
(9b)

when $ξ1−ξ2>d$.

We reformulate the above equations of motion using VO-FC. We start by defining the function s(t) that allows capturing the characteristic feature of the underlying contact mechanism (i.e., the gap width) as a function of the system's response. The function is defined as
$s(t)=ξ1(t)−ξ2(t)−d$
(10)
Note that the function $s(t)=ξ1(t)−ξ2(t)−d$ changes sign when the contact either closes or opens, hence it appears to be a reasonable choice to track the physical response of the contact. As previously mentioned, this function provides the link between selected physical features of the system and the VO $β(t)$ of the operator and it could potentially be defined in many different ways. Hence, this is one of the possible choices for such function. The variable-order $β(t)$ can now be defined according to Eq. (5) as
$β(t)=exp(−s0(ξ1(t)−ξ2(t)−d))$
(11)

The scaling factor s0 can be calculated as follows. We define the infimum of the function s(t) in the neighborhood of a point where the function switches sign as $ε=inf(|s(t)|)$. Now, $∀ ε>0$ and $δ>0 ∃ s0$ such that $e−s0ε<δ$; δ can be chosen to be infinitesimally small. Note that s0 depends on both the ε and δ parameters. We have graphically demonstrated the significance of the parameters δ and ε in Fig. 1.

From a general perspective, the parameter δ provides the characteristic time (i.e., the temporal resolution) for the switch to occur. In other terms, it determines the temporal range within which the switch is bound to occur. The parameter $ε=inf(s(t))$ physically represents the smallest possible change in the function s(t) following a change in the sign of s(t) that can be accurately detected. In other terms, the parameter ε controls the spatial resolution of the function s(t) within the interval defined by δ. For problems involving rapid changes in properties, such as contact problems, the switch should occur over a length scale much smaller than the characteristic spatial scale of the problem. The characteristic length scale is clearly dependent on the specific problem, hence the need for the scaling factor s0.

From a more practical perspective, the numerical value of the parameter s0 can be estimated as follows. First, we choose the parameter δ which provides the characteristic time (i.e., the temporal resolution) for the switch to occur. Note that δ should be greater than or equal to the time-step $(Δt)$ used to numerically simulate the system; the most accurate value being $δ=Δt$. As an example, for the simulations presented in the following for the contact problem, we chose $Δt=10−3$ s; hence, $δ=10−3$. Then, we define the parameter $ε=inf(|s(t)|)$, which physically represents the smallest possible change in the relative displacement of the two masses following a change in the status of the contact. As previously discussed, this parameter is related to the characteristic spatial scale of the problem. For contact problems, the initial gap distance d can be chosen as a possible characteristic spatial scale. Then, the spatial resolution for the switch can be set by requesting $ε≤P0d$, where P0 is a user-defined percentage value. At this point, it is easy to see that if tk is the instant when the contact status changes (i.e., $s(tk)=0$), then at the time-step $t̃k=tk+Δt$ following directly the change in the contact status, we have $s(t̃k)=ξ1(t̃k)−ξ2(t̃k)−d=inf(s(t))≠0$. Given the definition of the spatial resolution, the switch in the status of the contact should occur within the interval $0. For the simulations presented later in this paragraph, we chose $P0=10−2$ which means that the switch occurs when the function s(t) exceeds a value equal to 1% the initial gap distance d. From the relation $e−s0ε<δ$, we obtain $s0>1.4×102$. For a more conservative approximation, the value was rounded to $s0=103 m−1$.

At this point, we have all the information to define Keq and Ceq as VO-RL operators. Using Eqs. (8) and (11), it can be shown that Keq and Ceq are

$Keq=0RLDtβ(t)K0={0ξ1−ξ2≤dK0ξ1−ξ2>d$
(12a)
$Ceq=0RLDtβ(t)C0={0ξ1−ξ2≤dC0ξ1−ξ2>d$
(12b)

where K0 and C0 are constant values of stiffness and damping. We use Keq and Ceq to combine the equations of motion (see Eq. (12)) of the coupled oscillators in a single set of equations
$[M100M2]{ξ1˙˙ξ2˙˙}+[Ceq+C1−Ceq−CeqCeq+C2]{ξ1˙ξ2˙}+[Keq+K1−Keq−KeqKeq+K2]{ξ1ξ2}={F1F2}$
(13)

where Keq and Ceq evolve according to Eq. (12) guided by the evolution of the VO $β(t)$. The nonlinearity associated with the contact between the two masses is now fully captured in the VO. The equations of motion of the coupled oscillators are simplified, at least formally, and can evolve from linear to nonlinear simply based on the response of the system (i.e., as a function of the distance between the two masses). This discussion highlights a very unique feature of VO operators, that is, their evolutionary nature.

Another particularly intriguing aspect of this formulation is that the operator itself is still linear; a property inherited by the intrinsic nature of VO-FC [19,22]. However, we emphasize that the linearity of the VO-RL operator depends on the specific functional variation of the order [19]. In the context of our study, the fractional operator acts on piecewise constant functions (i.e., the contact stiffness and contact damping). The system consists of two oscillators and a set of spring-damper (K0 and C0) inserted between them. In this specific case, the application of the principle of superposition would correspond to adding an additional set of spring-damper (say, $K0′$ and $C0′$) between the two oscillators. For this class of operators, it is easy to show that
$0RLDtβ(x1,x2)[α1ϕ1+α2ϕ2]=α10RLDtβ(x1,x2)ϕ1+α20RLDtβ(x1,x2)ϕ2$
(14)

where $ϕ1$ (and $ϕ2$) is either K0 (and $K0′)$ or C0 (and $C0′)$.

Based on the above formulation, the equivalent VO-coupled oscillators system is schematically shown in Fig. 2(b). Note that the fractional Eq. (13) is an exact reformulation of the integer-order Eq. (9). Thus, we expect that the solutions obtained from both formulations were exactly equivalent (clearly, assuming the same initial and forcing conditions).

In order to validate this approach, we performed numerical simulations using the fractional-order system represented by Eq. (13) and compared the results with the solution obtained from the original integer-order system (Eq. (9)). Both systems of equations were solved numerically by using finite difference methods. The parameters for the simulations were chosen as $M1=M2=1$ kg, $C0=0, C1=0.1$ N·s/m, $C2=0.4$ N·s/m, $K0=2$ N/m, $K1=K2=1$ N/m, and d =5 m. Initial conditions were set at $ξ1(0)=0, ξ1˙(0)=−1$ m/s, and $ξ2(0)=0, ξ2˙(0)=0$. As previously discussed, the value of s0 used in the simulations was $s0=103 m−1$.

In order to showcase the capability of the VO operator applied to contact problems, we identified three different regimes, each corresponding to different values of the external forcing load: (1) (case a—linear) the gap does not close, hence the contact does not occur and the masses are effectively decoupled at all times; (2) (case b—weakly nonlinear) the gap closes aperiodically, and the overall response is quasi-periodic; (3) (case c—strongly nonlinear) the contact closes periodically hence giving rise to a highly nonlinear response. For case (a), we set $F1(t)=−3 sin(0.3t)$ N and $F2=0$; for case (b), $F1(t)=−10 sin(0.3t)$ N and $F2=0$; and for case (c), $F1(t)=5 sin(0.8t)$ N and $F2=−5 sin(0.8t)$ N. The results are shown in Fig. 3 in terms of the time history of the displacement of the two masses. The initial transient is very short and the results obtained in all the cases converge quickly to the steady-state. We denote the results obtained from the fractional-order model as $ξ1β$ and $ξ2β$ in order to differentiate them from the results of the integer-order model. From a general perspective, the results are in an excellent agreement with each other, and the nonlinearity in the contact problem is now fully captured in the order of the VO-RL operators (as evident from the switch in Keq).

Fig. 3
Fig. 3
Close modal

In case (a), the dead-zone between the masses #1 and #2 does not close, and hence, the oscillators are decoupled at all times. Thus, the response of the system to the sinusoidal forcing is linear as clearly seen in Fig. 3(a). As the amplitude of the loading increases (case (b)), the masses come in contact aperiodically and the response of the system becomes nonlinear but still quasi-periodic in nature; under these conditions, the system exhibits features typical of a weakly nonlinear dynamics. If the forcing amplitude is further increased, the masses come in contact periodically and the response becomes highly nonlinear (Fig. 3(c)). Also, note that each plot reports the time history of the spring constant Keq which clearly switches between the values 0 and K0 depending on the occurrence of contact.

We make a few remarks concerning the numerical scheme adopted to solve the VO fractional differential equations. From a general standpoint, it is not possible to obtain an analytical expression for the VO-RL derivative of a (nonconstant) function. Thus, the use of VO-RL derivatives in fractional differential equations requires numerical methods specialized for fractional-order derivatives. In this regard, there are several methods (e.g., finite difference, finite elements) that have been developed in recent years and that can be used for this purpose [24,42]. However, in the context of this study, we focused primarily on the use of VO-RL derivatives of constant functions. Under these specific conditions, it is possible to obtain an analytical expression (see Eq. (8)) of the derivative. It follows that specialized numerical methods are not required for this study. Indeed, the same second-order finite difference numerical scheme was used for both integer- and fractional-order models. As a result, both the computational time and the numerical scheme complexities are exactly equivalent for both models.

In addition, the comparison between the results obtained from both the integer- and the fractional-order formulations is exact (i.e., the error is exactly zero) at every time instant. This comment should be interpreted in the following way. Given that (1) both formulations have been numerically solved using the same finite difference scheme, and that (2) the VO-RL allows an exact definition of the switch, then the error between the two solutions is numerically zero. At the same time, both solutions to the initial systems of differential equations were obtained using a central-difference scheme in time. In this regard, both solutions are approximate. The error corresponding to the central-difference scheme is $O(Δt2)$, where $Δt$ indicates the step used for the temporal discretization. As a result, the error in the numerical results obtained via the VO-RL formulation is also $O(Δt2)$.

Despite the simplicity of this discrete two-body system, one can envision how this VO-FC formulation can be easily extended to the analysis of systems with N interacting bodies. The classical integer-order formulation of such a system will require a maximum of $2N−1$ separate equations that are N × N dimensional depending on the contacts allowed between individual bodies, thus resulting in a large number of equations. The VO-FC formulation, instead, allows simulating the entire system response via a single N × N dimensional equation where the contact parameters are VO-RL operators. Clearly, the VO-FC formulation has the potential to greatly simplify the analysis of nonlinear systems involving contact between many interacting bodies.

### 4.2 Contact Dynamics for Flexible Beams.

We extend the VO-RL formalism presented in Sec. 4.1 to simulate the dynamics of a cantilever beam undergoing a possible contact at its free end. We will use an approach exactly equivalent to that illustrated in the previous paragraph for the discrete system. More specifically, the beam is clamped at one end and separated by a dead-zone of amplitude Δ from a finite stiffness boundary (modeled via a spring) at the other end. When the beam is subject to an external excitation, its free end can clear the gap and close the contact with the spring element, as shown in Fig. 4(a). Similar to the previous problem involving coupled oscillators, when the contact closes the beam and the spring couple in a nonlinear fashion. The response of this system is governed by

$EI∂4η∂ξ4+ρA∂2η∂t2=FT(ξ,t)$
(15a)
$η|0,t=0 ∂η∂ξ|0,t=0 EI∂2η∂ξ2|L,t=0$
(15b)
$EI∂3η∂ξ3|L,t={0 η(L,t)>−ΔK0η(L,t) η(L,t)≤−Δ$
(15c)

where $η(ξ,t)$ indicates the transverse displacement of the beam, $FT(ξ,t)$ is the applied transverse load, while E, I, L, and A are the Young's modulus, the area moment of inertia, the length, and the cross-sectional area of the beam, respectively. K0 is the stiffness of the spring located at the contact point.

Fig. 4
Fig. 4
Close modal
By following the approach outlined in Sec. 4.1, we reformulate the equations of motion for the beam using VO operators. Consider the following expression for the VO $β(t)$:
$β(t)=exp [s0(η(L,t)+Δ)]$
(16)
where s0 is a scaling factor as discussed in the previous problem, $s(t)=−(η(L,t)+Δ)$ is the function that captures the main features of the contact and that changes its sign depending on the status of the contact. When the VO $β(t)$ is used within a VO-RL operator applied to the spring stiffness K0, it results in a nonlinear term that can capture the contact between the tip of the beam and the spring. Mathematically, using Eqs. (8) and (11), it can be shown that
$Keq=0RLDtβ(t)K0={0 η(L,t)>−ΔK0 η(L,t)≤−Δ$
(17)

We use the above VO-RL term Keq to reformulate the equation of motion of the beam in the following manner:

$EI∂4η∂ξ4+ρA∂2η∂t2=FT(ξ,t)$
(18a)
$η|0,t=0 ∂η∂ξ|0,t=0 EI∂2η∂ξ2|L,t=0$
(18b)
$EI∂3η∂ξ3|L,t=Keqη(L,t)$
(18c)

where the VO-RL term Keq evolves according to Eq. (17) guided by the VO $β(t)$. Clearly, the nonlinearity associated with the contact between the beam and the spring is now fully captured in the VO-RL term. Based on the above formulation, the equivalent VO system is schematically shown in Fig. 4(b). We emphasize that the modified VO Eq. (18) is an exact reformulation of the integer-order Eq. (15). In order to validate the VO-RL formulation, we solve the fractional model numerically and compare with the response from a nonlinear integer-order model.

To numerically evaluate the response of the beam, we discretize it into N elements having uniform length and mass lumped at the ends, as shown in Fig. 4(c). To ensure compatibility, it is essential that both the transverse displacement and the rotation (i.e., the first derivative of the transverse displacement) are continuous across these lumped elements. Thus, the degrees-of-freedom of each element are the transverse displacements and rotations of the ends of the element, denoted as $(η1,η2)$ and $(∂η1/∂ξ,∂η2/∂ξ)$, respectively. These degrees-of-freedom can be [43] interpolated using C1 continuous Hermitian shape functions to obtain the following element mass and the stiffness matrices
$Me=ρAle420[156−22le5413le−22le4le2−13le−3le254−13le15622le13le−3le222le4le2]$
(19a)
$Ke=2EIle3[6−3le−6−3le−3le2le23lele2−63le63le−3lele23le2le2]$
(19b)
where le is the length of each discretized element. The element mass and stiffness matrices are then assembled to generate the following equation of motion for the entire beam
$[M]{η¨}+[K]{η}={FT}$
(20)
where $[M]$ and $[K]$ are the global mass and stiffness matrices of the beam, respectively, and ${FT}$ is the force vector. The vector ${η}$ contains both the displacement and the rotation degrees-of-freedom. The variable stiffness at the end of the beam is accounted for by adding the VO-RL term Keq to the diagonal position of the global stiffness matrix $[K]$ corresponding to the transverse degree-of-freedom of the tip of the beam. Since the VO-RL term Keq evolves with time, the stiffness matrix $[K]$ evolves with time as well. In order to obtain the time response of the coupled system, the spatially discretized Eq. (20) is further discretized in time. By using a uniform time discretization, we obtain the following set of algebraic equations describing the response of the coupled beam-spring system
${η}m+1=2{η}m−{η}m−1+Δt2[M]−1[FTm−[K]m{η}m]$
(21)

where $Δt$ is the size of the time-step and the superscript m denotes the time-step such that at the step m time is $t=mΔt$. Note that this numerical approach can be applied to both the integer and the fractional-order models because the fractional-order term representing the contact is expressed entirely analytically. Hence, as in the previous example, there is no need for specialized numerical methods for the solution of the fractional-order problem.

By using the above described numerical model, the numerical response of the coupled system was obtained assuming the following parameters: EI =1 N·m2, ρA =1 kg/m, L =1 m, $K0=1$ N/m, $Δ=0.03$ m. The scaling parameter s0 was calculated according to the same procedure highlighted in the previous paragraph and it was found to be $s0=105m−1$ for a $P0=0.01$. The beam was discretized into N =25 elements. The initial condition of the beam was set at $η(ξ,0)=0$ and $η˙(ξ,0)=0$. The external transverse loading was chosen as $FT(ξ,t)=[0.1×δ(ξ−L)sin(t)]$ N and applied only at the tip of the beam. The results are presented in Fig. 5(a) in terms of the time history of the displacement of the tip of the beam. We have denoted the results obtained from the modified fractional equation of motion as $ηβ$ in order to differentiate it from the results obtained via the integer-order equation of motion.

Fig. 5
Fig. 5
Close modal

Additionally, we have also compared the response of the entire beam under an external loading $FT(ξ,t)=[δ(ξ−L)sin(4t)]$ N in order to highlight the effect of the contact on the curvature of the beam. The parameters used for this simulation are the same given above, except for Δ which is taken as 0.1 m. We obtained the response of the entire beam at two instants of time: (a) when the beam does not couple with the spring (denoted as $η1β$), and (b) when the beam couples with the spring (denoted as $η2β$). These two conditions occur at the time instants t =0.75 s and t =2.99 s, respectively. As evident from the results presented in Fig. 5(b), the curvature of the beam changes upon contact with the spring located at its tip.

As emphasized earlier, the match between the fractional- and the integer-order results is exact (i.e., their relative error is zero) at all times; however, both solutions are approximated due to the use of numerical approaches. The nonlinearity in the contact problem is now fully captured in the order of the VO-RL operator Keq which is evident from the switch in Keq (see Fig. 5(a)). A few closing comments on this problem. Although we have presented the results for the case of a concentrated load, the same formulation applies identically to distributed transverse loads. Also, the formulation can be extended to the study of beams whose finite stiffness contact can occur everywhere along the length and at multiple locations, simultaneously. This latter case could be seen as the counterpart of the more classical problem of a beam coming in contact with an elastic foundation.

### 4.3 Application to Nonlinear Hysteretic Behavior.

In this section, we extend the VO-RL formalism and illustrate how it can be used to model hysteretic behavior (associated with constitutive relations) in lumped parameter systems. Similar to the contact problem discussed above where the variable-order $β(t)$ was crafted to change the stiffness upon contact, VO operators can be crafted to capture the irreversible nonlinear evolution associated with hysteretic constitutive relations. Consider a lumped spring–mass–damper system (see Fig. 6(a)) such that the constitutive equation of the spring varies from linear to nonlinear irreversible response depending on the total elongation of the spring. Since the spring's behavior is assumed irreversible, the VO operator must be designed to capture the occurrence of a residual elongation in the spring. The equation of motion for the hysteretic behavior is generalized to the form

$Mξ¨+Cξ˙+[K−(0RLDtβ1(t)K̃(t)) (0RLDtβ2(t)K)]×[ξ−0RLDtβ3(t)γ(t)]=F(t)$
(22a)

where
$β1(t)=exp(−s1(|ξ(t)|−d0)︸s1(t))$
(22b)
$β2(t)=exp(−s2sgn(F˙(t))︸s2(t))$
(22c)
$β3(t)=exp(s3sgn(F˙(t))︸−s3(t))$
(22d)
$K̃(t)≜K̃(|F(t)|−FY)$
(22e)
$γ(t)=sup(ξ(t))−FsK$
(22f)
Fig. 6
Fig. 6
Close modal

In the above equations M, C, and K are the mass, damping, and linear elastic stiffness coefficients of the lumped system and $ξ(t)$ denotes the displacement of the oscillator. The system is subject to an external force F(t), and $Fs=F(ts)$ such that $ξ(ts)=sup(ξ(t))$. The coefficients $s1, s2$, and s3 (in the VO $β1(t), β2(t)$, and $β3(t)$) are scaling factors. The role of these scaling factors and the calculation approach is exactly identical to the one exemplified in previous paragraphs. The quantity d0 describes the linear response limit of the spring, and $FY=Kd0$ is the force applied on the spring at this limit. $K̃(t)$ is defined to be a function of $|F(t)|−FY$.

In order to understand the VO governing Eq. (22), the role of the different VO-RL operators must be discussed. The VO-RL operator $0RLDtβ1(t)K̃(t)$ accounts for the change in the stiffness of the system dictated by the order variation $β1(t)$. When the function $s1(t)=|ξ(t)|−d0$, which defines the order $β1(t)$, changes its sign (i.e., when $|ξ(t)|>d0$) the system enters the nonlinear regime. The VO-RL operator $0RLDtβ2(t)K$ accounts for the effect of the direction of loading on the stiffness via the function $s2(t)=sgn(F˙(t))$ in the order $β2(t)$. This means that the net stiffness of the system given by $K(t,ξ,F)=K−(0RLDtβ1(t)K̃(t)) (0RLDtβ2(t)K)$ fully captures both nonlinear and direction-dependent effects (i.e., irreversibility). Finally, $0RLDtβ3(t)γ(t)$ accounts for the residual elongation in the spring during unloading. Mathematically, it allows the spring to have a memory of the residual elongation via $β3(t)$. The VO-RL operators described above can be combined to capture the hysteretic behavior in the following fashion:
$Mξ¨+Cξ˙+Kξ︸I.Linear response−[0RLDtβ1(t)K̃(t) 0RLDtβ2(t)K]ξ︸II.Captures transition to nonlinear response−[K−(0RLDtβ1(t)K̃(t)) (0RLDtβ2(t)K)]×[0RLDtβ3(t)γ(t)]︸III.Captures residual elongation and unloading response=F(t)$
(23)

For clarity, we have marked these different regimes in Fig. 6(b) that presents the numerical results obtained to validate the VO formulation.

In order to verify the response of the system given in Eq. (22), we performed a numerical simulation in which the system was loaded quasi-statically with increasing quasi-static force F(t), and the released. The specific parameters used in the simulation were M =1 kg, C =1 N·s/m, K =1 N/m, $d0=5$ m, and $si=103$ m−1; i =1, 2, 3. We set $K̃(t)$ as
$K̃(t)=|F(t)|−54$
(24)
$[K−(0RLDtβ1(t)(|F0|−54)) (0RLDtβ2(t)K)]×[ξ0−0RLDtβ3(t)(supξ0−F0K)]=F0$
(25)

where ξ0 is the static displacement due to the quasi-static force $F(t)=F0$. We emphasize that the choice of $K̃(t)$ in this study is arbitrary but it can be easily connected to the specific properties of the material being modeled by the spring element.

In order to demonstrate the capability of the VO formulation, we identified three different cases and loaded the system quasi-statically using three different force amplitudes, each one followed by an unloading phase. In case (a), $0≤F0≤4$ N such that the spring exhibits linear response. In case (b), $0≤F0≤6$ N such that the linear limit of the spring is exceeded and the nonlinear irreversible behavior arises, hence giving rise to a residual elongation. Finally, in case (c), the system is loaded to a larger maximum force $0≤F0≤7$ N such that it accumulates a larger residual elongation. The simulations were obtained using Eq. (25) and the results are expressed in terms of the load–displacement curve in Fig. 6. Further, we have compared the results obtained from the VO equation of motion (denoted by $F0β$) with the corresponding integer-order equation of motion (denoted by F0) which incorporates the same hysteretic behavior. For the sake of completeness, we provide below the stiffness variation (in N/m) used to simulate the integer-order equation of motion
$K(t,ξ,F)={1 ξ(t)≤d0 ∩ F˙(t)≥09−|F(t)|4 ξ(t)>d0 ∩ F˙(t)≥01 F˙(t)<0$
(26)

As evident from Fig. 6, the match is exact. This result should not surprise because the fractional equation of motion is an exact reformulation of the integer-order equation of motion. We emphasize that all the characteristic of the hysteretic behavior have been captured via the VO operators given in Eq. (22). By visually inspecting the numerical results, it is evident that the spring responds in a linear fashion for $|F0| and in a nonlinear irreversible fashion for $|F0|>FY$. Note that the stiffness of the spring also depends on the loading direction as expected in a hysteretic system. Further, the spring remembers the residual elongation and the magnitude of the residual elongation depends on the maximum force exerted on the spring. All the above features are well known in hysteretic systems, and the numerical results show that the VO m odel capture this behavior very closely.

## 5 General Remarks

We emphasize that, despite using the same functional form for the governing equations, different laws for the order variations can result in largely different dynamics. This is mostly due to the impact of the order law on the memory of the fractional-order operator and the subsequent impact on the resulting dynamics of the system. We note that this is different from the order-memory of the different classes of variable-order operators discussed in Ref. [19]. In the hysteresis example, the VO operator retains a memory of the residual elongation during the unloading cycle. This is unlike the VO operator used in modeling the contact problem where the operator changes instantaneously depending on whether the contact has occurred or not. In this latter case, the change of the operator is instantaneous and does not have any memory of the previous states. Despite the simplicity of the discrete models presented in this study, it is immediate to envision how a similar modeling approach could be implemented for complex structural problems involving, for example, bistable dynamics, discontinuities, contacts, and variable stiffness joints.

## 6 Conclusions

This work explored the use of variable-order fractional operators for the analysis of selected discrete and continuous nonlinear dynamical systems. The most remarkable result consists of the observed evolutionary nature of the VO operator which allows a natural mathematical treatment of dynamical systems transitioning from linear to highly nonlinear behavior without requiring a reformulation of the underlying governing equations. In other terms, thanks to the ability to update the order of the differ-integral operator as a function of either dependent or independent variables, the governing equations can evolve seamlessly and in real-time from linear to nonlinear simply based on the instantaneous response of the system. In this context, it was discussed how an apparently unsettling property of the RL operator (that is the nonvanishing derivative of a constant when using a finite lower terminal of integration) has very useful implications to model systems with strong and evolving nonsmooth nonlinearities. We developed a physics-driven strategy to exploit this unique behavior of Riemann–Liouville derivatives and to define VO operators capable of capturing complex transitions in nonlinear dynamical systems. Specific applications to the modeling of either contact dynamics or hysteretic behavior were addressed. Results demonstrated that the VO operator is able to capture either the contact status or the transition from linear to hysteretic response. Furthermore, despite the equations can describe systems exhibiting a marked nonlinear dynamics, the operators remain linear. This latter aspect suggests that many operations and algorithms developed for systems of linear differential equations could still be applicable to this class of equations. We envision that the variable-order formulation will provide a solid foundation to build powerful computational tools for the analysis of both discrete and continuous nonlinear systems.

## Acknowledgment

The authors gratefully acknowledge the financial support of the National Science Foundation under the grants DCSD #1825837 and MOMS #1761423, and the Defense Advanced Research Project Agency under grant #D19AP00052. The content and information presented in this paper do not necessarily reflect the position or the policy of the government. The material is approved for public release; distribution is unlimited.

## Funding Data

• DARPA (Grant No. D19AP00052; Funder ID: 10.13039/100000185).

• National Science Foundation (Grant Nos. 1825837 and 1761423; Funder ID: 10.13039/100000001).

## References

1.
Mainardi
,
F.
,
1996
, “
Fractional Relaxation-Oscillation and Fractional Diffusion-Wave Phenomena
,”
Chaos, Solitons Fractals
,
7
(
9
), pp.
1461
1477
.10.1016/0960-0779(95)00125-5
2.
Fellah
,
Z. E. A.
,
Fellah
,
M.
,
Sebaa
,
N.
,
Lauriks
,
W.
, and
Depollier
,
C.
,
2006
, “
Measuring Flow Resistivity of Porous Materials at Low Frequencies Range Via Acoustic Transmitted Waves
,”
J. Acoust. Soc. Am.
,
119
(
4
), pp.
1926
1928
.10.1121/1.2179749
3.
Hollkamp
,
J. P.
,
Sen
,
M.
, and
Semperlotti
,
F.
,
2019
, “
Analysis of Dispersion and Propagation Properties in a Periodic Rod Using a Space-Fractional Wave Equation
,”
J. Sound Vib.
,
441
, pp.
204
220
.10.1016/j.jsv.2018.10.051
4.
Buonocore
,
S.
,
Sen
,
M.
, and
Semperlotti
,
F.
,
2019
, “
Occurrence of Anomalous Diffusion and Non-Local Response in Highly-Scattering Acoustic Periodic Media
,”
New J. Phys.
,
21
(
3
), p.
033011
.10.1088/1367-2630/aafb7d
5.
Bagley
,
R. L.
, 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
6.
Bagley
,
R. L.
, and
Torvik
,
P. J.
,
1985
, “
Fractional Calculus in the Transient Analysis of Viscoelastically Damped Structures
,”
AIAA J.
,
23
(
6
), pp.
918
925
.10.2514/3.9007
7.
Chatterjee
,
A.
,
2005
, “
Statistical Origins of Fractional Derivatives in Viscoelasticity
,”
J. Sound Vib.
,
284
(
3–5
), pp.
1239
1245
.10.1016/j.jsv.2004.09.019
8.
Carpinteri
,
A.
,
Cornetti
,
P.
, and
Sapora
,
A.
,
2014
, “
Nonlocal Elasticity: An Approach Based on Fractional Calculus
,”
Meccanica
,
49
(
11
), pp.
2551
2569
.10.1007/s11012-014-0044-5
9.
Sumelka
,
W.
, and
Blaszczyk
,
T.
,
2014
, “
Fractional Continua for Linear Elasticity
,”
Arch. Mech.
,
66
(
3
), pp.
147
172
.https://am.ippt.pan.pl/index.php/am/article/view/v66p147
10.
Patnaik
,
S.
,
Sidhardh
,
S.
, and
Semperlotti
,
F.
,
2020
, “
A Ritz-Based Finite Element Method for a Fractional-Order Boundary Value Problem of Nonlocal Elasticity
,” arXiv preprint
arXiv:2001.06885
.https://arxiv.org/abs/2001.06885
11.
Patnaik
,
S.
,
Sidhardh
,
S.
, and
Semperlotti
,
F.
,
2020
, “
Fractional-Order Models for the Static and Dynamic Analysis of Nonlocal Plates
,” arXiv preprint
arXiv:2002.10244
.https://arxiv.org/abs/2002.10244
12.
Hollkamp
,
J. P.
,
Sen
,
M.
, and
Semperlotti
,
F.
,
2018
, “
Model-Order Reduction of Lumped Parameter Systems Via Fractional Calculus
,”
J. Sound Vib.
,
419
, pp.
526
543
.10.1016/j.jsv.2018.01.011
13.
Hollkamp
,
J. P.
, and
Semperlotti
,
F.
,
2020
, “
Application of Fractional Order Operators to the Simulation of Ducts With Acoustic Black Hole Terminations
,”
J. Sound Vib.
,
465
, p.
115035
.10.1016/j.jsv.2019.115035
14.
Podlubny
,
I.
,
1998
,
Fractional Differential Equations: An Introduction to Fractional Derivatives, Fractional Differential Equations, to Methods of Their Solution and Some of Their Applications
, Vol.
198
,
Elsevier
, San Diego, CA.
15.
Richard
,
H.
,
2014
,
Fractional Calculus: An Introduction for Physicists
,
World Scientific
,
Singapore
.
16.
Glöckle
,
W. G.
, and
Nonnenmacher
,
T. F.
,
1995
, “
A Fractional Calculus Approach to Self-Similar Protein Dynamics
,”
Biophys. J.
,
68
(
1
), pp.
46
53
.10.1016/S0006-3495(95)80157-8
17.
Samko
,
S. G.
, and
Ross
,
B.
,
1993
, “
Integration and Differentiation to a Variable Fractional Order
,”
Integr. Transforms Spec. Funct.
,
1
(
4
), pp.
277
300
.10.1080/10652469308819027
18.
Samko
,
S. G.
,
1995
, “
Fractional Integration and Differentiation of Variable Order
,”
Anal. Math.
,
21
(
3
), pp.
213
236
.10.1007/BF01911126
19.
Lorenzo
,
C. F.
, and
Hartley
,
T. T.
,
2002
, “
Variable Order and Distributed Order Fractional Operators
,”
Nonlinear Dyn.
,
29
(
1/4
), pp.
57
98
.10.1023/A:1016586905654
20.
Sun
,
H.
,
Chen
,
W.
, and
Chen
,
Y.
,
2009
, “
Variable-Order Fractional Differential Operators in Anomalous Diffusion Modeling
,”
Phys. A
,
388
(
21
), pp.
4586
4592
.10.1016/j.physa.2009.07.024
21.
Chen
,
W.
,
Zhang
,
J.
, and
Zhang
,
J.
,
2013
, “
A Variable-Order Time-Fractional Derivative Model for Chloride Ions Sub-Diffusion in Concrete Structures
,”
Fractional Calculus Appl. Anal.
,
16
(
1
), pp.
76
92
.10.2478/s13540-013-0006-y
22.
Coimbra
,
C. F.
,
2003
, “
Mechanics With Variable-Order Differential Operators
,”
Ann. Phys.
,
12
(
1112
), pp.
692
703
.10.1002/andp.200310032
23.
Diaz
,
G.
, and
Coimbra
,
C.
,
2009
, “
Nonlinear Dynamics and Control of a Variable Order Oscillator With Application to the Van Der Pol Equation
,”
Nonlinear Dyn.
,
56
(
1–2
), pp.
145
157
.10.1007/s11071-008-9385-8
24.
Patnaik
,
S.
,
Hollkamp
,
J. P.
, and
Semperlotti
,
F.
,
2020
, “
Applications of Variable-Order Fractional Operators: A Review
,”
Proc. R. Soc. A
,
476
(
2234
), p.
20190498
.10.1098/rspa.2019.0498
25.
Lorenzo
,
C. F.
, and
Hartley
,
T. T.
,
2007
, “
Initialization, Conceptualization, and Application in the Generalized (Fractional) Calculus
,”
Crit. Biomed. Eng.
,
35
(
6
), pp.
447
553
.10.1615/CritRevBiomedEng.v35.i6.10
26.
Chunduru
,
S. P.
,
Kim
,
M. J.
, and
Mirman
,
C.
,
2011
, “
Failure Analysis of Railroad Couplers of Aar Type e
,”
Eng. Failure Anal.
,
18
(
1
), pp.
374
385
.10.1016/j.engfailanal.2010.09.016
27.
Wu
,
Q.
,
Luo
,
S.
,
Xu
,
Z.
, and
Ma
,
W.
,
2013
, “
Coupler Jackknifing and Derailments of Locomotives on Tangent Track
,”
Veh. Syst. Dyn.
,
51
(
11
), pp.
1784
1800
.10.1080/00423114.2013.830184
28.
Wilson
,
J. F.
, and
Callis
,
E. G.
,
2004
, “
The Dynamics of Loosely Jointed Structures
,”
Int. J. Non-Linear Mech.
,
39
(
3
), pp.
503
514
.10.1016/S0020-7462(02)00219-6
29.
Semperlotti
,
F.
, and
Conlon
,
S. C.
,
2010
, “
Structural Damage Identification in Plates Via Nonlinear Structural Intensity Maps
,”
J. Acoust. Soc. Am.
,
127
(
2
), pp.
EL48
EL53
.10.1121/1.3290175
30.
Semperlotti
,
F.
,
Wang
,
K.
, and
Smith
,
E.
,
2009
, “
Localization of a Breathing Crack Using Super-Harmonic Signals Due to System Nonlinearity
,”
AIAA J.
,
47
(
9
), pp.
2076
–20
86
.10.2514/1.38947
31.
Kim
,
G.-W.
,
Johnson
,
D. R.
,
Semperlotti
,
F.
, and
Wang
,
K. W.
,
2011
, “
Localization of Breathing Cracks Using Combination Tone Nonlinear Response
,”
Smart Mater. Struct.
,
20
(
5
), p.
055014
.10.1088/0964-1726/20/5/055014
32.
Chondros
,
T. G.
,
Dimarogonas
,
A. D.
, and
Yao
,
J.
,
2001
, “
Vibration of a Beam With a Breathing Crack
,”
J. Sound Vib.
,
239
(
1
), pp.
57
67
.10.1006/jsvi.2000.3156
33.
Zhang
,
M.
,
Shen
,
Y.
,
Xiao
,
L.
, and
Qu
,
W.
,
2017
, “
Application of Subharmonic Resonance for the Detection of Bolted Joint Looseness
,”
Nonlinear Dyn.
,
88
(
3
), pp.
1643
1653
.10.1007/s11071-017-3336-1
34.
Semperlotti
,
F.
, and
Conlon
,
S. C.
,
2010
, “
Nonlinear Structural Surface Intensity: An Application of Contact Acoustic Nonlinearity to Power Flow Based Damage Detection
,”
Appl. Phys. Lett.
,
97
(
14
), p.
141911
.10.1063/1.3491801
35.
Lamberti
,
A.
, and
Semperlotti
,
F.
,
2013
, “
Detecting Closing Delaminations in Laminated Composite Plates Using Nonlinear Structural Intensity and Time Reversal Mirrors
,”
Smart Mater. Struct.
,
22
(
12
), p.
125006
.10.1088/0964-1726/22/12/125006
36.
Srinivas
,
K.
,
Rao
,
C. L.
, and
Parameswaran
,
V.
,
2019
, “
Experimental Investigation of Rate Sensitive Mechanical Response of Pure Polyurea
,”
Dynamic Behavior of Materials
, Vol.
1
,
Springer
,
Cham, Switzerland
, pp.
173
179
.
37.
Zhang
,
J.
,
Lin
,
Z.
,
Wong
,
A.
,
Kikuchi
,
N.
,
Li
,
V.
,
Yee
,
A.
, and
Nusholtz
,
G.
,
1997
, “
Constitutive Modeling and Material Characterization of Polymeric Foams
,”
ASME J. Eng. Mater. Technol.
,
119
(
3
), pp.
284
291
.10.1115/1.2812258
38.
Hopkinson
,
B.
, and
Williams
,
T. G.
,
1912
, “
The Elastic Hysteresis of Steel
,”
Proc. R. Soc. London. Ser. A
,
87
(
598
), pp.
502
511
.10.1098/rspa.1912.0104
39.
Baber
,
T. T.
, and
Noori
,
M. N.
,
1986
, “
Modeling General Hysteresis Behavior and Random Vibration Application
,”
ASME J. Vib., Acoust., Stress, Reliab. Des.
,
108
(
4
), pp.
411
420
.10.1115/1.3269364
40.
Bhattacharjee
,
A.
, and
Chatterjee
,
A.
,
2013
, “
Dissipation in the Bouc–Wen Model: Small Amplitude, Large Amplitude and Two-Frequency Forcing
,”
J. Sound Vib.
,
332
(
7
), pp.
1807
1819
.10.1016/j.jsv.2012.10.026
41.
Moon
,
F. C.
, and
Kalmár-Nagy
,
T.
,
2001
, “
Nonlinear Models for Complex Dynamics in Cutting Materials
,”
Philos. Trans. R. Soc. London. Ser. A
,
359
(
1781
), pp.
695
711
.10.1098/rsta.2000.0751
42.
Li
,
C.
, and
Zeng
,
F.
,
2015
,
Numerical Methods for Fractional Calculus
,
Chapman and Hall/CRC
,
London
.
43.
Reddy
,
J. N.
,
1993
,
An Introduction to the Finite Element Method
,
McGraw-Hill
,
New York
.