A Variable-Order Fractional Constitutive Model to Characterize the Rate-Dependent Mechanical Behavior of Soft Materials

: Building an accurate constitutive model for soft materials is essential for better under-standing its rate-dependent deformation characteristics and improving the design of soft material devices. To establish a concise constitutive model with few parameters and clear physical meaning, a variable-order fractional model is proposed to accurately describe and predict the rate-dependent mechanical behavior of soft materials. In this work, the discrete variable-order fractional operator enables the predicted stress response to be entirely consistent with the whole stress history and the fractional order’s path-dependent values. The proposed model is further implemented in a numerical form and applied to predict several typical soft materials’ tensile and compressive deformation behavior. Our research indicates that the proposed variable-order fractional constitutive model is capable of predicting the nonlinear rate-dependent mechanical behavior of soft materials with high accuracy and more convinced reliability in comparison with the existing fractional models, where the fractional order contains a constant initial order to depict the initial elastic response and a linear variable-order function to account for the strain hardening behavior after acrossing the yield point.


Introduction
Soft material covers a wide range of materials, including polymers, hydrogels, elastomers et al., it has broad prospects of application in emerging industries and technologies such as soft display [1], flexible microelectronics [2] and soft robots [3] due to their excellent stretchability, lightweight and high energy density.Practically, the rate-dependent viscoelastic behavior will significantly affect their performance in service [4].For instance, due to the vast range of technical applications in the aerospace and automotive sectors, as well as in the field of fluid transportation, Animasaun [5] confirmed that the viscoelastic characteristics of polymeric soft materials have received extensive study.In this regard, characterizing the viscoelastic response of soft materials is of great importance for their applications.It is crucial to develop accurate constitutive models to describe their mechanical response to design soft materials with desired properties.
For most soft materials, their internal structures have microscopic polyphaser characteristics.Hepburn [6] reported that the macromolecular long chains of polyurethane elastomer soft materials are intertwined and distributed roundaboutly.The long chains of molecules are divided into two regions: the hard chain segment and the soft one.Under external load, the release of the wrapped soft chain in the deformation process will cause a change in the density of the hard chain and the soft chain area, leading to stress softening, and the soft material is easy to produce macroscopic deformation [7,8].Furthermore, the movement of the molecular chain inside the rubber soft material will also be hindered and affected by the surrounding molecules.Its mechanical properties in the deformation process are a complex behavior of the mixture of elasticity and viscosity, which makes it exhibit rate-dependent property in the process of loading and unloading, also known as time correlation.Given the complexity of the internal structure of soft materials, the deformation of soft materials is entirely different from that of traditional linear elastic materials, which will produce complex and variable mechanical responses under different environmental conditions.This challenges the establishment of constitutive models for soft materials [9,10].
As described above, the stress-strain behavior of soft materials exhibits obvious rate dependency or is a typical time-dependent mechanical behavior, and its stress response is not only related to the strain but also the strain rate of deformation.Qu et al. [11] studied the rate-dependent tensile mechanical response of VHB947 soft material under pure shear deformation mode.The results showed that the acrylic elastomeric soft material had obvious tensile rate dependency.Wissler et al. [12] used the Prony series theory to model the time-dependent viscoelastic behavior of VHB4910 elastomeric soft materials.Foo et al. [13] and Zhao et al. [14] used a standard linear solid rheological model satisfying the thermodynamic consistency to describe the rate-dependent mechanical behavior of elastomeric soft materials.Their model could well capture the rate-dependent behavior of soft materials.However, the internal variables are introduced in the modeling process, where the multiplicative decomposition of the deformation tensor is adopted to modify the traditional linear superposition process.This treatment significantly increases the complexity of the proposed constitutive model and the number of governing equations.There still lacks a proper constitutive model of soft materials with concise form and few model parameters.
In the last decades, fractional calculus has convinced to be a powerful mathematical approach to modeling the time-dependent mechanical behavior of viscoelastic material.The superiority lies in the remarkable reduction in model parameters required to develop a constitutive model involving fractional operators.Such models with the constant fractional order are among the most fashionable ones in the study of fractional viscoelasticity, including soils [15], polymers [16,17], and construction materials [18].However, the rate-dependent mechanical behavior, as well as its significant nonlinearities under extensive deformation regimes derived from those soft materials like rubbers [19,20], hydrogels [21-23] and dielectric elastomers [24,25], both featuring memory effects whose behavior varies over time, which exhibits highly complex material properties and poses a significant challenge with regards to constitutive modeling of them.Recently, the original point of building fractional operators by enabling the order of fractional calculus to be a function of the state variable, i.e., space or time, has increasingly arisen widespread attention.Such operators involving variable fractional orders are proved to be a promising tool to depict the complicated physical phenomena relating to memory effects.Numerous operators with variable fractional order have been suggested in the literature like physics [26], anomalous diffusion [27][28][29], signal process [30], solid mechanics [31] and transport process in inhomogeneous media [32], in which the fractional viscoelasticity is one of the most fashionable application domains.Meng et al. [16] constructed a fractional operator with variable order to model the nonlinear viscoelastic response of polymers associated with the evolution of material properties.Xiang et al. [4] investigated the creep behavior of natural fiber polymer composites (NFPCs) material based on variable-order fractional operators.Additionally, extensive rheological studies of viscoelastic materials are reported by [33][34][35].In these variable-order fractional definitions mentioned above, a typically handling variable order extension of classical fractional derivatives is acquired by substituting an order function α(t) for the constant fractional order α, and the performance of these models was obtained by simply fitting it to the experimental data.This treatment of the fractional constitutive model with variable-order makes the description of its physical process unclear.
As Garrappa et al. [36] reported, the mathematical process of these operators obtained by simply replacing the constant order is problematic.Besides, the general solution of a variable-order fractional operator cannot be formulated by directly substituting a given order function for the constant fractional order.Di Paola et al. [37] pointed out that the existing modeling process is not derived from a robust mathematical framework and cannot be consistent with the classical Boltzmann principle, given viscoelasticity.To establish a simple constitutive model with clear physical meaning to accurately describe the rate-dependent mechanical behavior of soft materials, in this paper, with the help of the variable-order fractional calculus in describing the time-dependent mechanical behavior of soft materials, we develop a nonlinear fractional constitutive model to predict the rate-dependent viscoelastic response of soft materials, the fractional operator with variable-order is numerically approximated by the Grünwald-Letnikov derivative and can be easily applied to construct the lower triangular strip matrix, this treatment dramatically simplifies the operation of the numerical algorithm.
The structure of this work is arranged as follows.The preliminary knowledge of fractional derivatives is introduced in Section 2, and Section 3 gives the scheme of the numerical discretization method based on the Grünwald-Letnikov derivative.The nonlinear variable-order fractional model is established and solved in Section 4, and the prediction results of the numerical calculation of the proposed model.Finally, the main conclusions are arranged in Section 5.

Preliminary Knowledge of Fractional Calculus
Generally, the fractional differential and integral operators are defined based on a natural extension of the order of integer operators.Herein, three commonly accepted fractional definitions are recalled, namely the Grünwald-Letnikov (GL), Riemann-Liouville (RL), and Caputo definitions [38][39][40].Firstly, the fractional integral D −α t of the function g(t) is defined as: where Γ(z) = ∞ 0 t z−1 e −t dt is the Euler-Gamma function.The Caputo derivative of the function f (t) with fractional order α > 0 satisfies where m − 1 < α < m, m ∈ Z + .By interchanging the processes of differentiation and integration in Equation ( 2), the so-called RL derivative of the function f (t) with order α > 0 can be obtained as follows Secondly, the GL derivative of the function f (t) with fractional order α > 0 reads where m − 1 < α < m, j = (t − t 0 )/h and α j denotes Newton's binomial coefficients with the form of If one considers a class of functions f (t) having m + 1 continuous derivatives for t ≥ 0, the RL definition is equivalent to the GL one, yielding While the Caputo and RL derivatives are not equivalent unless the initial condition satisfies

Numerical Discretization Method Based on the GL Derivative
From a different point of view, considering the classical backward finite differences scheme, the first derivative of the function g(t) can be formulated as [41,42] where t k is the generic time and k is the steps of discretized variable range, h is the temporal step between two subsequent values of variable t, i.e., h = t k − t k−1 .Then, the N equations expressed by Equation ( 8) can be grouped in the following matrix, giving where B N is a lower triangular strip matrix.By generalizing Equation ( 9), one can obtain the j-th derivative with the form of [43,44] N B (1) Correspondingly, the j-th integral proves to be governed by the inverse operator, yielding where the matrix B (j) N and its inverse form B can be built by the following simple rules: By substituting the derivative order with a non-integer value α > 0 in Equation ( 12), the GL fractional derivative and integral can be written as where ∇ α and ∇ −α denotes the GL derivative and integral, respectively, and N are the lower triangular strip matrix related to the fractional differential and integral operator, which the following matrix can evaluate: in which the element ω (α) j can be recursively constructed as following: ω It can be observed that to evaluate the fractional derivative of a function f (t) at a specific point t n , all the precedent values of the function have to be considered.Besides, we have where β (α) j can be obtained by substituting (−α) into ω (α) j .As indicated above, the GL derivative may be the correct way of studying the fractional derivative since all the other definitions can be deduced from it.In addition, it might be the most suitable approach for numerical computations due to its recursive coefficients and discretization scheme.This will be utilized to calculate the following nonlinear variableorder fractional constitutive model.

Model Development
For the sake of interpretation, the spring-pot is considered a reference to model the stress-stretch response of soft materials due to its mechanical property naturally transforming from the pure solid-like behavior to that of ideal fluid as the fractional order varies from 0 to 1.The variable-order fractional constitutive relation of the spring-pot is defined as [16,[45][46][47] where is the stress, ε is the strain, and E(MPa) is the material parameter considered as Young's modulus, θ(s) rep-resents the viscous parameter with the form of η/E, commonly defined as the relaxation time and η is the viscosity coefficient, D α(t) t represents the Caputo variable-order fractional derivative.From the numerical perspective, we discrete Equation ( 16) in time using the uniform grids with 0 < α(t) < 1, that is Meanwhile, the order of fractional derivative in Equation ( 16) is generally defined as a piecewise function of discretized time increments, then the Caputo derivative of strain can be approximated by the following relation [38,44] where c α n k are binomial coefficients defined with variable fractional order.Compared to the fractional operator with constant order, the elements of the lower triangular strip matrix B α n n−1 are extended as (also refer to [37,48]) Accordingly, its counterpart with that of fractional integral A α n n−1 can be formulated by the ensuing relation Follows from Equation ( 17), the stress response described by the discrete fractional model can be formulated as Considering the following formula For the sake of simplicity, the above stress response can be expressed under the zero initial condition From Equation (22), it is observed that the formula underlines the property of historical dependence and non-locality of fractional calculus, namely, the stress value σ n at t = t n is dependent on all the previous stress history of t = 0, 1, • • • , n − 1 and the strain values ε n at t n .Furthermore, this treatment dramatically simplifies the calculation of the fractional model, which benefits from the recursion formula of matrix coefficients.In the following work, the performance of the numerical stress solution will be testified by adopting several material tests.

Validity of the Proposed Model under Tensile Tests 4.2.1. Results for Nitrile Rubber
In this subsection, to testify to the applicability of the discrete fractional model, a series of uniaxial tensile tests were conducted on the nitrile rubber at stretch rates of 20, 50, 200, and 500/min using an electronic universal-testing machine (WDW-3050, Changchun Testing Machine Research Institute, Changchun, China) in our lab.All tests were performed in a temperature-controlled room (70 • C).The samples were kept in a thermostat for at least 30 min to eliminate the residual stress and thermal history before each set of subsequent experiments.The obtained tensile behaviors of described rubber samples are shown in Figure 1 (empty blue circles).It can be observed that the four stress-strain curves at different stretch rates all share a similar trend: the stress grows with the increasing applied strain, showing a convex shape with intervals of strain-softening, and it also increases distinctly with increasing stretch rate for a given strain, this implies that the described rubber material exhibits a typical rate-dependent behavior.
ing applied strain, showing a convex shape with intervals of strain-softening, and it also increases distinctly with increasing stretch rate for a given strain, this implies that the described rubber material exhibits a typical rate-dependent behavior.
To determine the optimal values of unknown model parameters E ,  , and n  , an optimization method is used by minimizing the residual at the discretized points goodness.Specifically, Moghaddam et al. [49] and Manish et al. [50] minimized  represent the discretized stress of the experimental point and calcu- lated the numerical stress solution of the proposed model.The Newton backward difference algorithm performed all the calculation steps in MATLAB (Version: R2016a).For the sake of simplicity, herein, the fractional order for each simulation is set to a constant considering the uniform mechanical properties reflected by the stress-strain curves with 1200 j = . Table 1 shows the derived optimal material parameters of the predicted stress response for the uniaxial tensile regime at four different stretch rates.Figure 1 shows the calibrated stress response against the experimental stress data of described rubber materials at four stretch rates.Based on the experimental and numerical results of the discrete variable fractional order model, it is inferred that the proposed model gives the perfect description of the rate-dependent mechanical responses of the described rubber material.Furthermore, Figure 2 displays the variance of obtained fractional orders plotted as a function of the stretch rate.It is observed that the order decreases with the increasing stretch rate.As an essential parameter in characterizing the material To determine the optimal values of unknown model parameters E, θ, and α n , an optimization method is used by minimizing the residual at the discretized points goodness.Specifically, Moghaddam et al. [49] and Manish et al. [50] minimized 2 where σ i and σ i represent the discretized stress of the experimental point and calculated the numerical stress solution of the proposed model.The Newton backward difference algorithm performed all the calculation steps in MATLAB (Version: R2016a).For the sake of simplicity, herein, the fractional order for each simulation is set to a constant considering the uniform mechanical properties reflected by the stress-strain curves with j = 1200.Table 1 shows the derived optimal material parameters of the predicted stress response for the uniaxial tensile regime at four different stretch rates.Figure 1 shows the calibrated stress response against the experimental stress data of described rubber materials at four stretch rates.Based on the experimental and numerical results of the discrete variable fractional order model, it is inferred that the proposed model gives the perfect description of the rate-dependent mechanical responses of the described rubber material.Furthermore, Figure 2 displays the variance of obtained fractional orders plotted as a function of the stretch rate.It is observed that the order decreases with the increasing stretch rate.As an essential parameter in characterizing the material behavior during the process of fractional modeling, it is of great importance to illustrate the physical significance of the fractional order.For the sake of brevity, the discussion by introducing the equivalent viscoelasticity between the time-varying viscosity model and the fractional viscoelastic model in terms of the physical meaning of the fractional order is available in the Appendix A. The results indicate that the fractional order α is commonly identified as an index of material properties, which changes continuously from Hooke's elasticity to Newtonian viscosity as α varies from zero to one.The described soft material thereby depicts the viscoelastic solid-like behavior during stretching, and the observed descending trend of order variation shown in Figure 2 vividly reflects that the rubber material exhibits less viscous behavior at a high loading rate and, conversely, becomes softer at a low stretch rate.Figure 1 shows the calibrated stress response against the experimental stress da described rubber materials at four stretch rates.Based on the experimental and nume results of the discrete variable fractional order model, it is inferred that the prop model gives the perfect description of the rate-dependent mechanical responses o described rubber material.Furthermore, Figure 2 displays the variance of obtained tional orders plotted as a function of the stretch rate.It is observed that the order decre with the increasing stretch rate.As an essential parameter in characterizing the mat behavior during the process of fractional modeling, it is of great importance to illus the physical significance of the fractional order.For the sake of brevity, the discussio introducing the equivalent viscoelasticity between the time-varying viscosity model the fractional viscoelastic model in terms of the physical meaning of the fractional o is available in the Appendix A. The results indicate that the fractional order  is c monly identified as an index of material properties, which changes continuously Hooke's elasticity to Newtonian viscosity as  varies from zero to one.The descr soft material thereby depicts the viscoelastic solid-like behavior during stretching, and observed descending trend of order variation shown in Figure 2 vividly reflects tha rubber material exhibits less viscous behavior at a high loading rate and, conversely comes softer at a low stretch rate.

Results for Soft Elastomer
The mechanical response of soft elastomer is conventionally studied in uniaxial deformation at constant stretch rates.This typical soft material can generate a fast and large deformation in response to external environmental factors.In this subsection, to further verify the performance of the proposed discrete fractional model with variable order, the uniaxial tensile test data of VHB 4910 soft elastomer [51] is adopted for calculation.
From the outset, a representative set of uniaxial tensile experimental data of VHB 4910 at the stretch rate of 0.065/s is employed to exemplify the application of the discrete variable fractional order model in characterizing the deformation process under tensile tests.As shown in Figure 3, the stress-stretch curve (empty blue circle) exhibits a typical sigmoidal diagram with intervals of strain hardening.This stark distinction from the previous observations in Figure 1 demonstrates a complex material property of mixed elasticity and viscosity.As seen in Figure 3, the initial stress response first goes through an elastic region.This stage can be described by a constant fractional order α 0 , as reported in Section 4.2.1, followed by a sudden increase during stretching.This implies that a mutation of material properties appears before and after crossing the yield point, and the material behavior exhibits more nonlinear.In Figure 4, it can be observed that the strain hardening behavior can be well depicted by a variable order function rather α 1 (λ) than a constant order, which is defined as a linear trend to decrease according to our previous work [45][46][47], namely, α 1 (λ) = −mλ + β 0 where β 0 is the initial value of variable fractional order and m denotes the rate of variable order change.As shown in Figure 3, it can be observed that the proposed fractional model provides an excellent description of the uniaxial tensile stress response of described soft material.All the calculated material parameters are listed in Table 2.The above analysis indicates that the nonlinear mechanical response of the described elastomer, including the initial elastic stress response and the followed strain hardening behavior, was well-controlled by the two-stage fractional order variation comprising a fixed order α 0 and variable orders α 1 (λ).To illustrate the impact of the fractional order, a parameter study is performed in terms of the fixed order α 0 and the initial order β 0 .As shown in Figure 5, the variance of the fractional order has a significant influence on the stress response of soft elastomer: increasing α 0 will result in significant softening effects, making the elastomer exhibit more viscous behavior and, conversely, the strain hardening becomes more pronounced as the material property changes with decreasing the initial order β 0 .According to the physical significance of the fractional order discussed in the Appendix A and the basic theory of fractional viscoelasticity, the "hardening" behavior of the described soft material can be characterized by the decreasing variable orders, which corresponds to the relaxation of polymer chains during the deformation process, namely, the faster-growing stress response will induce a dramatically declining of the related variable fractional orders.To further demonstrate the applicability of the discrete fractional model to the tensile mechanical response of soft materials, the uniaxial stress-stretch experimental data [51] of soft elastomers at diverse loading rates are predicted.As described above, for each data set of tensile stress-stretch response, the modeling approach would be illustrated by the variance of the fractional order comprising an initial constant and a consequent variable.For the sake of comparison, the existing two fractional modeling methods are adopted to simulate the described stress response, where the order function of the former is piecewise-constant (CO model) [52] and that of the latter is assumed as piecewise-linear (VO model) [16,53].The simulation results of these two fractional models and the prediction of our discrete fractional model are plotted in Figure 6a.All the fitted model parameters are listed in Tables 3 and 4 for the CO and VO models, respectively.The calculated order curves of our model are displayed in Figure 6b.
For the sake of comparison, the existing two fractional modeling methods are adopted simulate the described stress response, where the order function of the former is pie wise-constant (CO model) [52] and that of the latter is assumed as piecewise-linear (V model) [16,53].The simulation results of these two fractional models and the predicti of our discrete fractional model are plotted in Figure 6a.All the fitted model paramete are listed in Tables 3 and 4 for the CO and VO models, respectively.The calculated ord curves of our model are displayed in Figure 6b.As shown in Figure 6a, the solid red line represents the predicted results of the proposed variable-order fractional constitutive model.The blue solids line and the black dash line denote the simulated results obtained by the CO and VO fractional models, respectively.As can be seen, the proposed variable-order fractional model provides an excellent description of the nonlinear rate-dependent mechanical behavior of the soft elastomers.In contrast, the simulation result of the CO fractional model exhibits large deviations from the experimental data, especially during the subsequent development stage of applied stretches.Otherwise, although the model with piecewise-linear orders (VO) can describe the stress response with acceptable accuracy, it is crucial to emphasize that the existing approach of the VO model is merely appreciable by the simplest fitting to the experimental data.And by contrast, the proposed variable-order fractional constitutive model is numerically discretized by the GL derivative, where the recursive matrix coefficients approximately substitute the fractional operator with variable order.This discretization scheme dramatically simplifies the operation of fractional calculus and makes the whole characteristics of the deformation process described by the proposed model much clear.As illustrated, all the precedent values of the function until the lower derivative limit are considered in numerical form.Namely, the predicted stress responses σ n at t = t n is entirely consistent with the whole stress history until t = t n−1 and with all the path-dependent values of variable order α k .Furthermore, the calculated evolutions of the variable fractional order over increasing stretches are plotted in Figure 6b.To characterize the deformation behaviors of the described soft elastomer, the fractional order includes two portions: an initial constant order for depicting the initial elastic response, and the declining order is utilized to account for the strain hardening behavior after crossing the yield point.The above analysis indicates that the proposed discrete fractional model is more suitable for depicting the nonlinear stress response of soft materials.

Validity of the Proposed Model under Compression
In this subsection, the proposed discrete variable fractional order model is applied to describe the compressive stress-strain behavior of aluminum foams with different porosities.The loading rate is 5 mm/min, as reported by [54].Our model's predicted mechanical responses under compression are shown in Figure 7a, with the calculated parameters listed in Table 5.
As displayed in Figure 7a, the compressive stress-strain curve of closed aluminum foams exhibits three obvious stages: linear deformation stage, plastic plateau stage, and densification stage.Furthermore, the yield stress increases at a faster rate with the decrease in porosity.The proposed variable-order fractional constitutive model is further used to describe the nonlinear mechanical behavior of aluminum foam soft materials under large plastic deformation.The simulated stress response is denoted by the solid red lines, as shown in Figure 7a.Based on the experimental and modeling observations in Figure 7a, our model's predictions agree with the test data of described aluminum foams.Moreover, as displayed in Figure 7b, the curves of fractional order are divided into two similar regions: a fixed order α 0 and a variable order function α 1 (λ) = −mλ + β 0 to consecutively illustrate the compressive mechanical behavior of described aluminum foams, where the initial purely elastic response is reflected by the constant order close to zero.The linear variable order trend subsequently depicts the densification phenomenon under compression to decrease.Additionally, the mutation of fractional order at λ 0 suggests that the mechanical property is mutated and departs from the linear elasticity at the yield point.Furthermore, as shown in Figure 7b, the fractional order increases along with the decrease in material porosities, suggesting that the aluminum foam exhibits more viscous behavior at lower porosity.This observation is also in line with the explanation of fractional order in characterizing the softness of described soft materials by the theory of fractional mechanics.
Fractal Fract.2022, 6, x FOR PEER REVIEW 12 of As shown in Figure 6a, the solid red line represents the predicted results of the p posed variable-order fractional constitutive model.The blue solids line and the black da line denote the simulated results obtained by the CO and VO fractional models, resp tively.As can be seen, the proposed variable-order fractional model provides an excell description of the nonlinear rate-dependent mechanical behavior of the soft elastome In contrast, the simulation result of the CO fractional model exhibits large deviations fr the experimental data, especially during the subsequent development stage of appl stretches.Otherwise, although the model with piecewise-linear orders (VO) can descr the stress response with acceptable accuracy, it is crucial to emphasize that the exist approach of the VO model is merely appreciable by the simplest fitting to the expe mental data.And by contrast, the proposed variable-order fractional constitutive mo is numerically discretized by the GL derivative, where the recursive matrix coefficie approximately substitute the fractional operator with variable order.This discretizat scheme dramatically simplifies the operation of fractional calculus and makes the wh characteristics of the deformation process described by the proposed model much cle As illustrated, all the precedent values of the function until the lower derivative limit considered in numerical form.Namely, the predicted stress responses n  at n tt = is tirely consistent with the whole stress history until and with all the path-depen ent values of variable order k  .
Furthermore, the calculated evolutions of the variable fractional order over incre ing stretches are plotted in Figure 6b.To characterize the deformation behaviors of described soft elastomer, the fractional order includes two portions: an initial const order for depicting the initial elastic response, and the declining order is utilized to count for the strain hardening behavior after crossing the yield point.The above analy indicates that the proposed discrete fractional model is more suitable for depicting nonlinear stress response of soft materials.

Validity of the Proposed Model under Compression
In this subsection, the proposed discrete variable fractional order model is applied describe the compressive stress-strain behavior of aluminum foams with different por ities.The loading rate is 5 mm/min, as reported by [54].Our model's predicted mechani responses under compression are shown in Figure 7a, with the calculated paramet listed in Table 5.

Conclusions
In this paper, we develop a nonlinear variable-order fractional constitutive model to predict the rate-dependence viscoelastic response of soft materials, in which the fractional operator with variable order is numerically approximated by the Grünwald-Letnikov derivative.Moreover, the physical meaning of the fractional order is also explored by introducing the principle of equivalent viscoelasticity between the fractional viscoelastic model and the time-varying viscosity model.The proposed discrete variable-order fractional model is applied to predict the mechanical response of several soft materials.The main conclusions can be summarized as follows.Firstly, the proposed variable-order fractional model can depict the tensile behavior of described soft materials, including the rate-dependent response and strain-hardening effect, which exhibits higher accuracy and more convinced reliability compared with the existing fractional models with piecewiseconstant and piecewise-linear orders.Secondly, the proposed model can also well predict the compressive stress-strain response.The fractional order contains a constant initial order, and a linear variable-order function is proved to be able to characterize the whole deformation process of soft materials, where the former is used to depict the initial elastic response and the latter is utilized to account for the strain hardening behavior after crossing the yield point.
coefficient mirrors the capability of a given viscoelastic material to resist deformation, where the solid-like material exhibits extremely high viscosity coefficients as claimed by Su et al. [56].As  increases, continuous attenuation of the resistance of flow contributes to the reduction of the elastic properties, making the material behavior approaches Newtonian fluid as  is close to 1.To summarize, the order  of fractional derivative is commonly identified as an index of material properties, which changes continuously from Hooke's elasticity to Newtonian viscosity as  varies from zero to one.

Fractal 17 Figure 3 .
Figure 3.Comparison of the experimental data of soft elastomer at stretch rate 0.065/s and the prediction of the discrete fractional model with variable order.

Figure 4 .Table 2 .Figure 3 .
Figure 4.The derived evolution of variable order.Table 2. Derived material parameters and error values of soft elastomer at stretch rate 0.065/s.

Figure 3 .
Figure 3.Comparison of the experimental data of soft elastomer at stretch rate 0.065/s and the prediction of the discrete fractional model with variable order.

Figure 4 .
Figure 4.The derived evolution of variable order.

Figure 4 .
Figure 4.The derived evolution of variable order.

Figure 3 .
Figure 3.Comparison of the experimental data of soft elastomer at stretch rate 0.065/s and the pr diction of the discrete fractional model with variable order.

Figure 4 .
Figure 4.The derived evolution of variable order.

Figure 5 .
Figure 5.The sensitivity of stress response to order parameters; (a) α 0 , (b) β 0 .(The experimental data was taken from the uniaxial tensile tests of VHB 4910 at a rate of 0.065/s, as reported in [51]).

Figure 6 .
Figure 6.(a) Comparison of the simulated results of different fractional models and the experimen data at stretch rates of 0.0065/s and 0.1613/s; (b) Derived order variance of the presented numeri stress response.(Experimental data were taken from [51]).

Figure 6 .
Figure 6.(a) Comparison of the simulated results of different fractional models and the experimental data at stretch rates of 0.0065/s and 0.1613/s; (b) Derived order variance of the presented numerical stress response.(Experimental data were taken from [51]).

Figure 7 .
Figure 7. (a) Comparison of the compressive experimental data of described aluminum foams w different porosities and the prediction of our model; (b) Derived order variance of the propo numerical stress response.(Experimental data were taken from [54].).

Figure 7 .
Figure 7. (a) Comparison of the compressive experimental data of described aluminum foams with different porosities and the prediction of our model; (b) Derived order variance of the proposed numerical stress response.(Experimental data were taken from [54].).

Figure A1 .
Figure A1.Variance of creep viscosity function over fractional order ranges from zero to one.

Table 1 .
Derived material parameters and error values of four uniaxial tests for nitrile rubber.

Table 2 .
Derived material parameters and error values of soft elastomer at stretch rate 0.065/s.

Table 2 .
Derived material parameters and error values of soft elastomer at stretch rate 0.065/s.

Table 2 .
Derived material parameters and error values of soft elastomer at stretch rate 0.065/s.

Table 3 .
Derived model parameters of constant-order (CO) fractional model.

Table 4 .
Derived model parameters of data fitting of the variable-order (VO) fractional model.

Table 3 .
Derived model parameters of constant-order (CO) fractional model.

Table 4 .
Derived model parameters of data fitting of the variable-order (VO) fractional model.

Table 5 .
Derived model parameters of predicted compressive stress response.