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

Abstract

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.

1. 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 FGM shells subjected to low-velocity impact. Rousseau and Tippur [22, 23, 24] studied the dynamic fracture of a three-point bending beam made of epoxy/glass phases under low-velocity impact. Cheng et al. [25] used a peridynamic model to investigate dynamic fracture of FGMs. Lindholm and Doshi [26] looked into a slender bar with free ends and obtained a solution that is synthesized from eigenfunctions by using the principle of virtual work. Karlsson et al. [27] developed a Green’s function approach for 1D transient wave propagation in composite materials. Santare et al. [28] used graded finite elements to simulate elastic wave propagation in graded materials. Chakraborty and Gopalakrishnan [18] developed a spectral element to study wave propagation behavior in FGM beams subjected to high-frequency impact loads.

The development of parallel computers and improvements in computer hardware has been of significant importance [29]. A dynamic event has a key characteristic that the incident pulse duration is very small in the order of microseconds that makes the frequency content of pulse very high (in the order of kHz). Thus, all the higher-order modes participate in the dynamic response, and the finite element mesh must be very fine enough to capture the small wavelengths [18]. This makes the system size enormously large. But in an explicit analysis, storage of large matrices is avoided with the added advantage of not requiring to solve linear algebraic equations [29]. Explicit methods are also conditionally stable, and so the time step size must be below a critical value that is dependent on the size of the smallest finite element [30]. Due to mesh size constraint, extremely small time steps (0.01–100 μs) are needed to solve explicit problems. This can be overcome by parallel computing. Krysl and Belytschko [30] investigated to parallelize an existing object-oriented code written in C where a Parallel Virtual Machine (PVM) was used for the necessary communication. Krysl and Bittnar [31] and Rao [32, 33] adopted the Message Passing Interface (MPI) standard for the implementation of the exchange algorithm, which is used in this chapter.

This chapter addresses the behavior of functionally graded materials under dynamic impact loading within the framework of linear elasticity using MPI-based parallel explicit algorithm. Numerical examples are presented that verify the dynamic explicit finite element code and demonstrate the dynamic response of graded materials. Finite element modeling and simulation can be a critical tool to help understand physics behind the dynamic events.

2. Finite element formulation for dynamic analysis

The governing equation for structural dynamics is given by [34]

$${\int }_{\mathrm{\Omega }}\left(\mathrm{\sigma }:\mathrm{\delta E}+\mathrm{\rho }\ddot{\mathrm{u}}\bullet \mathrm{\delta u}\right)\mathrm{d\Omega }-{\int }_{{\Gamma }_{\mathit{ext}}}{\mathrm{T}}_{\mathrm{ext}}\bullet \mathrm{\delta u}\mathrm{d\Gamma }=0$$E1

where Ω is domain volume, Γ is the boundary with a normal vector n, u is the displacement vector, T is the traction vector, σ is the Cauchy stress tensor, $\ddot{\mathrm{u}}$ is the acceleration vector, ρ is the mass density, ${\Gamma }_{\mathit{ext}}$ represents the boundary surface on which external traction ${\mathrm{T}}_{\mathrm{ext}}$ is applied, and E is the Green strain tensor. For a finite element of volume V and surface area S, the work balance becomes

$$\begin{array}{l}\int {\left\{\text{δu}\right\}}^{\text{T}}\left\{\text{F}\right\}\text{dV}+\int {\left\{\text{δu}\right\}}^{\text{T}}\{\varnothing \}\text{dS}+\sum _{\text{i}=1}^{\text{n}}{\left\{\text{δu}\right\}}_{\text{i}}^{\text{T}}{\left\{\text{p}\right\}}_{\text{i}}\hfill \\ =\int \left({\left\{\delta u\right\}}^T\rho \left\{\ddot{u}\right\}+{\left\{\delta u\right\}}^{T}c\left\{\dot{u}\right\}+ {\left\{\delta \epsilon \right\}}^T\left\{\sigma \right\}\right)dV\hfill \end{array}$$E2

where {F} and {Φ} represent prescribed body forces and surface tractions, ${\left\{p\right\}}_i$ and ${\left\{\mathit{\delta u}\right\}}_i\left({\left\{\mathit{\delta \epsilon }\right\}}_i\right)$ represent prescribed concentrated loads and virtual displacements and strains at a total of n points, and c is a damping coefficient. The finite element discretization provides

$$\left\{u\right\}=\left[N\right]\left\{d\right\};\left\{\dot{u}\right\}=\left[N\right]\left\{\dot{d}\right\};\left\{\ddot{u}\right\}=\left[N\right]\left\{\ddot{d}\right\};\left\{\epsilon \right\}=\left[B\right]\left\{d\right\}$$E3

in which N are the shape functions, the nodal displacement vectors $\left\{d\right\}$ are functions of time, and [B] is the matrix of shape function derivatives. The equation of motion is

$$\left[M\right]\left\{\ddot{D}\right\}+\left[C\right]\left\{\dot{D}\right\}+\left[K\right]\left\{D\right\}=\left\{F\right\}$$E4

where the consistent mass (m), damping (c), and stiffness (k) matrices are given by

$$\left[M\right]=\int \varrho \left(\mathit{x}\right){\left[N\right]}^T\left[N\right]\mathit{dV}\left[C\right]=\int c\left(\mathit{x}\right){\left[N\right]}^T\left[N\right]\mathit{dV}\left[K\right]=\int {\left[B\right]}^T\left[D\left(\mathit{x}\right)\right]\left[B\right]\mathit{dV}$$E5

where mass density, damping coefficient, and constitutive matrices are functions of spatial positions.

3. Explicit method, stability, and mass lumping

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]

$$\left[M\right]{\left\{\ddot{u}\right\}}_n+\left[C\right]{\left\{\dot{u}\right\}}_n+\left[K\right]{\left\{u\right\}}_n={\left\{F^{\mathit{ext}}\right\}}_n$$E6

Discretization in time is accomplished using finite difference approximations of time derivatives. Explicit algorithms use a difference expression of the general form:

$${\left\{\mathrm{u}\right\}}_{\mathrm{n}+1}=\mathrm{f}\left({\left\{\mathrm{u}\right\}}_{\mathrm{n}},{\left\{\dot{\mathrm{u}}\right\}}_{\mathrm{n}},{\left\{\ddot{\mathrm{u}}\right\}}_{\mathrm{n}},{\left\{\mathrm{u}\right\}}_{\mathrm{n}-1},\dots \right)$$E7

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

$$∆t\le l_e/C_d$$E8

where $l_e$ is the shortest distance between two nodes in the mesh and the dilatational wave speed C_d is expressed as

$$C_d^2=\frac{E\left(x\right)\left(1-\nu \left(x\right)\right)}{\left(1+\nu \left(x\right)\right)\left(1-2\nu \left(x\right)\right)\rho \left(x\right)}$$E9

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:

$$m_{\mathit{ii}}^l=m_{\mathit{ii}}\left(\sum _{j=1}^n{m}_{\mathit{ij}}\right)/\left(\sum _{i=1}^n{m}_{\mathit{ii}}\right)$$E10

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.

4. 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 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.

Figure 1.

Time-integration procedure in the parallel execution of the explicit Newmark-β method [37] using MPI standard [36].

5. 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].

Figure 2.

Epoxy/glass beam subjected to velocity impact: (a) geometry and boundary conditions of the three-point bending specimen and (b) one-quarter model with symmetric boundary conditions. Stress values are retrieved at P(0, 0.2 W). (c) 3D FE mesh of the one-quarter model. The FE mesh contains 14,085 15-node wedge elements and 44,827 nodes. (d) Zoom of the FE mesh near the left bottom corner, i.e., at x = 0.

Three material gradation cases are considered for the dynamic analysis of the three-point bending beam: a homogeneous beam (Homog, E₂ = E₁), a beam stiffer at the impacted surface (StiffTop, E₂ > E₁), and a beam softer at the impacted surface (StiffBot, E₂ < E₁), 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³ and 2000 kg/cm³. 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.,

$$\frac{E}{\rho }=\frac{1}{W}{\int }_0^W\frac{E\left(y\right)}{\rho \left(y\right)}\mathit{dy}$$E11

Poisson’s ratio is 0.33 for all the three cases. The average dilatational wave speed (C_d(avg)) is defined as

$$C_d\left(\mathit{avg}\right)=\frac{1}{W}{\int }_0^W{C}_d\left(y\right)\mathit{dy}$$E12

where C_d is calculated assuming plane stress behavior given by

$$C_d^2=\frac{E\left(y\right)}{\left(1-{\nu }^2\left(y\right)\right)\rho \left(y\right)}$$E13

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′ = 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 3.

Stress history σ_x at location P (0, 0.2 W) for homogeneous and graded beams with linearly varying E and ρ, subjected to velocity impact (V = 5 m/s) at the top.

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 4.

Stress history σ_y at location P (0, 0.2 W) for homogeneous and graded beams with linearly varying E and ρ, subjected to velocity impact (V = 5 m/s) at the top.

Figure 5 shows the 2D contour plot of the stress component σ_x at a specific time instant of 90 μs (t′ ∼ 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.

Figure 5.

Stress contour (σ_x MPa) at t = 90 μs due to impact velocity of 5 m/s. (a) StiffBot beam (E₂ < E₁), (b) homogeneous beam (E₂ = E₁), and (c) StiffTop beam (E₂ > E₁).

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.

Figure 6.

Stress history σ_x at location P(0, 0.2 W) for homogeneous and graded beams subjected to two different velocities.

Figure 7.

Stress history σ_y at location P(0, 0.2 W) for homogeneous and graded beams subjected to two different velocities.

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).

Figure 8.

Stress history σ_x at location P(0, 0.2 W) for homogeneous and graded beams subjected to impact velocity of 10 m/s as a line load and point load (see the insert). Thick and thin lines correspond to line load and point load, respectively. Solid, dashed, and dash-dot lines correspond to StiffTop, Homog, and StiffBot beams.

Figure 9.

Stress history σ_y at location P(0, 0.2 W) for homogeneous and graded beams subjected to impact velocity of 10 m/s as a line load and point load (see the insert). Thick and thin lines correspond to line load and point load, respectively. Solid, dashed, and dash-dot lines correspond to StiffTop, Homog, and StiffBot beams.

6. 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.

Figure 10.

(a) Schematic of the 3D bar and (b) step load.

Material	E (GPa)	ν	ρ (kg/m^3)	Cd (m/s)
Steel	210	0.31	7800	6109
Alumina	390	0.22	3950	10,617

Table 1.

Material properties of steel and alumina.

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]

$$u\left(x,t\right)=\frac{P_0}{\mathit{EA}}x-\frac{8P_{0}L}{{\pi }^2\mathit{EA}}\sum _{n=1}^{\infty }\frac{{\left(-1\right)}^{n-1}}{{\left(2n-1\right)}^2}sin\left(\left(2n-1\right)\frac{\mathit{\pi x}}{2L}\right)cos\left(\left(2n-1\right)\frac{\mathit{\pi c}}{2L}t\right)$$E14

where P₀, 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)². 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.

Figure 11.

Normalized longitudinal stress history at three locations (see the insert) of the homogeneous bar subjected to step loading using 1D analytical solution given in Eq. (14).

Figure 12.

Normalized longitudinal stress history at three locations (see the insert) of the homogeneous bar subjected to step loading obtained using finite element analysis.

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.

Figure 13.

Normalized longitudinal stress history of six points (see the insert) on a graded bar (z = 0.025 m) subjected to step loading. Solid and dashed lines indicate points at x = 0 and x = 0.5 L, respectively. Thick, intermediate-thick, and thin lines indicate alumina-rich side, midplane, and steel-rich side points, respectively.

Figure 14.

Normalized longitudinal stress history of six points (see the insert) on a graded bar (z = 0 m) subjected to step loading. Solid and dashed lines indicate points at x = 0 and x = 0.5 L, respectively. Thick, intermediate-thick, and thin lines indicate alumina-rich side, midplane, and steel-rich side points, respectively.

7. 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.