Analysis of Strain-Hardening Viscoplastic Wide Sheets Subject to Bending under Tension

: The present paper provides an accurate solution for ﬁnite plane strain bending under tension of a rigid/plastic sheet using a general material model of a strain-hardening viscoplastic material. In particular, no restriction is imposed on the dependence of the yield stress on the equivalent strain and the equivalent strain rate. A special numerical procedure is necessary to solve a non-standard ordinary differential equation resulting from the analytic treatment of the boundary value problem. A numerical example illustrates the general solution assuming that the tensile force vanishes. This numerical solution demonstrates a signiﬁcant effect of the parameter that controls the loading speed on the bending moment and the through-thickness distribution of stresses.


Introduction
Bending is frequently used as a test for determining various material properties. A four-point bend apparatus was used in [1] for getting stress-strain curves in the tension and compression of three materials: beryllium, cast iron, and copper. No significant difference between these curves and the corresponding curves from the conventional uniaxial tension and compression tests was found. Papers [2,3] applied the same method for determining the stress-strain curves for annealed copper and cement-based composites, respectively. The theoretical solution proposed in [1] was modified in [4]. Then, the method was used for determining the stress-strain curves for pure magnesium and S45C steel. Paper [5] studied the elastic/plastic behavior and the failure of CLARE laminates in bending experimentally. An elastic/plastic material model was proposed using these experimental data. Lightweight-aggregate concrete beams were tested in bending to failure in [6]. This work emphasized the location of the neutral axis at failure. Adhesively bonded bending specimens were employed in [7] for determining the bilinear elastic/plastic shear stress-strain behavior of acrylic adhesives. Hybrid sandwich structures were tested in bending in [8] to construct failure mode maps.
Experimental data should be usually supplemented with theoretical solutions for identifying material models. Most semi-analytic solutions have dealt with strain hardening materials [9][10][11][12]. Paper [12] also accounted for plastic anisotropy and tension/compression asymmetry. The representation of strain or work hardening in the strain reversal region has been simplified in these works. Paper [13] developed a semi-analytic method to treat this region without any simplification. This method was combined with various constitutive equations to describe pure bending and bending under tension of wide sheets [14]. Under certain conditions, viscoplastic or strain-hardening viscoplastic constitutive equations are required to analyze bending processes adequately [15][16][17]. Many material models of this type are available in the literature [18][19][20][21][22][23][24]. For identifying constitutive equations, it is desirable to have a theoretical solution for a general material model. The motivation of the present paper is to provide a semi-analytic solution of the same level of complexity as available solutions but for quite a general material model. It is believed that such a solution is useful for identifying constitutive equations. The solution is based on the general method proposed in [13].
The solution found is facilitated by using the equivalent strain instead of the space coordinate as an independent variable. Similar changes of independent variables have been successfully used in several works [25][26][27]. In particular, the process of bending was analyzed in [26]. However, in all these cases, it has been assumed that the yield stress depends on the single kinematic quantity, the equivalent plastic strain. The novelty of the present solution is that the yield stress is an arbitrary function of the equivalent strain and the equivalent strain rate. An important property of the new solution is that process parameters depend on the loading speed. In addition, the solution shows that some assumptions used in simplified formulations are not valid and corresponding numerical results may lead to the misinterpretation of experimental data.

Statement of the Problem
What is considered is the process of bending under tension due to large strains under plane strain conditions in the classical formulation proposed in [28], where a model of rigid perfectly plastic material was adopted. A schematic diagram of the process is shown in Figure 1. The cross section of a sheet is rectangular at the initial instant. The initial thickness of the sheet is H. As the deformation proceeds, the initially straight lines A 1 B 1 and C 1 D 1 become circular arcs AB and CD, respectively. The initially straight lines B 1 C 1 and A 1 D 1 remain straight (lines BC and AD, respectively). The radii of circular arcs AB and CD are denoted as r AB and r CD , respectively. A Cartesian coordinate system (x, y) is introduced. Its x-axis coincides with the axis of symmetry of the process. Its origin coincides with the intersection of the outer radius of the sheet and the axis of symmetry. Due to symmetry, it is sufficient to consider the domain y ≥ 0. The orientation of line BC relative to the x-axis is θ 0 . The actual distribution of stresses over BC is replaced with the bending moment M and the tensile force F. Both are per unit length. The surface AB is traction-free. Equilibrium requires that the tensile force is compensated with pressure over the surface CD. It is assumed that this pressure is uniformly distributed. equations, it is desirable to have a theoretical solution for a general material model. The motivation of the present paper is to provide a semi-analytic solution of the same level of complexity as available solutions but for quite a general material model. It is believed that such a solution is useful for identifying constitutive equations. The solution is based on the general method proposed in [13]. The solution found is facilitated by using the equivalent strain instead of the space coordinate as an independent variable. Similar changes of independent variables have been successfully used in several works [25][26][27]. In particular, the process of bending was analyzed in [26]. However, in all these cases, it has been assumed that the yield stress depends on the single kinematic quantity, the equivalent plastic strain. The novelty of the present solution is that the yield stress is an arbitrary function of the equivalent strain and the equivalent strain rate. An important property of the new solution is that process parameters depend on the loading speed. In addition, the solution shows that some assumptions used in simplified formulations are not valid and corresponding numerical results may lead to the misinterpretation of experimental data.

Statement of the Problem
What is considered is the process of bending under tension due to large strains under plane strain conditions in the classical formulation proposed in [28], where a model of rigid perfectly plastic material was adopted. A schematic diagram of the process is shown in Figure 1. The cross section of a sheet is rectangular at the initial instant. The initial thickness of the sheet is H. As the deformation proceeds, the initially straight lines 1 x y is introduced. Its x-axis coincides with the axis of symmetry of the process. Its origin coincides with the intersection of the outer radius of the sheet and the axis of symmetry. Due to symmetry, it is sufficient to consider the domain 0 y  . The orientation of line BC relative to the x-axis is 0  . The actual distribution of stresses over BC is replaced with the bending moment M and the tensile force F. Both are per unit length. The surface AB is traction-free. Equilibrium requires that the tensile force is compensated with pressure over the surface CD . It is assumed that this pressure is uniformly distributed.
(a) (b)  Then, The constitutive equations include an isotropic pressure-independent yield criterion (i.e., the linear invariant of the stress tensor does not affect plastic yielding) and its associated flow rule. The elastic portion of the strain tensor is neglected. The present paper assumes that the tensile yield stress is a function of the equivalent strain, ε eq , and the equivalent strain rate, ξ eq , obtained by fitting any experimental data. Then, the plane-strain yield criterion can be represented as Here σ 1 and σ 2 are the principal stresses, σ 0 is a reference stress, and Φ ε eq , ξ eq is an arbitrary function of its arguments. The associated flow rule is equivalent to the equation of incompressibility and the condition that the stress and strain rate tensors are coaxial. In the case under consideration, the equivalent strain rate is determined as where ξ 1 and ξ 2 are the principal strain rates. The equivalent plastic strain is given by here d/dt denotes the convected derivative.

Kinematics
The approach for solving plane-strain bending problems proposed in [13] consists of two major steps. The first step describes the kinematics of the process. This step is the same for all isotropic incompressible materials. The second step concerns the determination of the stresses, bending moment, and tensile force. This step depends on the constitutive equations chosen. The present section outlines the basic relations from the first step derived in [13]. These relations are necessary for the second step.
Let (ζ, η) be the Lagrangian coordinate system satisfying the conditions x = Hζ and y = Hη at the initial instant. Then, ζ = 0 on AB and ζ = −1 on CD throughout the process of deformation ( Figure 1). The transformation equations between the (ζ, η)− and (x, y)− coordinate systems are Here a is a time-like variable and s is a function of a. At the initial instant, a = 0. The function s(a) should be found from the solution. This function must satisfy the condition at a = 0. The coordinate lines of the Lagrangian coordinate system coincide with principal strain rate trajectories [13]. Then, the coordinate lines of the Lagrangian coordinate system coincide with principal stress trajectories for the coaxial models. Therefore, the contour of the cross-section is free of shear stresses throughout the process of deformation. It is convenient to introduce a cylindrical coordinate system (r, θ) by the transformation equations: x Equations (5) and (7) combine to give Therefore, the radial, σ r , and circumferential, σ θ , stresses are the principal stresses. In what follows, it is assumed that σ r ≡ σ 1 and σ θ ≡ σ 2 . Accordingly, ξ r ≡ ξ 1 and ξ θ ≡ ξ 2 . Here ξ r is the radial strain rate and ξ θ is the circumferential strain rate. It follows from (8) that It is straightforward to find the principal strain rates using (5). In particular, Equations (3) and (10) combine to give The equivalent strain rate vanishes at the neutral line. Therefore, it follows from (11) that the ζ− coordinate of the neutral line is The neutral line moves through the material as the deformation proceeds. It is known that it moves towards CD in the case of strain-hardening materials. Assume that it is true in the case of the constitutive equations under consideration. Specific calculations below will verify this assumption. The initial condition for Equation (4) is at a = 0. Let ζ n0 be the value of ζ n at a = 0. It is necessary to consider three regions for determining the equivalent strain. These regions are: 0 ≥ ζ ≥ ζ n0 (Region 1), ζ n0 ≥ ζ ≥ ζ n (Region 2), and ζ n ≥ ζ ≥ 0 (Region 3). In Region 1, ζ − ζ n = ζ + ds/da ≥ 0 throughout the process of deformation. Therefore, one can substitute (11) into (4) and immediately integrate the resulting equation using (13) and (6) to get In Region 3, ζ − ζ n = ζ + ds/da ≤ 0 throughout the process of deformation. Therefore, one can substitute (11) into (4) and immediately integrate the resulting equation using (13) and (6) to get In Region 2, ζ − ζ n = ζ + ds/da ≥ 0 instantly. However, strain reversal occurs in this region. Consider any ζ− curve of Region 2. Let a c be the value of a at which this curve coincides with the neutral line and s c is the corresponding value of s. Evidently, both a c and s c are functions of ζ. Then, the equivalent strain in Region 2 is given by

Stress Solution
The only stress equilibrium equation that is not satisfied automatically in the cylindrical coordinate system is ∂σ r ∂r Metals 2022, 12, 118
Equations (24) and (25) combine to give Let ε 1 be the value of ε eq on surface AB (Figure 1). It follows from (14) that Since σ r = 0 on surface AB, the distribution of the radial stress in region 1 is determined from (26) and (27) as Here µ is a dummy variable of integration and Λ 1 ε eq denotes the function Λ ε eq when β = +1 in (23). The value of ε eq at ζ = ζ n0 is determined from (14) as Then, the radial stress at ζ = ζ n0 is found from (28) in the form Let ε 3 be the value of ε eq on surface CD (Figure 1). It follows from (15) that Using (9), one can represent (1) as where f = F/(σ 0 H). It follows from (26), (31), and (32) that the distribution of the radial stress in region 3 is Here Λ 3 ε eq denotes the function Λ ε eq when β = −1 in (23). The value of ε eq at ζ = ζ n is determined from (12) and (15) as Then, the radial stress at ζ = ζ n is found from (33) in the form In region 2, Equation (20) becomes Here the function Λ 2 (ζ) is determined by replacing ε eq and χ in Φ ε eq , ξ 0 χ with the right-hand sides of Equations (16) and (21), respectively. The radial stress must be continuous at ζ = ζ n and ζ = ζ n0 . Therefore, using (35), one can represent the solution of Equation (36) as Equations (30) and (37) combine to give This equation connects a, s, and ds/da. Therefore, it can be considered an ordinary differential equation to find s as a function of a. Its solution must satisfy the condition (6).
Once Equation (38) has been solved, the distribution of the radial stress in all three regions is readily determined from (28), (33), and (37). Equations (28) and (33) should be used in conjunction with Equations (14) and (15), respectively. The distribution of the circumferential stress is found from (18). The bending moment is given by Here h is the current thickness of the sheet, h = r AB − r CD (Figure 1). Using (8), one can rewrite (39) as Note that m = 1 in the case of rigid perfectly plastic material if f = 0 [28].

Numerical Method
It has been noted that Equation (38) is an ordinary differential equation. Several well-established methods are available for solving such equations numerically (see, for example, [29]). However, the form of the differential equation in Equation (38) is nonstandard. Therefore, a special numerical method should be developed to solve this equation. The method below is a modification of the method developed in [30] for solving a system of partial differential equations.
Equation (38), together with Equation (6), constitutes the initial value problem. However, a difficulty is that the equation provides no information at the initial instant. In particular, each term of this equation vanishes at a = 0 since ε eq = 0 everywhere. It is a consequence of rigid plasticity. Such difficulty does not arise in elastic/plastic solutions (see, for example, [31]). Below, the difficulty noted is resolved by finding the state of stress at the initial instant without using Equation (38) and developing a general iterative procedure.
Assume that s = s i and ds/da = q i at a = a i , and s i and q i are known. Let ∆a be a small increment of a. Then, a i+1 = a i + ∆a, and the value of s at a = a i+1 can be approximated as Here q i+1 is the value of the derivative ds/da at a = a i+1 , and its value is unknown. Using (27), (29), (31), (34), and (41), one can express the limits of integration involved in (38) in terms of q i+1 and known quantities. Thus Equation (38) determines q i+1 . The procedure above can be repeated for a = a i+2 = a i+1 + ∆a and subsequent values of a.
It remains to find s and ds/da at a = 0 to start the iterative procedure above. The required value of s is given by (6). The initial cross-section of the sheet is shown in Figure 1a. The through-thickness stress vanishes at this instant. The equivalent strain vanishes everywhere. Equation (21) becomes Here q 0 is the value of the derivative ds/da at a = 0. Equations (2) and (42) combine to give The tensile force at the initial instant is determined as Substituting (43) into (44) gives This equation determines q 0 . In particular, q 0 = 1/2 if f = 0 (pure bending). The dimensionless bending moment can be found using (40), where σ θ should be replaced with σ η .

Numerical Example
The numerical results below are for pure bending. Several functions Φ ε eq , ξ eq are available in the literature [18][19][20][21][22][23][24]. Since bending of strain-hardening materials, including quantitative results, has been investigated in many papers, the numerical example below emphasizes the effect of loading speed on process parameters.
In most cases, the function Φ ε eq , ξ eq is represented as a product of two functions. One of these functions depends on the equivalent strain and the other on the equivalent strain rate. The illustrative example in this section assumes that Φ ε eq , ξ eq is the product of Swift's strain-hardening law and Herschel-Bulkley's viscoplastic law. Then, In all calculations, λ = 0.73, ν = 0.17, and n = 0.2. The characteristic strain rate is involved in the solution only through γ. In particular, one can rewrite (46) as Φ ε eq , χξ 0 = 1 + λε eq ν (1 + χ n ).
The numerical solution has been found using the procedure described in the previous section. It demonstrates the effect of γ on the through-thickness distribution of stresses and the bending moment. Figures 2 and 3 depict the through-thickness distributions of the radial and circumferential stresses, respectively, at several process stages identified by the inner radius of the sheet. These distributions correspond to γ = 0.1. In these figures,  Figure 6 demonstrates the effect of  on the dimensionless bending moment. The broken line corresponds to the rate-independent material model (i.e., the material obeys Swift's strain-hardening law). The bending moment increases as  increases, which is in agreement with physical expectations. The increase in the bending moment is sharper for smaller values of  . It is worthy of note that the bending moment is nearly constant for a   Figure 6 demonstrates the effect of  on the dimensionless bending moment. The broken line corresponds to the rate-independent material model (i.e., the material obeys Swift's strain-hardening law). The bending moment increases as  increases, which is in agreement with physical expectations. The increase in the bending moment is sharper for smaller values of  . It is worthy of note that the bending moment is nearly constant for a given value of  .  Thus Z = 0 corresponds to the inner radius of the sheet and Z = h/H to its outer radius. Using (8) and (9), one transforms Equation (48) to At the final stage of the process, the inner radius of the sheet is rather small (smaller than the initial thickness of the sheet). Figures 4 and 5 show the same distributions at γ = 1. The circumferential stress is discontinuous at the neutral line. Figures 3 and 5 show that the neutral line moves to the inner surface of the sheet as the deformation proceeds. Thus the assumption made in Section 3 has been confirmed. The magnitude of the radial stress is quite high near the neutral line, except at the very beginning of the process (see Table 1). Therefore, the assumption used in simplified solutions [32,33] that the radial stress is negligible is not valid. The rightmost points of the curves in Figures 2-5 correspond to Z = h/H. Therefore, it is seen from these figures that the change in the thickness is negligible.         Figure 6 demonstrates the effect of γ on the dimensionless bending moment. The broken line corresponds to the rate-independent material model (i.e., the material obeys Swift's strain-hardening law). The bending moment increases as γ increases, which is in agreement with physical expectations. The increase in the bending moment is sharper for smaller values of γ. It is worthy of note that the bending moment is nearly constant for a given value of γ.

Conclusions
An accurate rigid plastic solution for a general strain-hardening viscoplastic material model has been found. The solution describes the plane strain process of bending under tension due to large strains. The elastic portion of the strain tensor has been neglected. However, the effect of elasticity on the accuracy of solutions is negligible except at the beginning of the process [34].
The numerical solution reduces to an ordinary differential equation in a non-standard form. A special numerical procedure is required to solve this equation. The procedure is described in Section 5. It is used to determine the stress field and the bending moment in the pure bending of a sheet made of a material obeying Equation (46). This numerical solution is presented and discussed in detail in Section 6.
The solution presented may be considered a starting point for designing experiments to identify constitutive equations. It is advantageous in this respect that the general solution is not restricted to a specific function or even a class of functions but is valid for an arbitrary dependence of the yield stress on the equivalent strain and the equivalent strain rate.
The solution can be used as a benchmark problem for verifying numerical codes. This necessary step applies to different constitutive equations, including equations used in metal forming applications [35][36][37].