An Entropy Dynamics Approach for Deriving and Applying Fractal and Fractional Order Viscoelasticity to Elastomers

Constitutive model development for polymeric materials traditionally starts with a free energy function that contains information about the internal energy of the bonds between atoms and the entropy governing heat transport and the configurations of the polymer network as a function of macroscopic deformation [1–3]. Information about the internal degrees-of-freedom that span quantum, molecular, mesoscale, and macroscales remains extraordinarily difficult to explicitly quantify and therefore techniques that use entropy to quantify measures of uncertainty are important to approximate internal forces and transport phenomena (e.g., heat transport, chemical diffusion, photochemistry) without complete knowledge across all scales. It was rather fortuitous that the pioneering work by Shannon’s theory of information [4] was useful in modeling thermodynamic entropy [5,6]. In this work, we use the tools from information theory to better understand the nonlinear, rate-dependent mechanics of elastomers that exhibit complexity that we encode as fractals or multifractals (random fractals) [7]. We do this by using entropy dynamics [8] and modifying the framework to include fractional constraints that restrict the multiscale material structure to move along fractal dimensions. We relate fractal and fractional order operators to fractal structure to construct models that more accurately predict nonlinear, rate-dependent deformation in elastomers guided by an entropy dynamic framework.

The field of fractals and fractional order operators to materials science is extensive and has grown significantly in the past several years. Seminal works in the broader field of fractional calculus [9,10] and fractal diffusion in materials [11,12] have motivated more recent efforts to understand how both fractal and fractional order operators provide new insights into complex material behavior [13–20]. Mandelbrot’s concept of fractals [21] provides a measure of 1D, 2D, and 3D geometry or time series characteristics that are otherwise difficult to quantify on a finite domain. For example, the circumference of an island can be infinitely long given a repeated (fractal) structure over all length scales. These measures are made finite on a fractal domain. In terms of time scales, fractal behavior manifest in terms of rates of change of some physical property such as velocity that may follow a fractal response as the kinetic behavior of a material is viewed on different time scales. This fractal mathematical property is defined over all length scales that serves as an approximate model of a material with finite bounds in space and time. Upper and lower limits of the fractal domain have been considered to address these finite length scales [22,23,15].

We focus on constitutive relations near equilibrium while accommodating dissipation characterized by viscoelasticity. Elastomers exhibit complexities due to the vast number of polymer network configurations that result in nonlinear deformation and thermal dissipation upon time-dependent deformation. These complexities across molecular to continuum scales remain challenging as error propagation limits prediction across a broad range of deformation states and long time periods. In the limiting case of Gaussian processes, approximations can be made to develop idealized constitutive relations that lead to the very well-known neo-Hookean model of hyperelasticity [2,3]. Beyond moderate strain levels, elastomers exhibit nonlinearities characterized by non-Gaussian processes. Phenomenological methods can accommodate these nonlinearities in the stress–stretch behavior such as the Mooney–Rivlin and Ogden models [2]. However, these models are limited by their phenomenological nature requiring parameter estimation from data that are not well informed by the polymer structure. Other modeling approaches accommodate polymer structure and microscopic effects of polymer entanglement, cross-linking, and nonaffine deformation [24,25]. These models make assumptions on the microscopic volume elements that are homogenized over a continuum volume. Although these models have been successful in predicting finite deformation in certain elastomers, homogenization across continuum scales is based on a limited number of microscopic factors that often ignore non-Gaussian statistical distributions, direct relations to multiscale structure, and connections to viscoelasticity. In many cases, fat-tailed probability distributions of polymer molecular motion becomes important, and the fractal nature of materials can provide inputs to more accurately reflect derivative operations in the continuum scale constitutive model. This can have a significant influence on macroscale constitutive behavior as extreme events in the tails of the probability densities can have a cascading influence on macroscopic thermodynamic and kinetic properties.

West and Grigolini have argued that the fractal structure is often best approximated by fractional calculus operators [9]. The pervasive nature of fractal structures in nature, including materials, clearly illustrates a unique opportunity to understand how to apply these operators to develop constitutive relations in complex materials. One striking example is the fractional form of the continuity equation or fractional conservation of mass [26]. It was shown that when fluid density can be homogenized as a power-law function in space (e.g., fractal function), the Caputo fractional Taylor expansion of two terms exactly represents the density’s power-law relation. From a physical perspective, a power-law density relation can be motivated by fluids within a porous media [26]. This facilitates formulating an exact representation of the continuity equation by using first-order fractional derivatives. This same power-law density function would otherwise require an infinite number of integer-order derivatives from the Taylor expansion to obtain the same accuracy as one fractional Caputo derivative. The compactness gained by a small number of fractional derivatives must be balanced by the nonlocal properties of the fractional order operator that requires an integral calculation for the Caputo or Riemann–Louiville version [27]. In other works, Mainardi et al. [13] have given solutions to space–time fractional order balance equations that follow a stretched exponential under certain limiting conditions in space and time. From a Bayesian perspective, the problem can be flipped around and cast as an inference of a probability density from sparse fractal data. In doing so, we illustrate how fractional order kinematic constraints lead to a stretched exponential density that contains the limiting case of a Gaussian density. It is known that a microscopic field distributed about this stretched exponential leads to the fractal space–time diffusion equation. The result highlights how collecting sparse information in space and time can be used to formulate fractal partial differential equations to identify and simulate this type of material behavior. Importantly, such results motivate understanding the similarities and differences in local fractal derivatives versus nonlocal fractional derivatives for material modeling. We investigate how these characteristics inferred from Bayesian statistics are used to guide the application of fractal or fractional order operators used in developing rate-dependent constitutive equations in elastomers.

Given the subtle but important distinctions between fractal and fractional order operators in material modeling, we evaluate assumptions that constrain elastomeric materials to follow certain trajectories that are better predicted by fractal or fractional order operators in space and time. We apply entropy dynamics [8] to better understand how our assumptions, cast as a set of constraints about polymer network displacements, provide insights on appropriate derivative operators. We illustrate how a fractional (power-law) constraint of particle displacements leads to a stretched exponential in space, which, in turn, suggests that a fractal derivative is also ideal to approximate deformation via a fractal deformation gradient. Time-dependent properties, such as viscoelasticity, require additional assumptions on the irreversible characteristics of entropy generation. The covariance properties of particle interactions are introduced in the entropy dynamics framework to provide guidance on time-dependent particle interactions. Under Gaussian approximations, the variance or covariance of particle motion is linear in time, which can be accurately approximated by integer-order time derivatives. If the covariance evolves over a power-law in time, we show that a fractal time derivative within the diffusion equation exactly models such behavior. We compare this relation to fractional order viscoelasticity assumptions through experimental validation on the uniaxial stress–stretch behavior of dielectric elastomers, VHB 4910 and 4949. We also offer insight on when fractal and fractional order derivatives will give the same or different material model predictions.

The entropy dynamics framework quantifies conditional probability densities of future material particle positions given their original undeformed positions. The particles are broken down into two different sets. One set of particles is defined to be observable or controllable by an external loading device in a way that particles move in an affine manner based on boundary conditions. These displacements are typically defined by transformations that map changes on a boundary to the bulk volume. The set of particles is denoted by their Eulerian or deformed position in a three-dimensional space as $xα(t)$ for α = 1, …, n particles in some representative volume element (Eulerian frame volume). A second set of unobservable or uncontrollable particles are denoted by $yα(t)$ for α = 1, …, m such that the total number of particles in the material is N = n + m. These unobservable particles contribute to heat, residual effects, damage, or other irreversibilities that cannot be controlled by an external loading device. We apply entropy dynamics to construct a model containing constraints on both the controllable $xα$ and uncontrollable $yα$ particle positions and their time-dependent properties. These manifest through quantifying the uncertainty of a continuum homogenization of the particles as a field mapped onto the Lagrangian (undeformed) configuration that we denote by X. This requires first quantifying a Bayesian posterior density as a joint probability of the continuum homogenized kinematics: x = x(X, t) and y = y(X, t). Once the Bayesian posterior density is quantified, this naturally leads to a maximization of the likelihood of future deformation x(X, t) and internal state changes y(X, t) based on thermodynamic functions that are obtained from the posterior density. The thermodynamic functions provide the information required to quantify fractal or fractional hyperelastic and viscoelastic constitutive relations. We then extend the kinematic relations to the fractal domain to accommodate complexities that follow power-law behavior in space and time.

In the following section, we first outline how we apply entropy dynamics to obtain a Bayesian posterior probability. This is followed by the introduction of conditional probability constraints on x = x(X, t) and y = y(X, t) as well as fractional constraints that lead to thermodynamic potentials and entropy generation expression that relate the fractal elastomer network structure to kinematics constrained to move along fractal paths. We then numerically evaluate the model and validate it against viscoelastic elastomer data. Discussions are then given about the importance of constraints and its interpretation with respect to material structure when developing material models based on fractional order and fractal order derivatives. Conclusions are given in Sec. 5.

The constraint that we introduce in Sec. 2.2 is to minimize the relative positions between neighboring Lagrangian points: Δx = x(X + ΔX, t) − X and Δy = y(X + ΔX, t) − Y, where we denote Y as the undeformed position of the unobservable deformed state y. If a constraint based on Δx or Δy is added to Shannon’s entropy, we find that maximizing the Shannon entropy tends to broaden the posterior of possible particle positions while a constraint on the squared magnitude of Δx or Δynarrows the posterior. The posterior becomes a Gaussian distribution rather than a uniform distribution, and the width of the Gaussian depends on knowledge of the covariance matrix, which is used to constrain the vector products, ΔxΔx, ΔyΔy, and ΔxΔy as described in the following section. The quantitative form of the Gaussian requires knowledge of covariance matrices that represent distributions of the self-interactions of ΔxΔx and ΔyΔy and the interactions ΔxΔy. These covariance matrices may come from theoretical conjectures, sparse observations from experiments, or high fidelity simulations. Entropy dynamics provides a tool to evaluate the uncertainties in terms of estimates of the covariance matrix, which may include time invariant and time-varying components. This provides advantages in distinguishing reversible, path-independent behavior (time-independent covariance) from irreversible, hysteretic behavior (time-dependent covariance). It also provides a way to connect information theory and information entropy with thermodynamic entropy and constitutive relations for non-Gaussian processes. This latter effect will be explored using fractional order constraints, as opposed to a quadratic constraint, to relate fractal structure to fractional order properties. This will impact both the reversible hyperelastic behavior and irreversible viscoelasticity.

We briefly illustrate how this Gaussian density leads to a neo-Hookean, viscoelastic constitutive model to highlight how constraints are used within the information theoretic approach. Importantly, our selection of the quadratic constraint restricts the posterior density to material displacements that are distributed about a Gaussian density. This neglects the possibility of nonlinear effects that may manifest as extreme events such as internal displacements during unraveling of polymer entanglement and random walks along a fractal polymer network, which are later described in Sec. 3.1.

To illustrate, we evaluate the characteristics of an ideal Gaussian polymer network given the Gaussian posterior that was obtained from the tensor product constraint in the cost function. Simplifications are made by assuming all material points are observable meaning all material points are distributed about Gaussian density and remain Gaussian for all deformation states considered. This means y can be neglected, and Σ_xy and Σ_yy are irrelevant. We note that the assumption of Gaussian behavior is applicable to all future material displacements relative to the undeformed state. This is inherent in the fixed constraint Σ_xx within Eq. (4). It is well known that the Gaussian network neglects significantly large deformation and is only applicable for stretch ratios before strain hardening occurs [3,24]. This limitation is addressed in Sec. 3.

In Sec. 3, we will change the tensor constraint from a quadratic to a fractional constraint to accommodate fat-tail probability distributions of particle displacements often exhibited by fractal structures. We will also introduce rate-dependent effects governing viscoelasticity. This provides a means to correlate non-Gaussian polymer network configurations with complex constitutive relations where rate-dependent polymer deformation may follow fractal space–time complexities.

Large elastomer deformation can lead to non-Gaussian polymer network distributions where softening followed by hardening in the stress–stretch behavior may be observed as the underlying polymer network deforms significantly. In this section, we extend the information theoretic framework that resulted in a neo-Hookean free energy in Sec. 2.2 to fractional order deformation constraints to accommodate nonlinear stress–stretch behavior often observed during large elastomer deformation.

To establish translational and rotational invariance, we restrict the radius of some sphere to a finite set that we define by the characteristic length l₀, which is used to define a representative continuum volume element of sufficient size to represent the fractal structure. If we assume that β = ν from (12), this relation becomes

We introduce viscoelasticity into this model by adjusting the constraints for particle motion as given by the last term in Eq. (16). This constraint on relative particle motion (⁠$dμ$⁠) is divided into observable (⁠$dμ(x)$⁠) and unobservable (⁠$dμ(y)$⁠) terms that are independent. The unobservable relative displacements contribute to the viscoelastic effect described in Sec. 3.3. This gives three separate constraints in the cost function with a fractional variance on the observable $(Σxxν2)$⁠, unobservable $(Σyyν2)$⁠, and their interactions $(Σxyν2)$ fractal relative displacements.

This section re-visits the integer Taylor series of a function and then compares it to the fractal Taylor and fractional series expansions space. To simplify the analysis, we illustrate its effect in terms of a general scalar function f(x). Its use is given here to demonstrate the utility of fractal derivatives in space and time for systems that follow power-law (fractal relationships). This has implications on predicting the kinematics of elastomers upon mechanical loading.

It is clear that the first two terms of the series provide a linear approximation of the function around x = x₀. For nonlinear functions, extrapolation in space and time limits accuracy, dependent upon the magnitude of the nonlinearity.

One significant difference between the fractal Taylor series and Caputo-based fractional Taylor series is their basis of expansions. As mentioned earlier, fractal series relies on $(xα−x0α)$ type of basis, while the fractional series uses (x − x₀)^α. This difference in basis results in a slightly different approximation of functions. More importantly, depending on the functions being approximated, one method can outperform the other as shown in Fig. 1. For example, two-term fractional Taylor series is exact for f(x) = q + (x − x₀)^p but not for $g(x)=q+xp−x0p$ (again p > 0), while the opposite holds true for the fractal two-term Taylor expansion. This difference in basis and approximation is explored by applying both operators in time domain to model viscoelasticity.

Our assumption of the fractal time derivative model has assumed that the covariance matrix governing the internal state increases in time proportional to a special power-law function. In the limiting case, α = 1, we obtain the classical diffusion process that follows a Gaussian density. We have generalized this to power-law relations with the limiting case of a Gaussian density when the time order is α = 1 and the spatial order is ν = 2. We give additional support for using a fractal time derivative by starting with a stretched exponential probability density in space and time. It is shown that when the internal state μ is homogenized over a stretched exponential, fractal time derivatives in space and time lead to the fractal diffusion equation. This extends results illustrated by Falconer [33] for the limiting case of a Gaussian distribution.

Stress–stretch experiments on the dielectric elastomer very high bond (VHB) 4910 and 4949 manufactured by 3M are used to validate the model. A series of uniaxial cyclic loading/unloading stress/stretch tests were conducted to observe the change in steady-state hysteretic stress behavior while varying the stretch rate. Details on the experimental methods can be found elsewhere [34]. These stress/stretch results are used to compare a fractal order linear viscoelastic model to the fractional order viscoelastic model discussed here.

The total stress can then be solved by combining Eqs. (47), (48), and either (49) or (50) depending on whether the fractal or fractional order derivative, respectively, is being used.

Using Bayesian statistical analysis, the parameter set, $θ=[η,α,f011,f012,f022,ν]$⁠, for both the fractional and fractal viscoelastic models can then be inferred, where η and α describe the viscoelastic behavior and $f011$⁠, $f012$⁠, $f022$⁠, and ν describe the hyperelastic model. The delayed rejection adaptive metropolis algorithm is employed based on the Markov Chain Monte Carlo sampling technique to construct posterior distributions [35–37]. All the chains look to have converged to a stable mean with approximately the same statistical distribution along any point in the sample space for about 2.5 × 10⁵ iterations for the fractional and fractal order models. Figure 2 shows the marginal posterior densities for both the fractional and fractal models for VHB 4910 when calibrated at the fastest stretch rate tested. All six parameters look to have a Gaussian distribution, and except for η, the parameter distributions overlap almost completely. The difference in the values for η can be explained by the difference in Eqs. (44) and (45). Similar plots for VHB 4949 are provided in Fig. 3.

Mean parameter values can be found in Table 1 for each tested stretch rate. Overall, each parameter set is similar to the one previously discussed where all the parameters except for η have similar values and η is different by the previously discussed factor for the fractional and fractal models.

Using the determined parameter values, a model fit for the VHB 4910 experimental data at the fastest speed tested can be seen in Figs. 4(b) and 5(b) along with calculated 95% credible and prediction intervals for the fractional and fractal order models, respectively. The same procedure is then used to obtain separate parameter sets for each tested rate. The calibrated model fits at 6.7 × 10⁻⁵ can be found in Fig. 4(a) for the fractional order model and Fig. 5(a) for the fractal order model. The fits between both models look to be almost identical at each tested rate. The fractional and fractal model fit on VHB 4949 data are provided in Figs. 6 and 7, respectively, and mean model parameters for both models are listed in Table 2.

To more accurately assess the fits of the calibrated plots, a sum-of-squares error between the experimental data and model using the mean parameter values can be determined. Tables 3 and 4 give the errors for each calibrated model fit at each stretch rate. By comparing the fractional and fractal errors, we see that the errors are within 0.1% of each other with the fractal order model slightly outperforming the fractional order model 50% of the time. We also note that the mean estimates for η and α converge to consistent values at higher stretch rates. This is because viscoelasticity is more dominant at high stretch rates. The inferred η and α at higher stretch rates are thus more representative of viscoelastic behavior. Moreover, fixed η and α from higher stretch rates can accurately predict hysteresis at lower stretch rates; see Ref. [14].

A new modeling framework for polymer mechanics has been developed using entropy dynamics. A cost functional that combined Shannon informational entropy with a fractional order constraint of relative particle displacements was introduced to model non-Gaussian deformation processes that often occur in materials that exhibit the fractal structure. The fractional order constraint leads to a Bayesian posterior density that is a function of observable and unobservable internal state deformation variables. Through the use of Boltzmann entropy, we obtain a thermodynamic entropy and entropy generation functions from the time-varying Bayesian posterior that accommodates nonlinearities associated with fractal or power-law characteristics in both space and time. The Bayesian posterior obtained from this cost function has the form of a stretched exponential, which accommodates non-Gaussian deformation for hyperelastic and viscoelastic behaviors. Such fat-tailed probabilities are often seen in complex systems that are fractal in nature. While materials have finite length scales and may not follow an exact power-law or fractal structure, Bayesian inference gives an estimate of uncertainty in the model parameters that result from the power-law assumption.

The current study applied fractional and fractal order operators to characterize the viscoelastic behavior of polymers under constant strain rate loading. Under such loading, local fractal time derivatives and nonlocal fractional order time derivatives were shown to be equivalent excluding a multiplicative constant. The two viscoelastic models were calibrated to experimental data from two different polymers and both illustrated equivalent estimates of viscoelasticity over a broad range of deformation rates, indicating their ability to characterize the materials’ viscoelastic behavior accurately. This new approach has important implications for providing stronger connections between complex multifractal structure and power-law hyperelastic and viscoelastic properties that can be used to design and optimize polymer-based materials and structures, as traditional linear viscoelastic models may not capture their nonlinear behavior. Further research could explore applying these models to a broader range of polymers and other loading conditions to assess distinctions between model predictions of data when using local fractal or nonlocal fractional order operators. The study highlights the potential of fractional and fractal order models to advance our understanding of complex material behavior.

W. O. and B. R. P. appreciate support from the Department of Defense Basic Research Program for HBCU/MI (W911NF-19-S-0009). E. S. acknowledges the support from the DOD SMART Scholarship. The research of S. M. was partially supported by the National Science Foundation grant DBI 2109990. Any opinions, findings, and conclusions or recommendations expressed in this publication are those of the authors and do not necessarily reflect the views of the funding sponsors.

There are no conflicts of interest.

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.