Three-Dimensional Finite Element Analysis for Nonhomogeneous Materials Using Parallel Explicit Algorithm

This chapter addresses the behavior of functionally graded solids under dynamic impact loading within the framework of linear elasticity using parallel explicit algorithm. Numerical examples are presented that verify the dynamic explicit finite element code and demonstrate the dynamic response of graded materials. A three-point bending beam made of epoxy and glass phases under low-velocity impact is studied. Bending stress history for beam with higher values of material properties at the loading edge is consistently higher than that of the homogeneous beam and the beam with lower values of material properties at the loading edge. Larger bending stresses for the former beam may indicate earlier crack initiation times, which were proven by experiments performed by other researchers. Wave propagation in a 3D bar is also investigated. Poisson ’ s ratio and thickness effects are observed in the dynamic behavior of the bar. Finite element modeling and simulation discussed herein can be a critical tool to help understand physics behind the dynamic events.


Introduction
Functionally graded materials (FGMs) are materials characterized by smooth variation in composition, microstructure, and properties. These materials have emerged with the need to enhance material performance and to meet specific functions and applications [1,2]. The material gradation concept has been utilized in various applications [3][4][5][6][7][8][9][10][11][12][13][14][15][16]. The concept of FGMs has been applied to thermal barrier structures and wear-and corrosion-resistant coatings and also used for joining dissimilar materials [17]. FGMs are subjected to harsh thermal/mechanical/dynamic environments. Thin-walled plates and shells, which are used in reactor vessels, turbines, and other machine parts, are susceptible to buckling failure, large deflections, or excessive stresses induced by thermomechanical loading. Functionally graded coatings on these structural elements may help reduce these failures [18].
The dynamic response of FGMs has been investigated [19]. Reddy and Chin [20] carried out nonlinear transient thermoelastic analysis of FGMs by using plate elements for moderately large rotations. Gong et al. [21] studied the elastic response of Dynamic response is evaluated at discrete instants of time separated by time increments Δt. At the nth time step, the equation of motion for linear dynamic problems is given by [34] M ½ € u f g n þ C ½ _ u f g n þ K ½ u f g n ¼ F ext f g n (6) Discretization in time is accomplished using finite difference approximations of time derivatives. Explicit algorithms use a difference expression of the general form: u f g nþ1 ¼ f u f g n ; _ u f g n ; € u f g n ; u f g nÀ1 ; … À Á which contains only historical information on its right-hand side. This difference expression is combined with the equation of motion at time step n.
The stability of the explicit finite element scheme is governed by the Courant condition [34], which provides a critical time step size beyond which the results become unstable. It is given by where l e is the shortest distance between two nodes in the mesh and the dilatational wave speed C d is expressed as Since C d varies due to material gradation, the desired time step for the finite element implementation is chosen based on the largest dilatational wave speed and is maintained constant throughout the simulation. The mass matrix used in the explicit method described above is typically lumped in order to avoid matrix inversion. The HRZ technique [35] is used to lump the mass matrix in which only the diagonal terms are retained and scaled so as to preserve the total mass: where m l is the lumped element mass matrix, m is the consistent mass matrix, and n is the number of degrees of freedom in the element.

Explicit parallel FEA using MPI
A master-slave approach in the parallel finite element code is used here to solve wave propagation-type problems using the Message Passing Interface (MPI) standard [36]. Figure 1 shows the procedure used in the parallel execution of the Time-integration procedure in the parallel execution of the explicit Newmark-β method [37] using MPI standard [36]. explicit Newmark-β method. The first step is to partition the structure (FE mesh) into number of sub-domains. The stiffness matrices are calculated for elements in the local processor. The stiffness matrix needs not be assembled in an explicit analysis as the algorithm is implemented on an element basis only. The parallel explicit time integration analysis of finite elements from onetime step to another consists in obtaining the unknown velocity and displacement of master and slaves, sending the effective force vector from slaves to the master for assembly and sending back the assembled force vector to the slaves for calculation of acceleration vector. The communication between the slaves and the master processor is done using MPI-based function calls as shown in the figure. The entire force vector at each processor is being sent to the master processor for assembly, and the master processor returns the assembled force vector to the slaves in the same manner. These steps are repeated for the rest of the time integration. Local field quantities like element stresses and strains can be calculated separately on each processor without any communication.

Example 1: 3D graded beam under velocity impact
A 3D three-point bending beam under velocity impact is studied to understand the influence of material gradation on the bending behavior in a 3D space. The beam is a real FGM system made of glass/epoxy phases. The dynamic fracture experiments of the linearly graded specimens have been conducted by Rousseau and Tippur [37][38][39][40].
In this paper, we investigate the 3D dynamic behavior of a beam. Using 3D explicit finite element code, we can simulate the 3D dynamic behavior of the beam. This 3D finite element study offers enhanced knowledge of the dynamic behavior of this material system and understanding of the stress field in homogeneous and graded materials, which help to predict fracture initiation times in various graded specimens.
Consider a three-point bending specimen under velocity impact load at the top as shown in Figure 2(a). Due to the symmetry of the geometry and the loading conditions, one-fourth of the beam is modeled as shown in Figure 2(b). Figure 2(c) and 2(d) shows the 3D FE mesh of the one-fourth model and a zoom of the FE mesh at the left bottom corner, respectively. The mesh is refined along the loaded edge with a uniform element size of 92.5 μm, and the stress history at point P(0, 0.2 W) is retrieved, where W is the beam depth. Point P is of significance because it corresponds to the location of the crack tip in the dynamic fracture analysis of the cracked beam [38,41].
Three material gradation cases are considered for the dynamic analysis of the three-point bending beam: a homogeneous beam (Homog, E 2 = E 1 ), a beam stiffer at the impacted surface (StiffTop, E 2 > E 1 ), and a beam softer at the impacted surface (StiffBot, E 2 < E 1 ), where subscripts 1 and 2 denote bottom and top surfaces of the beam, respectively. For the graded beam, Young's modulus (E) is varied linearly between 4 GPa and 12 GPa, and the mass density (ρ) is varied linearly between 1000 kg/cm 3 and 2000 kg/cm 3 . These material properties correspond to values used in [42]. For the homogeneous beam, the mass density (ρ) is taken as the average of the graded beam counterparts, and Young's modulus (E) is calculated such that the equivalent E/ρ value equals that of the graded specimen in an average sense, i.e., Poisson's ratio is 0.33 for all the three cases. The average dilatational wave speed (C d(avg) ) is defined as where C d is calculated assuming plane stress behavior given by The dilatational wave speed (C d ) is calculated assuming plane stress behavior to compare results with those obtained by Zhang and Paulino [42]. The difference in C d(avg) for each material gradation is marginal with 2421.5 m/s for the homogeneous beam and 2418.4 m/s for the graded beams. Hence, the dilatational wave speed is assumed constant as 2421.5 m/s, and the results are normalized based on this average dilatational wave speed. Upon impact on the top surface of the beam, compressive stress waves are generated and propagate toward the bottom surface and reflect back as tensile waves from the bottom surface.
A detailed stress history analysis is performed to obtain the influence of material gradation on the bending behavior of the beam. Figure 3 shows the bending stress (σ x ) history of point P due to impact velocity of 5 m/s for the three cases mentioned above with normalized time t 0 = t*C d /W as the x-axis and stress (MPa) as the y-axis. We see that material gradation leads to considerable change in the stress field at point P. The stress wave takes normalized time of about 0.7 to reach the point P, which is different from 0.8 in Zhang and Paulino [42]. This is due to 3D effects. In a 3D medium, stress waves travel at different speeds when compared to a 2D medium. At this time, the point P undergoes compressive stress until the normalized time of around 1.1. This compressive stress is due to the Poisson's ratio (ν) effect. In a separate simulation of the beam under similar loading conditions, the Poisson's ratio of the beam was made zero, and we found that the point P underwent tensile stress throughout the simulation (result not shown). From this instant, i.e., after normalized time of 1.1, the point P undergoes monotonically increasing tensile stress. The maximum tensile stress is consistently attained in the StiffTop beam, and the minimum in the StiffBot beam and the Homog beam undergoes stresses in between these two. This type of behavior may be intuitively observed as the material properties (E, ρ) of a StiffTop beam at point P are lower than the StiffBot and Homog beams and hence may undergo higher stresses. We also found that for the line load, the stress component σ x along the thickness (z) direction at point P(y = 0.2 W) does not vary considerably. Figure 4 shows the stress history of component σ y at point P for different material gradations due to impact velocity of 5 m/s. We see that there is a considerable effect of the gradation on this stress behavior, too. The first stress wave approaches point P at normalized time of 0.7 due to the compressive velocity loading. At this point, there is a sudden increase in σ y stress magnitude (compressive) until a normalized time of 1.25. The stress wave reaches the bottom surface of the beam and reflects as tensile waves at which point there is a gradual decrease in the compressive stress until the normalized time of 2.0. Between normalized times 2.0 and 2.7, the stress behavior does not change considerably, although there are small variations, which may be because of the discrete wave reflections from the edges that occur in finite configurations, to form a plateau in the stress field. The stress wave, after reflecting from the top surface, becomes compressive and leads to an increase in compressive stress between normalized times 2.7 and 3.2. Between the normalized times of 3.2 and 6, the stress behavior tends to repeat itself as seen between 0.7 and 3.2, although with larger variations of increase and decrease. The results plotted in Figures 13 and 14 are consistent with the plane stress simulation results obtained by Zhang and Paulino [42] and Rousseau and Tippur [38]. Figure 5 shows the 2D contour plot of the stress component σ x at a specific time instant of 90 μs (t 0 $ 5.9) for the three beams subjected to the impact velocity of 5 m/s. We firstly see that the stress pattern is considerably different among the beams considered. Further, at regions close to the loading, compressive stress is concentrated for all the beams, although the size of the compressive zone differs considerably. The StiffTop beam has the largest compressive region at the loading region followed by the Homog beam and the StiffBot beam. At regions close to the loading for the StiffTop beam, we see that the material properties are relatively higher than the other cases and hence the compressive load gets transferred to nearby regions. For the StiffBot case, the compressive region is localized in a small region because of the compliance at the top. Also, we see that the compressive region extends the most in the x direction in the same order as mentioned above. We observe that the size of the tensile zone near the bottom of the beam (i.e., x = 0, y = 0) varies in the same order as the compressive zone at the loading region. The neutral axis location varies for each of the graded beams considered which may be the reason behind the larger tensile region being developed for StiffTop beam when compared to the other two cases. For the dynamic loading considered, the neutral axis shifts in time; nevertheless, the overall stress pattern will be similar to those shown in Figure 5. With respect to crack initiation in such graded beams (assuming that the crack lies parallel to the applied loading), we know that the σ x stress component dominates the most. We see that the StiffTop beam is consistently subjected to larger tensile stresses than other two cases at point P and will lead to earlier crack initiation. Figures 6 and 7 show the σ x and σ y plots for the different graded beams subjected to 5 m/s and 10 m/s velocities of impact. We see that the stress behavior increases linearly (gets doubled almost) for V = 10 m/s when compared to the case of V = 5 m/s for all the gradation cases.
We have so far investigated the 2D dynamic behavior of a beam under the line load using 3D finite element method. Using 3D explicit finite element code, we can simulate the 3D dynamic behavior of the beam. We apply the impact point load of V = 10 m/s at a central node at the top edge of the beam. Note that this cannot be simulated using 2D finite elements in Zhang and Paulino [41]. The stress history results are obtained at point P and compared against those shown in Figure 7.  Figures 8 and 9 show the comparison of σ x and σ y , respectively, for the three beams under impact velocity of V = 10 m/s applied throughout the thickness (line load) and at one central node (point load). We see that the stresses for the latter case (point load) are much lower than the former case (line load). Maximum tensile stress (σ x ) at point P due to point load is experienced by the StiffTop beam followed by the Homog beam and the StiffBot beam, similar to the line load case. The stress values at point P across the z (thickness) direction do not change noticeably (result not shown).

Example 2: wave propagation in 3D bar
Wave propagation in a fixed-free 3D bar with homogeneous and graded materials in the vertical direction (y) is simulated to illustrate the relevance of material gradation and also to verify the parallel explicit 3D finite element code. Consider a  fixed-free bar in Figure 10(a). The square bar is of length L = 1 m and height H = 0.05 m. A step loading is used (Figure 10(b)) at the free end of the bar. The material considered is steel for the homogeneous bar and steel and alumina for the graded bar. The properties of steel and alumina are given in Table 1. For the graded bar, the lower surface of the bar is made of alumina, and the upper part is made of steel, and the material properties vary linearly from alumina to steel.  Explicit time integration method is used for analyzing wave propagation type of problems since it avoids storage of large matrices and does not require the solution of global systems of equations, but is limited by its conditional stability [43]. The time step size has to be smaller than a critical value which is directly dependent on the largest frequency of the FE discretization (smallest element). Higher-order modes participate in the bar response which in turn requires refined element size. To capture the transient response, the time step size needs to be small enough in order to also track the change of applied loads, and adequate element size can be estimated according to the Courant condition [34]. Considering material gradation in the y direction, the FE mesh is discretized into 300 Â 15 Â 15 quads, each quad divided into four 15-node wedge elements, totaling 270,000 elements and 725,836 nodes. The 1D analytical solution for the motion of a bar with step force loading at the free end is obtained using a finite sine transform method given by [44] where P 0 , E, A, L, and c represent the applied force magnitude, the Young's modulus, cross-sectional area, length, and dilatational wave speed of the bar, respectively. The first term on the right side of Eq. (14) is the static solution, whereas the series represents a superposition of normal modes of amplitudes inversely proportional to (2n À 1) 2 . The axial stress can be easily obtained by differentiating displacement Eq. (14).
Results are presented for the 3D wave propagation in a bar under step loading considering homogeneous and graded material properties. Figures 11 and 12 show  the normalized longitudinal stress at three locations along the x-axis obtained using the exact solution and the finite element method, respectively. For the exact solution, we see small oscillation of the stress solution at x = 0 at t = 1 which is due to finite series of n used in Eq. (14). Poisson's ratio (ν) is taken as zero for the present finite element simulation. From the figures, we see that the overall stress behavior of the finite element simulation is similar to the exact solution, although there are  some oscillations. These oscillations are seen at times when there is a sudden increase or decrease in stress due to wave reaching a certain location in the bar. This oscillation effect may be due to discrete wave reflections that occur in finite configurations and also due to the finite element discretization. Figures 13 and 14 show the normalized longitudinal stress at six locations (see the insert) for the functionally graded bar along z = 0.025 m and z = 0 m,  respectively. The stress history at x = L is not plotted as the behavior is similar to the homogeneous case. We see that the stress wave gets distorted in time because of the material gradation. The stress wave at x = 0.5 L and mostly at x = 0 differs highly across the thickness direction due to difference in wave speeds. As expected, the alumina side undergoes higher stresses than the steel side, more so at the fixed end than at other x locations.

Conclusions
Dynamic and wave propagation behavior of 3D functionally graded continua is investigated using explicit finite element formulations. Three numerical examples are presented. The first example is on dynamic analysis of a three-point bending beam made of epoxy and glass phases under low-velocity impact. Material gradation considerably affects the dynamic stress behavior of the beam. Tensile stress is maximum for StiffTop beam at the imaginary crack-tip location for stress component (σ x ) indicating that crack initiation will occur earlier for this type of beam which was also verified by experiments conducted by Rousseau and Tippur [40].
The second example is on 3D wave propagation behavior of a homogeneous and functionally graded bar under sinusoidal and step loading. Stress results for the homogeneous bar with zero Poisson's ratio obtained using the finite element method match very well with the exact solution. Poisson's ratio and material constraints in a 3D continuum affect the stress behavior for both homogeneous and functionally graded bars. In the functionally graded bar, the alumina side consistently experienced more stress than the steel side. Finite element modeling and simulation discussed herein can be a critical tool to help understand physics behind the dynamic events.