Modeling Contacts and Hysteretic Behavior in Discrete Systems Via Variable-Order Fractional Operators

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 [1–4], formulation of constitutive equations for viscoelastic materials [5–7], nonlocal elasticity [8–11], 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. [17–19] 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).

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.

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 a₀) 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 a₀ 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 a₀ to a [25 ]. Further details on the initialization procedure and its importance can be found inRef. [25 ].

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.

where $m−1≤β0<m$⁠.

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 s₀ for the different physical problems is reported in each specific section.

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

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.

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.

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 [30–32], 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 [36–38]. 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.

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 (K₀) and damper (C₀) are engaged only when $ξ1(t)−ξ2(t)>d$⁠.

when $ξ1−ξ2>d$⁠.

The scaling factor s₀ 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 s₀ 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 s₀.

From a more practical perspective, the numerical value of the parameter s₀ 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 P₀ is a user-defined percentage value. At this point, it is easy to see that if t_k 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<s(t̃k)≤P0d$⁠. 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 K_eq and C_eq as VO-RL operators. Using Eqs. (8) and (11), it can be shown that K_eq and C_eq are

where K_eq and C_eq 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.

where $ϕ1$ (and $ϕ2$⁠) is either K₀ (and $K0′)$ or C₀ (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 s₀ 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 K_eq).

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 K_eq which clearly switches between the values 0 and K₀ 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.

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

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. K₀ is the stiffness of the spring located at the contact point.

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

where the VO-RL term K_eq 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.

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·m², ρA = 1 kg/m, L = 1 m, $K0=1$ N/m, $Δ=0.03$ m. The scaling parameter s₀ 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.

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 K_eq which is evident from the switch in K_eq (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.

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

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 s₃ (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 d₀ 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$⁠.

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

where ξ₀ 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.

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

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.

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.

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.

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

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