A Flory’s only framework for rate-dependent viscoelastoplasticity at large strains

Abstract

Flory's multiplicative decomposition of the deformation gradient is a powerful tool for developing analytical and computational solutions for hyperelastic materials. This approach simplifies the analysis by separating the deformation into distinct components, isochoric and volumetric. This allows for material models that follow the growth condition, preventing unrealistic material self-intersection. Despite being widely used in elasticity, its use as a basis for theoretical developments in plasticity and viscosity is associated with the Kröner–Lee multiplicative decomposition. As a novelty, in this study it is observed that the Flory’s decomposition can be interpreted as the generation of isochoric and volumetric spaces. From this space identification, an additive decomposition of strains can be used, resulting in an alternative framework for the computational modeling of viscoelastoplastic materials. In this framework, the plastic flow direction and the viscous strain rate are taken inside the isochoric Flory’s spaces, simplifying the resulting constitutive models. In particular, a viscoelastoplastic model with isotropic and kinematic hardenings and viscous parameters dependent on the rate and intensity of strains is developed. To demonstrate the effectiveness of the proposed framework, the author developed a specific model suitable for computer simulations using the finite element method. The proposed approach is then validated by comparing representative examples—simulated with an in-house code—to literature experimental data.

Introduction

The behavior of many solids depends on how fast they deform, not just on the amount of deformation. This "time rate dependence" becomes crucial even without considering acceleration. Therefore, static analysis can be performed when both the material's response is time-independent and inertial forces are negligible.

A material is considered viscoelastic when a dissipative phenomenon takes place at any level of stress and rate dependence is present; on the other hand, a material is considered elastoplastic if there is a yielding stress under which elastic deformations take place and there is not rate dependence.

A material can possess both viscoelastic and elastoplastic characteristics simultaneously, which is the focus of this study. It can be in the viscoelastic regime for low stress levels (under yielding) or in the viscoelastoplastic regime for high stress levels (over yielding). The present study proposes a new viscoelastoplastic framework to model materials developing large strains including variable kinematic and isotropic hardenings and viscosity parameters dependent on the level of strains and strains’ rate.

To start the discussion on modeling large-strain behavior of viscoelastic, viscoplastic and viscoelastoplastic materials, it is worth noting the valuable reviews provided by [1] and a previous study of the author [2]. These references offer a comprehensive overview of the development in this field summarized as follows. That the framework given by Perzina [3] and Duvaut–Lions [4] fits into the overstress models category and can represent viscoelastoplastic materials. In this case, following a good description given by [1], these models do not use typical Kuhn–Tucker equations to solve the plastic multiplier (or equivalent plastic strain), but calculate a temporal rate of equivalent (visco-) plastic strain through a direct equation that relates overstress and viscosity. These models are widely used in studies involving small strains [5, 6–7]. The Perzyna's strategy was also extended to model large strains by references [8, 9–10]. As it is considered an alternative formulation for large strains, it is important to mention (following the description given by reference [1]) that Perzyna's strategy is not extendable to the inviscid case (elastoplasticity), resulting in unstable algorithms in quasi-inviscid situations.

However, references [1] and the review carried out by the author in [2] do not indicate the existence of a preferred thermodynamically consistent formulation for treating viscoplasticity. For example, reference [1] only indicates the so-called consistency models introduced by [11] and explored by other authors such as [12, 13–14]. This strategy introduces stresses’ time dependence directly on the yielding surface, transferring viscosity to the evolution of the equivalent plastic strain. By the other hand, reference [2] indicates that, when using the Kröner–Lee multiplicative decomposition [15, 16–17], the thermodynamically consistent elastoplasticity—separated from viscoelasticity—is already a well-established subject.

One can say that the development of thermodynamically consistent formulations that combine viscoelasticity and viscoplasticity resulting in a unified viscoelastoplastic framework is still an open and broad field, taking as a starting point the works of [10, 18]. In this strategy, the Kröner–Lee [19, 20] multiplicative decomposition can be done in parallel (leading to multiple visco/elastic and elasto/plastic behaviors) or in series (leading to elastic/viscous/plastic behavior). It is opportune to mention that this framework is quite complicated from the algebraic and interpretative point of view and a considerable number of studies also combine the Flory’s decomposition with the Kröner–Lee one [21, 22, 23–24].

As an alternative to the consistency viscoplastic model, reference [1] proposed the concept of continuum elastic corrector rate to deal with the viscoelastoplastic problem at small strains, indicating that its extension to finite strains is possible, if the elastoplastic strategy presented by [25] is followed. The constitutive model proposed by [1] is based on the maximum dissipation principle and concentrates the viscous effect on the evolution of the equivalent plastic strain, imposing viscosity on the consistency condition (Kuhn–Tucker). As reference [1] considers viscosity in the plastic portion of the model, it does not directly consider the possibility of developing viscous stress in the elastic phase. Using the viscoelastoplastic serial arrangement model proposed by [18], reference [2] develops a model that considers viscosity in both elastic and plastic parts. The model was verified using experimental results for polytetrafluoroethylene presented by [26].

Another alternative elastic–viscoplastic framework follows the decomposition proposed in by Brünig [27] based on the mixed-invariant metric [28, 29], making a multiplicative decomposition in logarithm strains measures. This formulation is applied to solve elastic–viscoplasticity in metals considering hydrostatic stresses’ sensitivity.

As large-strain viscoelastoplasticity directly involves large-strain elastoplasticity, some comments regarding this subject must be given. As described by [30], elastoplastic finite strain continuous models can be divided into two large groups; the first—in its initial stage—uses the additive decomposition of strain rates at the current configuration (Eulerian space) in elastic and plastic parts, and searches for coherent constitutive laws to calculate objective stress rates. This framework is usually called rate-type model [31] and, sometimes, hypoelastic model. The last denomination is applied when the elastic range is considered infinitesimal [30], but results nonobjective for large rotation applications. A large amount of works are dedicated to this strategy, since pioneer studies as [32], through important contributions that improve the performance of the technique, as for example, [33, 34–35]. In this group, one may include the self-consistent Eulerian formulation [36], in which logarithm strains are used, rendering significant improvements. A recent study [37] details the rate-type models of elastoplastic materials, it explores—in a strong mathematical language—most aspects of this framework and includes the proposition of a rate-type model that does not use any strain decomposition.

The second group of elastoplastic finite strain continuous models is called hyperelastic-based plasticity [30]. Models of this group are developed applying the multiplicative decomposition of Kröner–Lee [19, 20], which splits the gradient deformation into plastic and elastic parts, defining an intermediate space where only elastic strains are developed, usually accompanied by residual stresses [38]. According to [30], the intermediate configuration is valid in a point-wise sense, i.e., in general, compatible unstressed configurations do not exist. The elastic part is modeled by hyperelastic potentials that guaranty objectivity in large rotations and strains. As mentioned by [39], this approach becomes widespread after the large computational evolution established by [40] and the volume conserving plastic flow proposed by [41, 42], adapted in a consistent way by [16].

Although some studies related with plasticity and viscoplasticity uses Flory’s decomposition [24, 26] to establish isochoric flows, they consider this decomposition linked with the Kröner–Lee framework. As a novelty, in the present study a new viscoelastoplastic framework (including a specific model) for isotropic materials (based only on the Flory’s decomposition) is presented. The proposed model allows additive decomposition of elastic, plastic and viscous components and results in very simple numerical expressions. In summary, the proposed model does not use hereditary integrals as the Perzyna-like models, it does not limit the viscous behavior in the plastic multiplier space as the consistency formulations, and it does not use the Kröner–Lee decomposition or the Brünig decomposition. In addition, the proposed model is objective for large rotations and strains.

It is important to mention that in hyperelasticity the use of Flory’s decomposition is usual, see [44, 45], for instance, and that in soft tissue viscoelastic applications, some studies use the Flory’s decomposition to model viscosity, see [46, 47] for instance.

As mentioned, the starting point of the proposed model is the Flory’s multiplicative [43] decomposition, which separates the continuum isochoric and volumetric behaviors. This separation allows the definition of viscous flow and plastic flow directly in the isochoric directions, which allows a direct calculation of the viscous (isochoric) strain rate and cumulated plastic stress. One may cite some previous author’s works [48, 49] that used the isochoric strain directions (achieved from Flory’s decomposition) to solve rate-independent plasticity and viscoelasticity. Here, both plastic and viscous models of [48, 49] are combined to result the proposed viscoelastoplastic formulation, and an Appendix is provided to show the thermodynamic consistency of dissipations.

To present this study, Sect. 2 describes the Flory’s decomposition, the elastic potentials used and what is understood here as isochoric and volumetric strain components (or space). Section 3 shows the Zener-type viscoelastic constitutive model that serves as the basis of the proposed viscoelastoplastic model. Section 4 presents the proposed viscoelastoplastic model, introducing the concept of cumulated plastic stress and updated isochoric plastic direction. Section 5 presents in a briefly way the variational form of the motion equation and a specific finite element strategy presented by the author in [50, 51]. Section 6 presents viscoelastic, elastoplastic and viscoelastoplastic examples in which theoretical benchmarks are used to validate the proposed model and experimental data are used to show its good behavior. Although all variables are described when first appears, Appendix A provides a list of variables. In Appendix B, a simple implementation guide is provided facilitating the implementation of the proposed framework. At the end of the study, conclusions and future perspectives are presented.

Flory’s multiplicative decomposition and basics

As the proposed model is based on the Flory’s [43] decomposition and on the understanding of isochoric directions and/or the definition of isochoric spaces, one starts with the hyperelastic specific strain energy function, from which these concepts naturally arise. It adopts the Green–Lagrange strain; thus, the specific strain energy can be directly written as a function of the right Cauchy–Green stretch tensor $\mathbf{C}={\mathbf{A}}^T·\mathbf{A}$, i.e., $W(\mathbf{C})$, in which $\mathbf{A}$ is the deformation gradient. In order to split the specific strain energy in its isochoric and volumetric parts, one needs to proceed with the Flory’s multiplicative decomposition [43] of $\mathbf{A}$ as follows:

The isochoric part of $\mathbf{A}$, called here $\overline{\mathbf{A}}$, is given by

1

$$\overline{\mathbf{A}}=J^{-1/3}\mathbf{A}\text{with}Det\left(\overline{\mathbf{A}}\right)=Det\left(\mathbf{A}\right){\left(J^{-1/3}\right)}^3=1$$and its volumetric part ($\widehat{\mathbf{A}}$) by

2

$$\widehat{\mathbf{A}}=J^{1/3}\mathbf{I}\text{with}Det\left(\widehat{\mathbf{A}}\right)=Det\left(\mathbf{A}\right){\left(J^{1/3}\right)}^3=J$$in which $J=Det(\mathbf{A})$, and $\mathbf{I}$ is the identity tensor of second order.

Using Eqs. (1) and (2), one recovers the deformation gradient by the product

3

$$\mathbf{A}=\widehat{\mathbf{A}}·\overline{\mathbf{A}}=\overline{\mathbf{A}}·\widehat{\mathbf{A}}$$in which “$·$” is a simple contraction.

From Eqs. (1) and (2), one writes the decomposition of the right Cauchy–Green stretch tensor as:

4

$$\overline{\mathbf{C}}={\overline{\mathbf{A}}}^T·\overline{\mathbf{A}}=J^{-2/3}\mathbf{C}$$

5

$$\widehat{\mathbf{C}}={\widehat{\mathbf{A}}}^T·\widehat{\mathbf{A}}=J^{2/3}\mathbf{I}$$

Identifying the Flory’s decomposition of the right Cauchy–Green stretch tensor by the following operations:

6

$$\mathbf{C}=\widehat{\mathbf{C}}·\overline{\mathbf{C}}=\overline{\mathbf{C}}·\widehat{\mathbf{C}}$$

7

$$Det\left(\overline{\mathbf{C}}\right)={\left(J^{-2/3}\right)}^3{J}^2=1Det\left(\widehat{\mathbf{C}}\right)={\left(J^{2/3}\right)}^3=J^2$$being $\overline{\mathbf{C}}$ isochoric and $\widehat{\mathbf{C}}$ volumetric. Using this decomposition, the strain energy density (scalar) can be written as:

8

$$\psi ={\psi }^{\mathit{vol}}(\widehat{\mathbf{C}})+{\psi }^{\mathit{iso}}(\overline{\mathbf{C}})\text{or}\psi ={\psi }^{\mathit{vol}}(J)+{\psi }^{\mathit{iso}}(\overline{\mathbf{C}})$$in which $\psi$ stands for strain energy density.

Using the Rivlin–Saunders [52] isochoric strain energy and the Hartmann–Neff [53] volumetric component, Eq. (8) turns into:

9

$$\psi ={\psi }^{\mathit{vol}}(J)+{\psi }^{iso1}({\overline{I}}_1)+{\psi }^{iso2}({\overline{I}}_2)$$in which

10

$${\psi }^{\mathit{vol}}(J)=\frac{K}{8n^2}\left(J^{2n}+J^{-2n}-2\right),{\psi }^{iso1}=\frac{G_1}{4}({\overline{I}}_1-3),{\psi }^{iso2}=\frac{G_2}{4}({\overline{I}}_2-3)$$being ${\overline{I}}_1$ and ${\overline{I}}_2$ the first and second invariants of the isochoric Cauchy–Green stretch tensor, Eq. (4), $K$ is the bulk modulus, $G_1=G_2=G$ is the elastic shear modulus, and $n$ is a growth parameter that keeps the convexity of the energy density. The upper index $\mathit{vol}$ stands for volumetric, $iso1$ for the first isochoric invariant and $iso2$ for the second isochoric invariant.

Using the Green–Lagrange $\mathbf{E}$ strain, given by:

11

$$\mathbf{E}=\frac{1}{2}(\mathbf{C}-\mathbf{I})=\frac{1}{2}({\mathbf{A}}^T·\mathbf{A}-\mathbf{I})$$as strain measure and, from Eq. (10), one calculates the elastic second Piola–Kirchhoff stress as:

12

$$\mathbf{S}=\frac{\partial \psi }{\partial \mathbf{E}}=\frac{\partial {\psi }^{\mathit{vol}}}{\partial \mathbf{E}}+\frac{\partial {\psi }^{iso1}}{\partial \mathbf{E}}+\frac{\partial {\psi }^{iso2}}{\partial \mathbf{E}}\text{or in compact form}\mathbf{S}={\mathbf{S}}^{\mathit{vol}}+{\mathbf{S}}^1+{\mathbf{S}}^2$$with

13

$${\mathbf{S}}^{\mathit{vol}}=\frac{\partial {\psi }^{\mathit{vol}}}{\partial \mathbf{E}}=\left\{\frac{K}{4n}\left(J^{2n-1}-J^{-(2n+1)}\right)\right\}{\mathfrak{E}}^{\mathit{vol}}\text{and}{\mathfrak{E}}^{\mathit{vol}}=\frac{\partial J}{\partial \mathbf{E}}=J{\mathbf{C}}^{-1}$$

14

$${\mathbf{S}}^1=\frac{G_1}{4}{\mathfrak{E}}^1=\frac{G}{4}{\mathfrak{E}}^1\text{and}{\mathfrak{E}}^1=\frac{\partial {\overline{I}}_1}{\partial \mathbf{E}}=-\frac{2}{3}{J}^{-2/3}Tr(\mathit{C}){\mathit{C}}^{-1}+2J^{-2/3}\mathit{I}$$

15

$${\mathbf{S}}^2=\frac{G_2}{4}{\mathfrak{E}}^2=\frac{G}{4}{\mathfrak{E}}^2\text{and}{\mathfrak{E}}^2=\frac{\partial {\overline{I}}_2}{\partial \mathbf{E}}=2J^{-4/3}\left(-\frac{2}{3}{\mathit{C}}^{-1}{I}_2+\left\{Tr(\mathit{C})\mathit{I}-{\mathit{C}}^t\right\}\right)$$

In Eqs. (14) and (15), tensors ${\mathfrak{E}}^1$ and ${\mathfrak{E}}^2$ are defined here as isochoric strains defining a non-unitary isochoric space base, because when applying the pushing forward [44, 45], the second Piola–Kirchhoff-like stresses ${\mathbf{S}}^1$ and ${\mathbf{S}}^2$ correspond to the deviatory Cauchy stress components, see references [48, 50] for a straightforward operational demonstration. An isochoric unitary space base arises when one uses tensors ${\mathfrak{E}}^1/\sqrt{{\mathfrak{E}}^1:{\mathfrak{E}}^1}$ and ${\mathfrak{E}}^2/\sqrt{{\mathfrak{E}}^2:{\mathfrak{E}}^2}$ to assemble other quantities. In addition, ${\mathbf{C}}^{-1}$ is the inverse of the Cauchy stretch tensor, $I_1$ is the first invariant of the Cauchy stretch and $I_2$ its second invariant.

Zener-like viscoelastic model: isochoric strain rate components:

Before introducing the proposed viscoelastic behavior, it is important to mention that after establishing the strain space separation in volumetric and isochoric components, viscosity (and plasticity) will develop only regarding isochoric components, i.e., the volumetric behavior is considered pure elastic. Using this separation, in Appendix D, it is shown the thermodynamic consistence of the proposed viscous strategy.

Two ways are followed to write the isochoric viscoelastic part of the proposed model using the isochoric directions [49]. One based directly on the isochoric strain rates ${\dot{\mathfrak{E}}}^i$ and other based on the time rate of the first and second invariants ${\dot{\overline{I}}}_i$ of the isochoric Cauchy–Green stretch tensor. In Fig. 1b, one can see the Zener-like model in which the invariants (scalars) are additively split. Additionally, thanks to the multiplicative Flory’s decomposition, the isochoric strains rates always fall in the same isochoric space. Thus differently of all consulted works that use Kröner–Lee decomposition, one can also split the isochoric strain rates in an additive way, see the schematic representation of Fig. 1a.

[See PDF for image]

Fig. 1

Zener-like isochoric viscoelastic models

Zener-like model: isochoric strain rate components

From figure Fig. 1a, one can state that the total isochoric Piola–Kirchhoff stress ${\mathit{S}}_T^i$ (following direction $i$) can be written as:

16

$${\mathbf{S}}_T^i=\mathbf{S}_{}^i+{\mathbf{S}}_v^i$$in which $\mathbf{S}_{}^i$ is the elastic part of the deviatory second Piola–Kirchhoff-like stress and ${\mathbf{S}}_v^i$ the viscous part, given by:

17

$$\mathbf{S}_{}^i=\frac{G}{4}\mathfrak{E}_{}^i{\mathbf{S}}_v^i=\frac{\overline{G}}{4}{\mathfrak{E}}_e^i=\frac{\eta (\bullet )}{4}{\dot{\mathfrak{E}}}_v^i$$in which the viscous parameter $\eta (\bullet )$ may depend upon the isochoric Green stretch tensor or upon the isochoric strain rate (or both). In Eq. (17), $\mathfrak{E}_{}^i$ is the total isochoric strain, ${\mathfrak{E}}_e^i$ is the elastic isochoric strain at the viscoelastic branch of the model, ${\dot{\mathfrak{E}}}_v^i$ is the time rate of the viscous part of isochoric strain at the viscoelastic branch of the model, and $\overline{G}$ is the shear elastic modulus at the viscoelastic branch of the model. To be concise, in the following equations the sign $(\bullet )$ is suppressed. From Eqs. (14) and (15), one can see that the isochoric strain rates falls in the same isochoric original space. It means that the time derivative of isochoric strains is also isochoric. Thus, remembering that upper dot means time derivative, thanks to the Flory’s decomposition, being ${\dot{\mathfrak{E}}}_e^i$ and ${\dot{\mathfrak{E}}}_v^i$ at the same constrained direction, one considers that:

18

$$\dot{\mathfrak{E}}_{}^i={\dot{\mathfrak{E}}}_e^i+{\dot{\mathfrak{E}}}_v^i$$

Using the second part of Eq. (17), one writes:

19

$${\dot{\mathfrak{E}}}_v^i=\frac{\overline{G}}{\eta }{\mathfrak{E}}_e^i$$

Substituting Eq. (19) into Eq. (18) results:

20

$$\dot{\mathfrak{E}}_{}^i={\dot{\mathfrak{E}}}_e^i+\frac{\overline{G}}{\eta }{\mathfrak{E}}_e^i$$

For simplicity, in this study the backward finite difference is adopted in Eq. (20)—initially over ${\dot{\mathfrak{E}}}_e^i$—achieving:

21

$${\left({\mathfrak{E}}_e^i\right)}_{s+1}=\left({\left({\dot{\mathfrak{E}}}^i\right)}_{s+1}+\frac{{\left({\mathfrak{E}}_e^i\right)}_s}{\mathrm{\Delta }t}\right)/\left(\frac{\overline{G}}{\eta }+\frac{1}{\mathrm{\Delta }t}\right)$$in which $s+1$ means current instant and $\mathrm{\Delta }t$ is a time step.

However, it is important to note that Eq. (20), after some mathematical steps, could be integrated in time by other more elaborated strategies as, for example, the exponential mapping, see references [16, 24, 54, 55], for instance.

Using (16) and (17), for the current instant $s+1$ results:

22

$${\left({\mathbf{S}}_T^i\right)}_{s+1}=\frac{G}{4}\left(\mathfrak{E}_{}^i\right)_{s+1}^{}+\frac{\overline{G}}{4}{\left({\mathfrak{E}}_e^i\right)}_{s+1}$$and using Eq. (21), one writes the proposed numerical equation for the total isochoric part of the second Piola–Kirchhoff stress as:

23

$${\left({\mathbf{S}}_T^i\right)}_{s+1}=\left(\mathbf{S}_{}^i+{\mathbf{S}}_v^i\right)_{s+1}^{}=\frac{G}{4}\left(\mathfrak{E}_{}^i\right)_{s+1}^{}+\frac{\overline{G}}{4}\left({\left({\dot{\mathfrak{E}}}^i\right)}_{s+1}+\frac{{\left({\mathfrak{E}}_e^i\right)}_s}{\mathrm{\Delta }t}\right)/\left(\frac{\overline{G}}{\eta }+\frac{1}{\mathrm{\Delta }t}\right)$$in which ${\left({\dot{\mathfrak{E}}}^i\right)}_{s+1}$ is also calculated by finite difference and ${\left({\mathfrak{E}}_e^i\right)}_s$ is the result of a previous time step.

Regarding time discretization, example 1 shows that volumetric and isochoric strains are independent, despite using ${\mathfrak{E}}_s^i$ and ${\left({\dot{\mathfrak{E}}}_e^i\right)}_s$ at time $t_{s+1}$. In Sect. 5.2 (solution process), the dependence of viscosity regarding strain rate and strain level is considered and some examples explore the suggested ones.

Zener-like model: strain invariants’ rate

Although Eq. (23) gives very good results for practical applications, this item shows an alternative version in which viscosity depends on scalar isochoric invariants’ rate, see Fig. 1b. The deviatoric viscous Piola–Kirchhoff stress component ${\mathbf{S}}_v^i$, Fig. 1b, becomes:

24

$${\mathbf{S}}_v^i=\frac{\overline{\eta }(\bullet )}{4}{\dot{\overline{I}}}_{\mathit{iv}}{\mathfrak{E}}^i\text{or}{\mathbf{S}}_v^i=\frac{\overline{G}}{4}\frac{\partial \overline{I}_{\mathit{ie}}^{}}{\partial \mathbf{E}}$$in which $\overline{\eta }$ is the viscous parameter of this model, and ${\overline{I}}_{\mathit{iv}}$ and ${\overline{I}}_{\mathit{ie}}$ are parts of the total isochoric strain invariant ${\overline{I}}_i$. Being the strain invariant a scalar value, it is fair to define its parts (not necessarily with real meaning) as:

25

$${\overline{I}}_i={\overline{I}}_{\mathit{ie}}+{\overline{I}}_{\mathit{iv}}=\alpha {\overline{I}}_i+(1-\alpha ){\overline{I}}_i$$in which ${\overline{I}}_{\mathit{ie}}=\alpha {\overline{I}}_i$ and ${\overline{I}}_{\mathit{iv}}=(1-\alpha ){\overline{I}}_i$. Scalar $\alpha$ is used to simplify understanding.

Using Eq. (25) into Eq. (24) results:

26

$${\mathbf{S}}_v^i=(1-\alpha )\frac{\overline{\eta }}{4}{\dot{\overline{I}}}_i{\mathfrak{E}}^i\text{and}{\mathbf{S}}_v^i=\frac{\overline{G}}{4}\alpha {\mathfrak{E}}^i$$

Thus, considering that both values are equal:

27

$$\alpha =\overline{\eta }{\dot{\overline{I}}}_i\left./\phantom{\overline{\eta }{\dot{\overline{I}}}_i\left(\overline{G}+\overline{\eta }{\dot{\overline{I}}}_i\right)}\right)\left(\overline{G}+\overline{\eta }{\dot{\overline{I}}}_i\right)$$

Returning $\alpha$ into the second expression of Eq. (26) results

28

$${\mathbf{S}}_v^i=\frac{\overline{G}}{4}\frac{\overline{\eta }{\dot{\overline{I}}}_i}{\left(\overline{G}+\overline{\eta }{\dot{\overline{I}}}_i\right)}{\mathfrak{E}}^i$$

The total isochoric Piola–Kirchhoff stress is achieved introducing Eqs. (28) and the first of (17) into Eq. (16), i.e.,

29

$${\mathbf{S}}_T^i=\frac{1}{4}\left[G+\frac{\overline{\eta }{\dot{\overline{I}}}_i}{\left(1+\left(\overline{\eta }{\dot{\overline{I}}}_i/\overline{G}\right)\right)}\right]{\mathfrak{E}}^i$$in which one also used finite difference to calculate ${\dot{\overline{I}}}_i$.

It is worth noting that, in this case, for a simple relaxation test one has ${\dot{\overline{I}}}_i=0$, thus it is necessary to calculate the viscous stress from:

30

$${\mathbf{S}}_v^i=\frac{\overline{\eta }}{4}{\dot{\overline{I}}}_{\mathit{vi}}{\mathfrak{E}}^i=\frac{\overline{G}}{4}\alpha {\mathfrak{E}}^i\text{or}\alpha =\frac{\overline{\eta }}{\overline{G}}{\dot{\overline{I}}}_{\mathit{iv}}$$

Introducing $\alpha$ into Eq. (25) in which one keeps $I_{\mathit{iv}}$ as an independent unknown, results

31

$${\overline{I}}_i=\frac{\overline{\eta }}{\overline{G}}{\dot{\overline{I}}}_{\mathit{iv}}{\overline{I}}_i+{\overline{I}}_{\mathit{iv}}$$and applying finite difference one writes:

32

$${\overline{I}}_{\mathit{iv}}^{s+1}={\overline{I}}_i\left(1+(\overline{\eta }/\overline{G}){\overline{I}}_{\mathit{iv}}^s/\mathrm{\Delta }t\right)/\left(1+(\overline{\eta }/\overline{G})\overline{I}_i^{}/\mathrm{\Delta }t\right)$$

From this result, one calculates the time rate ${\dot{\overline{I}}}_{\mathit{iv}}^{s+1}=({\overline{I}}_{\mathit{iv}}^{s+1}-{\overline{I}}_{\mathit{iv}}^s)/\mathrm{\Delta }t$ and introduces in the first of Eq. (30) to compute ${\mathbf{S}}_v^i$. At the beginning of an instantaneous imposed strain ($s=0$), one has ${\overline{I}}_{\mathit{iv}}^0=(1-{\alpha }_0){\overline{I}}_i$, being $\overline{I}_i^{}$ known and ${\alpha }_0$ calculated by Eq. (27) as ${\dot{\overline{I}}}_v\ne 0$ at instant zero.

Viscoelastoplastic models

The proposed viscoelastoplastic models are achieved simply by adding a plastic “brake” with a hardening spring at the elastic branch of the viscoelastic models, i.e., the elastic stress ${\mathbf{S}}^i$ calculated for both viscoelastic models by the first of Eqs. (17) passes to be calculated as a rate-independent (consistent) elastoplastic stress, see Fig. 2a, b for a schematic representation. The plastic framework used here has been recently developed by the author in [48] and is summarized here for clarity.

[See PDF for image]

Fig. 2

Zener-like isochoric viscoelastoplastic models

In Fig. 2, $H$ represents isotropic or cinematic hardening and $\overline{\tau }$ the shear elastic limit, to be explained in what follows.

Defining the elastic limits: ductile material in Flory’s framework

As usual the von-Mises strength criterion,

33

$$J_2<\frac{{\mathit{\sigma }}_y^2}{3}\text{or}\sqrt{\frac{3}{2}{\mathit{\sigma }}^{\mathit{des}}:{\mathit{\sigma }}^{\mathit{des}}}<{\sigma }_y=2\overline{\tau }$$is used to limit the elastic behavior. In Eq. (33), $J_2=\left({\mathit{\sigma }}^{\mathit{des}}:{\mathit{\sigma }}^{\mathit{des}}\right)/2$ is the second invariant of the deviatory Cauchy stress ${\mathit{\sigma }}^{\mathit{des}}$ and $\overline{\tau }$ is the yielding shear stress limit usually considered equivalent to half of the uniaxial elastic limit ${\sigma }_y$. Observing the first of Eqs. (17), the comments after Eq. (10) that defines $G_1=G_2=G$ and Eq. (12), one may split the von-Mises yielding criterion into two Lagrangian components, as follows:

34

$$\frac{3}{2}{\mathbf{S}}^1:{\mathbf{S}}^1-{\overline{\tau }}^2<0\text{and}\frac{3}{2}{\mathbf{S}}^2:{\mathbf{S}}^2-{\overline{\tau }}^2<0$$

In which ${\overline{\tau }}_1={\overline{\tau }}_2=\overline{\tau }$; however in general cases, when $G_1\ne G_2$ the stress limits can be calibrated. Using Eqs. (14) and (15), one writes Eq. (34) in a unified manner as:

35

$$\frac{3}{2}\frac{G^2}{16}{\mathfrak{E}}^i:{\mathfrak{E}}^i-{\overline{\tau }}^2<0$$in which $i$ represents the Lagrangian isochoric strain component (no summation is implied). Equation (35) is the yielding limit before any plastic evolution occurs. It is interesting to note that from Eqs. (35) and (36), two Lagrangian plastic flow directions arise instead of a single Eulerian direction established when applying associative (or nonassociative) flow rules on Eq. (33).

Alternative elastoplastic definitions

The assumed elastoplastic definitions are given in details by the author in [48] and here the necessary ones to present the plastic evolution algorithm are presented. In order to simplify expressions, the isochoric direction symbol $i$ is omitted, see Eq. (35), remembering that all equations are valid for both isochoric directions. The first and more important definition is the cumulated plastic stress ${\mathbf{S}}^p$ that should always be in the isochoric direction, i.e.,

36

$${\mathbf{S}}^p=\frac{G}{4}{\lambda }^{\zeta }\frac{\mathfrak{E}}{\sqrt{\mathfrak{E}:\mathfrak{E}}}$$in which, ${\lambda }^{\zeta }$ is defined as the equivalent scalar plastic strain and $\mathfrak{E}/\sqrt{\mathfrak{E}:\mathfrak{E}}$ is the isochoric unitary direction, the same presented in Eqs. (14) and (15), used here without indices, as the operations that follow are similar for both isochoric directions. The evolution of ${\lambda }^{\zeta }$ is given by:

37

$${\lambda }^{\zeta }={\left({\lambda }^{\zeta }\right)}^{\mathit{ac}}+\zeta \mathrm{\Delta }\lambda$$in which the superscript $_{}^{\mathit{ac}}$ represents cumulated, $\zeta$ is the sign of the kinematic plastic evolution defined by Eq. (40), and $\mathrm{\Delta }\lambda$ is the plastic multiplier.

The kinematic hardening $H^c$ is assumed constant at iterations of the numerical solution process, but may vary between time steps as a function of ${\lambda }^{\zeta }$. To shorten equations, one defines the dimensionless hardening:

38

$$\beta ({\lambda }^{\zeta })=H^c({\lambda }^{\zeta })/(G/4)$$

Thus, the internal variable that controls the change of the origin of the yielding function due to kinematic hardening becomes:

39

$$q=q^{\mathit{ac}}+\beta ({\lambda }^{\zeta })\zeta \mathrm{\Delta }\lambda$$

From the internal variable $q$, the scalar evolution signal (see Eq. (37)) is defined as:

40

$$\zeta =sign\left(\sqrt{\mathfrak{E}:\mathfrak{E}}-{\lambda }^{\zeta }-q\right)$$

And the back stress tensor results:

41

$$\mathit{\chi }=\frac{G}{4}q\frac{\mathfrak{E}}{\sqrt{\mathfrak{E}:\mathfrak{E}}}$$

Similarly, the isotropic hardening $H^i$ is considered and its dimensionless counterpart is:

42

$$\varpi (\lambda )=H^i(\lambda )/(G/4)$$

The isotropic hardening (constant at iterations) also depends upon the cumulated plastic multiplier $\lambda$, for which evolution is given by:

43

$$\lambda ={\lambda }^{\mathit{ac}}+\mathrm{\Delta }\lambda$$

The isotropic hardening internal variable is defined by:

44

$$\kappa ={\kappa }^{\mathit{ac}}+\varpi (\lambda )\mathrm{\Delta }\lambda$$

Elastoplastic objective function and plastic evolution

Adopting the symbol ${\mathit{S}}^e$ for any ${\mathit{S}}^1$ or ${\mathit{S}}^2$ calculated by Eqs. (14) or (15) and using the definitions given by Eqs. (36) through (44), the von-Mises stress value that should be used in Eqs. (34) is:

45

$${\mathbf{S}}^{\mathit{vM}}={\mathbf{S}}^{\mathit{EP}}-\mathit{\chi }={\mathbf{S}}^e-{\mathbf{S}}^p-\mathit{\chi }\text{or}{\mathbf{S}}^{\mathit{vM}}=\frac{G}{4}\mathfrak{E}-\frac{G}{4}{\lambda }^{\zeta }\frac{\mathfrak{E}}{\sqrt{\mathfrak{E}:\mathfrak{E}}}-\frac{G}{4}q\frac{\mathfrak{E}}{\sqrt{\mathfrak{E}:\mathfrak{E}}}=\frac{G}{4}\left(1-\frac{{\lambda }^{\zeta }}{\sqrt{\mathfrak{E}:\mathfrak{E}}}-\frac{q}{\sqrt{\mathfrak{E}:\mathfrak{E}}}\right)\mathfrak{E}$$in which the superscript $_{}^{\mathit{vM}}$ means von-Mises and $_{}^{\mathit{EP}}$ elastoplastic.

Thus, one calculates:

46

$${\mathbf{S}}^{\mathit{vM}}:{\mathbf{S}}^{\mathit{vM}}=\frac{G^2}{16}{\left(1-\frac{{\lambda }^{\zeta }}{\sqrt{\mathfrak{E}:\mathfrak{E}}}-\frac{q}{\sqrt{\mathfrak{E}:\mathfrak{E}}}\right)}^2\mathfrak{E}:\mathfrak{E}$$

Substituting Eq. (46) into Eq. (35), the yielding (objective) function becomes:

47

$$f=\frac{3}{2}\frac{G^2}{16}{\left(1-\frac{{\lambda }^{\zeta }}{\sqrt{\mathfrak{E}:\mathfrak{E}}}-\frac{q}{\sqrt{\mathfrak{E}:\mathfrak{E}}}\right)}^2\mathfrak{E}:\mathfrak{E}-{\left(\sqrt{\frac{3}{2}}\frac{G}{4}\kappa +\overline{\tau }\right)}^2<0$$in which both growing (isotropic hardening) and translation (kinematic hardening) of the yielding surface are considered.

As usual, when $f<0$ the regime is elastic and there is no plastic evolution. But, when a trial of $f$ is greater than or equal to 0, i.e.,

48

$$f^{\mathit{tr}}=\frac{3}{2}\frac{G^2}{16}{\left(1-\frac{{\lambda }^{ac\zeta }}{\sqrt{{\mathfrak{E}}^{\mathit{tr}}:{\mathfrak{E}}^{\mathit{tr}}}}-\frac{q^{\mathit{ac}}}{\sqrt{{\mathfrak{E}}^{\mathit{tr}}:{\mathfrak{E}}^{\mathit{tr}}}}\right)}^2{\mathfrak{E}}^{\mathit{tr}}:{\mathfrak{E}}^{\mathit{tr}}-{\left(\sqrt{\frac{3}{2}}\frac{G}{4}{\kappa }^{\mathit{ac}}+\overline{\tau }\right)}^2\ge 0$$a plastic evolution takes place and $\mathrm{\Delta }\lambda$ can be calculated by the following equality:

49

$$\frac{3}{2}\frac{G^2}{16}{\left(1-\frac{({\lambda }^{ac\zeta }+q^{\mathit{ac}})}{\sqrt{{\mathfrak{E}}^{\mathit{tr}}:{\mathfrak{E}}^{\mathit{tr}}}}-(1+\beta )\zeta \frac{\mathrm{\Delta }\lambda }{\sqrt{{\mathfrak{E}}^{\mathit{tr}}:{\mathfrak{E}}^{\mathit{tr}}}}\right)}^2{\mathfrak{E}}^{\mathit{tr}}:{\mathfrak{E}}^{\mathit{tr}}={\left(\sqrt{\frac{3}{2}}\frac{G}{4}({\kappa }^{\mathit{ac}}+\varpi \mathrm{\Delta }\lambda )+\overline{\tau }\right)}^2$$in which the evolution Eqs. (37), (39), (43) and (44) are used.

As hardenings are considered constant at iterations, Eq. (49), the plastic multiplier is the smallest value of:

50

$$\mathrm{\Delta }\lambda =\left[\left(\sqrt{{\mathfrak{E}}^{\mathit{tr}}:{\mathfrak{E}}^{\mathit{tr}}}-({\lambda }^{ac\zeta }+q^{\mathit{ac}})-{\kappa }^{\mathit{ac}}\right)-\overline{\tau }/\left(\sqrt{\frac{3}{2}}\frac{G}{4}\right)\right]/\left(\zeta (1+\beta )+\eta \right)$$and

51

$$\mathrm{\Delta }\lambda =\left[\left(\sqrt{{\mathfrak{E}}^{\mathit{tr}}:{\mathfrak{E}}^{\mathit{tr}}}-({\lambda }^{ac\zeta }+q^{\mathit{ac}})+{\kappa }^{\mathit{ac}}\right)+\overline{\tau }/\left(\sqrt{\frac{3}{2}}\frac{G}{4}\right)\right]/\left(\zeta (1+\beta )-\eta \right)$$

Finally, the elastoplastic part of the viscoelastoplastic algorithm constitutive tensor is:

52

$$\mathfrak{C}={\mathfrak{C}}^{\mathit{vol}}+\frac{G}{4}(1-\mathrm{\Delta }\lambda )\frac{{\partial }^2({\overline{I}}_1+{\overline{I}}_2)}{\partial \mathit{E}\otimes \partial \mathit{E}}={\mathfrak{C}}^{\mathit{vol}}+{\mathfrak{C}}^{\mathit{ep}}$$in which ${\mathfrak{C}}^{\mathit{vol}}$ is the elastic volumetric part. Thus, see Eq. (45), the stresses components at the elastoplastic part of the proposed viscoelastoplastic model are given by

53

$${\mathbf{S}}_{\mathit{EP}}^i=\frac{G}{4}\left[{\mathfrak{E}}^i-{\lambda }_i^{\zeta }\frac{{\mathfrak{E}}^i}{\sqrt{{\mathfrak{E}}^i:{\mathfrak{E}}^i}}\right]$$in which one recovers components $i$ in the isochoric elastoplastic ($\mathit{EP}$) stress representation.

Viscoelastoplastic stress expressions

Substituting into Eq. (23) the elastic stresses of the viscoelastic models by their elastoplastic counterparts given by Eq. (53), Eq. (23) turns into

54

$${\left({\mathbf{S}}_T^i\right)}_{s+1}=\frac{G}{4}{\left[{\mathfrak{E}}^i-{\lambda }_i^{\zeta }\frac{{\mathfrak{E}}^i}{\sqrt{{\mathfrak{E}}^i:{\mathfrak{E}}^i}}\right]}_{s+1}+\frac{\overline{G}}{4}\left({\left({\dot{\mathfrak{E}}}^i\right)}_{s+1}+\frac{{\left({\mathfrak{E}}_e^i\right)}_s}{\mathrm{\Delta }t}\right)/\left(\frac{\overline{G}}{\eta }+\frac{1}{\mathrm{\Delta }t}\right)$$for the model based on isochoric strain rate components, Fig. 2a.

Making the same for the model based on the time rate of isochoric invariants, Fig. 2b, Eq. (29) turns into

55

$${\mathbf{S}}_T^i=\frac{1}{4}\left[G\left[1-\frac{{\lambda }_i^{\zeta }}{\sqrt{{\mathfrak{E}}^i:{\mathfrak{E}}^i}}\right]+\frac{\overline{\eta }{\dot{\overline{I}}}_i}{\left(1+\left(\overline{\eta }{\dot{\overline{I}}}_i/\overline{G}\right)\right)}\right]{\mathfrak{E}}^i$$

The contribution of these isochoric viscoelastoplastic stresses in the global equilibrium of solids is discussed in the next section, mainly in Eqs. (56) and (62).

Finite element, solution technique and specific models

Preliminaries

One adopts a total Lagrangian finite element based in positions [50, 51] that has as starting point the following weak form of equilibrium equation (achieved from the virtual work principle):

56

$$\delta W=\int _{V_0}^{}{\rho }_0\overrightarrow{\ddot{y}}·\delta \overrightarrow{y}d{V}_0-\int _{V_0}^{}\overrightarrow{b}_{}^0·\delta \overrightarrow{y}d{V}_0-\int _{S_0}^{}{\overrightarrow{p}}_0·\delta \overrightarrow{y}d{S}_0+\int _{V_0}^{}{\mathbf{S}}^{\mathit{cpl}}:\delta \mathbf{E}d{V}_0=0$$in which $\delta \overrightarrow{y}$ is a position variation, ${\rho }_0$ is the Lagrangian density, ${\overrightarrow{b}}_0$ is the volumetric force, ${\overrightarrow{p}}_0$ is the surface force, ${\mathbf{S}}^{\mathit{cpl}}={\mathbf{S}}^{\mathit{vol}}+{\mathbf{S}}_T^i$ is the complete current viscoelastoplastic second Piola–Kirchhoff stress considering both volumetric and isochoric parts—see Eqs. (16), (54) and (55)—and $\delta \mathbf{E}$ is a variation of the Green–Lagrange strain $\mathbf{E}$.

In Fig. 3, one can see, as an illustration, the deformation function approximated by a prismatic finite element that has cubic approximation in the triangular base and linear approximation in its thickness.

[See PDF for image]

Fig. 3

Deformation function by positional FEM

Observing Fig. 3, the initial and current mappings are written as:

57

$$\overrightarrow{f}_{}^0=\overrightarrow{x}={\varphi }_{\alpha }\left(\overrightarrow{\xi }\right){\overrightarrow{X}}_{\alpha }\text{and}\overrightarrow{f}_{}^1=\overrightarrow{y}={\varphi }_{\alpha }\left(\overrightarrow{\xi }\right)Y_{\alpha }$$being $\overrightarrow{x}$ and $\overrightarrow{y}$ initial and current domain positions, ${\varphi }_a$ shape functions related to node $\alpha$ and ${\overrightarrow{X}}_{\alpha }$ and ${\overrightarrow{Y}}_{\alpha }$ nodal positions. One should note that there is a summation in $\alpha$ in Eq. (57). Observing Fig. 3, the deformation results as the composition:

58

$$\overrightarrow{y}(\overrightarrow{x})=\overrightarrow{f}(\overrightarrow{x})={\overrightarrow{f}}^1(\overrightarrow{\xi })\circ {\left({\overrightarrow{f}}^0(\overrightarrow{\xi })\right)}^{-1}$$and its gradient is written as:

59

$$\mathbf{A}={\mathbf{A}}^1{\left({\mathbf{A}}^0\right)}^{-1}$$numerically representing the chain rule. The mappings’ gradients are given by:

60

$${\mathbf{A}}^0=\frac{\partial {\overrightarrow{f}}^0}{\partial \overrightarrow{\xi }}\text{and}{\mathbf{A}}^1=\frac{\partial {\overrightarrow{f}}^1}{\partial \overrightarrow{\xi }}$$

Knowing that the deformation gradient $\mathbf{A}$ is a function of nodal positions the variation of the Green–Lagrange strain is written as:

61

$$\delta \mathbf{E}=\frac{\partial \mathbf{E}}{\partial \overrightarrow{Y}}·\delta \overrightarrow{Y}$$

Thus, the approximate form of Eq. (56) is:

62

$$\left(\int _{V_0}^{}{\rho }_0\overrightarrow{\varphi }\otimes \overrightarrow{\varphi }d{V}_0·\overrightarrow{\ddot{Y}}-\int _{V_0}^{}\overrightarrow{\varphi }\otimes \overrightarrow{\varphi }d{V}_0·{\overrightarrow{B}}^0-\int _{S_0}^{}\overrightarrow{\phi }\otimes \overrightarrow{\phi }d{S}_0·\overrightarrow{P}+\int _{V_0}^{}\left({\mathbf{S}}_T^1+{\mathbf{S}}_T^2+{\mathbf{S}}^{\mathit{vol}}\right):\frac{\partial \mathbf{E}}{\partial \overrightarrow{Y}}d{V}_0\right)·\delta \overrightarrow{Y}=0$$in which domain and surface forces are also approximated by shape functions $\overrightarrow{\varphi }$ and $\overrightarrow{\phi }$. Remembering that one can eliminate $\delta \overrightarrow{Y}$ in Eq. (62) as it is arbitrary one writes:

63

$$\mathbf{M}·\overrightarrow{\ddot{Y}}-{\overrightarrow{F}}^{\mathit{ext}}+{\overrightarrow{F}}^{\mathit{int}}(\overrightarrow{Y})=\overrightarrow{0}$$in which ${\overrightarrow{F}}^{\mathit{ext}}$ are conservative (in this study) external forces, $\mathbf{M}$ is the mass matrix (constant), and ${\overrightarrow{F}}^{\text{int}}$ are internal viscoelastoplastic forces that depends on positions, strain rates and plastic evolution. More steps to achieve the FEM formulation can be consulted in [50, 51].

Solution procedure

The dynamic equilibrium system of Eq. (63) is solved along time using the Newmark $\beta$ method, which is adequate for constant mass matrix [51]. For a current instant $t_{s+1}$, the approximated movement equations become:

64

$$\overrightarrow{g}\left({\overrightarrow{Y}}_{s+1}\right)={\overrightarrow{F}}_{s+1}^{el.vol}+{\overrightarrow{F}}_{s+1}^{vep.iso}+\frac{\mathbf{M}}{\beta \mathrm{\Delta }{t}^2}{\overrightarrow{Y}}_{s+1}-\mathbf{M}·{\overrightarrow{Q}}_s+\frac{\gamma \mathbf{C}}{\beta \mathrm{\Delta }t}·{\overrightarrow{Y}}_{s+1}+\mathbf{C}·{\overrightarrow{R}}_s-\gamma \mathrm{\Delta }t\mathbf{C}·{\overrightarrow{Q}}_s-{\overrightarrow{F}}_{s+1}^{\mathit{ext}}(t)=\overrightarrow{0}$$where $\gamma$ and $\beta$ are well-known Newmark’s parameters. For completeness, a damping matrix $\mathbf{C}$ is introduced. The internal force is divided into a volumetric elastic ${\overrightarrow{F}}_{s+1}^{el.vol}$ and isochoric viscoelastoplastic ${\overrightarrow{F}}_{s+1}^{vp.iso}$. The quantities

65

$${\overrightarrow{Q}}_s=\left(\frac{{\overrightarrow{Y}}_s}{\beta \mathrm{\Delta }{t}^2}+\frac{{\dot{\overrightarrow{Y}}}_s}{\beta \mathrm{\Delta }t_{}^{}}+\left(\frac{1}{2\beta }-1\right){\overrightarrow{\ddot{Y}}}_s\right)\text{and}{\overrightarrow{R}}_s=\left[{\dot{\overrightarrow{Y}}}_s+\mathrm{\Delta }t\left(1-\gamma \right){\ddot{\overrightarrow{Y}}}_s\right]$$are part of the Newmark’s procedure.

The current position ${\overrightarrow{Y}}_{s+1}$ is understood as a trial position and in Eq. (64) the residuum $\overrightarrow{g}({\overrightarrow{Y}}_{s+1})$ is null only if this trial is the correct solution. Making a Taylor’s expansion of $\overrightarrow{g}$ (truncated at first order) results the Newton–Raphson’s method:

66a

$$\overrightarrow{0}=\overrightarrow{g}(\overrightarrow{Y}_{s+1}^{})\cong \overrightarrow{g}\left({\overrightarrow{Y}}_{s+1}^0\right)+\mathrm{\nabla }\overrightarrow{g}\left({\overrightarrow{Y}}_{s+1}^0\right)·\mathrm{\Delta }\overrightarrow{Y}$$in which

66b

$$\mathbf{H}=\mathrm{\nabla }\overrightarrow{g}\left({\overrightarrow{Y}}_{s+1}\right)={\left(\frac{\partial {\overrightarrow{F}}^{el.vol}}{\partial \overrightarrow{Y}}\right|}_{s+1}+{\left(\frac{\partial {\overrightarrow{F}}^{vep.iso}}{\partial \overrightarrow{Y}}\right|}_{s+1}+\left(\frac{\mathbf{M}}{\beta \mathrm{\Delta }{t}^2}+\frac{\gamma \mathbf{C}}{\beta \mathrm{\Delta }t}\right)={\mathbf{H}}^{\mathit{el}\mathbf{.}\mathit{v}\mathit{o}\mathit{l}}+{\mathbf{H}}^{vep.iso}+{\mathbf{H}}^{\mathit{din}}$$is the Hessian matrix with its parts:.

67

$${\mathbf{H}}^{\mathit{el}\mathbf{.}\mathit{v}\mathit{o}\mathit{l}}=\frac{\partial {\overrightarrow{F}}^{el.vol}}{\partial \overrightarrow{Y}}=\int _{V_0}^{}\frac{\partial \mathbf{E}}{\partial \overrightarrow{Y}}:\frac{\partial {\mathbf{S}}^{\mathit{vol}}}{\partial \mathbf{E}}:\frac{\partial \mathbf{E}}{\partial \overrightarrow{Y}}+{\mathbf{S}}^{\mathit{vol}}:\frac{{\partial }^2\mathbf{E}}{\partial \overrightarrow{Y}\otimes \partial \overrightarrow{Y}}d{V}_0$$

68

$${\mathbf{H}}^{\mathit{vp}\mathbf{.}\mathit{i}\mathit{s}\mathit{o}}=\frac{\partial {\overrightarrow{F}}^{vep.iso}}{\partial \overrightarrow{Y}}=\int _{V_0}^{}\frac{\partial \mathbf{E}}{\partial \overrightarrow{Y}}:\frac{\partial {\mathbf{S}}_T}{\partial \mathbf{E}}:\frac{\partial \mathbf{E}}{\partial \overrightarrow{Y}}+{\mathbf{S}}_T:\frac{{\partial }^2\mathbf{E}}{\partial \overrightarrow{Y}\otimes \partial \overrightarrow{Y}}d{V}_0$$

At the same time one is solving Eq. (66), plastic evolution may take place in Eq. (48).

The correction $\mathrm{\Delta }\overrightarrow{Y}$ (calculated at Eq. (66)) is used to improve the trial solution by ${\overrightarrow{Y}}_{s+1}={\overrightarrow{Y}}_{s+1}+\mathrm{\Delta }\overrightarrow{Y}$ until $∥\mathrm{\Delta }\overrightarrow{Y}∥/∥\overrightarrow{X}∥\le Tolerance$ [50, 51].

The values $\partial \mathbf{E}/\partial \overrightarrow{Y}$ and ${\partial }^2\mathbf{E}/(\partial \overrightarrow{Y}\otimes \partial \overrightarrow{Y})$ are well known, see, for example, the author’s studies [48, 50]. In next items, the necessary values $\partial {\mathbf{S}}^{el.vol}/\partial \mathbf{E}$ and $\partial ({\mathbf{S}}_T^1+{\mathbf{S}}_T^2)/\partial \mathbf{E}$ are provided.

(a) Elastic volumetric part $\partial {\mathbf{S}}^{\mathit{vol}}/\partial \mathbf{E}$

69

$${\mathfrak{C}}^{\mathit{vol}}=\frac{\partial {\mathbf{S}}^{el.vol}}{\partial \mathbf{E}}=\frac{K}{4n}\left\{\left((2n-1)J^{2(n-1)}+(2n+1)J^{-2(n+1)}\right){\mathfrak{E}}^{\mathit{vol}}\otimes {\mathfrak{E}}^{\mathit{vol}}+\left(J^{2n-1}-J^{-(2n+1)}\right)\frac{\partial {\mathfrak{E}}^{\mathit{vol}}}{\partial \mathbf{E}}\right\}$$with

70

$$\frac{\partial {\mathfrak{E}}_{j\gamma }^{\mathit{vol}}}{\partial E_{\mathit{oz}}}=J\left\{D_{j\gamma }{D}_{\mathit{oz}}-2D_{\mathit{jo}}{D}_{z\gamma }\right\}$$in which $\mathbf{D}={\mathbf{C}}^{-1}$.

b) Viscoelastoplastic isochoric part$\partial ({\mathbf{S}}_T^1+{\mathbf{S}}_T^2)/\partial \mathbf{E}$.

The isochoric components depend on the choice of the viscosity strategy and if the viscous parameter has some dependence on isochoric invariants or on the isochoric strains rates. Remembering that ${\mathbf{S}}_T={\mathbf{S}}_T^1+{\mathbf{S}}_T^2$ and that ${\eta }_i$ and ${\overline{\eta }}_i$ are presented without indices and arguments one has:

(b1) Viscoelastoplastic model 1—strains rate components with viscous parameter depending on strain invariants and viscous rate components.

Using Eqs. (52) and (54), after some algebraic manipulations, see Appendix C, one finds the final expression of the isochoric viscoelastoplastic part of model 1:

71

$${\left(\frac{\partial {\mathbf{S}}_v^i}{\partial \mathbf{E}}\right)}_{s+1}={\mathfrak{C}}^{\mathit{ep}}+\left\{\frac{\overline{G}}{4\mathrm{\Delta }t}\left(\frac{\partial \mathfrak{E}_{}^i}{\partial \mathbf{E}}\right)_{s+1}^{}+\overline{G}\left({\left(\eta \right)}^{-2}\frac{\partial \eta }{\partial \mathbf{E}}\right)\otimes {\left({\mathbf{S}}_v^i\right)}_{s+1}\right\}/\left(\frac{\overline{G}}{\eta }+\frac{1}{\mathrm{\Delta }t}\right)={\mathfrak{C}}^{\mathit{ep}}+{\mathfrak{C}}^{\mathit{vis}}$$

When considering $\eta$ as a function of ${\dot{\mathfrak{E}}}^i$ and ${\overline{I}}_i$ in the following form:

72

$$\eta ={\eta }^a({\dot{\mathfrak{E}}}_v^i){\eta }^b({\overline{I}}_i)$$one finds (also see Appendix C)

73

$$\frac{\partial \eta }{\partial \mathbf{E}}={\eta }^b\left[\overline{G}\frac{\partial {\eta }^a}{\partial {\dot{\mathfrak{E}}}_v^i}:\frac{\partial {\mathfrak{E}}^i}{\partial \mathbf{E}}-\left({\dot{\mathfrak{E}}}_v^i:\partial {\eta }^a/\partial {\dot{\mathfrak{E}}}_v^i\right){\eta }^a\frac{\partial {\eta }^b}{\partial {\overline{I}}_i}{\mathfrak{E}}^i\right]/\left(\overline{G}\mathrm{\Delta }t+\eta +\left({\dot{\mathfrak{E}}}_v^i:\partial {\eta }^a/\partial {\dot{\mathfrak{E}}}_v^i\right){\eta }^b\right)+{\eta }^a\frac{\partial {\eta }^b}{\partial {\overline{I}}_i}{\mathfrak{E}}^i$$

A simpler alternative that does not interfere significantly in the convergence rate is:

74

$$\frac{\partial \eta }{\partial \mathbf{E}}={\eta }^b\overline{G}\left[{\dot{\mathfrak{E}}}_v^i:\partial {\eta }^a/\partial {\dot{\mathfrak{E}}}_v^i\right]/\left[\overline{G}\mathrm{\Delta }t+{\eta }^a+\left({\dot{\mathfrak{E}}}_v^i:\partial {\eta }^a/\partial {\dot{\mathfrak{E}}}_v^i\right)\right]+{\eta }^a\frac{\partial {\eta }^b}{\partial {\overline{I}}_i}{\mathfrak{E}}^i$$(b2) Viscoelastoplastic model 2—strain invariant rates with viscous parameter depending on strain invariants.

In this case, one assumes $\overline{\eta }({\overline{I}}_i)$, and using Eqs. (52) and (55), one is able to write the final expression of the isochoric viscoelastoplastic part of model 2:

75

$${\left(\frac{\partial {\mathbf{S}}_T^i}{\partial E}\right|}_{s+1}={\mathfrak{C}}^{\mathit{ep}}+\left\{\frac{\overline{G}}{4}\left[Q{\mathfrak{E}}_{s+1}^i\otimes {\mathfrak{E}}_{s+1}^i+Z\frac{\partial {\mathfrak{E}}_{s+1}^i}{\partial \mathbf{E}}\right]-Q{\mathfrak{E}}_{s+1}^i\otimes {\left({\mathbf{S}}_v^i\right)}_{s+1}\right\}\left./\phantom{\left\{\frac{\overline{G}}{4}\left[Q{\mathfrak{E}}_{s+1}^i\otimes {\mathfrak{E}}_{s+1}^i+Z\frac{\partial {\mathfrak{E}}_{s+1}^i}{\partial \mathbf{E}}\right]-Q{\mathfrak{E}}_{s+1}^i\otimes {\left({\mathbf{S}}_v^i\right)}_{s+1}\right\}\left(1+Z\right)}\right)\left(1+Z\right)={\mathfrak{C}}^{\mathit{ep}}+{\mathfrak{C}}^{\mathit{vis}}$$with auxiliary values given as:

76

$$Z=\left(\overline{\eta }{\dot{\overline{I}}}_i/\overline{G}\right)\frac{\partial Z}{\partial \mathbf{E}}=\frac{1}{\overline{G}}\left(\frac{d\overline{\eta }}{d{\overline{I}}_i}{\left({\dot{\overline{I}}}_i\right)}_{s+1}+\frac{\overline{\eta }}{\mathrm{\Delta }t}\right){\mathfrak{E}}_{s+1}^i=Q{\mathfrak{E}}_{s+1}^i$$

As model 1 is more general, see Eqs. (24) through (32), it is not introduced more general viscous parameter dependence for model 2.

Specific adopted viscosities

In order to model the behavior of the tested materials, viscous parameters $\eta$ and $\overline{\eta }$ should be more elaborated than simply constant. Thus, the following sets are used in the general expressions presented above:

Set 1:

77

$${\eta }^a=\eta _{}^{\ast }+\upsilon {\left({\dot{\mathfrak{E}}}_v^i:{\dot{\mathfrak{E}}}_v^i+rv\right)}^{\ell }\text{and}{\eta }^b=1$$

Set 2:

78

$${\eta }^a=\tilde{\eta }{e}^{\tilde{\upsilon }{\dot{\mathfrak{E}}}_v^i:{\dot{\mathfrak{E}}}_v^i}\text{and}{\eta }^b=1$$

Set 3:

79

$${\eta }^a=\eta {\left({\dot{\mathfrak{E}}}_v^i:{\dot{\mathfrak{E}}}_v^i+rv\right)}^{\ell }\text{and}{\eta }^b=\epsilon +{({\overline{I}}_1-\rho )}^{\overline{\upsilon }-1}$$

${\eta }^a$ and ${\eta }^b$ are used in Eq. (72) for model 1 and ${\eta }^a=\overline{\eta }$ is used for model 2.

Finally, the derivatives of the isochoric strains regarding the Green–Lagrange strain are given by:

80

$$\frac{{\partial }^2{\overline{I}}_1}{\partial E_{\mathit{ij}}\partial E_{k\ell }}=\frac{\partial {\mathfrak{E}}_{\mathit{ij}}^{iso1}}{\partial E_{k\ell }}=\frac{4}{3}{J}^{-2/3}\left\{\frac{1}{3}\left\{D_{\mathit{ij}}{D}_{k\ell }+3D_{\mathit{ik}}{D}_{\ell j}\right\}{I}_1-D_{\mathit{ij}}{\delta }_{k\ell }-D_{k\ell }{\delta }_{\mathit{ij}}\right\}$$

81

$$\frac{{\partial }^2{\overline{I}}_2}{\partial E_{\mathit{ij}}\partial E_{k\ell }}=\frac{\partial {\mathfrak{E}}_{\mathit{ij}}^{iso2}}{\partial E_{k\ell }}=\frac{8}{3}{J}^{-4/3}\left\{\begin{array}{c}\left[\frac{2}{3}{D}_{\mathit{ij}}{D}_{k\ell }+D_{\mathit{ik}}{D}_{\ell j}\right]I_2-\left[C_{\mathit{zz}}\left(D_{\mathit{ij}}{\delta }_{k\ell }+D_{k\ell }{\delta }_{\mathit{ij}}\right)\right]+\\ +D_{\mathit{ij}}{C}_{k\ell }+D_{k\ell }{C}_{\mathit{ij}}+\frac{3}{2}\left[{\delta }_{\mathit{ij}}{\delta }_{k\ell }-{\delta }_{\mathit{jk}}{\delta }_{i\ell }\right]\\ \end{array}\right\}$$in which $\mathbf{D}={\mathbf{C}}^{-1}$. For completeness, Eqs. (80) and (81) are given in index notation.

Examples

This section presents examples that validate and clarify some properties of the formulation, such as: (1) The independence of volumetric loads regarding the viscous isochoric behavior. (2) Validate the viscoelastic Zener-like models regarding instantaneous and time dependent responses including rate-dependent viscosity. (3) Verify the proposed elastoplastic model for cases independent of strain rates. (4) Verify the objectivity of the proposed elastoplastic framework—the same as the viscous one. (5) Confirms the proposed formulation with experimental results, showing the applicability of the technique.

Testing orthogonality of the numerical solution

This example is presented in order to verify that the isochoric viscosity (of the proposed framework) has no influence in the volumetric change and vice versa. For this purpose, two linear prismatic finite elements are used to discretize a unitary cube $(1\text{mm}\times 1\text{mm})$ with properties $G=10\text{MPa}$, $\overline{G}=8\text{MPa}$ and $K=1000\text{MPa}$. The boundary conditions and the adopted discretization can be seen in Fig. 4. An instantaneous hydrostatic stress is applied by modifying Eq. (13), i.e.,$${\mathbf{S}}^{\mathit{vol}}=\left\{\frac{K}{4n}\left(J^{2n-1}-J^{-(2n+1)}\right)-1000\text{MPa}\right\}{\mathfrak{E}}^{\mathit{vol}}$$

[See PDF for image]

Fig. 4

Geometry and discretization

No yielding takes place and an instantaneous displacement occurs for all proposed models, resulting in a final and constant (static analysis) volume $V=4015{\text{mm}}^3$.

Thus, when an isolated volumetric change is prescribed—using the proposed framework—no viscous stress takes place.

Viscoelastic rate-dependent response of damping rubbers—experimental validation

A recent study performed by [56] investigates the properties of high damping rubber (HDR), highlighting the significant influence of strain rate on viscous properties. The researchers employed a Kröner–Lee-based approach considering the dependence of viscosity on the rate of isochoric strain.

In this example, the elastic part of the proposed model 3.1 together with the specific viscous model called Set 1 (${\eta }^b=1$) is adopted. In order to do so, one uses $\overline{\tau }=100\text{MPa}$, i.e., avoiding plasticity. To validate the proposed model capability, an experimental shear relaxation test extracted from reference [56] is replicated for high damping rubber (HDR) using the same discretization scheme as the previous example. The following parameters are used $K=1000\text{MPa}$, $G=0.5\text{MPa}$, $\overline{G}=0.6\text{MPa}$, $\eta _{}^{\ast }=10.0\text{MPa}\text{s}$, $\ell =-0.5$, $\upsilon =1.8\text{MPa}{\text{s}}^{2\ell +1}=1.8\text{MPa}$ and $rv=0.0004{\text{s}}^{-2}$.

As the constitutive model used by [56] is not similar to the ones presented here, its physical properties are not listed.

Figure 5 compares the Cauchy stress predicted by the presented numerical model for three imposed shear strains with the experimental data from [56]. The simulations used a time step of $\mathrm{\Delta }t=1.0s$.

[See PDF for image]

Fig. 5

Comparing numerical results for shear relaxation (values of $\gamma$ in caption) and experimental data from [56], histories separated by $200s$

The intermediate shear strain level is used to calibrate the adopted constants and results are practically equal to the experimental data. Despite its relative simplicity, the proposed model achieves good agreement with the ones presented by the numerical results of reference [56]. In order to analyze the sensitivity of the polynomial model Set 1 regarding parameter $\upsilon$, one keeps $\eta _{}^{\ast }=10.0\text{MPa}\text{s}$ constant in Fig. 6a, while to analyze the sensibility of $\eta _{}^{\ast }$, one keeps $\upsilon =0.0$ in Fig. 6b. For both cases, it is adopted: $rv=0.0003{\text{s}}^{-2}$, $\ell =-0.5$ and $\gamma =2.50$.

[See PDF for image]

Fig. 6

Sensitivity SET 1 regarding parameters $\upsilon$ and $\eta _{}^{\ast }$

Applying an instantaneous distortion $\gamma =1.0$ Fig. 7 shows the sensitivity of exponential model SET 2 regarding $\tilde{\upsilon }$ keeping $\tilde{\eta }=100\text{MPa}\text{s}$ constant and the behavior of the polynomial model SET 1 regarding $\ell$ keeping $\eta _{}^{\ast }=10.0\text{MPa}\text{s}$ and $\upsilon =10.8\text{MPa}{\text{s}}^{2\ell +1}$ constants.

[See PDF for image]

Fig. 7

Sensitivity of exponential and polynomial models

The influence of parameter $\ell$ in the viscosity behavior of SET 1 is summarized as: When adopting $\ell <0$ viscosity becomes relatively low for high strain rates and relatively high for low strain rates. For $\ell >0$, the viscosity will be relatively high for high strain rates and relatively low for low strain rates. The exponential model is less versatile than the polynomial one; thus in next viscoelastoplastic analyses, only polynomial models are used.

Rate-independent elastoplastic model

In this example, the same unitary cube is used to model an elastoplastic material. In this case, considering the polynomial rate-dependent viscoelastoplastic model (SET 1), one adopts the following material properties: ${\overline{\sigma }}_y/2=\overline{\tau }=0.4\text{MPa}$, $K=800\text{MPa}$, $G_1=G_2=G=24\text{MPa}$, $\overline{G}=0$, $\eta _{}^{\ast }=0$, $\upsilon =0$, $rv=0$ and $\ell =0$. In reality $\overline{G}$ and $\eta _{}^{\ast }$ are very small values to avoid dividing by zero into calculations. Three different constant hardenings are adopted and their values are given in the captions of Fig. 8. Position control is applied in the third direction for both tensioning and compressing situations; in addition a cyclic loop using kinematic hardening to show the energy consistence of the plastic flow.

[See PDF for image]

Fig. 8

Cyclic behavior of an elastoplastic material for different constant hardenings

Tensioning phases:

• First phase: 200 equally spaced steps from ${\lambda }_3=1.0$ until ${\lambda }_3=1.4$

• Second phase: 100 equally spaced steps from ${\lambda }_3=1.4$ until ${\lambda }_3=1.2$

• Third phase: 200 equally spaced steps from ${\lambda }_3=1.2$ until ${\lambda }_3=1.6$

Compressing phases:

• First phase: 200 equally spaced steps from ${\lambda }_3=1.0$ until ${\lambda }_3=0.6$

• Second phase: 100 equally spaced steps from ${\lambda }_3=0.6$ until ${\lambda }_3=0.8$

• Third phase: 200 equally spaced steps from ${\lambda }_3=0.8$ until ${\lambda }_3=0.4$

For both tensioning and compressing cases, Fig. 8 shows the Cauchy stress ${\sigma }_{33}$ as a function of $\left({\lambda }_3-1\right)$. It is important to mention that the arrows present in Fig. 8 represent the sense of the applied strain control for each phase.

As one can see the “geometric hardening” is very pronounced at the compressive situation; however, both cases perfectly close the loading cycle, indicating the energy consistence of the model. Figure 9 presents the initial and final deformed configurations for both cases with $H^c=0.125\text{MPa}$.

[See PDF for image]

Fig. 9

Final position for compression ${\lambda }_3=0.4$ and tension ${\lambda }_3=1.6$, $H^c=0.125MPa$

To show that the formulation proposed in this study is objective regarding large rotations, the same cube with $H^c=0.125\text{MPa}$ is analyzed using other boundary conditions. In this case, all degrees of freedom are restricted and an isochoric strain (distortion $0<\gamma <0.5$) is applied in two ways.

In the first way, the cube is kept in its original orientation and a movement in direction $x_1$ is applied at points with coordinates $x_2=1$. In the second way, at the same time that the distortion is applied, the cube is rotated in the plane $x_1{x}_2$ from 0 to − 45 degrees. Figure 10a shows three snapshots without rotation the axes, and Fig. 10b shows the same three snapshots with the imposed rigid body rotation. A total of 1000 position increments were used to run the example in both situations. Figure 10 also shows the values of the Cauchy stress component ${\sigma }_{12}$ for the presented snapshots. Obviously, the shear stresses are different for both cases.

[See PDF for image]

Fig. 10

Shear stresses—${\sigma }_{12}(MPa)$—for applied distortion, with and without rigid body rotation

Figure 11a shows a graph with the three nonzero Cauchy stress components of the problem. The letter R is used to indicate the case with rigid body rotation and the letter H indicates the case without rigid body rotation. Finally, the Cauchy stresses—obtained in the analysis that includes rigid body rotation—are rotated in a way that their reference is changed according to the current orientation of the body. The resulting values indicated by the symbol RR are compared in Fig. 11b with the stresses of case H—non-rotated case. As can be seen, the responses presented in Fig. 11b are coincident, revealing the objectivity of the proposed model.

[See PDF for image]

Fig. 11

Cauchy stresses as a function of applied distortion (with and without rigid rotation)

Pinched cylinder with rigid end diaphragms

In this example, the well-known benchmark “Pinched cylinder with rigid end diaphragms” is solved in its elastoplastic version, see [57, 58–59], for instance. Its geometry, loading and boundary conditions are depicted in Fig. 12a.

[See PDF for image]

Fig. 12

Problem definition—loading boundary conditions and discretization

Only one-eighth of the problem is discretized, see Fig. 13b, due to the symmetry of the problem. The mesh has 20 divisions (triangular elements) along the circumferential direction and 12 divisions along the longitudinal axis. It has only one division along thickness (with cubic approximation—triangular) and each division along longitudinal axis is depicted by three subdivisions due to a particularity of the in-house post-processor (prismatic element). This is a dense discretization, called 20 × 12 prismatic elements with cubic approximation, resulting in 27,084 DOF.

[See PDF for image]

Fig. 13

Final position of analysis, a vertical displacement, b cumulated plastic multiplier ${\lambda }_1$, c Cumulated plastic multiplier ${\lambda }_2$

The physical properties, directly extracted from reference [57], are: $E=3000$, $\nu =0.3$, ${\sigma }_y=24.3$ and isotropic hardening $H^i=300$. Thus, one finds the properties to be used in the proposed model, i.e., $K=2500$, $G=1153\text{.85}$ and $\overline{\tau }=12\text{.15}$. In this case, one can adopt directly $H_1^i=H_2^i=37.5$ with null kinematic hardening. No viscous parameters are used. In Fig. 13, the final position of the analysis (only 1/8 of the body) is presented. Over the deformed mesh of Fig. 13 one can see the vertical displacement (a), the cumulated plastic multiplier in direction ${\overline{I}}_1$ (b) and in direction ${\overline{I}}_2$ (c).

In Fig. 14, the force versus displacement behavior achieved by the proposed model is compared with the ones achieved by [57] using a mesh of 16 × 16 SHB20 implemented in ABACUS and the solution given by [59].

[See PDF for image]

Fig. 14

Force versus displacement

As one can see, the results are in good agreement with those presented by references. It is interesting to mention that references’ results also present depressions for displacement control, not perceptible here as only selected points are taken. It is interesting noting that the two elements that contain the loaded nodes are considered elastic in the proposed model. Equal steps of $0.25$ was adopted when using displacement control, limiting the displacement step to less than $10\%$ of the shell thickness. For force control, $100$ steps of $23.667$ were sufficient to achieve convergence for $tol=6.5x{10}^{-4}$ in positions. It is interesting to mention that when force control is used, snaps through do not appear in Fig. 14. From the above results, it is possible to conclude that the proposed model can be used as an alternative approach for large-strain plasticity associated with large displacements and rotations.

General viscoelastoplastic test—high strain rate with strain rate and isochoric invariant dependence

This example shows the possibility of the proposed model in reproducing the experimental response under fast loading of a complex material, for example, the polytetrafluoroethylene. The response regarding long-term loads is shown in example 6.6, indicating ways for future implementations to be carried out by researchers who will use the Flory's multiplicative decomposition framework to model specific materials (which is not the objective of the present study).

The present example makes extreme demands on the developed viscoelastoplastic formulation, as it refers to the experimental test of polytetrafluoroethylene (room temperature) at very different strain rates. Four longitudinal engineering strain rates are considered: $\dot{\epsilon }={10}^{-4}{s}^{-1}$, $\dot{\epsilon }={10}^{-2}{s}^{-1}$, $\dot{\epsilon }=1.0s^{-1}$ e $\dot{\epsilon }={10}^3{s}^{-1}$.

It should be noted that the highest applied rate is ${10}^7$ times larger than the lowest strain rate, revealing how extreme is the performed test. The experimental test was carried out by [26] using a servo-hydraulic testing machine (MTS) controlling the strain rate. All necessary actions to reduce friction and discount strain of the MTS were followed by [26]. Being an experimental work, reference [26] presented a uniaxial viscoelastoplastic model in additive decomposition, with the aim of modeling only the simple uniaxial behavior of the material.

Reference [2] presented a model based on the Kröner–Lee decomposition, but did not present results for the more aggressive strain rate ($\dot{\epsilon }={10}^3{\text{s}}^{-1}$). As models are different, only the constants adopted in this study are presented as follows:

• Elastic constants: $K=1820\text{MPa}$, $G=G_1=G_2=98\text{MPa}$, $\overline{G}={\overline{G}}_1={\overline{G}}_2=98\text{MPa}$

• Elastic limit:${\overline{\tau }}_1={\overline{\tau }}_2=\tau =3.41\text{MPa}$

• Plastic constants: Kinematic and isotropic hardenings are constant at a time interval and equal for both isochoric directions, the adopted values are given in Table 1.

• Viscous constants for SET 3: $\eta =32\text{MPa}{\text{s}}^{4.2}$, $\ell =-0.45$, $rv=0.01{\text{s}}^{-2}$, $\epsilon =0$, $\rho =2.9$ and $\overline{\upsilon }=1.6$

Table 1. Kinematic and Isotropic hardening for plastic multipliers’ intervals

Plastic multipliers’ intervals $\lambda ={\lambda }_1={\lambda }_2$	$H^C=H_1^C=H_2^C$ (MPa)	$H^I=H_1^I=H_2^I$ (MPa)
$0\le \lambda <0.15$	$7.414$	$0.726$
$0.15\le \lambda <0.25$	$6.403$	$0.627$
$0.25\le \lambda <0.50$	$5.729$	$0.561$
$0.50\le \lambda <0.70$	$3.370$	$0.330$
$0.70\le \lambda <0.90$	$4.044$	$0.396$
$0.90\le \lambda$	$2.696$	$0.426$

The original experimental specimen [26] is cylindrical with a diameter of $25.4\text{mm}$ and length of $30.5\text{mm}$ and was compressed following the above-mentioned strain rates in its longitudinal direction. After carrying out static and dynamic analysis with density $2200{\text{kg/m}}^3$, it is verified that there is no significant influence of inertia for this dimension; thus, a cubic numerical test specimen with side of $20\text{mm}$ is chosen and the static analysis is performed. Different time steps were used for each strain rate as follows: $\mathrm{\Delta }t=1.0\times {10}^0\text{s}$ for $\dot{\epsilon }={10}^{-4}{\text{s}}^{-1}$, $\mathrm{\Delta }t=1.0\times {10}^{-2}\text{s}$ for $\dot{\epsilon }={10}^{-2}{\text{s}}^{-1}$, $\mathrm{\Delta }t=1.0\times {10}^{-4}\text{s}$ for $\dot{\epsilon }=1.0{\text{s}}^{-1}$ and $\mathrm{\Delta }t=1.0\times {10}^{-7}\text{s}$ for $\dot{\epsilon }={10}^3{\text{s}}^{-1}$.

A linear spatial approximation was used for the numerical solutions—model 3.1 (SET 3)—with a discretization similar to that presented in the first example—Fig. 4.

In Fig. 15, the results obtained by the present formulation are compared with the experimental results of [26] in a situation of monotonically increasing strain. A dynamic result is included for an additional numerical cubic specimen with side $20\text{mm}$ under $\dot{\epsilon }={10}^3{\text{s}}^{-1}$ showing the influence of inertia (that appears fir this strain rate) at the beginning of the test that is propagated to the global behavior of the solid throughout the analysis. It is important to note that, for the dynamic analysis, no restrictions are imposed in the transverse directions of the specimen (no sliding planes).

[See PDF for image]

Fig. 15

Different rates—increasing strain—comparison with experimental results from [26]

Figure 16 shows the result when imposing strain with constant strain rate $\dot{\epsilon }={10}^{-2}{\text{s}}^{-1}$ punctuated with periods of relaxation lasting 36 s. The proposed formulation result is compared with the experimental result given by [26]. The relaxation constant stretches values are: at the loading phase ${\lambda }_3=-0.0832$ and ${\lambda }_3=-0.1832$ and at the unloading phase ${\lambda }_3=-0.1832$. The maximum stretch is of ${\lambda }_3=0.2832$ and the minimum stretch (where the strain-controlled stress cancels out) is of ${\lambda }_3=-0.1712$. From this last stretch value, the numerical test specimen was allowed to conform freely, a procedure similar to for the experimental test by [26]. It is observed that the stress reductions in the 36 s relaxation periods are smaller than the experimental results because the long-term behavior of the material. This long-term behavior is presented in example 6.6.

[See PDF for image]

Fig. 16

Cyclic loading with short-term relaxation—comparison with experimental results from [26]

Considering the complexity of the analysis, the results of the proposed model fit very well with the behavior of the tested material; in addition, it is confirmed that the dimensions of the experimental specimen chosen by reference [26] are suitable even for the most extreme strain rate.

Long-term response of polytetrafluoroethylene:

It is not the objective of this study to include more complex models in the originally developed viscoelastoplastic model. However, Fig. 17 shows how the viscoelastoplastic model should be like for the solution combining short- and long-term behaviors. In Fig. 17, index “$a$” indicates the short term and index “$b$” the long-term branch of the model. Due to the fact that the proposed formulation works in the isochoric space through the Flory’s decomposition, it is possible to solve the general expressions of the model (not carried out in this study) by the equality of the total stress and the allowed additive decomposition in the isochoric space.

[See PDF for image]

Fig. 17

Unified model for short- and long-term behaviors of this specific material

In this study, to illustrate the behavior of this specific material, a series association of three prismatic finite elements is made, two sets of two elements with length of 19.9 mm and two central elements with length of 0.2 mm as shown in Fig. 18a, resulting in a prismatic specimen with base 2 mm × 2 mm and length 4 mm. In Fig. 18, elements type “a” and type “b” are indicated, alluding to the physical properties of the model presented in Fig. 17. Figure 18a also shows the tiny central element with low transverse elastic and viscous moduli (called type “c”) responsible for allowing “free” transversal relative movement between elements “a” and “b.” Thanks to the Hartmann–Neff [53] volumetric strain energy density, Eq. (10), elements “c” do not invert. In Fig. 18b, just as an illustration, the maximum deformation of the specimen is shown.

[See PDF for image]

Fig. 18

Discretization, material types and deformed shape

For material “a,” it is adopted the same parameters used in example 6.5, that is:

• Elastic constants: $K_a=1820\text{MPa}$, $G_a=G_1^a=G_2^a=98\text{MPa}$,$\overline{G}_{}^a={\overline{G}}_1^a={\overline{G}}_2^a=98\text{MPa}$

• Elastic limits:${\overline{\tau }}_1^a={\overline{\tau }}_2^a=\tau _{}^a=3.41\text{MPa}$

• Plastic constants: Kinematic and isotropic hardenings are constant at a time interval and equal for both isochoric directions, the adopted values are given in Table 1, considering material “a.”

• Viscous constants (SET 3): $\eta =32\text{MPa}{\text{s}}^{4.2}$, ${\ell }_a=-0.45$, $r{v}_a=0.01{\text{s}}^{-2}$, ${\epsilon }_a=0$, ${\rho }_a=2.9$ and ${\overline{\upsilon }}_a=1.6$

For material “b,” it is adopted:

• Elastic constants: $K_b=1820\text{MPa}$, $G_b=G_1^b=G_2^b=250\text{MPa}$,$\overline{G}_{}^b={\overline{G}}_1^b={\overline{G}}_2^b=3800\text{MPa}$

• Plastic constants: ${\overline{\tau }}_1^b={\overline{\tau }}_2^b=\tau _{}^b=9.40\text{MPa}$ with constant hardenings $(H_1^C=H_2^C=H^C)$ and $(H_1^I=H_2^I=H_{}^I)$: $H_b^C=404.4\text{MPa}$ and $H_b^I=39.6\text{MPa}$$\forall \lambda$.

• Viscous constants (SET 3): $\eta =120\text{GPa}{\text{s}}^{4.2}$, ${\ell }_b=0.0$, $r{v}_b=0.01{\text{s}}^{-2}$, ${\epsilon }_b=0$, ${\rho }_b=2.9$ and ${\overline{\upsilon }}_b=0.3$

The connecting element “$c$” has viscoelastic behavior with properties:

• Elastic constants: $K_c=1820\text{MPa}$, $G_c=G_1^c=G_2^c=9.9\text{kPa}$,$\overline{G}_{}^c={\overline{G}}_1^c={\overline{G}}_2^c=3.8\text{kPa}$

• Viscous constants (SET 3): ${\eta }_c=0.32\text{MPa}\text{s}$, ${\ell }_c=0.0$, $r{v}_c=0.01{\text{s}}^{-2}$, ${\epsilon }_c=0$, ${\rho }_c=2.9$ and ${\overline{\upsilon }}_c=0.3$

Two situations were tested and compared with the experimental results from reference [26]. The first is the cyclic loading considered in the previous example (Fig. 16) in which relaxation time intervals of 36 s are used, because the model applied there did not include the long-term part “b.” Now the loading and unloading paths continue with the same strain rate $\dot{\epsilon }=0.01{\text{s}}^{-1}$, while relaxation takes 6 h. The results are compared with the experimental ones in Fig. 19. As can be seen, the reduction in stress levels in the relaxation periods is much more pronounced, reproducing the experimental behavior throughout the loading process with greater precision.

[See PDF for image]

Fig. 19

Cyclic result with long-term relaxation—comparison with experimental result from [26]

In the second situation (creep), the applied stresses are controlled and the results are compared with the experimental response. For all load increasing path, the same engineering stress rate $0.160833\text{MPa/s}$ is employed with a time step of $\mathrm{\Delta }t=0.1\text{s}$ and durations $60\text{s}$, $32.1\text{s}$ and $16.5\text{s}$, respectively. After each load increase, constant load levels of 6 h are imposed with values $9.65\text{MPa}$, $14.89\text{MPa}$ and $17.42\text{MPa}$. For the constant loads, a time step of $\mathrm{\Delta }t=21.6\text{s}$ is assumed. Figure 20 presents the engineering strain as a function of time for the entire analysis, showing that the numerical results adequately approximate the material experimental behavior.

[See PDF for image]

Fig. 20

Long-term creep results—comparison with experimental results from [26]

It is worth noting that, differently to the previous examples, the strain indicated in Figs. 19 and 20 corresponds to the total displacement divided by the length of element type “a” and that the engineering stress in Fig. 19 was calculated by dividing the total reactive force by the initial cross-sectional area of the specimen (0.02 m × 0.02 m). It is observed that results of Figs. 19 and 20 are close to the experimental data, even using the numerical device of Fig. 18. It is also commented that no further calibration was sought for this example since it is not a complete constitutive model, just an indicative of how one constructs it using the proposed Flory’s framework.

Conclusions

In this study, a viscoelastoplastic framework for finite strains based on the Flory’s decomposition without using Kröner–Lee multiplicative decomposition was presented. From Flory’s decomposition, the isochoric and volumetric directions of strains are defined as a tensor basis in which directions of plastic flow and viscous dissipation are originally written. Using this premise, a plastic evolution was established with constant isotropic and kinematic hardenings throughout numerical iterations, but which may vary along the load or time steps. Using the stated tensor base, viscous dissipation is introduced in two ways. The first directly applies the isochoric strain rate and the second uses invariant strain rate multiplied by isochoric strain directions. As the viscous dissipation based on strain rate proved to be more general, it is extended to represent rate-dependent viscous relations. The hyperelastic, plastic and viscous models were put together and implemented in an in-house total Lagrangian finite element computational code, assembling the viscoelastoplastic framework. Numerical tests and comparisons with experimental results prove the capabilities of the proposed framework. The use of isochoric directions (or tensor base) makes possible the additive decomposition at large strains, opening the possibility of developing other finite strain viscoelastoplastic models for specific materials. Further improvements include nonlinear hardening and updated FEM Lagrangian implementation.

Acknowledgements

This research has been supported by the São Paulo Research Foundation, Brazil—Grant #2020/05393-4.

Humberto Breves Coda

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

References

1. Nguyen, K; Amores, VJ; Sanz, MA; Montáns, FJ. Thermodynamically consistent nonlinear viscoplastic formulation with exact solution for the linear case and well-conditioned recovery of the inviscid one. Comput Mech; 2021; 67, pp. 1349-1373.4248602 [DOI: https://dx.doi.org/10.1007/s00466-021-01991-w] 1476.74013

2. Carvalho, PRP; Coda, HB; Sanches, RAK. A large strain thermodynamically-based viscoelastic–viscoplastic model with application to finite element analysis of polytetrafluoroethylene (PTFE). Eur J Mech A Solids; 2023; 97, 4509435 [DOI: https://dx.doi.org/10.1016/j.euromechsol.2022.104850] 1503.74017 104850.

3. Perzyna, P. Fundamental problems in viscoplasticity. Adv Appl Mech; 1966; 9, C pp. 243-377. [DOI: https://dx.doi.org/10.1016/S0065-2156(08)70009-7] 1266.74009

4. Duvaut, G; Lions, J. Les Inequations en Mecanique et en Physique; 1972; Paris, Dunod: 0298.73001

5. Caggiano, A; Martinelli, E; Said Schicchi, D; Etse, G. A modified Duvaut-Lions zero-thickness interface model for simulating the rate-dependent bond behavior of FRP-concrete joints. Compos Part B Eng; 2018; 149, April pp. 260-267. [DOI: https://dx.doi.org/10.1016/j.compositesb.2018.05.010] 1440.74374

6. Chaboche, JL. Constitutive equations for cyclic plasticity and cyclic viscoplasticity. Int J Plast; 1989; 5, May pp. 247-302. [DOI: https://dx.doi.org/10.1016/0749-6419(89)90015-6] 0695.73001

7. Cormeau, I. Numerical stability in quasi-static elasto/viscoplasticity. Int J Numer Methods Eng; 1975; 9, 1 pp. 109-127.436739 [DOI: https://dx.doi.org/10.1002/nme.1620090110] 0293.73022

8. Ibrahimbegović, A; Chorfi, L. Viscoplasticity model at finite deformations with combined isotropic and kinematic hardening. Comput Struct; 2000; 77, 5 pp. 509-525. [DOI: https://dx.doi.org/10.1016/S0045-7949(99)00232-1] 1151.74320

9. Kowalczyk-Gajewska, K; Pieczyska, EA; Golasinski, K; Maj, M; Kuramoto, S; Furutab, T. A finite strain elastic-viscoplastic model of Gum Metal. Int J Plast; 2019; 119, 2018 pp. 85-101. [DOI: https://dx.doi.org/10.1016/j.ijplas.2019.02.017]

10. Nedjar, B. Frameworks for finite strain viscoelastic-plasticity based on multiplicative decompositions. Part I: Continuum formulations. Comput Methods Appl Mech Eng; 2002; 191, 15–16 pp. 1541-1562. [DOI: https://dx.doi.org/10.1016/S0045-7825(01)00337-1] 1141.74317

11. Wang, WM; Sluys, LJ; De Borst, R. Viscoplasticity for instabilities due to strain softening and strain-rate softening. Int J Numer Methods Eng; 1997; 40, 20 pp. 3839-3864. [DOI: https://dx.doi.org/10.1002/(SICI)1097-0207(19971030)40:20<3839::AID-NME245>3.0.CO;2-6] 0974.74511

12. Carosio, A; Willam, K; Etse, G. On the consistency of viscoplastic formulations. Int J Solids Struct; 2000; 37, 48 pp. 7349-7369. [DOI: https://dx.doi.org/10.1016/S0020-7683(00)00202-X] 1056.74512

13. Heeres, OM; Suiker, ASJ; De Borst, R. A comparison between the Perzyna viscoplastic model and the consistency viscoplastic model. Euro J Mech A/Solids; 2002; 21, 1 pp. 1-12. [DOI: https://dx.doi.org/10.1016/S0997-7538(01)01188-3] 1005.74013

14. Ristinmaa, M; Ottosen, NS. Consequences of dynamic yield surface in viscoplasticity. Int J Solids Struct; 2000; 37, 33 pp. 4601-4622. [DOI: https://dx.doi.org/10.1016/S0020-7683(99)00158-4] 0982.74014

15. Haupt, P. On the concept of an intermediate configuration and its application to a representation of viscoelastic-plastic material behavior. Int J Plast; 1985; 1, 4 pp. 303-316.792657 [DOI: https://dx.doi.org/10.1016/0749-6419(85)90017-8] 0632.73042

16. Simo, J. Algorithms for static and dynamic multiplicative plasticity that preserve the classical return mapping schemes of the infinitesimal theory. Comput Methods Appl Mech Eng; 1992; 99, 1 pp. 61-112.1182873 [DOI: https://dx.doi.org/10.1016/0045-7825(92)90123-2] 0764.73089

17. Khan, AS; Huang, S. Continuum theory of plasticity; 1995; New York, WileyInterscience: 440. [DOI: https://dx.doi.org/10.1021/jp108806v] 0856.73002

18. Nedjar, B. Frameworks for finite strain viscoelastic-plasticity based on multiplicative decompositions. Part II: computational aspects. Comput Methods Appl Mech Eng; 2002; 191, 15 pp. 1563-1593. [DOI: https://dx.doi.org/10.1016/S0045-7825(01)00336-X] 1141.74318

19. Kröner, E. Allgemeine kontinuumstheorie der versetzungen und eigenspannungen. Arch Ration Mech Anal; 1960; 4, 4 pp. 273-334.1269780090.17601

20. Lee, EH. Elastic–plastic deformation at finite strains. J Appl Mech; 1969; 36, 1 pp. 1-6.3331838 [DOI: https://dx.doi.org/10.1115/1.3564580] 0179.55603

21. Dettmer, W; Reese, S. On the theoretical and numerical modelling of Armstrong–Frederick kinematic hardening in the finite strain regime. Comput Methods Appl Mech Eng; 2004; 193, 1 pp. 87-116. [DOI: https://dx.doi.org/10.1016/j.cma.2003.09.005] 1063.74020

22. Garino, CG; Vairo, MR; Fagés, SA; Mirasso, A; Ponthot, J-P. Numerical simulation of finite strain viscoplastic problems. J Comput Appl Math; 2013; 246, pp. 174-184.3018106 [DOI: https://dx.doi.org/10.1016/j.cam.2012.10.008] 1426.74057

23. Petiteau, J-C; Verron, E; Othman, R; Le Sourne, H; Sigrist, J-F; Barras, G. Large strain rate-dependent response of elastomers at different strain rates: convolution integral vs. internal variable formulations. Mech Time-Dependent Mater; 2013; 17, 3 pp. 349-367. [DOI: https://dx.doi.org/10.1007/s11043-012-9188-7]

24. Zopf, C; Garcia, MA; Kaliske, M. A continuum mechanical approach to model asphalt. Int J Pavement Eng; 2015; 16, 2 pp. 105-124. [DOI: https://dx.doi.org/10.1080/10298436.2014.927065] 1425.74052

25. Latorre, M; Montáns, F. A new class of plastic flow evolution equations for anisotropic multiplicative elastoplasticity based on the notion of a corrector elastic strain rate. Appl Math Modell; 2018; 55, pp. 716-740.3739883 [DOI: https://dx.doi.org/10.1016/j.apm.2017.11.003] 1480.74031

26. Khan, AS; Zhang, H. Finite deformation of a polymer: Experiments and modeling. Int J Plast; 2001; 17, pp. 1167-1188. [DOI: https://dx.doi.org/10.1016/S0749-6419(00)00073-5] 1097.74501

27. Brünig, M. Numerical analysis and large strain elastic–viscoplastic behavior of hydrostatic stress-sensitive metals. Int J Solids Struct; 2001; 38, pp. 635-656. [DOI: https://dx.doi.org/10.1016/S0020-7683(00)00092-5] 0984.74016

28. Lehmann, T. Some remarks on the decomposition of deformations and mechanical work. Int J Eng Sci; 1982; 20, pp. 281-288. [DOI: https://dx.doi.org/10.1016/0020-7225(82)90023-4] 0472.73034

29. Lehmann, T. Some theoretical considerations and experimental results concerning elastic-plastic stress-strain relations. Arch Appl Mech; 1982; 52, pp. 391-403.0496.73034

30. Brepols, T; Vladimirov, IN; Reese, S. Numerical comparison of isotropic hypo- and hyperelastic-based plasticity models with application to industrial forming processes. Int J Plast; 2014; 63, pp. 18-48. [DOI: https://dx.doi.org/10.1016/j.ijplas.2014.06.003] 1507.74208

31. Sansour, C; Bednarczyk, H. A study on rate-type constitutive equations and the existence of a free energy function. Acta Mech; 1993; 100, pp. 205-221.1238791 [DOI: https://dx.doi.org/10.1007/BF01174790] 0777.73003

32. Argyris, JH; Kleiber, M. Incremental formulation in nonlinear mechanics and large strain elasto-plasticity—natural approach. Part 1. Comput Methods Appl Mech Eng; 1977; 11, 2 pp. 215-247. [DOI: https://dx.doi.org/10.1016/0045-7825(77)90060-3] 0384.73053

33. Atluri, SN. On constitutive relations at finite strain: hypo-elasticity and elasto-plasticity with isotropic or kinematic hardening. Comput Methods Appl Mech Eng; 1984; 43, 2 pp. 137-171. [DOI: https://dx.doi.org/10.1016/0045-7825(84)90002-1] 0571.73001

34. Hughes, TJR; Winget, J. Finite rotation effects in numerical integration of rate constitutive equations arising in large-deformation analysis. Int J Numer Methods Eng; 1980; 15, 12 pp. 1862-1867.601532 [DOI: https://dx.doi.org/10.1002/nme.1620151210] 0463.73081

35. Kojić, M; Bathe, KJ. Studies of finite element procedures—stress solution of a closed elastic strain path with stretching and shearing using the updated Lagrangian Jaumann formulation. Comput Struct; 1987; 26, 1–2 pp. 175-179. [DOI: https://dx.doi.org/10.1016/0045-7949(87)90247-1] 0609.73074

36. Bruhns, OT; Xiao, H; Meyers, A. Self-consistent Eulerian rate type elasto-plasticity models based upon the logarithmic stress rate. Int J Plast; 1999; 15, 5 pp. 479-520. [DOI: https://dx.doi.org/10.1016/S0749-6419(99)00003-0] 1036.74010

37. Giorgi, C; Morro, A. A thermodynamic approach to rate-type models of elastic-plastic materials. J Elast; 2021; 147, pp. 113-148.4368613 [DOI: https://dx.doi.org/10.1007/s10659-021-09871-3] 1489.74002

38. Mandel J (1971) Plasticite classique et viscoplasticity. In: CISM course 97. Springer, Udine

39. Zhang, M; Montán, FJ. A simple formulation for large-strain cyclic hyperelasto-plasticity using elastic correctors. Theory and algorithmic implementation. Int J Plast; 2019; 113, pp. 185-217. [DOI: https://dx.doi.org/10.1016/j.ijplas.2018.09.013]

40. Simo, J; Ortiz, M. A unified approach to finite deformation elastoplastic analysis based on the use of hyperelastic constitutive equations. Comput Methods Appl Mech Eng; 1985; 49, pp. 221-245. [DOI: https://dx.doi.org/10.1016/0045-7825(85)90061-1] 0566.73035

41. Weber, G; Anand, L. Finite deformation constitutive equations and a time integration procedure for isotropic, hyperelastic–viscoplastic solids. Comput Methods Appl Mech Eng; 1990; 79, pp. 173-202. [DOI: https://dx.doi.org/10.1016/0045-7825(90)90131-5] 0731.73031

42. Eterovic, AL; Bathe, KJ. A hyperelastic-based large strain elasto-plastic constitutive formulation with combined isotropic-kinematic hardening using the logarithmic stress and strain measures. Int J Numer Methods Eng; 1990; 30, pp. 1099-1114. [DOI: https://dx.doi.org/10.1002/nme.1620300602] 0714.73035

43. Flory, PJ. Thermodynamic relations for high elastic materials. Trans Faraday Soc; 1961; 57, pp. 829-838.128117 [DOI: https://dx.doi.org/10.1039/tf9615700829] 1380.81315

44. Ogden, RW. Nearly isochoric elastic deformations: application to rubberlike solids. J Mech Phys Solids; 1978; 26, pp. 37-57.505905 [DOI: https://dx.doi.org/10.1016/0022-5096(78)90012-1] 0377.73044

45. Holzapfel, GA. Nonlinear solid mechanics. A continuum approach for engineering; 2000; Chichester, Wiley: 0980.74001

46. Holzapfel, GA; Gasser, TC; Stadler, M. A structural model for the viscoelastic behavior of arterial walls: continuum formulation and finite element analysis. Eur J Mech A/Solids; 2002; 21, 2002 pp. 441-463. [DOI: https://dx.doi.org/10.1016/S0997-7538(01)01206-2] 1100.74597

47. Holzapfel, GA; Gasser, TC. A viscoelastic model for fiber-reinforced composites at finite strains: continuum basis, computational aspects and applications. Comput Methods Appl Mech Eng; 2001; 190, pp. 4379-4403. [DOI: https://dx.doi.org/10.1016/S0045-7825(00)00323-6] 1146.74342

48. Coda, HB. A finite strain elastoplastic model based on Flory's decomposition and 3D FEM applications. Comput Mech; 2022; 69, pp. 245-266.4373425 [DOI: https://dx.doi.org/10.1007/s00466-021-02092-4] 07492668

49. Hayashi, EY; Coda, HB. Alternative finite strain viscoelastic models: constant and strain rate-dependent viscosity. Acta Mech; 2024; 235, pp. 3699-3719.4752952 [DOI: https://dx.doi.org/10.1007/s00707-024-03914-1] 07866249

50. Coda, HB; Sanches, RAK. Unified solid–fluid Lagrangian FEM model derived from hyperelastic considerations. Acta Mecanica; 2022; 233, 7 pp. 2653-2685.4448975 [DOI: https://dx.doi.org/10.1007/s00707-022-03237-z] 1493.74107

51. Coda, HB; Paccola, RR. An alternative positional FEM formulation for geometrically non-linear analysis of shells: curved triangular isoparametric elements. Comput Mech; 2007; 40, pp. 185-200. [DOI: https://dx.doi.org/10.1007/s00466-006-0094-1] 1178.74163

52. Rivlin, R; Saunders, D. Large elastic deformations of isotropic materials VII. Experiments on the deformation of rubber. Philos Trans R Soc Lond Ser A; 1951; 243, pp. 251-288. [DOI: https://dx.doi.org/10.1098/rsta.1951.0004] 0042.42505

53. Hartmann, S; Neff, P. Polyconvexity of generalized polynomial-type hyperelastic strain energy functions for near-incompressibility. Int J Solids Struct; 2003; 40, pp. 2767-2791.2003053 [DOI: https://dx.doi.org/10.1016/S0020-7683(03)00086-6] 1051.74539

54. Simo J, Hughes T (2000) Computational inelasticity. In: Interdisciplinary applied mathematics. Springer, New York

55. Chaboche, JL. A review of some plasticity and viscoplasticity constitutive theories. Int J Plast; 2008; 24, pp. 1642-1693. [DOI: https://dx.doi.org/10.1016/j.ijplas.2008.03.009] 1142.74012

56. Amin, AFMS; Lion, A; Sekita, S; Okui, Y. Nonlinear dependence of viscosity in modeling the rate-dependent response of natural and high damping rubbers in compression and shear: experimental identification and numerical verification. Int J Plast; 2006; 22, pp. 1610-1657. [DOI: https://dx.doi.org/10.1016/j.ijplas.2005.09.005] 1146.74310

57. Wang, P; Chalal, H; Abed-Meraim, F. Quadratic solid–shell elements for nonlinear structural analysis and sheet metal forming simulation. Comput Mech; 2017; 59, pp. 161-186.3603805 [DOI: https://dx.doi.org/10.1007/s00466-016-1341-8] 1293.74415

58. Wriggers, P; Eberlein, R; Reese, S. A comparison of three-dimensional continuum and shell elements for finite plasticity. Int J Solids Struct; 1996; 33, pp. 3309-3326. [DOI: https://dx.doi.org/10.1016/0020-7683(95)00262-6] 0913.73071

59. Eberlein, R; Wriggers, P. Finite element concepts for finite elastoplastic strains and isotropic stress response in shells: theoretical and computational analysis. Comput Methods Appl Mech Eng; 1999; 171, pp. 243-279. [DOI: https://dx.doi.org/10.1016/S0045-7825(98)00212-6] 0957.74047