Finite element implementation of a certain class of elasto-viscoplastic constitutive models

In this paper a selected type of elasto-viscoplastic constitutive equations is considered. The viscoplastic model which was proposed by Marquis is generalized by allowing multiple terms describing the isotropic and the kinematic hardening. Furthermore, the presented model formulation enables one to use an arbitrary equivalent plastic strain function to describe the isotropic hardening behavior. What is more, the backstress evolution equation which was proposed by Marquis was modified so that any equivalent plastic strain function can be used in the recovery term now. The general form of the generalized constitutive model obtained this way was subsequently implemented into the finite element method (FEM). For that purpose, the radial-return mapping algorithm was utilized. The consistent tangent operator was derived for the considered class of models and is presented in the paper. A developed user material subroutine (UMAT) which allows one to use the viscoplastic models under consideration in the FEM program CalculiX is attached in the appendix section. A number of numerical simulations were conducted in order to verify the performance of the developed UMAT code.


Introduction
The classical small strain plasticity theory finds a wide range of applications in the engineering analysis. This theory is based on the Huber-von Mises-Hencky (HMH) yield criterion and is most commonly utilized to capture the mechanical properties of metals, e.g., [22]. One of the simplest versions of this theory assumes the so-called isotropic hardening rule, i.e., the expansion of the yield surface in the stress space which occurs during the plastic straining. Linear, piecewise linear and nonlinear relations are used to describe the evolution of the yield stress. A slightly more sophisticated approach assumes the usage of kinematic hardening rule, i.e., the translation of the yield surface which occurs during straining. The so-called backstress variable which defines the translation of the yield surface center is, in the simplest case, taken to evolve according to the linear equation proposed by Prager [17]. Furthermore, a group of mixed-hardening models exists which combine the two aforementioned behaviors. The constitutive models listed above are available in practically every software which is utilized for engineering analysis, e.g., [6,9,11,15]. In the case of cyclic loadings, the nonlinear kinematic hardening rules, such as the Armstrong-Frederick (A-F) rule [1], allow to achieve a substantially better agreement between the model predictions and the experimental results.
In order to capture such effects as strain rate sensitivity, creep or stress relaxation, a flow plasticity model must be generalized by formulating the so-called viscoplastic model. One of the classical formulations in this field is the constitutive model by Perzyna [16] and presents itself an extension of the flow plasticity model with isotropic hardening. A viscoplastic model assuming the mixed-hardening rule was proposed by Chaboche [5]. This model is a generalization of the previously proposed Chaboche-Rousselier elastoplastic model [3,4] and is more suitable for predicting material response when subjected to cyclic loadings. One possible extension of the Chaboche viscoplastic model was introduced in the work by Benallal and Marquis [2]. The proposed model assumes utilizing a nonlinear kinematic hardening rule which was obtained by modifying the A-F equation. A function of the equivalent plastic strain was added to the so-called recovery term in the backstress evolution equation. In the original formulation by Marquis, this multiplier was taken in the form of an exponential law. Furthermore, it was assumed that the kinematic hardening behavior is captured by a single backstress variable, whereas the isotropic hardening is governed by the nonlinear rule developed by Voce [21]. The viscous effects follow a nonlinear power law.
In this study a generalization of the viscoplastic model by Marquis is proposed. It is assumed that the kinematic hardening can be described by any defined number of backstress variables, whereas the backstress evolution law proposed by Marquis is extended by allowing the usage of any selected function of the equivalent plastic strain as the additional multiplier in the evolution law's recovery term. Furthermore, the model is generalized by allowing the isotropic hardening behavior to be described by arbitrary function of the equivalent plastic strain.
In the case of more advanced constitutive equations, the analytical solutions to boundary value problems are available for a limited number of cases only. Thus, in order to obtain solutions it is necessary to utilize numerical methods such as the finite element method (FEM), for instance. For that purpose, a numerical calculations algorithm has to be developed for the considered constitutive model.
In the work by Zienkiewicz and Cormeau [24], an algorithm for the finite element (FE) implementation of elasto-viscoplastic constitutive equations was proposed. Different variants of the associated and non-associated flow rule were considered. The cases of HMH and Tresca yield conditions with the assumption of ideal plasticity were analyzed. Furthermore, in one of the proposed model formulations the HMH yield condition with isotropic hardening was used. The aforementioned constitutive equations were utilized to solve several test problems using the introduced model integration algorithm. Moreover, an extended version of the constitutive model was proposed which included additional terms accounting for the viscoelastic effects. These studies were continued as reported in [25] by analyzing other viscoplastic model formulations which would utilize the Drucker-Prager or the Coulomb-Mohr yield conditions. A review of the elastoplastic constitutive models along with their integration using the radial-return mapping algorithm was presented in [7] with the main focus on the models utilizing the isotropic hardening. The Newton-Raphson (N-R) method was used for solving the nonlinear equations which result from the algorithm. The case of nonlinear stress-strain relation in the elastic range was considered as well. In [8] a modified version of the radial-return algorithm was proposed which is correct for any defined number of internals variables which simulate the loading history effects. The nonlinear equation system which follows from the applied model discretization algorithm was solved using the N-R method. The proposed methodology was utilized to implement into FEM various elastoplastic and elasto-viscoplastic constitutive models. In [10] the FE implementation of two alternative viscoplastic models was discussed, i.e., the Perzyna and the so-called consistency model. A comparative study was conducted.
In the paper by Yang and Feng [23], a viscoplastic model was considered in which additive decomposition of the plastic strain was assumed, i.e., into the steady state and the viscoplastic parts. A numerical integration algorithm was developed for the constitutive model under consideration which was subsequently used to implement the model into the FE analysis program ABAQUS.
Simo and Taylor [19] used the radial-return mapping algorithm to derive the consistent tangent operators for elastoplastic models with mixed hardening behavior. In particular, the nonlinear isotropic hardening was assumed, whereas the kinematic hardening was simulated with a modified linear Prager rule. In the previous work by Suchocki [20], the FE implementation of the cyclic elastoplasticity was discussed. In particular, some generalizations of the Chaboche-Rousselier elastoplastic model were analyzed. Kullig and Wippler [14] developed a numerical integration algorithm for the viscoplastic Chaboche model including a static recovery term in the backstress evolution equation.
In this work the FE implementation of the proposed Marquis viscoplastic model generalization is discussed. The considered constitutive equation makes use of the modified A-F rule to capture the kinematic hardening behavior. Implementing this particular model formulation into FEM has not been thoroughly discussed in the literature yet. The radial-return mapping algorithm was used to integrate the considered constitutive model. A general formula for the fourth-order consistent tangent operator was derived and is presented in the paper. To the best of author's knowledge, the specific forms of the tangent operator and the stress update procedure obtained for the considered class of viscoplastic constitutive models have not been reported in the literature yet. The derived operator along with the relations for the stress were utilized to implement the model into the non-commercial, open-source FE program CalculiX. This was achieved by taking advantage of UMAT (UserMATerial) subroutine option which is provided by CalculiX, c.f. [9]. In order to verify the performance of the developed UMAT subroutine, numerous FE simulations were performed. Selected results are included in this work. In the case of processes with the homogeneous stress and strain fields, such as uniaxial tension/compression (UT/UC), equibiaxial tension/compression (BT/BC) and simple shear (SS), the developed discrete version of the model was used to derive equation systems which describe the aforementioned deformations. The obtained equations were subsequently utilized to write proper programs in Scilab [18]. These programs allow one to calculate the material's stress response for a given deformation process and strain history. Thus, the homogeneous deformations previously analyzed in CalculiX using the developed UMAT code were simulated again in Scilab with the prepared scripts. The results produced using these two different methods were compared to verify the performance of the UMAT subroutine. A perfect agreement was found between the results generated by CalculiX and those obtained in Scilab. Additionally, some results achieved when analyzing stress relaxation, creep and nonhomogeneous deformations are included. A notched flat bar subjected to cyclic biaxial tension/compression and a thick-walled tube loaded by internal pressure were analyzed. The developed UMAT subroutine's performance has been tested for different types of finite elements. The code is attached in the appendix section.

Basic notions
The total stress tensor σ is taken to be the sum of the volumetric stress p and the stress deviator s, i.e., with e = tr(ε e ) and e e = ε e − 1 3 tr(ε e )1, where ε e is the small elastic strain tensor, 1 is the identity tensor, whereas tr(•) denotes the trace operation. B and μ are the bulk and the shear elastic moduli, respectively. It is assumed that the total small strain tensor can be additively decomposed into the elastic and the plastic components, i.e., ε = ε e + ε p .
Thus, in order to calculate the stress for a given strain it is, according to Eqs. (1) and (2) necessary to know the amount of plastic strain. The plastic strain tensor is determined utilizing the associated flow rule with the assumption of plastic incompressibility, i.e., tr(ε p ) = 0. The following form of the yield surface is assumed which takes into account both the isotropic and the kinematic hardening behaviors: where J 2 (s − X) = 3 2 (s − X) · (s − X) with X being the so-called backstress tensor. The initial yield stress is denoted as k, whereas R(ē p ) is the isotropic hardening function withē p being the equivalent plastic strain defined as:ē The flow rule associated with the yield surface given by Eq. (3) takes the form: withn being the normalized effective stress. The effective plastic strain rate is assumed to be given by the following relation:ė where • are Macaulay brackets while K and m are the viscoplasticity constants while ζ = 1 [1/s] is the viscosity parameter. In general, the isotropic hardening (or softening) behavior can be governed by any chosen function R(ē p ). However, this function is commonly assumed in the form of a series of terms proposed by Voce [21], i.e., with Q i and b i (i = 1, 2, . . . , N ) being the material parameters. In the original formulation presented in [2], a single backstress tensor X was used. However, the total backstress can be taken to be the sum of multiple components, i.e., In [3][4][5] the backstress variables were assumed to follow the nonlinear evolution law developed by Armstrong-Frederick, c.f. [1], i.e.,Ẋ where C i and γ i (i = 1, 2, . . . , M) are material parameters. In the paper by Benallal and Marquis [2], a modified version of the Armstrong-Frederick (A-F) evolution law was used. The modification concerns introducing a plastic strain function ϕ(ē p ) in the equation's "recovery" term, i.e., with the following specific form of ϕ(ē p ) used in [2]: where ϕ ∞ ∈ 0; 1 andω are material constants. This concept can be generalized for the total backstress being a sum of M components as given in Eqs. (8) and (9), i.e., with and ϕ ∞ i ∈ 0; 1 . In the following derivations the general case is considered, i.e., the forms of the functions R(ē p ) and ϕ(ē p ) remain unspecified.

Model discretization
For the purpose of discretizing the constitutive model described in the previous section, the radial-return mapping algorithm is utilized. For each computational step n + 1, the strain increment is initially assumed to be purely elastic. Thus, the deviatoric stress is calculated as: with the stress predictor s is satisfied, then no yielding takes place and so s n+1 = s pr n+1 andē p n+1 =ē p n . If the condition is not satisfied the plastic flow occurs. The equivalent plastic strain is given asē p n+1 =ē p n + ē p where the plastic strain increment ē p has to be determined so that the other variables, i.e., s n+1 , X n+1 and e p n+1 , can be updated. To this end, Eq. (12) is discretized using the backward Euler method, c.f. [14,20]: with Subtracting X n+1 from Eq. (14) and substituting e p = 3 2 ē pn after some manipulations leads to: or, after using Eq. (5) 2 : During the plastic flow, the discretized form of Eq. (6) is given as: Equation (19) can be solved for J 2 (s n+1 − X n+1 ), i.e., Thus, It follows from Eq. (21) 1 that: Thus, the following scalar, nonlinear equation is obtained which can be utilized to compute the effective plastic strain increment ē p : Equation (23) can be solved iteratively using the Newton-Raphson method. Thus, during every step of the iterative process a correction term is calculated: which is subsequently used to update the equivalent plastic strain increment, i.e., ē p i+1 = ē p i + c i+1 . The derivative of r ( ē p ) with respect to ē p takes the form: where and After computing ē p , the backstress components can be updated according to Eq. (15) and subsequently the total deviatoric stress can be calculated, i.e., The presented integration algorithm is used to determine the consistent tangent operator which is required in order to implement the considered model into FEM.

Consistent tangent operator
Taking a variation of the deviatoric stress with respect to all quantities leads to: It is seen that the variations δē p and δn n+1 are necessary in order to determine the formula for the consistent tangent operator. The term δē p can be found by taking the variation on Eq. (23), i.e., The variations δh( ē p ), δ J 2 (Z) and δα( ē p ) need to be found so that Eq. (31) could be solved for δē p . The variation δh( ē p ) is determined by taking a variation on Eq. (20), thus, In order to find δ J 2 (Z), a variation must be taken on Eq. (17) 2 which results in: with Thus, where the following notation is introduced: Taking the variation on Eq. (21) 2 leads to the formula for δα( ē p ): After inserting Eqs. (32), (35) and (37) into Eq. (31), it can be solved for δē p ; thus, where the so-called effective hardening modulus has been introduced: It can be shown thatn n+1 = Z J 2 (Z) , e.g., [20]. Taking a variation on this relation δn n+1 can be found, i.e., Substituting Eqs. (38) and (40) into Eq. (30), after some rearrangements and taking into account the volumetric stress component, leads to a linear relation between the variations of stress and strain, i.e., with the consistent viscoplastic tangent operator given as: where I is the fourth order identity tensor, whereas the introduced effective Lamé constants are given by the following relations: where the Lamé parameter λ = B − 2 3 μ.

User material (UMAT) subroutine for cyclic elasto-viscoplasticity
Many commercial and non-commercial FE packages provide some interfaces enabling the users to code their own constitutive laws, e.g., [6,9,11,15]. The open-source FE analysis program CalculiX offers two options for writing a user subroutine UMAT (UserMATerial), i.e., the ABAQUS UMAT or the CalculiX native UMAT interface, c.f. [9]. Usually, the native interface proves to be more efficient than ABAQUS UMAT used under CalculiX. Thus, it was decided to code the constitutive model under study using CalculiX UMAT. Some details regarding the interaction of UMAT subroutine with CalculiX can be found in [9] or [20]. The FE solver provided by CalculiX assumes that the tangent operator is symmetric with 21 independent components only. Thus, the determined operator given by Eq. (42) was symmetrized for the purpose of implementing the considered viscoplastic model via CalculiX UMAT. It has been demonstrated that performing symmetrization does not affect the accuracy of computations and has very little influence on their robustness, e.g., [12,13,20].

Exemplary problems
In order to verify the performance of the constitutive model under consideration, a number of FE simulations have been conducted. These included simulations of processes with homogeneous stress and strain fields, i.e., uniaxial tension/compression (UT/UC), equibiaxial tension/compression (BT/BC) and simple shear (SS). The discretized version of the model which is discussed in Sect. 3 was used do derive the process equations for the aforementioned homogeneous deformations. These equations were subsequently used to write Scilab programs [18]. Each program can be utilized to calculate the stress response for a given strain history. The results generated in Scilab were used to verify the computations performed in CalculiX with the use of the written UMAT code. In addition, some more complicated problems were analyzed, i.e., stress relaxation test, creep test, notched flat bar in equibiaxial tension/compression and thick-walled tube loaded by internal pressure.
The material parameters gathered in Table 1 were utilized for the conducted simulations. Thus, a single backstress variable (M = 1) and the function ϕ(ē p ) form as given by Eq. (11) were used, whereas the isotropic hardening behavior was simulated using a single Voce term, i.e., N = 1 in Eq. (7). The Lamé parameters are calculated as: λ = Eν (1+ν) (1−2ν) , μ = E 2(1+ν) , whereas the bulk modulus B = E 3(1−2ν) . The results presented in the following paragraphs were obtained for C3D8 elements. 1 The developed UMAT procedure has been tested for other types of finite elements as well. Kinematic hardening parameter 10 -

Uniaxial tension/compression (UT/UC)
A one-dimensional stress state is considered, i.e., In each increment the axial stress predictor σ pr n+1 is calculated for the given axial strain increment ε a : Subsequently, the yield criterion is checked. If σ pr n+1 − X n < k + R(ē p n ), the material response is elastic and σ n+1 = σ pr n+1 is assumed. Otherwise, the plastic flow occurs and the plastic strain increment ē p must be determined. This is done by solving Eq. (23) which in this particular case takes the form: with h( ē p ) being specified by Eq. (20), whereas: The fsolve function offered by Scilab was used for solving Eq. (47). After obtaining ē p , the backstress values are updated, i.e., where sgn(•) denotes the signum function, while w i (i = 1, 2, . . . , M) is given by Eq. (16). Finally, the total axial stress is updated: The computations according to Eqs. (46-50) were performed in Scilab for every single increment of the performed simulations.
In CalculiX the UT/UC FE simulations were conducted on a 1 mm×1 mm×1 mm cube meshed with a single C3D8 element. The considered deformation along with the applied boundary conditions can be seen in Fig. 1a. A kinematic loading in the form of a displacement δ in the direction "1" (Fig. 1a,e) was applied to the cube's ABC D face. Furthermore, a zero displacement in direction "1" was set on the face E FG H, a zero displacement in direction "2" was set on the face AE H D, and zero displacement in direction "3" was set on face AB F E. In the first approach a ramp tension test was considered with the maximum displacement value δ = 0.005 mm. The kinematic excitation was applied with different deformation rates, i.e., 1 × 10 −6 s −1 , 1 × 10 −2 s −1 and 1 s −1 . The simulations were performed in CalculiX utilizing the developed UMAT code and later repeated in Scilab using the program described above. An excellent agreement between the results obtained in CalculiX and those produced by Scilab was observed (Fig. 2).
In the second approach a cyclic deformation of the cube was analyzed. The kinematic excitation was defined using a saw function with the amplitude of δ a = 0.005 mm and the period T = 0.2 s. A preloading to the maximum value of δ a was defined which was followed by 10 full deformation cycles. Again, for the verification purposes, the considered deformation process was analyzed both in CalculiX and Scilab. As previously, the obtained results were found to be in an excellent agreement (Fig. 3). Both simulations were later repeated for the cube meshed with 125 C3D8 elements with the same results as previously.

Equibiaxial tension/compression (BT/BC)
A two-dimensional stress state is considered, i.e., consequently: The total and plastic strain tensors have the following components: with ε a and ε l being the axial and lateral strains, respectively. The axial stress in the increment n + 1 can be expressed as follows: If the condition σ pr n+1 − X n < k + R(ē p n ) is satisfied, the predictor stress value is taken to be the total updated axial stress, i.e., σ n+1 = σ pr n+1 . Otherwise, the equivalent plastic strain increment ē p needs to be determined by solving the following equation: where and where w i (i = 1, 2, . . . , M) is given by Eq. (16). After solving Eq. (55), the backstress variables can be updated, i.e., and subsequently the total axial stress value can be computed: The calculations according to Eqs. (54-59) are performed in every increment of the analysis.
The FE simulation in CalculiX was performed on a 1 mm×1 mm×1 mm cube meshed with a single C3D8 element. The boundary conditions which were applied to the cube are depicted in Fig. 1b. A kinematic excitation was assumed in the form of a displacement δ set on two faces, i.e., the displacement in direction "1" was applied to the cube's face ABC D, whereas the displacement in direction "2" was applied to the face B FGC (Fig. 1b,e). A zero displacement in direction "1" was defined on the face E FG H, while a zero displacement in direction "2" was set on the face AE H D, and finally, a zero displacement in direction "3" was defined on the face AB F E. The time history of the displacement δ was assumed in the form of a saw function with the amplitude δ a = 0.005 mm and the period T = 20 s. Consequently, the strain rate in each of the straining axes was equal to 1 × 10 −3 s −1 . The initial prestraining and 10 subsequent cycles were analyzed. Again, the process equations derived for the discretized version of the model were used to write a Scilab program which calculates the material stress response for a given strain history. This program was utilized to verify the results obtained with CalculiX. An excellent agreement was found (Fig. 4). The considered FE simulations were subsequently repeated for the cube meshed with 125 C3D8 elements with the same results.

Simple shear (SS)
In the case of simple shear (SS) process, the stress tensors are assumed to have the following components: with τ and X τ being the shear stress and the shear backstress component, respectively. For the analyzed process, the following scalar equations can be obtained: with τ pr n+1 being the shear predictor stress in the increment n + 1 while ε p s is the plastic shear strain increment, whereas ε s is the total shear strain increment. If √ 3 τ pr n+1 − X τ n < k+ R(ē p n ) it is assumed that τ n+1 = τ pr n+1 . Otherwise, the plastic flow occurs and the equivalent plastic strain increment ē p must be determined. This is done by solving the nonlinear equation: where while w i (i = 1, 2, . . . , M) is given by Eq. (16). After calculating the plastic strain increment, the stresses can be updated, i.e., and The computations listed in Eqs. (61-65) are performed for every analysis step. As before, the analyzed deformation process was simulated in CalculiX for a cube meshed with a single C3D8 element. The deformed shape and the applied boundary conditions can be seen in Fig. 1c, d with γ = 2ε s being the shear angle. A kinematic excitation was assumed in the form of the displacement δ in the direction "1" applied to the cube's face B FGC where the displacements in the directions "2" and "3" were set to zero. Due to the small deformations, it is assumed that δ/1 = tg (γ ) ≈ γ . The displacements in all directions were defined as zero on the face AE H D. A cyclic deformation was considered with the amplitude δ a = 0.01 mm and the period T = 20 s which resulted in the shear strain rate 1 × 10 −3 s −1 . A preloading and 10 subsequent cycles were analyzed. Again, the set of equations derived for the discretized version of the model was utilized to write a Scilab program which was used to verify the results obtained in CalculiX. It can be seen in Fig. 5 that once again excellent agreement was found.

Stress relaxation
A one-dimensional stress relaxation process was simulated. Again, the FE simulation was performed in Cal-culiX using a cube meshed with a single C3D8 element. The boundary conditions as described in Sect. 6.1 were used with a different strain history which acted as an excitation (Fig. 6a). The material stress response can be seen in Fig. 6b. In order to validate the obtained FE results, a Scilab program utilizing Eqs. (46-50) was used to simulate the considered deformation process again. It can be seen in Fig. 6b that the results obtained in CalculiX and Scilab are in a perfect agreement.

Creep
A creep test simulation was performed. In CalculiX, as before, a cube meshed with a C3D8 element was used. The defined boundary conditions are the same as described in Sect. 6.1 with the only difference being a stress excitation applied to the face ABCD instead of a displacement. The assumed stress history and the strain response can be seen in Fig. 7a and b, respectively. For the purpose of validating the obtained FE results, the axial strain history computed in CalculiX was used as an input for two different Scilab programs in order to check if the correct stress history will be recreated by each one of them. The first program utilizes the radial-return mapping algorithm which for the one-dimensional case leads to Eqs. (46-50). The second Scilab program partially takes advantage of the analytical formulas. In order to check the performance of the developed UMAT code, a notched flat bar was considered which was subjected to cyclic biaxial tension/compression. The used geometry along with the dimensions, the applied boundary conditions and the coordinate system can be seen in Fig. 8a. The bar was meshed with C3D8 elements with three elements defined along the bar's thickness which equaled to 0.1 mm (Fig. 8b). The kinematic excitation was defined on bar's outer faces in the form of displacement histories δ 1 and δ 2 (Fig. 8a). The displacement δ 1 was set as a ramp elongation until its maximum value of 0.0025 mm which is reached for the time instant t = 10 s. The displacement δ 2 was defined with a saw function with the amplitude δ 2a = 0.005 mm and the period T = 40 s. The total of 10 cycles were performed. The final contour maps obtained for the displacement magnitude, the Huber-von Mises-Hencky (HMH) equivalent stress and the equivalent total strain are shown in Fig. 9.    Thick-walled tube loaded with internal cyclic pressure: a displacement magnitude, b HMH equivalent stress. c equivalent total strain 6.7 Thick-walled tube loaded with internal pressure As another problem a thick-walled tube loaded with internal pressure was considered. The assumed geometry along with the boundary conditions and the coordinate system are depicted in Fig. 10a. The pressure loading was defined using a saw function with the pressure limits p max = 470 MPa, p min = 0 MPa and the period T = 20 s. The total of 10 cycles were analyzed. The tube was assumed to be 100 mm thick and meshed with C3D8 elements (Fig. 10b). The contour plots showing the final distributions of the displacement magnitude, the HMH stress and the equivalent total strain can be seen in Fig. 11.

Conclusions
The viscoplastic model proposed in the work by Benallal and Marquis [2] was analyzed. This model presents itself a certain generalization of the constitutive model which was introduced by Chaboche [3][4][5]. In this study the model by Benallal and Marquis was further generalized by assuming the isotropic hardening behavior to be governed by any selected function of the equivalent plastic strain R(ē p ), c.f. Eq. (3). Moreover, the strain function in the "recovery" term of the backstress evolution law was taken to be any chosen function of the plastic strain ϕ(ē p ), see Eq. (12). The model was also extended by allowing multiple backstress variables to be utilized.
The obtained generalized viscoplastic constitutive equation was subsequently discretized by using the radial-return mapping algorithm. A general formula for the consistent tangent operator was derived. According to the author's best knowledge, this formula has not been presented in the literature before. The discrete model was implemented into the non-commercial, open-source FE program CalculiX by writing a UMAT subroutine which is included in the appendix section.
In order to verify the performance of the written UMAT code, proper equation sets were derived from the discretized version of the model which describe selected homogeneous deformation processes. The obtained equations were in the next step used to write Scilab scripts which computed the stress responses for the deformation histories in question. The simulation results obtained in Scilab were found to be in an excellent agreement to those obtained by performing the FE simulations in CalculiX using the developed UMAT subroutine. What is more, the UMAT code was used to solve several other boundary value problems such as notched flat bar in biaxial tension/compression and thick-walled tube loaded with internal pressure.