An Effective Computational Approach for the Numerical Simulation of Elasto-/Viscoplastic Solid Materials

Abstract

The effects of different loading programs on the elasto-/viscoplastic behavior of rate-sensitive materials are analyzed with specific numerical examples. An appropriate solution scheme and a consistent tangent operator are applied which are capable of being adopted for general computational procedures. Numerical computations and results are reported which illustrate the rate-dependence of the elasto-/viscoplastic constitutive model in use. In the numerical analysis the loading is applied by increasing the pressure and accordingly a nondimensional loading program parameter is introduced. In the numerical results the significance of the loading program is thus emphasized with reference to the nonlinear response of the elasto-/viscoplastic material behavior of solids.

1. Introduction

The variational formulation of boundary value problems governed by linear and nonlinear operators is a problem which is of interest in many fields of theoretical mechanics. A comprehensive presentation is reported, for example, by Vainberg [1] and by Oden and Reddy [2]. Nevertheless in recent times an increasing interest has been devoted in the literature to the statement of variational principles able to extend to elastoplasticity the variational formulations holding in linear elasticity and to provide valid mathematical bases to the development of computational algorithms. These efforts are shown by precursor works due to Capurso [3] and to Capurso and Maier [4]. More recently contributions may be quoted due to Simo et al. [5], Simo [6], Maier and Novati [7], and Reddy and Martin [8].

The nonlinear and nonsmooth model problem in elasto-/viscoplasticity is governed by nonlinear field equations and internal and external constraint equations ruled by multivalued operators. In such cases the recourse to the tools of subdifferential calculus appears to be advantageous since it allows presenting a sound variational formulation of the evolutive elasto-/viscoplastic problem in a general and compact framework; see, for example, De Angelis [9].

In the present paper we focus attention on the elastic/viscoplastic model problem and the related integration procedure. In elastoplasticity the system of variational inequalities is associated with a class of return mapping algorithms based on the generalized midpoint rule. Application of this operator split methodology is based on an elastic prediction and a plastic correction phase; see, for example, Nagtegaal [10], Ortiz and Popov [11], and Simo and Taylor [12]. Considerations on the stability of the generalized midpoint rule integration algorithms are reported by Ortiz and Popov [11] and by Simo [6] for elastoplasticity and by Hughes and Taylor [13] and Simo and Govindjee [14] for viscoplasticity. The accuracy analysis and the use of generalized midpoint rule algorithms have been analyzed in detail by Ortiz and Popov [11]. The need for a connection between the local algorithmic integration of the constitutive equations and the global structural problem was first indicated by Nagtegaal [10] and then developed by Simo and Taylor [12] by the introduction of the consistent tangent operator which results in a quadratic rate of asymptotic convergence for the computational procedure.

The derivation of a class of return mapping algorithms associated with the system of variational inequalities in viscoplasticity is not a minor problem with respect to the rate-independent one. In this regard Zienkiewicz and Cormeau [15] discussed integration procedures and considered time step restrictions for the Euler forward difference method in quasi-static elasto/viscoplasticity. Hughes and Taylor [13] reconsidered the application of implicit methods by the use of an algorithmic procedure which requires the inversion of a compliance matrix.

Integration algorithms for viscoplastic models involving nonsmooth yield surfaces are reported by Simo et al. [5], while stability properties of algorithms are investigated by Simo [6] and Simo and Govindjee [14]. Significant attempts of integration procedures for viscoplastic models are also found in Ju [16], where a closed form expression of an algorithmic compliance tensor is derived. For the solution of stiff equations arising in low-rate-sensitive materials Perić [17] proposed a perturbation method, while the linearization of the algorithm is performed by a systematic application of the directional derivative formula. A procedure for general yield criteria was presented by Alfano et al. [18]. Models including nonlinear kinematic hardening behavior and integration procedures applied for complex material models were analyzed in detail by Chaboche and Cailletaud [19].

A face-based smoothed finite element method was discussed by Nguyen-Thoi et al. [20] to improve the accuracy and convergence rate of the existing standard finite element method for the solid mechanics problems. The method was further extended to more complicated viscoelastoplastic analyses of 3D solids using the von-Mises yield function and the Prandtl-Reuss flow rule. A node-based smoothed finite element method was also proposed by Nguyen-Thoi et al. [21] for the solid mechanics problems. The system stiffness matrix was computed using the smoothed strains over the smoothing domains associated with nodes of the element mesh. The material behavior included perfect viscoelastoplasticity and viscoelastoplasticity with isotropic hardening and linear kinematic hardening. A dual formulation with displacements and stresses as the main variables was performed. In the numerical procedure the stress variables were eliminated and the problem became only displacement-dependent. An edge-based smoothed finite element method using triangular elements was also recently proposed by Nguyen-Thoi et al. [22] to improve the accuracy and convergence rate of the existing standard finite element method for the elastic solid mechanics problems. In the paper the edge-based smoothed finite element method was extended to more complicated viscoelastoplastic analyses using the von-Mises yield function and the Prandtl-Reuss flow rule. A scaled boundary polygon formulation was furtherly proposed by Ooi et al. [23] to model elastoplastic material responses in structures. The polygons showed flexible mesh generation and accurate solutions capabilities, especially for problems with cracks and notches. Within this approach standard finite element procedures were used to formulate the stiffness matrix and the residual load vector. The nonlinear material constitutive matrix and the internal stresses were approximated locally in each polygon by a polynomial function.

In the present paper we refer to the class of elasto-/viscoplastic material behavior often denoted in the literature by rate-sensitive materials; see, for example, Naghdi and Murch [24], Perzyna [25], and Skrzypek and Hetnarski [26], where it is assumed that viscous effects are exhibited beyond the elastic range. The numerical resolution of the evolutive problem for the adopted constitutive model is performed by following a procedure initially suggested by Simo [27] and by extending it with the use of a tangent operator which encompasses the case of yield functions of arbitrary type. Another characteristical feature of the presented procedure is that it applies to different elasto-/viscoplastic constitutive models by suitably specializing the flow function of the adopted constitutive model.

In the numerical analysis the loading is performed by increasing the pressure and accordingly a nondimensional loading program parameter is introduced in order to account for the effects of the loading rate. In the numerical results the significance of the loading program and the significance of the rate of the prescribed loading are thus emphasized with reference to the nonlinear response of the elasto-/viscoplastic material behavior of solids. Numerical applications and computational results are finally reported for a typical benchmark problem in order to illustrate the effectiveness and robustness of the adopted integration procedure.

2. Continuum Problem of Evolution

Let $\Omega \subset {\mathfrak{R}}^n$ , 1 ≤ n ≤ 3, be the reference configuration of body $\mathcal{B}$ with particles identified by a position vector x ∈ Ω relative to a Cartesian coordinate system. We indicate with $\mathcal{T}\subset {\mathfrak{R}}_{+}$ the time interval of interest and with V the space of displacements, D the strain space, and S the dual stress space. Let us indicate with $\mathbf{u}:\Omega \times \mathcal{T}\to \mathbf{V}$ the displacement vector and with $\mathit{\sigma }:\Omega \times \mathcal{T}\to \mathbf{S}$ the stress tensor. The compatible strain tensor is defined as $\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathit{\epsilon }(\mathbf{u})={\nabla }^s(\mathbf{u}):\Omega \times \mathcal{T}\to \mathbf{D}$ , where ∇ ^s is the symmetric part of the gradient.

In the analysis we consider quasi-static deformations and viscous effects are assumed to show beyond the elastic range; see, for example, Naghdi and Murch [24]. Thus we refer to the class of rate-sensitive materials; see, for example, Skrzypek and Hetnarski [26], and for a survey account see Lemaitre and Chaboche [28]. We assume an additive decomposition of total strain ∊ into the elastic and inelastic parts. Accordingly, we denote by ∊^vp = ∊ − ∊ ^e the viscoplastic strain, where combined viscous and plastic effects are represented.

A dual pair of kinematic $\mathit{\alpha }=({\mathit{\alpha }}_{\text{kin}},{\alpha }_{\text{iso}})\in \mathbf{X}\times \mathfrak{R}$ and static $\mathbf{\chi }=({\mathbf{\chi }}_{\text{kin}},{\chi }_{\text{iso}})\in {\mathbf{X}}^{\prime }\times \mathfrak{R}$ internal variables are considered in order to represent hardening behavior, where, respectively, ${\alpha }_{\text{iso}}\in \mathfrak{R}$ and ${\chi }_{\text{iso}}\in \mathfrak{R}$ model isotropic hardening and α_kin ∈ X and χ_kin ∈ X′ model kinematic hardening, with X and X′ being dual spaces. Static and kinematic internal variables are associated with the relation χ = Hα, where H = diag ⁡ [H_kin, H_iso] is the hardening matrix.

The convex elastic domain $\mathcal{E}$ is assumed to be characterized by means of a convex yield function f(σ, χ_kin, χ_iso) so that

$\mathcal{E}\stackrel{\scriptscriptstyle\mathrm{def}}{=}\left\{\left(\mathit{\sigma },{\mathbf{\chi }}_{\text{kin}},{\chi }_{\text{iso}}\right)\in \mathbf{S}\times {\mathbf{X}}^{\prime }\times \mathfrak{R}:f\left(\mathit{\sigma },{\mathbf{\chi }}_{\text{kin}},{\chi }_{\text{iso}}\right)\le \mathrm{0}\right\}.$

(1)

A general class of viscoplastic hardening materials arises when the yield function is expressed in the following form:

$f\left(\mathit{\sigma },{\mathbf{\chi }}_{\text{kin}},{\chi }_{\text{iso}}\right)=F\left(\mathit{\sigma }-{\mathbf{\chi }}_{\text{kin}}\right)-\kappa \left({\chi }_{\text{iso}}\right),$

(2)

where κ(χ_iso) is a material hardening parameter.

In the sequel we resort to the compact notation introduced by Halphen and Nguyen [29] for the generalized standard material model, so that strains and kinematic internal variables, as well as the corresponding dual ones, are collected in suitably defined generalized variables:

$\begin{array}{c}\mathbf{E}=(\mathit{\epsilon },\mathbf{o}),\hspace{1em}\hspace{1em}{\mathbf{E}}^e=({\mathit{\epsilon }}^e,\mathit{\alpha }),\hspace{1em}\\ {\mathbf{E}}^{\text{vp}}=({\mathit{\epsilon }}^{\text{vp}},-\mathit{\alpha }),\hspace{1em}\hspace{1em}\mathbf{\Sigma }=(\mathit{\sigma },\mathbf{\chi }).\end{array}$

(3)

The generalized kinematic and static variables are, respectively, defined in the product spaces $\stackrel{\sim }{\mathbf{D}}=\mathbf{D}\times \mathbf{X}\times \mathfrak{R}$ and $\stackrel{\sim }{\mathbf{S}}=\mathbf{S}\times {\mathbf{X}}^{\prime }\times \mathfrak{R}$ .

By indicating with Π*(Γ) a viscoplastic convex potential, the evolutive problem in viscoplasticity is expressed as optimality condition of the convex potential function:

${\mathcal{L}}^{\text{vp}}\left(\mathbf{\Gamma }\right)\stackrel{\scriptscriptstyle\mathrm{def}}{=}-〈\mathbf{\Gamma },{\dot{\mathbf{E}}}^{\text{vp}}〉+{\Pi }^{*}\left(\mathbf{\Gamma }\right),$

(4)

where the symbol 〈·, ·〉 represents a nondegenerate bilinear form acting on dual spaces, $\mathbf{\Gamma }=(\mathit{\tau },\mathbf{q})\in \stackrel{\sim }{\mathbf{S}}$ indicates a generic generalized stress, and the relevant value at the solution is denoted by Σ = (σ, χ).

The optimality condition for the potential thus implies

$\mathrm{0}\in {\left[\partial {\mathcal{L}}^{\text{vp}}\left(\mathbf{\Gamma }\right)\right]}_{\left(\mathbf{\Sigma }\right)}⟺{\dot{\mathbf{E}}}^{\text{vp}}\in \partial {\Pi }^{*}\left(\mathbf{\Sigma }\right),$

(5)

which represents the evolutive law of the viscoplastic strain and of the kinematic internal variables. The above equation is expressed in subdifferential form (see, e.g., Hiriart-Urruty and Lemaréchal [30]), so that the multivaluedness of the viscoplastic evolutive problem is suitably taken into account.

The conjugate of the viscoplastic potential is defined as the viscoplastic dissipation and it is expressed by

${\mathcal{D}}^{\text{vp}}\left({\dot{\mathbf{E}}}^{\text{vp}}\right)=\underset{\mathbf{\Gamma }\in \tilde{\mathbf{S}}}{\sup }\left\{〈\mathbf{\Gamma },{\dot{\mathbf{E}}}^{\text{vp}}〉-{\Pi }^{*}\left(\mathbf{\Gamma }\right)\right\}.$

(6)

Consequently the evolutive law (5)₂ may be expressed in the equivalent inverse form:

$\mathbf{\Sigma }\in \partial {\mathcal{D}}^{\text{vp}}\left({\dot{\mathbf{E}}}^{\text{vp}}\right),$

(7)

and both (5)₂ and (7) are also equivalently expressed in the form of Fenchel's equality:

${\Pi }^{*}\left(\mathbf{\Sigma }\right)+{\mathcal{D}}^{\text{vp}}\left({\dot{\mathbf{E}}}^{\text{vp}}\right)=〈\mathbf{\Sigma },{\dot{\mathbf{E}}}^{\text{vp}}〉.$

(8)

By considering the classical penalty regularization procedure (Yosida [31]) the constrained optimization problem related to the plastic model is phrased as an unconstrained problem. Let g⁺ (x) be a penalty function of the constraint f(Γ) satisfying the conditions: g⁺ (x) is continuous in [0, ∞), g⁺ (x) ≥ 0 and convex in [0, ∞), g⁺ (x) = 0 if and only if x ≤ 0. In these assumptions a regularized form of the potential function has the expression

${\mathcal{L}}_{\eta }^{\text{vp}}\left(\mathbf{\Gamma }\right)\stackrel{\scriptscriptstyle\mathrm{def}}{=}-〈\mathbf{\Gamma },{\dot{\mathbf{E}}}^{\text{vp}}〉+\frac{\mathrm{1}}{\eta }{g}^{+}\left(f\left(\mathbf{\Gamma }\right)\right),$

(9)

where η > 0 is a penalty parameter, which in viscoplasticity has the meaning of a viscosity coefficient.

The viscoplastic dissipation is thus expressed in the regularized form:

${\mathcal{D}}_{\eta }^{\text{vp}}\left({\dot{\mathbf{E}}}^{\text{vp}}\right)=\underset{\mathbf{\Gamma }\in \tilde{\mathbf{S}}}{\sup }\left\{〈\mathbf{\Gamma },{\dot{\mathbf{E}}}^{\text{vp}}〉-\frac{\mathrm{1}}{\eta }{g}^{+}\left(f\left(\mathbf{\Gamma }\right)\right)\right\}.$

(10)

Accordingly the viscoplastic problem is considered in the literature as a penalty regularization of the plastic problem and the solution Σ_η of the regularized problem tends to the solution Σ of the constrained plastic problem for η → 0⁺ (Luenberger [32]).

Different expressions of the viscoplastic evolutive equations are obtained by specializing the penalty function. For example, when the penalty function g⁺ is chosen to be in the form

$g^{+}\stackrel{\scriptscriptstyle\mathrm{def}}{=}\left\{\begin{array}{cc}\frac{1}{2}{x}^{\mathrm{2}}\hfill & \text{for}\hspace{0.17em}\hspace{0.17em}x\ge \mathrm{0}\hfill \\ 0\hfill & \text{for}\hspace{0.17em}\hspace{0.17em}x<\mathrm{0},\hfill \end{array}\right)$

(11)

the derivative is given by dg⁺ (x)/dx = 〈x〉, where the MacAuley bracket 〈·〉 is defined as 〈x〉 = (x + |x|)/2. Nonlinear viscous effects are modelled by assuming a flow function Φ(x) such that dg⁺ (x)/dx = 〈Φ(x)〉. Accordingly, by recalling expression (9), the optimality condition of the unconstrained problem yields in the regularized form

$\mathrm{0}\in {\left[\partial {\mathcal{L}}_{\eta }^{\text{vp}}\left(\mathbf{\Gamma }\right)\right]}_{\left(\mathbf{\Sigma }\right)}⟺{\dot{\mathbf{E}}}^{\text{vp}}\in \frac{\mathrm{1}}{\eta }〈\Phi \left(f\left(\mathbf{\Sigma }\right)\right)〉\partial f\left(\mathbf{\Sigma }\right),$

(12)

which represents the evolutive law for the Perzyna [25] viscoplastic constitutive model expressed in subdifferential form. A standard choice of the flow function for linear viscous effects is Φ(f) = f. Other proposed expressions of the flow function for nonlinear viscous effects are reported, for example, by Skrzypek and Hetnarski [26].

The interpretation of the evolutive problem in viscoplasticity as optimality condition of a convex optimization problem appears to be convenient to supply the viscoplastic structural problem with a complete variational formulation (see, e.g., De Angelis [9, 33]) and it is ideally suited for the development of numerical algorithms in finite element applications.

3. Numerical Solution Procedure

Let the time interval of interest $\mathcal{T}\subset {\mathfrak{R}}_{+}$ be specified as [0, T]. At time t _n ∈ [0, T] the total and viscoplastic strain fields and the internal variables are assumed to be known. For a prescribed increment of the displacement field Δu, the basic problem is to update the unknown fields at time t_{n + 1} ∈ [0, T] in a manner consistent with the evolutive equation (12)₂, which for convenience is reported below in components:

$\begin{array}{c}{\dot{\mathit{\epsilon }}}^{\text{vp}}=\frac{\mathrm{1}}{\eta }〈\Phi \left(f\left(\mathit{\sigma },\mathbf{\chi }\right)\right)〉\frac{\partial f\left(\mathit{\sigma },\mathbf{\chi }\right)}{\partial \mathit{\sigma }},\\ -{\dot{\mathit{\alpha }}}_{\text{kin}}=\frac{\mathrm{1}}{\eta }〈\Phi \left(f\left(\mathit{\sigma },\mathbf{\chi }\right)\right)〉\frac{\partial f\left(\mathit{\sigma },\mathbf{\chi }\right)}{\partial {\mathbf{\chi }}_{\text{kin}}},\\ -{\dot{\alpha }}_{\text{iso}}=\frac{\mathrm{1}}{\eta }〈\Phi \left(f\left(\mathit{\sigma },\mathbf{\chi }\right)\right)〉\frac{\partial f\left(\mathit{\sigma },\mathbf{\chi }\right)}{\partial {\chi }_{\text{iso}}}.\end{array}$

(13)

A von-Mises yield criterion with linear hardening is considered in the form

$f\left(\mathit{\sigma },{\mathbf{\chi }}_{\text{kin}},{\chi }_{\text{iso}}\right)=‖\mathit{s}-{\mathbf{\chi }}_{\text{kin}}‖-\kappa \left({\chi }_{\text{iso}}\right)=‖\mathit{\xi }‖-\sqrt{\frac{2}{3}}\left({\sigma }_{yo}+{\chi }_{\text{iso}}\right)\le \mathrm{0},$

(14)

where s ≝ dev σ is the stress deviator, ξ ≝ s − χ_kin is the relative stress, $\kappa ({\chi }_{\text{iso}})=\sqrt{2/3}({\sigma }_{yo}+{\chi }_{\text{iso}})$ represents the current radius of the yield surface in the deviatoric plane, and σ _yo denotes the uniaxial yield stress of the virgin material. For simplicity of notation we also set β ≝ χ_kin for the back stress.

In the assumption of linear hardening behavior, the static internal variable related to isotropic hardening is specified as χ_iso = H_isoα_iso, where the dual kinematic internal variable α_iso is represented by the equivalent viscoplastic strain ${\stackrel{-}{e}}^{\text{vp}}\stackrel{\scriptscriptstyle\mathrm{def}}{=}{\int }_{\mathrm{0}}^t\sqrt{2/3}‖{\dot{\mathit{\epsilon }}}^{\text{vp}}‖\text{d}t$ .

By setting the normal to the yield surface as n ≝ d _σ f = ξ/|ξ|, the evolutive laws may thus be expressed in components by

$\begin{array}{c}{\dot{\mathit{\epsilon }}}^{\text{vp}}=\frac{\mathrm{1}}{\eta }〈\Phi \left(f\left(\mathit{\sigma },\mathit{\chi }\right)\right)〉\frac{\mathit{\xi }}{‖\mathit{\xi }‖}=\frac{\mathrm{1}}{\eta }〈\Phi \left(f(\mathit{\sigma },\mathit{\chi })\right)〉\mathbf{n},\\ \dot{\mathit{\beta }}=\frac{2}{3}{H}_{\text{kin}}{\dot{\mathit{\epsilon }}}^{\text{vp}}=\frac{\mathrm{1}}{\eta }〈\Phi \left(f\left(\mathit{\sigma },\mathit{\chi }\right)\right)〉\frac{2}{3}{H}_{\text{kin}}\frac{\mathit{\xi }}{‖\mathit{\xi }‖}=\frac{\mathrm{1}}{\eta }〈\Phi \left(f(\mathit{\sigma },\mathit{\chi })\right)〉\frac{2}{3}{H}_{\text{kin}}\mathbf{n}\mathbf{,}\\ {\dot{\stackrel{-}{e}}}^{\text{vp}}=\frac{\mathrm{1}}{\eta }〈\Phi \left(f(\mathit{\sigma },\mathit{\chi })\right)〉\sqrt{\frac{2}{3}}.\end{array}$

(15)

A fully implicit integration scheme yields the following discrete forms of the evolutive equations:

$\begin{array}{c}{\mathit{\epsilon }}_{n+\mathrm{1}}^{\text{vp}}={\mathit{\epsilon }}_n^{\text{vp}}+{\gamma }_{n+\mathrm{1}}^{\text{vp}}{\mathbf{n}}_{n+\mathrm{1}},\hspace{1em}\hspace{1em}{\mathit{\beta }}_{n+\mathrm{1}}={\mathit{\beta }}_n+\frac{2}{3}{H}_{\text{kin}}{\gamma }_{n+\mathrm{1}}^{\text{vp}}{\mathbf{n}}_{n+\mathrm{1}},\\ {\stackrel{-}{e}}_{n+\mathrm{1}}^{\text{vp}}={\stackrel{-}{e}}_n^{\text{vp}}+\sqrt{\frac{2}{3}}{\gamma }_{n+\mathrm{1}}^{\text{vp}},\hspace{1em}\hspace{1em}{\kappa }_{n+\mathrm{1}}={\kappa }_n+\frac{2}{3}{H}_{\text{iso}}{\gamma }_{n+\mathrm{1}}^{\text{vp}},\end{array}$

(16)

where n_{n + 1} = d _σ f_{n + 1} = ξ_{n + 1}/|ξ_{n + 1}| and a time step dependent viscoplastic multiplier has been introduced as

${\gamma }_{n+\mathrm{1}}^{\text{vp}}=\frac{\mathrm{1}}{\eta }〈\Phi \left(f({\mathit{\sigma }}_{n+\mathrm{1}},{\mathit{\beta }}_{n+\mathrm{1}},{\kappa }_{n+\mathrm{1}})\right)〉\Delta t,$

(17)

with Δt = t_{n + 1} − t _n .

An elastic prediction-plastic correction scheme yields the following trial values:

${\mathit{s}}_{n+\mathrm{1}}^t=\mathrm{2}G\left({\mathit{\nu }}_{n+\mathrm{1}}-{\mathit{\nu }}_n^{\text{vp}}\right),\hspace{1em}\hspace{1em}{\mathit{\xi }}_{n+\mathrm{1}}^t={\mathit{s}}_{n+\mathrm{1}}^t-{\mathit{\beta }}_n,$

(18)

where ν ≝ dev ∊ is the strain deviator and G denotes the shear modulus. The stress deviator and the relative stress are accordingly expressed at time t_{n + 1} by

$\begin{array}{c}{\mathit{s}}_{n+\mathrm{1}}={\mathbf{s}}_{n+\mathrm{1}}^t-\mathrm{2}G{\gamma }_{n+\mathrm{1}}^{\text{vp}}{\mathbf{n}}_{n+\mathrm{1}},\\ {\mathit{\xi }}_{n+\mathrm{1}}={\mathit{\xi }}_{n+\mathrm{1}}^t-\left(\mathrm{2}G+\frac{2}{3}{H}_{\text{kin}}\right){\gamma }_{n+\mathrm{1}}^{\text{vp}}{\mathbf{n}}_{n+\mathrm{1}},\end{array}$

(19)

whereas the trial value of the yield function is given by

$f_{n+\mathrm{1}}^t\stackrel{\scriptscriptstyle\mathrm{def}}{=}‖{\mathit{\xi }}_{n+\mathrm{1}}^t‖-{\kappa }_n.$

(20)

If f_{n + 1} ^t ≤ 0, the step is elastic and the values of the unknowns at the end of the time step are set equal to the trial values. If f_{n + 1} ^t > 0, the viscoplastic multiplier γ_{n + 1}^vp and the normal n_{n + 1} are to be evaluated. By inverting relation (17) it follows that

$f_{n+\mathrm{1}}={\Phi }^{-\mathrm{1}}\left(\frac{\eta }{\Delta t}{\gamma }_{n+\mathrm{1}}^{\text{vp}}\right).$

(21)

By taking the dot product of relation (19)₂ with n_{n + 1} and recalling (14) and the above expressions of f_{n + 1} ^t and f_{n + 1} the following nonlinear equation is obtained:

$\left(‖{\mathit{\xi }}_{n+\mathrm{1}}^t‖-{\kappa }_n\right)-\left(\mathrm{2}G+\frac{2}{3}{H}_{\text{kin}}+\frac{2}{3}{H}_{\text{iso}}\right){\gamma }_{n+\mathrm{1}}^{\text{vp}}={\Phi }^{-\mathrm{1}}\left(\frac{\eta }{\Delta t}{\gamma }_{n+\mathrm{1}}^{\text{vp}}\right),$

(22)

which can be solved for γ_{n + 1}^vp by a Newton iteration scheme. In addition, at the end of the time step it is as follows:

${\mathbf{n}}_{n+\mathrm{1}}=\frac{{\mathit{\xi }}_{n+\mathrm{1}}^t}{‖{\mathit{\xi }}_{n+\mathrm{1}}^t‖}.$

(23)

The values of the variables are consequently updated according to (16) and (19), while the stress tensor is evaluated as

${\mathit{\sigma }}_{n+\mathrm{1}}={\mathbf{C}}^{\text{vp}}\left({\mathit{\epsilon }}_{n+\mathrm{1}}-{\mathit{\epsilon }}_{n+\mathrm{1}}^{\text{vp}}\right).$

(24)

In the adopted numerical solution procedure a viscoplastic tangent operator C^vp has been utilized [18]. This operator holds for arbitrary yield criteria and for different viscoplastic constitutive models. Thus an enhancement relative to other procedures is obtained; see, for example, [5]. For more details, see also De Angelis [33].

4. Computational Results

We consider the problem of an infinitely long internally pressurized thick-walled cylinder subject to an increasing internal pressure; see, for example, Ortiz et al. [34]. The boundary value problem is assumed to be in a state of plane-strain. The cylinder is characterized by an internal radius of 100 mm and an external radius of 200 mm. In the example 4-node bilinear isoparametric quadrilateral elements have been used and, for symmetry reasons, only 1/4 of the cylinder is considered. The present numerical algorithmic scheme has been implemented into the finite element analysis program (FEAP); see Zienkiewicz et al. [35] and Taylor [36]. The adopted finite element mesh consists of 221 nodes and 192 elements.

The adopted material properties are elastic modulus E = 2.1·10⁵ MPa, Poisson's ratio ν = 0.3, yield limit σ _yo = 240 MPa, and hardening moduli H_kin = H_iso = 0 MPa. In the analysis J₂ material behavior is assumed by a viscoplastic constitutive law of the Perzyna-type. Loading is performed by linearly increasing the internal pressure p up to the final maximum value p _o = 192 MPa which results, for the rate-independent perfectly plastic behavior, in a full plasticization of the cylinder.

For modeling the elasto-/viscoplastic material behavior and the related rate-dependent effects, a nondimensional loading program parameter τ is introduced as

$\tau =\frac{t_R}{{\sigma }_{yo}}\frac{\Delta p}{\Delta t},$

(25)

which takes into account the loading rate Δp/Δt and the intrinsic properties of the material by means of the relaxation time t _R = η/2G and the yield stress σ _yo .

During the loading process the plastic strains appear at first in correspondance with the internal radius of the cylinder. By increasing the internal pressure the plastic interface evolves from the cylindrical surface corresponding to the internal radius to the cylindrical surface corresponding to the external radius.

The pressure versus radial displacement curve is plotted in Figure 1 for a value of the loading program parameter τ = 0, which corresponds to an elastic perfectly plastic material behavior. In Figure 2 for a static loading program corresponding to τ = 0 the maximum equivalent plastic strain is plotted as a function of the radial displacement. In Figure 3 the maximum equivalent plastic strain is plotted as a function of the internal pressure for τ = 0.

Figure 1:

Pressure versus radial displacement for a loading program parameter τ = 0.

Figure 2:

Maximum equivalent plastic strain versus radial displacement for a loading program parameter τ = 0.

Figure 3:

Maximum equivalent plastic strain versus internal pressure for a loading program parameter τ = 0.

The effectiveness of the present numerical approach is clearly shown by comparing the above results, related to a loading rate parameter τ = 0, with the analytical solution of the thick-walled cylinder subject to increasing internal pressure which is provided by Prager and Hodge Jr. [37] for the elastic perfectly plastic behavior. The excellent correspondance between the herein reported numerical results and the analytical results present in the literature illustrates the robustness and soundness of the proposed computational approach. This correspondance between numerical results and analytical solution validates the proposed numerical procedure.

For an elasto-/viscoplastic material behavior the pressure versus radial displacement curves are plotted in Figure 4 for different values of the loading program parameter τ at constant material properties, that is, for different values of the prescribed loading rate.

Figure 4:

Pressure versus radial displacement for different values of the loading program parameter τ.

The rate-independent elastoplastic behavior is correctly recovered for τ = 0, with a static imposition of the load, whereas the elasto-/viscoplastic material behavior is described by the curves characterized by nonnull values of τ, corresponding to nonnull values of the prescribed loading rate.

For a higher prescribed loading rate it may be observed in Figure 4 that a higher pressure is required for the yielding of the material relative to the case of static loading; see, for example, Skrzypek and Hetnarski [26] and Lemaitre and Chaboche [28].

In Figure 5 the contour plot of the equivalent plastic strain is illustrated by showing the evolution of the plastic interface for an increasing value of the internal pressure p = 0.98p _o and a loading program parameter τ = 0.

Figure 5:

Contour plot of the equivalent plastic strain showing the evolution of the plastic interface for an increasing value of the internal pressure p = 0.98p _o and a loading program parameter τ = 0.

The solution strategy is well suited to be applied to different yield criteria and different constitutive models with the appropriate modifications, thus enhancing other procedures; see, for example, [5, 14].

5. Conclusions

In the present paper a procedure for the numerical analysis of elasto-/viscoplastic material behavior has been illustrated. An evaluation of the effects of the different loading programs on the nonlinear analysis of elasto-/viscoplastic materials is reported. An appropriate solution strategy is adopted which is well suited to be applied for general cases with the convenient modifications on the procedure.

In the present paper, differently from other previous analyses, for example, De Angelis et al. [38] and De Angelis and Cancellara [39], the loading is applied by increasing the pressure rather than by increasing the prescribed displacements. Accordingly, a suitable nondimensional loading program parameter has been introduced. The present paper represents a natural evolution of the analysis presented by De Angelis and Cancellara [40], where the effects of the loading procedures in the nonlinear inelastic behavior of solid materials have also been introduced.

An appropriate computational algorithmic scheme has been presented for the simulation of the elasto-/viscoplastic material behavior of solids. The proposed computational approach properly reduces to the inviscid limit for null viscosity parameter. A suitable solution algorithm and a consistent tangent operator have been applied which can be adopted for general computational procedures and different viscoplastic constitutive models by suitably specifying the flow function of the elasto-/viscoplastic constitutive model in use. The consequences of different loading programs on the nonlinear behavior of elasto-/viscoplastic materials have been analyzed in detail with the appropriate numerical examples.

Numerical computations and results have been finally reported so that the rate-dependency of the elasto-/viscoplastic constitutive model has been described and the influence of the loading rate on the nonlinear elasto-/viscoplastic material behavior has been illustrated with specific numerical examples.

Conflict of Interests

The author declares that there is no conflict of interests regarding the publication of this paper.