The variation in elastic properties inherent to FGMs is implemented at the element level by making use of user subroutines within the commercial finite element package ABAQUS. The graded elements described in
Section 2 can be readily implemented by using a USDFLD subroutine, for a Gauss points-based implementation, or a UFIELD subroutine, for a nodal-based graded element. In addition, as discussed above, the temperature can be used as an auxiliary field to effectively implement a generalized isoparametric graded element. The user must provide the material properties as a function of a user defined field (or temperature). Then, a suitable field (or temperature distribution) is defined to match the material property variation desired. A direct comparison between the two approaches in terms of computational time is hindered by their different implementations; nevertheless, the sampling of material properties at integration points or nodes is achieved at a negligible computational cost. The performance of different types of graded elements will be benchmarked by considering a Gauss point-based implementation (
Section 2.1) and, via temperature dependent properties, a generalized isoparametric approach (
Section 2.2). For the sake of simplicity, we will consider bi-dimensional problems and four quadrilateral element types (see
Figure 3): linear elements with reduced integration (Q4R), linear elements with full integration (Q4), quadratic elements with reduced integration (Q8R), and quadratic elements with full integration (Q8).
Three plane problems for which analytical solutions can be obtained will be addressed, as shown in
Figure 4. Young’s modulus will be varied along the
x-direction and Poisson’s ratio will be assumed to be constant. We assume plane stress conditions. As sketched in
Figure 4, the three boundary value problems considered involve a functionally graded plate being subjected to (i) uniform displacement perpendicular to the material gradient direction, (ii) uniform traction perpendicular to the material gradient direction, and (iii) uniform traction in the direction parallel to material gradation.
In all cases, we assume that the Young’s modulus varies exponentially as
with
and
being material constants. This choice is motivated by the existence of analytical solutions for functionally graded solids exhibiting an exponential variation of the elastic properties; see, for example, Refs. [
13,
22]. Many other functions have been employed in the literature (see Ref. [
23] for a review), but our choice is appropriate for our aim: comparing on equal footing different graded finite element implementations. Consistent units are throughout the manuscript and, therefore, units will be omitted subsequently. We consider a width
, a total height of
, and choose
and
so as to vary
E gradually in the
x-direction from
to
.
3.1. Uniform Displacement Perpendicular to the Material Gradient Direction
Consider first the case of an FGM plate subjected to a remote strain
, where
denotes the displacement in the remote boundary and
h denotes half the height of the plate. A constant Poisson’s ratio of
throughout the plate is assumed. The relevant stress component is given by
This analytical solution is compared in
Figure 5 with the finite element results obtained for the element types and graded element formulations described above. A remote displacement of
is prescribed in all cases. As shown in the insets of
Figure 5, a uniform mesh of 4 by 12 elements is employed.
Results reveal differences between the different types of graded elements. Consider first the case of the linear element with full integration Q4—see
Figure 5a. The Gauss integration point-based approach accurately captures the material gradient depicted by the analytical solution. In fact, the numerical result is exact at the integration points since the displacement field is linear; a single Q4 element will suffice to capture the FGM response. However, the generalized isoparametric formulation (via temperature-dependent properties) exhibits a step-type variation with constant stress in each element. This behaviour is inherently related to how ABAQUS interpolates nodal temperature values. As many other finite element packages, ABAQUS interpolates nodal temperatures with shape functions that are one order lower than those used for the displacements, so as to obtain an equivalent distribution of mechanical and thermal strains. An average value of the temperature in the nodes is passed to the integration points when using linear elements and a linear variation is assumed in quadratic elements. In agreement with expectations, the use of linear elements with reduced integration (Q4R, see
Figure 5b) exhibits the response inherent to homogeneous elements for both cases. However, the constant value of
attained in each element depends on the implementation approach. The integration point-based scheme computes the exact
at the element centroid, where
E is sampled. On the other hand, the nodal-based approach averages nodal temperatures, introducing a source of error when the elastic properties vary in a nonlinear manner. The use of quadratic elements leads to a good agreement with the analytical solution for both schemes, although differences are observed (
Figure 5c,d). We plot the error obtained with both schemes for the element Q8 in
Figure 6. It is evident that the integration point-based approach reproduces the analytical result more accurately.
Further insight is gained by reproducing the analysis with a single element in the
x-direction—see
Figure 7. Inspection of
Figure 7a reveals substantial differences between Gauss points-based and nodal-based implementations. The former reproduces precisely the analytical stress distribution, with the exact result being obtained at the integration points. Contrarily, the approximation is much poorer when the material gradient is implemented via the temperature. Differences are particularly significant at the left side of the specimen, where the error is on the order of 50%.
Remarkably,
negative stresses are predicted for
. These non-physical compressive stresses arise as a consequence of the particularities of ABAQUS’ criterion for interpolating nodal temperatures, which does not correspond to the linear interpolation outlined in
Figure 2. In ABAQUS, the nodal temperature values are multiplied by certain weights, such that the temperature
T at an integration point
i is given by
Here,
is the temperature in node
j,
the weight associated with the nodal temperature
j and integration point
i, and
n and
m respectively denote the total number of nodes and integration points. The specific values of
depend on a number of numerical considerations and can be easily obtained by means of a one-element model. This criterion is motivated by numerical convergence in thermomechanical problems, as it smoothens localized temperature peaks. The resulting variation in the elastic properties within the element is shown in
Figure 7b. Differences with a direct linear interpolation are evident. The weighting procedure implemented in ABAQUS brings non-physical values of
E, but it shows a better agreement with the material gradation profile at the Gauss integration points. In turn, this better approximation of
reduces the error in the computation of the stresses, as quantified in
Figure 7c as a function of the number of elements. The log–log plot of
Figure 7c shows that the weighted interpolation of ABAQUS exhibits a smaller maximum error in the computation of
at the Gauss points, as well as a faster convergence rate. Consequently, the error intrinsic to a temperature-based graded element is magnified if a standard linear interpolation is used and, therefore, the conclusions of the present study are even more relevant to finite element codes that employ non-weighted interpolations of nodal temperatures. Recall that the integration point-based scheme presented in
Section 2.1 captures the analytical solution at the Gauss points with a single element.
3.2. Uniform Traction Perpendicular to the Material Gradient Direction
Consider now the case where the remote load is prescribed as a traction perpendicular to the elastic gradient—see
Figure 4c. The Dirichlet boundary conditions of the problem read
In the case of a plate with infinite height, the only non-zero component of the Cauchy stress tensor is
. Following Refs. [
13,
22], a membrane resultant
N along the
line can be defined as a function of the remote stress
and the width,
The compatibility condition
requires the strain component to be of the form
and, consequently, one can readily obtain the stress field by considering the exponential elastic modulus variation assumed (
12) and making use of Hooke’s law as
The coefficients
A and
B are obtained by solving
such that
The displacement solution can be readily obtained by making use of the strain-displacement relations and applying the boundary conditions (
15) and (
16),
The analytical solution in the middle line
is compared with the numerical predictions for
. Results are shown in
Figure 8 for a uniform mesh of 48 plane stress elements.
Differences between the integration point-based scheme and the nodal/temperature-based implementation are particularly significant for the case of linear elements with full integration (Q4,
Figure 8a). Sampling the material gradient directly at the Gauss points leads to a good agreement with the analytical solution. However, using temperature-based properties renders the homogeneous element solution. Both the analytical and Gauss point-based solutions show an increasing
along
x within those elements close to the left edge. Contrarily, the inverse response is observed when using temperature-dependent properties, as the strain field decreases with
x and
E is constant element-wise. The Q4 element predicts in all cases a linear variation of
within each element for the quadratic displacement solution under consideration. On the other hand, identical predictions between graded element schemes are obtained when using linear elements with reduced integration (Q4R,
Figure 8b). This is unlike the case of a prescribed displacement (see
Figure 5b), as the error in the approximation of the strain field
(exact at the integration points only for the Gauss points-based scheme) is compensated. Quadratic elements (Q8 and Q8R) introduce an element-wise variation of
E in both approaches and, therefore, differences appear to be smaller than in linear elements—see
Figure 8c,d. The error in the approximation is shown in
Figure 9 for the Q8 element case. When using full integration, the Gauss point-based approach outperforms the temperature-based, generalized isoparameteric graded element.
Further insight is gained by analysing a very coarse mesh with a single element in the
x-direction; results are shown in
Figure 10. Regarding the stress (
Figure 10a), the prediction obtained from the integration point-based implementation of graded elements agrees well with the analytical solution, being exact at the Gauss points (symbols). On the other hand, the approximation via a nodal/temperature-based approach introduces a significant source of error (larger than 20% in all the integration points). Moreover, when the load is prescribed as a traction, the approximation of the material gradient also influences the strain field, see Equations (
8) and (
9). As shown in
Figure 10b, a better approximation is attained if the material properties are sampled directly at the Gauss points. Both schemes differ at the edges with the analytical solution for a plate of infinite height.
3.3. Uniform Traction Parallel to the Material Gradient Direction
The last case study involves a functionally graded plate subjected to traction in the
x-direction, parallel to the elastic gradient, see
Figure 4d. Under those conditions, the normal stress component equals the applied stress
if Poisson’s ratio is made equal to zero,
. The results obtained for a uniform mesh of 75 plane stress elements are shown in
Figure 11.
Contrarily to what has been observed so far, the Q4 results (
Figure 11a) show that the nodal/temperature-based approach outperforms the integration point-based counterpart. Hooke’s law requires the strains to vary according to an inverse exponential distribution to obtain a constant stress for an exponentially varying
E. However, linear elements predict a constant strain field and, consequently, an effectively homogeneous element will predict a constant stress. This trend is inverted for the case of a quadratic element with full integration (Q8). As shown in
Figure 11b, again the use of a Gauss point-based graded element formulation approximates the analytical solution better. The results pertaining to reduced integration elements (Q4R and Q8R), not shown for brevity, reveal a perfect agreement with the analytical solution in all cases. Thus, reduced integration improves precision in this specific case study as the resulting stress is constant—either because there is a single integration point (Q4R), leading to a constant
and
E, or because both
and
E are element-wise linear (Q8R).