A Theory of Elastic/Plastic Plane Strain Pure Bending of FGM Sheets at Large Strain

An efficient analytical/numerical method has been developed and programmed to predict the distribution of residual stresses and springback in plane strain pure bending of functionally graded sheets at large strain, followed by unloading. The solution is facilitated by using a Lagrangian coordinate system. The study is concentrated on a power law through thickness distribution of material properties. However, the general method can be used in conjunction with any other through thickness distributions assuming that plastic yielding initiates at one of the surfaces of the sheet. Effects of material properties on the distribution of residual stresses are investigated.


Introduction
Structures made of functionally graded materials (FGM) are advantageous for many applications. A difficulty with theoretical analysis and design is that structures made of FGM are classified by a much greater number of parameters than similar structures made of homogeneous materials. For this reason, it is desirable to perform parametric studies by analytic or semi-analytic methods as much as possible. A review of results related to the analysis of FGM and published before 2007 is presented in [1]. This review focuses on structures with through-thickness variation of material properties. Analytic solutions derived in [1][2][3][4][5] belong to this class of FGM as well. In [2][3][4], elastic and elastic/plastic spherical vessels subjected to various loading conditions are considered. Thermo-elastic simply supported and clamped circular plates are studied in [5]. Many analytic and semi-analytic solutions are available for FGM discs and cylinders assuming that material properties vary in the radial direction but are independent of the circumferential and axial directions. Purely elastic solutions for a hollow disc or cylinder subjected to internal or/and external pressure are derived in [6][7][8]. An axisymmetric thermo-elastic solution for a hollow cylinder subjected quite a general system of thermo-mechanical loading is presented in [9]. It is assumed that the temperature varies along the radial coordinate. A plane strain analytic elastic/plastic solution for pressurized tubes is found in [10]. The solution is based on the Tresca yield criterion. Many solutions are proposed for functionally graded solid and hollow rotating discs. Purely elastic solutions for solid discs of constant thickness are given in [11,12], a purely elastic solution for a hollow disc of variable thickness in [13], a purely elastic solution for hollow polar orthotropic discs in [14], and a solution for hollow cylinders using the theory of electrothermoelasticity in [15]. An elastic perfectly plastic stress solution for hollow discs is derived in [16] using the von Mises yield criterion.

Basic Equations
The process of plane strain pure bending is illustrated in Figure 1. The approach proposed in [18] for solving the corresponding boundary value problem is based on the following transformation equations: where (x, y) is an Eulerian-Cartesian coordinate system and (ζ, η) is a Lagrangian coordinate system. Without loss of generality, it is possible to assume that the origin of the Cartesian coordinate system is located at the intersection of the axis of symmetry of the process and the outer surface AB and that the x-axis coincides with the axis of symmetry. The Lagrangian coordinate system is chosen such that ζ = x/H and η = y/H at the initial instant where H is the initial thickness of the sheet. It is evident from these relations and the geometry in Figure 1 that ζ = 0 on AB, ζ = −1 on CD, η = L/H on CB and η = −L/H on AD throughout the process of deformation. Here, L is the initial width of the sheet. In Equation (1), a is a time-like variable. In particular, a = 0 at the initial instant. In Equation (1), s is a function of a. This function should be found from the stress solution and therefore depends on constitutive equations. The condition in Equation (2) is satisfied if at a = 0. It is possible to verify by inspection that the mapping in Equation (1) satisfies the equation of incompressibility. Moreover, this mapping transforms initially straight lines A 1 B 1 and C 1 D 1 into circular arcs AB and CD and initially straight lines C 1 B 1 and A 1 D 1 into circular arcs CB and AD after any amount of deformation ( Figure 1). Furthermore, coordinate curves of the Lagrangian coordinate system coincide with trajectories of the principal strain rates and, for coaxial models, with trajectories of the principal stresses. Thus, the shear stress vanishes in the Lagrangian coordinates. In particular, the contour ABCD is free of shear stresses. Let σ ζ and σ η be the physical stress components referred to the Lagrangian coordinates. The stress solution should satisfy the boundary conditions for ζ = −1 and ζ = 0. The only non-trivial equilibrium equation in the Lagrangian coordinates has been derived in [18] as The initial plane strain yield criterion of the functionally graded sheet is supposed to be where σ 0 is a material constant and β(x/H) is an arbitrary function of its argument. It is assumed that material properties are not affected by plastic deformation. Therefore, Equation (6) can be rewritten in the form In this case, the yield locus is invariant along the motion. The importance of this property of material models has been emphasized in [25]. Let τ ζ and τ η be the deviatoric portions of σ ζ and σ η , respectively. Since the material is incompressible, τ ζ + τ η = 0 under plane strain conditions. Then, the yield criterion in Equation (7) is equivalent to Hooke's law generalized on functionally graded materials reads It has been taken into account here that Poisson's ratio is equal to 1/2 for incompressible materials. In addition, ε e ζ and ε e η are the total strain components in elastic regions and the elastic portions of the total strain components in plastic regions referred to the Lagrangian coordinate system, G 0 is a material constant and g(ζ) is an arbitrary function of its argument. Geometric parameters shown in Figure 1 depend on a and are expressed as [18] Once s has been found as a function of a, these parameters are immediate from Equation (10).

Stress Solution at Loading
It is assumed that the functions β(ζ) and g(ζ) involved in Equations (7) and (9) are such that plastic yielding can only initiate at ζ = 0 or ζ = −1. This assumption can be verified using the purely elastic solution with no difficulty. At the very beginning of the process, the entire sheet is elastic. As deformation proceeds, one of the following three cases arises: (i) plastic yielding initiates at the surface ζ = −1; (ii) plastic yielding initiates at the surface ζ = 0; and (iii) plastic yielding initiates simultaneously at the surfaces ζ = −1 and ζ = 0. These cases should be treated separately. In the following, ζ 1 is the elastic/plastic boundary between the plastic region that propagates from the surface ζ = 0 and the elastic region and ζ 2 is the elastic/plastic boundary between the plastic region that propagates from the surface ζ = −1 and the elastic region. It is evident that both ζ 1 and ζ 2 depend on a. The general structure of the solution with two plastic regions is illustrated in Figure 2. Let M be the bending moment. Then, its dimensionless representation is in terms of the Lagrangian coordinates given by [18]  In the elastic region, the whole strain is elastic. Therefore, it follows from Equation (1) that the principal logarithmic strains are Since σ ζ − σ η = τ ζ − τ η , Equations (5) and (9) combine to give ∂σ ζ ∂ζ Eliminating the strain components in Equation (13) Integrating this equation with respect to ζ and using the boundary condition in Equation (4) at where k = σ 0 /(3G 0 ) and χ is a dummy variable of integration. The expression for σ η in Equation (15) has been derived using the identity σ η = σ ζ − τ ζ + τ η , and Equations (9) and (12). In the case of the purely elastic solution, Equation (15) must satisfy the boundary condition in Equation (4) at ζ = −1. Then, the equation for the function s(a) is Using Equation (15), in which s should be eliminated by means of the solution of Equation (16), and the yield criterion in Equation (8), it is possible to determine which of the three cases mentioned above occurs for given material properties. Simultaneously, the value of a at which plastic yielding initiates is determined. This value of a is denoted as a e . In the following, it is assumed that a ≥ a e . It is now necessary to consider Cases (i), (ii) and (iii) separately.
Case (i). There are two regions. A plastic region occupies the domain −1 ≤ ζ ≤ ζ 2 and an elastic region the domain ζ 2 ≤ ζ ≤ 0. Equation (15) is valid in the elastic region. However, the function s(a) is not determined from Equation (16). It is reasonable to assume that σ η < σ ζ in the plastic region. Therefore, the yield criterion in Equation (7) becomes Substituting Equation (17) into Equation (5) and integrating yields the dependence of the stress σ ζ on ζ. Using Equation (17) again provides the dependence of the stress σ η on ζ. As a result, It is evident that this solution satisfies the boundary condition in Equation (4) at ζ = −1. Both σ ζ and σ η should be continuous across ζ = ζ 2 . Consequently, τ ζ is continuous across ζ = ζ 2 . The stress τ ζ on the elastic side of the elastic/plastic boundary is determined from Equation (15) and on the plastic side from Equation (8). Then, the condition of continuity of τ ζ across the surface ζ = ζ 2 is represented as Solving this equation for s yields Using Equations (15) and (18), the condition of continuity of σ ζ across the surface ζ = ζ 2 is represented as In this equation, s can be eliminated by means of Equation (20). The resulting equation should be solved numerically to find ζ 2 as a function of a. Then, s as a function of a is readily found from Equation (20). The yield criterion should be checked in the elastic region using the solution in Equation (15). The calculation should be stopped when the yield condition is satisfied at one point of the elastic region. Denote the corresponding value of a as a 2 .
In Case (i), Equation (11) becomes In the first integrand, σ η /σ 0 should be eliminated by means of Equation (18) and in the second by means of Equation (15).
Case (ii). There are two regions. A plastic region occupies the domain ζ 1 ≤ ζ ≤ 0 and an elastic region the domain −1 ≤ ζ ≤ ζ 1 . The elastic solution in Equation (15) satisfies the boundary condition in Equation (4) at ζ = 0. Therefore, it is convenient to rewrite this solution as The elastic solution in this form satisfies the boundary condition in Equation (4) at ζ = −1. It is reasonable to assume that σ η > σ ζ in the plastic region. Therefore, the yield criterion in Equation (7) becomes Substituting Equation (24) into Equation (5) and integrating yields the dependence of the stress σ ζ on ζ. Using Equation (24) again provides the dependence of the stress σ η on ζ. As a result, It is evident that this solution satisfies the boundary condition in Equation (4) at ζ = 0. Both σ ζ and σ η should be continuous across ζ = ζ 1 . Consequently, τ ζ is continuous across ζ = ζ 1 . The stress τ ζ on the elastic side of the elastic/plastic boundary is determined from Equation (23) and on the plastic side from Equation (8). Then, the condition of continuity of τ ζ across the surface ζ = ζ 1 is represented as Solving this equation for s yields Using Equations (23) and (25), the condition of continuity of σ ζ across the surface ζ = ζ 1 is represented as In this equation, s can be eliminated by means of Equation (27). The resulting equation should be solved numerically to find ζ 1 as a function of a. Then, s as a function of a is readily found from Equation (27). The yield criterion should be checked in the elastic region using the solution in Equation (23). The calculation should be stopped when the yield condition is satisfied at one point of the elastic region. Denote the corresponding value of a as a 1 .
In Case (ii), Equation (11) becomes In the first integrand, σ η /σ 0 should be eliminated by means of Equation (23) and in the second by means of Equation (25). Case (iii). In this case, there are two plastic regions, −1 ≤ ζ ≤ ζ 2 and ζ 1 ≤ ζ ≤ 0, and one elastic region, ζ 1 ≤ ζ ≤ ζ 2 . At the beginning of this stage of the process, a = a 1 and ζ 2 = −1 or a = a 2 and ζ 1 = 0. Let σ n1 be the value of σ ζ at ζ = ζ 1 and σ n2 be the value of σ ζ at ζ = ζ 2 . Then, the elastic solution in Equation (15) can be rewritten as It follows from this solution that The solution in Equation (18) is valid in the plastic region −1 ≤ ζ ≤ ζ 2 and the solution in Equation (25) in the plastic region ζ 1 ≤ ζ ≤ 0 . Then, and Equations (20) and (27) are valid. Therefore, and Equations (31)-(33) combine to give Eliminating in this equation s by means of Equation (20) or Equation (27) and then a by means of Equation (35) supplies the equation to find ζ 1 as a function of ζ 2 (or ζ 2 as a function of ζ 1 ). Then, a as a function of ζ 1 (or ζ 2 ) is found from Equation (35) and s as a function of ζ 1 (or ζ 2 ) from Equation (20) or (27). The distribution of the stresses is determined from Equation (30) with the use of Equations (32) and (33) in the elastic region, from Equation (18) in the region −1 ≤ ζ ≤ ζ 2 and from Equation (25) in the region ζ 1 ≤ ζ ≤ 0.
In Case (iii), Equation (11) becomes In the first integrand, σ η /σ 0 should be eliminated by means of Equation (18), in the second by means of Equation (30) and the third by means of Equation (25). As usual, it is necessary to verify that the yield criterion is not violated in the elastic region.

Unloading
It is assumed that unloading is purely elastic. This assumption should be verified a posteriori. At this stage of the process, the strains can be considered as infinitesimal. Let a f and s f be the values of a and s, respectively, at the end of loading. These values are known from the solution given in the previous section. Using Equation (10), the values of r CD and r AB at the end of loading, r f CD and r f AB , are determined as It is convenient to introduce a polar coordinate system (r, θ) with the origin at x = −H √ s f /a f and y = 0 (point O 1 in Figure 1). The coordinate curves of this coordinate system coincide with the coordinate curves of the (ζ, η)-coordinate system. Therefore, σ ζ = σ r and σ η = σ θ where σ r and σ θ are the normal stresses in the polar coordinate system. Moreover, r = R f H at ζ = 0 and r = r f H at ζ = −1. The equilibrium equation for the increment of the stresses, ∆σ ζ and ∆σ η , in the polar coordinate system can be written as where ρ = r/H. Since σ ζ = 0 at ζ = 0 and ζ = −1 at any stage of the process, the increment of this stress should satisfy the conditions for ζ = 0 and ζ = −1.
The displacement components from the configuration corresponding to the end of loading in the polar coordinate system are supposed to be where U 0 and V 0 are dimensionless constants. Using Equation (41), the increment of the normal strains in the polar coordinate system is determined as The increment of the deviatoric stresses is found from Equation (42) and the Hooke's law (Equation (9)) where the stresses and strains should be replaced with the corresponding increments. Then, Using this solution, the right hand side of Equation (36) can be rewritten as The Lagrangian coordinate ζ at the end of loading is expressed in terms of ρ as [18] It is evident that this solution satisfies the boundary condition in Equation (40) at ζ = −1 (or ρ = r f ). The other boundary conditions in Equations (40) and (46) combine to yield Solving this equation for V 0 results in Using Equations (43) and (46), it is possible to represent the distribution of ∆σ η as The constant V 0 can be eliminated in Equations (46) and (49) by means of Equation (48). It is then obvious that both ∆σ ζ and ∆σ η are proportional to U 0 . The distribution of the residual stresses follows from Equations (46) and (49) in the form σ res (50) As before, ζ should be eliminated by means of Equation (45) and V 0 by means of Equation (48). The constant U 0 remains to be found. To this end, it is necessary to use the condition that the bending moment vanishes at the end of unloading. Using Equations (11) and (45), this condition can be represented as This equation should be solved for U 0 numerically. Then, Equation (50) supplies the distribution of the residual stresses. To verify that the solution given in Section 4 is valid, this distribution should be substituted into the yield criterion in Equation (7) where σ ζ and σ η should be replaced with σ res ζ and σ res η , respectively. The left-hand side of Equation (7) should be less than or equal to 2/ √ 3 σ 0 β(ζ) in the range −1 ≤ ζ ≤ 0 .

Numerical Examples
Several numerical examples are presented in this section, based on the analytical solutions developed in the previous sections. Our chosen modulus gradient function is g(ζ) = 1 + (G 1 /G 0 − 1)(−ζ) N , and yield stress gradient function β(ζ) = 1 + (σ 1 /σ 0 − 1)(−ζ) N . The power law exponent N controls the functional distribution of material properties along the thickness coordinate ζ. The power law distributions in modulus and yield stress with the same N have been proposed in the literature [26,27]. The material parameters used in our numerical calculations are listed in Table 1.

Homogeneous Sheet under Bending
When the homogenous sheet is under bending, both edges will simultaneously develop plastic zones. Figure 3a shows the movement of the two elastic-plastic boundaries toward the centerline of the sheet, as deformation magnitude increases. The deformation magnitude is measured by parameter a. The applied bending moment is a function of a, as shown in Figure 3b. As an illustration of the developed analytical solutions in previous sections, Figure 4 shows the stress distributions along the sheet under two different deformation magnitudes. As can be seen, the plastic zones increase with a for σ η , while σ ζ remains in elastic regime. The reason for σ η is not perfectly horizontal in the plastic zone is because our numerical codes do not allow N set equal to zero, hence a very small N is chosen, as shown in Table 1. After unloading, Figure 5a,b shows the residual stress distributions under two different as. Larger a increases the magnitude of residual stresses after unloading. Moreover, the residual stress σ ζ is zero at the left and right edges, as indicated by the red short-dashed line.

FGM Sheet Belonged to Case (i) under Bending
In Case (i), the left edge (ζ = −1) of the sheet has smaller yield stress, hence a plastic zone will start on the left edge first. Figure 6 shows the relationship between the applied bending moment and deformation magnitude a with the gradient function exponent N = 1 and 3. As can be seen, larger m is required for N = 3 than that for N = 1, as a increases. Under given deformation magnitudes, Figure 7 shows the stress distributions in the Case (i) FGM under plastic deformation. Larger plastic zone is developed at the left edge as deformation increases. After unloading, residual stress distributions are shown in Figures 8 and 9 for the N = 1 and N = 3 FGM, respectively. Residual stresses are more predominant at the left edge. In addition, the residual stress σ ζ is zero at the left and right edges, as indicated by the red short-dashed line.

FGM Sheet Belonged to Case (ii) under Bending
In Case (ii), the right edge (ζ = 0) of the sheet has smaller yield stress, hence a plastic zone will start on the right edge first. Figure 10 shows the relationship between the applied bending moment and deformation magnitude a with the gradient function exponent N = 1 and 3. As can be seen, larger m is required for N = 1 than that for N = 3, as a increases. Under given deformation magnitudes, Figure 11 shows the stress distributions in the Case (ii) FGM under plastic deformation. Larger plastic zone is developed at the right edge as deformation increases. After unloading, residual stress distributions are shown in Figures 12 and 13 for the N = 1 and N = 3 FGM, respectively. Residual stresses are more predominant at the right edge. The residual stress σ ζ is zero at the left and right edges, as indicated by the red short-dashed line. The results from the illustrative examples solved here may serve as benchmark solutions for data obtained from numerical or experimental methods.

Conclusions
Efficient analytical and numerical methods and procedures have been developed and programmed to predict the distribution of stresses in a sheet of incompressible material subject to plane strain pure bending at large strain and then the distribution of residual stresses after unloading. Springback is also predicted. It has been assumed that the sheet is made of functionally graded material. The general theory has been developed for an arbitrary through thickness distribution of material properties assuming that the initiation of plastic yielding occurs at one of the surfaces of the sheet. This assumption can be verified using the purely elastic solution (Equation (15)) and the yield criterion (Equation (7)). It is possible to use the general solutions (Equations (15), (18) and (25)) even if the assumption is not satisfied but constants of integration should be added. Then, these general solutions should be combined to satisfy the boundary conditions and the conditions at elastic/plastic boundaries. An illustrative example is concentrated on power law distributions of material properties. Using the numerical code developed in this work enables the effect of parameters involved in these laws to be predicted effectively. The calculated examples show the analytical solutions derived here can systematically treat the plastic problems of the homogenous or functionally graded sheet. The magnitudes of applied moment may be strongly influenced by the power law exponent as deformation increases, which provides an effective way to design the functionally graded sheets.
The method employed to derive the solution in this paper can be extended to cyclic loading. This new solution may be useful for the interpretation of experimental data from the reverse bending fatigue test (for example, [28]).