Evolution of internal variables in an expanding hollow cylinder at large plastic strains

An efficient method for calculating the evolution of internal variables in an expanding hollow cylinder of rigid/plastic material is proposed. The conventional constitutive equations for rigid plastic, hardening material are supplemented with quite an arbitrary set of evolution laws for internal variables assuming that the material is incompressible. No restriction is imposed on the hardening law. The problem is solved in Lagrangian coordinates. This significantly facilitates a numerical treatment of the problem. In particular, the initial/boundary value problem is reduced to a system of equations in characteristic coordinates. A finite difference scheme is used for solving these equations. An illustrative example is presented assuming that the internal variables are the equivalent plastic strain and a damage parameter.

model capable of such predictions is based on the conventional system of equations of plasticity theory (i.e. a yield criterion and its associated flow rule) and evolution equations for internal variables. These equations should be solved simultaneously. Commonly used internal variables are the equivalent strain, various damage parameters (Lemaitre 1985;Chandrakanth and Pandey 1993; Besson 2010 among many others), dislocation density (Ganapathysubramanian and Zabaras 2004;Lin and Dean 2005;He et al. 2012;Hore et al. 2015) and many other parameters that characterize the microstructure of material.
It is shown in the present paper that for a generic set of evolution equations for internal variables the initial/boundary value problem for an expanding hollow cylinder is reduced to a hyperbolic system of equations. The characteristic curves of this system are coordinate lines of a Lagrangian coordinate system. Therefore, a very high accuracy of numerical solutions can be easily achieved using a finite difference method and these solutions can serve as simple benchmark tests for numerical packages that deal with the prediction of the evolution of internal variables in metal forming processes. A necessity of such tests has been pointed out by Roberts et al. (1992). Note that even in the case of linear elasticity numerical results often depend on a particular method adopted to solve the problem (Helsing and Jonsson 2002). Therefore, reliable benchmark tests are important for large strain plasticity models that include internal variables.
The method proposed is an extension of the method developed in Alexandrov and Jeng (2011). In this paper, the evolution of damage in an expanding/contracting hollow sphere has been investigated.

Statement of the problem
Consider a long thick-walled tube of initial outer radius b 0 and initial inner radius a 0 . The tube is subject of uniform pressure p 0 applied over its inner radius. It is assumed that the magnitude of p 0 is high enough to initiate plastic yielding everywhere in the tube. Since the model adopted is rigid plastic, the magnitude of p 0 is determined from the solution. The current inner and outer radii of the tube will be denoted by a c and b c , respectively. It is convenient to introduce a cylindrical coordinate system (r, θ , z) with its z-axis coinciding with the axis of symmetry of the tube. The state of strain is two-dimensional in this coordinate system (ε z = 0). Here ε z is the axial strain. The initial/boundary value problem is axisymmetric and its solution is independent of θ. The circumferential velocity vanishes everywhere. The normal stresses in the cylindrical coordinate system are the principal stresses. The stress boundary condition is for r = b c . Here σ r is the radial stress (σ θ will stand for the circumferential stress). Models of rate-independent plasticity will be considered in the present paper. Therefore, it is possible, with no loss of generality, to arbitrarily prescribe the radial velocity at r = b c . It is convenient to put for r = b c . Here u is the radial velocity and u 0 is constant.
According to a widely used model of hardening plasticity the yield stress depends on the equivalent strain, ε eq , and the latter is defined by the following equation Here t is the time, d dt denotes the convected derivative and ξ eq is the equivalent strain rate. In the case under consideration, the latter is expressed in terms of the physical components of the strain rate tensor in the cylindrical coordinate system as where ξ r and ξ θ are the radial and circumferential strain rates, respectively. It has been taken into account here that the axial strain rate vanishes everywhere since ε z = 0. The constitutive equations include the yield criterion, its associated flow rule and internal variable evolution equations. The von Mises yield condition modified to account for the effect of internal variables on the yield stress is written as where σ eq is the equivalent stress, σ 0 is a reference stress, α p are internal variables other than ε eq (1 ≤ p ≤ k), � ε eq is an arbitrary function of its argument satisfying the conditions �(0) = 1 and d� dε eq ≥ 0 for all ε eq . In the case under consideration, the equivalent stress is defined by It has been taken into account here that σ θ > σ r in the case of tube expansion. The associated flow rule reduces to Here is a non-negative multiplier. The magnitude of is not important for applications. It is evident that ξ θ > 0 in the boundary value problem under consideration. Therefore, the inequality > 0 is satisfied and Eq. (7)  where σ is the hydrostatic stress and F is quite an arbitrary function of its arguments.
(3) dε eq dt = ξ eq . (4) In the case under consideration the hydrostatic stress is defined as The initial conditions to Eqs. (3) and (10)  It is convenient to introduce the following dimensionless quantities

Analytic treatment
By definition, u 0 = db c dt and u = dr dt. Therefore, using Eq. (15) Since ξ r = ∂u ∂r and ξ θ = u/ r, Eq. (8) becomes ∂u ∂r + u/ r = 0. Integrating this equation and using the boundary condition (2) lead to u = u 0 b c r −1 or using Eq. (15) to Equations (16) and (17) Integrating this equation and using the initial condition (12) yield Equations (5) and (6) combine to give Using Eqs. (15), (18), (19) and (22), Eq. (14) can be transformed to The equivalent strain in this equation should be eliminated by means of Eq. (21). Therefore, the right hand side of Eq. (23) is a function of b, R and α p (1 ≤ p ≤ k). In the Lagrangian coordinates Eq. (10) is Here Eqs. (9), (15), (17) and (18) have been used. The first argument of the function F is determined from Eqs. (5), (6) and (11) as The equivalent strain in Eqs. (24) and (25) should be eliminated by means of Eq. (21). Therefore, the right hand side of Eq. (24) is a function of b, R, σ r and α p (1 ≤ p ≤ k). Thus Eqs. (23) and (24) constitute the system of equations for σ r and α p in the Lagrangian coordinates. It is evident that the characteristic curves of this system are R = constant and b = constant. The boundary condition (1) becomes for R = 1. Using this condition Eq. (24) can be integrated along the characteristic curve R = 1. In particular, substituting Eq. (26) into Eqs. (21) and (25) and, then, the resulting expressions for σ r and ε eq into Eq. (24) yield The initial conditions to these equations follow from Eq. (13) in the form for b = 1. It is evident that Eq. (27) can be solved for any given function F. Analogously, Eq. (23) can be integrated along the characteristic curve b = 1. In particular, using Eqs. (12) and (13)

Integrating this equation and using the boundary condition (26) yield
The solution of Eqs. (27) and (29) facilitate a numerical treatment of Eqs. (23) and (24).

Numerical treatment
A finite difference method is adopted to solve Eqs. (23) and (24). N + 1 nodes are chosen on the spatial coordinate and T + 1 nodes on the time-like coordinate. Therefore, the mesh used is The value of β has been introduced in Eq. (15) (24) the derivatives ∂σ r ∂R and ∂α p ∂b can be found at the point j = N and i = 2, ∂σ r ∂R | N ,2 and ∂α p ∂b | N ,2 , 1 ≤ p ≤ k. Then, the derivatives within the intervals 1 − R ≤ R ≤ 1 and 1 ≤ b ≤ 1 + b are approximated by Replacing the derivatives in Eq. (31) by the values given in Eq. (32) a more accurate approximation of σ r and α p at the point j = N and i = 2 is determined. It is obvious that this procedure can be extended to N − 1 ≥ j ≥ 1. As a result, the distribution of σ r and α p along the line i = 2 is obtained. The line i = 3 can be treated in a similar manner since the solution of Eq. (27) is available for any b. This procedure can be extended to any number i.
where M is the Ramberg-Osgood hardening exponent, and q and A are material constants. Then, Eq. (27) becomes Integrating this equation and using the initial condition (28)

Conclusions
The initial/boundary value problem for calculating the evolution of internal variables in an expanding hollow tube at large strains has been reduced to simple equations in characteristic coordinates. An advantage of the approach proposed is that an arbitrary hardening law and an arbitrary set of internal variable evolution equations are included in the formulation. The equilibrium equation has been integrated analytically at b = 1 [see Eq. (29)]. The internal variable evolution equations can be integrated analytically or, in the most general case, reduced to numerical integration at R = 1 independently of the general solution of the initial/boundary value problem [see Eq. (27)]. These solutions have been used in the numerical scheme proposed. They can also be used in conjunction with any other numerical scheme to increase the accuracy of calculation and to reveal possible errors in numerical solutions. The illustrative example deals with the evolution of damage. It is seen from Figs. 1, 2 and 3 that the initiation of ductile fracture occurs at the inner radius of the tube in all the cases considered. It is worth noting here that the fracture initiation in an expanding thick ring may occur at mid-annulus (Tomkins and Atkins 1981). It is evident that the damaged model adopted in the present paper is not appropriate in such cases. However, the general solution found and the general numerical scheme developed are independent of the specific form of the functions and F. The present general solution can be used to test various function involved in Eqs. (5) and (10) to find functions and F that are in agreement with the experimental result presented by Tomkins and Atkins (1981).
The accuracy of the numerical solution is rather high. Therefore, the final result is useful for verifying more general numerical codes, which can be used to solve arbitrary initial/boundary value problems.