An Exact Axisymmetric Solution in Anisotropic Plasticity

A hollow cylinder of incompressible material obeying Hill’s orthotropic quadratic yield criterion and its associated flow rule is contracted on a rigid cylinder inserted in its hole. Friction occurs at the contact surface between the hollow and solid cylinders. An axisymmetric boundary value problem for the flow of the material is formulated and solved, and the solution is in closed form. A numerical technique is only necessary for evaluating ordinary integrals. The solution may exhibit singular behavior in the vicinity of the friction surface. The exact asymptotic representation of the solution shows that some strain rate components and the plastic work rate approach infinity in the friction surface’s vicinity. The effect of plastic anisotropy on the solution’s behavior is discussed.


Introduction
One can rarely find an analytical solution to an axisymmetric problem in plasticity. Such solutions exist for the one-dimensional classical problems, such as expanding tubes and spherical shells [1]. A solution for the flow of rigid perfectly plastic material through an infinite conical channel has been presented in [2]. This solution is valid for an arbitrary isotropic pressure-independent yield criterion. Several papers have generalized the solution [2]. In particular, a solution for the flow of rigid/linear-hardening materials has been given in [3]. This solution has been further generalized on an arbitrary nonlinear strainhardening law in [4]. Paper [5] has dealt with the radial flow of rigid/linear-hardening materials between two conical walls. The solution [2] has been used in [6,7] for constructing a solution for the flow of multilayered material. The same approach in conjunction with the solutions given in [4,5] has been adopted in [6] to include strain hardening into consideration. The radial flow of viscoplastic solids has been studied in [8], where solutions for other rate-dependent materials have been discussed. Solutions with non-zero circumferential velocity have been derived in [9,10]. In both papers, the material response is rigid perfectly plastic. Paper [11] has generalized the solution [2] on the double shearing model of ideal granular material. The double shearing model is described in this paper as well. It is worthy of note that of the solutions above, only the rigid perfectly plastic solutions are exact. The other solutions involve more approximations or have more limitations than the solution [2].
Another series of exact solutions has started from the solution given in [12]. Considered is a rigid/perfectly plastic hollow cylinder with a rigid solid cylinder inserted in its hole. The rigid/perfectly plastic cylinder is subject to uniform contraction over its outer surface. The constitutive equations comprise Tresca's yield criterion and its associated flow rule. In [13] the solution [12] has been extended to the rigid/perfectly plastic material model based on the yield criterion proposed in [14]. Papers [15] and [16] have generalized the solution [12] on a model of viscous material and the double shearing model, respectively. generalized the solution [12] on a model of viscous material and the double shearing model, respectively.
The importance of analytical solutions is at least twofold. Firstly, such solutions can be used as benchmark problems for verifying numerical codes, which a necessary step before using such codes [17,18]. Secondly, analytical solutions are useful for understanding the qualitative behavior of general solutions. Of particular interest for rigid/plastic models is the maximum friction surfaces where the velocity field may be singular [19]. This singularity causes difficulty with numerical solutions. In particular, finite element solutions based on typical shape functions do not converge [20,21].
One of the present paper's objectives is to extend the solution [12] to a model of anisotropic plasticity. It is known that even mild plastic anisotropy may significantly affect some features of isotropic solutions [22]. An example of such an effect has been presented, for example, in [23]. It is assumed that the material obeys Hill's quadratic yield criterion and its associate flow rule [1]. This model is widely used until now [24][25][26][27]. It is known that plane strain solutions for this model may be singular [28]. This singularity is associated with envelopes of characteristics. The equations for axisymmetric flow are not hyperbolic. However, isolated characteristic surfaces in ideal plasticity may exist even if the equations are not hyperbolic [29]. Another objective of the present paper is to demonstrate it for the model adopted using a particular solution.

Statement of the Problem
A rigid solid cylinder is inserted into the hole of a rigid/plastic hollow cylinder ( Figure 1). The radius of the rigid cylinder and the radius of the hollow cylinder's hole is denoted by 0 a , and the outer radius of the hollow cylinder by 0 b . The length of both the cylinders is 2L. The hollow cylinder is subject to uniform contraction over its outer surface. Cylindrical coordinates ( ) ,, rz  are taken, with the z-axis coinciding with the axis of symmetry of each cylinder. The plane 0 z = is a plane of symmetry of the boundary value problem. The solution is independent of  . The boundary conditions are the same as those in [12]. In particular, the boundary conditions at 0 ra = and 0 rb = are satisfied exactly, except for the condition on the radial stress. The exact boundary conditions at 0 z = and zL = are replaced with integral boundary conditions. Since 0 z = is the plane of symmetry for the flow, it is sufficient to find the solution in the region 0 zL  .  The non-zero components of the stress tensor in the cylindrical coordinates are denoted by σ rr , σ θθ , σ zz , and σ rz ; and the non-zero components of the velocity vector by u r and u z . Since the inner cylinder is rigid, for r = a 0 . The other velocity boundary condition that is exactly satisfied is for r = b 0 . Here U > 0. The exact velocity boundary condition at the plane of symmetry is replaced with for z = 0. The stress boundary conditions that are exactly satisfied are for r = b 0 and Prandtl's friction law at r = a 0 . This law depends on the constitutive equations and will be formulated below. The integral stress boundary condition proposed in [12] reads At the end of the hollow cylinder, an axial stress is required to deform it according to the boundary conditions above. Once the boundary value problem has been solved, the average value of this stress is calculated as This stress is involved in the method for predicting the brittle fracture of fibers in composites proposed in [12].
The material of the hollow cylinder is plastically orthotropic. The principal axes of anisotropy coincide with coordinate curves of the cylindrical coordinate system. The elastic portion of strain is neglected. The constitutive equations comprise the yield criterion proposed in [1] and its associated flow rule. In the case under consideration, the yield criterion reads Here where R is the tensile yield stress in the radial direction, Θ is the tensile yield stress in the circumferential direction, Z is the tensile yield stress in the axial direction, and 1/ √ 2M is the shear yield stress with respect to the r-, z-axes. All of them are supposed to be constant. The flow rule is Here λ is a non-negative multiplier. Prandtl's friction law for isotropic materials postulates that the friction stress τ f is a fraction of the shear yield stress. Therefore, it is reasonable in the case under consideration to represent this law as In what follows, it is convenient to use the nondimensional quantities:

Analytic Solution
The solution provided in [12] suggests that one looks for velocity solutions of the form Here η(ρ) and µ(ρ) are arbitrary functions only of ρ, and A is constant. Equation (8) results in the equation of incompressibility ξ ρρ + ξ θθ + ξ zz = 0. In terms of the velocity components, this equation becomes ∂u r /∂r + u r /r + ∂u z /∂z = 0. Or, using Equation (10), Substituting Equation (11) into Equation (12) one gets This equation can be immediately integrated to give Here A 0 is constant. Equations (1), (2), (10), (11), and (14) combine to give Substituting Equation (15) into Equation (14) and the resulting equation into (11) leads to Then, using Equation (10), one can find the non-zero components of the strain rate tensor referred to the cylindrical coordinates as b 0 In a generic meridian plane, the slope ϕ of the principal axes of strain rate with respect to the r-axis is given by Equations (17) and (18) combine to give Symmetry 2021, 13, 825

of 15
The slope ϕ can also be expressed in terms of the stress components using Equations (8) and (18). As a result, It is convenient to introduce the new stress variables p 1 , p 2 , and p 3 as It is evident that Equations (20) and (21) combine to give It follows from Equations (8) and (21) (10) and (16), one transforms the latter equation to One can solve Equations (22) and (24) for p 3 and p 2 to get In Equations (24) and (25), t and w are functions of ρ defined as Equations (23) and (25) combine to give Substituting Equation (21) into Equation (7) and using Equations (25) and (27), one arrives at the following equation for p 1 Assuming that the components of the deviatoric stress tensor are independent of z, the equilibrium equations reduce to ∂σ ∂ρ Here σ is the hydrostatic stress, τ rr = σ rr − σ, and τ θθ = σ θθ − σ. The third normal component of the deviatoric stress tensor is τ zz = σ zz − σ. The Equations in (29) where B is constant and σ 0 (ρ) is a function of ρ. Substituting Equation (30) into Equation (29), one gets dσ 0 dρ + dτ rr dρ + τ rr − τ θθ ρ = 0 and dσ rz dρ where C is constant. The boundary conditions in Equations (4) and (9) serve for determining B and C. As a result, using Equation (10), Then, Equation (32) becomes Equations (27) and (34) combine to give Eliminating p 1 between Equations (28) and (35), one arrives at the following equation for tan 2ϕ Substituting tan 2ϕ back into Equation (35) leads to It is seen from Equation (26) that t > 0, t → ∞ as ρ → a , and dt/dρ < 0 for all ρ. Moreover, w → 1 as t → ∞ (or ρ → a ) and dw/dt < 0 for all t. Summarizing the relations above, one can conclude that the direction of the friction stress ( Figure 1) demands that dµ/dρ > 0. Then, Equations (19) and (36) combine to give w ≥ 1 and t > 0. (38) Here, Equation (38) has been taken into account. It follows from Equation (39) that Here, µ 0 is constant. Its value is found applying Equation (3). Using Equations (10), (11) and (40), one gets The integrals in Equations (40) and (41) should be evaluated numerically. Equation (21) is equivalent to These equations can be solved for the components of the deviatoric stress tensor. As a result, Substituting Equation (43) into the first equation in Equation (31), one finds Equations (19), (35), and the inequality dµ/dρ > 0 demands that p 1 < 0. Then, the lower sign should be chosen in Equation (37). Using this equation, one can rewrite Equation (44) as Integrating and using Equation (43) gives Equation (5) serves for determining D. The radial stress is found from Equations (30), (33) and (46) as Substituting Equation (47) into Equation (5) and using Equation (10) gives This equation completes the solution. Summarizing the solution above, one finds the velocity field from Equations (16), (40) and (41). Equations (34), (47), and (48) supply the shear and radial stresses. The other normal stresses follow from Equations (30) and (33), and the equations σ θθ = σ + τ θθ and σ zz = σ + τ zz . The deviatoric stress components involved in these equations are determined from Equations (25), (37) and (43). Having found the axial stress, one calculates the value of q from Equation (6).
The plastic work rate W can be calculated as [1]: where σ and ξ are the equivalent stress and strain rate, respectively, which are given by

Asymptotic Analysis
Many rigid/plastic solutions are singular in the vicinity of maximum friction surfaces [19,28]. In the case under consideration, the maximum friction surface is ρ = a if m = 1. Equation (39) at m = 1 becomes It is seen from Equation (26) that t → ∞ as ρ → a . However, Therefore, the function ρ 2 1 − a 2 2 − a 2 1 − ρ 2 2 controls the behavior of the righthand side of Equation (52) as ρ → a . Since as ρ → a , it is evident from Equation (52) that as ρ → a . Then, it follows from Equation (11) that the shear strain rate approaches infinity in the vicinity of the maximum friction surface. In particular, as ρ → a . The normal strain rates referred to the cylindrical coordinate system are not singular. Using Equations (53) and (54), one can rewrite Equation (45) as as ρ → a . Integrating as ρ → a . It follows from Equations (25) and (37) in which the lower sign is chosen, (53), and (54) that as ρ → a . If F = H, Equations (43) and (59) combine to give as ρ → a . Comparing Equations (58) and (60), one concludes that as ρ → a . At any ζ, it follows from Equations (30) and (62) that as ρ → a . Equations (58), (60), and (63) show that the absolute values of the derivatives ∂σ θθ /∂ρ and ∂σ zz /∂ρ approach infinity as ρ → a but the absolute value of the derivative ∂σ rr /∂ρ does not. The asymptotic Expansions (56), (60), (61) and (63) show that traditional finite elements are not capable of solving this problem. It has been already demonstrated in [21] for an isotropic model that calculation does not converge. These asymptotic expansions can be used in conjunction with the extended finite element method [38].
Equations (51) and (56) show that the equivalent strain rate approaches infinity near the bi-material interface. This strain rate controls the evolution of many material properties. Therefore, Equations (51) and (56) predict a high gradient of such properties near the bimaterial interface. This phenomenon has been known for a long time (see, for example, [39]). In particular, the generation of microstructure in the vicinity of bi-material interfaces has been studied in [40,41]. Nevertheless, the asymptotic expansions found cannot be directly used in conjunction with conventional evolution equations for material properties because ξ → ∞ as ρ → a . An approach to overcome this difficulty for isotropic materials has been proposed in [42]. Equation (56) suggests that a similar approach can be developed for anisotropic materials.

Numerical Example
The numerical example provided in this section evaluates the effect of the hollow cylinder's material anisotropy on the stress and strain fields. It was assumed that the hollow cylinder is made from aluminum alloy AA5086. Due to various parameters of the manufacturing process, the same material may exhibit different anisotropic properties characterized by the symmetry of mechanical properties with respect to the principal axes of anisotropy [43][44][45]. The following types of anisotropy were considered: isotropic, transversely isotropic (mechanical properties are identical in any direction in the zθ-coordinate surface but differ from the properties in the other coordinate surfaces) and orthotropic materials. The AA5086 alloy's mechanical properties involved in yield criterion (7) are available in the literature [46]. The mechanical properties responsible for plastic anisotropy are summarized in Table 1. The solution provided in Section 3 has been applied to calculate the distribution of the velocities, strain rates, stresses, and plastic work rate at m = 1. It has also been assumed that a 2 = 1/2 and L/b 0 = 10. The first equation in Equation (16) supplies the radial distribution of the radial velocity at any value of a. This simple formula does not require a graphical illustration. No material properties affect the radial velocity. The second equation in Equation (16) shows that the function µ(ρ) completely controls the influence of plastic anisotropy on the axial velocity. Figure 2 illustrates this function. Using this Figure, one can visualize the radial distribution of the axial velocity with ease by adding the first term on the second equation's right-hand side in Equation (16). This term is constant at any value of ζ. It is seen from Figure 2 that the influence of plastic anisotropy on the axial velocity is not significant. However, one can see an interesting qualitative feature of the solution that the magnitude of this velocity in the vicinity of the friction surface is largest for transversely isotropic material, even though this material model is intermediate with respect to the other two models. Another qualitative feature of the solution is that the magnitude of the axial velocity in the cylinder of transversely isotropic material is the largest in the vicinity of the friction surface but is the smallest at the outside radius, as compared to the other two models. The curves in Figure 2 are in qualitative agreement with Equation (55). The solution provided in Section 3 has been applied to calculate the distribution of the velocities, strain rates, stresses, and plastic work rate at 1 m = . It has also been assumed that 2 12 a = and 0 10 Lb = . The first equation in Equation (16) supplies the radial distribution of the radial velocity at any value of a. This simple formula does not require a graphical illustration. No material properties affect the radial velocity. The second equation in Equation (16) shows that the function ( )  completely controls the influence of plastic anisotropy on the axial velocity. Figure 2 illustrates this function. Using this Figure, one can visualize the radial distribution of the axial velocity with ease by adding the first term on the second equation's right-hand side in Equation (16). This term is constant at any value of  . It is seen from Figure 2 that the influence of plastic anisotropy on the axial velocity is not significant. However, one can see an interesting qualitative feature of the solution that the magnitude of this velocity in the vicinity of the friction surface is largest for transversely isotropic material, even though this material model is intermediate with respect to the other two models. Another qualitative feature of the solution is that the magnitude of the axial velocity in the cylinder of transversely isotropic material is the largest in the vicinity of the friction surface but is the smallest at the outside radius, as compared to the other two models. The curves in Figure 2 are in qualitative agreement with Equation (55). The strain rate tensor's normal components are given by the simple formulae presented in Equation (17). Therefore, no graphical illustration of these components is required. The last equation in Equation (17) shows that the function ( )  completely controls the influence of plastic anisotropy on the shear strain rate component. The variation of this component with the dimensionless radius is depicted in Figure 3. The effect of plastic anisotropy is not significant. As in the case of the axial velocity, the curve corresponding to transversely isotropic material does not lie between the curves corresponding to the isotropic and orthotropic materials. The magnitude of the shear strain The strain rate tensor's normal components are given by the simple formulae presented in Equation (17). Therefore, no graphical illustration of these components is required. The last equation in Equation (17) shows that the function µ(ρ) completely controls the influence of plastic anisotropy on the shear strain rate component. The variation of this component with the dimensionless radius is depicted in Figure 3. The effect of plastic anisotropy is not significant. As in the case of the axial velocity, the curve corresponding to transversely isotropic material does not lie between the curves corresponding to the isotropic and orthotropic materials. The magnitude of the shear strain rate approaches infinity in the vicinity of the friction surface in accordance with Equation (56).   The simple formula in Equation (34) gives the solution for the shear stress. The effect of plastic anisotropy vanishes if this stress is normalized by √ M. Therefore, no graphical illustration of this stress is required. The most significant influence of plastic anisotropy is observed on τ rr , τ θθ , and σ (Figure 4). The change of the type of anisotropy affects these distributions quantitively. In particular, the isotropic model results in the largest absolute value of τ rr , the smallest absolute value of τ θθ , and the intermediate value of σ. The transversely isotropic model results in the smallest absolute value of τ rr , the largest absolute value of τ θθ , and the largest absolute value of σ. The orthotropic model results in the intermediate absolute values of both τ rr and τ θθ , and the smallest absolute value of σ. The effect of the type of plastic anisotropy on the deviatoric stress τ zz is negligible. The difference in the value of σ between the orthotropic and transversely isotropic models is about 10%. This result is important for predicting the brittle fracture of fibers in composites using the method proposed in [12]. It is seen from Figure 4 that the tangents to all the curves corresponding τ rr , τ zz , and σ tend to a vertical line in the vicinity of the inner radius. This behavior of the curves is in accordance with Equations (60) and (63). The curves corresponding to τ θθ reveal the same feature for the orthotropic and transversely isotropic models but not for the isotropic model. The latter agrees with Equation (61) since F = H = G in the case.
The simple formula in Equation (34) gives the solution for the shear stress. The effect of plastic anisotropy vanishes if this stress is normalized by M . Therefore, no graphical illustration of this stress is required. The most significant influence of plastic anisotropy is observed on rr  ,   , and  ( Figure 4). The change of the type of anisotropy affects these distributions quantitively. In particular, the isotropic model results in the largest absolute value of  is negligible. The difference in the value of  between the orthotropic and transversely isotropic models is about 10%. This result is important for predicting the brittle fracture of fibers in composites using the method proposed in [12]. It is seen from Figure  4 that the tangents to all the curves corresponding rr  , zz  , and  tend to a vertical line in the vicinity of the inner radius. This behavior of the curves is in accordance with Equations (60) and (63). The curves corresponding to   reveal the same feature for the orthotropic and transversely isotropic models but not for the isotropic model. The latter agrees with Equation (61) since F H G == in the case. The radial distribution of the plastic work rate calculated using Equations (49)-(51) is depicted in Figure 5. Since the stresses and normal strain rates are bounded, it is seen from these equations and Equation (56) that This asymptotic representation of the solution is visible in Figure 5.  The radial distribution of the plastic work rate calculated using Equations (49)-(51) is depicted in Figure 5. Since the stresses and normal strain rates are bounded, it is seen from these equations and Equation (56)

Prediction of the Brittle Fracture of Fibers in Composites
The solution above can be used for predicting the failure of composite sheets made of ductile matrix and brittle fibers subject to tension in the direction of fibers using the approach suggested in [12]. According to this approach, the failure occurs by plastic flow

Prediction of the Brittle Fracture of Fibers in Composites
The solution above can be used for predicting the failure of composite sheets made of ductile matrix and brittle fibers subject to tension in the direction of fibers using the approach suggested in [12]. According to this approach, the failure occurs by plastic flow of the matrix if L < L c and by brittle fracture of the fiber if L > L c . The critical length L c should be found from the solution. The criterion for finding this length is that both failure mechanisms occur simultaneously. Equation (6) at L = L c supplies the tensile load q c at which the fiber breaks. The parameters from Section 5 have been used in the calculation. In addition, the mean tensile stress at which the fiber breaks, T, is required. Its magnitude has been varied in the range typical for the fibrous composites [47].
The change of the type of anisotropy affects the ultimate tensile strength quantitively ( Figure 6). In particular, the orthotropic model results in the largest values of q c on the range of T, the smallest values of q c correspond to the transversely isotropic model; the isotropic model results in the intermediate values. The lower the difference between the ultimate stress of the fiber and the yield stress of the matrix, the greater the effect of the plastic properties anisotropy on the brittle fracture of the fiber.

Prediction of the Brittle Fracture of Fibers in Composites
The solution above can be used for predicting the failure of composite sheets made of ductile matrix and brittle fibers subject to tension in the direction of fibers using the approach suggested in [12]. According to this approach, the failure occurs by plastic flow of the matrix if In addition, the mean tensile stress at which the fiber breaks, T, is required. Its magnitude has been varied in the range typical for the fibrous composites [47].
The change of the type of anisotropy affects the ultimate tensile strength quantitively ( Figure 6). In particular, the orthotropic model results in the largest values of

Conclusions
An exact axisymmetric solution has been found for the stress and velocity fields in a hollow rigid/plastic anisotropic cylinder contracted a solid rigid cylinder. Three types of Figure 6. Dependence of the tensile load at which the fiber breaks on its ultimate stress: --isotropic material; ----orthotropic material; ···· -transversely isotropic material.

Conclusions
An exact axisymmetric solution has been found for the stress and velocity fields in a hollow rigid/plastic anisotropic cylinder contracted a solid rigid cylinder. Three types of anisotropy have been adopted for obtaining quantitative results: orthotropic material, transversely isotropic material, and isotropic material. The following conclusions have been reached from this theoretical analysis: 1. The type of anisotropy has a negligible effect on kinematics variables and the plastic work rate. This feature is attributed to the boundary conditions that specify the radial velocity at both the inner and outer radii of the rigid/plastic cylinder and the equation of incompressibility valid for all the models adopted.
2. The type of anisotropy has a significant effect on the radial and circumferential deviatoric stress and hydrostatic stress. The type of anisotropy has a negligible effect on the axial deviatoric stress. A consequence of such a solution's behavior is that the type of anisotropy has a significant effect on the axial stress involved in Equation (6).
3. If the friction factor is equal to unity, then the solution becomes singular in the vicinity of the inner radius of the rigid/plastic cylinder. In particular, the shear strain rate and plastic work rate approach infinity near this radius according to an inverse square root rule. This behavior may cause difficulties with solving other problems using numerical methods.
4. The solution found can be used in conjunction with the approach proposed in [12] for predicting the brittle fracture of fibers in composites.