A Novel Approach to Elasto-Plastic and Plastic Finite Element Analysis of Functionally Graded Beams under Transverse Pressure Loading

Article history: Received August 15, 2017 Revised October 2, 2017 Accepted October 20, 2017 In this study a new numerical-analytical method for elasto-plastic and plastic modelling and simulating of FG beams is presented. The functionally graded (FG) beam composed of ten layers through the thickness and it is assumed that the mechanical properties of the beam vary through the thickness described by a simple power law distribution in terms of the volume fractions of constituents. The beam is assumed to be under transverse pressure load. In this paper a new method is presented based on linearization of the nonlinear part of the stress-strain curve of the material of the layers of the FG beam and using the elastic relations for bending analysis of beams. Numerical results for functionally graded beam are given and results of this paper for homogeneous beam are compared with other methods and good agreement is obtained between them. In addition, the effects of material properties on the stress field through the thickness of the FG beam are determined and discussed. Keyword: Elastic–plastic stress analysis Functionally graded beams FGM Simple power law Transverse pressure loads


Introduction
Nowadays, beams are widely used in different industries especially in civil structures. They can be made of different materials depending on their application. Meanwhile, the complex case regarding the ingredients, are beams made of functionally graded materials. In these beams, mechanical properties are varied along its thickness as a function of thickness of the beam. Functionally materials were invented by Japanese scientists as advances composite materials [1]. The use of these materials is being developed in various industries.
Less work has been done in elasto-plastic analysis of beams and plates made of advanced materials. Akbulut et al. [2] investigated elasto-plastic stress analysis of composite plates with metal context under plane loads. Sayman et al. [3] performed elasto-plastic stress analysis of a one-side clamped composite beam reinforced with thermoplastic fibers under a uniform distributed load. Sayman et al. [4] proposed an analytical method for elasto-plastic stress analysis of one-side clamped thermoplastic composite beams under bending moment. Štok and Halilovicˇ EE studied elastoplastic bending of beams with rectangular cross section area. The functional solutions derived under the assumption of at most quadratic bending moment distribution enable fully analytical tracing of the elasto-plastic state evolution in structural components by monotonic and proportional application of loads to a beam structure. Wen and Zeng [5] presented an incremental secant stiffness formulation for materially non-linear analysis of planar beam structures under monotonically increasing external loads. To describe the elasto-plastic behaviour of a typical beam member, a set of non-dimensional plasticity coefficients are introduced to progressively deteriorate the elastic stiffness properties over an incremental load history. Singha et al. [6] investigated the nonlinear behaviors of functionally graded material (FGM) plates under transverse distributed load using a high precision plate bending finite element. Material properties of the plate were assumed to be graded in the thickness direction according to a simple power-law distribution in terms of volume fractions of the constituents. This paper studies a new numerical-analytic method for elasto-plastic and plastic modeling of FG beams under transverse pressure loads. This method is based on linearization of plastic zone of the beam's material and applying analytical relationships of elastic case. In other words, nonlinear region in the stress-strain diagram of the material is converted into small linear regions and required parameters are calculated using correlations in linear analysis and method of width changing of beam layers. The FE-based software ABAQUS is implemented for simulations. Axial stress values in different layers of FG beam and displacement of its center are presented.

Properties of FGMs
A FG beam of length a and thickness h, composed of a mixture of ceramics and metals, as shown in figure 1, is considered. It is assumed the top surface of the beam (z = h/2) is ceramic, whereas the bottom surface (z = -h/2) is metal-rich. The mechanical properties of this beam such as the Young's modulus (E), the shear modulus (G) and the poison's ratio are functions of depth (z). In other words, these mechanical properties are considered to be varied functionally through the thickness of FG beam. In this study the simple power law distribution function is used for description of mechanical properties. The functional relationship between P and z for ceramic and metal FG beams is given by [1], denotes a typical material property (E,G,υ ). 1 P and 2 P denote values of the variables at top and bottom surfaces of beam, respectively. n is a variable parameter. The working range of n can be taken as 3 1 -3, as any magnitude outside this range will provide an FGM having too much of a phase.

Governing equations for elastic bending of beams
In this section, the fundamental equations of the elastic bending and deflection analysis of pressure-loaded functionally graded beams are briefly outlined. Position of the neutral line is determined after determining the equal cross area of FG beam, By using ( ) n z the above integrals can be easily calculated and the position of neutral line is obtained. Crosssectional moment of inertia is determined as, Where, c is the distance from neutral line. Radius of curvature and displacement of the simply supported beam with transverse uniform distributed load are obtained by [8], Where, R is the neutral line curvature radius and q is the transverse uniform distributed load on the beam. To investigate the plastic case, passing of elasto-plastic boundary at loading stage must be studied; i.e. the material begins to yielding and changes its condition from elastic to plastic. For this purpose, a criterion should be used for determining the yield onset of materials.

Yield criterion
In this paper, Von-mises yield criterion is used. Initial yielding occurs in the material when [9] Y o x = σ (9)

Equation of equal cross area of beam
As it is said, the mechanical properties of FG beam vary continuously through the thickness. In these patterns, it is necessary to determine the equal cross area of beam based on the Young's modulus of various layers of beam. The equal cross area of beam is calculated based on the Young's modulus of the bottom layer of beam and equation of equal cross area as, ( )

Elasto-plastic analysis procedure
In this section the mathematical code is written using MATHEMATICA and it is presented as a flowchart. This code does the elasto-plastic analysis using the elastic bending relation of beams and the stress-strain curvature.
In this study, the nonlinear part of the stress-strain curvature is converted to small linear increments, as shown in figure 3. In composite and FG beams, it is necessary to determine the equal cross area of beam based on the smallest Young's modulus of elasticity of the layers of beam using the Eqs. (3),(4), firstly. Then the uniformly transverse load is exerted step by step on beam and in each step, the stress of each layer in the middle of beam and deflection are calculated based on the Eqs. (4), (5), respectively. In each step the mathematical code compare the calculated stress and the yield stress of each layer of beam. If the calculated stress of a given layer is smaller than its yield stress, the code continues the procedure and goes the later step. But if the calculated stress of a given layer of beam is equal or lager than its yield stress, the code uses the small linear increment of the stress-strain curve of that layer and using the slip of that increment, the equal cross area of beam recalculated and other calculation is done based on the new cross area. In other words, this code based on the magnitude of load and the magnitude of stress of each layer of beam at its middle, recalculates and improves the equal cross area and calculates stresses and deflections based on the elastic relations. This procedure is regularly done step by step by increasing of load. Results are depicted in the end of each step but the latest results due to complete magnitude of loading are significant. Flowchart of this mathematical code is completely presented in figure 4.

Properties of materials of FG beam
As previously showed, the analysis of FG beam is applied for type of ceramic and metal combination. The set of materials considered is alumina and aluminum. In all cases, the lower layer of the beam is assumed to be pure metal Calculate the axial strain of each layer of beam using the elastic relations.
Is the calculated stress for each layer larger than the yield stress of it?
Put the Young's modulus of elasticity equal with the calculated amount from simple power law distribution for each layer.
Put the Young's modulus of elasticity equal with (Calculated stress / calculated strain) for each layer.
Recalculate equal cross area parameters (neutral line position, moment of inertia and the magnitude of area of beam) using obtained result until here.
Calculate the axial stresses in each layer of beam using elastic relations. Is the calculated stress for each layer larger than the yield stress of it?
Is TEMP equal with 1? Temp=1 Is the calculated stress for each layer larger than the yield stress of it?
Put the Young's modulus of elasticity equal with the calculated amount from simple power law distribution for each layer.
Put the Young's modulus of elasticity equal with (Calculated stress / calculated strain) for each layer.
Increase the magnitude of load in the amount of one step-tome and put the stress of each layer of beam in new load equal with calculated stress until here.
Is the magnitude of load equal with the final load? Let TEMP be zero.
Calculate the center deflection of FG beam using elastic relations and print the amount of center deflection and axial stresses for each layer in the middle of the beam.

Yes No
No Yes (aluminum) and the upper layer is assumed to be pure ceramic (alumina). The properties of alumina and aluminum are represented in table 1. In elasto-plastic and plastic analysis, other mechanical properties such as the yield and ultimate stresses of the materials of FG beam are required. In table 2 these properties of alumina and aluminum are indicated. In this section some assumption are taken as • Compressive yield stress and compressive ultimate stress of aluminum are equal with tensile yield and ultimate stresses. • Tensile ultimate stress of alumina is equal with tensile yield stress of it.

Yes
• Compressive ultimate stress of alumina is equal with compressive yield stress of it.
Diagram of stress-strain for aluminum used in this section is as following:  Figure 6 shows the elasto-plastic displacement diagram at the center of FG beam with power law distribution function for n 0.5 = bases on the applied load. As seen in these figures, the results of the ABAQUS simulation and numerical-analytical method are approximately corresponding to each other until the deformations are linear, but by onset of elasto-plastic deformations they are separated. This is due to many factors, including how to enter the layer properties in software, the number of time intervals in the software and the rate of load increase in each interval. Also, it can be observed that the rate of displacement of the beam center greatly increases after the plastic zone onset. = bases on the applied load. The behavior of this diagram is same as previous one, however, due to varying yield stress and the different number of layers entering the plastic zone, displacement values and slope of diagrams will be different for both cases. Furthermore, the results of ABAQUS and the presented numerical-analytical method provide an acceptable agreement with each other. The elasto-plastic displacement diagram at the center of FG beam with power law distribution function for n 2 = bases on the applied load is represented in figure 8. In this case, 13 layers of the beam thickness are entered into the plastic zone, but the slope of displacement-time diagram is not changed. This is due to the properties variation in a way that effect of lower layers of neutral line in beam stiffness is less compared to upper layers.

Figure 8. Elasto-plastic displacement of FG beam center for n 2 =
The elasto-plastic displacement diagram at the center of FG beam with power law distribution function for n = ∞ bases on the applied load is shown in figure 9. It is evident that the results of ABAQUS and the presented numericalanalytical method provide an acceptable agreement with each other. The elasto-plastic displacement diagram at the center of FG beam with power law distribution function for n = ∞ bases on the applied load is shown in figure 10 using both analytical and numerical-analytical methods. It shows the validity of results of the analytical and numerical-analytical method compared by ones available in the literature [9].

Axial stress diagrams
Unlike the diagrams provided for large and small deformations which all were plotted at an applied load, in this section diagrams of each case are presented for specific load. Also, due to the lack of plastic zone in the stress-strain diagram of aluminum, the axial stress diagram is not considered for the beam with power law distribution for 0 = n . Figure 11 shows the elasto-plastic axial stress diagram at the center of FG beam with power law distribution function for n 0.

Conclusion
In this paper a new numerical-analytical method is presented for simulation and analysis of beam made of functionally graded materials at elasto-plastic and plastic cases. It is assumed that the mechanical properties of the beam vary through the thickness described by a simple power law distribution in terms of the volume fractions of constituents. The displacement rate of FG beam with power law distribution function is reduced for n 0.5, 1, 2 = by increasing n. By increasing n, the ceramic layers decreases, and the average beam modulus of elasticity decreases and causes displacement increase of beam center at non-plastic deformations. But in this case the reason of displacement rate reduction is due to fact that in elasto-plastic analysis in addition to variation of elasticity modulus, variation of layers yield stress are major factors affecting the behavior of the beam. It is observed that by increasing n the beam layers are yielded at higher loads and this is also due to the variation of yield stress in the beam layers. The number of yielded layers for a specific load (-1.8 MN/m) is reduced by increasing n, and consequently the beam displacement is reduced. The first layer plastic deformation is from internal layers of beam for n 0 = , n 1 = cases, whereas for other cases the onset of plastic zone occurs at the outer layers of beam like homogenous beam. This is due to three main reasons: variation of elasticity modulus, yield stress of layers and distance of each layer from neutral line. The obtained results for a specific case (homogenous beam) are in good agreement with analytical results.