Finite Plane Strain Bending under Tension of Isotropic and Kinematic Hardening Sheets

The present paper provides a semianalytic solution for finite plane strain bending under tension of an incompressible elastic/plastic sheet using a material model that combines isotropic and kinematic hardening. A numerical treatment is only necessary to solve transcendental equations and evaluate ordinary integrals. An arbitrary function of the equivalent plastic strain controls isotropic hardening, and Prager’s law describes kinematic hardening. In general, the sheet consists of one elastic and two plastic regions. The solution is valid if the size of each plastic region increases. Parameters involved in the constitutive equations determine which of the plastic regions reaches its maximum size. The thickness of the elastic region is quite narrow when the present solution breaks down. Elastic unloading is also considered. A numerical example illustrates the general solution assuming that the tensile force is given, including pure bending as a particular case. This numerical solution demonstrates a significant effect of the parameter involved in Prager’s law on the bending moment and the distribution of stresses at loading, but a small effect on the distribution of residual stresses after unloading. This parameter also affects the range of validity of the solution that predicts purely elastic unloading.


Introduction
Finite plane strain bending with and without tension is one of the classical problems of plasticity. All sheet-metal processes incorporate some bending, which demonstrates the practical importance of this problem. A brief overview of such processes was presented in [1]. Moreover, the process of bending is often used as a test for determining material properties, which requires a theoretical description of the process. In particular, cyclic bending tests were employed in [2,3]. A four-point bending test was carried out in [4] for investigating the large-strain elastic constitutive behavior of thin unidirectional composites. Anisotropic hardening of metal sheets was studied in [5]. An in-plane bending test was used in [6] for determining the uniaxial large strain work-hardening behavior of highstrength steel. A pure bending test was employed in [7] to study of the strength differential effect of DP98 steel sheets.
The first rigorous rigid plastic solution was found in [8] for the model of isotropic perfectly plastic material. The work in [9] extended this solution for bending with no tension to rigid plastic, work-hardening material. However, this solution assumes an artificial strain distribution between the current neutral line and the original center fiber. The study in [10] further extended the solution for pure plane strain bending given in [8] study in [10] further extended the solution for pure plane strain bending given in [8] to rigid plastic-laminated sheet metals. Anisotropic material properties were taken into account in [11]. Deformation theories of plasticity were employed to analyze plane strain bending at large strains in [12,13]. The work in [14] dealt with finite plastic bending of a compressible elastoplastic strip with combined isotropic and kinematic hardening at a given stretch normal to the bending plate. The flow theory of plasticity was adopted. Since the elastic portion of the strain tensor was not neglected in [14], the solution was much more complicated than those in [8][9][10][11]. It was shown in [15] that the assumption of incompressibility retained the simplicity of the rigid plastic solutions above for the corresponding elastic/plastic solutions. In particular, the Kirchhoff stress, which is involved in the equation that connects the elastic stretching and an objective stress rate, is equal to the Cauchy stress, and all objective stress rates reduce to the convected derivative. The authors of [15] adopted an inverse method that assumed transformation equations that connect Eulerian and Lagrangian coordinates, and showed that these equations described the geometry of plane strain bending at large strains. These equations can be used in conjunction with many constitutive equations [16,17]. The basic geometric assumptions used in [15] were similar to or even less restrictive than those accepted in papers published in journals devoted to metal forming [18,19].
The present paper extends the solution given in [17] to finite plane strain bending under tension of isotropic and kinematic hardening sheets. Prager's law [20] is used to describe kinematic hardening. Isotropic hardening is controlled by an arbitrary function of the equivalent plastic strain. Elastic unloading is also included in the solution.

Constitutive Equations
A wide elastic/plastic sheet is subject to bending under tension by the bending moment M and the tensile force F under plane strain conditions. The initial cross-section of the sheet is rectangular (A1B1C1D1 in Figure 1a). The initial thickness and width of the sheet are H and 2L, respectively. The Eulerian Cartesian coordinate system ( ) , x y is chosen such that its x-axis coincides with the axis of symmetry of the cross-section throughout the process of deformation. The origin of the coordinate system is situated at point O, which is the intersection of the axis of symmetry and surface AB at any deformation stage ( Figure 1b). The surface AB coincides with the plane A1B1 at the initial instant.  The commonly used Eulerian rate-type theory of elastoplasticity is adopted. This theory is presented, for example, in [21]. The Cauchy stress and Hencky strain are adopted throughout the paper. The material is assumed to be incompressible. Therefore, Poisson's ratio is 1/2. Let τ 1 and τ 2 be the principal deviatoric stresses in the planes of flow. It will be seen later that all objective rates reduce to the convected derivative. Then, the elastic portions of the principal strain rates, ξ e 1 and ξ e 2 , are given by: where the superimposed dot denotes the convected derivative and G is the shear modulus of elasticity. The plane strain criterion is assumed to be: where α 1 and α 2 are the principal back stresses, ε p eq is the equivalent plastic strain, σ 0 is the yield stress in tension at ε p eq = 0, Φ ε p eq is an arbitrary function of its argument satisfying the conditions Φ(0) = 1 and dΦ ε p eq /dε p eq ≥ 0 for all ε p eq . The following equation defines the equivalent strain: where ξ p 1 and ξ p 2 are the plastic portions of the principal strain rates. The overdot in (3) denotes the convected derivative. In what follows, Prager's law for the back stresses is accepted [20]. Then: where K is a constitutive parameter. The plastic flow rule associated with the yield criterion (2) leads to: ξ The total strain rate tensor is the sum of the elastic and plastic strain rate tensors. Then: Equations (1), (5) and (6) combine to give the equation of incompressibility:

Strain Analysis
For a class of constitutive equations, strain analysis can be carried out independently of the stress equations. The compatibility of the strain distribution found with the stress equations should be verified a posteriori. This approach for pure plane strain bending has been proposed in [15] and extended to plane strain bending under tension in [17]. These papers have not dealt with kinematic hardening. However, one can try to combine the strain analysis completed there and kinematic hardening laws. For the reason of readability, the main results of [15,17] are presented in this section. These results are also necessary for finding the stress solution in the next section.
In addition to the Cartesian coordinate system (x, y), one can introduce the Lagrangian coordinate system (ζ, η) such that: x = Hζ and y = Hη at the initial instant. Then, it is seen from Figure 1 that ζ = 0 on AB, ζ = −1 on CD, η = −L/H on AD, and η = L/H on CB throughout the process of deformation. The mapping between the (x, y) and (ζ, η) coordinate systems proposed in [15] is: where a is a time-like variable such that a = 0 at the initial instant and s is a function of a. This function should be found from the stress solution. The condition (8) is satisfied if: at a = 0. It is possible to verify by inspection that the mapping given by (9) and (10) satisfy Equation (7). Moreover, AB and CD are circular arcs, CB and AD are straight, and the coordinate curves of the (ζ, η)− coordinate system are principal strain rate and principal strain trajectories. Then, as noted in the previous section, all objective rates reduce to the convected derivative. The normal strains referred to (ζ, η)− coordinate system are determined from (9) as: These strains are the principal strains. The equations in (4), (5) and (6) can be immediately integrated to get: and: Equations (3) and (5) combine to give As before, one can integrate these equations to get: The inequalities should be satisfied throughout the process of deformation. Using (9), one can readily find that: where r AB is the radius of arc AB, r CD is the radius of arc CD, θ 0 is the orientation of line CB relative to the x-axis, and h is the current thickness of the sheet (Figure 1). The methodology of combining the strain analysis above and constitutive equations is provided in Appendix A.

Stress Analysis
The constitutive equations formulated in Section 2 result in the coaxiality of stress and strain rate tensors. Therefore, the coordinate curves of the Lagrangian coordinate system coincide with principal stress trajectories. In particular, the surfaces AB, CD, CB and AD are free of shear stresses. Moreover, the normal stresses referred to the (ζ, η)− coordinate system, σ ζ and σ η , are the principal stresses. Using (12) and (13), one can rewrite the yield criterion in Equation (2) as: It is assumed that the principal stresses are independent of η. Then, the only stress equilibrium equation that is not identically satisfied in the Lagrangian coordinates reads: Equilibrium demands that some pressure P applies over the surface CD if F = 0 ( Figure 1b). It is assumed that this pressure is uniformly distributed. Then, F = r CD P. Using (16), one determines: where p = P/σ 0 and f = F/(σ 0 H).
The dimensionless bending moment, m, is determined as [17]: where H/h should be eliminated by means of (16).

Purely Elastic Solution (Stage 1 of the Process)
The purely elastic solution is independent of the yield criterion. Therefore, the solution given in [17] is valid. In particular, the general solution in any elastic region is: where α = σ 0 /G and D 0 is constant. In deriving this result, it has been taken into account that τ 1 = τ ζ , τ 2 = τ η , τ ζ = σ ζ − σ, τ η = σ η − σ, and σ = (σ ζ + σ η )/2. In particular, Equation (1) after integrating becomes: The purely elastic solution should satisfy the boundary conditions (20) and (21). As a result, it follows from (23) that: and (25) show that plastic yielding initiates at ζ = 0. The corresponding values of a and s are denoted as a e and s e , respectively. One can find from (17) and (23) that 4s e = exp α/ √ 3 . Replacing s with s e in the second equation in (25) and solving the resulting equation for a provides a e . In what follows, it is assumed that a ≥ a e .

Elastic/Plastic Solution with One Plastic Region (Stage 2 of the Process)
There are one elastic region, −1 ≤ ζ ≤ ζ 1 , and one plastic region, ζ 1 ≤ ζ ≤ 0, during this stage of the process. Thus ζ 1 is the elastic/plastic boundary. It is evident that ζ 1 is a function of a and ζ 1 = 0 at a = a e . The solution (23) is valid in the elastic region, but (25) is not. The boundary condition (21) and the solution (23) combine to give: It follows from (23) and (26) that: at ζ = ζ 1 .
The distribution of σ ζ and σ η in the plastic region is found from (17) and (18). In this region, σ η > σ ζ and ξ p η > 0. Therefore, Equation (15) becomes: Moreover, using (28), one transforms Equation (17) to: Equations (18) and (29) combine to give: where β = K/σ 0 . The elastic strains are determined from (11), (14) and (28) as: Equations (24), (29) and (31) combine to give: It follows from this equation that where Φ ε p eq = dΦ/dε p eq . One replaces differentiation with respect to ζ with differentiation with respect to ε p eq in (30) using (33) to transform (30) to Let ε 1 be the value of ε p eq at ζ = 0. Then, the boundary condition (20) becomes: for ε p eq = ε 1 . Both σ ζ and σ η should be continuous across the elastic/plastic boundary. Moreover, ε p eq = 0 at this boundary. Therefore, the solution of Equation (34) satisfying the condition that σ ζ is continuous is where Ψ ε p eq is the antiderivative of Φ ε p eq and p 2 should be eliminated using (27). The stress σ η is determined from (29) and (36) as: The continuity of σ η across the elastic/plastic boundary is equivalent to the continuity of σ η − σ ζ . The latter combined with (27) and (29) at ε p eq = 0 leads to Using the boundary condition (35), one finds from (36) that It follows from (32) that Using this equation, one eliminates s in the expression for p 2 given in (27). As a result, Equation (39) connects a and ε 1 . This equation should be solved numerically to find the dependence of ε 1 on a. Then, s, ζ 1 , and D 0 as functions of a are readily determined from (40), (38) and (26). The equations in (23) supply the distribution of the stresses in the elastic region. Equations (32), (36) and (37) provide the distribution of the stresses in the plastic region in parametric form with ε p eq being the parameter. The bending moment is determined from (22) by numerical integration.
This stage of the process ends when the yield criterion starts to satisfy at ζ = −1. It follows from (23) and (17) Since s as a function of a has already been found, this equation allows one to determine the value of a corresponding to plastic yielding initiation at ζ = −1. This value is denoted as a 1 .

Elastic/plastic Solution with Two Plastic Regions (Stage 3 of the Process)
There are one elastic region, ζ 2 ≤ ζ ≤ ζ 1 , and two plastic regions, ζ 1 ≤ ζ ≤ 0 and −1 ≤ ζ ≤ ζ 2 , at a > a 1 . Thus ζ 1 and ζ 2 are the elastic/plastic boundaries. Both are functions of a and ζ 2 = −1 at a = a 1 . The solution (23) is valid in the elastic region, but (25) and (26) are not. Moreover, the general solution in the plastic region ζ 1 ≤ ζ ≤ 0 provided in Section 4.2 is also valid. Consider the plastic region −1 ≤ ζ ≤ ζ 2 . In this region, σ η < σ ζ and ξ p ζ > 0. Therefore, Equation (15) becomes: Moreover, using (42), one transforms Equation (17) to: Equations (18) and (43) combine to give: The elastic strains are determined from (11), (14) and (42) as: Equations (24), (43) and (45) combine to give: It follows from this equation that: One replaces differentiation with respect to ζ with differentiation with respect to ε p eq in (44) using (47) to transform (44) to: Let ε 2 be the value of ε p eq at ζ = −1. Then, the boundary condition (21) becomes: for ε p eq = ε 2 . Both σ ζ and σ η should be continuous across the elastic/plastic boundary ζ = ζ 2 . Moreover, ε p eq = 0 at this boundary. Therefore, the solution of Equation (48) satisfying the condition that σ ζ is continuous is: where p 1 is determined from (23) as: The stress σ η is found from (43) and (50) as: The continuity of σ η across the elastic/plastic boundary is equivalent to the continuity of σ η − σ ζ . The latter combined with (23) and (43) at ε p eq = 0 leads to: Using the boundary condition (49), one finds from (50) that: It follows from (46) that: The first equation in (23) supplies the value of p 2 involved in (36) and (37) in the form: Here Equation (38) has been taken into account. Equations (51) and (53) combine to give: It is seen from (56) and (57) that p 1 = p 2 during this stage of the process. Therefore, one can replace p 2 with p 1 in (36), (37) and (39). Then, one eliminates p 1 between (39) and (54) to get: Equations (40) and (55) are solved for a and s as: Using (59), one eliminates a and s in (58). The resulting equation connects ε 2 and ε 1 . This equation should be solved numerically to find the dependence of ε 1 on ε 2 . Then, a, s, ζ 1 , ζ 2 , p 1 , and D 0 as functions of ε 2 are readily determined from (59), (38), (53), (54), and (57), respectively. It is evident that it is also possible to express ε 1 , ε 2 , s, ζ 1 , ζ 2 , p 1 , and D 0 as functions of a using this solution. The equations in (23) supply the distribution of the stresses in the elastic region. Equations (32), (36) and (37) provide the distribution of the stresses in the plastic region ζ 1 ≤ ζ ≤ 0 in parametric form with ε p eq being the parameter. Equations (46), (50) and (52) provide the distribution of the stresses in the plastic region −1 ≤ ζ ≤ ζ 2 in parametric form with ε p eq being the parameter. The bending moment is determined from (22) by numerical integration.
The solution above is valid if the size of each plastic region increases as the deformation proceeds. The condition is equivalent to: These inequalities should be verified numerically. It is worthy to note that the process' analysis beyond the instant when one of the conditions in (60) is violated is possible but requires additional constitutive equations that describe reversed plastic yielding.

Residual Stresses
Consider elastic unloading, assuming that the strains and the displacements are small at this stage. Therefore, variations of the shape can be neglected. Denote r CD = R 0 and r AB = R 1 at the end of loading when a = a l , s = s l , p = p l , h = h l and m = m l . All these values are calculated using the solution given in the previous section. The distribution of σ ζ and σ η corresponding to a = a l is denoted as σ (l) ζ and σ (l) η . It is convenient to introduce a polar coordinate system (r, θ) by the following transformation equations: It is evident from these equations that the coordinate curves of the (r, θ)-coordinate system coincide with the coordinate curves of the (ζ, η)-coordinate system. Therefore, σ (l) ζ = σ r and σ (l) η = σ θ . Using (61), one can represent the solution given in the previous section in the polar coordinate system with no difficulty.
The general solution for the increment of the principal stresses from the configuration corresponding to a = a l is independent of the solution at loading. Therefore, one may adopt the solution provided in [17]. According to this solution: where ρ 0 = R 0 /R 1 . The radial distribution of the residual stresses can be found from the equations σ res r = σ res ζ = σ r + σ * r and σ res θ = σ θ + σ * θ . The solution in this section is valid if no reverse yielding occurs. Using Equation (17), the corresponding condition can be written in the form: In this equation, ε p eq and ε p ζ are understood to be calculated at a = a l .

Numerical Example
Even though the solution found is semianalytic, its full parametric analysis is not feasible due to the significant number of parameters that affect the solution. Since the solution for isotropic hardening laws has already been discussed in [17], the numerical example below focuses on kinematic hardening. The numerical integration and solution of transcendental equations have been performed using the corresponding build-in commands in Mathematica (version 11.3) [22].
Swift's law describes isotropic hardening. Then: It is assumed that n = 0.25 and ε 0 = 0.222. In all calculations, α = 3 × 10 −3 . The parameter β varies in the range 0 ≤ β ≤ 10. The material with no kinematic hardening is obtained at β = 0. The through-thickness distribution of the principal stresses and the bending moment are calculated according to the procedure described in Section 4.
The inequalities in (60) control the range of validity of the solution at loading. Figure 2 demonstrates the variation of ζ 1 and ζ 2 with H/r CD at β = 5 for several values of f. The convex down lines correspond to ζ 1 and the convex up lines to ζ 2 . Stars denote the extremum points. It is evident that the derivatives with respect to a and H/r CD have the same sign and vanish simultaneously. It is seen from Figure 2 that the solution may break down because (i) dζ 1 /da = 0 and dζ 2 /da > 0, or (ii) dζ 2 /da = 0 and dζ 1 /da < 0, or (iii) dζ 1 /da = 0 and dζ 2 /da = 0 at a certain value of H/r CD . These conditions are taken into account in the numerical example below.   Since ζ = ζ 1 and ζ = ζ 2 are the elastic/plastic boundaries, Figure 2 illustrates the elastic and plastic regions' thickness in the Lagrangian coordinates. In particular, the thickness of the plastic region adjacent to the outside surface is |ζ 1 |, the thickness of the plastic region adjacent to the inside surface is (1 + ζ 2 ), and the thickness of the elastic region is (ζ 1 − ζ 2 ). It is seen from this figure that the thickness of the elastic region decreases quickly at the very beginning of the process and, as a result, becomes small after a small amount of deformation. Figures 3 and 4 depict the variation of the dimensionless bending moment with H/r CD at f = 0 (pure bending) and f = 0.2, respectively. Each figure shows several curves corresponding to different values of β. It is evident that m = 0 and H/r CD = 0 at the initial instant. Figures 3 and 4 reveal a significant effect of β on the dimensionless bending moment. In general, the material with no kinematic hardening requires the least moment for the deformation to proceed.         This value of corresponds r CD to the third stage of the process. In these figures, X is the dimensionless distance from the surface ζ = −1 defined as: The distribution of the stress σ ζ at f = 0 and f = 0.2 is depicted in Figures 5 and 7, respectively. Figures 6 and 8 show the distribution of the stress σ η at the same set of parameters. It is seen from these figures that the thickness of the elastic region is quite narrow at this stage of the process. The effect of β on the magnitude of the stress σ η is not significant in some vicinity of the elastic/plastic boundaries and increases as the distance from these boundaries increases. The effect is most pronounced and significant at the surfaces ζ = −1 and ζ = 0 and in the elastic region's central part. The effect of β on the magnitude of the stress σ ζ is most significant in the sheet's central part. It is not surprising since the magnitude of this stress is prescribed at ζ = −1 and ζ = 0 by the boundary conditions.           σ sharply increases as X increases in the region that is purely elastic at the loading stage. This stress significantly decreases as X increases in the two other regions and vanishes at two points. One of these points lies in the region that is purely elastic at the loading stage, and the other in the region 1 0 ζ < ζ < . The latter region is plastic at the loading stage.  ζ with X at f = 0 and f = 0.2 is depicted in Figures 9 and 11, respectively. Figures 10 and 12 show the through-thickness distribution of the stress σ res η . The magnitude of the stress σ res ζ is quite small. The magnitude of the stress σ res η sharply increases as X increases in the region that is purely elastic at the loading stage. This stress significantly decreases as X increases in the two other regions and vanishes at two points. One of these points lies in the region that is purely elastic at the loading stage, and the other in the region ζ 1 < ζ < 0. The latter region is plastic at the loading stage.           The distribution of the residual stresses does not satisfy the inequality (64) if β is large enough. It is evident that one can represent Λ involved in (64) as a function of X. This function is depicted in Figure 13 for f = 0 and Figure 14 for f = 0.2. It is seen from Figure 14 that Λ practically vanishes at one point if β = 6. Therefore, the solution for the residual stresses is not valid for β > 6. For this reason, the curves for β = 8 and β = 10 are not shown in Figures 11 and 12. For consistency, such curves are not depicted in Figures 9 and 10.

Discussion
A semianalytic solution describing the process of finite plane strain bending under tension of an elastic/plastic sheet has been derived assuming combined isotropic and kinematic hardening. The material is supposed to be incompressible. The solution has been found by an inverse method starting from the transformation equations between Eulerian and Lagrangian coordinates given in Equation (9). In general, the process consists of three stages. The sheet is purely elastic during the first stage. This stage ends when a plastic region starts to spread from the surface AB to the surface CD (Figure 1). The second stage of the process ends when another plastic region starts to spread from the surface CD to the surface AB (Figure 1). The solution is valid if the inequalities in Equation (60) are satisfied. Figure 2 shows that these solution nonexistence conditions depend on process parameters. The process of unloading has also been considered assuming that no reverse plastic region appears. Equation (64) controls the latter condition.
The solution is semianalytic but includes several parameters and an arbitrary function Φ ε p eq , its derivative and antiderivative. Therefore, the complete parametric analysis of the solution is not feasible, but it is straightforward to find any quantity of interest for any given set of parameters and the function Φ ε p eq . In particular, the derivative and antiderivative of Φ ε p eq can be found in terms of elementary functions for all widely used isotropic hardening laws.
One of the limiting cases of the process considered in the present paper is tension with no bending. In this case, there is one elastic region at the beginning of the process. This region is replaced with a plastic region at a certain value of the tensile force. The distribution of stresses and strains is uniform in this limiting case. If the bending moment is small, then Stage 3 of the process never occurs. The solutions given in Sections 4.1 and 4.2 are valid. As the deformation proceeds, the elastic/plastic boundary may reach the inside surface of the sheet, ζ 1 = −1. If the deformation further proceeds, the plastic region occupies the entire sheet.
The numerical example has been provided for the function Φ ε p eq shown in Equation (65). Figures 3-14 demonstrate that the parameter β involved in Prager's law significantly affects the bending moment and the through-thickness distribution of the stresses at loading but slightly affects the through-thickness distribution of the residual stresses. Of course, this conclusion is only valid for the parameters used in the numerical solution.
The solution found in conjunction with experimental research can be used for identifying constitutive models. Various bending tests have already been employed for this purpose [2][3][4][5][6][7]. An advantage of the present solution is that it is given in closed form, except simple numerical techniques for calculating the bending moment and solving transcendental equations, for arbitrary function Φ ε p eq involved in Equation (2). One can vary this function for improving the agreement between the theory and experiment.
The solution can be used as a benchmark problem for verifying numerical codes, which a necessary step before using such codes [23,24].

Conclusions
The following conclusions were reached from the theoretical analysis of plane strain bending under tension of a wide sheet: (1) The constitutive equations comprising combined isotropic and kinematic hardening laws and the normality rule permit a semianalytic solution. No restriction, other than the conventional restrictions from the general theory of plasticity, is imposed on the isotropic hardening law. Prager's law describes kinematic hardening. Numerical techniques are only necessary for solving transcendental equations and evaluating ordinary integrals. (2) An inverse method is used for finding the solution. Its key point is transformation equations that connect Eulerian and Lagrangian coordinates. These equations describe the geometry of the process and are independent of the stress equations. (3) The distribution of residual stresses is also given, assuming that unloading is purely elastic. (4) An applied aspect of the solution is that it can be used in conjunction with experimental data for identifying material properties. An advantage of the solution in this respect is that it is valid for the arbitrary isotropic hardening law. (5) One can use the solution for verifying numerical codes, which is a necessary step for their subsequent use.