Analysis of Stress Concentration in Functionally Graded Plates with Linearly Increasing Young’s Modulus

Abstract

In this article, the strain and stress analyses of functionally graded plates with circular holes that are subject to a uniaxial far-field traction load are analytically considered. The Young’s modulus is assumed to vary linearly along the radial direction around the hole. The adoption of such a type of inhomogeneity variation can be justified as follows. Firstly, and among all the possible variations of stiffness, the linear one is indeed the simplest inhomogeneity distribution. Surprisingly however, according to our knowledge extent, the associated elastic fields were not yet addressed in the literature. Secondly, a linearly varying stiffness could reasonably imply a remarkable advantage from a technological point of view. In fact, unlike nonlinearly varying stiffness plates, manufacturing routes are only required to handle constant variations throughout the radial domain. After recalling the basic equations for plane stress elasticity, the displacement, strain, and stress fields around the hole were numerically tackled and discussed for different stiffness ratios. A comparison was also carried out with other Young’s modulus distributions that have been commonly employed in the literature.

1. Introduction

The stress evaluation in plates endowed with geometrical discontinuities constitutes one of the major problems in mechanical engineering. Kirsch [1] analytically assessed the stress concentration factor at the rim of a circular hole under the assumption of a homogeneous, linearly elastic and isotropic material under a far-field, in-plane traction load. He found that the hoop stress, at the rim of the hole, comes out to be identically three times the applied load when the dimension of the hole is very small compared to the dimensions of the plate. Since then, the interest in reducing such a factor by abandoning the above assumptions has been considerably increasing. Analytical and numerical approaches have yielded results for both infinite and finite plates. For instance, an overview of the analytic methods describing the elastic stresses in two-dimensional flat anisotropic plates with notches is presented in [2]. In [3], the authors report studies on the stress concentration arising near to holes and different cutouts when the panels are made of fiber composite materials. The authors in [4] derived closed-form solutions for circularly notched composite laminates that were subject to tensile loading. The mechanical behavior of composite finite laminates endowed with open holes was assessed by analytical methods in [5]. In [6], the authors analyzed the stress concentration at the edge of a circular hole in sheets with different widths by means of the finite element method.

On the other hand, the increasing interest in the so-called functionally graded materials has been gaining considerable interest by researchers, since the adoption of such a material concept could lead to remarkable stress and strain enhancements in several mechanical models. For instance, the vibrations of thin, microstructured, functionally graded plates were studied in [7,8] by means of the so-called tolerance modeling method. In particular, free vibration frequencies bands were analyzed, and the effect of various distribution functions on the plate properties was shown. In [9], the thermomechanical stresses occurring in functionally graded rotating hollow disks were analytically computed and numerically optimized for a particular gradation law along the radial direction of the disk. Lightweight and economical designs for functionally graded rotating disks of a fixed mass based on the limitation of their angular speeds before material yielding can be found in [10]. The mass values, strains, and stresses of rotating disks subject to thermal loads were numerically enhanced within a multiobjective optimization study in [11]. In [12], the thermal and mechanical stresses occurring in functionally graded pressurized cylindrical vessels with spherical domes were computed and optimized using a commercial finite-element-based software. The flexural, free vibration, and buckling responses of two-dimensional functionally graded curved beams whose material properties in both directions were varied according to power law distributions were numerically studied in [13]. In [14], the maximization of the sound transmission loss across a functionally graded material in the form of a cylindrical shell was conducted using a genetic algorithm and compared to the responses associated with several traditional materials.

As far as plates with geometrical discontinuities (e.g., a circular hole) are concerned, the literature offers substantial works where the stress analysis has been carried out either analytically or numerically—once the homogeneity assumption is dropped. For instance, the effect of the material inhomogeneity on the stress field due to circular holes was studied in [15] by means of the finite element method. In particular, a considerable reduction in the stress concentration factor was obtained by allowing the Young’s modulus to vary in a power law or exponential fashion. In [16], the stress concentration of an infinite plate subjected to uniform biaxial tension and pure shear was analytically obtained by considering an exponential variation of Young’s modulus throughout the radius. In [17], closed-form solutions for the stress fields of functionally graded plates subject to uniform far-field tensile traction were derived. The elastic response of a functionally graded annular ring inserted around the hole of a homogeneous plate was derived analytically in [18,19] under different far-field loading conditions. The stress concentration factor measured for exponentially graded plates with circular holes and subject to biaxial loads was studied in [20] by means of the boundary element method. The main results of these works can be summarized by the following two statements:

• The stress concentration factor decreases only when the stiffness progressively increases away from the circular hole;
• The effect of the variation in the Poisson’s ratio on the stress distribution is negligible.

These two findings entail a remarkable rule of thumb for designers and stress analysts.

Despite the enhanced stress results arising from the relaxation of the homogeneity assumption, it is surprising that studies on the stress concentration behavior in functionally graded plates with a circular hole whose Young’s modulus linearly increases are not considered in the literature. Although a number of works that adopted a linear (or piecewise linear) stiffness distribution can be found in the literature for different mechanical models, e.g., [21,22,23,24], such a distribution has been so far neglected as far as stress concentration problems are concerned. Supposedly, a linearly increasing stiffness is the simplest variation form among all of the possible variations and could present an advantage from a manufacturing viewpoint with respect to its nonlinear counterparts. According to [25], this could be one of the desired characteristics when manufacturing such materials, since steep variations in the stiffness along the radial domain could be, reasonably, more costly and more difficult to obtain.

The article is organized as follows. Section 2 recalls the governing equations for the elastic behavior of the functionally graded plate. A boundary value problem is obtained and solved numerically in Section 3. The numerical solutions for the displacement, strain, and stress fields are presented in graphical form and discussed in Section 4. Finally, conclusions are drawn in Section 5.

2. Governing Equations

Consider a locally isotropic, linearly elastic, and functionally graded infinite plate with a circular hole of radius a and that is subject to a far-field uniaxial load ${\sigma }_0$, as shown in Figure 1a, where the generic point P is described by the polar coordinate system $(r,\theta )$, whose origin is at the center of the circular hole.

Let SQ and MN denote the horizontal and vertical lines, respectively, associated with the polar angles $\theta =0$ and $\theta =\pi /2$, respectively. Moreover, let the inhomogeneity be described by a radially linear variation in the Young’s modulus, of slope e, which is given by

$$E\left(r\right)=E_a+e(r-a) ,$$

(1)

where $E_a$ denotes the Young’s modulus at $r=a$.

By letting the thickness of the plate be sufficiently small to the point that the stress state results are two-dimensional (plane stress condition) and considering no body forces, the equilibrium equations in terms of the plane stresses in the polar coordinates read as [26]

$$\left\{\begin{array}{c}{\displaystyle \frac{\partial {\sigma }_r(r,\theta )}{\partial r}+\frac{1}{r}\frac{\partial {\tau }_{r\theta }(r,\theta )}{\partial \theta }+\frac{{\sigma }_r(r,\theta )-{\sigma }_{\theta }(r,\theta )}{r}=0 ,}\hfill \\ {\displaystyle \frac{\partial {\tau }_{r\theta }(r,\theta )}{\partial r}+\frac{1}{r}\frac{\partial {\sigma }_{\theta }(r,\theta )}{\partial \theta }+\frac{2}{r}{\tau }_{r\theta }(r,\theta )=0 ,}\hfill \end{array}\right.$$

(2)

where ${\sigma }_r$, ${\sigma }_{\theta }$, and ${\tau }_{r\theta }$ are, respectively, the radial, hoop, and shear stresses, whose orientations are shown in Figure 1b. These elastic stresses are, in turn, linked to the corresponding strains by the plane stress constitutive relations, namely [26],

$$\left\{\begin{array}{c}{\displaystyle {\sigma }_r(r,\theta )=\frac{E\left(r\right)}{1-{\nu }^2}\left({\epsilon }_r(r,\theta )+\nu {\epsilon }_{\theta }(r,\theta )\right) ,}\hfill \\ {\displaystyle {\sigma }_{\theta }(r,\theta )=\frac{E\left(r\right)}{1-{\nu }^2}\left({\epsilon }_{\theta }(r,\theta )+\nu {\epsilon }_r(r,\theta )\right) ,}\hfill \\ {\displaystyle {\tau }_{r\theta }(r,\theta )=\frac{E\left(r\right)}{2(1+\nu )}{\gamma }_{r\theta }(r,\theta ) ,}\hfill \end{array}\right.$$

(3)

where $\nu$ is Poisson’s ratio (which is supposed to be constant for simplicity) and ${\epsilon }_r$, ${\epsilon }_{\theta }$, and ${\gamma }_{r\theta }$ are the radial, hoop, and shear strains, respectively. From Figure 1a, the boundary conditions can be described as follows:

• At the hole periphery, both the radial and shear stresses are zero for all values of $\theta$;
• At a sufficiently large distance from the hole, the radial and shear stresses on the MN approach zero;
• At a sufficiently large distance from the hole, the radial and shear stresses on the SQ tend to ${\sigma }_0$ and to zero, respectively.

Furthermore, strains are linked to the plane displacement field through the following kinematic relations in polar coordinates [26]:

$$\left\{\begin{array}{c}{\displaystyle {\epsilon }_r(r,\theta )=\frac{\partial u(r,\theta )}{\partial r} ,}\hfill \\ {\displaystyle {\epsilon }_{\theta }(r,\theta )=\frac{u(r,\theta )}{r}+\frac{1}{r}\frac{\partial v(r,\theta )}{\partial \theta } ,}\hfill \\ {\displaystyle {\gamma }_{r\theta }(r,\theta )=\frac{1}{r}\frac{\partial u(r,\theta )}{\partial \theta }+\frac{\partial v(r,\theta )}{\partial r}-\frac{v(r,\theta )}{r} ,}\hfill \end{array}\right.$$

(4)

where u and v are the radial and hoop components of the displacement field, respectively. It is worthwhile to mention that these displacement functions have to be sagaciously chosen so that the elastic strains obey the following compatibility equation [26]:

$$\frac{{\partial }^2{\epsilon }_{\theta }(r,\theta )}{\partial r^2}+\frac{1}{r^2}\frac{{\partial }^2{\epsilon }_r(r,\theta )}{\partial {\theta }^2}+\frac{2}{r}\frac{\partial {\epsilon }_{\theta }(r,\theta )}{\partial r}-\frac{1}{r}\frac{\partial {\epsilon }_r(r,\theta )}{\partial r}=\frac{1}{r}\frac{{\partial }^2{\gamma }_{r\theta }(r,\theta )}{\partial r\partial \theta }+\frac{1}{r^2}\frac{\partial {\gamma }_{r\theta }(r,\theta )}{\partial \theta } .$$

(5)

If the Young’s modulus is assumed to be constant in Equations (1)–(4) (namely, homogeneous plates), the elastic stresses can be determined analytically. More precisely, by splitting the elastic problem into a biaxial (axisymmetric) problem and a pure shear problem, the superposition principle and the introduction of the Airy stress function lead to the well-known Kirsch stress field [26]:

$$\left\{\begin{array}{c}{\displaystyle {\sigma }_r(r,\theta )=\frac{{\sigma }_0}{2}\left(1-\frac{a^2}{r^2}\right)+\frac{{\sigma }_0}{2}\left(1-\frac{4a^2}{r^2}+\frac{3a^4}{r^4}\right)\cos 2\theta ,}\hfill \\ {\displaystyle {\sigma }_{\theta }(r,\theta )=\frac{{\sigma }_0}{2}\left(1+\frac{a^2}{r^2}\right)-\frac{{\sigma }_0}{2}\left(1+\frac{3a^4}{r^4}\right)\cos 2\theta ,}\hfill \\ {\displaystyle {\tau }_{r\theta }(r,\theta )=-\frac{{\sigma }_0}{2}\left(1+\frac{2a^2}{r^2}-\frac{3a^4}{r^4}\right)\sin 2\theta .}\hfill \end{array}\right.$$

(6)

Equation (6) forms one of the few closed-form solutions in the theory of elasticity, whose expressions are twofold. Firstly, it demonstrates that if r and $\theta$ tend to a and $\pi /2$, respectively, the stress concentration factor, namely, the ratio between the hoop stress and the applied load, is identically three. Another important consequence consists of the derivation of the expression for the far-field radial and shear stresses at a sufficiently large distance from the hole by making r tend to $+\infty$, namely,

$$\left\{\begin{array}{c}{\displaystyle \underset{r\to +\infty }{\lim }{\sigma }_r(r,\theta )=\frac{{\sigma }_0}{2}(1+\cos 2\theta ) ,}\hfill \\ {\displaystyle \underset{r\to +\infty }{\lim }{\tau }_{r\theta }(r,\theta )=-\frac{{\sigma }_0}{2}\sin 2\theta .}\hfill \end{array}\right.$$

(7)

Equation (7) not only transcribes the boundary conditions as a function of $\theta$, but also allows for a transformation of the rectangular domain into a circular one as long as a limit radius, hereinafter denoted by A, is provided (see Figure 1c). The latter can be defined by referring to a fixed threshold of the discrepancy of the stress concentration factor associated with finite and infinite homogeneous plates. For instance, in [27], the solution to the problem for radially inhomogeneous material was carried out with a limit radius that was ten times the radius of the hole. Such a value, for homogeneous plates, yields a discrepancy in the stress concentration factor of less than 3%. In the present study, we took as the limit radius the double of such a value.

Next, the analytical derivation for the elastic displacement field, and therefore the strains and stresses from Equations (3) and (4), is shown by adopting a set of compatible displacement functions, in the sense of Equation (5), which is borrowed from the generalized Michell solution in the polar coordinates [28]. Three coupled differential equations are obtained, whose analytical treatment is cumbersome, thus enforcing one to resort to a numerical tool when considering the resulting boundary value problem (BVP).

3. Solution of the Problem

Let the set of radial and hoop displacements be a function of r and $\theta$, which is given by [27]:

$$\left\{\begin{array}{c}u(r,\theta )=\overline{u}\left(r\right)+\widehat{u}\left(r\right)\cos 2\theta ,\hfill \\ v(r,\theta )=\widehat{v}\left(r\right)\sin 2\theta .\hfill \end{array}\right.$$

(8)

It is easy to show that the radial displacements along the MN and SQ are given by $\overline{u}\left(r\right)-\widehat{u}\left(r\right)$ and $\overline{u}\left(r\right)+\widehat{u}\left(r\right)$, respectively. Moreover, it is required that $\overline{u}\left(A\right)-\widehat{u}\left(A\right)$ be strictly negative, thus exhibiting the well-known Poisson’s effect. Based on the kinematic relations (4), the elastic strains read as follows:

$$\left\{\begin{array}{c}{\displaystyle {\epsilon }_r={\overline{u}}^{\prime }+{\widehat{u}}^{\prime }\cos 2\theta ,}\hfill \\ {\displaystyle {\epsilon }_{\theta }=\frac{\overline{u}}{r}+\left(\frac{\widehat{u}}{r}+2\frac{\widehat{v}}{r}\right)\cos 2\theta ,}\hfill \\ {\displaystyle {\gamma }_{r\theta }=-\left(\frac{\widehat{v}}{r}-{\widehat{v}}^{\prime }+2\frac{\widehat{u}}{r}\right)\sin 2\theta ,}\hfill \end{array}\right.$$

(9)

where the dependence on r and $\theta$ are omitted for the sake of light notation, and a prime denotes a first derivative with respect to r. It is worth noting that Equation (9) satisfies the compatibility condition (5) identically, thus suggesting promising candidates for the description of the displacement field. Hence, the elastic stress field reads as follows:

$$\left\{\begin{array}{c}{\displaystyle {\sigma }_r=\frac{E}{1-{\nu }^2}\left[{\overline{u}}^{\prime }+\nu \frac{\overline{u}}{r}+\left({\widehat{u}}^{\prime }+\nu \frac{\widehat{u}}{r}+2\nu \frac{\widehat{v}}{r}\right)\cos 2\theta \right] ,}\hfill \\ {\displaystyle {\sigma }_{\theta }=\frac{E}{1-{\nu }^2}\left[\frac{\overline{u}}{r}+\nu {\overline{u}}^{\prime }+\left(\frac{\widehat{u}}{r}+2\frac{\widehat{v}}{r}+\nu {\widehat{u}}^{\prime }\right)\cos 2\theta \right] ,}\hfill \\ {\displaystyle {\tau }_{r\theta }=-\frac{E}{2(1+\nu )}\left(\frac{\widehat{v}}{r}-{\widehat{v}}^{\prime }+\frac{2\widehat{u}}{r}\right)\sin 2\theta ,}\hfill \end{array}\right.$$

(10)

whose substitution in the equilibrium relations of (2) yields the following three second-order linear differential and coupled equations:

$$\left\{\begin{array}{c}{\widehat{u}}^{″}+{\alpha }_1{\widehat{u}}^{\prime }+{\alpha }_2{\widehat{v}}^{\prime }+{\alpha }_3\widehat{u}+{\alpha }_4\widehat{v}=0 ,\hfill \\ {\overline{u}}^{″}+{\alpha }_1{\overline{u}}^{\prime }+{\alpha }_6\overline{u}=0 ,\hfill \\ {\widehat{v}}^{″}+{\alpha }_1{\widehat{v}}^{\prime }+{\alpha }_7{\widehat{u}}^{\prime }+{\alpha }_5\widehat{v}+{\alpha }_8\widehat{u}=0 ,\hfill \end{array}\right.$$

(11)

where ${\alpha }_i$ ($i=1,2,\dots ,8$) are functions that depend on the distribution of the mechanical properties along the radial direction. Precisely, they are given by the following:

$$\left\{\begin{array}{c}{\alpha }_1=\frac{r{E}^{\prime }+E}{rE} ,\hfill \\ {\alpha }_2=\frac{\nu +1}{r} ,\hfill \\ {\alpha }_3=-\frac{3-2\nu }{r^2}+\frac{\nu E^{\prime }}{rE} ,\hfill \\ {\alpha }_4=-\frac{3-\nu }{r^2}+\frac{2\nu E^{\prime }}{rE} ,\hfill \\ {\alpha }_5=-\frac{E^{\prime }}{rE}+\frac{9+\nu }{(\nu -1)r^2} ,\hfill \\ {\alpha }_6=\frac{\nu E^{\prime }}{rE}-\frac{1}{r^2} ,\hfill \\ {\alpha }_7=\frac{2\nu +2}{(\nu -1)r} ,\hfill \\ {\alpha }_8=-\frac{2E^{\prime }}{rE}+\frac{6-2\nu }{(\nu -1)r^2}\hfill \end{array}\right.$$

(12)

Finally, as derived from Equation (10), the boundary conditions at the hole periphery, as well as at the limit radius, can be written in terms of $\overline{u}$, $\widehat{u}$, and $\widehat{v}$ as follows:

$$\left\{\begin{array}{c}{\displaystyle a{\overline{u}}^{\prime }\left(a\right)+\nu \overline{u}\left(a\right)=0 ,}\hfill \\ {\displaystyle a{\widehat{u}}^{\prime }\left(a\right)+\nu \widehat{u}\left(a\right)+2\nu \widehat{v}\left(a\right)=0 ,}\hfill \\ {\displaystyle a\widehat{v}\left(a\right)+a{\widehat{v}}^{\prime }\left(a\right)-2\widehat{u}\left(a\right)=0 ,}\hfill \\ {\displaystyle A{\overline{u}}^{\prime }\left(A\right)+\nu \overline{u}\left(A\right)-\frac{{\sigma }_0}{2}\frac{1-{\nu }^2}{E_A}=0 ,}\hfill \\ {\displaystyle A{\widehat{u}}^{\prime }\left(A\right)+\nu \widehat{u}\left(A\right)+2\nu \widehat{v}\left(A\right)-\frac{{\sigma }_0}{2}\frac{1-{\nu }^2}{E_A}=0 ,}\hfill \\ {\displaystyle A\widehat{v}\left(A\right)+A{\widehat{v}}^{\prime }\left(A\right)-2\widehat{u}\left(A\right)-\frac{{\sigma }_0(1+\nu )A}{E_A}=0 ,}\hfill \end{array}\right.$$

(13)

where $E_A$ is the Young’s modulus at the limit radius, which is evaluated by substituting $r=A$ in Equation (1).

The system of differential equations in Equation (11) and the boundary conditions (13) form a BVP, whose solution can be carried out numerically using a dedicated solver [29]. More precisely, the differential equations in Equation (11) are integrated on a mesh formed by N points, which divides the interval of integration into subintervals. A mesh convergence study led us to take $N=200$. The collocation polynomial approximating the solutions provides C¹ continuous solutions. The solver determines a numerical solution by solving a global system of algebraic equations resulting from the boundary conditions (13), as well as the collocation conditions imposed on all the subintervals. The solver then estimates the error of the numerical solution on each subinterval; if the solution does not satisfy the tolerance criteria, the solver adapts the mesh and repeats the process.

4. Numerical Results

In this section, a numerical solution for the above mentioned BVP is presented and shown in graphical form. For simplicity, let the Young’s modulus linear distributions be described by the ratio $E_A/E_a$. Consider four instances of this ratio: One of them is associated with the value $E_A/E_a=1$ referring to the homogeneous case for the mere sake of comparison, while the others are 5, 10, and 20. Unless explicitly specified, let the dashed lines refer to numerical solutions associated with the homogeneous case. The derivation of the stresses in units of MPa requires one to express the stiffness in the same unit, the radius of the hole a, the limit radius A in mm, and the slope e in Equation (1), which is given by

$$e=\frac{E_A-E_a}{A-a} ,$$

(14)

in MPa/mm. However, displacements, strains, and stresses are hereinafter normalized with respect to $A{\sigma }_0/E_a$, ${\sigma }_0/E_a$, and ${\sigma }_0$, respectively.

For clarity of exposition, the radial displacement fields along the MN and SQ are shown in Figure 2a,b respectively, which show the quantities $\overline{u}-\widehat{u}$ and $\overline{u}+\widehat{u}$, respectively. It is worth noting that the radial displacement along the MN was negative along the lines and up to the limit radius (whose values are highlighted by the scatter points in Figure 2a), thus confirming the Poisson’s effect when the plate is subject to a far-field traction load. As far as the displacement along the SQ at the limit radius is concerned, it is not difficult to notice that it decreased as long as the stiffness ratio $E_A/E_a$ reached higher values, which is in agreement with the physical idea of the corresponding continuously hardened plates.

Figure 3a,b show the normalized radial and hoop stresses, respectively, along the MN line.

It is not difficult to notice that the boundary conditions of the radial stress were satisfied at the hole periphery and at the limit radius. Furthermore, the peak value was progressively reduced as the stiffness ratio increased, and the prefixed radial stress at the limit radius was reached from the left with higher slopes. As far as the hoop stresses along the MN line are concerned, Figure 3b highlights two important aspects. The first one consists of the remarkable reduction in the maximum stress with respect to the homogeneous case. In fact, the hoop stresses associated with linear Young’s modulus distributions showed less stress concentration factors at the rim of the hole. To give some numerical values, the stress concentration factor decreased by approximately 39%, 58%, and 73% for the three considered stiffness ratios, respectively. The second aspect regards the stress uniformity along the MN line. In particular, the stress concentration factor decreased at the hole periphery, while a maximum value was recorded at some point along the radial direction. This inevitable stress rise is a consequence of the conservation of the integral of the hoop stress resulting from the equilibrium along the vertical MN line. Nevertheless, the rate of change of the maximum value of the hoop stress was considerably lower than the rate of change of the stress concentration factor at the rim of the hole, thus yielding lower variance values for the stress profiles along the MN line.

5. Discussion

The above considerations shed the light on the existence of a particular value of the stiffness ratio, say $\overline{E_A/E_a}$, beyond which the maximum value for the hoop stress does not occur at the rim of the hole, but elsewhere (see Figure 4).

For the considered limit radius, such a value turned out to be 11. For halved and doubled limit radii, the $\overline{E_A/E_a}$ was 9 and 13, respectively, and the discrepancy between the corresponding peak stress values was less than 5%. Figure 4 shows both the hoop stress at the rim of the circular hole and its maximum value as functions of $E_A/E_a$. In agreement with Figure 3b, it is easy to show that both of these stresses were identical up to the aforementioned value, namely, for a stiffness ratio close to 11, after which the maximum value for the hoop stress seemed to increase very marginally.

For the sake of validation of these numerical results, a numerical model for the plate has been developed by means of the finite element methods (FEM). Due to symmetrical load and geometrical considerations, the numerical model consisted of the quarter of the plate and employed plane stress elements. Symmetric boundary conditions and the uniaxial load have been correctly applied to the model. The radial direction has been discretized into radial strips, each of which is isotropic and homogeneous and has the same mechanical properties. Adjacent layers present different properties such that the resulting piecewise constant variation mimics the linear Young’s modulus distribution. Maximum hoop stresses and the stress concentration factors have been computed and reported as scatter points in Figure 4, which shows a remarkable fit with the numerical results obtained by the method described in Section 3.

Eventually, some remarks on the worthiness of the considered Young’s modulus distribution will be addressed. This can be assessed by fixing an instance for the stiffness ratio, e.g., $E_A/E_a=10$, and comparing its elastic performance with that associated with different laws, which are largely employed in the literature (Figure 5a). Sound comparison criteria could include the maximum hoop stress along the radial domain and determing whether this value occurs at the hole periphery. For instance, for the specified inhomogeneity parameters $\beta$ and n in the generalized power laws considered in [17], namely,

$$E\left(r\right)=E_A\left[1+\beta {\left(\frac{r}{a}\right)}^n\right] ,$$

(15)

the hoop stress at the rim of the hole was 40% of the applied load, while the peak stress was 73% more than the applied load (see the dotted line in Figure 5b). In [18], elastic stresses occurring in a functionally graded ring that was perfectly bonded to the circular hole were derived analytically when the stiffness in the ring was described by the power law,

$$E\left(r\right)=E_A{\left(\frac{r}{b}\right)}^m ,$$

(16)

where b is the radius of the ring (see the dahsed–dotted lines in Figure 5a). For fixed geometry and stiffness ratios, the grading factor m in Equation (16) can be computed as

$$m=\frac{\log \left(E_a/E_A\right)}{\log (a/b)} .$$

(17)

The hoop stress associated with two rings whose radii are three and five times the radius of the hole is shown in Figure 5b (dashed–dotted lines). Also here, in both cases, the maximum hoop stresses, whose values were, respectively, 54% and 36% more than the applied load, did not occur at the hole periphery. On the other hand, the maximum hoop stress for the linearly varying stiffness was 26% more than the applied load (solid line in Figure 5b), and this value decreased when the stiffness ratio assumed higher values (as shown in Figure 3b).

It is worth mentioning that most of the works available in the literature prefix the Young’s modulus distribution a priori and address, either analytically or numerically, the associated elastic problem. For instance, the solutions in the present work, as well as in [17,18], were obtained by heuristically choosing the parameters describing the stiffness throughout the plate. In fact, a proper choice of the parameters $\beta$ and n in (15) and $b/a$ and m in (16) led to enhanced values for the stress concentration factor, while in the present work, it was shown that there exists a slope e of the stiffness distribution that yields comparable results. However, we believe that a more intriguing approach should consist in finding the Young’s modulus distribution without prefixing its functional form beforehand. This approach could be introduced within an optimization problem aimed at reducing the stress concentration factor as much as possible. After discretizing the domain into a finite number of, e.g., radials strips, possible decision variables could include the Young’s modulus in each strip. Hence, the results of the present study can be useful while developing numerical codes, as these typically employ linear initial guesses before seeking the optimum solution.

6. Conclusions

The stress behavior of functionally graded plates with a circular hole and that are subject to a far-field uniaxial traction load has been addressed in this article. The Young’s modulus was allowed to vary linearly along the radial domain, while the Poisson’s coefficient was kept constant. A set of coupled linear differential equations describing the displacement, strain, and stress fields was analytically derived and numerically tackled for a few linear stiffness distributions that were associated with different slopes. Numerical solutions have been shown in graphical forms and discussed. Specifically, along the vertical MN line, the following concluding remarks are highlighted:

• Similar to the homogeneous case, a linearly varying stiffness distribution along the radial domain makes each radial strip become shortened and stretched along the radial and circumferential direction, respectively. These strain measures decrease in absolute value as the stiffness ratio increases.
• The peak value for the radial stress progressively decreases as the stiffness ratio increases.
• The stress concentration factor remarkably decreases as the stiffness ratio increases, thus confirming one of the most established findings in the literature.
• There exists a stiffness value beyond which the maximum hoop stress does not occur at the rim of the hole.

These findings have been validated by the finite element method and compared to the elastic performance associated with other stiffness distributions that have been commonly employed in the literature.