Finite Pure Plane Strain Bending of Inhomogeneous Anisotropic Sheets

Abstract

The present paper concerns the general solution for finite plane strain pure bending of incompressible, orthotropic sheets. In contrast to available solutions, the new solution is valid for inhomogeneous distributions of plastic properties. The solution is semi-analytic. A numerical treatment is only necessary for solving transcendent equations and evaluating ordinary integrals. The solution’s starting point is a transformation between Eulerian and Lagrangian coordinates that is valid for a wide class of constitutive equations. The symmetric distribution relative to the center line of the sheet is separately treated where it is advantageous. It is shown that this type of symmetry simplifies the solution. Hill’s quadratic yield criterion is adopted. Both elastic/plastic and rigid/plastic solutions are derived. Elastic unloading is also considered, and it is shown that reverse plastic yielding occurs at a relatively large inside radius. An illustrative example uses real experimental data. The distribution of plastic properties is symmetric in this example. It is shown that the difference between the elastic/plastic and rigid/plastic solutions is negligible, except at the very beginning of the process. However, the rigid/plastic solution is much simpler and, therefore, can be recommended for practical use at large strains, including calculating the residual stresses.

1. Introduction

Sheet metal forming processes usually include bending. A brief review of typical sheet metal forming processes that incorporate bending is provided in [1]. The process of bending is also an essential test for identifying material properties, for example [2,3,4,5,6,7,8,9,10]. Theoretical analyses of the bending process are necessary for interpreting test results.

An exact rigid perfectly plastic solution for finite pure plane strain bending and plane strain bending under tension of sheets has been found in [11]. This solution has been adopted in [12] for deriving a closed form expression for strain at any fiber. This paper has concluded that the basic assumptions made in [11] are plausible. The solution [11] has been extended to many material models. Of primary interest for the present paper are solutions for anisotropic materials and inhomogeneous sheets. A brief review of such solutions is given below.

Elastic and plastic anisotropy is a typical property of many metallic and non-metallic materials [13,14,15,16]. A solution for pure plane strain bending of anisotropic sheets has been derived in [17]. The model adopted incorporates the Bauschinger effect and strain hardening. Plastic anisotropy is described by Hill’s quadratic yield criterion [11]. It has been found that the effect of plastic anisotropy on the bending moment is rather significant. Papers [18,19] present an analysis of elastic/plastic bending of orthotropic sheets. Elastic properties are isotropic, and plastic yielding obeys Hill’s quadratic yield criterion [11]. A more detailed review of solutions for plane strain bending of anisotropic sheets is provided in [20]. An extension of the solutions above has been proposed in [21] where tension/compression asymmetry of plastic properties has been considered. All these solutions are for homogeneous sheets.

Several solutions are available for piecewise homogenous sheets. The pure bending of bonded laminated metals under plane strain conditions has been considered in [22,23]. Different rigid plastic material models have been adopted in these papers. Elastic properties have been taken into consideration in [24]. This paper emphasizes the prediction of springback. The distribution of residual stresses in bilayer sheets after bending has been found in [25]. A simplified solution for bending and subsequent unloading of bilayer sheets has been proposed in [26].

Different methods have been employed to get the solutions above. A unified approach for analyzing finite pure plane strain bending has been developed in [27]. The approach applies to a broad class of incompressible materials. In particular, the corresponding solutions for anisotropic and bilayer sheets have been found in [28] and [29], respectively.

In many cases, the through-thickness distribution of material properties in sheets is non-uniform but, in contrast to [22,23,24,25,26], is described by continuous functions [30]. Using the approach [27], the pure bending of isotropic elastic/plastic functionally graded sheets has been analyzed in [31]. The present paper extends this solution to anisotropic sheets. In addition, a rigid/plastic model is considered. The general solution involves quite an arbitrary through-thickness distribution of plastic properties. The numerical example is based on the experimental data presented in [32]. In the case of the rigid/plastic solution, the prediction of residual stresses can be made using the methodology proposed in [33] and further developed and discussed in [34]. According to this methodology, for computing residual stresses in a deformation process, the elasticity is neglected during the process’s loading phase. It is shown that the solution for purely elastic unloading is not valid even for a relatively large inside radius.

2. Statement of the Problem

A metallic sheet is bent under plane strain conditions by two couples M, as shown in Figure 1. The initial shape of the sheet in the planes of flow is a rectangular of thickness H and width 2L.

Curves AB and CD are circular arcs, and AD and BC are straight throughout the process of deformation (except the initial instant when AB and CD are straight). It has been shown in [27] that the following equations describe the transformation of the initial shape into the shape after any amount of strain:

$$\frac{x}{H}=\sqrt{\frac{\zeta }{a}+\frac{s}{a^2}}\cos \left(2a\eta \right)-\frac{\sqrt{s}}{a}\hspace{1em}\mathrm{and}\hspace{1em}\frac{y}{H}=\sqrt{\frac{\zeta }{a}+\frac{s}{a^2}}\sin \left(2a\eta \right).$$

(1)

Here $\left(x,\hspace{0.17em}y\right)$ are Eulerian Cartesian coordinates and $\left(\zeta ,\hspace{0.17em}\eta \right)$ are Lagrangian coordinates. Additionally, a is a time-like variable such that $a=0$ at the initial instant and s is a function of a. The latter should be found from the solution. The x-axis is an axis of symmetry of the process. The Cartesian coordinate system’s origin is situated at the intersection of the axis of symmetry and curve AB. At the initial instant,

$$x=\zeta H\hspace{1em}\mathrm{and}\hspace{1em}y=\eta H.$$

(2)

This condition is satisfied if

$$s=\frac{1}{4}$$

(3)

at $a=0$. It is convenient to introduce a plane polar coordinate system $\left(r,\hspace{0.17em}\theta \right)$ with the origin at $x=-H\sqrt{s}/a$ and $y=0$. It follows from (1) that

$$\frac{r}{H}=\frac{\sqrt{\zeta a+s}}{a}\hspace{1em}\mathrm{and}\hspace{1em}\theta =2a\eta .$$

(4)

It is seen from (2) and Figure 1 that $\zeta =0$ on curve AB and $\zeta =-1$ on curve CD throughout the process of deformation. The radii of circular arcs AB and CD are determined from (4) as

$$\frac{r_{AB}}{H}=\frac{\sqrt{s}}{a} \mathrm{and}\hspace{1em}\frac{r_{CD}}{H}=\frac{\sqrt{s-a}}{a},$$

(5)

respectively. The transformation (1) satisfies the equation of incompressibility.

Let ${\sigma }_{\zeta }$ and ${\sigma }_{\eta }$ be the normal stresses referred to the Lagrangian coordinate system. This coordinate system is orthogonal, and the normal stresses ${\sigma }_{\zeta }$ and ${\sigma }_{\eta }$ are the principal stresses. The principal axes of anisotropy coincide with the x—and y—axes at the initial instant. Then, according to the model proposed in [35], the principal anisotropy axes coincide with the $\zeta -$ and $\eta -$ coordinate curves throughout the process of deformation. Under plane strain conditions, the anisotropic yield criterion proposed in [11] reads

$$\left|{\sigma }_{\eta }-{\sigma }_{\zeta }\right|=2T\sqrt{1-c}$$

(6)

Here T is the shear yield stress with respect to the $\zeta -$ and $\eta -$ coordinate curves, and c is expressed through the yield stresses in the principal anisotropy axes’ directions [11]. Both T and c depend on $\zeta$. It is convenient to represent T as $T=T_0\omega \left(\zeta \right)$ where $T_0$ is constant with the dimensions of stress. Then, Equation (6) becomes

$$\left|{\sigma }_{\eta }-{\sigma }_{\zeta }\right|=T_{0}W\left(\zeta \right)$$

(7)

where

$$W\left(\zeta \right)=2\omega \left(\zeta \right)\sqrt{1-c}.$$

(8)

The transformation equations between the $\left(x,\hspace{0.17em}y\right)-$ and $\left(\zeta ,\hspace{0.17em}\eta \right)-$ coordinate systems in (1) satisfy the flow rule associated with the yield criterion (7).

Let ${\xi }_{\zeta }$ and ${\xi }_{\eta }$ be the principal components of the total strain rate tensor. This tensor is the sum of the elastic and plastic strain rate tensors. Then,

$${\xi }_{\zeta }={\xi }_{\zeta }^e+{\xi }_{\zeta }^p\hspace{1em}\mathrm{and}\hspace{1em}{\xi }_{\eta }={\xi }_{\eta }^e+{\xi }_{\eta }^p.$$

(9)

The superscript e denotes the elastic portion of the strain rate components and the superscript p the plastic portion. The elastic portion is related to the stress components as

$${\dot{\sigma }}_{\zeta }-\dot{\sigma }=2G{\xi }_{\zeta }^e\hspace{1em}\mathrm{and}\hspace{1em}{\dot{\sigma }}_{\eta }-\dot{\sigma }=2G{\xi }_{\eta }^e.$$

(10)

Here the superimposed dot denotes the convected derivative, $\sigma$ is the hydrostatic stress and G is the shear modulus of elasticity. Since the material is incompressible, $\sigma =\left({\sigma }_{\zeta }+{\sigma }_{\eta }\right)/2$. Then, Equation (10) becomes

$${\dot{\sigma }}_{\zeta }-{\dot{\sigma }}_{\eta }=4G{\xi }_{\zeta }^e\hspace{1em}\mathrm{and}\hspace{1em}{\dot{\sigma }}_{\eta }-{\dot{\sigma }}_{\zeta }=4G{\xi }_{\eta }^e$$

(11)

The only stress equilibrium equation which is not identically satisfied in the polar coordinate system is

$$\frac{\partial {\sigma }_r}{\partial r}+\frac{{\sigma }_r-{\sigma }_{\theta }}{r}=0.$$

(12)

It is evident from (4) that ${\sigma }_r={\sigma }_{\zeta }$ and ${\sigma }_{\theta }={\sigma }_{\eta }$. Therefore, Equation (12) becomes

The boundary conditions are

$$\frac{\partial {\sigma }_{\zeta }}{\partial r}+\frac{{\sigma }_{\zeta }-{\sigma }_{\eta }}{r}=0.$$

(13)

$${\sigma }_{\zeta }=0$$

(14)

for $r=r_{AB}$ (or $\zeta =0$) and $r=r_{CD}$ (or $\zeta =-1$).

The bending moment is determined as

$$M={\displaystyle \underset{r_{CD}}{\overset{r_{AB}}{\int }}{\sigma }_{\eta }rdr}.$$

(15)

3. General Solution

Using (1), one can immediately find the total principal strains as

$${\epsilon }_{\zeta }=-\frac{1}{2}\ln \left[4\left(\zeta a+s\right)\right]\hspace{1em}\mathrm{and}\hspace{1em}{\epsilon }_{\eta }=\frac{1}{2}\ln \left[4\left(\zeta a+s\right)\right].$$

(16)

Using (4), one transforms Equation (13) to

$$\frac{\partial {\sigma }_{\zeta }}{\partial \zeta }+\frac{a\left({\sigma }_{\zeta }-{\sigma }_{\eta }\right)}{2\left(\zeta a+s\right)}=0$$

(17)

and Equation (15) to

$$m=\frac{4M}{T_0{H}^2}=\frac{2}{a}{\displaystyle \underset{-1}{\overset{0}{\int }}\frac{{\sigma }_{\eta }}{T_0}d\zeta }.$$

(18)

It is worthy of note that $m=1$ if $W\left(\zeta \right)=1$ in (7), as follows from [11].

3.1. Purely Elastic Solution

In the case of the purely elastic solution, ${\xi }_{\zeta }^e={\xi }_{\zeta }$ and ${\xi }_{\eta }^e={\xi }_{\eta }$ in (11). The solution of Equations (11) and (17) supplemented with the equation of strain compatibility is [36]

$$\frac{{\sigma }_{\zeta }}{T_0}=\frac{1}{2k}{\ln }^2\left[4\left(\zeta a+s\right)\right]+C,\hspace{1em}\frac{{\sigma }_{\eta }}{T_0}=\frac{1}{2k}{\ln }^2\left[4\left(\zeta a+s\right)\right]+\frac{2}{k}\ln \left[4\left(\zeta a+s\right)\right]+C.$$

(19)

Here $k=T_0/G$ and C is constant.

3.2. Solution in Plastic Regions Where ${\sigma }_{\eta }>{\sigma }_{\zeta }$

In this case, the yield criterion (7) becomes

$${\sigma }_{\eta }-{\sigma }_{\zeta }=T_{0}W\left(\zeta \right).$$

(20)

Substituting this equation into (17) leads to

$$\frac{\partial {\sigma }_{\zeta }}{T_0\partial \zeta }-\frac{aW\left(\zeta \right)}{2\left(\zeta a+s\right)}=0.$$

(21)

One integrates this equation to get

$$\frac{{\sigma }_{\zeta }}{T_0}=\frac{a}{2}{\displaystyle \underset{{\zeta }^{\left(1\right)}}{\overset{\zeta }{\int }}\frac{W\left(\mu \right)}{\left(\mu a+s\right)}d\mu }+\frac{{\sigma }_{\zeta }^{\left(1\right)}}{T_0}.$$

(22)

Here $\mu$ is a dummy variable of integration and

$${\sigma }_{\zeta }={\sigma }_{\zeta }^{\left(1\right)}.$$

(23)

at $\zeta ={\zeta }^{\left(1\right)}$. Equations (20) and (22) combine to give

$$\frac{{\sigma }_{\eta }}{T_0}=\frac{a}{2}{\displaystyle \underset{{\zeta }^{\left(1\right)}}{\overset{\zeta }{\int }}\frac{W\left(\mu \right)}{\left(\mu a+s\right)}d\mu }+W\left(\zeta \right)+\frac{{\sigma }_{\zeta }^{\left(1\right)}}{T_0}.$$

(24)

3.3. Solution in Plastic Regions Where ${\sigma }_{\eta }<{\sigma }_{\zeta }$

In this case, the yield criterion (7) becomes

$${\sigma }_{\zeta }-{\sigma }_{\eta }=T_{0}W\left(\zeta \right).$$

(25)

Substituting this equation into (17) leads to

$$\frac{\partial {\sigma }_{\zeta }}{T_0\partial \zeta }+\frac{aW\left(\zeta \right)}{2\left(\zeta a+s\right)}=0.$$

(26)

One integrates this equation to get

$$\frac{{\sigma }_{\zeta }}{T_0}=\frac{a}{2}{\displaystyle \underset{\zeta }{\overset{{\zeta }^{\left(2\right)}}{\int }}\frac{W\left(\mu \right)}{\left(\mu a+s\right)}d\mu }+\frac{{\sigma }_{\zeta }^{\left(2\right)}}{T_0},$$

(27)

where

$${\sigma }_{\zeta }={\sigma }_{\zeta }^{\left(2\right)}$$

(28)

at $\zeta ={\zeta }^{\left(2\right)}$. Equations (25) and (27) combine to give

$$\frac{{\sigma }_{\eta }}{T_0}=\frac{a}{2}{\displaystyle \underset{\zeta }{\overset{{\zeta }^{\left(2\right)}}{\int }}\frac{W\left(\mu \right)}{\left(\mu a+s\right)}d\mu }-W\left(\zeta \right)+\frac{{\sigma }_{\zeta }^{\left(2\right)}}{T_0}.$$

(29)

4. Initiation of Plastic Yielding

The entire sheet is elastic at the beginning of the process. In this case, the solution (19) is valid in the range $-1\le \zeta \le 0$. Therefore, this solution should satisfy the boundary conditions in Equation (14). Then,

$$\frac{1}{2k}{\ln }^2\left(4s\right)+C=0\hspace{1em}\mathrm{and}\hspace{1em}\frac{1}{2k}{\ln }^2\left[4\left(s-a\right)\right]+C=0.$$

(30)

Eliminating C between these equations yields

$$16s\left(s-a\right)=1.$$

(31)

Solving this equation for s and then using any of the equations in (30), one gets

$$s=\frac{2a+\sqrt{4a^2+1}}{4},\hspace{1em}C=-\frac{1}{2k}{\ln }^2\left(2a+\sqrt{4a^2+1}\right).$$

(32)

Substituting (32) into (19) supplies the principal stresses’ through-thickness distribution at any value of a. This solution is valid if the yield criterion is not violated in the range $-1\le \zeta \le 0$. Equations (7) and (19) combine to give

$$-\frac{k}{2}\le \Phi \left(\zeta \right)\le \frac{k}{2}$$

(33)

where

$$\Phi \left(\zeta \right)=\frac{\ln \left[4\left(\zeta a+s\right)\right]}{W\left(\zeta \right)}.$$

(34)

If $\Phi \left(\zeta \right)$ is a monotonically increasing or decreasing function of its argument, then the yield criterion may violate at $\zeta =0$, or $\zeta =-1$, or $\zeta =0$ and $\zeta =-1$ simultaneously. The corresponding conditions are

$$\frac{\ln \left(4s\right)}{W\left(0\right)}=\pm \frac{k}{2},\hspace{1em}\mathrm{or}\hspace{1em}\frac{\ln \left[4\left(s-a\right)\right]}{W\left(-1\right)}=\mp \frac{k}{2},\hspace{1em}\mathrm{or}\hspace{1em}\frac{\ln \left(4s\right)}{W\left(0\right)}=\pm \frac{k}{2}\hspace{1em}\mathrm{and}\hspace{1em}\frac{\ln \left[4\left(s-a\right)\right]}{W\left(-1\right)}=\mp \frac{k}{2}.$$

(35)

Here the upper sign corresponds to monotonically increasing functions, and the lower sign to monotonically decreasing functions.

Let $a_e$ and $s_e$ be the values of a and s, respectively, that correspond to plastic yielding initiation. These values are readily found from (32) and (35) if the function $W\left(\zeta \right)$ is prescribed. In particular, if the distribution of material properties is symmetric relative to the surface $\zeta =-1/2$, then $W\left(0\right)=W\left(-1\right)$ and it follows from the third case in (35) that $16s_e\left(s_e-a_e\right)=1$. The latter equation coincides with (31). Therefore, if the distribution of material properties is symmetric relative to the surface $\zeta =-1/2$, then the initiation of plastic yielding occurs at $\zeta =0$ and $\zeta =-1$ simultaneously.

If the $\Phi \left(\zeta \right)$ has a local minimum or maximum, then one should consider the possibility of the initiation of plastic yielding at such points, in addition to the points $\zeta =0$ and $\zeta =-1$. The initiation of plastic yielding at $\zeta ={\zeta }_m$ where $-1<{\zeta }_m<0$ occurs if

$${\Phi }^{\prime }\left({\zeta }_m\right)=0\hspace{1em}\mathrm{and}\hspace{1em}\left|{\sigma }_{\eta }-{\sigma }_{\zeta }\right|=T_{0}W\left({\zeta }_m\right).$$

(36)

Here ${\Phi }^{\prime }\left(\zeta \right)\equiv d\Phi /d\zeta$. Using (19) and (34), one transforms Equation (36) to

$$aW\left({\zeta }_m^{\left(1\right)}\right)-\left({\zeta }_m^{\left(1\right)}a+s\right)W^{\prime }\left({\zeta }_m^{\left(1\right)}\right)\ln \left[4\left({\zeta }_m^{\left(1\right)}a+s\right)\right]=0\hspace{1em}\mathrm{and}\hspace{1em}\ln \left[4\left({\zeta }_m^{\left(1\right)}a+s\right)\right]=\frac{k}{2}W\left({\zeta }_m^{\left(1\right)}\right)$$

(37)

or

$$aW\left({\zeta }_m^{\left(2\right)}\right)-\left({\zeta }_m^{\left(2\right)}a+s\right)W^{\prime }\left({\zeta }_m^{\left(2\right)}\right)\ln \left[4\left({\zeta }_m^{\left(2\right)}a+s\right)\right]=0\hspace{1em}\mathrm{and}\hspace{1em}\ln \left[4\left({\zeta }_m^{\left(2\right)}a+s\right)\right]=-\frac{k}{2}W\left({\zeta }_m^{\left(2\right)}\right)$$

(38)

Here $W^{\prime }\left(\zeta \right)\equiv dW/d\zeta$. The equations in (37) combine to give

$$a_e^{\left(1\right)}=\frac{k}{8}{W}^{\prime }\left({\zeta }_m^{\left(1\right)}\right)\exp \left[\frac{k}{2}W\left({\zeta }_m^{\left(1\right)}\right)\right].$$

(39)

Then, Equation (32) supplies $s=s_e^{\left(1\right)}$ as a function of ${\zeta }_m^{\left(1\right)}$. One can replace a and s in the second equation in (37) with $a_e^{\left(1\right)}$ and $s_e^{\left(1\right)}$, respectively, to arrive at the equation for ${\zeta }_m^{\left(1\right)}$. This equation should be solved numerically. Then, $a_e^{\left(1\right)}$ is readily determined from (39).

One applies the same procedure to (38) for determining $a_e^{\left(2\right)}$, $s_e^{\left(2\right)}$, and ${\zeta }_e^{\left(2\right)}$. In particular,

$$a_e^{\left(2\right)}=-\frac{k}{8}{W}^{\prime }\left({\zeta }_m^{\left(2\right)}\right)\exp \left[-\frac{k}{2}W\left({\zeta }_m^{\left(2\right)}\right)\right].$$

(40)

If $-1<{\zeta }_e^{\left(1\right)}<0$ and $-1<{\zeta }_e^{\left(2\right)}<0$, then $a_e=\min \left\{a_e^{\left(1\right)},a_e^{\left(2\right)}\right\}$.

In what follows, it is assumed that $a\ge a_e$.

5. Rigid Plastic Solution

For many applications, it is possible to assume that the inelastic behavior is dominant and to neglect the elastic response [37,38]. In this case, the general solution given in Section 3 simplifies. There are two plastic regions throughout the process of deformation. The inequality ${\sigma }_{\zeta }>{\sigma }_{\eta }$ is valid in one of these regions, and the inequality ${\sigma }_{\zeta }<{\sigma }_{\eta }$ in the other. The neutral line separates the plastic regions. The stress ${\sigma }_{\eta }$ is discontinuous across the neutral line.

Let $\zeta ={\zeta }_n$ be the neutral line. It has been shown in [27] that

$${\zeta }_n=-\frac{ds}{da}.$$

(41)

The solution (27) is valid in the region $-1\le \zeta \le {\zeta }_n$. This solution should satisfy the boundary condition (14) at $\zeta =-1$. Then, ${\zeta }^{\left(2\right)}=-1$ and ${\sigma }_{\zeta }^{\left(2\right)}=0$ in (27), and this equation becomes

$$\frac{{\sigma }_{\zeta }}{T_0}=-\frac{a}{2}{\displaystyle \underset{-1}{\overset{\zeta }{\int }}\frac{W\left(\mu \right)}{\left(\mu a+s\right)}d\mu }.$$

(42)

Consequently, Equation (29) becomes

$$\frac{{\sigma }_{\eta }}{T_0}=-\frac{a}{2}{\displaystyle \underset{-1}{\overset{\zeta }{\int }}\frac{W\left(\mu \right)}{\left(\mu a+s\right)}d\mu }-W\left(\zeta \right).$$

(43)

The solution (22) is valid in the region ${\zeta }_n\le \zeta \le 0$. This solution should satisfy the boundary condition (14) at $\zeta =0$. Then, ${\zeta }^{\left(1\right)}=0$ and ${\sigma }_{\zeta }^{\left(1\right)}=0$ in (22), and this equation becomes

$$\frac{{\sigma }_{\zeta }}{T_0}=\frac{a}{2}{\displaystyle \underset{0}{\overset{\zeta }{\int }}\frac{W\left(\mu \right)}{\left(\mu a+s\right)}d\mu }.$$

(44)

Consequently, Equation (24) becomes

$$\frac{{\sigma }_{\eta }}{T_0}=\frac{a}{2}{\displaystyle \underset{0}{\overset{\zeta }{\int }}\frac{W\left(\mu \right)}{\left(\mu a+s\right)}d\mu }+W\left(\zeta \right).$$

(45)

The stress ${\sigma }_{\zeta }$ must be continuous across the neutral line. Then, Equations (41), (42) and (44) combine to give

$${\displaystyle \underset{-ds/da}{\overset{0}{\int }}\frac{W\left(\mu \right)}{\left(\mu a+s\right)}d\mu }={\displaystyle \underset{-1}{\overset{-ds/da}{\int }}\frac{W\left(\mu \right)}{\left(\mu a+s\right)}d\mu }.$$

(46)

Since W is a known function of its argument, this equation is an ordinary differential equation for determining s as a function of a. However, its form is non-standard. To develop a numerical method for solving (46), it is advantageous that one considers the initial instant separately. The distribution of the stress ${\sigma }_{\eta }$ at the initial instant is illustrated in Figure 2. The stress ${\sigma }_{\zeta }$ vanishes everywhere.

Then, it follows from (7) that

$${\sigma }_{\eta }=\{\begin{array}{c}{T}_{0}W\left(\zeta \right) \mathrm{in} \mathrm{the} \mathrm{region} {\zeta }_n\le \zeta \le 0\\ -T_{0}W\left(\zeta \right) \mathrm{in} \mathrm{the} \mathrm{region} -1\le \zeta \le {\zeta }_n.\end{array}$$

(47)

In the case of pure bending, equilibrium demands

$${\displaystyle \underset{-1}{\overset{0}{\int }}{\sigma }_{\eta }d\zeta }=0.$$

(48)

Substituting (47) into (48) and using (41) yields

$${\displaystyle \underset{-ds/da}{\overset{0}{\int }}W\left(\zeta \right)d\zeta }-{\displaystyle \underset{-1}{\overset{-ds/da}{\int }}W\left(\zeta \right)d\zeta }=0.$$

(49)

This equation should be solved for $ds/da$. In general, a numerical method should be used. However, in some cases, an analytic solution is available. For example, if W is constant then $ds/da=1/2$ at the initial instant. If the distribution of material properties is symmetric relative to the central plane of the sheet, then $W_Z\left(Z\right)$ is an even function of Z where $Z=\zeta +1/2$ and $W_Z\left(Z\right)=W\left(\zeta \right)=W\left(Z-1/2\right)$. In this case, Equation (49) becomes

$${\displaystyle \underset{-ds/da+1/2}{\overset{1/2}{\int }}{W}_Z\left(Z\right)dZ}-{\displaystyle \underset{-1/2}{\overset{-ds/da+1/2}{\int }}W\left(\zeta \right)d\zeta }=0.$$

(50)

Integrating gives

$$\Omega \left(1/2\right)-2\Omega \left(-ds/da+1/2\right)+\Omega \left(-1/2\right)=0$$

(51)

where $\Omega \left(Z\right)$ is the anti-derivative of $W_Z\left(Z\right)$. Therefore, $\Omega \left(Z\right)$ is an odd function of Z and $\Omega \left(1/2\right)+\Omega \left(-1/2\right)=0$. Then, it follows from (51) that

$$\frac{ds}{da}=\frac{1}{2}$$

(52)

at the initial instant.

Using (47), the bending moment at the initial instant is determined as

$$M=T_0{\displaystyle \underset{-ds/da}{\overset{0}{\int }}W\left(\zeta \right)\left(\zeta +\frac{ds}{da}\right)d\zeta }+T_0{\displaystyle \underset{-1}{\overset{-ds/da}{\int }}W\left(\zeta \right)\left(\zeta +\frac{ds}{da}\right)d\zeta }.$$

(53)

If the distribution of material properties is symmetric relative to the central plane of the sheet, then Equation (53) becomes

$$M=T_0{\displaystyle \underset{0}{\overset{1/2}{\int }}{W}_Z\left(Z\right)ZdZ}-T_0{\displaystyle \underset{-1/2}{\overset{0}{\int }}{W}_Z\left(Z\right)ZdZ}=2T_0{\displaystyle \underset{0}{\overset{1/2}{\int }}{W}_Z\left(Z\right)ZdZ}.$$

(54)

Here Equation (52) has been taken into account.

Assume that one needs to find the solution of (46) at $a=a_f$. Let the interval $0\le a\le a_f$ be subdivided into an arbitrary number of small segments by the points $a_i$. The values of s and $ds/da$ at $a=a_i$ are denoted as $s_i$ and $s_i^{\prime }$, respectively. It is convenient to choose these points a constant distance $\Delta a$ apart. Let the value of $a_i$, $s_i$, and $s_i^{\prime }$ be known. Then, $a_{i+1}=a_i+\Delta a$. The value of $s_{i+1}$ can be approximated as

$$s_{i+1}=s_i+\frac{\left(s_i^{\prime }+s_{i+1}^{\prime }\right)}{2}\Delta a.$$

(55)

Equation (46) at $a=a_{i+1}$ becomes

$${\displaystyle \underset{-s_{i+1}^{\prime }}{\overset{0}{\int }}\frac{W\left(\mu \right)}{\left(\mu a_{i+1}+s_{i+1}\right)}d\mu }={\displaystyle \underset{-1}{\overset{-s_{i+1}^{\prime }}{\int }}\frac{W\left(\mu \right)}{\left(\mu a_{i+1}+s_{i+1}\right)}d\mu }.$$

(56)

One can eliminate $s_{i+1}$ in (56) using (55). The resulting equation contains one unknown, $s_{i+1}^{\prime }$. This equation should be solved numerically. Then, $s_{i+1}$ is readily determined from (55). Having found $s_{i+1}$ and $s_{i+1}^{\prime }$, one can apply the procedure above to find the solution at $a_{i+1}$. This procedure should be repeated until $a_{i+1}$ becomes equal to $a_f$. It remains to find the input data for the first step. It is evident that $a_0=0$ and $a_1=\Delta a$. It follows from (3) that $s_0=1/4$. The solution of (49) supplies $s_0^{\prime }$. However, if the distribution of material properties is symmetric relative to the central plane of the sheet, then $s_0^{\prime }=1/2$, as follows from (52).

6. Unloading

Consider purely elastic unloading. Variations of the shape are neglected during this stage of the process. At the end of loading, $a=a_l$, $s=s_l$, and $m=m_l$. The corresponding values of the inside and outside radii are denoted as $r_{CD}=R_0$ and $r_{AB}=R_1$, respectively. The distribution of ${\sigma }_{\zeta }$ and ${\sigma }_{\eta }$ at $a=a_l$ is denoted as ${\sigma }_{\zeta }^{\left(l\right)}$ and ${\sigma }_{\eta }^{\left(l\right)}$. All the quantities introduced above can be calculated using the solution given in the previous sections.

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 for the plane strain bending under tension provided in [36] assuming the tensile force vanishes. In our nomenclature, this solution reads

$$\begin{array}{l}\frac{\Delta {\sigma }_{\zeta }}{2G}=V_0\ln \left(\frac{r}{R_0}\right)-U_0\left[1-{\left(\frac{R_0}{r}\right)}^2\right],\hspace{1em}\hspace{1em}\frac{\Delta {\sigma }_{\eta }}{2G}=V_0\left[1+\ln \left(\frac{r}{R_0}\right)\right]+U_0\left[1+{\left(\frac{R_0}{r}\right)}^2\right],\\ U_0=-\frac{{\rho }_0^2\ln {\rho }_0{m}_{l}k}{2\left[{\left({\rho }_0^2-1\right)}^2-4{\rho }_0^2{\ln }^2{\rho }_0\right]}\frac{H^2}{R_0^2},\hspace{1em}{V}_0=-\frac{{\rho }_0^2\left(1-{\rho }_0^2\right)m_{l}k}{2\left[{\left({\rho }_0^2-1\right)}^2-4{\rho }_0^2{\ln }^2{\rho }_0\right]}\frac{H^2}{R_0^2}\end{array}$$

(57)

where ${\rho }_0=R_0/R_1$. Equation (4), in which a should be replaced with $a_l$ and s with $s_l$, supplies the dependence of r on $\zeta$. The distribution of residual stresses is given by

$${\sigma }_{\zeta }^{res}={\sigma }_{\zeta }^{\left(l\right)}+\Delta {\sigma }_{\zeta }\hspace{1em}\mathrm{and}\hspace{1em}{\sigma }_{\eta }^{res}={\sigma }_{\eta }^{\left(l\right)}+\Delta {\sigma }_{\eta }.$$

(58)

After calculating the residual stresses, it is necessary to verify that the yield criterion is not violated in the range $-1\le \zeta \le 0$. The corresponding condition follows from (7) in the form

$$\left|{\sigma }_{\eta }^{res}-{\sigma }_{\zeta }^{res}\right|\le T_{0}W\left(\zeta \right).$$

(59)

7. Practical Example

The through-thickness distribution of the coefficients involved in Hill’s quadratic yield criterion [11] has been experimentally determined in rolled sheets of Al-Mg-Si alloy in [32]. The distribution is symmetric relative the central plane of the sheets. In our nomenclature,

Table 1. Through-thickness distribution of W.

Surface (ζ = 0 and ζ = −1)	ζ = −1/4 and ζ = −3/4	Center (ζ = −1/2)
0.644	0.659	0.6

represents the results from [32]. It is seen from this table that the function W is non-monotonic in each half of the sheet. This discrete function is approximated to the following continuous function:

$$W\left(\zeta \right)=0.644-0.837\zeta -4.888{\zeta }^2-8.103{\zeta }^3-4.051{\zeta }^4.$$

(60)

It is worthy to note that the experimental data were first approximated by an even function of $\zeta +1/2$ and then Equation (60) was derived. The numerical solution has been found using the general solutions described in Section 3, Section 4, Section 5 and Section 6 and (60). In all calculations, $k=0.001$.

It has been found that $a_e\approx 6\cdot {10}^{-5}$. Two plastic regions initiate in the vicinity of the inside and outside surface almost simultaneously and quickly propagate to the corresponding stress-free surface. This stage of the process is very short and is not significant for bending at large strains. The solution with three regions (two plastic regions and an elastic region between them) starts at $a=a_p\approx 2.7\cdot {10}^{-4}$. At $a\ge a_p$, the plastic region adjacent to the surface $\zeta =0$ (plastic region 1) occupies the domain ${\zeta }_1\le \zeta \le 0$ and the plastic region adjacent to the surface $\zeta =-1$ (plastic region 2) occupies the domain $-1\le \zeta \le {\zeta }_2$. The elastic region occupies the domain ${\zeta }_2\le \zeta \le {\zeta }_1$. Thus $\zeta ={\zeta }_1$ and $\zeta ={\zeta }_2$ are the elastic/plastic boundaries. Both ${\zeta }_1$ and ${\zeta }_2$ depend on a.

The distribution of ${\sigma }_{\zeta }$ and ${\sigma }_{\eta }$ in plastic region 1 follows from (22) and (24), and in plastic region 2 from (27) and (29). Using the boundary conditions in (14), one finds that ${\sigma }_{\zeta }^{\left(1\right)}=0$, ${\zeta }^{\left(1\right)}=0$,${\sigma }_{\zeta }^{\left(2\right)}=0$, and ${\zeta }^{\left(2\right)}=-1$. The stresses ${\sigma }_{\zeta }$ and ${\sigma }_{\eta }$ must be continuous across the elastic/plastic boundaries. This requirement is equivalent to the requirement that ${\sigma }_{\zeta }$ and ${\sigma }_{\zeta }-{\sigma }_{\eta }$ are continuous across the elastic/plastic boundaries. Using (19), (20), (22), (24), (25), (27) and (29), one can represent the latter as

$$\begin{array}{c}{\ln }^2\left[4\left({\zeta }_{1}a+s\right)\right]+2kC=ak{\displaystyle \underset{0}{\overset{{\zeta }_1}{\int }}\frac{W\left(\mu \right)}{\left(\mu a+s\right)}d\mu ,}\\ {\ln }^2\left[4\left({\zeta }_{2}a+s\right)\right]+2kC=ak{\displaystyle \underset{{\zeta }_2}{\overset{-1}{\int }}\frac{W\left(\mu \right)}{\left(\mu a+s\right)}d\mu ,}\\ 2\ln \left[4\left({\zeta }_{1}a+s\right)\right]=kW\left({\zeta }_1\right),\hspace{1em}2\ln \left[4\left({\zeta }_{2}a+s\right)\right]=-kW\left({\zeta }_2\right).\end{array}$$

(61)

One can eliminate C between the first two equations. In the resulting equation, $\ln \left[4\left({\zeta }_{1}a+s\right)\right]$ and $\ln \left[4\left({\zeta }_{2}a+s\right)\right]$ can be eliminated using the third and fourth equations. Then,

$${\displaystyle \underset{0}{\overset{{\zeta }_1}{\int }}\frac{W\left(\mu \right)}{\left(\mu a+s\right)}d\mu }-{\displaystyle \underset{{\zeta }_2}{\overset{-1}{\int }}\frac{W\left(\mu \right)}{\left(\mu a+s\right)}d\mu }=\frac{k}{4a}\left[W^2\left({\zeta }_1\right)-W^2\left({\zeta }_2\right)\right].$$

(62)

Moreover, the third and fourth equations in (61) can be solved for a and s. As a result,

$$\begin{array}{c}a=\frac{1}{4\left({\zeta }_1-{\zeta }_2\right)}\left\{\exp \left[\frac{k}{2}W\left({\zeta }_1\right)\right]-\exp \left[-\frac{k}{2}W\left({\zeta }_2\right)\right]\right\},\\ s=\frac{1}{4\left({\zeta }_2-{\zeta }_1\right)}\left\{{\zeta }_2\exp \left[\frac{k}{2}W\left({\zeta }_1\right)\right]-{\zeta }_1\exp \left[-\frac{k}{2}W\left({\zeta }_2\right)\right]\right\}.\end{array}$$

(63)

One can now eliminate a and s in (62) using (63). The resulting equation contains the two unknowns, ${\zeta }_1$ and ${\zeta }_2$. A numerical solution of this equation supplies the dependence of ${\zeta }_1$ on ${\zeta }_2$. Having found this dependence, a and s are readily calculated from (63) giving ${\zeta }_1$, ${\zeta }_2$, and s as functions of a in implicit form. Then, C is determined from the first equation in (61). It has been found that ${\zeta }_1\approx -0.189$ and ${\zeta }_2\approx -0.811$ at $a=a_p$. It is more informative to use $H/r_{CD}$ as a time-like variable instead of a. Equation (5) allows one to replace a with $H/r_{CD}$ with no difficulty. The variation of ${\zeta }_1$ and ${\zeta }_2$ with $H/r_{CD}$ is depicted in Figure 3. It is seen from this figure that the thickness of the elastic region decreases with $H/r_{CD}$ very quickly. In particular, this thickness is less than 0.8% of the sheet’s initial thickness when $H/r_{CD}=0.04$ (e.g., the inside radius of the sheet is 25 times larger than the sheet’s initial thickness). The through-thickness distribution of ${\sigma }_{\zeta }$ and ${\sigma }_{\eta }$ at $H/r_{CD}=0.04$ is shown in Figure 4 and Figure 5, respectively. In these figures, X is the dimensionless distance from the inside surface of the sheet defined as

$$X=\frac{r-r_{CD}}{H}.$$

(64)

Using the numerical procedure described in Section 5, s and $ds/da$ have been found as functions of a. Then, Equation (5) has been used to replace a with $H/r_{CD}$. The dependence of ${\zeta }_n=-ds/da$ on $H/r_{CD}$ in the range $0<H/r_{CD}\le 0.04$ is depicted in Figure 3. It is seen from this figure that the neutral line found from the rigid/plastic solution is located between the two elastic/plastic boundaries found from the elastic/plastic solution.

The dependence of ${\zeta }_n=-ds/da$ on $H/r_{CD}$ in the range $0<H/r_{CD}\le 40$ is shown in Figure 6. It is seen from this figure that the location of the neutral line changes very quickly at the beginning of the process but gradually at large $H/r_{CD}$. The through-thickness distribution of ${\sigma }_{\zeta }$ and ${\sigma }_{\eta }$ at $H/r_{CD}=0.04$ is shown in Figure 4 and Figure 5, respectively. It is seen from Figure 4 that the difference between the elastic/plastic and rigid/plastic solutions is insignificant. The difference is most considerable near the neutral line in the rigid/plastic solution, and the rigid/plastic solution predicts a slightly higher value of $\left|{\sigma }_{\zeta }\right|$ than the elastic/plastic one. The difference between the elastic/plastic and rigid/plastic solutions is invisible in Figure 5. The through-thickness distribution of ${\sigma }_{\zeta }$ and ${\sigma }_{\eta }$ at several stages of the process found by means of the rigid/plastic solution is shown in Figure 7 and Figure 8, respectively. The maximum value of $\left|{\sigma }_{\zeta }\right|$ significantly increases as the deformation proceeds.

The distribution of ${\sigma }_{\eta }$ has a weak local maximum in the vicinity of the outside surface, which may affect bendability. Having found the stress solution, one can calculate the bending moment using (18). The variation of the dimensionless bending moment with $H/r_{CD}$ is depicted in Figure 9. It is seen from this figure that the dimensionless bending moment changes very quickly at the beginning of the process but gradually at large $H/r_{CD}$. Its value attains a local minimum around $H/r_{CD}=13$.

The procedure described in Section 6 has been applied to calculate the distribution of residual stresses. The methodology proposed in [34] has been adopted. Figure 4 and Figure 5 justify the validity of this methodology. These calculations have shown that inequality (59) is not satisfied with relatively large values of $H/r_{CD}$, requiring consideration of reversed yielding. In particular, even if $H/r_{CD}=0.04$, inequality (59) is not satisfied at a narrow region near the neutral line.

8. Conclusions

A general semi-analytic solution for finite pure plane strain bending of a sheet made of incompressible material has been found. A distinguished feature of this solution is that the sheet is plastically anisotropic and inhomogenous. No restriction of the through-thickness distribution of plastic properties is imposed. Both elastic/plastic and rigid/plastic models have been considered. It has been shown that the difference between the elastic/plastic and rigid/plastic solutions is negligible, except at the very beginning of the process. Since the rigid/plastic solution is much simpler than the elastic/plastic one, it is recommended to use rigid/plastic models at large strains. It is possible even if residual stresses should be calculated. In this case, the methodology proposed in [34] can be adopted. The numerical example uses real material properties provided in [32]. The through-thickness distribution of these properties is symmetric relative to the center line of the original sheet. The through-thickness stress distributions are illustrated in Figure 4, Figure 5, Figure 7 and Figure 8. It is seen from Figure 5 and Equation (18) that the difference between the bending moment found using the elastic/plastic and rigid/plastic solutions is negligible if $H/r_{CD}\ge 0.04$, at least.

An advantage of the general solution is that it is valid for an arbitrary distribution of material properties. Therefore, in conjunction with experimental data, this solution can be readily used for identifying these properties.

It is crucial to predict springback in bending followed by unloading accurately [39]. The solution for purely elastic unloading is not valid for all cases considered, and it is necessary to consider the appearance of reversed yielding. The latter will be the subject of a subsequent investigation.