Sublinear drag regime at mesoscopic scales in viscoelastic materials

Stressed soft materials commonly present viscoelastic signatures in the form of power-law or exponential decay. Although exponential responses are the most common, power-law time dependencies arise peculiarly in complex soft materials such as living cells. Understanding the microscale mechanisms that drive rheologic behaviors at the macroscale shall be transformative in fields such as material design and bioengineering. Using an elastic network model of macromolecules immersed in a viscous fluid, we numerically reproduce those characteristic viscoelastic relaxations and show how the microscopic interactions determine the rheologic response. The macromolecules, represented by particles in the network, interact with neighbors through a spring constant k and with fluid through a non-linear drag regime. The dissipative force is given by γvα, where v is the particle’s velocity, and γ and α are mesoscopic parameters. Physically, the sublinear regime of the drag forces is related to micro-deformations of the macromolecules, while α ≥ 1 represents rigid cases. We obtain exponential or power-law relaxations or a transitional behavior between them by changing k, γ, and α. We find that exponential decays are indeed the most common behavior. However, power laws may arise when forces between the macromolecules and the fluid are sublinear. Our findings show that in materials not too soft not too elastic, the rheological responses are entirely controlled by α in the sublinear regime. More specifically, power-law responses arise for 0.3 ⪅ α ⪅ 0.45, while exponential responses for small and large values of α, namely, 0.0 ⪅ α ⪅ 0.2 and 0.55 ⪅ α ⪅ 1.0.

Stressed soft materials commonly present viscoelastic signatures in the form of power-law or exponential decay.Understanding the origins of such rheologic behaviors is crucial to find proper technological applications.Using an elastic network model of macromolecules immersed in a viscous fluid, we numerically reproduce those characteristic viscoelastic relaxations and show how the microscopic interactions determine the rheologic response.We find that exponential relaxations are indeed the most common behavior.However, power laws may arise when drag forces between the macromolecules and the fluid are sublinear, which is related to micro-deformations of the macromolecules.
Purely elastic and purely viscous behaviors are limited cases of constitutive equations of materials [1].Actual substances may deform and flow, but one of these attributes usually dominates the other, depending on the applied conditions.This solid-liquid duality has teased researchers since at least the 19 th century.Back then, pioneers such as James Maxwell and Ludwig Boltzmann proposed analytical models based on series and parallel associations of springs and dashpots to explain the peculiar characteristics observed in silk, glass fibers, and steel wires [2,3].The effective response of such early models invariably presents exponential relaxation decays, regardless of how springs and dashpots are connected.However, these simple approaches are only suited for some viscoelastic materials nowadays.
In modern society, soft matter is ubiquitous and broadly accessible.The emergence of such complex materials has triggered new theoretical models and the improvement of proper experimental techniques to explain and control their viscoelastic properties [4,5].Nanoindentation methods, such as Atomic Force Microscopy, have become essential to characterize viscoelastic features at micro and nanometer scales by probing materials with nano-sized indenters [6].The characterization of viscoelastic materials attempts to determine the relaxation function that possesses both qualitative and quantitative information.
Exponential and power-law relaxation functions are the two major types of experimentally probed responses.Polyacrylamide gels [7,8] and aqueous solutions of cationic surfactants [9], for instance, present exponentiallike responses with a relaxation time for the material to achieve a new equilibrium configuration.On the other hand, living cells [10], microgel dispersions [11], soft glassy materials [12], and hydrogels [13] present a timeinvariant power-law-like behavior.As observed in elastic materials [14,15], macroscopic physical parameters are intrinsically connected to their microscopic interactions and structures [16][17][18].
Power laws and exponentials arise in many physical phenomena having a deep origin in their dynamic pro-cesses.For instance, in non-additive entropy systems, many physical variables are described by power-law distributions instead of the traditional exponential functions in the counterpart entropy [19,20].Exponential and power-law canonical distributions emerge naturally regarding whether the heat capacity of the heat bath is constant or diverges [21].Moreover, power laws are associated with emergence phenomena where exponents display scaling behaviors as they approach criticality [22].Systems with precisely the same critical exponents belong to the same universality class, and a small set of universality classes describes almost all material phase transitions.
One of the challenges in material science is linking the physical mechanisms at microscopic scales to macroscopic functional behavior.This approach is especially relevant for soft matter because properties on the molecular scale are linked to conformational and compositional fluctuations on the nanometer and micrometer scale and, in addition, span many orders of magnitude in length [23,24].Soft matter holds rich structures and various interactions at the mesoscale, where thermal energy per unit volume is negligible, in contrast with the high energy density stored in atomic bonds of crystalline structures [25].While exponential materials can be modeled by an association of springs and dashpots, such as the so-called standard linear solid model, power-law materials are usually obtained by fractional rheology [26,27] or glassy rheology models [28,29].However, these models cannot explain the connection between macroscopic responses and their elastic and viscous components.
We design a model of viscoelastic materials composed of an immersed elastic network of macromolecules to study how mesoscopic interactions influence macroscopic rheological behavior in soft materials [30].We assume non-linear hydrodynamic drag forces act between the macromolecule and the fluid, where the contribution of elastic and viscous interactions are controlled at the mesoscopic level.By changing the physical parameters of elastic and drag forces, we obtain materials with exponential or power-law relaxations or then an intermediary behavior for responses not clearly characterized.Our results show that exponential behavior is, in fact, the most common regime of deformation, being described by the standard linear solid model.Power-law responses are exceptional outcomes for a particular range of sublinear drag forces.
Our model consists of N spherical particles of diameter σ and mass m arranged in a face-centered cubic (FCC) lattice with dimensions x × y × z given by σH sin(π/3) × σH sin(π/3) × σH, where H is the number of layers in z direction, as shown in Fig. 1(a).Every particle interacts with its twelve nearest neighbors through an elastic potential with an effective spring constant k.The spring network is immersed in a viscous fluid, where drag forces act on moving particles.In this coarse-grained approach, the particles represent macromolecules commonly found in suspended polymer chains, colloidal aggregations, and other load-bearing structures of soft matter.
We perform computational indentation assays to probe this network's effective viscoelastic properties.Firstly, a rigid spherical indenter presses down the network at a constant rate.This loading stage is done during a time τ l until a maximum indentation depth δ max is achieved.After that, called dwell stage, the indenter stays still while the network rearranges towards a minimum energy con-figuration.Particles in the bottom layer are not allowed to move along the z-axis (where deformation is applied) but are free to slide horizontally.To avoid finite-size effects, we limit the maximum indentation to less than 10% of the network height [31], and thus we apply δ max ≈ σ.
The equation of motion of the i th particle, at position r i , is given by the following equation [32,33] where U i is the interaction potential of particle i due to other particles and the indenter, given by ( The summation in the first part runs over the neighbors, where r ij = | r i − r j | and are, respectively, the distance and the equilibrium distance between the centers of particles i and j.The last term in Eq. ( 2) represents a hard-core potential applied only to those particles in contact with the indenter, where is an energy parameter, and r is = | r i − r s | is the distance between the center of the particle and the indenter.The exponent ξ must be large enough to keep the stiffness of the indenter.The last term of Eq. ( 1) represents a generalized drag force acting oppositely to the particle velocity, v i = v i vi , with magnitude given by γv α i , where γ and α are related to the particle geometry and the fluid properties in which the particles are immersed.Dissipation vanishes for γ = 0, leading to purely elastic networks where our model reproduces the well-known Hertz behavior for mechanical contacts [30].However, when local friction becomes relevant, γ > 0, the model may produce distinguished behaviors of viscoelastic materials.Notice that α = 1 and α = 2 represent typical values for drag forces acting on rigid structures.The linear regime, known as Stoke's law, arises for small Reynolds numbers when viscous forces dominate over inertial forces, where γ is proportional to the medium's viscosity and the particle's diameter.On the other hand, the quadratic drag is dominant for large Reynolds numbers.In this case, γ is proportional to the medium's density and the cross-sectional area of the particle.
The physical origins of the sublinear regime (α < 1), however, are entirely different and are related to the deformability/adaptability of bodies subjected to drag forces [34].Many living beings present sublinear drag behaviors, where deformability is a survival strategy to protect their fragile structures under hydrodynamic conditions.Characterizing the drag exponents in deformable systems is recurrent in botany, aerodynamics, and hydrodynamics [35][36][37][38].Typical drag exponents for algae are as small as 0.34 [39].On the other hand, tulip and willow oak trees are more rigid and present exponents close to unity, α = 0.92 and α = 0.94 [36], respectively.The limiting case of α = 0 corresponds to constant frictional forces, regardless of the particle velocity.In our model, micro deformations are associated with a new molecular arrangement of the macromolecules.The numerical solution of Eq. ( 1) is performed through molecular dynamics simulations [40,41] with periodic boundary conditions applied to the horizontal plane, and time integration is done with the velocity Verlet algorithm with time step dt.See [42] for the parameters used.The force is computed as the sum of all collisions on the indenter at each time step.Each set of k, γ, and α represents a specific material and defines the macroscopic viscoelastic responses.To investigate how microscopic properties lead to the rheological behavior of the entire network, we perform simulations varying the spring constant between 100 and 1000, the drag constant between 10 and 100, and the drag exponent between 0.0 and 1.0, totalizing 2100 different networks.
Figure 1(b) shows typical force curves for γ = 80 and k = 800 and different values of the drag exponent.F (t) is the computational analogous of an indentation assay used to characterize macroscopic responses of actual materials.In the loading stage, the indenter slowly deforms the network until the contact force reaches its maximum value, which is inversely proportional to α.In the dwell stage, when the indenter stays still, F (t) relaxes until part of the mechanical energy is lost through friction.The timescale of such dissipation depends on α in a highly complex fashion.For α = 2, dissipation is fast enough so that the response behavior is qualitatively identical to a perfectly elastic network.For α < 1, drag forces slowly remove energy from the network leading to a viscoelastic decay.
In a continuum approach, the analytical timedependent contact force of a viscoelastic sample indented by a spherical punch follows the convolution integral [30] where R(t) is the relaxation function and δ(t) the indentation depth history.The contact force in Eq. ( 3) is normalized by the constant 4 , where ν is the Poisson coefficient.Typical relaxation behaviors for viscoelastic materials are given by where R P (t) is the relaxation model for power-law materials with an exponent β and R E (t) for exponential materials with relaxation time τ .In both models, E ∞ is the elastic modulus at considerable times when the material is completely relaxed, and ∆E is the difference between the maximum value, at t = τ l , and E ∞ .We assume that materials become perfectly elastic after some time, although this long-time elasticity plateau is only sometimes observed in actual materials [43].However, micro-deformations are allowed only due to the movement of the pseudo-atoms, which is represented by the sublinear drag regime.Solving Eq. (3) with Eqs.(4) in the dwell stage leads, respectively, to where the constant a ≈ ∆E Γ(1−β) is obtained by expanding the incomplete beta function, Γ(1−β) is the gamma function, and b = 3τ 2τ l ∆E and c = 3 , where erfi( τ l τ ) is the imaginary error function of x.
Fitting the dwell part of the force curve allows us to map the macroscopic mechanical properties (E ∞ , ∆E, τ , β) with the mesoscopic parameters (k, γ and α) as done previously for α = 1 [30].Here we focus on finding the qualitative rheological behavior rather than describing relations among parameters.Figures 1(c) and (d) show the same force curves as in (b) for α = 0.45 and α = 0.85, respectively, where the curves are fitted both with F P (t) and F E (t).Clearly, the α = 0.45 case is better fitted with the power-law model, while α = 0.85 with an exponential.
The mean-square error determines the goodness-of-fit parameter between the obtained numerical force curve F (t) and the analytical force models given in Eqs. 5 and 6, where N p is the number of points in the force curves and M stands for exponential (E) or power-law (P ) model type.This statistical index, calculated for each combination of k, γ, α, is used in a K-means method as an unsupervised clustering strategy [44,45] to classify the computational force curves as either an exponential or a power-law behavior, or even a transitional regime that cannot be clearly classified.In this machine learning process, the classification of the material takes into account not only the values of χ E and χ P but also the trends and distributions in the phase space of χ E and χ P .Figure 2(a) shows the graph visualization of this clustering process.The exponential behavior is indeed the most common one, as shown in the big cluster in Fig. 2(a), but small groups of power-law materials are also observed.The parallel coordinates plot in Fig. 2(b) summarizes the rheological outcomes.Each line passes through every combination of k, γ, and α ending up in one of the three boxes representing the response behavior of the network.
To understand why small clusters of power-law materials form in Fig. 2(a), we must investigate the impact of considering sublinear drag regimes.Fig. 3 shows the normalized probability P (α) of finding PL, EXP, or TR behavior for a given α.In panel (a), P (α) is computed for all values of k and γ considered here.The three distributions strongly overlap, making classification a challenging exercise.In panel (b), on the other hand, we remove small values of k and γ.For networks not too soft, k ≥ 800, and not too elastic, γ ≥ 80, those distributions split for different regions of α, and the drag exponent becomes the central controller to characterize the macroscopic behavior.Exponential behaviors are found mainly for α 0.2 and α 0.55, while the relaxation is a power law for α between 0. as an exponential or a power law.
For those cases where the network presents a power-law behavior, we show in Fig. 4 the relationship between the macroscopic relaxation exponent β and the drag exponent.For each value of α, there is a small dispersion distribution of β whose mean range is between 1.05 and 1.35.Actual materials exhibiting power-law relaxation usually present smaller exponents and exhibit structural disorder and metastability [27][28][29]46].Computational investigations of disordered two-dimensional networks obtained relaxation exponents between 0.5 and 0.75, depending on the network arrangement [18].Our simulations exhibit macroscopic relaxation exponents above 1.0, mainly because our viscoelastic solid model is structurally ordered and stable, restricting the viscoelastic responses to faster relaxation regimes.Our results clearly show that structural disorder is not a mandatory ingredient for powerlaw-like viscoelasticity and can only change the exponent.
In conclusion, we show why viscoelastic materials present power-law or exponential relaxation responses.Using molecular dynamics simulations, we perform numerical indentations onto a network of interacting macromolecules immersed in a fluid and reproduce typical vis- coelastic signatures as those found experimentally.The macromolecules interact with the fluid through a nonlinear drag regime given by γv α , where γ and α are mesoscopic parameters, and v is the particle's velocity.For each set of the interacting parameters, we classify the macroscopic viscoelasticity using an unsupervised clustering algorithm according to the type of relaxation of the deformed network, namely, exponential or power-law relaxations or transitional behavior between them.While exponential behaviors are predominant, power laws may arise in the sublinear regime.In fact, the drag exponent alone may explain the macroscopic viscoelastic relaxation for materials not too elastic nor too soft.More specifically, power-law responses are found for 0.3 α 0.45, while exponential responses for 0.0 α 0.2 and 0.55 α 1.0.
The authors acknowledge the financial support from the Brazilian agencies CNPq, CAPES, and FUNCAP.

FIG. 1 :
FIG. 1: (a) Viscoelastic model made by a face-centered cubic (FCC) lattice of macromolecules immersed in a viscous fluid (not shown in the image).Each spherical particle of diameter σ represents a collection of molecules interacting with an effective spring constant k with its twelve closest neighbors.A rigid spherical punch of diameter σs indents the network from top to bottom.H is the number of layers in the vertical direction.(b) Computational force curves for γ = 80 and k = 800 and several values of α exhibiting different viscoelastic relaxations.The cases for α = 0.45 and α = 0.85 are shown in panels (c) and (d), respectively, where the dashed and the dotted lines represent fit with the exponential and power-law model of Eqs. 5 and 6.

FIG. 2 :
FIG. 2: (a) Graph visualization of the clustering process using K-means method.Each dot represents a computational experiment for a given combination of k, γ, and α, and the colors represent different relaxation outcomes, namely, power-law (PL) in blue, exponential (EXP) in yellow and transitional behavior (TR) in red.(b) Parallel plot showing how different values of the mesoscopic network parameters lead to different types of viscoelastic relaxation.Each line crosses a combination of k, γ, α, and the corresponding rheological behavior.

FIG. 3 :
FIG. 3: The probability that a given α leads to a power law (blue bars), an exponential (yellow bars), or a transitional response (red bars).The solid lines show the probability density computed by the KDE (kernel density estimation) algorithm.Panel (a) shows results for the whole range of k and γ, while panel (b) considers only networks with k (k ≥ 800) and γ (γ ≥ 80).

FIG. 4 :
FIG.4:For those networks classified as a power-law material in Figure2(b), we show the relationship between the relaxation exponent β and the mesoscopic drag exponent γ.