Viscoelastic multiscaling in immersed networks

Rheological responses are the most relevant features to describe soft matter. So far, such constitutive relations are still not well understood in terms of small scale properties, although this knowledge would help the design of synthetic and bio-materials. Here, we investigate, computational and analytically, how mesoscopic-scale interactions influence the macroscopic behavior of viscoelastic materials. We design a coarse-grained approach where the local elastic and viscous contributions can be controlled. Applying molecular dynamics simulations, we mimic real indentation assays. When elastic forces are dominant, our model reproduces the hertzian behavior. However, when friction increases, it restores the Standard Linear Solid model. We show how the response parameters depend on the microscopic elastic and viscous contributions. Moreover, our findings also suggest that the relaxation times, obtained in relaxation and oscillatory experiments, obey a universal behavior in viscoelastic materials.

One big challenge in science and engineering is the socalled multiscale modeling, namely, how constitutive relations of a material depend on its smaller scale interaction and composition [1,2].Such emergence problem, as described by P. W. Anderson [3], is even more remarkable in soft matter, where mesoscopic structures combine the atomistic and macroscopic frameworks and are responsible to the material elasticity and viscosity [4,5].In fact, the link among features of different scales forming the matter is highly nontrivial and subjected to intense research [6].The seminal work of H. Hertz about mechanical contacts, for instance, is still used nowadays to measure elastic properties of materials by analyzing how samples are deformed under applied stresses [7].Although it is one of the most used models to investigate the stiffness of a material, it fails to relate the macroscopic behavior with its inner parts.Besides, well-known analytical approaches such as the Maxwell, Kelvin-Voigt, and Standard Linear Solid models apply circuit analogies to propose a simple manner of how elastic and viscous terms mix up in viscoelastic materials [8].Again, these rheological models neglect any downscaling analysis.
Currently, many sophisticated methods have been developed to study the viscoelasticity of materials.Nanoindentation experiments, for instance, such as those with Atomic Force Microscopy, are extensively applied to investigate the mechanical properties at micrometer scale.Besides condensed matter, such technique has also been used to investigate soft materials [9][10][11].The way the sample responses to external stresses depends on its elastic and viscous terms.However, the complex structure and composition in soft matter systems make it challenging to investigate accurately the distribution of applied stresses through their interior.Many soft matter systems hold large colloidal aggregations or long (linear or cross-linked) polymeric chains, which have an essential mechanical role [4].In living cells, for instance, the cytoskeleton and the cytoplasmic fluid are constantly exchanging momentum with the extracellular surroundings [12,13].The knowledge of why these organelles influence the cell stiffness differently in healthy and sick cells, may lead to treatments for several diseases [14].
To investigate the link between mesoscopic and macroscopic features, we employ a molecular coarse-grained approach to reproduce rheological behavior of soft matter.We design a viscoelastic material composed of a particle-spring network immersed in a viscous medium so that the contribution of elastic and viscous interactions can be controlled at the mesoscopic level.This immersed network model is compatible with suspended polymer chains, colloidal aggregations, and other loadbearing structures as commonly found in soft materials [4].Although real materials unlikely present purely quadratic potentials, for small deviations, molecular interactions displaying a potential well can be described appropriately by a spring-like interaction.For instance, the equivalent spring constant of the Lennard-Jones interatomic potential can be calculated as 72 0 / 2 0 , where 0 is the height of the potential well and 0 the equilibrium distance between molecules [15].Moreover, the computational model we develop here is easily changed to support a more realistic potential in order to study more complex materials.We also consider the spring network as a regular lattice.Although complex structures may influence rheological properties of a material, simple geometries have been successfully applied to a large number of problems in mesoscopic models of condensed [16,17] and soft materials [18][19][20].
Our computational model consists of N spherical particles of diameter σ and mass m arranged in a face-centered cubic (FCC) lattice with a height given by H and a bottom and top plane given by 20σ × 20σsin(π/3).Every particle interacts with its 12 nearest neighbors by a spring potential.This elastic network is immersed in a medium of viscosity µ.Viscous effects are taken into account by considering friction between the particles and the medium.The network is indented by a hard sphere of diameter σ s , as shown in Fig. 1.This spherical indenter moves down with constant speed until a maximum indentation depth, δ 0 , is achieved.To avoid non-linearities, we consider δ 0 = σ, i.e., the maximum indentation is equal to the thickness of one layer of particles.Once in contact, the indenter applies a stress onto the particles.The bottom layer of the network is in contact to a hard substrate, where particles cannot move down but are free to slide horizontally.
The equation of motion of the i th particle is given by the following Langevin-like equation [21,22] where r i and v i are the position and velocity vectors of particle i, respectively.The term γ v i , known as Stoke's law, describes the drag force acting on spherical bodies, where the coefficient of friction is given by γ = 3πσµ [23].
The potential U i of particle i due to the neighboring particles and the indenter is given by where r ij = | r i − r j | is the distance and the equilibrium distance between the centers of neighboring particles i and j. k ij is the spring constant of the bond between these two particles.The last term of Eq. 2 applies only to surface particles, when they are in contact with the indenter.r is = | r i − r s | is the distance between the center of particle i and the center of the indenter, is a constant of energy, and α regulates the hardness of the indenter.The particles and the indenter exclusively interact through this hard-sphere potential with a high value of α.See in Ref. [24] the constants used in this work.
We solve the equations of motion in (1) through molecular dynamics simulations with a time integration done with the Velocity Verlet algorithm, and periodic boundary conditions applied to the horizontal plane [25,26].The contact force F is the sum of all collisions on the indenter, computed at each time as the indenter slowly presses down the sample.F increases with the indentation depth, δ, since more collisions occur on the indenter.These collisions cause a fluctuation in F , but an equilibrium state is reached in around 10 5 time steps.After equilibrium, we perform an additional 2 × 10 5 time steps to average the quantities of interest.
Before viscoelastic materials, we explore the elasticity in networks without local friction (γ = 0).Initially, we study homogenous networks, where all bonds have the same spring constant, k ij = k, and later we investigate the role of heterogeneities of local elasticity in the global stiffness.For simplicity, we consider heterogeneous networks made of only two types of bonds with spring constants given by k 1 and k 2 , respectively.Moreover, this binary network is built in two fashions depending on the distribution of k 1 and k 2 , namely, the double-layered and the binary random network.In the first, two homogenous networks, with spring constants given by k 1 and k 2 , respectively, are deposited one on the top of the other (see the inset in Fig. 2(a)).In the second, k ij is randomly assigned to k 1 or k 2 according to probabilities φ 1 and φ 2 , respectively (see the inset in Fig. 2(b)).In both cases, φ 1 and φ 2 also stand to the fraction of each bond and φ 1 + φ 2 = 1.For the limit cases of φ 1 = 0 and φ 1 = 1, both the double-layered and the binary random networks become homogeneous with spring constant k 2 and k 1 , respectively.
The distribution of stress into an elastic sample, due to an indentation δ, can be tricky [27], however the force experienced by a spherical indenter is well described by the Hertz model where E is the effective elastic modulus of the material [28] (see Fig. S1 in the Supplementary Material).The contact forces computed in elastic networks perfectly agree with the Hertz model, where E can be obtained by fitting the data with Eq. 3. The results for homogeneous elastic networks are shown in Fig. S2 in the Supplementary Material for different values of k.As expected for this simple network, the effective modulus of elasticity scales linearly with k.In addition, E is proportional to the sample height in the form E ∝ H −ζ , where ζ depends on the indenter geometry.Spherical and flat indenters, for instance, give ζ = 1/3 and 1, respectively.In fact, in the latter case, E represents the Young's modulus of the material and can be found analytically with springs in series, while in the other case, E seems to take into account, besides Young's modulus, the shear modulus of the material.The effective elastic modulus as a function of φ 1 is shown in Fig. 2(a), for double-layered networks, and in In order to study the contribution of viscous effects (γ > 0) in the macroscopic rheological behavior, we perform relaxation and oscillatory assays in homogenous viscoelastic networks.In relaxation experiments, the indenter is initially located above the network and moving down into the top surface.Before touching the sample, the force is zero.When the indenter touches the surface particles, at time t 1 = 0, the force rises until the maximum indentation depth, δ 0 , is achieved, at t 3 .Figure 3 shows the time evolution of the normalized contact force, f , for k = 500 and several values of γ.The loading time, τ l = t 3 − t 1 , is the time the indenter pushes the sample.After reaching δ 0 , the indenter stops moving, but the sample particles may continue to move subjected to friction (the dwell stage).For elastic γ = 0, the particles immediately rearrange themself to an equilibrium state and thus the force becomes nearly constant, except for some noise due to the perpetual undamped movement of the particles.On the other hand, in viscoelastic networks, this noise quickly vanishes while the force decreases continuously to the value of the undamped case, at t 4 , since particles need time to relax into the equilibrium state.For long times, networks with the same elastic properties experience the same contact force, regardless of the viscous contribution.
Our numerical results are compatible with the Standard Linear Solid (SLS) model composed of a Maxwell material (i.e., a spring of elastic modulus E 0 in series with a dashpot of viscosity η) in parallel to another spring of elastic modulus E ∞ , see Fig. 3.Such analytical model successfully describes rheological behaviors of a large number of soft materials, such as polymers [29], soft gels [30], and living cells [31,32].The relaxation function of this model is know as R(t) = E ∞ + E 0 e −t/τ , where the relaxation time is given by τ = η/E 0 .Initially, the effective modulus of elasticity is given by the sum of E ∞ and E 0 , but, for long times, the influence of E 0 vanishes, leaving E ∞ to dominate the elasticity of the material, regardless of the relaxation time.Notice that, E ∞ corresponds to the elastic modulus computed in elastic networks.The analytical force curves in the SLS model with a spherical indenter, f (t), can be obtained separately in the loading, f l (t), and dwell stages, f d (t), as follows where erf (t) is the error function and the constants are given by See more details in Section I of the Supplementary Material.
Combining our numerical results with Eqs. 4 allows us to find the relations between macroscopic (τ and E 0 ) In oscillatory experiments, instead of stoping the indenter after the maximum indentation is achieved, a sinusoidal strain is imposed with the indentation depth following the function, δ(t) = δ 0 + δ a sin(ωt), where δ a (δ a = σ) is the amplitude and ω the angular frequency of oscillation.The response force F also presents a periodic behavior but with a time delay that depends on viscoelastic properties, F = F 0 + F a sin(ωt + λ), where F 0 is the force at δ 0 , F a is the amplitude of the force and λ is the phase lag between force and indentation.Figure S3(a) and (b), in the Supplementary Material, show δ(t) and F (t), respectively, for a sample with k = 500 and three values of γ.The phase lag vanishes for γ = 0, as it is expected for elastic materials, but increases with γ, for viscoelastic materials.
The tangent of λ, defined as the ratio between the loss modulus, G = δa Fa sin(λ), and the storage modulus, G = δa Fa cos(λ), represents a quantitative way to define the fluidity of viscoelastic materials, i.e., tan(λ) < 1 leads to solid-like materials while tan(λ) > 1 to fluid-like materials.High values of ω favor the fluidity of the material, as shown in panel of Fig. 5 for several values of γ.only observe viscoelastic solids in all cases of k and γ considered in this work.Fluid-like behavior should be observed for very high values of ω, which requires a time step too small compared to the time of the experiment.Moreover, viscoelastic materials may present different relaxation times regarding of how the sample is deformed.
In this work, we obtained τ and τ o in relaxation and oscillatory conditions, respectively.Standard Linear Solids under oscillatory strain leads to the following behavior [8] tan which is used to fit those numerical data in Fig. 5 in order to obtain τ o , where A is a constant of the material.The main graph in Fig. 5 shows that τ and τ o depend on each other by the following exponential relation, τ = 0.466(e 0.08τo − 1).Each point in the graph of Fig. 5(b) represents a different material, with different k and γ, but they all collapse on that relation, suggesting some kind of universal behavior.
In conclusion, we studied macroscopic rheological behaviors in viscoelastic soft materials in terms of their microscopic elastic and viscous constituents.Different indentation experiments were performed to investigated how the local stiffness and the coefficient of friction in- duce different macroscopic behaviors.Elastic networks recovers the Hertz model for mechanical contact, where the effective modulus of elasticity, E, is obtained.In homogeneous elastic networks, we found that E ∝ H −ζ , where H is the network height and the exponent depends on the indenter geometry.For spherical and flat indenters, for instance, ζ = 1/3 and ζ = 1, respectively.Heterogeneous elastic networks, for simplicity, are assumed to be binary networks, where each bond holds one of only two possible values, k 1 or k 2 .When k 1 and k 2 are separated forming layers, E follows a cubic function with φ 1 , however, when k 1 and k 2 are randomly distributed, the higher order terms vanish and E changes linearly with φ 1 .Viscoelastic networks are compatible with the analytical Standard Linear Solid (SLS) model, which is described by two elastic moduli, E 0 and E ∞ , and a relaxation time, τ .Therefore, our work applies to all those materials that can be described by the Hertz and SLS models, such as polymers, gels and living cells.In oscillatory experiments, we also computed the phase lag between force and indentation, λ, and the relaxation time, τ o .We showed how all these macroscopic properties (E 0 , E ∞ , τ , λ, and τ o ) depend on the microscopic parameters.In addition, the interplay between different relaxation times, τ and τ 0 , obtained in relaxation and oscillatory experiments, follows an exponential behavior, regardless of the mesoscopic interactions.All viscoelastic networks simulated here collapsed on this curve, suggesting a universal behavior in viscoelastic materials.We acknowledge financial support from the Brazilian agencies CNPq, CAPES, FUNCAP.This work was also supported by the Serrapilheira Institute (grant number Serra-1709-18453).

FIG. 1 :
FIG. 1: Particle-spring network arranged in an FCC (facecentered cubic) lattice with height H and immersed in a viscous medium (not shown in the figure).Particles are bonded to their 12 nearest neighbors.A spherical indenter of diameter σs is used to press down the top surface of the network.

Fig. 2 (
Fig. 2(b), for binary random networks, for k 2 = 500 and different values of k 1 .E decreases with φ 1 for k 1 < k 2 , but increases with φ 1 for k 1 > k 2 , regardless the distribution of k 1 and k 2 .In double-layered networks, E is fitted by the cubic function E = 570.75+ aφ + bφ 2 + cφ 3 , shown in colored solid lines.On the other hand, when k 1 and k 2 are randomly distributed over the network, E changes according to the straight line E = 570.75+ dφ, shown in colored solid lines.The constants a, b, c, and d depend on k 1 .Interestingly, the random distribution of local stiffness vanishes those polynomial higher order found in (a).

FIG. 2 :
FIG. 2: The effective modulus of elasticity, E, as a function of the fraction of k1-bonds in heterogeneous elastic networks, φ1, for k2 = 500 and several values of k1.The double-layered network is in shown in (a) while the binary random network in (b).The inset of each graph shows how k1 and k2 are distributed in the network.The black dashed lines in each graph represent the elastic modulus of homogeneous networks, for k1 = k2 = 500, with E = 570.75.Cubic functions and straight lines, shown in colored solid lines, are used to fit each set of points in (a) and (b), respectively.

FIG. 3 :
FIG. 3: Relaxation experiments in viscoelastic networks, for k = 500 and different values of γ.The contact force increases in the loading stage and decreases in the dwell stage.The force is normalized by geometric constants (see Eq. 4 in the Supplementary Material).Our numerical results are compatible with the Standard Linear Solid (SLS) model given in Eqs. 4. A schematic illustration of the SLS model is shown in the inset, with a dashpot of viscosity η and two springs of elastic moduli E0 and E∞.

FIG. 4 :
FIG. 4: The macroscopic properties τ and E0, in (a) and (b), respectively, as a function of microscopic parameters, k and γ, in viscoelastic networks under relaxation assays.

FIG. 5 :
FIG.5:In oscillatory experiments, the tangent of the phase lag, λ, as a function of ω, is shown in the panel for k = 500 and different values of γ.Our results recover in a good agreement the behavior of tan(λ) given in Eq. 5, where the relaxation time in oscillatory experiments, τo, is obtained.The relaxation times, τ and τo, obtained in different experiments of deformation, depend on each other as follows, τ = 0.466(e 0.08τo − 1) (shown in black solid line), regardless of k and γ, as shown in the main graph.