A constitutive modeling approach to simulate time‐dependent phenomena in rubber‐like materials

Industrial rubber materials and similar rubber‐like materials exhibit complex mechanical behavior including nonlinear stress‐strain behavior up to large strains, hysteresis, Mullins effect, permanent set, as well as time‐dependent effects, like relaxation and rate dependency. The modeling of all these phenomena is a challenging task. One constitutive model that is capable to capture many effects in the large strain regime is the so‐called model of rubber phenomenology (MORPH). However, it is not able to capture any of the mentioned time‐dependent phenomena. Thus, this contribution deals with the extension of this model to simulate relaxation and rate dependency of rubber‐like materials. Thereby, a basic modeling approach is proposed, where only a slight change in the original system of tensorial equations is necessary. This approach is further extended by applying not only a constant relaxation time, but also loading history‐depended relaxation properties. For numerical simulations, a time discretization has to be applied to the constitutive equations, which is briefly shown. Finally, the extended MORPH model is applied to two different materials: an industrial rubber material and a polyurethane based adhesive. In both cases, the extended MORPH model is adjusted to uniaxial tension test data by parameter identification procedures. It will be shown that basic time‐dependent phenomena of rubber‐like materials can be captured well with the proposed extensions.


INTRODUCTION
Due to their outstanding mechanical properties, rubber-like materials, like filled rubbers or thermoplastic elastomers, are widely used in many technical applications.Thereby, their mechanical behavior is characterized by a comparably low stiffness and high flexibility.Thus, these materials can sustain very large strains.However, they also reveal comparably complex phenomena within their mechanical behavior.Main aspects are the highly nonlinear stress-strain behavior accompanied with hysteresis effects, permanent set and the Mullins effect.Beyond that, viscoelastic properties are observed as well.This is expressed, for example, in relaxation and a certain degree of rate dependency.Concerning rate dependency, a distinction has to be made between different specific classes of rubber-like materials.For example, thermoplastic elastomers show both relaxation and rate dependency in a significant degree [1].Filled rubbers, however, only show a weakly pronounced rate dependency over a wide frequency range, while still having significant relaxation properties [2][3][4].
Due to its complex behavior, constitutive modeling of rubber-like materials is a challenging task.Indeed, many different modeling approaches with different levels of complexity exist, where models are chosen in such a way that the main effects of the application to be analyzed are captured.One of those research directions is the study of viscoelastic properties of rubber-like materials, which have been investigated already for many decades [5,6].However, due to the complex phenomena, research on rubber viscoelasticity is still ongoing, which can, for example, be seen by the works [2,3,[7][8][9][10][11].Thereby, one basic procedure is to use an appropriate model to capture time-independent effects and to extend this model by several Maxwell elements [7,10].In other approaches, the viscoelastic properties are not only described by constant relaxation parameters, but by loading history-dependent functions.One promising approach is the application of stress-dependent relaxation times [2,12].Another approach is to formulate an ordinary differential equation (ODE) for the evolution of the relaxation time, where the evolution depends on the complete strain loading history [3,13].In the literature [3], a basic study was conducted on the time-dependent effects of filled rubber (significant relaxation, low rate dependency).The material behavior was explained by a so-called time-rescaling invariant property, which means that the same stress response is obtained in different strain loading programs, if these strain loading programs are only scaled in time.
In this study, the so-called model of rubber phenomenology (MORPH) [14] is regarded in detail.This model is especially suitable for large strain analyses, where inelastic effects play a crucial role, which has already been proven in several applications (e.g.[15,16]).However, in its current state it is not possible to simulate any time-dependent effects.The time-independent MORPH model and its extension to represent time-dependent phenomena is presented in Section 2.Moreover, a short note on time discretization for numerical simulation is given in Section 2.4.Finally, in Section 3 the extended model is applied to different materials: a filled rubber material [2] and a polyurethane based adhesive [17] to show the capability of this approach to capture the basic time-dependent properties of rubber-like materials.

Notation
The MORPH model and its extensions are formulated in the framework of large deformations.Thereby, the local deformation is described by the deformation gradient  =   ̄ ⊗ ̄ , with   being its coefficients with respect to a Cartesian basis vector system ̄ = { ̄ , ̄ , ̄ }.The underscores denote the order of tensorial quantities.Based on this, the local volume ratio , the left Cauchy-Green tensor , and the isochoric part of the left Cauchy-Green tensor b are defined by where det[ ] is the determinant and ( )  stands for the transpose of a second order tensor.Moreover, the velocity gradient , the deformation rate tensor  and the spin tensor  are defined by Here, • ( ) represents the Lagrangian (material) time derivative, and the symbols ( ) −1 , ( )  and ( )  denote the inverse, the symmetric part and the antisymmetric part, respectively, of a second order tensor.The Zaremba-Jaumann rate of a second order tensor is given by *  = •  − ⋅  +  ⋅  . ( With this definition, the evolution of b can be formulated by the Here,  ′ =  − This norm is applied to define the following scalar equivalent measures at time point : Thereby,   () is employed to represent a scalar measure of the absolute value of the current deformation, and *   () represents a scalar measure of the absolute value of the current deformation rate.

Time-independent MORPH model
The MORPH model is described by the following set of constitutive equations [16].First, the total Kirchhoff stress tensor  (also known as weighted Cauchy stress tensor) is defined by a sum of three different terms as follows Therein, the first term is a Neo-Hookean like stress contribution with a loading history-dependent material function (   ) (see explanation below), and    () is the maximum value of the equivalent deformation measure   () within the time span 0 ≤  ≤  of a considered deformation process: The second term in Equation ( 7) is the so-called additional stress   , and the last term is a volumetric contribution with the bulk modulus  to capture nearly incompressible behavior.The additional stress   is governed by the following ODE where (   ) is a second loading history-dependent material function, and   is the limiting stress defined by Here, exp(⋅) is the exponential tensor function and (   ) is a third loading history-dependent material function.The loading history-dependent material functions appearing in Equations ( 7), ( 9) and ( 10) are defined as follows: Altogether, the MORPH model contains nine material parameters: parameters  1 −  8 appearing in Equations ( 10) and ( 11) as well as the bulk modulus  appearing in Equation ( 7).

Extension of the MORPH model by time-dependent phenomena
One typical approach to extend any time-independent model by time dependency is to add further stress contributions that describe viscoelastic behavior.In the case of solid materials, typically one or more Maxwell elements are added in parallel to the original model (also known as Prony series).Thereby, the stress contributions of the original model and all Maxwell elements sum up to a total stress.Moreover, the constitutive equations of the original model remain unchanged.
In the framework of large deformations, there exist numerous approaches to describe the phenomenological behavior of a Maxwell element (see, e.g., [18]).One approach can, for example, be given by the ODE: where  is the shear modulus and   is the relaxation time.Thus, the last term of this ODE is responsible for the relaxation and rate dependency of this model.In this study, another basic approach is applied to incorporate time dependent effects into the MORPH model.Thereby, the already existing set of constitutive equations of the original MORPH model is modified.In particular, the ODE for the evolution of the additional stresses   given in Equation ( 9) is extended by a relaxation term according to Equation (12).Thereby, the relaxation properties are not defined by a constant relaxation time   , but by a scalar material function (…) ≥ 0, which enables the consideration of process-depending relaxation properties.Thus, a general form of the proposed ODE for the additional stress reads: There are different possibilities to define the function (… ).
Thereby, *   is employed to represent the norm of the strain rate and  0 ≥ 0 is the initial value for the ODE.In the following, two different specific functions are regarded.In the first version, an ODE for time-rescaling invariant behavior is adopted from the literature [3].Thereby, the basic structure of the ansatz is exactly recovered, but a minor change concerning the norm of the strain rate is applied, since only a one-dimensional model was regarded in the literature [3].
and   ( = 0) =  0 Here,  1 ,  2 and  3 are material parameters.Moreover, the function (   ) is applied as pre-factor to the relaxation term.The additional factor  0 is introduced to keep the physical unit of  unchanged.In Section 3.2, this version is applied to the polyurethane based adhesive, which shows distinct rate dependency.

Time stepping algorithm
In order to solve the constitutive equations of the MORPH model numerically, an appropriate time stepping algorithm has to be applied.To this end, a representative time step Δ =  +1 −   > 0 is considered.Quantities at the time point   , which are assumed to be completely known, are denoted by  ( ).Quantities of the new time point  +1 are denoted by +1 ( ).Thus, the aim is to calculate the total stress +1  = ( +1 ) according to Equation (7).Scalar and tensorial Lagrangian time derivatives occurring in the constitutive equations are approximated by +1 ḟ( +1 ) ≈ ( +1  −  )∕Δ , +1 •  ≈ It is assumed that the deformation gradient +1  is known at time point  +1 within the local stress calculation algorithm.
Thus, by virtue of the approximation (17) 2 it is possible to directly calculate +1 •  as well as all other kinematic quantities, like +1 , +1 * b, +1  or +1  (according to their definitions given in Section 2.1).From this, also the limiting stress +1   and two terms of Equation ( 7) directly follow.Special regard, instead, has to be applied to the additional stress +1   , for which an ODE has to be solved numerically.To derive its update formula, first the definition of the Zaremba-Jaumann rate ( 3) is substituted into ODE (13).After rearranging terms, one obtains This equation is solved numerically by fixed-point iteration, which converged fast and did not show any stability issues in the current study.If any issues are encountered in future studies, an application of Newton's algorithm can be considered to improve the numerical behavior.
For the calculation of Equation ( 19), the current value +1  = ( +1 *   ,  +1 ) has to be calculated first by solving the corresponding ODE.In the case that the general ODE ( 14) is formulated for the current time point and the approximation ( 17) is substituted, one obtains the general update formula: The derived update formulas for the specific model versions ( 15) and ( 16) are as follows.
F I G U R E 1 Comparison of experiment and simulation for cyclic loading with intermediate relaxation periods for a filled rubber material (left and middle) and identified material parameters (right).Experimental data were acquired from the literature [2].

APPLICATION EXAMPLES
In this section, the proposed modeling approach is applied to two different materials by adjusting the model versions to respective experimental data resulting from uniaxial tension tests.With () being the stretch in loading direction, the gradient of these experiments is defined by: Thereby, the stretch in loading direction () is employed as time resolved input data   = (  ), with  = 1 … , where is the number of time steps.The additional lateral stretch   () occurring in Equation ( 22) is calculated algorithmically for each time step by the help of the stress boundary conditions   () =   () = 0, where   and   are coefficients of the first Piola-Kirchhoff stress tensor  =  −1 ⋅ .Thus, the engineering stress in loading direction   () is calculated and compared with the experimental data.
The stress calculation is embedded in an optimization procedure to adjust the model to the experimental data and thereby identify the included material parameters.The stress calculation and the optimization are implemented in Matlab scripts, where the optimization is based on the Matlab function lsqnonlin.

Filled rubber data
In a first example, uniaxial tension test data of a filled rubber material were acquired from the literature [2].The regarded data includes a displacement driven uniaxial tension test, which consists of periods of cyclic loading with constant loading rate and periods of constant strain (relaxation periods).Within the cyclic loading periods, different maximum strains are approached.For each maximum strain value, loading/unloading is repeated four times.Within the simulation, the timerescaling invariant approach I according to Equation ( 15) is used.Therein, the bulk modulus and the initial value for the function  were set to  = 2000 MPa and  0 = 0, respectively.The results of the parameter identification are depicted in Figure 1 (including a list of the identified material parameters).Therein, the stresses of experiment and simulation are plotted against the time and the stretch, respectively.It can be observed that the basic nonlinear stress-strain behavior can be captured very well, which already is a feature of the original MORPH model.However, with the proposed extension, also relaxation can be captured, which can be seen in the relaxation periods in Figure 1 (e.g.periods around  = 60 s or  = 80 s).In these periods, the original MORPH model would show constant stress values.

Polyurethane based adhesive
In a second application example, the polyurethane based adhesive Teroson PU 1510 [17,19] is considered in its fully cured state.For this material, uniaxial tension tests with different loading programs were conducted and the extended MORPH (approach II according to Equation ( 16)) was adjusted to these data.The tests include the following three loading scenarios.
(i) Cyclic loading at different strain rates: Cyclic loading was conducted up to different maximum stretches (1.25 and 1.5) with three repeated load cycles for each maximum stretch level.The cross-head speeds were 10 mm/min and 400 mm/min.The comparison of the stress-strain curves of experiments and simulations is depicted in Figure 2. It can be observed that different strain rates lead to different stress responses, which is captured by the extended MORPH (approach II).(ii) Loading and unloading with intermediate relaxation periods (virgin material): In this loading scenario, first monotonic loading with cross-head speed 100 mm/s was applied with intermediate relaxation periods at stretches 1.2, 1.4 and 1.6.Subsequently, monotonic unloading with the same cross-head speed and additional relaxation periods at stretches 1.4 and 1.2 was conducted.This scenario was applied to a virgin material (i.e., a material without any pre-loading history).The comparison of the experiment and the according simulation is shown in Figure 3.It can be observed that relaxation can be captured very well, both in loading and unloading of a virgin material.(iii) Loading and unloading with intermediate relaxation periods (pre-conditioned material): In this loading scenario, first a cyclic loading with increasing maximum stretches (1.2, 1.4 and 1.6) was applied (similar to loading scenario (i)) to bring the material to a pre-conditioned state.Subsequently, monotonic loading with intermediate relaxation periods according to loading scenario (ii) was applied.Thereby, the relaxation behavior in a stationary loading cycle can be analyzed.The corresponding comparison between experiment and simulation is depicted in Figure 4. Here, it can be seen that the relaxation behavior is also captured very well for the case of a stationary loading cycle (i.e., in a pre-conditioned material).

I G U R E 4
Monotonic loading with intermediate relaxation periods (scenario (iii)) for a pre-conditioned material -comparison of experiment and simulation.
Note that the data of all load scenarios were used in one optimization.Thus, all corresponding simulations were done with the same set of material parameters.Thereby,  0 = 9.5 ⋅ 10 −3 s −1 and  = 8000 MPa were pre-set.All other identified parameters are listed in Equation (23).

CONCLUSION
In this study, an approach is proposed to extend the MORPH model in such a way that time-dependent phenomena can be captured.In a first step, the original constitutive equations of the MORPH model were extended by a relaxation term.Thus, the system of tensorial equations is only slightly modified and its size remains unchanged.Additionally, two different specific approaches were applied: a time rescaling invariant approach I, which is especially suitable to simulate the timerescaling invariant behavior of filled rubber materials, and a rate-dependent approach II, which shows a good fitting quality when applied to a polyurethane based adhesive.In these cases, time-dependent phenomena could be included by extending the model by only two (approach I) and three (approach II) additional material parameters, respectively.In future work, more detailed analyses of the approaches should be conducted.This especially includes the validation on the basis of more experimental data and the finite element implementation of the extended MORPH model to enable the simulations of real components.

A C K N O W L E D G M E N T S
Open access funding enabled and organized by Projekt DEAL.

R E F E R E N C E S
If, for example, a constant relaxation time should be used, then this function is just a constant factor reading  = 1∕  .In the case that stress-dependent relaxation properties should be included, the function could be formulated as an algebraic equation as function of the current stress norm, that is () = (||()||).If the relaxation properties depend on the complete strain loading history, then the function must be governed by an integral equation or an ODE.A general form of such an ODE reads Ṁ =   ( *   , ) with ( = 0) =  0 .

F I G U R E 2
Cyclic loading at different strain rates (experiments vs. simulations).FI G U R E 3Monotonic loading with intermediate relaxation periods (scenario (ii)) for a virgin material -comparison of experiment and simulation.
Thus, the one-dimensional strain rate is replaced by *  1 and  2 are material parameters.This version is applied to a filled rubber material in Section 3.1.In a second model version, another ODE for the function  is chosen to capture strain rate effects as well.