Calculation of stress in FGM beams

Content of the paper is oriented to calculation of elastic normal stress in the Functionally Graded Beams (FGM). Spatial variation of material properties is considered in the lateral, transversal and longitudinal direction of the straight beam. The displacements and internal forces are calculated using our new FGM finite beam element. Heterogeneous material properties are homogenized by extended mixture rules, laminate theory and reference volume element (RVE). Obtained results by our approach are evaluated and compared with the ones obtained by the 3D solid finite elements.


Introduction
Many papers deal with static analysis of the FGM single 2D beams with transverse variation of material properties only [1][2][3][4][5].Less attention is paid to longitudinal and lateral variation of material properties.However, the authors did not find papers which consequently deal with all the longitudinal, transversal and lateral variation of material properties by single beams or beam structures built of such FGM.
The presented contribution is the continuation of our previous work dealing with the derivation of a 3D FGM beam finite element with longitudinal varying effective material properties, which is suitable for analysis of beam structures made of spatially varying FGM.Effect of the axial and shear forces and Winkler elastic foundation is included, as well.Homogenization of the spatial varying material properties in the real FGM beam and the calculation of effective parameters are done by the multilayer or direct integration method (MLM or DIM) [6].If only transversal and lateral variations of material properties are considered in the real FGM beam, longitudinally constant effective material properties arise from the homogenization.This method can also be used in the homogenization of multilayer beams with discontinuous variation of material properties in transversal and lateral direction.Numerical experiments are performed to calculate the elastostatic response of chosen spatial FGM beam structures with rectangular cross-sections with symmetrically lateral, transversal and longitudinal variations of material properties.Expressions for normal stress calculation in the cross-section caused by the axial force as well by the bending moments are established.For torsional loading only the twist angle and internal torsional moment is evaluated.The solution results are discussed and compared to those obtained by means of very fine 3D -solid finite element mesh.

3D FGM beam finite element with spatially varying material properties
Let us consider a 3D straight finite beam element (Timoshenko beam theory and Saint-Venant torsion theory) of doubly symmetric cross-section -Fig. 1.The nodal degrees of freedom at node i are: the displacements u i , v i , w i in the local axis direction x, y, z and the cross-sectional area rotations - M .The first derivative with respect to x of the relevant variable is denoted with an apostrophe "´".
are the elastic foundation modules (the torsional elastic foundation is not considered).The effective homogenized and longitudinally varying stiffness reads: is the effective elasticity modulus for axial loading), is the effective elasticity modulus for bending about axis y), is the flexural stiffness about the axis z, (  is the effective elasticity modulus for bending about axis z), is the effective shear modulus and sm y k is the average shear correction factor in y -direction), is the effective shear modulus and sm y k is the average shear correction factor in z -direction), is the effective torsional stiffness  is the torsional elasticity modulus and T I is the torsion constant).
Mixture rules are one of the methods for micromechanical modeling of heterogeneous materials.Extended mixture rules [7] are based on the assumption that the constituents volume fractions (formally only denoted here as fibres -f and matrix -m) continuously vary as polynomial The condition has to be fulfilled.The appropriated material property distribution in the real FGM beam (Fig. 2a.) then reads Here, are the spatial distributions of material properties of the FGM constituents.The extended mixture rule (1) can be analogically used for FGM material made of more than two constituents.The assumption of a polynomial variation of the constituent's volume fractions and material properties enables an easier establishing of the main appropriated field equations and allows the modeling of many common realizable variations.
In literature and in practical applications, mostly the one directional variation of the FGM properties is considered.There, an exponential law for transversal variation of the constituents volume fractions is often presented, e.g. in [8,9,10] and in references therein.
For the FGM beams the transversal variation (continuously or discontinuously, symmetrically or asymmetrically) has been mainly considered.The homogenization of such material properties variation is relatively simple.If the material properties vary only with respect to the longitudinal direction, the homogenization is frequently not needed since there are new FGM beam and link finite elements established that consider such variations in a very accurate and effective way, e.g. in [6,11].A more complicated case is, if the material properties vary in three directions -namely in transversal, lateral and longitudinal direction of the FGM real beam and the torsion is included as well.
In this contribution, the homogenization techniques for spatially varying (continuously or discontinuously and symmetrically in transversal and lateral direction, and continuously in longitudinal direction) material properties of FGM beams of selected doubly-symmetric cross-sections are considered.The expressions are proposed for the derivation of effective elasticity modules for axial loading, the transversal and lateral bending, the shear modules for transversal and lateral shear and for uniform torsion by the extended mixture rules (EMR) and the multilayer method (MLM).
Let us consider a two nodded 3D straight beam element with double symmetric crosssectional area A (Fig. 2).The composite material of this beam arises from mixing two components.The continuous polynomial spatial variation of the elasticity moduli and mass density can be caused by continuous polynomial spatial variation of both the volume ) and material properties of the FGM constituents ( In our case the elasticity modulus E(x,y,z), the Poisson ratio for the real beam have been calculated by expression (1).The FGM shear modulus can be calculated by expression: If the constituents Poisson's ratio are approximately of the same value and the constituent volume fractions variation is not strong, then the FGM shear modulus can by calculated using a simplification [12]: where  is an average value of the function 

Fig. 2. FGM beam with rectangular cross-section
Homogenization of the spatially varying material properties (the reference volume is the volume of the whole beam) are done in two steps.In the first step, the real beam (Fig. 2a) is transformed into a multilayer beam (Fig. 2b).Material properties of the layers are calculated with the EMR [13].We assume that each layer cross-section at position x has constant material properties.They are calculated as an average value from their values at the boundaries of the respective layer cross-section.Polynomial variation of these properties appears in the longitudinal direction of the layer.Sufficient accuracy of the proposed substitution of the continuous transversal and lateral variation of material properties by the layer-wise constant distribution of material properties is reached if the division to layers is fine enough.In the second step, the effective longitudinal material properties of the homogenized beam are derived using the MLM [14].These homogenized material properties are constant through the beam's height and depth but they vary continuously along the longitudinal beam axis.Accordingly, the beam finite element equation is established for the homogenized beam (Fig. 2c) in order to calculate the primary beam unknowns -the displacements and internal forces and critical buckling force in our case.The stress has to be calculated on the real beam and it is the main goal of this paper that will be presented in the next chapter.

Normal stress calculation in the FGM rectangular crosssectional area
Stress calculation depends of the cross-sectional area form.In the next we focus on the calculation of normal stress distribution on rectangular cross section.Calculation of the shear and torsional stress in the FGM beams is more complicated in our approach.We deal with this task very active and the results will be published in our future publication.The rectangular cross-section at position x with doubly symmetric variation of material properties is considered -Fig.3.At the centroid acts the normal force   are the quadratic moments of the rectangular cross-section.

Fig. 3. To normal stress calculation
The resultant normal stress at point (x,y,z) is the sum: There, is the normal stress caused by normal force and is the effective normal stress in homogenized cross-section while constant effective normal strain over the whole cross-sectional area is assumed.Here

 
x E NH L is the effective elasticity modulus for tension-compression [7].
If the effective bending stiffness about the axis z is , then the distribution of the effective bending strain is is the effective elasticity modulus for bending about z-axis [7].
Analogically, if the effective bending stiffness about the axis y is


, then the distribution of the effective bending strain is and the effective bending stress is is the effective elasticity modulus for bending about z-axis [7].It is not simply obtaining functional distribution of the normal stress.For discontinuous variation of the material properties in transversal and lateral direction, the stress has to be calculated at discrete points of crosssection.

Numerical experiments
In the chapter, the results of numerical investigations concerning the normal stress calculation in the FGM beams with spatially varying material properties are presented.

Normal stress calculation in the FGM cantilever beam
The cantilever FGM beam is considered (as shown in Fig. 4), which is loaded at its free end by forces 10    The TiC volume fraction varies in the y and z direction linearly and symmetrically according to the x-y and x-z planes: -the core of the beam is made from pure Al6061-TO and linearly varies to the edges that are made from pure TiC.Constant effective material properties are considered in the local x -direction of the beam.162.6 MPa have been calculated [15].
The FGM beam, clamped at the node i, has been studied by the elastic-static and buckling analysis.All the calculations were done with our 3D FGM beam finite element (NFE) which we have implemented into the code MATHEMATICA [16].Additionally, the effect of axial force was considered.It has to be pointed out that the entire structure is discretized using only one our finite element.The critical buckling force is 28 10.
In the elasto-static analysis the axial force (tension and compression) N N II  have been chosen as a part of the critical buckling force II Ki N .The same problem has been solved using a very fine mesh -28889 of SOLID186 elements of the FEM program ANSYS [17].The results of ANSYS as well as the results of the NFE are presented in Table 1.The average relative difference  [%] between displacements calculated by our method and the ANSYS solution has been evaluated.As is shown in Table 1, the accuracy of our results is excellent compared to the ANSYS solution.The bending moments about the y and z -axis, the torsional moment, the axial, transversal and shear forces are calculated.To save space, the distribution of the bending moment about y-axis is shown in Fig. 5 only.As expected, the axial force x F influence distribution of the bending moment significantly.The high efficiency of our method (NFE) is obvious since our results are evaluated using only one finite beam element compared to the large number of 41924 elements used in the continuum mesh (ANSYS).

Displacements
According the expressions in chapter 2, the normal stress caused by axial force are calculated.Because of discontinuous variation of elasticity modulus in real beam, the cross section is discretized on layers with constant elasticity modulus -Fig.6.

Fig. 6. Compression normal stress distribution caused by the axial force
The maximal axial force stress are in the layer number 10 and minimal in layer 1 -see Table 2. Accuracy of the calculation approach depends of the number of the layers.As shown in Table 2, for then layers is the difference between our results (NFE) and the ones obtained by very fine mesh of 3D-ssolid FE (ANSYS) very small.By very fine division to layers a continuous distribution of the stress can be obtained.Distribution of the bending normal stress is shown in Fig. 7 and Table 3.

Normal stress calculation in the FGM cantilever beam on elastic foundation
The cantilever FGM beam on varying Winkler foundation is considered (as shown in Fig. 8).Its rectangular cross-section is constant with height h = 0.005 m and width b = 0.01 m.The length of the beam is L = 0.1 m.

Fig. 8. Cantilever beam on elastic foundation with spatially varying material properties
The beam is made of a mixture of two components: Aluminum Al6061-TO and Titanium Carbide TiC, their constant constituent's material properties are given in Case I.The Aluminum volume fraction, in this case, varies linearly and symmetrically according to the x-y and x-z planes: At node i is GPa.The FGM cantilever beam has been studied by elastic-static analysis.All the calculations were done with our 3D FGM beam finite element (NFE).Additionally, the effect of axial force was considered.It has to be pointed out that the entire structure is discretized using only one herein proposed finite element.The cantilever FGM beam resting on varying vertical Winkler elastic foundation are used [12].The displacements at node j are evaluated using the only one new FGM beam finite element (NFE).The same problem is solved using a very fine mesh -23015 of SOLID186 elements of the FEM program ANSYS [18 34].The results of ANSYS as well as the results of the NFE are presented in Table 4.The average relative difference  [%] between displacements calculated by our method and the ANSYS solution is evaluated.The comparison of the vertical beam deflection curve with and without elastic foundation is shown in Fig. 10.

Fig. 10. Vertical beam deflection curve with and without elastic foundation
The bending moments about the y and z -axis for case without elastic foundation are shown in Fig. 11 and Fig. 12, respectively.The Fig. 13 and Fig. 14 show the transversal force in y and z -axis.The comparison of the bending moments for the case F x = -2 kN calculated by our approach and by ANSYS are compared in Table 5.In Table 6 the comparison of the resultant normal stress, caused by axial, transversal and lateral forces, in the clamped cross-section (x = 0) calculated by our approach (Fig. 15) and by ANSYS is shown.[MPa] at position x = 0 calculated by our approach In Fig. 16, the resultant normal stress, caused by axial, transversal and lateral forces (marked in blue flags), in the clamped cross-section is shown that was calculated by ANSYS (using a very fine mesh -23015 of SOLID186).As can be seen, a very good agreement of both solution method has been obtained at all marked points excluding the corners.As is well known, the solutions with 3D solid finite elements produce in the sharp corners incorrect stress first of all by the very fine meshes.In the nearby points of the sharp corners is the match of the results very good (see the details in Fig. 16).

Conclusions
In this contribution, our new 3D beam finite element for elasto-static and buckling analysis of the FGM single beam and beam structures has been used for calculation of the deformation and internal forces and moments in FGM beams with spatial variation of material properties.Effect of varying Winkler elastic foundations and shear force deformation effect is taken into account.The effect of axial force has been taken into account for the flexural loading.For an elastic-static analysis, the expressions for

Fig. 1 .
Fig. 1.The local internal variables and loads rectangular cross-section is constant with height h = 0.005 m and width b = 0.01 m.The length of the beam is L = 0.1 m.The local coordinates system is denoted by the axis x, y, z.

Fig. 4 .
Fig. 4. Cantilever FGM beam with spatial variation of material properties Material of the beam consists of two components: Aluminum Al6061-TO -denoted with index m and Titanium carbide TiC -denoted with index f.The material properties of the components are assumed to be constant and their values are: Aluminum Al6061-TOthe elasticity modulus 0 69.E m  GPa, the Poisson's ratio 33 0. m   ; Titanium carbide TiC -the elasticity modulus 0 480.E f  GPa, the Poisson's ratio 20 0. f   .

Fig. 7 . 3 .
Fig. 7. Resultant biaxial bending normal stress calculated by ANSYS Table 3. Comparing of the biaxial bending normal stress at chosen position in clamped cross-section Position  [kPa] NFE

Table 2 .
Axial force normal stress at the clamped cross-section

Table 4 .
Displacements at node j with and without elastic foundation