The development of the formalism of movable cellular automata for modeling the nonlinear mechanical behavior of viscoelastic materials

The paper is devoted to the development of the formalism of the computational method of discrete elements (DEM) for describing the mechanical behavior of consolidated viscoelastic materials. We considered an advanced implementation of DEM, namely, the method of movable cellular automata (MCA). A feature of this implementation of DEM is the use of a generalized many-body formulation of the relations for the forces of element-element interaction. 3D numerical models of viscoelastic material with a spectrum of relaxation times (Kelvin and Maxwell models, the standard model of elastomers, and others) were developed within the formalism of MCA. The correctness of the developed discrete element formalism and its applicability for modeling the processes of deformation and fracture of viscoelastic materials under dynamic loading are shown using the standard model of elastomers as an example. The relevance of the results is determined by the prospects for the further development of DEM and its application to study and predict the mechanical response of viscoelastic materials of various nature under dynamic loading including contact problems.


Introduction
Discrete Element Methods (DEM) are reliable and effective tools for numerical analysis of mechanical behavior of particulate and consolidated materials under complex loading [1][2][3][4].
The key advantage of DEM in comparison with the methods of continual mechanics is the possibility of direct simulation of the complex processes of fracture, friction and wear (including the formation and evolution of the so-called third body) [5][6][7], the flow of granular media, etc. Despite the wide range of potential applications of the method, nowadays classic DEM implementations are used mainly for modeling brittle materials and granular media. This limitation is due to the simplified consideration of a discrete element as a pseudo-rigid body. The term "pseudo-rigid" means only local deformation of the element in the contact patch with the neighboring element. At the same time, the volume of the element remains unchanged, while other contact spots of the element "do not feel" the deformation in the considered contact. Mathematically, this simplification is formulated in terms of the pair-wise element-element interaction potentials/forces. The limited possibilities of the pair approximation do not allow the adequate description of the processes of deformation and fracture of dense granular and consolidated materials with complex rheological properties (including elastic-plastic and viscoelastic materials, permeable fluid-saturated materials, etc.). Therefore, at present, the most common tool for modeling the mechanical behavior of such consolidated materials are numerical methods based on continuum mechanics [8][9][10][11].
In recent years, there have been different efforts to solve this problem by introducing the deformability of discrete elements. We can mention the work of Brodu et al [12], Rojek et al. [13,14], Celigueta et al. [15] and others. In these works, the idea of the deformability of an element (that is, changes in its volume and shape during deformation of contacts) is implemented in various forms. However, the developed models are applicable for describing only granular media consisting of linearly elastic powder particles (the powder rheology corresponds to the simplest Kelvin viscoelastic model).
Recently, a more general approach has been proposed. It is based on the consideration of a discrete element as uniformly (homogeneously) deformable. This approach uses the generalized formulation of relations for the forces of element-element interaction in the many-body form [16,17]. Within the framework of this formulation, similarly, to the method of an embedded atom, the force of the central interaction of two elements is represented as a superposition of the pair-wise and volume-dependent components. The developed approach makes it possible to model both the consolidated and granular materials and to implement the complex rheological models of "structural units" (discrete elements). This implementation of the method is called the movable cellular automaton method (MCA). This name originated historically since the MCA method was first developed to simulate shock-induced exothermal chemical reactions in intermetallic powder mixtures [18]. To date, macroscopic models of elastic-plastic materials of various nature (metals, rocks, etc.) [16,17,19] have been implemented in the MCA method within the framework of the many-body approximation. At the same time, a wide class of technical and natural materials with viscoelastic rheological properties remained outside the scope of consideration. Examples of materials of this class are rubber, biological tissues (bone and soft tissues), and even some rocks. The implementation of viscoelasticity models makes it possible to apply the method of discrete elements for the numerical study of many topical problems, including the study of the laws of friction and wear in technical and natural tribological pairs. Moreover, hybrid techniques can be developed on the basis of viscoelastic discrete-element-based models to simulate unsteady processes in viscoelastic porous fluid-saturated materials (one of the variants of such hybrid technique was recently developed by the authors as applied to porous brittle materials [19]).
The goal of this paper is to construct a mathematical formalism of a 3D macroscopic discrete-element model of viscoelastic materials. The method of moving cellular automata was used as a computational basis for such a model.

Model description
In the framework of the DEM, the modeled material is represented by an ensemble of elements (material areas) of finite size and a given initial shape, which may change during loading. The size and volume of elements are determined based on specified packing [17]. When modeling a consolidated material, it is assumed that the elements interact through flat faces similar to flat-joint model [20]. The initial value of the area of the face is determined by the linear dimensions of the elements and packaging type.
The MCA method belongs to the group of explicit DEM, which is based on the principle of local force balance (force balance is determined for each discrete element rather than for the whole system). The MCA method uses the common approximation of equivalent disks (2D problem statement) or spheres (3D statement) [21,22]. This approximation means that the kinematics of a discrete element is assumed to correspond to the kinematics of the disk/sphere with the same mass and volume values as the element. The equation of the dynamics of a discrete element in this formulation is written in the form of Newton-Euler: ( where S ij is the area of the surface of interaction between elements i and j.
The key advantages of the MCA method are the use of the approximation of a homogeneously deformable element and the realization of this approximation by means of a many-body formulation of the relations for the element-element interaction forces.
In the framework of the approximation of a homogeneously deformable element/automaton, the stress state of the automaton is characterized by the averaged stress tensor: where R i is the radius of the equivalent disk/sphere, 0 ij S is the initial value of the contact area in the pair of the automaton i with the neighbor j, 0 i V is the initial value of automaton volume, , = x,y,z (XYZ is laboratory coordinate system), is projection of the unit vector connecting the centers of automata i and j onto the axis .
In recent papers [16,17], it has been shown by the example of isotropic elastic-plastic materials that the relations for the specific values of the forces of the central ( ij ) and tangential ( ij ) interaction of automata should be formulated in the form corresponding to the governing relations for the simulated material (that is, the relations connecting the diagonal/non-diagonal components of the stress and strain tensors or stress rate and strain rate tensors). Note that this is a many-body form of relations for  ij and  ij . It includes the invariants of the averaged stress tensor (3). In this case, the mechanical behavior of an ensemble of movable cellular automata rigorously corresponds to the mechanical response of the simulated material.
In the present work, we apply the developed formalism to construct a discrete-element model of viscoelastic materials. Below is a description of the numerical implementation of one of the most common viscoelastic rheological models, namely the standard model of elastomers. Note that this model is used to describe the mechanical behavior of viscoelastic materials of various nature, including rubber, bone and soft tissue.
The standard model of elastomers is a three-parameter model. It is a superposition of Kelvin element (spring with Young's modulus E K ) and Maxwell element (spring with Young's modulus E M , connected in series with a damper characterized by dynamic viscosity  M ). Specific values of resistance force of these elements are added (= K + M ), and the strains are the same (= К = М ). A scheme of the model in the simplest 1D problem statement is shown in Figure 1. In the three-dimensional problem statement, the defining relations between stresses and strains for the standard model are written in the following form (hypoelastic notation): Here,   are components of the Cauchy strain tensor, G K and G M are shear moduli of Kelvin and Maxwell elements,   are Kronecker symbols,  mean is mean strain, is mean stress.
As noted above, expressions for calculating the specific central force of the response of automaton to the mechanical action of neighboring automaton ( ij ) should be constructed by means of reformulation of the defining relations of the simulated material for the diagonal components of the stress and strain tensors. Similarly, when constructing expressions for specific tangential (shear) response forces ( ij ), we use the defining relations that connect the off-diagonal components of the stress and strain tensors. As applied to the viscoelastic model (the standard model of elastomers), the expressions for  ij and  ij are formulated in incremental form as follows: Here, the symbol  denotes the increment of the corresponding parameter per integration step, t is the integration step of the numerical scheme, is the increment of the normalized (by the radius of the equivalent sphere) distance from the center of mass of the automaton i to the central point of the interaction surface with the neighboring automaton j,  ij is the shear angle of the automaton i in the pair i-j [16,17], Here: The "cur" and "pre" indices here denote the belonging of the parameter values to the current and previous step of integrating the equations of motion, r ij is the change in the distance between the centers of the automata over a time step t, sh ij l  is the relative shear displacement of the automata of the pair i-j (it is calculated taking into account possible rotations of the automata [17,22]). The pairs of equations (8)-(9) are solved for the unknown increments of strains  ij and  ji (and also for  ij and  ji ). This allows us to then calculate the current specific values of the central and tangential interaction forces   Local fracture in MCA method (as well as in the conventional implementations of DEM) is modeled by means of change of the state of the pair of cellular automata (from "linked" state of the pair to "unlinked" state) [16,17,22]. This transition is governed by specified fracture criterion for the pair. In the study we used two-parametric failure criterion of Drucker and Prager in the following form: . (14) where  eq and  mean are equivalent and mean stresses at the surface of the interaction of linked automata, b= UCS / UTS , σ UCS is the uniaxial compressive strength of the material, σ UTS is the uniaxial tensile strength. These two uniaxial strengths of the material are input parameters of the fracture criterion. Failure criterion is applied to the pair of linked movable cellular automata. Details of numerical implementation of the criterion can be found in [16,17].

Verification of the model
The developed particle-based model of elastomers was verified by the comparing results of numerical simulations with well-known analytical solutions. For the standard model, there are analytical estimates of the dependences of the components of the complex elastic modulus on the frequency of periodic (alternating) uniaxial loading of a viscoelastic sample. When a deformation, varying harmonically with frequency , is applied to a sample of a viscoelastic material, then after a transient process the stress (specific force of the sample response to punch) will change according to a periodic law with the same frequency but with a phase lag : where E and E are real (elastic) and imaginary (loss) components of the complex Young's modulus of the material [23]. For the standard model, the steady-state values of the real and imaginary moduli under the condition of cyclic loading with fixed  0 are determined by the frequency of oscillation  and relaxation time T M as follows: Verification of the numerical model was performed on cylindrical samples. The sample (Figure 2a) was subjected to harmonic loading by periodically displacing the automata of its end surfaces along the cylinder axis (cyclic compression-tension) with a frequency . A time equal to 20 periods of oscillations was assumed as a characteristic time of transition of the system to the mode of steady-state deformation. This is sufficient for both "low frequency" loading (T M <0.1) and for "high" frequencies (T M >10). Real (elastic) and imaginary (loss) moduli were derived from the loading curve - by the standard method [23].  Figure 2b shows an example of a numerically obtained dependences of the real and imaginary parts and the absolute value of the complex Young's modulus versus frequency in logarithmic coordinates. The figure shows that the simulation results with high accuracy coincide with the analytical solution (16)- (17). Also, note some difference between numerical points and analytical solution at high loading frequencies (T M >5). This difference slightly increases with the further increase of . The reason for this is the fact that the analytical relations do not take into account the inertial characteristics and the final dimensions of the sample. At high cycling frequencies, a complex wave pattern is formed and continuously changes in the sample. It causes some difference in the registered dynamic elastic characteristics from those predicted analytically.
The frequency dependence of the dynamic values of the real and imaginary parts of the elastic modulus of a viscoelastic material is the result of its non-linear behavior under loading.
In the case of monotonic uniaxial loading with  =const, the nonlinearity of behavior appears in the logistic dependence of the effective (secant) elastic modulus on the applied strain rate  . Here, the secant Young's modulus E sec is determined as the ratio of the magnitude of the axial stress  to the magnitude of the corresponding axial strain . The value of E sec depends on both  , and . To study the character of this relationship, we performed a computer simulation of the uniaxial loading of the sample, shown in Figure 2a, with different strain rates  . One can see that in the process of deformation there is a decrease in the instantaneous (tangent) Young's modulus from the upper limit E (at the initial phase) to the lower limit E K . This is a consequence of the shear stress relaxation process. As the strain rate decreases, the nonlinearity of the material behavior increases.  Figure 3b shows the dependence of the effective (secant) Young's modulus of the sample E sec of the studied model material on the applied strain rate  . In all simulations, the secant modulus was determined at the moment the specimen fracture. The dependence is logistic and is approximated with good accuracy by sigmoid (18) with an exponent m2. The "midpoint" point of this dependence is a material parameter, which depends on the combination of material strength, elastic constants, and relaxation time T M . Note that the numerically derived logistic character of the secant Young's modulus dependence on strain rate coincides with an analytical estimate and is qualitatively confirmed by the results of experimental studies of bone tissues [24].
The results indicate the correctness of the developed numerical model and its applicability for studying the behavior (including fracture) of viscoelastic materials described by the standard model of elastomers under complex loading conditions including dynamic impacts.

Summary
In the paper, we proposed an approach to the construction of numerical models of viscoelastic materials in the framework of the method of homogeneously deformable discrete elements. It is important to note that in order to correctly simulate the nonlinear behavior of consolidated viscoelastic materials, it is necessary to use many-body formulations of relations for element-element interaction forces (this is realized in the framework of the movable cellular automaton method considered here). Although the paper presents a numerical implementation of the simplest viscoelastic model (standard model of elastomers), more complex models with a wide range of relaxation times can be implemented using Prony series in the same way. Moreover, many-body formulation of element-element interaction makes it possible to implement complex rheological models including viscoelasticity in the integral form with relaxation core functions R(t) [25], combined visco-elastic-plastic models [26][27][28] and coupled thermomechanical models (with the temperature dependence of dynamic viscosity [29]).
Some comments have to be done in connection with the shown logistic dependence of effective elastic modulus on strain rate under monotonic loading. A similar (sigmoidal) type of dependence of the effective Young's modulus on strain rate takes place for fluidsaturated elastic-brittle materials [19,30]. Indeed, in both cases, relaxation processes take place in the deformation process and reduce the magnitude of local stresses. In the first case (viscoelastic material) they are associated with the motion of defects and/or reconfiguration of the system of molecules. In the second system, the mechanism of stress relaxation is the transport of interstitial fluid in pore space.
Note that it is possible to build a "unified" (general) curve in terms of effective elastic modulus and the dimensionless parameter (Darcy number). Darcy number combines the characteristics of porosity and permeability of skeleton, fluid viscosity, pore pressure gradient and fluid flow distance (all these characteristics determine the possibility of interstitial fluid redistribution in pore space and corresponding stress relaxation in solid skeleton). The revealed generality of the strain rate dependences for viscoelastic and porous fluid-saturated material suggests that for viscoelastic materials there should be a similar dimensionless parameter. The identification of such a parameter is the goal of further research.
Moreover, the effective mechanical characteristics of fluid-saturated viscoelastic materials should also have a sigmoid dependence on some dimensionless parameter, which is a combination of the Darcy number and the corresponding dimensionless parameter of the viscoelastic skeleton. To carry out such an analysis, it is promising to develop a discrete element based coupled model of fluid-saturated viscoelastic materials by analogy with the implementation of the model of porous fluid-saturated elastic-plastic materials [19].
This research was supported by the Russian Science Foundation (Project 17-11-01232).